Weakly interacting Bose gas with two-body losses

We study the many-body dynamics of weakly interacting Bose gases with two-particle losses. We show that both the two-body interactions and losses in atomic gases may be tuned by controlling the inelastic scattering process between atoms by an optical Feshbach resonance. Interestingly, the low-energy behavior of the scattering amplitude is governed by a single parameter, i.e. the complex $s$-wave scattering length $a_c$. The many-body dynamics are thus described by a Lindblad master equation with complex scattering length. We solve this equation by applying the Bogoliubov approximation in analogy to the closed systems. Various peculiar dynamical properties are discovered, some of them may be regarded as the dissipative counterparts of the celebrated results in closed Bose gases. For example, we show that the next-order correction to the mean-field particle decay rate is to the order of $|n a_c^{3}|^{1/2}$, which is an analogy of the Lee-Huang-Yang correction of Bose gases. It is also found that there exists a dynamical symmetry of symplectic group Sp$(4,\mathbb{C})$ in the quadratic Bogoliubov master equation, which is an analogy of the SU(1,1) dynamical symmetry of the corresponding closed system. We further confirmed the validity of the Bogoliubov approximation by comparing its results with a full numerical calculation in a double-well toy model. Generalizations of other alternative approaches such as the dissipative version of the Gross-Pitaevskii equation and hydrodynamic theory are also discussed in the last.


Introduction
Interacting Bose gas is a central topic in the fields of ultracold atoms and condensed matter physics, for it represents a paradigm of quantum many-body physics.Theoretically, if one focuses on the low-energy physics, the interactions between bosons can always be simplified by a zero-range one with a real s-wave scattering length a s that reproduces the same low-energy scattering phase shift [1], despite that the actual potentials are usually very complicated.Various many-body effects can then be studied theoretically with the help of the zero-range model [2,3].More importantly, experimental techniques such as Feshbach resonance can adjust a along the real-axis in cold atomic gases [4], which allows us to demonstrate various many-body effects, from the Lee-Huang-Yang correction [5][6][7][8][9] for positive scattering length to the Bose nova effect [10][11][12] for negative scattering length.
Recently, with the development of theoretical and experimental methods [13][14][15][16][17][18][19][20][21][22][23][24][25][26][27][28][29], much more attention has been paid to behaviors of cold atom systems with dissipation.For bosonic and fermionic models, single particle and two-body dissipation processes such as particle pump and loss can be realized [30][31][32][33][34][35][36].To better understand the many-body physics in open systems, it is then useful to introduce a zero-range model for Bose gases with two-body losses.In a previous work [37], the authors have shown that with a proper renormalization or regularization approach, the boson-boson interactions and losses can be effectively described by a complex contact interaction parameterized by a complex s-wave scattering length a c .The physics of the dissipative Bose gases can then be regarded as an extension of a s from the real axis to the complex plane.
In this paper, we focus on the dissipative dynamics of weakly interacting Bose gas with two-body losses.Firstly, we show that complex scattering length a c can be effectively adjusted in the lower complex plane through optical Feshbach resonances.Then we study the manybody dynamics in the weakly interacting and dissipating region by introducing the Bogoliubov approximation [38].Within the Bogoliubov approximation, the evolution process is governed by a quadratic time-dependent Lindblad equation.By numerically solving a toy model, we verify the accuracy of this approximation in open systems.Furthermore, we find within this approximation, the superoperators in the master equation form a closed algebra, which corresponds to a symplectic Sp(4, ) dynamical symmetry group, and helps to derive an exact solution for the evolution of density matrix.Our analysis also reveals that an n-mode bosonic system governed by a quadratic Lindbladian always has Sp(2n, ) dynamical symmetry.Finally, within the framework of the Keldysh path-integral method, we obtain a generalized non-Hermitian Gross-Pitaevskii equation which reproduces the same excitation spectrum as the Lindblad equation and leads to a set of dissipative hydrodynamic equations in the long wavelength limit.
The paper is organized as follows.In Section.2, we show that the complex scattering length a c in an atomic gas can be tuned across the lower half complex plane via optical Feshbach resonances.Then in Section. 3 we introduce the single-channel master equation and the Bogoliubov approximation, which we use to solve the many-body properties of the dissipative Bose gases.We verify the validity of the Bogoliubov approximation by numerically solving a toy model in Section.4, and discuss the dynamical symmetry of the Bogoliubov Lindbladian in Section. 5.In Section.6, we derive the dissipative version of the Gross-Pitaevskii equation and hydrodynamic theory using the Keldysh path integral formalism for open systems.

Tuning complex scattering lengths in experiments
In the scattering theory, it is known that the inelastic collisions between particles will cause two-body losses and eventually lead to a complex s-wave scattering length, which contains all the information of the low-energy scattering amplitude [1,[39][40][41].This suggests that one might control the two-particle interaction and losses experimentally by tuning the inelastic scattering process between atoms in a cold atomic gas.Experimentally, this tuning can be realized through the optical Feshbach resonance technique [42].As shown in Fig. 1 (a), in the experiment, an external laser is applied to the atomic gas, and the frequency of the optical field is tuned close to a transition between two ground state atoms and a highly excited molecular state (i.e. a two-body bound state consists of a ground state atom and a highly excited atom).When the excited molecule spontaneously decays and emits a photon, the two atoms will be kicked out of the system because of the huge recoil momentum, thus causing two-particle losses in the atomic gas.
Theoretically, this process might be captured by a detailed calculation of the scattering amplitude for a finite-range multi-channel model, which shows that the complex scattering length is given by [39][40][41] Here ν is the detuning of the laser field, γ is the decay rate of the excited molecule, a bg is the (real) background s-wave scattering length, E is the scattering energy, and Γ (I) stands for the width of the resonance and is proportional to the intensity I of the laser field.
In this chapter, we provide a simplified zero-range model to reproduce this formula and discuss the domains on the complex plane where a c might reach using optical Feshbach resonances.To construct the model, we first introduce the bosonic field operator b ( d) for the ground state atoms (excited molecules), and the Hamiltonian of the system may be written as (ħ h is set to unity in this paper) , where g is the interacting strength between atoms in the open channel, ε k = k 2 2m and ξ k = k 2 4m +ν, m is the atom mass, Ω is the volume.The two channels are coupled by laser light with strength α.
The spontaneous radiation of the bound-state molecule can be characterized by the following Lindblad master equation [43] For a two-body process, the dynamics of the corresponding density matrix (in two-particle Fock space) is equivalent to the evolution under a non-hermitian Hamiltonian Ĥ2-body = Ĥofr − i γ We can obtain the two-body scattering matrix T 2 for this Ĥeff at relative kinetic energy By matching the derived T-matrix with the usual low-energy expansion formula 4π m 1 a c (E) + i mE −1 [1], we obtain the renormalize relation for the complex scattering length a c , where Λ is the momentum cut-off and a c (E) can be further written as 2π 2 +mgΛ .The above formula is exactly eq.(1) if we take the limit Λ → ∞, and write Γ (I) = m|α re |2 4πa bg . 1or a dilute gas at extremely low temperature, the kinetic energy E is negligible and the many-body effect can be well captured by the zero energy scattering length a c (E = 0).A general a c (0) can be achieved via changing the detuning ν and density I. To be more clear, as shown in Fig. 1(b), at fixed Γ (I), the trajectory of a −1 c (E = 0) form a circle tangent to the realaxis at a −1 bg with radius γa bg in the upper half complex when tuning ν ∈ (−∞, +∞).Therefore, as we gradually increase Γ (I), this circle can sweep across the entire upper half-plane.
Here we comment on the difference between tuning complex scattering length and real scattering length.The real scattering length can also be controlled in a magnetic-induced Feshbach resonance by varying the strength of the magnetic field.However, to achieve a general complex scattering length, we need at least two adjustable parameters.Thus we choose the optical resonance model rather than the magnetic Feshbach resonance.

Dissipative Bose gas
For non-dissipative atomic gases, it is known that the interaction between particles may be described by a single-channel contact potential, provided that the Feshbach resonance is "wide", i.e. that the occupation in the closed molecular channel is small compared to the open scattering channel [41].The Hermitian Hamiltonian of a Bose gas may be written as where âk is the annihilation operator for the atoms with momentum k, ε k is the kinetic energy, g is the coupling constant, Ω is the system volume.
The single-channel Hamiltonian greatly simplifies the calculation of the many-body properties in closed Bose gases.Thus it is desirable to construct a single-channel model to describe the many-body dynamics of dissipative atomic gases with complex scattering lengths.In the previous work [37], the authors have studied this problem and constructed a renormalizable single-channel model for systems across a "wide" optical Feshbach resonance.In this work, we focus on the many-body dynamics of this model, and the study of systems across "narrow" optical Feshbach resonances will be pursued in the future.In the single-channel model, the open system dynamics are governed by the master equation [37,43] with the non-Hermitian effective Hamiltonian and the recycling term This master equation describes the dynamics of Bose atoms, denoted by â, as they interact with each other through a bare coupling constant g b while decaying to the environment (i.e.no longer confined by the trapping potential) via a two-body losses process characterized by a bare coupling constant γ b .To regularize the contact interaction, we can use the renormalization relation and define g − iγ ≡ 4πa c m as renormalized complex coupling constant.

Bogoliubov approximation
To solve the spectrum of isolated weakly interacting Bose gas, we can use the Bogoliubov approximation and diagonalize the quadratic Bogoliubov Hamiltonian.Similarly, we can generalize Bogoliubov's transformation in the open system.Starting from a condensate initial state and assuming the condensate part is always much larger than quantum depletion during the time evolution, we can trace out the condensed part by replacing all the â0 , â † 0 with the square root of the zero-momentum atom number N − k̸ =0 â † k âk , then we obtain an approximate master equation for the reduced density matrix ρ′ as and where N is the total atom number and n = N Ω is the density.Here we use renormalized parameters g, γ to replace the bare parameters g b , γ b like in conventional Bogoliubov Hamiltonian [38].Different from the time-independent Bogoliubov Hamiltonian in a closed system, Lindbladian after Bogoliubov transformation is time-dependent due to the two-body losses.The density n is decreasing as a function of time t, at mean-field level, we have The conventional Bogoliubov Hamiltonian is the linear combination of three operators As a result, for an isolated system, the Hamiltonian has SU(1,1) dynamical symmetry, and the Heisenberg equations of these three operators are closed.When introducing the two-body losses, the Heisenberg equations for these three operators are no longer closed.However, the evolution of their expectation values Here, we focused on the short-time dynamics where γnt ≪ 1 such that we may approximate n as a constant at a fixed time t.We can diagonalize the 3×3 matrix in eq. ( 12) to obtain the eigenvalue −2iξ k,i , where the quasi-steady state eigenvalue is A phase boundary between stable and unstable BEC can be obtained by this quasi-steady state eigenvalue [37].

Lee-Huang-Yang correction
Furthermore, we will give a detailed derivation for the Lee-Huang-Yang correction [47,48] for this open system.Consider quasi-steady-state eigenvalue, up to O (γnt), simple solution of eq. ( 12) can be obtained constant vectors B, C, D can be determined by the initial value of A, and We see that A never reaches the steady value A s because the asymptotic solution eq. ( 13) only works for time interval 0 ≤ t ≪ ħ h/γn.Nevertheless, it is still useful to calculate some physical properties for this steady state, since this helps to verify the renormalization relation given in eq. ( 9).When A reaches this steady value, we see that It is then straightforward to show that This leads to From this equation, we see that the first term corresponds to the mean-field decay which leads to n ≃ n(1 + 2γnt) −1 as mentioned previously.After inserting the steady value, we have We immediately find that the summation diverges for large momentum.Similar divergence occurs in the calculation for ground state energy of Bose gas in a closed system [49,50].This divergence arises from the fact that the renormalized interaction is only valid for small momenta.To cure this divergence, one can introduce an intermediate momentum cutoff and then consider an effective interaction with second-order processes in the mean-field energy component.For our case, the second-order effective interaction parameter g and γ is given by We thus obtain Replacing γ in the mean-field part of eq. ( 17) by γ, we have with c θ a constant that only depends on the argument of the complex scattering length a c , Here, θ is defined as where a c is expressed as two real parameters: a c = a r + ia i .We note that the n|a c | 3 term in eq. ( 20) is an analog of the celebrated Lee-Huang-Yang correction to the ground state energy of a weakly interacting Bose gas [47,48], With the help of renormalization relation, we can eliminate the divergence eq. ( 16) and obtain a physical result of particle loss rate.
It is well known that, for Bose gas without a two-particle loss (a i = 0), the system is dynamically unstable when a r < 0. This instability is reflected in eq. ( 21), which becomes ill-defined for a negative real scattering length.Similarly, eq. ( 20) reflects certain dynamic instability in open systems.One can check that the coefficient c θ becomes ill-defined for θ ∈ [5π/6, π], suggesting that our approach is not valid in this regime.This is because the quantum depletion (Bogoliubov modes) grows rapidly in this regime, which invalidates the assumption that â0 ≈ â † 0 ≈ N 0 [37].

A toy model
Our calculation is based on the assumption that the Bogoliubov approximation is always valid during the dynamical evolution.It is impossible to verify this approximation numerically in the thermodynamic limit due to the exponential growth of the Hilbert space dimension in the quantum many-body system.To numerically compare the dynamics process generated by the original Lindbladian and the Bogoliubov Lindbladian, we introduce a toy model that can capture the interacting and dissipation features of the open bosonic system.

Numerical model
As shown in Fig. 2, we consider a double-well model, the annihilation operator for the bosonic mode in the right (left) well is denoted by br ( bl ).The evolution of this system is governed by the master equation where the Hamiltonian is , quasi-steady state eigenvalue has a positive imaginary part, which represents the exponential growth of non-condensate particle number, the system is in the unstable phase.When , the number of the minority particle is always much smaller than the condensate number during the time evolution, the system is in the stable phase.
with t the hopping strength between the two wells, g the on-site interaction strength.The Lindblad operators for the onsite two-body losses are Γr = br br , Γl = bl bl , with γ the decay rate.
Corresponding to the Bose gas model, we consider the "zero momentum" mode â0 = 1 2 ( br + bl ) as majority part while the "non-zero momentum" mode â1 = 1 2 ( br − bl ) as small depletion part.We can shift the energy for the mode â0 to zero without loss of generality and the mode â1 has a detuning ε = 2t from the mode â0 .Omitting the interacting or losses terms only including minority part, such as â † 1 â † 1 â1 â1 , the master equation is then given by where the effective Hamiltonian Ĥeff a for â0 , â1 is and the recycling term in eq. ( 24) are similar to Lindbladian of Bose gas with two-body losses, In this model, the effective Hamiltonian Ĥeff a conserves the particle number n = â † 0 â0 + â † 1 â1 while the recycling term always annihilates two particles, thus the master equation can be decomposed to a series of hierarchy equations for ρi,j [37] where ρi,j = Pi ρ Pj , with Pi the projection operator for the i-particle Fock space.Then for the dynamical evolution of eq. ( 24) with initial density matrix a pure state of particle number N , only the projected density matrix ρN,N , ρN−2,N−2 , ρN−4,N−4 ... is involved and the numerical simulation is performed by calculating these differential equations.Now considering the Bogoliubov approximation, we replace the â0 and â † 0 with the square root of the condensate atoms number n − â † 1 â1 , then we can obtain the approximated Lindbladian as  24); dashed: approximated Lindbladian eq. ( 27)).The initial condition is all particles condensate at groundstate, n 0 = N .(a),(b)System at stable phase with coupling strength g N = 0.52ε and dissipation strength γN = 0.20ε.(c),(d) System at stable phase with g N = −0.52εand γN = 0.02ε.
where the Hamiltonian ĤB dw is and time-dependent particle number is given by n = N (1 + 2γN t) −1 .Similar to eq. ( 12) in the Bose gas case, the expectation value of the three SU(1,1) operators with A dw = (〈 Â1 dw 〉 , 〈 Â2 dw 〉 , 〈 Â3 dw 〉) T .The three eigenvalues −2iξ i for the 3×3 matrix in eq. ( 29) are When the imaginary part of ξ 1 is larger than zero, the corresponding eigenvalue becomes a positive number, thus the particle number n1 grows exponentially in short time leading to an unstable dynamic.The phase boundary2 for this stability is given by as shown in Fig. 2(b).

Numerical results
We now check the accuracy of the Bogoliubov approximation by comparing the numerical simulation result from eq. ( 24) and eq.( 27).The evolution begins with an initial state with N = 100 particles occupying the mode â0 .As shown in Fig. 3, the two evolution results almost coincide in the stable phase.While in the unstable phase, due to the increase of the quantum depletion, the particle number n1 predicted by the Bogoliubov approximation eq. ( 27) deviates from the exact time evolution governed by eq. ( 24) at long time.This result confirms that our Bogoliubov approximation for the Lindblad master equation is valid as long as the ratio between the density of quantum depletion to the total density is much smaller than the unitary.

Dynamical symmetry of Lindbladian
In this section, we discuss the dynamical symmetry of quadratic Lindbladians in detail, which can greatly simplify the calculation of the master eq. ( 10).It is worth mentioning that a time-independent quadratic Lindbladian may be solved explicitly using the third quantization method [51,52].However, this cannot be directly applied to our problem because of the timedependent n(t) contained in the master eq. ( 10).Fortunately, in our Bogoliubov Lindbladian, it can be shown that, besides the obvious U(1) and Z 2 symmetry, the algebraic structures of the superoperators defined on the density matrix space also contain a hidden symplectic Sp(4, ) dynamical symmetry.This dynamical symmetry provides a relatively simple way [53,54] to construct exact solutions to the time-dependent Lindbladian algebraically.

Symmetry
The Lindbladian L B in eq. ( 10) has two obvious symmetries, i.e. an extended U(1) symmetry and a Z 2 symmetry.To see the U(1) symmetry, one only needs to verify that the Lindbladian is invariant under transformation âk → âk e iφ , â † k → â † k e −iφ .The corresponding conserved superoperator is thus Q = [ N , •].While for the Z 2 symmetry, it can be verified by realizing L B is invariant under âk → −â k or â † k → −â † k .However, these two symmetries alone are not enough to obtain an analytical solution to the master eq. ( 10).To solve the quench dynamics governed by this time-dependent master equation, we need to generalize the concept of dynamical symmetry to open systems, i.e. to find a closed algebra formed by superoperators that contain the Lindbladian L B itself.

Closed algebra
To find this algebra, we note that L B can be expressed in the superoperator-formalism as Here Lk is a superoperator act only on modes k and −k, which is a linear combination of seven quadratic superoperators.We label these superoperators by hk i , i = 1, 2, . . .7. They are defined by Among these seven superoperators, hi k , i = 1, 2, . . .6 represent the (anti-)commutators between the (anti-)Hermitian parts of the effective Hamiltonian and the density matrix ρ, and the last one represents the recycling term of Lk .
However, only these seven superoperators cannot form a closed algebra, i.e. their commutators are not necessarily linear combinations of themselves.For example, we have which is a superoperator that can not be written in the form of It is thus necessary to introduce three more superoperators hi k , i = 8, 9, 10 to close the algebra.They are All of the superoperators hi k , i = 1, 2, . . ., 10 now form a closed algebra, 3 whose commutation relations are listed in Table .1. Hence these superoperators are generators of the dynamical symmetry group of Lindbladian.In the following of this section, because only particles at momentum k and −k are coupled, we will omit the label k of superoperators.Below we will prove this algebra can map to the algebra of two coupled harmonic oscillators and the corresponding group is isomorphic to Sp(4, ).Before that, we will formally write the exact solution of LB .

Exact solution
Based on this closed algebra structure, the dynamical problem of time-dependent Lindbladian can be solved exactly [53,54].First of all, formally solve the master equation ρ as Λ(t, t 0 ) is a dynamical map from time t 0 to time t.In general, this mapping is a semi-group which satisfies Λ(t 2 , t 0 ) = Λ(t 2 , t 1 ) Λ(t 1 , t 0 ) but don't satisfy the unitary condition under Lindblad time evolution.Thanks for the closed algebra, the Lindbladian can be written as a linear combination of the dynamical symmetry group generators L = 7 l=1 u l (t) hl hence we can parametrize the dynamical map by these ten generators, 3 In fact, the structure of the master equation can further restrict these 10 elements into seven: H0 Substitute this equation back to the master equation, we can get where u l (t) are time-dependent parameters defined by Hamiltonian.Solving the dynamics of the time-evolution operator is equivalent to solving the complex functions g i (t).Right multiply the inverse of Λ(t, t 0 ) at both sides of eq. ( 35) and write the explicit expression of time derivative, we can obtain Using Baker-Campbell-Hausdorf formula [55] and commutation relation in Table .1, we can formally write L.H.S of eq. ( 36) as: where the η mn are analytic functions of g.Considering the linear independent property of ten generators, we can obtain a set of coupled first-order differential equations.In consequence, by solving these coupled differential equations for g(t), we will get the exact solution of the dynamics governed by time-dependent Lindbladian.In practice, obtaining an analytical solution for a general time dependence u l (t) is difficult.However, eq. ( 37) offers a method to numerically obtain the time evolution of the density matrix, which is similar to the evolution of a Gaussian state under a quadratic Hamiltonian [46,56].
In the following, we show that the superoperators hi form an algebra of sp(4, ) by mapping them to a coupled 2-mode harmonic oscillator.We further generalize this result to an n-mode quadratic Lindblad bosonic system, in which the superoperators form a closed algebra of sp(2n, ).

Map to harmonic oscillators
Even though the commutation relations between hi k seem complicated as shown in Table .1, it is worth noting that they only consist of quadratic superoperators.A natural question is then whether there is an isolated system that has the same algebra structure.Indeed, we find that if we consider two coupled harmonic oscillators with ladder operators a and b, all the superoperators hi can be mapped to linear combinations of quadratic forms of a, a † , b, b † .In Table .2, we give the explicit correspondence of this mapping, and it is not difficult to show that the mapping preserves commutation relations, i.e. it is indeed an isomorphism.The isomorphism greatly the structure of the algebra.As an example, we will show in the following that the quadratics of ladder operators â and b form the algebra of symplectic group Sp(4, ).

Sp(2n, ) dynamical symmetry group
To show this, we prove a general result, i.e. the quadratic forms of ladder operators a n in an n-mode harmonic oscillators (n = 1, 2, . ..) have symplectic Sp(2n, ) dynamical symmetry.A Table 2: Operators in coupled harmonic oscillator system and superoperators in dissipative Bose gas system after Bogoliubov approximation with momentum k.These 10 operators and superoperators have the same commutation relation in Table .1 and they can form the algebra sp(4, ).

Harmonic Oscillator Bogoliubov Lindbladian
To make the coefficient matrices unique, we require B T = B and C T = C.Note that the coefficient matrix may be written as with M satisfying which is the generator for Lie group Sp(2n, ).
We now define new operators which are related by b i = Ω i j b j and b i = Ω i j b j with matrices We thus have the following useful commutation relations With all these definitions, the generic quadratic form in Eq. ( 38) can then be expressed as with the upper (lower) index representing the column (row) index.The above formula gives a one-to-one mapping between the quadratic form and the Lie algebra sp(2n, ).More importantly, it can be proved (Appendix A) which implies that the mapping is an isomorphism.Now we prove n-mode quadratic bosonic Hamiltonian has an Sp(2n, ) dynamical symmetry, meanwhile, in the last subsection we proved our Bogoliubov Lindbladian is isomorphic to coupled 2-mode harmonic oscillator.In consequence, we can conclude that Bogoliubov Lindbladian has an Sp(4, ) dynamical symmetry.
Furthermore, our results are not only restricted to Bogoliubov Lindbladian.Because quadratic operators in 2n-mode coupled harmonic oscillator have the same commutation relation with the superoperators in n-mode quadratic Lindbladian, we can conclude that quadratic Lindbladian which is constituted by n-mode bosons has Sp(2n, ) dynamical symmetry.This result will be useful for the further analytical study of open systems.

Dissipative Gross-Pitaevskii equation and hydrodynamic theory
In closed systems, weakly interacting Bose gas can also be treated using other theoretical approaches such as the Gross-Pitaevskii equation and the hydrodynamic theory [49,50,57].In this section, we generalize these two descriptions to open systems with weak two-body losses.We derive the dissipative versions of the Gross-Pitaevskii equation and hydrodynamic equations based on the Keldysh path integral formalism.

Keldysh formalism
We start from the master equation ∂ t ρ = L ρ of interacting bosons subject to two-body losses in real space, where ψ(r) is the bosonic annihilation operator at position r, and is Hermitian Hamiltonian of interacting bosons.The recycling term J ρ is given by Using the keldysh path-integral representation, we introduce fields ψ ± and ψ± living respectively on the time contour C + and C − , where time runs from −∞ to +∞ and then back to −∞ in the closed-contour C + ∪ C − .Then we can write partition function Z = Tr(ρ(t)) of eq. ( 48), where we omit the initial condition and the Lindblad-Keldysh action [58] is where time runs from −∞ to +∞, and s η = 0 for η = +, s η = 1 for η = −.

Saddle-point approximation
We consider the case that almost all particles occupy the ground state energy level, which means N ≈ N 0 ≫ 1.In the closed system, we can assume the system is always in a coherent state, and the condensate wavefunction can be found by the saddle-point equation.Similarly, here we consider the losses process to be weak and slow so that we can still take the coherent state assumption.As a result, we can take the saddle-point equation in an open system, then we obtain where the complex interacting strength g c = g − iγ.

Non-Hermitian Gross-Pitaevskii equation
The regularization of Keldysh action requires the relation ψ + = ψ − , ψ+ = ψ− [59], combining with the saddle-point equations eq. ( 54), we will obtain the Gross-Pitaevskii equation under the two-body losses (ψ + = ψ), This equation substitutes the coupling parameter g to the complex version g − iγ in the conventional Gross-Pitaevskii equation.However, due to the particle loss, the solution of this dissipative Gross-Pitaevskii equation becomes quite different from the conventional solution for a closed system.For example, when V = 0, it is easy to show that we have a plane-wave-like solution given by ψ 0 (r) = n 0 e −i g c 2γ log(1+2γn 0 t) .
And in finite momentum we can also obtain an exact solution, ψ p (r) = n 0 e i(p•r−p 2 t/2) e −i g c 2γ log(1+2γn 0 t) .
Then we can solve the excitation spectrum of this dissipative Gross-Pitaevskii equation by inserting small perturbation ψ = ψ 0 + δψ and its complex conjugate and expanding to the linear order in δψ: To solve this equation, we may apply a unitary transformation δφ = e i g 2γ log(1+2γn 0 t) δψ which leads to a simplified equation in the "rotating wave frame", with n = n 0 (1 + 2γn 0 t) −1 .Now if we consider perturbations with momentum k and naively treat n(t) as time-invariant, we may find the eigenfrequencies of the above equations are which coincide with the values we previously obtained using the Bogoliubov approximation.
As a result, we can also get the phase diagram by solving this dissipative Gross-Pitaevskii equation.

Hydrodynamic equation
In the closed system, fruitful physical consequences can be obtained by solving the hydrodynamic theory of interacting BEC, such as superfluidity, anisotropic expansion, and low-energy modes in a harmonic trap [49,50,57].It is then interesting to construct the hydrodynamic equation with two-body loss.Starting from the time-dependent dissipative Gross-Pitaevskii equation eq. ( 55), we can derive the hydrodynamic equations which govern the dynamics of a dissipative fluid.By decomposing ψ = ρe iθ , the real part and imaginary part give the hydrodynamic equations: The first equation is called the continuity equation, which reflects the conservation of particle numbers.Now two-body losses in interacting BEC bring new term γρ 2 in the continuity equation, which breaks the particle number conservation.And the second equation, the Newton equation, is same as the conventional one.For a uniform system with V (r) = 0, the uniform solution is ρ(r, t) = ρ 0 1 + 2γtρ 0 , (59) which shows that condensate particles always decay with time and the vacuum is the only true steady state of this system.With the new dissipation term γρ 2 , finding the solutions for the general case is challenging, we will leave this for future research.

Conclusion
In summary, we systematically study the many-body dynamics of weakly interacting Bose gases with two-body losses.It is shown that both the two-body interactions and losses in a cold atomic gas may be described by a complex scattering length a c , which may be controlled via tuning an external laser field.We generalize Bogoliubov approximation to open systems and verify the validity of this approximation by numerical simulating a toy model that has a similar structure.Based on this time-dependent Bogoliubov Lindbladian, we study the quench problem and prove this system has a Sp(4, ) dynamical symmetry, which is crucial for the exact calculation of quench dynamics.Furthermore, we show a general n-mode quadratic Lindbladian of the bosonic system has a dynamical symmetry of Sp(2n, ), which is useful for the analytical understanding of the dissipative open system.On the other hand, we also generalize the Gross-Pitaevskii equation and hydrodynamic theory to dissipative Bose gases.Hydrodynamic equations have rich interesting solutions in closed systems, we only discuss the solution without external potential in this paper, it would be significant to generalize these solutions to open systems and study the stable property of these results.Finally, from a closed interacting bosonic model to open system, other types dissipation are admitted such as single body loss and pump, particle number form dissipation.It is also interesting to discuss the different properties of these open systems.We hope our work can be helpful for further understanding the non-equilibrium dissipative dynamics in cold atom systems.

Figure 1 :
Figure 1: (a)The schematic for the optical Feshbach resonance model.In the open channel, a pair of atoms b is interacting via strength g.They are coupled to a molecule state d in the closed channel by laser light at detuning ν.The molecular state spontaneously radiates at rate γ.(b) The range of a −1 c (0).For fixed Γ (I), the trajectory of a −1 c (0) forms a circle tangent to the real-axis at a −1 bg with radius Γ (I) γa bg in the upper half complex plane when changing ν within (−∞, +∞).Then the radius of this circle will increase as we gradually turn on Γ (I).