Path-Integral Formulation of Bosonic Markovian Open Quantum Dynamics with Monte Carlo stochastic trajectories using the Glauber-Sudarshan P, Wigner, and Husimi Q Functions and Hybrids

Abstract

The Monte Carlo (MC) trajectory sampling of stochastic differential equations (SDEs) based on the quasiprobabilities, such as the Glauber-Sudarshan P, Wigner, and Husimi Q functions, enables us to investigate bosonic open quantum many-body dynamics described by the Gorini-Kossakowski-Sudarshan-Lindblad (GKSL) equation. In this method, the MC samplings for the initial distribution and stochastic noises incorporate quantum fluctuations, and thus, we can go beyond the mean-field approximation. However, description using SDEs is possible only when the corresponding Fokker-Planck equation has a positive-semidefinite diffusion matrix. In this work, we analytically derive the SDEs for arbitrary Hamiltonian and jump operators based on the path-integral formula, independently of the derivation of the Fokker-Planck equation (FPE). In the course of the derivation, we formulate the path-integral representation of the GKSL equation by using the $s$-ordered quasiprobability, which systematically describes the aforementioned quasiprobabilities by changing the real parameter $s$. The essential point of this derivation is that we employ the Hubbard-Stratonovich (HS) transformation in the path integral, and its application is not always feasible. We find that the feasible condition of the HS transformation is identical to the positive-semidefiniteness condition of the diffusion matrix in the FPE. In the benchmark calculations, we confirm that the MC simulations of the obtained SDEs well reproduce the exact dynamics of physical quantities and non-equal time correlation functions of numerically solvable models, including the Bose-Hubbard model. This work clarifies the applicability of the approximation and gives systematic and simplified procedures to obtain the SDEs to be numerically solved.

Figures (5)

• Figure 1: Schematic images of (a) the path-integral representation \ref{['eq:path-integral representaiton continuous']} and (b)-(c) the approximations of the GKSL equation in the phase space. (b) Within the first order of quantum fluctuations, the GKSL equation is approximated into the generalized Liouville equation (GLE) GerlichSteeb, where each point distributed by the initial $\vec{s}$-ordered quasiprobability distribution function $W_{\vec{s}}(\vec{\alpha}_0,\vec{\alpha}^*_0,t_0)$ follows the classical equation of motion (CEOM), the equation of motion of the classical path. (c) The effects of the second order of the quantum fluctuations are incorporated into the classical path as Gaussian noises, where each point follows the stochastic differential equation (SDE). Here, the GKSL equation is approximated into the Fokker-Planck equation (FPE).
• Figure 2: Relaxation dynamics of two-site non-interacting atoms obeying the GKSL equation \ref{['eq:GKSL equation Model 1']} starting from the pure coherent state $\hat{\rho}(0) = \ket{\alpha_{{\rm I}1},\alpha_{{\rm I}2}}\bra{\alpha_{{\rm I}1},\alpha_{{\rm I}2}}$, where $\alpha_{{\rm I}1} = \sqrt{N_{{\rm I}1}}e^{i\pi/8}$ and $\alpha_{{\rm I}2} = \sqrt{N_{{\rm I}2}}e^{i\pi/4}$ with $N_{{\rm I}1} = 8$ and $N_{{\rm I}2} = 2$. Shown are (a) the difference of the remaining fractions of atoms at each sites $n_{12}$, (b) the correlation of atoms at different sites $C_{12}$, and real (c) and imaginary (d) parts of the non-equal time correlation function $G_{12}(t,0)$, which are defined by Eq. \ref{['eq:physical quantities']}. In each panel, we compare the numerically exact result (Exact) obtained by directly solving the GKSL equation \ref{['eq:GKSL equation Model 1']} and the ones of the second-order approximation using the Glauber-Sudarshan P (P:$2$nd), Wigner (W:$2$nd), and Husimi Q function (Q:$2$nd). We choose $\mu/(\hbar\gamma) = 1$ and $J/(\hbar\gamma) = 1$, and take 1000 samples for the initial conditions and 100 samples for the stochastic processes in the second-order approximation. The insets depict the difference between the results of each approximation and the numerically exact one. The shaded regions represent the standard error with respect to the sampling over the initial states.
• Figure 3: Relaxation dynamics of two-site Bose-Hubbard model obeying the GKSL equation \ref{['eq:GKSL equation Model2']} starting from the pure coherent state $\hat{\rho}(0) = \ket{\alpha_{{\rm I}1},\alpha_{{\rm I}2}}\bra{\alpha_{{\rm I}1},\alpha_{{\rm I}2}}$, where $\alpha_{{\rm I}1} = \sqrt{N_{{\rm I}1}}e^{i\pi/8}$ and $\alpha_{{\rm I}2} = \sqrt{N_{{\rm I}2}}e^{i\pi/4}$ with $N_{{\rm I}1} = N_{{\rm I}2} = 5$. The quantities and notations are the same as those in Fig. \ref{['fig:Model1']} except that we calculate the non-equal time correlation function $\bar{G}_{12}(t,0)$ instead of $G_{12}(t,0)$, and use the first-order approximation for the Glauber-Sudarshan P function. The other parameters are $\mu/(\hbar N_{\rm I}\gamma) = 1$, $J/(\hbar N_{\rm I}\gamma) = 1$, $U_{11}/(\hbar\gamma_3) = U_{22}/(\hbar\gamma_3) = 0.2$, and $\gamma_1/(N_{\rm I}\gamma_3) = \gamma_2/(N_{\rm I}\gamma_3) = 0.6$, where $N_{\rm I} = N_{{\rm I}1} + N_{{\rm I}2} = 10$, and we take 1000 samples for the initial conditions and 100 samples for the stochastic processes (W:$2$nd and Q:$2$nd). The insets depict the difference between the results of the approximations (P:$1$st, W:$2$nd, and Q:$2$nd) and the numerically exact one (Exact).
• Figure 4: Relaxation dynamics of two-component BEC obeying the GKSL equation \ref{['eq:GKSL equation Model3']} starting from the pure coherent state $\hat{\rho}(0) = \ket{\alpha_{{\rm I}1},\alpha_{{\rm I}2}}\bra{\alpha_{{\rm I}1},\alpha_{{\rm I}2}}$, where $\alpha_{{\rm I}1} = \sqrt{N_{{\rm I}1}}e^{i\pi/8}$ and $\alpha_{{\rm I}2} = \sqrt{N_{{\rm I}2}}e^{i\pi/4}$ with $N_{{\rm I}1} = 8$ and $N_{{\rm I}2} = 2$ at $\mu/(\hbar N_{\rm I}\gamma_3) = 1$, $U_{11}/(\hbar\gamma_3) = U_{22}/(\hbar\gamma_3) = U_{12}/(\hbar\gamma_3) = 1$, $\gamma_1/(N_{\rm I}\gamma_3) = 0.1$, and $\gamma_2/(N_{\rm I}\gamma_3) = 1$, where $N_{\rm I} = N_{{\rm I}1} + N_{{\rm I}2} = 10$. The quantities and notations are the same as those in Fig. \ref{['fig:Model2']} except that we use the first-order approximation using the Husimi Q function. The insets depict the difference between the results of the second-order approximation using the Wigner function (W:$2$nd) and the numerically exact one (Exact).
• Figure 5: Relaxation dynamics of two-site Bose-Hubbard model obeying the GKSL equation \ref{['eq:GKSL equation Model4']} starting from the pure coherent state $\hat{\rho}(0) = \ket{\alpha_{{\rm I}1},\alpha_{{\rm I}2}}\bra{\alpha_{{\rm I}1},\alpha_{{\rm I}2}}$, where $\alpha_{{\rm I}1} = \sqrt{N_{{\rm I}1}}e^{i\pi/8}$ and $\alpha_{{\rm I}2} = \sqrt{N_{{\rm I}2}}e^{i\pi/4}$ with $N_{{\rm I}1} = 8$ and $N_{{\rm I}2} = 2$ The other parameters are $\mu/(\hbar N_{\rm I}\gamma_3) = 1$, $J/(\hbar N_{\rm I}\gamma_3) = 0.4$, $U_{11}/(\hbar\gamma_3) = 1$, $U_{22}/(\hbar\gamma_3) = 0.2$, $\gamma_1/(N_{\rm I}\gamma_3) = 0.25$, and $\gamma_2/(N_{\rm I}\gamma_3) = 1.0$, where $N_{\rm I} = N_{{\rm I}1} + N_{{\rm I}2} = 10$. The notations in this figure are the same as those in Figs. \ref{['fig:Model3']} except that we use the second-order approximation for a hybrid of the Wigner and Husimi Q functions. The insets depict the difference between the results of the second-order approximations (W:$2$nd and W+Q:$2$nd) and the numerically exact one (Exact).