text:
\section{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~\cite{landau2013quantum}, 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~\cite{blatt1991theoretical,huang1957quantum}. More importantly, experimental techniques such as Feshbach resonance can adjust $a$ along the real-axis in cold atomic gases~\cite{inouye1998observation}, which allows us to demonstrate various many-body effects, from the Lee-Huang-Yang correction~\cite{altmeyer2007precision,shin2008realization,papp2008bragg,navon2010equation,navon2011dynamics} for positive scattering length to the Bose  nova effect~\cite{donley2001dynamics,cornish2006formation,strecker2002formation} for negative scattering length.

 Recently, with the development of theoretical and experimental methods~\cite{ritsch2013cold,kasprzak2006bose,mabuchi2002cavity,oztop2012excitations,deng2010exciton,garcia2009dissipation,yamamoto2019theory,he2020universal,bouchoule2020effect,yamamoto2021collective,bouchoule2021breakdown,rossini2021strong,rosso2022one,nakagawa2021exact,bouchoule2021losses}, 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~\cite{gerbier2010heating,bouganne2020anomalous,tindall2019heating,kantian2009atomic,syassen2008strong,tomita2017observation,sponselee2018dynamics}. 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~\cite{wang2021complex}, 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~\cite{bogoliubov1947theory}. 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.~\ref{sec2}, 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.~\ref{sec3} 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.~\ref{sec4}, and discuss the dynamical symmetry of the Bogoliubov Lindbladian in Section.~\ref{sec5}. In Section.~\ref{sec6}, we derive the dissipative version of the Gross-Pitaevskii equation and hydrodynamic theory using the Keldysh path integral formalism for open systems.

\section{Tuning complex scattering lengths in experiments}
\label{sec2}
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~\cite{landau2013quantum,bohn1996semianalytic,bohn1999semianalytic,chin2010feshbach}. 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~\cite{thalhammer2005inducing}. As shown in Fig.~\ref{sm} (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
\begin{align}
    a_c(E)=a_{\text{bg}}\left(1+\frac{\Gamma(I)}{E-\nu-\Gamma(I)+i\gamma/2}\right).\label{a_c_general}
\end{align}

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)
\begin{equation}
\begin{aligned}
\hat{H}_{\text{ofr}} =
&\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'},\mathbf{p}}\hat{b}^{\dagger}_{\mathbf{k}+\mathbf{p}}\hat{b}^{\dagger}_{\mathbf{k'}-\mathbf{p}}\hat{b}_{\mathbf{k'}}\hat{b}_{\mathbf{k}}\\  + &\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),
\nonumber
\end{aligned}
\end{equation}
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$.

\begin{figure}[ht]
	\includegraphics[width=1\linewidth]{sm}
	\caption{(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)$.} \label{sm}
\end{figure}
The spontaneous radiation of the bound-state molecule can be characterized by the following Lindblad master equation~\cite{breuer2002theory}
\begin{equation}
\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}.
\nonumber
\end{equation}
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
\begin{equation}
\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}}.
\end{equation}
We can obtain the two-body scattering matrix $T_{2}$ for this $\hat{H}_\text{eff}$ at relative kinetic energy $E$
\begin{equation}
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}. \nonumber
\end{equation}
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}$~\cite{landau2013quantum}, we obtain the renormalize relation for the complex scattering length $a_c$,
\begin{equation}
\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},
\end{equation}
where $\Lambda$ is the momentum cut-off and $a_{c}(E)$ can be further written as
\begin{equation}
a_c(E)=a_\text{bg}+\frac{m}{4\pi}\frac{|\alpha_{\text{re}}|^2}{E-\nu_{\text{re}}+i\gamma/2}
\end{equation}\label{a_c_derived}
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.~\eqref{a_c_general} if we take the limit $\Lambda \to \infty$, and write $\Gamma(I)=\frac{m|\alpha_{\text{re}}|^2}{4\pi a_{\text{bg}}}$~\cite{gamma}.

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.~\ref{sm}(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.

\section{dissipative Bose  gas}
\label{sec3}

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~\cite{chin2010feshbach}. The Hermitian Hamiltonian of a Bose  gas may be written as
\begin{align}
		\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}',\mathbf{p}}\hat{a}^\dagger_{\mathbf{k+p}}\hat{a}^\dagger_{\mathbf{k}'-\mathbf{p}}\hat{a}_{\mathbf{k}'}\hat{a}_\mathbf{k},
\end{align}
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~\cite{wang2021complex}, 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~\cite{breuer2002theory,wang2021complex}
 \begin{equation}
    \partial_{t}\hat{\rho} =-i( \hat{H}_\text{eff}\hat{\rho}-\hat{\rho}\hat{H}_\text{eff}^\dagger) + \mathcal{J}\hat{\rho},  \label{lind1}
\end{equation}
with the non-Hermitian effective Hamiltonian
\begin{align}
		\hat{H}_\text{eff}=H-\frac{i\gamma_b}{2\Omega}\sum_{\mathbf{k},\mathbf{k}',\mathbf{p}}\hat{a}^\dagger_{\mathbf{k+p}}\hat{a}^\dagger_{\mathbf{k}'-\mathbf{p}}\hat{a}_{\mathbf{k}'}\hat{a}_\mathbf{k},
	\end{align}
	and the recycling term
	\begin{align}
		\mathcal{J}\hat{\rho}=\frac{\gamma_b}{\Omega}\sum_{\mathbf{k},\mathbf{k}',\mathbf{p}}\hat{a}_{\mathbf{k}'}\hat{a}_\mathbf{k}\hat{\rho}\hat{a}^\dagger_{\mathbf{k'-p}}\hat{a}^\dagger_{\mathbf{k}+\mathbf{p}},
	\end{align}
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~\cite{wang2021complex}
\begin{equation}
    \frac{1}{g_b-i\gamma_b} = \frac{1}{g-i\gamma}-\frac{1}{\Omega}\sum_{\textbf{k}}\frac{1}{2\epsilon_{\textbf{k}}},
    \label{renormalization}
\end{equation}
and define $g-i\gamma \equiv \frac{4\pi a_c}{m}$ as renormalized complex coupling constant.

\emph{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}\ne \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}'$ as $\partial_t \hat{\rho}' \approx \mathcal{L}_B(\hat{\rho}' )$ with
\begin{equation}
	 \mathcal{L}_B(\hat{\rho}' ) = -i [\hat{H}_B,\hat{\rho}' ]-2\gamma n \sum_{\textbf{k}\neq0} \{ \hat{a}_{\textbf{k}}^\dagger \hat{a}_{\textbf{k}},\hat{\rho}' \}+4\gamma n \sum_{\textbf{k}\neq0}\hat{a}_{\textbf{k}}\hat{\rho}' \hat{a}_{\textbf{k}}^\dagger  \label{ME}
\end{equation}
and
\begin{equation}
\hat{H}_B = \frac{gnN}{2}+\sum_{\textbf{k}\neq0}((\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.),\nonumber
\end{equation}
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~\cite{bogoliubov1947theory}. 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~\cite{chen2020many,cheng2021many,lv2020s},
\begin{equation}
[\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}}. \label{su11}
\end{equation}
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
\begin{align}
	\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).\label{matrix_eq}
\end{align}

Here, we focused on the short-time dynamics where $\gamma nt\ll1$ such that we may approximate $n$ as a constant at a fixed time t. We can diagonalize the $3 \times 3$ matrix in eq.~\eqref{matrix_eq} to obtain the eigenvalue $-2i\xi_{\textbf{k},i}$, where the quasi-steady state eigenvalue is
\begin{equation}
    \xi_{\mathbf{k},0} = -2i\gamma n,
    \quad \xi_{\mathbf{k},(1,2)}=-2i\gamma n  \pm \sqrt{\epsilon_\mathbf{k}^2+2g_0n\epsilon_\mathbf{k}-\gamma_0^2n^2}.
    \nonumber
\end{equation}
And a phase boundary between stable and unstable BEC can be obtained by this quasi-steady state eigenvalue~\cite{wang2021complex}.

\emph{Lee-Huang-Yang correction-} Furthermore, we will give a detailed derivation for the Lee-Huang-Yang correction~\cite{lee1957eigenvalues,lee1957many} for this open system. Consider quasi-steady-state eigenvalue, up to $O\left(\gamma nt\right)$, simple solution of eq.~\eqref{matrix_eq} can be obtained
\begin{align}
	\mathbf{A}(t)=\mathbf{A}_\text{s}+e^{-4\gamma nt}\left(\mathbf{B}+\mathbf{C}\cos 2\xi_\mathbf{k}t+\mathbf{D}\sin2\xi_\mathbf{k}t\right)\label{asy_sol
\end{align}
constant vectors $\mathbf{B,C,D}$ can be determined by the initial value of $\mathbf{A}$, and
\begin{equation}
\begin{aligned}
	\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^2n^2}
	\\-\frac{gn\epsilon_\mathbf{k}+g^2n^2+2\gamma^2n^2}{2\epsilon_\mathbf{k}^2+4gn\epsilon_\mathbf{k}+6\gamma^2n^2}
	\\-\frac{\gamma n(\epsilon_\mathbf{k}-gn)}{2\epsilon_\mathbf{k}^2+4gn\epsilon_\mathbf{k}+6\gamma^2n^2}
	\end{array}\right).
	\label{steady_value}
\end{aligned}
\end{equation}
We see that $\mathbf{A}$ never reaches the steady value $\mathbf{A}_\text{s}$ because the asymptotic solution eq.~\eqref{asy_sol} 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.~\eqref{renormalization}.

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
\begin{align}
	\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_0a_0\right)=\text{tr}\left(\hat{\rho}\mathcal{L}'(a_0^\dagger a_0)\right),\nonumber
\end{align}
with $\mathcal{L}'$ defined by
\begin{equation}
\begin{aligned}
	\mathcal{L}'(\hat{O})\equiv i\left[\hat{H},\hat{O}\right]-\frac{\gamma_b}{2V}\sum_{\mathbf{k,k',p}}\left\{a^\dagger_{\mathbf{k+p}}a^\dagger_{\mathbf{k'-p}}a_\mathbf{k'}a_\mathbf{k},\hat{O}\right\}
	\\+\frac{\gamma_b}{2V}\sum_{\mathbf{k,k',p}}a^\dagger_{\mathbf{k+p}}a^\dagger_{\mathbf{k'-p}}\hat{O}a_\mathbf{k'}a_\mathbf{k}.\nonumber
\end{aligned}
\end{equation}
It is then straightforward to show that
\begin{align}
	\mathcal{L}'(a^\dagger_0 a_0)\simeq -2\gamma_b nN-2\sum_{\mathbf{k}\neq0}\left(\gamma_b nA_1^\mathbf{k}+g_b nA_2^\mathbf{k}\right).
\end{align}
This leads to
\begin{align}
	\frac{dN_0}{dt}=-2\gamma_b nN-2\sum_{\mathbf{k}\neq0}\left(\gamma_b n \langle A_1^{\mathbf{k}}\rangle+g_bn\langle A_2^{\mathbf{k}}\rangle\right).\label{summation}
\end{align}
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
\begin{align}
	\frac{dN_0}{dt}=-2\gamma_b nN+2\gamma n\sum_\mathbf{k}\frac{gn\epsilon_\mathbf{k}+\gamma^2n^2}{\epsilon_\mathbf{k}^2+2gn\epsilon_\mathbf{k}+3\gamma^2n^2}.
\end{align}
We immediately find that the summation diverges for large momentum. To cure this divergence, we use again the renormalization relation eq.~\eqref{renormalization}, which gives
\begin{equation}
\begin{aligned}
	g_b-i\gamma_b&=\frac{1}{\frac{1}{g-i\gamma}-\frac{1}{V}\sum_{\mathbf{k}}\frac{1}{2\epsilon_\mathbf{k}}}\\
	&=(g-i\gamma)\left(1+\frac{g-i\gamma}{V}\sum_\mathbf{k}\frac{1}{2\epsilon_\mathbf{k}}+\ldots\right).
\end{aligned}
\end{equation}
We thus obtain
\begin{align}
	\gamma_b\simeq\gamma+\frac{g\gamma}{V}\sum_\mathbf{k}\frac{1}{\epsilon_\mathbf{k}}.\label{renormalization_2nd}
\end{align}
Using this expansion, we have
\begin{equation}
\begin{aligned}
	&\frac{dN}{dt}=\frac{dN_0}{dt}\\
	&=-2\gamma nN+2\gamma n\sum_\mathbf{k}\left(\frac{gn\epsilon_\mathbf{k}+\gamma^2n^2}{\epsilon_\mathbf{k}^2+2gn\epsilon_\mathbf{k}+3\gamma^2n^2}-\frac{g_0n}{\epsilon_\mathbf{k}}\right)\\
	&=-2\gamma nN\left[1+c_\theta(n|a_c|^3)^{1/2}\right],\label{lhy_losses}
\end{aligned}
\end{equation}
with $c_\theta$ a constant that only depends on the argument of the complex scattering length $a_c$,
\begin{align}
	c_\theta=4\sqrt{2\pi}\left(\frac{\cos2\theta}{\sqrt{\cos(\theta-\pi/3)}}+2\cos\theta\sqrt{\cos(\theta-\pi/3)}\right).\nonumber
\end{align}
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.~\eqref{lhy_losses} is an analog of the celebrated Lee-Huang-Yang correction to the ground state energy of a weakly interacting Bose  gas~\cite{lee1957many,lee1957eigenvalues},
\begin{align}
	E_0=\frac{g_0nN}{2}\left[1+\frac{128}{15\pi^{1/2}}(na_r^3)^{1/2}+O(n\log n)\right].
\end{align}

With the help of renormalization relation, we can eliminate the divergence eq.~\eqref{summation} and obtain a physical result of particle losses rate.

\section{A Toy model}
\label{sec4}
 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.
\begin{figure}
\begin{center}
\includegraphics[width=0.5\textwidth]{phase_fig.pdf}
\end{center}
\caption{(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{g n}{\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{g n}{\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.
}
\label{toymodel}
\end{figure}

\emph{Numerical model-}  As shown in Fig.~\ref{toymodel}, 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
\begin{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),
\end{equation}
where the Hamiltonian is
\begin{equation}
    \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}),
\end{equation}
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
\begin{equation}
   \partial_{t}\hat{\rho} = -i(\hat{H}_a^{\text{eff}}\hat{\rho}-\hat{\rho}\hat{H}_a^{\text{eff}\dagger})+\mathcal{J}\hat{\rho},
    \label{Ldw}
\end{equation}
where the effective Hamiltonian $\hat{H}_a^{\text{eff}}$ for $\hat{a}_{0} ,\hat{a}_{1} $is
\begin{equation}
\begin{aligned}
    \hat{H}_a^{\text{eff}} = &\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\\
    &+\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),
    \end{aligned}
\end{equation}
and  the recycling term in eq.~\eqref{Ldw} are similar to Lindbladian of Bose  gas with two-body losses,
\begin{equation}
\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).\nonumber
\end{equation}

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}$~\cite{wang2021complex}
\begin{equation}
\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},
\end{equation}
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.~\eqref{Ldw} 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
\begin{equation}
\begin{aligned}
 \hat{\mathcal{L}}_{\text{dw}}^B\hat{\rho} &= -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} \}\\
 &+\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),
 \end{aligned}\label{LB2}
 \end{equation}
where the Hamiltonian $\hat{H}_{\text{dw}}^B$ is
\begin{equation}
\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),
\end{equation}
and time-dependent particle number is given by $n = N(1+2\gamma N t)^{-1}$. Similar to eq.~\eqref{matrix_eq} 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,
\begin{align}
	\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),\label{matrix_eqdw}
\end{align}
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.~\eqref{matrix_eqdw} are
\begin{equation}
    \xi_0 = -2i\gamma n,\quad \xi_{(1,2)}= -2i\gamma n\pm \sqrt{\epsilon^2+2g n \epsilon-\gamma^2 n^2}.
\end{equation}
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{g n}{\epsilon} = -\frac{1+3 (\frac{\gamma n}{\epsilon})^2}{2}$ as shown in Fig.~\ref{toymodel}(b).

\begin{figure}
\begin{center}
\includegraphics[width=0.48\textwidth]{stable_unstable}
\end{center}
\caption{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.~\eqref{Ldw}; dashed: approximated Lindbladian eq.~\eqref{LB2}). 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$.}
\label{results}
\end{figure}

 \emph{Numerical results -}
 We now check the accuracy of the Bogoliubov approximation by comparing the numerical simulation result from eq.~\eqref{Ldw} and eq.~\eqref{LB2}. The evolution begins with an initial state with $N=100$ particles occupying the mode $\hat{a}_0$. As shown in Fig.~\ref{results}, 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.~\eqref{LB2} deviates from the exact time evolution governed by eq.~\eqref{Ldw} 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.

 \section{Dynamical symmetry of Lindbladian}
 \label{sec5}

In this section, we discuss the dynamical symmetry of quadratic Lindbladians in detail, which can greatly simplify the calculation of the master eq.~\eqref{ME}. It is worth mentioning that a time-independent quadratic Lindbladian may be solved explicitly using the third quantization method~\cite{prosen2008third,prosen2010quantization}. However, this cannot be directly applied to our problem because of the time-dependent $n(t)$ contained in the master eq.~\eqref{ME}. 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~\cite{wei1963lie,scopa2019exact} to construct exact solutions to the time-dependent Lindbladian algebraically.

  \emph{Symmetry-} The Lindbladian $\mathcal{L}_B$ in eq.~\eqref{ME} 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.~\eqref{ME}. 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.

 \emph{Closed algebra-} To find this algebra, we note that $\mathcal{L}_B$ can be expressed in the superoperator-formalism as
\begin{equation}
\mathcal{L}_B(\hat{\rho}' ) = \left( \sum_{\mathbf{k},k_{z}\ge 0}\tilde{\mathcal{L}}_{\mathbf{k}}\right)\circ \hat{\rho}'.
\end{equation}

Here $\tilde{\mathcal{L}}_{\mathbf{k}}$ is a superoperator act only on modes $\textbf{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,\ldots7$. They are defined by
\begin{align}
&\tilde{h}^{1}_{\mathbf{k}} \circ \hat{\rho}'= (\hat{a}_{\mathbf{k}}^{\dagger}\hat{a}_{\mathbf{k}}+\hat{a}_{-\mathbf{k}}\hat{a}_{-\mathbf{k}}^{\dagger}) \hat{\rho}',\quad \tilde{h}^{2}_{\mathbf{k}} \circ \hat{\rho}'= 2\hat{a}_{\mathbf{k}}^{\dagger}\hat{a}_{-\mathbf{k}}^{\dagger} \hat{\rho}',&\nonumber \\
&\tilde{h}^{3}_{\mathbf{k}} \circ \hat{\rho}'= 2\hat{a}_{\mathbf{k}}\hat{a}_{-\mathbf{k}} \hat{\rho}',\quad \tilde{h}^{4}_{\mathbf{k}} \circ \hat{\rho}'= \hat{\rho}'(\hat{a}_{\mathbf{k}}^{\dagger}\hat{a}_{\mathbf{k}}+\hat{a}_{-\mathbf{k}}\hat{a}_{-\mathbf{k}}^{\dagger}),&\nonumber\\
&\tilde{h}^{5}_{\mathbf{k}} \circ \hat{\rho}'= 2\hat{\rho}'\hat{a}_{\mathbf{k}}^\dagger\hat{a}_{-\mathbf{k}}^\dagger ,\quad \tilde{h}^{6}_{\mathbf{k}} \circ \hat{\rho}'= 2\hat{\rho}'\hat{a}_{\mathbf{k}}\hat{a}_{\mathbf{k}},\nonumber\\
&\tilde{h}^{7}_{\mathbf{k}} \circ \hat{\rho}'= \hat{a}_{\mathbf{k}}\hat{\rho}'\hat{a}_{\mathbf{k}}^\dagger+\hat{a}_{-\mathbf{k}}\hat{\rho}'\hat{a}_{-\mathbf{k}}^\dagger.&\nonumber
\end{align}
Among these seven superoperators, $\tilde{h}_\mathbf{k}^i,\ i=1,2,\ldots6$ 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
 \begin{equation}
 [\tilde{h}^{5}_{\mathbf{k}} ,\tilde{h}^{7}_{\mathbf{k}}]\circ \hat{\rho}' = -2(\hat{a}_{\mathbf{k}}\hat{\rho}' \hat{a}_{-\mathbf{k}}^\dagger+\hat{a}_{-\mathbf{k}} \hat{\rho}' \hat{a}_{\mathbf{k}}^\dagger).
 \end{equation}
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
 \begin{equation}
 \begin{aligned}
  &\tilde{h}^{8}_{\mathbf{k}} \circ \hat{\rho}'= \hat{a}_{-\mathbf{k}}^\dagger \hat{\rho}'\hat{a}_{\mathbf{k}}^\dagger+\hat{a}_{\mathbf{k}}^\dagger \hat{\rho}'\hat{a}_{-\mathbf{k}}^\dagger,\\
    &\tilde{h}^{9}_{\mathbf{k}} \circ \hat{\rho}'= \hat{a}_{\mathbf{k}}\hat{\rho}' \hat{a}_{-\mathbf{k}}^\dagger+\hat{a}_{-\mathbf{k}} \hat{\rho}' \hat{a}_{\mathbf{k}}^\dagger,\\
    &\tilde{h}^{10}_{\mathbf{k}} \circ \hat{\rho}'= \hat{a}_{\mathbf{k}}^\dagger\hat{\rho}' \hat{a}_{\mathbf{k}}+\hat{a}_{-\mathbf{k}}^\dagger \hat{\rho}' \hat{a}_{\mathbf{k}}.
  \end{aligned}\nonumber
 \end{equation}

 All of the superoperators $\tilde{h}_\mathbf{k}^i,\ i=1,2,\ldots,10$ now form a closed algebra~\cite{seven_element}, whose commutation relations are listed in Table.~\ref{tab1}. Hence these superoperators are generators of the dynamical symmetry group of Lindbladian. In the following of this section, because only particles at momentum $\textbf{k}$ and $-\textbf{k}$ are coupled, we will omit the label $\textbf{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$.

 \begin{table}[b]
\centering
\begin{tabular}{ l | llllllllll}
\hline
 & $\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}]$ \\ \hline
 $[\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 \\ \hline
\end{tabular}
\caption{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.}
 \label{tab1}
\end{table}

 \emph{Exact solution-}  Based on this closed algebra structure, the dynamical problem of time-dependent  Lindbladian can be solved exactly~\cite{wei1963lie,scopa2019exact}. 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
 \begin{equation}
 \hat{\rho}(t) = \tilde{\Lambda}(t,t_0) \circ \hat{\rho}(t_0).
 \end{equation}
 $\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,
 \begin{equation}
  \tilde{\Lambda}(t,t_0) =  \prod \limits_{i=1}\limits^{10}e^{g_i(t)\tilde{h}^i}.
 \end{equation}
 Substitute this equation back to the master equation, we can get
  \begin{equation}
  \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  \label{ge}\limits_{j=1}\limits^{10}e^{g_j(t)\tilde{h}^j},
 \end{equation}
 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.~\eqref{ge} and write the explicit expression of time derivative, we can obtain
 \begin{equation}
 \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}. \label{sol_evo}
  \end{equation}
  Using Baker-Campbell-Hausdorf formula~\cite{hall2013lie} and commutation relation in Table.~\ref{tab1}, we can formally write L.H.S of eq.~\eqref{sol_evo} as:
  \begin{equation}
  \sum_{mn} \partial_t g_m(t) \eta_{mn} \tilde{h}^n = \sum_{l=1}^{7}u_{l}(t)\tilde{h}^{l},
  \end{equation}
  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})$.

  \emph{Map to harmonic oscillators-}
  Even though the commutation relations between $\tilde{h}_\mathbf{k}^i$ seem complicated as shown in Table.~\ref{tab1}, 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.~\ref{tb2}, 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})$.

 \begin{table}[]
\centering
\begin{tabular}{|c|c|c|}
\hline
 &Harmonic Oscillator   &Bogoliubov Lindbladian\\  \hline
 $h^1$&  $( a^\dagger  a+ a  a^\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$&  $  a  a$
 &$(2\hat a_\textbf{k} \hat a_{-\textbf{k}})\circ$\\
  $h^4$&  $( b^\dagger  b+ b  b^\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$&  $ b  b$  &$\circ (2\hat a_\textbf{k} \hat a_{-\textbf{k}})$\\
 $h^7$&  $ a  b$
 &$\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$&  $ a  b^\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}}$\\ \hline
\end{tabular}
\caption{Operators in coupled harmonic  oscillator system and superoperators in dissipative Bose  gas system after Bogoliubov approximation with momentum $\textbf{k}$. These 10 operators and superoperators have the same commutation relation in Table.~\ref{tab1} and they can form the algebra $sp(4,\mathbb{C})$.}
\label{tb2}
\end{table}

 Sp\emph{$(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
\begin{align}
	\hat{M
	&=\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).\label{quadratic_1}\nonumber
\end{align}

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
\begin{align}
	\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}
						{}&I\\
						-I&{}
			\end{array}\right)\equiv M\Omega
\end{align}
with $M$ satisfying
\begin{align}
	\Omega M+M^\text{T}\Omega=0,\label{Lie_alg_prop}
\end{align}
which is the generator for Lie group Sp$(2n,\mathbb{C})$.

We now define new operators
\begin{align}
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.\\
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.
\end{align}
which are related by $b_i=\Omega_{ij}b^j$ and $b^i=\Omega^{ij}b_j$ with matrices
\begin{align}
	\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.\\
	\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.
\end{align}

We thus have following useful commutation relations
\begin{align}
	&[b_i,b^j]=\delta_i^j,\qquad\ \  [b^i,b_j]=-\delta^i_j,\\
	&[b^i,b^j]=\Omega^{ij},\qquad [b_i,b_j]=-\Omega_{ij}.
\end{align}

With all these definitions, the generic quadratic form in Eq.~\eqref{quadratic_1} can then be expressed as
\begin{align}
	\hat{M}=\frac{1}{2}b^iM_i^jb_j
\end{align}
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
\begin{align}
	[\hat{M},\hat{N}]=\frac{1}{2}b^i[M,N]_i^jb_j,\label{Lie_bracket}
\end{align}
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.

 \section{DISSIPATIVE Gross-Pitaevskii EQUATION and HYDRODYNAMIC THEORY}
  \label{sec6}
 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~\cite{dalfovo1999theory,pethick2008bose,zhai2021ultracold}. 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.

  \emph{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,
\begin{equation}
\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, \label{master3}
\end{equation}
where $\hat\psi(\textbf{r})$ is the bosonic annihilation operator at position $\textbf{r}$, and
\begin{align}
    \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})
\end{align}
is Hermitian Hamiltonian of interacting bosons.

The recycling term $\mathcal{J}\hat{\rho}$ is given by
\begin{equation}
\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}).
\end{equation}
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.~\eqref{master3},
\begin{equation}
Z = \int \mathcal{D}[\psi_+,\bar\psi_+,\psi_-,\bar\psi_-]e^{iS[\psi]},
\end{equation}
where we omit the initial condition and the Lindblad-Keldysh action~\cite{sieberer2016keldysh} is
\begin{equation}
\begin{aligned}
S &= \int dt d\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})))\\
&+\frac{i}{2}\gamma\int dt d\textbf{r} \sum_{\eta=\pm} \bar\psi_\eta(\textbf{r}) \bar\psi_\eta(\textbf{r}) \psi_\eta(\textbf{r}) \psi_\eta(\textbf{r})\\
&-i\gamma \int dt d\textbf{r} \quad \psi_+(\textbf{r}) \psi_+(\textbf{r})  \bar\psi_-(\textbf{r})  \bar\psi_-(\textbf{r}),
\end{aligned}
\end{equation}
where time runs from $-\infty$ to +$\infty$, and $s_\eta = 0$ for $\eta = +$, $s_\eta = 1$ for $\eta = -$.

\emph{Saddle-point approximation-}
We consider case that almost all particles occupies the ground state energy level, which means $N\approx N_0 \gg1$. 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,
\begin{equation}
\frac{\delta S}{\delta  \bar\psi_{\pm}} = 0,  \frac{\delta S}{\delta  \psi_{\pm}} = 0,
\end{equation}

then we obtain
\begin{equation}
 \begin{aligned}
&i\partial_t \psi_+ -(-\frac{\nabla^2}{2m}\psi_++g_c \bar\psi_+ \psi_+ \psi_+)=0\\
-&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\\
-&i\partial_t \psi_- +(-\frac{\nabla^2}{2m}\psi_-+g_c^* \bar\psi_- \psi_-\psi_-)-2i\gamma \psi_+\psi_+\bar\psi_-=0\\
&i\partial_t \bar\psi_-+(-\frac{\nabla^2}{2m}\bar\psi_-+g_c^* \bar\psi_-\psi_-\bar\psi_-)=0,
\label{saddle}
\end{aligned}
\end{equation}
where the complex interacting strength $g_c = g-i\gamma$.

 \emph{Non-Hermitian Gross-Pitaevskii equation -}
 The regularization of Keldysh action requires the relation $\psi_+=\psi_-,\bar\psi_+=\bar\psi_-$~\cite{kamenev2011field}, combining with the saddle-point equations eq.~\eqref{saddle}, we will obtain the Gross-Pitaevskii equation under the two-body losses $(\psi_+ = \psi)$,
\begin{equation}
i\partial_t\psi-\left(-\frac{\nabla^2}{2m}\psi+g_c|\psi|^2\psi+V\psi\right)=0 .\label{nGP}
\end{equation}
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
\begin{equation}
\psi_0(\textbf{r}) = \sqrt{n_0} e^{-i \frac{g_c}{2\gamma}\text{log}(1+2\gamma n_0t)}. \label{sol1}
\end{equation}
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^2t/2)}e^{-i \frac{g_c}{2\gamma}\text{log}(1+2\gamma n_0t)}$.

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$:
\begin{equation}
\begin{aligned}
&i\partial_t \delta \psi = -\frac{1}{2} \nabla^2 \delta\psi+2 g_c |\psi_0|^2\delta\psi+(g-i\gamma)\psi_0^2
\delta\psi^*,\\
&i\partial_t \delta \psi^* = \frac{1}{2} \nabla^2 \delta\psi^*-2 g_c^* |\psi_0|^2\delta\psi^*-(g+i\gamma)\psi_0^{*2}
\delta\psi.\nonumber
\end{aligned}
\end{equation}

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'',
\begin{equation}
\begin{aligned}
&i\partial_t \delta \phi = -\frac{1}{2} \nabla^2 \delta\phi+(2 g_c-g) n\delta\phi+g_c n \delta\phi^*
,\\
&i\partial_t \delta \phi^* = \frac{1}{2} \nabla^2 \delta\phi^*-(2 g_c^* -g)n \delta\psi^*-g_c^* n \delta\phi,
\end{aligned}\nonumber
\end{equation}
with $n = n_0(1+2\gamma n_0t)^{-1}$.
Now if we consider perturbations with momentum $\textbf{k}$ and naively treat $n(t)$ as time invariant, we may find the
eigen-frequencies of the above equations are
\begin{equation}
\omega_\textbf{k} = -2i\gamma n\pm\sqrt{\epsilon_\textbf{k}^2+2gn\epsilon_\textbf{k}-\gamma^2n^2}
\end{equation}
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.

\emph{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~\cite{pethick2008bose,dalfovo1999theory,zhai2021ultracold}. It is then interesting to construct the hydrodynamic equation with two-body loss. Starting from time-dependent dissipative Gross-Pitaevskii equation eq.~\eqref{nGP}, 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:
\begin{equation}
\begin{aligned}
&\partial_t \rho +\nabla(\rho\textbf{v})+2\gamma \rho^2 = 0 \\
m\partial_t\textbf{v} = -\nabla(&\frac{1}{2m}\frac{1}{\sqrt{\rho}}\nabla^2\sqrt{\rho}+\frac{1}{2}m v^2+V(\textbf{r})+g\rho).
\end{aligned}
\end{equation}
 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
\begin{equation}
\rho(\mathbf{r},t) = \frac{\rho_0}{1+2\gamma t \rho_0},
\end{equation}
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.

\section{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.

\section{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.

\bibliographystyle{apsrev4-1_our_style}

redpajama_set_name: RedPajamaArXiv


text: Q: Capitalizing educational Departments In a sentence, is a list of academic departments capitalized? For example, "the departments of education, business, and psychology." Or should it look like this "the Departments of Education, Business, and Psychology."

redpajama_set_name: RedPajamaStackExchange