# Weakly interacting Bose gas with two-body losses

Chang Liu  Institute for Advanced Study, Tsinghua University, Beijing 100084,
China  Zhe-Yu Shi  State Key Laboratory of Precision Spectroscopy, East China
Normal University, Shanghai 200062, China  Ce Wang<EMAIL_ADDRESS>School
of Physics Science and Engineering, Tongji University, Shanghai 200092, China

###### Abstract

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
$|na_{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 the
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.

##  I 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_{\text{s}}$ that reproduces the same low-
energy scattering phase shift Landau and Lifshitz (2013), despite the fact
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 Blatt and Weisskopf (1991); Huang and Yang (1957). More importantly,
experimental techniques such as Feshbach resonance can adjust $a$ along the
real-axis in cold atomic gases Inouye  _et al._ (1998), which allows us to
demonstrate various many-body effects, from the Lee-Huang-Yang correction
Altmeyer  _et al._ (2007); Shin  _et al._ (2008); Papp  _et al._ (2008); Navon
_et al._ (2010, 2011) for positive scattering length to the Bose nova effect
Donley  _et al._ (2001); Cornish  _et al._ (2006); Strecker  _et al._ (2002)
for negative scattering length.

Recently, with the development of theoretical and experimental methods Ritsch
_et al._ (2013); Kasprzak  _et al._ (2006); Mabuchi and Doherty (2002); Öztop
_et al._ (2012); Deng  _et al._ (2010); García-Ripoll  _et al._ (2009);
Yamamoto  _et al._ (2019); He  _et al._ (2020); Bouchoule  _et al._ (2020);
Yamamoto  _et al._ (2021); Bouchoule and Dubail (2021); Rossini  _et al._
(2021); Rosso  _et al._ (2022); Nakagawa  _et al._ (2021); Bouchoule  _et al._
(2021), 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
Gerbier and Castin (2010); Bouganne  _et al._ (2020); Tindall  _et al._
(2019); Kantian  _et al._ (2009); Syassen  _et al._ (2008); Tomita  _et al._
(2017); Sponselee  _et al._ (2018). 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 Wang  _et al._ (2021), 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_{\text{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 many-body dynamics in the weakly
interacting and dissipating region by introducing the Bogoliubov approximation
Bogoliubov (1947). 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 an symplectic
Sp$(4,\mathbb{C})$ 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,\mathbb{C})$ 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. II, 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. III 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. IV, and discuss the dynamical symmetry of the
Bogoliubov Lindbladian in Section. V. In Section. VI, we derive the
dissipative version of the Gross-Pitaevskii equation and hydrodynamic theory
using the Keldysh path integral formalism for open systems.

##  II Tuning complex scattering lengths in experiments

In the scattering theory, it is known that the two-boy losses will cause
inelastic collisions between particles and eventually lead to a complex
$s$-wave scattering length, which contains all the information of the low-
energy scattering amplitude Landau and Lifshitz (2013); Bohn and Julienne
(1996, 1999); Chin  _et al._ (2010). 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
Thalhammer  _et al._ (2005). 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

$\displaystyle
a_{c}(E)=a_{\text{bg}}\left(1+\frac{\Gamma(I)}{E-\nu-\Gamma(I)+i\gamma/2}\right).$
(1)

Here $\nu$ is the detuning of the laser field, $\gamma$ is the decay rate of
the excited molecule, $a_{\text{bg}}$ is the (real) back ground $s$-wave
scattering length, $E$ is the scattering energy, and $\Gamma(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 $\hat{b}$ ($\hat{d}$) for the ground state atoms
(excited molecules), and the Hamiltonian of the system may be written as
($\hbar$ is set to unity in this paper)

$\displaystyle\hat{H}_{\text{ofr}}=$
$\displaystyle\sum_{\mathbf{k}}\left(\epsilon_{\mathbf{k}}\hat{b}_{\mathbf{k}}^{\dagger}\hat{b}_{\mathbf{k}}+\xi_{\mathbf{k}}\hat{d}_{\mathbf{k}}^{\dagger}\hat{d}_{\mathbf{k}}\right)+\frac{g}{2\Omega}\sum_{\mathbf{k},\mathbf{k^{\prime}},\mathbf{p}}\hat{b}^{\dagger}_{\mathbf{k}+\mathbf{p}}\hat{b}^{\dagger}_{\mathbf{k^{\prime}}-\mathbf{p}}\hat{b}_{\mathbf{k^{\prime}}}\hat{b}_{\mathbf{k}}$
$\displaystyle+$
$\displaystyle\frac{1}{\sqrt{\Omega}}\sum_{\mathbf{k},\mathbf{p}}\left(\alpha\hat{d}_{\mathbf{p}}^{\dagger}\hat{b}_{\frac{\mathbf{p}}{2}-\mathbf{k}}\hat{b}_{\frac{\mathbf{p}}{2}+\mathbf{k}}+\alpha^{*}\hat{d}_{\mathbf{p}}\hat{b}_{\frac{\mathbf{p}}{2}-\mathbf{k}}^{\dagger}\hat{b}_{\frac{\mathbf{p}}{2}+\mathbf{k}}^{\dagger}\right),$

where $g$ is the interacting strength between atoms in the open channel,
$\epsilon_{\mathbf{k}}=\frac{k^{2}}{2m}$ and
$\xi_{\mathbf{k}}=\frac{k^{2}}{4m}+\nu$, $m$ is the atom mass, $\Omega$ is the
volume. The two channels are coupled by laser light with strength $\alpha$.

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 $\nu$. The molecular state spontaneously radiates at rate $\gamma$.
(b) The range of $a_{c}^{-1}(0)$. For fixed $\Gamma(I)$, the trajectory of
$a_{c}^{-1}(0)$ forms a circle tangent to the real-axis at
$a_{\text{bg}}^{-1}$ with radius $\frac{\Gamma(I)}{\gamma a_{\text{bg}}}$ in
the upper half complex plane when changing $\nu$ within $(-\infty,+\infty)$.
Then the radius of this circle will increase as we gradually turn on
$\Gamma(I)$.

The spontaneous radiation of the bound-state molecule can be characterized by
the following Lindblad master equation Breuer  _et al._ (2002)

$\partial_{t}\hat{\rho}=\mathcal{L}(\hat{\rho})=-i[\hat{H}_{\text{ofr}},\hat{\rho}]-\frac{\gamma}{2}\sum_{\mathbf{k}}\\{\hat{d}^{\dagger}_{\mathbf{k}}\hat{d}_{\mathbf{k}},\hat{\rho}\\}+\gamma\sum_{\mathbf{k}}\hat{d}_{\mathbf{k}}\hat{\rho}\hat{d}_{\mathbf{k}}^{\dagger}.$

For 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

$\hat{H}_{\text{2-body}}=\hat{H}_{\text{ofr}}-i\frac{\gamma}{2}\sum_{\mathbf{k}}\hat{d}_{\mathbf{k}}^{\dagger}\hat{d}_{\mathbf{k}}.$
(2)

We can obtain the two-body scattering matrix $T_{2}$ for this
$\hat{H}_{\text{eff}}$ at relative kinetic energy $E$

$T_{2}(E)=\left(\left(g+\frac{|\alpha|^{2}}{E-\nu+i\gamma/2}\right)^{-1}-\frac{1}{\Omega}\sum_{k}\frac{1}{E-k^{2}/m}\right)^{-1}.$

By matching the derived T-matrix with the usual low-energy expansion formula
$\frac{4\pi}{m}\left(\frac{1}{a_{c}(E)}+i\sqrt{mE}\right)^{-1}$ Landau and
Lifshitz (2013), we obtain the renormalize relation for the complex scattering
length $a_{c}$,

$\frac{m}{4\pi
a_{c}(E)}=\left(g+\frac{|\alpha|^{2}}{E-\nu+i\gamma/2}\right)^{-1}+\frac{1}{\Omega}\sum_{|\mathbf{k}|<\Lambda}\frac{m}{k^{2}},$
(3)

where $\Lambda$ is the momentum cut-off and $a_{c}(E)$ can be further written
as

$a_{c}(E)=a_{\text{bg}}+\frac{m}{4\pi}\frac{|\alpha_{\text{re}}|^{2}}{E-\nu_{\text{re}}+i\gamma/2}$
(4)

with $\frac{m}{4\pi a_{\text{bg}}}=\frac{1}{g}+\frac{m\Lambda}{2\pi^{2}}$,
$|\alpha_{\text{re}}|^{2}=\frac{|\alpha|^{2}}{(1+mg\Lambda/(2\pi^{2}))^{2}}$,
$\nu_{\text{re}}=\nu-\frac{m\Lambda|\alpha|^{2}}{2\pi^{2}+mg\Lambda}$.

The above formula is exactly eq. (1) if we take the limit $\Lambda\to\infty$,
and write $\Gamma(I)=\frac{m|\alpha_{\text{re}}|^{2}}{4\pi a_{\text{bg}}}$ gam
.

For 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 $\nu$ and density $I$. To be more clear, as shown in
Fig. 1(b), at fixed $\Gamma(I)$, the trajectory of $a_{c}^{-1}(E=0)$ form a
circle tangent to the real-axis at $a_{\text{bg}}^{-1}$ with radius
$\frac{\Gamma(I)}{\gamma a_{\text{bg}}}$ in the upper half complex when tuning
$\nu\in(-\infty,+\infty)$. Therefore, as we gradually increase $\Gamma(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 magnetic field.
However, to achieve a general complex scattering length, we need at least two
adjustable parameters. Thus we choose optical resonance model rather than the
magnetic Feshbach resonance.

##  III 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 Chin  _et
al._ (2010). The Hermitian Hamiltonian of a Bose gas may be written as

$\displaystyle\hat{H}=\sum_{\mathbf{k}}\epsilon_{\mathbf{k}}\hat{a}_{\mathbf{k}}^{\dagger}\hat{a}_{\mathbf{k}}+\frac{g}{2\Omega}\sum_{\mathbf{k},\mathbf{k}^{\prime},\mathbf{p}}\hat{a}^{\dagger}_{\mathbf{k+p}}\hat{a}^{\dagger}_{\mathbf{k}^{\prime}-\mathbf{p}}\hat{a}_{\mathbf{k}^{\prime}}\hat{a}_{\mathbf{k}},$
(5)

where $\hat{a}_{\mathbf{k}}$ is the annihilation operator for the atoms with
momentum $\mathbf{k}$, $\epsilon_{\mathbf{k}}$ is the kinetic energy, $g$ is
the coupling constant, $\Omega$ 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 Wang  _et al._
(2021), 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
on 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 Breuer  _et al._ (2002); Wang  _et al._ (2021)

$\partial_{t}\hat{\rho}=-i(\hat{H}_{\text{eff}}\hat{\rho}-\hat{\rho}\hat{H}_{\text{eff}}^{\dagger})+\mathcal{J}\hat{\rho},$
(6)

with the non-Hermitian effective Hamiltonian

$\displaystyle\hat{H}_{\text{eff}}=H-\frac{i\gamma_{b}}{2\Omega}\sum_{\mathbf{k},\mathbf{k}^{\prime},\mathbf{p}}\hat{a}^{\dagger}_{\mathbf{k+p}}\hat{a}^{\dagger}_{\mathbf{k}^{\prime}-\mathbf{p}}\hat{a}_{\mathbf{k}^{\prime}}\hat{a}_{\mathbf{k}},$
(7)

and the recycling term

$\displaystyle\mathcal{J}\hat{\rho}=\frac{\gamma_{b}}{\Omega}\sum_{\mathbf{k},\mathbf{k}^{\prime},\mathbf{p}}\hat{a}_{\mathbf{k}^{\prime}}\hat{a}_{\mathbf{k}}\hat{\rho}\hat{a}^{\dagger}_{\mathbf{k^{\prime}-p}}\hat{a}^{\dagger}_{\mathbf{k}+\mathbf{p}},$
(8)

This master equation describes Bose atoms $\hat{a}$ interacting with each
other at a bare coupling constant $g_{b}$ while decaying to the environment at
a bare coupling constant $\gamma_{b}$ though a two-particle losses process. To
regularize the contact interaction, we can use the renormalization relation
Wang  _et al._ (2021)

$\frac{1}{g_{b}-i\gamma_{b}}=\frac{1}{g-i\gamma}-\frac{1}{\Omega}\sum_{\textbf{k}}\frac{1}{2\epsilon_{\textbf{k}}},$
(9)

and define $g-i\gamma\equiv\frac{4\pi 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
$\hat{a}_{\mathbf{0}},\hat{a}_{\mathbf{0}}^{\dagger}$ with the square root of
the zero-momentum atom number
$N-\sum_{\mathbf{k}\neq\mathbf{0}}\hat{a}_{\mathbf{k}}^{\dagger}\hat{a}_{\mathbf{k}}$,
then we obtain an approximate master equation for the reduced density matrix
$\hat{\rho}^{\prime}$ as
$\partial_{t}\hat{\rho}^{\prime}\approx\mathcal{L}_{B}(\hat{\rho}^{\prime})$
with

$\mathcal{L}_{B}(\hat{\rho}^{\prime})=-i[\hat{H}_{B},\hat{\rho}^{\prime}]-2\gamma
n\sum_{\textbf{k}\neq
0}\\{\hat{a}_{\textbf{k}}^{\dagger}\hat{a}_{\textbf{k}},\hat{\rho}^{\prime}\\}+4\gamma
n\sum_{\textbf{k}\neq
0}\hat{a}_{\textbf{k}}\hat{\rho}^{\prime}\hat{a}_{\textbf{k}}^{\dagger}$ (10)

and

$\hat{H}_{B}=\frac{gnN}{2}+\sum_{\textbf{k}\neq
0}((\epsilon_{k}+gn)\hat{a}_{\textbf{k}}^{\dagger}\hat{a}_{\textbf{k}}+\frac{gn-i\gamma
n}{2}\hat{a}_{\textbf{k}}^{\dagger}\hat{a}_{-\textbf{k}}^{\dagger}+h.c.),$

where $N$ is the total atom number and $n=\frac{N}{\Omega}$ is the density.
Here we use renormalized parameters $g,\gamma$ to replace the bare parameters
$g_{b},\gamma_{b}$ like in conventional Bogoliubov Hamiltonian Bogoliubov
(1947). 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 $n(t)=n_{0}(1+2\gamma n_{0}t)^{-1}$.

The conventional Bogoliubov Hamiltonian is the linear combination of three
operators
$\hat{A}^{0}_{\textbf{k}}=\frac{1}{2}(\hat{N}_{\textbf{k}}+\hat{N}_{-\textbf{k}}+1)$,
$\hat{A}^{1}_{\textbf{k}}=\frac{1}{2}(\hat{a}_{\textbf{k}}^{\dagger}\hat{a}_{-\textbf{k}}^{\dagger}+h.c.)$
and
$\hat{A}^{2}_{\textbf{k}}=\frac{1}{2i}(\hat{a}_{\textbf{k}}^{\dagger}\hat{a}_{-\textbf{k}}^{\dagger}-h.c.)$,
these operators form a closed SU(1,1) algebra Chen  _et al._ (2020); Cheng and
Shi (2021); Lv  _et al._ (2020),

$[\hat{A}^{1}_{\textbf{k}},\hat{A}^{2}_{\textbf{k}}]=-i\hat{A}^{0}_{\textbf{k}},[\hat{A}^{2}_{\textbf{k}},\hat{A}^{0}_{\textbf{k}}]=i\hat{A}^{1}_{\textbf{k}},[\hat{A}^{0}_{\textbf{k}},\hat{A}^{1}_{\textbf{k}}]=i\hat{A}^{2}_{\textbf{k}}.$
(11)

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 $\mathbf{A}=(\langle A^{\mathbf{k}}_{0}\rangle,\langle
A^{\mathbf{k}}_{2}\rangle,\langle A^{\mathbf{k}}_{2}\rangle)^{\text{T}}$ are
still closed, satisfying

$\displaystyle\dot{\mathbf{A}}=-2\left(\begin{array}[]{ccc}2\gamma n&\gamma
n&-gn\\\ \gamma n&2\gamma n&\epsilon_{\mathbf{k}}+gn\\\
-gn&-\epsilon_{\mathbf{k}}-gn&2\gamma
n\end{array}\right)\mathbf{A}+\left(\begin{array}[]{c}2\gamma n\\\ 0\\\
0\end{array}\right).$ (18)

Here, we focused on the short-time dynamics where $\gamma nt\ll 1$ such that
we may approximate $n$ as a constant at a fixed time t. We can diagonalize the
$3\times 3$ matrix in eq. (18) to obtain the eigenvalue
$-2i\xi_{\textbf{k},i}$, where the quasi-steady state eigenvalue is

$\xi_{\mathbf{k},0}=-2i\gamma n,\quad\xi_{\mathbf{k},(1,2)}=-2i\gamma
n\pm\sqrt{\epsilon_{\mathbf{k}}^{2}+2g_{0}n\epsilon_{\mathbf{k}}-\gamma_{0}^{2}n^{2}}.$

And a phase boundary between stable and unstable BEC can be obtained by this
quasi-steady state eigenvalue Wang  _et al._ (2021).

_Lee-Huang-Yang correction-_ Furthermore, we will give a detailed derivation
for the Lee-Huang-Yang correction Lee  _et al._ (1957); Lee and Yang (1957)
for this open system. Consider quasi-steady-state eigenvalue, up to
$O\left(\gamma nt\right)$, simple solution of eq. (18) can be obtained

$\displaystyle\mathbf{A}(t)=\mathbf{A}_{\text{s}}+e^{-4\gamma
nt}\left(\mathbf{B}+\mathbf{C}\cos 2\xi_{\mathbf{k}}t+\mathbf{D}\sin
2\xi_{\mathbf{k}}t\right)$ (19)

constant vectors $\mathbf{B,C,D}$ can be determined by the initial value of
$\mathbf{A}$, and

$\displaystyle\mathbf{A}^{\mathbf{T}}_{\text{s}}=\left(\begin{array}[]{c}\frac{1}{2}+\frac{(g^{2}+\gamma^{2})n^{2}}{2\epsilon_{\mathbf{k}}^{2}+4gn\epsilon_{\mathbf{k}}+6\gamma^{2}n^{2}}\\\
-\frac{gn\epsilon_{\mathbf{k}}+g^{2}n^{2}+2\gamma^{2}n^{2}}{2\epsilon_{\mathbf{k}}^{2}+4gn\epsilon_{\mathbf{k}}+6\gamma^{2}n^{2}}\\\
-\frac{\gamma
n(\epsilon_{\mathbf{k}}-gn)}{2\epsilon_{\mathbf{k}}^{2}+4gn\epsilon_{\mathbf{k}}+6\gamma^{2}n^{2}}\end{array}\right).$
(20)

We see that $\mathbf{A}$ never reaches the steady value
$\mathbf{A}_{\text{s}}$ because the asymptotic solution eq. (19) only works
for time interval $0\leq t\ll\hbar/\gamma 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 $\mathbf{A}$ reaches this steady value, we see that
$n_{\mathbf{k}}+n_{\mathbf{-k}}=\langle A_{0}^{\mathbf{k}}\rangle-\frac{1}{2}$
also becomes steady. Thus we get the total particle losses rate
$\frac{dN}{dt}=\frac{dN_{0}}{dt}$. And we may calculate $dN_{0}/dt$ by

$\displaystyle\frac{d{N}_{0}}{dt}=\partial_{t}\text{tr}\left(\hat{\rho}a_{0}^{\dagger}a_{0}\right)=\text{tr}\left(\mathcal{L}(\hat{\rho})a^{\dagger}_{0}a_{0}\right)=\text{tr}\left(\hat{\rho}\mathcal{L}^{\prime}(a_{0}^{\dagger}a_{0})\right),$

with $\mathcal{L}^{\prime}$ defined by

$\displaystyle\mathcal{L}^{\prime}(\hat{O})\equiv
i\left[\hat{H},\hat{O}\right]-\frac{\gamma_{b}}{2V}\sum_{\mathbf{k,k^{\prime},p}}\left\\{a^{\dagger}_{\mathbf{k+p}}a^{\dagger}_{\mathbf{k^{\prime}-p}}a_{\mathbf{k^{\prime}}}a_{\mathbf{k}},\hat{O}\right\\}$
$\displaystyle+\frac{\gamma_{b}}{2V}\sum_{\mathbf{k,k^{\prime},p}}a^{\dagger}_{\mathbf{k+p}}a^{\dagger}_{\mathbf{k^{\prime}-p}}\hat{O}a_{\mathbf{k^{\prime}}}a_{\mathbf{k}}.$

It is then straightforward to show that

$\displaystyle\mathcal{L}^{\prime}(a^{\dagger}_{0}a_{0})\simeq-2\gamma_{b}nN-2\sum_{\mathbf{k}\neq
0}\left(\gamma_{b}nA_{1}^{\mathbf{k}}+g_{b}nA_{2}^{\mathbf{k}}\right).$ (21)

This leads to

$\displaystyle\frac{dN_{0}}{dt}=-2\gamma_{b}nN-2\sum_{\mathbf{k}\neq
0}\left(\gamma_{b}n\langle A_{1}^{\mathbf{k}}\rangle+g_{b}n\langle
A_{2}^{\mathbf{k}}\rangle\right).$ (22)

From this equation, we see that the first term corresponds to the mean-field
decay which leads to $n\simeq n(1+2\gamma nt)^{-1}$ as mentioned previously.
After insert in the steady value, we have

$\displaystyle\frac{dN_{0}}{dt}=-2\gamma_{b}nN+2\gamma
n\sum_{\mathbf{k}}\frac{gn\epsilon_{\mathbf{k}}+\gamma^{2}n^{2}}{\epsilon_{\mathbf{k}}^{2}+2gn\epsilon_{\mathbf{k}}+3\gamma^{2}n^{2}}.$
(23)

We immediately find that the summation diverges for large momentum. To cure
this divergence, we use again the renormalization relation eq. (9), which
gives

$\displaystyle g_{b}-i\gamma_{b}$
$\displaystyle=\frac{1}{\frac{1}{g-i\gamma}-\frac{1}{V}\sum_{\mathbf{k}}\frac{1}{2\epsilon_{\mathbf{k}}}}$
(24)
$\displaystyle=(g-i\gamma)\left(1+\frac{g-i\gamma}{V}\sum_{\mathbf{k}}\frac{1}{2\epsilon_{\mathbf{k}}}+\ldots\right).$

We thus obtain

$\displaystyle\gamma_{b}\simeq\gamma+\frac{g\gamma}{V}\sum_{\mathbf{k}}\frac{1}{\epsilon_{\mathbf{k}}}.$
(25)

Using this expansion, we have

$\displaystyle\frac{dN}{dt}=\frac{dN_{0}}{dt}$ (26) $\displaystyle=-2\gamma
nN+2\gamma
n\sum_{\mathbf{k}}\left(\frac{gn\epsilon_{\mathbf{k}}+\gamma^{2}n^{2}}{\epsilon_{\mathbf{k}}^{2}+2gn\epsilon_{\mathbf{k}}+3\gamma^{2}n^{2}}-\frac{g_{0}n}{\epsilon_{\mathbf{k}}}\right)$
$\displaystyle=-2\gamma nN\left[1+c_{\theta}(n|a_{c}|^{3})^{1/2}\right],$

with $c_{\theta}$ a constant that only depends on the argument of the complex
scattering length $a_{c}$,

$\displaystyle c_{\theta}=4\sqrt{2\pi}\left(\frac{\cos
2\theta}{\sqrt{\cos(\theta-\pi/3)}}+2\cos\theta\sqrt{\cos(\theta-\pi/3)}\right).$

Here $\theta$ is defined as
$\theta=-\arg(a_{c})=-\arg(a_{r}+ia_{i})\in[0,\pi]$.

We note that the $n|a_{c}|^{3}$ term in eq. (26) is an analog of the
celebrated Lee-Huang-Yang correction to the ground state energy of a weakly
interacting Bose gas Lee and Yang (1957); Lee  _et al._ (1957),

$\displaystyle
E_{0}=\frac{g_{0}nN}{2}\left[1+\frac{128}{15\pi^{1/2}}(na_{r}^{3})^{1/2}+O(n\log
n)\right].$ (27)

With the help of renormalization relation, we can eliminate the divergence eq.
(22) and obtain a physical result of particle losses rate.

##  IV 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 which can
capture the interacting and dissipation features of the open bosonic system.

Figure 2: (a) Schematic of double-well toy model with interacting and two-body
losses, $\epsilon$ represent the energy detuning. (b) The phase diagram of
double-well model on $\frac{\gamma n}{\epsilon}-\frac{gn}{\epsilon}$ plane.
When $\frac{gn}{\epsilon}<-\frac{1+3(\frac{\gamma n}{\epsilon})^{2}}{2}$,
quasi-steady state eigenvalue has positive imaginary part, which represents
the exponential growth of non-condensate particle number, the system is in the
unstable phase. When $\frac{gn}{\epsilon}>-\frac{1+3(\frac{\gamma
n}{\epsilon})^{2}}{2}$, number of the minority particle is always much smaller
than condensate number during the time evolution, the system is in the stable
phase.

_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 $\hat{b}_{r}(\hat{b}_{l})$. The evolution of this system is governed by the
master equation

$\partial_{t}\hat{\rho}=-i[\hat{H}_{\text{dw}},\hat{\rho}]-\gamma\sum_{i=r,l}\left(\\{\hat{\Gamma}_{i}^{\dagger}\hat{\Gamma}_{i},\rho\\}+2\hat{\Gamma}_{i}\rho\hat{\Gamma}^{\dagger}_{i}\right),$
(28)

where the Hamiltonian is

$\hat{H}_{\text{dw}}=-t(\hat{b}_{r}^{\dagger}\hat{b}_{l}+\hat{b}_{l}^{\dagger}\hat{b}_{r})+g(\hat{b}_{r}^{\dagger}\hat{b}_{r}^{\dagger}\hat{b}_{r}\hat{b}_{r}+\hat{b}_{l}^{\dagger}\hat{b}_{l}^{\dagger}\hat{b}_{l}\hat{b}_{l}),$
(29)

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
$\hat{\Gamma}_{r}=\hat{b}_{r}\hat{b}_{r},\hat{\Gamma}_{l}=\hat{b}_{l}\hat{b}_{l}$,
with $\gamma$ the decay rate.

Corresponding to the Bose gas model, we consider the “zero momentum” mode
$\hat{a}_{0}=\frac{1}{\sqrt{2}}(\hat{b}_{r}+\hat{b}_{l})$ as majority part
while the “non-zero momentum” mode
$\hat{a}_{1}=\frac{1}{\sqrt{2}}(\hat{b}_{r}-\hat{b}_{l})$ as small depletion
part. We can shift the energy for the mode $\hat{a}_{0}$ to zero without loss
of generality and the mode $\hat{a}_{1}$ has a detuning $\epsilon=2t$ from the
mode $\hat{a}_{0}$. Omitting the interacting or losses terms only including
minority part, such as
$\hat{a}_{1}^{\dagger}\hat{a}_{1}^{\dagger}\hat{a}_{1}\hat{a}_{1}$, the master
equation is then given by

$\partial_{t}\hat{\rho}=-i(\hat{H}_{a}^{\text{eff}}\hat{\rho}-\hat{\rho}\hat{H}_{a}^{\text{eff}\dagger})+\mathcal{J}\hat{\rho},$
(30)

where the effective Hamiltonian $\hat{H}_{a}^{\text{eff}}$ for
$\hat{a}_{0},\hat{a}_{1}$is

$\displaystyle\hat{H}_{a}^{\text{eff}}=$
$\displaystyle\epsilon\hat{a}_{1}^{\dagger}\hat{a}_{1}+\frac{g-i\gamma}{2}(\hat{a}_{0}^{\dagger}\hat{a}_{0}^{\dagger}\hat{a}_{0}\hat{a}_{0}+\hat{a}_{0}^{\dagger}\hat{a}_{0}^{\dagger}\hat{a}_{1}\hat{a}_{1}$
(31)
$\displaystyle+\hat{a}_{1}^{\dagger}\hat{a}_{1}^{\dagger}\hat{a}_{0}\hat{a}_{0}+4\hat{a}_{0}^{\dagger}\hat{a}_{0}\hat{a}_{1}^{\dagger}\hat{a}_{1}),$

and the recycling term in eq. (30) are similar to Lindbladian of Bose gas with
two-body losses,

$\mathcal{J}\hat{\rho}=\gamma(\hat{a}_{0}\hat{a}_{0}\hat{\rho}\hat{a}_{0}^{\dagger}\hat{a}_{0}^{\dagger}+\hat{a}_{1}\hat{a}_{1}\hat{\rho}\hat{a}_{0}^{\dagger}\hat{a}_{0}^{\dagger}+\hat{a}_{0}\hat{a}_{0}\hat{\rho}\hat{a}_{1}^{\dagger}\hat{a}_{1}^{\dagger}+4\hat{a}_{0}\hat{a}_{1}\hat{\rho}\hat{a}_{0}^{\dagger}\hat{a}_{1}^{\dagger}).$

In this model, the effective Hamiltonian $\hat{H}_{a}^{\text{eff}}$ conserves
the particle number
$\hat{n}=\hat{a}_{0}^{\dagger}\hat{a}_{0}+\hat{a}_{1}^{\dagger}\hat{a}_{1}$
while the recycling term always annihilates two particles, thus the master
equation can be decomposed to a series of hierarchy equations for
$\hat{\rho}_{i,j}$ Wang  _et al._ (2021)

$\partial_{t}\hat{\rho}_{i,j}=-i(\hat{H}_{a}^{\text{eff}}\hat{\rho}_{i,j}-\hat{\rho}_{i,j}\hat{H}_{a}^{\text{eff}\dagger})+\mathcal{J}\hat{\rho}_{i+2,j+2},$
(32)

where $\hat{\rho}_{i,j}=\hat{P}_{i}\hat{\rho}\hat{P}_{j}$, with $\hat{P}_{i}$
the projection operator for the $i$-particle Fock space. Then for the
dynamical evolution of eq. (30) with initial density matrix a pure state of
particle number $N$, only the projected density matrix $\hat{\rho}_{N,N}$,
$\hat{\rho}_{N-2,N-2}$,$\hat{\rho}_{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 $\hat{a}_{0}$ and
$\hat{a}_{0}^{\dagger}$ with the square root of the condensate atoms number
$n-\hat{a}_{1}^{\dagger}\hat{a}_{1}$, then we can obtain the approximated
Lindbladian as

$\displaystyle\hat{\mathcal{L}}_{\text{dw}}^{B}\hat{\rho}$
$\displaystyle=-i[\hat{H}_{\text{dw}}^{B},\hat{\rho}]-\frac{\gamma}{2}\\{n\hat{a}_{1}\hat{a}_{1}+n\hat{a}_{1}^{\dagger}\hat{a}_{1}^{\dagger}+4n\hat{a}_{1}^{\dagger}\hat{a}_{1},\hat{\rho}\\}$
(33)
$\displaystyle+\gamma(n\hat{a}_{1}\hat{a}_{1}\hat{\rho}+n\hat{\rho}\hat{a}_{1}^{\dagger}\hat{a}_{1}^{\dagger}+4n\hat{a}_{1}\hat{\rho}\hat{a}_{1}^{\dagger}),$

where the Hamiltonian $\hat{H}_{\text{dw}}^{B}$ is

$\hat{H}_{\text{dw}}^{B}=\epsilon\hat{a}_{1}^{\dagger}\hat{a}_{1}+\frac{g}{2}(n^{2}+n\hat{a}_{1}\hat{a}_{1}+n\hat{a}_{1}^{\dagger}\hat{a}_{1}^{\dagger}+2n\hat{a}_{1}^{\dagger}\hat{a}_{1}),$
(34)

and time-dependent particle number is given by $n=N(1+2\gamma Nt)^{-1}$.
Similar to eq. (18) in the Bose gas case, the expectation value of the three
SU(1,1) operators
$\hat{A}_{\text{dw}}^{1}=(\hat{a}_{1}^{\dagger}\hat{a}_{1}+\hat{a}_{1}\hat{a}_{1}^{\dagger})/2,\hat{A}_{\text{dw}}^{2}=(\hat{a}_{1}^{\dagger}\hat{a}_{1}^{\dagger}+\hat{a}_{1}\hat{a}_{1})/2$
and
$\hat{A}_{\text{dw}}^{3}=(\hat{a}_{1}^{\dagger}\hat{a}_{1}^{\dagger}-\hat{a}_{1}\hat{a}_{1})/2i$
form a closed set of differential equations,

$\displaystyle\dot{\mathbf{A}}_{\text{dw}}=-2\left(\begin{array}[]{ccc}2\gamma
n&\gamma n&-gn\\\ \gamma n&2\gamma n&\epsilon+gn\\\ -gn&-\epsilon-gn&2\gamma
n\end{array}\right)\mathbf{A}_{\text{dw}}+\left(\begin{array}[]{c}\gamma n\\\
0\\\ 0\end{array}\right),$ (41)

with
$\mathbf{A}_{\text{dw}}=(\langle\hat{A}_{\text{dw}}^{1}\rangle,\langle\hat{A}_{\text{dw}}^{2}\rangle,\langle\hat{A}_{\text{dw}}^{3}\rangle)^{\text{T}}$.
The three eigenvalues $-2i\xi_{i}$ for the $3\times 3$ matrix in eq. (41) are

$\xi_{0}=-2i\gamma n,\quad\xi_{(1,2)}=-2i\gamma
n\pm\sqrt{\epsilon^{2}+2gn\epsilon-\gamma^{2}n^{2}}.$ (42)

When the imaginary part of $\xi_{1}$ is larger than zero, the corresponding
eigenvalue becomes positive number, thus the particle number $\hat{n}_{1}$
grows exponentially in short time leading to an unstable dynamics. The phase
boundary for this stability is given by
$\frac{gn}{\epsilon}=-\frac{1+3(\frac{\gamma n}{\epsilon})^{2}}{2}$ as shown
in Fig. 2(b).

Figure 3: Quench dynamics under time evolution using two methods with total
particle number $N=100$ and energy in unit of $\epsilon$. (solid: exact time
evolution governed by Lindbladian eq. (30); dashed: approximated Lindbladian
eq. (33)). Initial condition is all particles condensate at groundstate,
$n_{0}=N$. (a),(b)System at stable phase with coupling strength
$gN=0.52\epsilon$ and dissipation strength $\gamma N=0.20\epsilon$. (c),(d)
System at stable phase with $gN=-0.52\epsilon$ and $\gamma N=0.02\epsilon$.

_Numerical results -_ We now check the accuracy of the Bogoliubov
approximation by comparing the numerical simulation result from eq. (30) and
eq. (33). The evolution begins with an initial state with $N=100$ particles
occupying the mode $\hat{a}_{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 $\hat{n}_{1}$
predicted by the Bogoliubov approximation eq. (33) deviates from the exact
time evolution governed by eq. (30) 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.

##  V 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 Prosen (2008); Prosen
and Seligman (2010). However, this cannot be directly applied to our problem
because of the time-dependent $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,\mathbb{C})$ dynamical symmetry. And this dynamical symmetry
provides a relatively simple way Wei and Norman (1963); Scopa  _et al._ (2019)
to construct exact solutions to the time-dependent Lindbladian algebraically.

_Symmetry-_ The Lindbladian $\mathcal{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 $a_{\textbf{k}}\rightarrow
a_{\textbf{k}}e^{i\phi},a_{\textbf{k}}^{\dagger}\rightarrow
a_{\textbf{k}}^{\dagger}e^{-i\phi}$. The corresponding conserved superoperator
is thus $\tilde{Q}=[\hat{N},\circ]$. While for the $Z_{2}$ symmetry, it can be
verified by realizing $\mathcal{L}_{B}$ is invariant under
$\hat{a}_{\textbf{k}}\rightarrow-\hat{a}_{\textbf{k}}$ or
$\hat{a}_{\textbf{k}}^{\dagger}\rightarrow-\hat{a}_{\textbf{k}}^{\dagger}$.

However, these two symmetries alone are not enough for obtaining 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 contains the Lindbladian $\mathcal{L}_{B}$ itself.

_Closed algebra-_ To find this algebra, we note that $\mathcal{L}_{B}$ can be
expressed in the superoperator-formalism as

$\mathcal{L}_{B}(\hat{\rho}^{\prime})=\left(\sum_{\mathbf{k},k_{z}\geq
0}\tilde{\mathcal{L}}_{\mathbf{k}}\right)\circ\hat{\rho}^{\prime}.$ (43)

Here $\tilde{\mathcal{L}}_{\mathbf{k}}$ is a superoperator act only on modes k
and $-\textbf{k}$, which is a linear combination of seven quadratic
superoperators. We label these superoperators by $\tilde{h}^{\mathbf{k}}_{i},\
i=1,2,\ldots 7$. They are defined by

$\displaystyle\tilde{h}^{1}_{\mathbf{k}}\circ\hat{\rho}^{\prime}=(\hat{a}_{\mathbf{k}}^{\dagger}\hat{a}_{\mathbf{k}}+\hat{a}_{-\mathbf{k}}\hat{a}_{-\mathbf{k}}^{\dagger})\hat{\rho}^{\prime},\quad\tilde{h}^{2}_{\mathbf{k}}\circ\hat{\rho}^{\prime}=2\hat{a}_{\mathbf{k}}^{\dagger}\hat{a}_{-\mathbf{k}}^{\dagger}\hat{\rho}^{\prime},$
$\displaystyle\tilde{h}^{3}_{\mathbf{k}}\circ\hat{\rho}^{\prime}=2\hat{a}_{\mathbf{k}}\hat{a}_{-\mathbf{k}}\hat{\rho}^{\prime},\quad\tilde{h}^{4}_{\mathbf{k}}\circ\hat{\rho}^{\prime}=\hat{\rho}^{\prime}(\hat{a}_{\mathbf{k}}^{\dagger}\hat{a}_{\mathbf{k}}+\hat{a}_{-\mathbf{k}}\hat{a}_{-\mathbf{k}}^{\dagger}),$
$\displaystyle\tilde{h}^{5}_{\mathbf{k}}\circ\hat{\rho}^{\prime}=2\hat{\rho}^{\prime}\hat{a}_{\mathbf{k}}^{\dagger}\hat{a}_{-\mathbf{k}}^{\dagger},\quad\tilde{h}^{6}_{\mathbf{k}}\circ\hat{\rho}^{\prime}=2\hat{\rho}^{\prime}\hat{a}_{\mathbf{k}}\hat{a}_{\mathbf{k}},$
$\displaystyle\tilde{h}^{7}_{\mathbf{k}}\circ\hat{\rho}^{\prime}=\hat{a}_{\mathbf{k}}\hat{\rho}^{\prime}\hat{a}_{\mathbf{k}}^{\dagger}+\hat{a}_{-\mathbf{k}}\hat{\rho}^{\prime}\hat{a}_{-\mathbf{k}}^{\dagger}.$

Among these seven superoperators, $\tilde{h}_{\mathbf{k}}^{i},\ i=1,2,\ldots
6$ represent the (anti-)commutators between the (anti-)Hermitian parts of the
effective Hamiltonian and the density matrix $\hat{\rho}$, and the last one
represents the recycling term of $\tilde{\mathcal{L}}_{\mathbf{k}}$.

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

$[\tilde{h}^{5}_{\mathbf{k}},\tilde{h}^{7}_{\mathbf{k}}]\circ\hat{\rho}^{\prime}=-2(\hat{a}_{\mathbf{k}}\hat{\rho}^{\prime}\hat{a}_{-\mathbf{k}}^{\dagger}+\hat{a}_{-\mathbf{k}}\hat{\rho}^{\prime}\hat{a}_{\mathbf{k}}^{\dagger}).$
(44)

which is an superoperator that can not be written in the form of
$\sum_{i=1}^{7}\tilde{h}_{\mathbf{k}}^{i}$.

It is thus necessary to introduce three more superoperators
$\tilde{h}_{\mathbf{k}}^{i},\ i=8,9,10$ to close the algebra. They are

$\displaystyle\tilde{h}^{8}_{\mathbf{k}}\circ\hat{\rho}^{\prime}=\hat{a}_{-\mathbf{k}}^{\dagger}\hat{\rho}^{\prime}\hat{a}_{\mathbf{k}}^{\dagger}+\hat{a}_{\mathbf{k}}^{\dagger}\hat{\rho}^{\prime}\hat{a}_{-\mathbf{k}}^{\dagger},$
$\displaystyle\tilde{h}^{9}_{\mathbf{k}}\circ\hat{\rho}^{\prime}=\hat{a}_{\mathbf{k}}\hat{\rho}^{\prime}\hat{a}_{-\mathbf{k}}^{\dagger}+\hat{a}_{-\mathbf{k}}\hat{\rho}^{\prime}\hat{a}_{\mathbf{k}}^{\dagger},$
$\displaystyle\tilde{h}^{10}_{\mathbf{k}}\circ\hat{\rho}^{\prime}=\hat{a}_{\mathbf{k}}^{\dagger}\hat{\rho}^{\prime}\hat{a}_{\mathbf{k}}+\hat{a}_{-\mathbf{k}}^{\dagger}\hat{\rho}^{\prime}\hat{a}_{\mathbf{k}}.$

All of the superoperators $\tilde{h}_{\mathbf{k}}^{i},\ i=1,2,\ldots,10$ now
form a closed algebra sev , 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 $-\textbf{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,\mathbb{C})$. Before that, we will formally write the exact solution of
$\tilde{\mathcal{L}}_{B}$.

| $\tilde{h}^{1}]$ | $\tilde{h}^{2}]$ | $\tilde{h}^{3}]$ | $\tilde{h}^{4}]$ | $\tilde{h}^{5}]$ | $\tilde{h}^{6}]$ | $\tilde{h}^{7}]$ | $\tilde{h}^{8}]$ | $\tilde{h}^{9}]$ | $\tilde{h}^{10}]$
---|---|---|---|---|---|---|---|---|---|---
$[\tilde{h}^{1},$ | 0 | $2\tilde{h}^{2}$ | $-2\tilde{h}^{3}$ | 0 | 0 | 0 | $-\tilde{h}^{7}$ | $\tilde{h}^{8}$ | $-\tilde{h}^{9}$ | $\tilde{h}^{10}$
$[\tilde{h}^{2},$ |  | 0 | $-4\tilde{h}^{1}$ | 0 | 0 | 0 | $-2\tilde{h}^{8}$ | 0 | $-2\tilde{h}^{10}$ | 0
$[\tilde{h}^{3},$ |  |  | 0 | 0 | 0 | 0 | 0 | $2\tilde{h}^{7}$ | 0 | $2\tilde{h}^{9}$
$[\tilde{h}^{4},$ |  |  |  | 0 | $2\tilde{h}^{5}$ | $-2\tilde{h}^{6}$ | $-\tilde{h}^{7}$ | $-\tilde{h}^{8}$ | $\tilde{h}^{9}$ | $\tilde{h}^{10}$
$[\tilde{h}^{5},$ |  |  |  |  | 0 | $-4\tilde{h}^{4}$ | $-2\tilde{h}^{9}$ | $-2\tilde{h}^{10}$ | 0 | 0
$[\tilde{h}^{6},$ |  |  |  |  |  | 0 | 0 | 0 | $2\tilde{h}^{7}$ | $2\tilde{h}^{8}$
$[\tilde{h}^{7},$ |  |  |  |  |  |  | 0 | $\tilde{h}^{6}$ | $\tilde{h}^{3}$ | $\tilde{h}^{1}+\tilde{h}^{4}$
$[\tilde{h}^{8},$ |  |  |  |  |  |  |  | 0 | $\tilde{h}^{1}-\tilde{h}^{4}$ | $\tilde{h}^{2}$
$[\tilde{h}^{9},$ |  |  |  |  |  |  |  |  | 0 | $\tilde{h}^{5}$
$[\tilde{h}^{10},$ |  |  |  |  |  |  |  |  |  | 0
Table 1: Commutation relation table of ten superoperators $\tilde{h}^{i}\
i=1,2,\ldots,10$. The superoperators
$\tilde{h}^{1},\tilde{h}^{2},\tilde{h}^{3}$ or $\tilde{h}^{4},\tilde{h}^{5},$
$\tilde{h}^{6}$ can form two $su(1,1)$ algebra, which is the result of the
SU(1,1) dynamical symmetry of Bogoliubov Hamiltonian in closed system.

_Exact solution-_ Based on this closed algebra structure, the dynamical
problem of time-dependent Lindbladian can be solved exactly Wei and Norman
(1963); Scopa  _et al._ (2019). First of all, formally solve the master
equation $\hat{\rho}$ as
$\partial_{t}\hat{\rho}=\tilde{\mathcal{L}}_{B}\circ\hat{\rho}$ , we can
denote the solution as

$\hat{\rho}(t)=\tilde{\Lambda}(t,t_{0})\circ\hat{\rho}(t_{0}).$ (45)

$\tilde{\Lambda}(t,t_{0})$ is a dynamical map from time $t_{0}$ to time $t$.
In general, this mapping is a semi-group which satisfies
$\tilde{\Lambda}(t_{2},t_{0})=\tilde{\Lambda}(t_{2},t_{1})\tilde{\Lambda}(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
$\tilde{\mathcal{L}}=\sum_{l=1}^{7}u_{l}(t)\tilde{h}^{l}$ hence we can
parametrise the dynamical map by these ten generators,

$\tilde{\Lambda}(t,t_{0})=\prod\limits_{i=1}\limits^{10}e^{g_{i}(t)\tilde{h}^{i}}.$
(46)

Substitute this equation back to the master equation, we can get

$\partial_{t}\prod\limits_{i=1}\limits^{10}e^{g_{i}(t)\tilde{h}^{i}}=\sum_{l=1}^{7}u_{l}(t)\tilde{h}^{l}\prod\limits_{j=1}\limits^{10}e^{g_{j}(t)\tilde{h}^{j}},$
(47)

where $u_{l}(t)$ are time-dependent parameters defined by Hamiltonian. And
solving dynamics of time-evolution operator is equivalent to solving the
complex functions $g_{i}(t)$. Right multiply the inverse of
$\tilde{\Lambda}(t,t_{0})$ at both side of eq. (47) and write the explicit
expression of time derivative, we can obtain

$\sum_{m}\partial_{t}g_{m}(t)\prod\limits_{i=1}\limits^{m}e^{g_{i}(t)\tilde{h}^{i}}\tilde{h}^{m}e^{-g_{i}(t)\tilde{h}^{i}}=\sum_{l=1}^{7}u_{l}(t)\tilde{h}^{l}.$
(48)

Using Baker-Campbell-Hausdorf formula Hall (2013) and commutation relation in
Table. 1, we can formally write L.H.S of eq. (48) as:

$\sum_{mn}\partial_{t}g_{m}(t)\eta_{mn}\tilde{h}^{n}=\sum_{l=1}^{7}u_{l}(t)\tilde{h}^{l},$
(49)

where the $\eta_{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 the following, we show that the superoperators $\tilde{h}^{i}$ form an
algebra of $sp(4,\mathbb{C})$ 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,\mathbb{C})$.

_Map to harmonic oscillators-_ Even though the commutation relations between
$\tilde{h}_{\mathbf{k}}^{i}$ seem complicated as shown in Table. 1, it is
worth noting that they only consist quadratic superoperators. A natural
question is then whether there is an isolated system which 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
$\tilde{h}^{i}$ can be mapped to linear combinations of quadratic forms of
${a},{a}^{\dagger},{b},{b}^{\dagger}$. 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 a example, we will show
in the following that the quadratics of ladder operators $\hat{a}$ and
$\hat{b}$ form the algebra of symplectic group Sp$(4,\mathbb{C})$.

| Harmonic Oscillator | Bogoliubov Lindbladian
---|---|---
$h^{1}$ | $(a^{\dagger}a+aa^{\dagger})/2$ | $(\hat{a}^{\dagger}_{\textbf{k}}\hat{a}_{\textbf{k}}+\hat{a}_{-\textbf{k}}\hat{a}_{-\textbf{k}}^{\dagger})\circ$
$h^{2}$ | $a^{\dagger}a^{\dagger}$ | $(2\hat{a}^{\dagger}_{\textbf{k}}\hat{a}_{-\textbf{k}}^{\dagger})\circ$
$h^{3}$ | $aa$ | $(2\hat{a}_{\textbf{k}}\hat{a}_{-\textbf{k}})\circ$
$h^{4}$ | $(b^{\dagger}b+bb^{\dagger})/2$ | $\circ(\hat{a}^{\dagger}_{\textbf{k}}\hat{a}_{\textbf{k}}+\hat{a}_{-\textbf{k}}\hat{a}_{-\textbf{k}}^{\dagger})$
$h^{5}$ | $b^{\dagger}b^{\dagger}$ | $\circ(2\hat{a}^{\dagger}_{\textbf{k}}\hat{a}_{-\textbf{k}}^{\dagger})$
$h^{6}$ | $bb$ | $\circ(2\hat{a}_{\textbf{k}}\hat{a}_{-\textbf{k}})$
$h^{7}$ | $ab$ | $\hat{a}_{\textbf{k}}\circ\hat{a}_{\textbf{k}}^{\dagger}+\hat{a}_{-\textbf{k}}\circ\hat{a}_{-\textbf{k}}^{\dagger}$
$h^{8}$ | $a^{\dagger}b$ | $\hat{a}_{-\textbf{k}}^{\dagger}\circ\hat{a}^{\dagger}_{\textbf{k}}+a_{\textbf{k}}^{\dagger}\circ a^{\dagger}_{-\textbf{k}}$
$h^{9}$ | $ab^{\dagger}$ | $\hat{a}_{\textbf{k}}\circ\hat{a}_{-\textbf{k}}^{\dagger}+\hat{a}_{-\textbf{k}}\circ\hat{a}_{\textbf{k}}^{\dagger}$
$h^{10}$ | $a^{\dagger}b^{\dagger}$ | $\hat{a}_{\textbf{k}}^{\dagger}\circ\hat{a}_{\textbf{k}}+\hat{a}_{-\textbf{k}}^{\dagger}\circ\hat{a}_{-\textbf{k}}$
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,\mathbb{C})$.

Sp _$(2n,\mathbb{C})$ 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,\ldots$) have symplectic
Sp$(2n,\mathbb{C})$ dynamical symmetry. A generic quadratic form is given by
is given by

$\displaystyle\hat{M}$
$\displaystyle=\frac{1}{2}\left(\begin{array}[]{cccccc}a_{1}^{\dagger},&\ldots,&a_{n}^{\dagger},&a_{1},&\ldots,&a_{n}\end{array}\right)\left(\begin{array}[]{cc}B&A\\\
A^{\text{T}}&C\end{array}\right)\left(\begin{array}[]{c}a_{1}^{\dagger}\\\
\vdots\\\ a_{n}^{\dagger}\\\ a_{1}\\\ \vdots\\\ a_{n}\end{array}\right).$ (59)

To make the coefficient matrices unique, we require $B^{\text{T}}=B$ and
$C^{\text{T}}=C$.

Note that the coefficient matrix may be written as

$\displaystyle\left(\begin{array}[]{cc}B&A\\\
A^{\text{T}}&C\end{array}\right)=\left(\begin{array}[]{cc}A&-B\\\
C&-A^{\text{T}}\end{array}\right)\left(\begin{array}[]{cc}{}\hfil&I\\\
-I&{}\hfil\end{array}\right)\equiv M\Omega$ (66)

with $M$ satisfying

$\displaystyle\Omega M+M^{\text{T}}\Omega=0,$ (67)

which is the generator for Lie group Sp$(2n,\mathbb{C})$.

We now define new operators

$\displaystyle
b_{i}=\left\\{\begin{array}[]{ll}a_{i},&i\in\\{1,2,\ldots,n\\}\\\
a_{i-n},&i\in\\{n+1,\ldots,2n\\}\\\ \end{array}\right.$ (70) $\displaystyle
b^{i}=\left\\{\begin{array}[]{ll}a_{i}^{\dagger},&i\in\\{1,2,\ldots,n\\}\\\
a_{i-n},&i\in\\{n+1,\ldots,2n\\}\end{array}\right.$ (73)

which are related by $b_{i}=\Omega_{ij}b^{j}$ and $b^{i}=\Omega^{ij}b_{j}$
with matrices

$\displaystyle\Omega_{ij}=\left\\{\begin{array}[]{ll}\delta_{i,j-n},&i\in\\{1,2,\ldots,n\\}\\\
-\delta_{i,j+n},&i\in\\{n+1,\ldots,2n\\}\\\ \end{array}\right.$ (76)
$\displaystyle\Omega^{ij}=\left\\{\begin{array}[]{ll}-\delta_{i,j-n},&i\in\\{1,2,\ldots,n\\}\\\
\delta_{i,j+n},&i\in\\{n+1,\ldots,2n\\}\\\ \end{array}\right.$ (79)

We thus have following useful commutation relations

$\displaystyle[b_{i},b^{j}]=\delta_{i}^{j},\qquad\ \
[b^{i},b_{j}]=-\delta^{i}_{j},$ (80)
$\displaystyle[b^{i},b^{j}]=\Omega^{ij},\qquad[b_{i},b_{j}]=-\Omega_{ij}.$
(81)

With all these definitions, the generic quadratic form in Eq. (59) can then be
expressed as

$\displaystyle\hat{M}=\frac{1}{2}b^{i}M_{i}^{j}b_{j}$ (82)

with the upper (lower) index represents the column (row) index. The above
formula gives a one-to-one mapping between the quadratic form and the Lie
algebra $sp(2n,\mathbb{C})$. More importantly, it can be proved

$\displaystyle[\hat{M},\hat{N}]=\frac{1}{2}b^{i}[M,N]_{i}^{j}b_{j},$ (83)

which implies that the mapping is an isomorphism.

Now we prove $n$-mode quadratic bosonic Hamiltonian has an Sp$(2n,\mathbb{C})$
dynamical symmetry, meanwhile, in 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,\mathbb{C})$ dynamical symmetry.

Furthermore, our results are not only restricted to Bogoliubov Lindbladian.
Because quadratic operators in 2$n$-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,\mathbb{C})$ dynamical symmetry. This result
will be useful for the further analytical study of open systems.

##  VI 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 Dalfovo  _et al._ (1999); Pethick and Smith (2008); Zhai
(2021). In this section, we generalize these two descriptions to open systems
with weak two-body losses. We derive the dissipative versions of Gross-
Pitaevskii equation and hydrodynamic equations based on the Keldysh path
integral formalism.

_Keldysh formalism-_ We start from the master equation
$\partial_{t}\hat{\rho}=\mathcal{L}\hat{\rho}$ of interacting bosons subject
to two-body losses in real space,

$\mathcal{L}\hat{\rho}=\frac{1}{i}[\hat{H},\hat{\rho}]-\frac{\gamma}{2}\int_{\textbf{r}}\\{\hat{\psi}^{\dagger}(\textbf{r})\hat{\psi}^{\dagger}(\textbf{r})\hat{\psi}(\textbf{r})\hat{\psi}(\textbf{r}),\hat{\rho}\\}+\mathcal{J}\hat{\rho},$
(84)

where $\hat{\psi}(\textbf{r})$ is the bosonic annihilation operator at
position r, and

$\displaystyle\hat{H}=\int_{\mathbf{r}}\hat{\psi}^{\dagger}(\mathbf{r})\left(-\frac{\nabla^{2}}{2m}\right)\hat{\psi}^{\dagger}(\mathbf{r})+\frac{g}{2}\int_{\mathbf{r}}\hat{\psi}^{\dagger}(\textbf{r})\hat{\psi}^{\dagger}(\textbf{r})\hat{\psi}(\textbf{r})\hat{\psi}(\textbf{r})$
(85)

is Hermitian Hamiltonian of interacting bosons.

The recycling term $\mathcal{J}\hat{\rho}$ is given by

$\mathcal{J}\hat{\rho}=\gamma\int_{\textbf{r}}\hat{\psi}(\textbf{r})\hat{\psi}(\textbf{r})\hat{\rho}\psi^{\dagger}(\textbf{r})\hat{\psi}^{\dagger}(\textbf{r}).$
(86)

Using the keldysh path-integral representation, we introduce fields
$\psi_{\pm}$ and $\bar{\psi}_{\pm}$ living respectively on the time contour
$\mathcal{C}_{+}$ and $\mathcal{C}_{-}$, where time runs from $-\infty$ to
$+\infty$ and then back to $+\infty$ in the closed-contour
$\mathcal{C}_{+}\cup\mathcal{C}_{-}$. Then we can write partition function
$Z=\text{Tr}(\rho(t))$ of eq. (84),

$Z=\int\mathcal{D}[\psi_{+},\bar{\psi}_{+},\psi_{-},\bar{\psi}_{-}]e^{iS[\psi]},$
(87)

where we omit the initial condition and the Lindblad-Keldysh action Sieberer
_et al._ (2016) is

$\displaystyle S$ $\displaystyle=\int
dtd\textbf{r}\sum_{\eta=\pm}(-1)^{s_{\eta}}(\bar{\psi}_{\eta}(\textbf{r})i\partial_{t}\psi_{\eta}(\textbf{r})-H(\bar{\psi}_{\eta}(\textbf{r}),\psi_{\eta}(\textbf{r})))$
(88) $\displaystyle+\frac{i}{2}\gamma\int
dtd\textbf{r}\sum_{\eta=\pm}\bar{\psi}_{\eta}(\textbf{r})\bar{\psi}_{\eta}(\textbf{r})\psi_{\eta}(\textbf{r})\psi_{\eta}(\textbf{r})$
$\displaystyle-i\gamma\int
dtd\textbf{r}\quad\psi_{+}(\textbf{r})\psi_{+}(\textbf{r})\bar{\psi}_{-}(\textbf{r})\bar{\psi}_{-}(\textbf{r}),$

where time runs from $-\infty$ to +$\infty$, and $s_{\eta}=0$ for $\eta=+$,
$s_{\eta}=1$ for $\eta=-$.

_Saddle-point approximation-_ We consider case that almost all particles
occupies the ground state energy level, which means $N\approx N_{0}\gg 1$. In
the closed system, we can assume the system is always in coherent state and
the condensate wavefunction can be found by saddle-point equation. Similarly,
here we consider losses process is weak and slow so that we can still take the
coherent state assumption. As a result, we can take the saddle-point equation
in open system,

$\frac{\delta S}{\delta\bar{\psi}_{\pm}}=0,\frac{\delta
S}{\delta\psi_{\pm}}=0,$ (89)

then we obtain

$\displaystyle
i\partial_{t}\psi_{+}-(-\frac{\nabla^{2}}{2m}\psi_{+}+g_{c}\bar{\psi}_{+}\psi_{+}\psi_{+})=0$
(90) $\displaystyle-$ $\displaystyle
i\partial_{t}\bar{\psi}_{+}-(-\frac{\nabla^{2}}{2m}\bar{\psi}_{+}+g_{c}\bar{\psi}_{+}\psi_{+}\bar{\psi}_{+})-2i\gamma\psi_{+}\bar{\psi}_{-}\bar{\psi}_{-}=0$
$\displaystyle-$ $\displaystyle
i\partial_{t}\psi_{-}+(-\frac{\nabla^{2}}{2m}\psi_{-}+g_{c}^{*}\bar{\psi}_{-}\psi_{-}\psi_{-})-2i\gamma\psi_{+}\psi_{+}\bar{\psi}_{-}=0$
$\displaystyle
i\partial_{t}\bar{\psi}_{-}+(-\frac{\nabla^{2}}{2m}\bar{\psi}_{-}+g_{c}^{*}\bar{\psi}_{-}\psi_{-}\bar{\psi}_{-})=0,$

where the complex interacting strength $g_{c}=g-i\gamma$.

_Non-Hermitian Gross-Pitaevskii equation -_ The regularization of Keldysh
action requires the relation $\psi_{+}=\psi_{-},\bar{\psi}_{+}=\bar{\psi}_{-}$
Kamenev (2011), combining with the saddle-point equations eq. (90), we will
obtain the Gross-Pitaevskii equation under the two-body losses
$(\psi_{+}=\psi)$,

$i\partial_{t}\psi-\left(-\frac{\nabla^{2}}{2m}\psi+g_{c}|\psi|^{2}\psi+V\psi\right)=0.$
(91)

This equation substitutes the coupling parameter $g$ to complex version
$g-i\gamma$ in the conventional Gross-Pitaevskii equation. However, due to the
particle loss, the solution of this dissipative Gross-Pitaevskii equation
becomes quiet different from the conventional solution for closed system. For
example, when $V=0$, it is easy to show that we have a plane-wave like
solution given by

$\psi_{0}(\textbf{r})=\sqrt{n_{0}}e^{-i\frac{g_{c}}{2\gamma}\text{log}(1+2\gamma
n_{0}t)}.$ (92)

And in finite momentum we can also obtain an exact solution,
$\psi_{\textbf{p}}(\textbf{r})=\sqrt{n_{0}}e^{i(\textbf{p}\cdot\textbf{r}-p^{2}t/2)}e^{-i\frac{g_{c}}{2\gamma}\text{log}(1+2\gamma
n_{0}t)}$.

Then we can solve the excitation spectrum of this dissipative Gross-Pitaevskii
equation by inserting small perturbation $\psi=\psi_{0}+\delta\psi$ and its
complex conjugate and expanding to the linear order in $\delta\psi$:

$\displaystyle
i\partial_{t}\delta\psi=-\frac{1}{2}\nabla^{2}\delta\psi+2g_{c}|\psi_{0}|^{2}\delta\psi+(g-i\gamma)\psi_{0}^{2}\delta\psi^{*},$
$\displaystyle
i\partial_{t}\delta\psi^{*}=\frac{1}{2}\nabla^{2}\delta\psi^{*}-2g_{c}^{*}|\psi_{0}|^{2}\delta\psi^{*}-(g+i\gamma)\psi_{0}^{*2}\delta\psi.$

To solve this equation, we may apply a unitary transformation
$\delta\phi=e^{i\frac{g}{2\gamma}}\log(1+2\gamma n_{0}t)\delta\psi$ which
leads to a simplified equation in the “rotating wave frame”,

$\displaystyle
i\partial_{t}\delta\phi=-\frac{1}{2}\nabla^{2}\delta\phi+(2g_{c}-g)n\delta\phi+g_{c}n\delta\phi^{*},$
$\displaystyle
i\partial_{t}\delta\phi^{*}=\frac{1}{2}\nabla^{2}\delta\phi^{*}-(2g_{c}^{*}-g)n\delta\psi^{*}-g_{c}^{*}n\delta\phi,$

with $n=n_{0}(1+2\gamma n_{0}t)^{-1}$. Now if we consider perturbations with
momentum k and naively treat $n(t)$ as time invariant, we may find the eigen-
frequencies of the above equations are

$\omega_{\textbf{k}}=-2i\gamma
n\pm\sqrt{\epsilon_{\textbf{k}}^{2}+2gn\epsilon_{\textbf{k}}-\gamma^{2}n^{2}}$
(93)

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 equations-_ 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
Pethick and Smith (2008); Dalfovo  _et al._ (1999); Zhai (2021). It is then
interesting to construct the hydrodynamic equation with two-body loss.
Starting from time-dependent dissipative Gross-Pitaevskii equation eq. (91),
we can derive the hydrodynamic equations which govern dynamics of a
dissipative fluid. By decomposing $\psi=\sqrt{\rho}e^{i\theta}$, the real part
and imaginary part give the hydrodynamic equations:

$\displaystyle\partial_{t}\rho+\nabla(\rho\textbf{v})+2\gamma\rho^{2}=0$ (94)
$\displaystyle m\partial_{t}\textbf{v}=-\nabla($
$\displaystyle\frac{1}{2m}\frac{1}{\sqrt{\rho}}\nabla^{2}\sqrt{\rho}+\frac{1}{2}mv^{2}+V(\textbf{r})+g\rho).$

The first equation is called the continuity equation, which reflects the
conservation of particle number. Now two-body losses in interacting BEC bring
new term $\gamma\rho^{2}$ in the continuity equation, which break the particle
number conservation. And the second equation, the Newton equation, is same as
the conventional one. For a uniform system with $V(\mathbf{r})=0$, the uniform
solution is

$\rho(\mathbf{r},t)=\frac{\rho_{0}}{1+2\gamma t\rho_{0}},$ (95)

which show that condensate particles are always decay with time and the vaccum
is the only true steady state of this system. With the new dissipation term
$\gamma\rho^{2}$, finding the solutions for general case is challenging, we
will leave this for future research.

##  VII Summary and Outlooks

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 system
and verify the validity of this approximation by numerical simulating a toy
model which has a similar structure. Based on this time-dependent Bogoliubov
Lindbladian, we study the quench problem and prove this system has a
Sp$(4,\mathbb{C})$ dynamical symmetry, which is crucial for the exact
calculation of quench dynamics. Furthermore, we show a general $n$-mode
quadratic Lindbladian of bosonic system has a dynamical symmetry of
Sp$(2n,\mathbb{C})$, which is useful for the analytical understanding of
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 system, we
only discuss the solution without external potential in this paper, it would
be significant to generalize these solutions to open system and study the
stable property of these results. Finally, from 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 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.

##  VIII Acknowledgements

We gratefully acknowledge fruitful discussions with Hui Zhai and Xin Chen.
This work is supported by NSFC under Grant No. 12004115, Program of Shanghai
Sailing Program Grant No. 20YF1411600.

## References

  * Landau and Lifshitz (2013) L. D. Landau and E. M. Lifshitz,  _Quantum mechanics: non-relativistic theory_ , Vol. 3 (Elsevier, 2013).
  * Blatt and Weisskopf (1991) J. M. Blatt and V. F. Weisskopf,  _Theoretical nuclear physics_ (Courier Corporation, 1991).
  * Huang and Yang (1957) K. Huang and C. N. Yang,  _Quantum-mechanical many-body problem with hard-sphere interaction_ , Physical review 105, 767 (1957).
  * Inouye  _et al._ (1998) S. Inouye, M. Andrews, J. Stenger, H.-J. Miesner, D. M. Stamper-Kurn, and W. Ketterle,  _Observation of Feshbach resonances in a Bose–Einstein condensate_ , Nature 392, 151 (1998).
  * Altmeyer  _et al._ (2007) A. Altmeyer, S. Riedl, C. Kohstall, M. Wright, R. Geursen, M. Bartenstein, C. Chin, J. H. Denschlag, and R. Grimm,  _Precision measurements of collective oscillations in the BEC-BCS crossover_ , Physical review letters 98, 040401 (2007).
  * Shin  _et al._ (2008) Y.-i. Shin, A. Schirotzek, C. H. Schunck, and W. Ketterle,  _Realization of a strongly interacting Bose-Fermi mixture from a two-component Fermi gas_ , Physical Review Letters 101, 070404 (2008).
  * Papp  _et al._ (2008) S. Papp, J. Pino, R. Wild, S. Ronen, C. E. Wieman, D. S. Jin, and E. A. Cornell,  _Bragg spectroscopy of a strongly interacting rb 85 bose-einstein condensate_ , Physical review letters 101, 135301 (2008).
  * Navon  _et al._ (2010) N. Navon, S. Nascimbene, F. Chevy, and C. Salomon,  _The equation of state of a low-temperature Fermi gas with tunable interactions_ , Science 328, 729 (2010).
  * Navon  _et al._ (2011) N. Navon, S. Piatecki, K. Günter, B. Rem, T. C. Nguyen, F. Chevy, W. Krauth, and C. Salomon,  _Dynamics and thermodynamics of the low-temperature strongly interacting Bose gas_ , Physical review letters 107, 135301 (2011).
  * Donley  _et al._ (2001) E. A. Donley, N. R. Claussen, S. L. Cornish, J. L. Roberts, E. A. Cornell, and C. E. Wieman,  _Dynamics of collapsing and exploding Bose–Einstein condensates_ , Nature 412, 295 (2001).
  * Cornish  _et al._ (2006) S. L. Cornish, S. T. Thompson, and C. E. Wieman,  _Formation of bright matter-wave solitons during the collapse of attractive Bose-Einstein condensates_ , Physical review letters 96, 170401 (2006).
  * Strecker  _et al._ (2002) K. E. Strecker, G. B. Partridge, A. G. Truscott, and R. G. Hulet,  _Formation and propagation of matter-wave soliton trains_ , Nature 417, 150 (2002).
  * Ritsch  _et al._ (2013) H. Ritsch, P. Domokos, F. Brennecke, and T. Esslinger,  _Cold atoms in cavity-generated dynamical optical potentials_ , Reviews of Modern Physics 85, 553 (2013).
  * Kasprzak  _et al._ (2006) J. Kasprzak, M. Richard, S. Kundermann, A. Baas, P. Jeambrun, J. M. J. Keeling, F. Marchetti, M. Szymańska, R. André, J. Staehli, _et al._ ,  _Bose–Einstein condensation of exciton polaritons_ , Nature 443, 409 (2006).
  * Mabuchi and Doherty (2002) H. Mabuchi and A. Doherty,  _Cavity quantum electrodynamics: coherence in context_ , Science 298, 1372 (2002).
  * Öztop  _et al._ (2012) B. Öztop, M. Bordyuh, Ö. E. Müstecaplıoğlu, and H. E. Türeci,  _Excitations of optically driven atomic condensate in a cavity: theory of photodetection measurements_ , New Journal of Physics 14, 085011 (2012).
  * Deng  _et al._ (2010) H. Deng, H. Haug, and Y. Yamamoto,  _Exciton-polariton bose-einstein condensation_ , Reviews of modern physics 82, 1489 (2010).
  * García-Ripoll  _et al._ (2009) J. J. García-Ripoll, S. Dürr, N. Syassen, D. M. Bauer, M. Lettner, G. Rempe, and J. I. Cirac,  _Dissipation-induced hard-core boson gas in an optical lattice_ , New Journal of Physics 11, 013053 (2009).
  * Yamamoto  _et al._ (2019) K. Yamamoto, M. Nakagawa, K. Adachi, K. Takasan, M. Ueda, and N. Kawakami,  _Theory of non-Hermitian fermionic superfluidity with a complex-valued interaction_ , Physical review letters 123, 123601 (2019).
  * He  _et al._ (2020) M. He, C. Lv, H.-Q. Lin, and Q. Zhou,  _Universal relations for ultracold reactive molecules_ , Science Advances 6, eabd4699 (2020).
  * Bouchoule  _et al._ (2020) I. Bouchoule, B. Doyon, and J. Dubail,  _The effect of atom losses on the distribution of rapidities in the one-dimensional Bose gas_ , SciPost Physics 9, 044 (2020).
  * Yamamoto  _et al._ (2021) K. Yamamoto, M. Nakagawa, N. Tsuji, M. Ueda, and N. Kawakami,  _Collective excitations and nonequilibrium phase transition in dissipative fermionic superfluids_ , Physical Review Letters 127, 055301 (2021).
  * Bouchoule and Dubail (2021) I. Bouchoule and J. Dubail,  _Breakdown of Tan’s relation in lossy one-dimensional Bose gases_ , Physical Review Letters 126, 160603 (2021).
  * Rossini  _et al._ (2021) D. Rossini, A. Ghermaoui, M. B. Aguilera, R. Vatré, R. Bouganne, J. Beugnon, F. Gerbier, and L. Mazza,  _Strong correlations in lossy one-dimensional quantum gases: from the quantum Zeno effect to the generalized Gibbs ensemble_ , Physical Review A 103, L060201 (2021).
  * Rosso  _et al._ (2022) L. Rosso, A. Biella, and L. Mazza,  _The one-dimensional Bose gas with strong two-body losses: the effect of the harmonic confinement_ , SciPost Physics 12, 044 (2022).
  * Nakagawa  _et al._ (2021) M. Nakagawa, N. Kawakami, and M. Ueda,  _Exact liouvillian spectrum of a one-dimensional dissipative hubbard model_ , Physical Review Letters 126, 110404 (2021).
  * Bouchoule  _et al._ (2021) I. Bouchoule, L. Dubois, and L.-P. Barbier,  _Losses in interacting quantum gases: Ultraviolet divergence and its regularization_ , Physical Review A 104, L031304 (2021).
  * Gerbier and Castin (2010) F. Gerbier and Y. Castin,  _Heating rates for an atom in a far-detuned optical lattice_ , Physical Review A 82, 013615 (2010).
  * Bouganne  _et al._ (2020) R. Bouganne, M. Bosch Aguilera, A. Ghermaoui, J. Beugnon, and F. Gerbier,  _Anomalous decay of coherence in a dissipative many-body system_ , Nature Physics 16, 21 (2020).
  * Tindall  _et al._ (2019) J. Tindall, B. Buča, J. R. Coulthard, and D. Jaksch,  _Heating-induced long-range $\eta$ pairing in the Hubbard model_, Physical Review Letters 123, 030603 (2019).
  * Kantian  _et al._ (2009) A. Kantian, M. Dalmonte, S. Diehl, W. Hofstetter, P. Zoller, and A. Daley,  _Atomic color superfluid via three-body loss_ , Physical review letters 103, 240401 (2009).
  * Syassen  _et al._ (2008) N. Syassen, D. M. Bauer, M. Lettner, T. Volz, D. Dietze, J. J. Garcia-Ripoll, J. I. Cirac, G. Rempe, and S. Durr,  _Strong dissipation inhibits losses and induces correlations in cold molecular gases_ , Science 320, 1329 (2008).
  * Tomita  _et al._ (2017) T. Tomita, S. Nakajima, I. Danshita, Y. Takasu, and Y. Takahashi,  _Observation of the Mott insulator to superfluid crossover of a driven-dissipative Bose-Hubbard system_ , Science advances 3, e1701513 (2017).
  * Sponselee  _et al._ (2018) K. Sponselee, L. Freystatzky, B. Abeln, M. Diem, B. Hundt, A. Kochanke, T. Ponath, B. Santra, L. Mathey, K. Sengstock, _et al._ ,  _Dynamics of ultracold quantum gases in the dissipative Fermi–Hubbard model_ , Quantum Science and Technology 4, 014002 (2018).
  * Wang  _et al._ (2021) C. Wang, C. Liu, and Z.-Y. Shi,  _Complex contact interaction for systems with short-range two-body loss_ , arXiv preprint arXiv:2109.11284 (2021).
  * Bogoliubov (1947) N. Bogoliubov,  _On the theory of superfluidity_ , J. Phys 11, 23 (1947).
  * Bohn and Julienne (1996) J. L. Bohn and P. S. Julienne,  _Semianalytic treatment of two-color photoassociation spectroscopy and control of cold atoms_ , Physical Review A 54, R4637 (1996).
  * Bohn and Julienne (1999) J. L. Bohn and P. S. Julienne,  _Semianalytic theory of laser-assisted resonant cold collisions_ , Physical Review A 60, 414 (1999).
  * Chin  _et al._ (2010) C. Chin, R. Grimm, P. Julienne, and E. Tiesinga,  _Feshbach resonances in ultracold gases_ , Rev. Mod. Phys. 82, 1225 (2010).
  * Thalhammer  _et al._ (2005) G. Thalhammer, M. Theis, K. Winkler, R. Grimm, and J. H. Denschlag,  _Inducing an optical Feshbach resonance via stimulated Raman coupling_ , Physical Review A 71, 033403 (2005).
  * Breuer  _et al._ (2002) H.-P. Breuer, F. Petruccione, _et al._ ,  _The theory of open quantum systems_ (Oxford University Press on Demand, 2002).
  * (42) This term could be regarded as $Gamma(I)$ because the square of the renormalized coupling strength $|\alpha_{\text{re}}|^{2}$ is also proportional to the laser intensity.
  * Chen  _et al._ (2020) Y.-Y. Chen, P. Zhang, W. Zheng, Z. Wu, and H. Zhai,  _Many-body echo_ , Physical Review A 102, 011301 (2020).
  * Cheng and Shi (2021) Y. Cheng and Z.-Y. Shi,  _Many-body dynamics with time-dependent interaction_ , Physical Review A 104, 023307 (2021).
  * Lv  _et al._ (2020) C. Lv, R. Zhang, and Q. Zhou,  _S U (1, 1) Echoes for Breathers in Quantum Gases_ , Physical Review Letters 125, 253002 (2020).
  * Lee  _et al._ (1957) T. D. Lee, K. Huang, and C. N. Yang,  _Eigenvalues and eigenfunctions of a Bose system of hard spheres and its low-temperature properties_ , Physical Review 106, 1135 (1957).
  * Lee and Yang (1957) T. Lee and C. Yang,  _Many-body problem in quantum mechanics and quantum statistical mechanics_ , Physical Review 105, 1119 (1957).
  * Prosen (2008) T. Prosen,  _Third quantization: a general method to solve master equations for quadratic open Fermi systems_ , New Journal of Physics 10, 043026 (2008).
  * Prosen and Seligman (2010) T. Prosen and T. H. Seligman,  _Quantization over boson operator spaces_ , Journal of Physics A: Mathematical and Theoretical 43, 392004 (2010).
  * Wei and Norman (1963) J. Wei and E. Norman,  _Lie algebraic solution of linear differential equations_ , Journal of Mathematical Physics 4, 575 (1963).
  * Scopa  _et al._ (2019) S. Scopa, G. T. Landi, A. Hammoumi, and D. Karevski,  _Exact solution of time-dependent Lindblad equations with closed algebras_ , Physical Review A 99, 022105 (2019).
  * (52) In fact, the structure of the master equation can further restrict these 10 elements into seven:$\tilde{\mathcal{H}}_{\mathbf{k}}^{0}=\tilde{h}^{1}_{\mathbf{k}}-\tilde{h}^{4}_{\mathbf{k}},\tilde{\mathcal{H}}_{\mathbf{k}}^{1}=\tilde{h}^{2}_{\mathbf{k}}-\tilde{h}^{6}_{\mathbf{k}},\tilde{\mathcal{H}}_{\mathbf{k}}^{2}=\tilde{h}^{3}_{\mathbf{k}}-\tilde{h}^{5}_{\mathbf{k}},\tilde{\mathcal{L}}_{\mathbf{k}}^{0}=\tilde{h}^{1}_{\mathbf{k}}+\tilde{h}^{4}_{\mathbf{k}}-2\tilde{h}^{7}_{\mathbf{k}},\tilde{\mathcal{L}}_{\mathbf{k}}^{1}=\tilde{h}^{2}_{\mathbf{k}}+\tilde{h}^{6}_{\mathbf{k}}-2\tilde{h}^{8}_{\mathbf{k}},\tilde{\mathcal{L}}_{\mathbf{k}}^{2}=\tilde{h}^{3}_{\mathbf{k}}+\tilde{h}^{5}_{\mathbf{k}}-2\tilde{h}^{9}_{\mathbf{k}},\tilde{\mathcal{L}}_{\mathbf{k}}^{3}=\tilde{h}^{1}_{\mathbf{k}}+\tilde{h}^{4}_{\mathbf{k}}-2\tilde{h}^{10}_{\mathbf{k}}$.
  * Hall (2013) B. C. Hall,  _Lie groups, Lie algebras, and representations_ , in  _Quantum Theory for Mathematicians_ (Springer, 2013) pp. 333–366.
  * Dalfovo  _et al._ (1999) F. Dalfovo, S. Giorgini, L. P. Pitaevskii, and S. Stringari,  _Theory of Bose-Einstein condensation in trapped gases_ , Reviews of modern physics 71, 463 (1999).
  * Pethick and Smith (2008) C. J. Pethick and H. Smith,  _Bose–Einstein condensation in dilute gases_ (Cambridge university press, 2008).
  * Zhai (2021) H. Zhai,  _Ultracold Atomic Physics_ (Cambridge University Press, 2021).
  * Sieberer  _et al._ (2016) L. M. Sieberer, M. Buchhold, and S. Diehl,  _Keldysh field theory for driven open quantum systems_ , Reports on Progress in Physics 79, 096001 (2016).
  * Kamenev (2011) A. Kamenev,  _Field theory of non-equilibrium systems_ (Cambridge University Press, 2011).