Stationary states in the many-particle description of Bose-Einstein condensates with balanced gain and loss

TL;DR

This paper investigates the existence and properties of stationary states in many-particle Bose-Einstein condensates with balanced gain and loss, extending mean-field PT-symmetric models to a full quantum description.

Contribution

It demonstrates the presence of non-oscillating trajectories and a non-equilibrium steady state in the many-particle framework, highlighting the effects of particle interactions.

Findings

Existence of two non-oscillating trajectories as many-particle counterparts of mean-field states.

Identification of a non-equilibrium steady state acting as an attractor.

Steady state remains a pure condensate under weak interactions with strong gain and loss.

Abstract

Bose-Einstein condensates with balanced gain and loss can support stationary states despite the exchange of particles with the environment. In the mean-field approximation this is described by the PT-symmetric Gross-Pitaevskii equation with real eigenvalues. In this work we study the role of stationary states in the appropriate many-particle description. It is shown that without particle interaction there exist two non-oscillating trajectories which can be interpreted as the many-particle equivalent of the stationary PT-symmetric mean-field states. Furthermore the system has a non-equilibrium steady state which acts as an attractor in the oscillating regime. This steady state is a pure condensate for strong gain and loss contributions if the interaction between the particles is sufficiently weak.

Figures

Equations

$$\frac{\mathrm{d}}{\mathrm{d}t}\hat{\rho }=-i[\hat{H},\hat{\rho }]+{\gamma }_{\mathrm{loss}}\mathcal{L}({\hat{a}}_1)\hat{\rho }+{\gamma }_{\mathrm{gain}}\mathcal{L}({\hat{a}}_2^{†})\hat{\rho },$$

$$\hat{H}=-J\left({\hat{a}}_1^{†}{\hat{a}}_2+{\hat{a}}_2^{†}{\hat{a}}_1\right)+\frac{U}{2}\left({\hat{a}}_1^{†}{\hat{a}}_1^{†}{\hat{a}}_1{\hat{a}}_1+{\hat{a}}_2^{†}{\hat{a}}_2^{†}{\hat{a}}_2{\hat{a}}_2\right),$$

$$\mathcal{L}({\hat{a}}_1)\hat{\rho }$$

$$\mathcal{L}({\hat{a}}_2^{†})\hat{\rho }$$

$${\gamma }_{\mathrm{loss}}=\frac{N_0+2}{N_0}{\gamma }_{\mathrm{gain}}\equiv \gamma ,$$

$$⟨{\hat{a}}_i^{†}{\hat{a}}_j{\hat{a}}_k^{†}{\hat{a}}_l{\hat{a}}_m^{†}{\hat{a}}_n⟩\approx ⟨{\hat{a}}_i^{†}{\hat{a}}_j{\hat{a}}_k^{†}{\hat{a}}_l⟩⟨{\hat{a}}_m^{†}{\hat{a}}_n⟩$$

$$+⟨{\hat{a}}_i^{†}{\hat{a}}_j{\hat{a}}_m^{†}{\hat{a}}_n⟩⟨{\hat{a}}_k^{†}{\hat{a}}_l⟩+⟨{\hat{a}}_k^{†}{\hat{a}}_l{\hat{a}}_m^{†}{\hat{a}}_n⟩⟨{\hat{a}}_i^{†}{\hat{a}}_j⟩$$

$$-2⟨{\hat{a}}_i^{†}{\hat{a}}_j⟩⟨{\hat{a}}_k^{†}{\hat{a}}_l⟩⟨{\hat{a}}_m^{†}{\hat{a}}_n⟩,$$

$${\hat{L}}_x$$

$${\hat{L}}_z$$

$${\Delta }_{jk}=⟨{\hat{A}}_j{\hat{A}}_k+{\hat{A}}_k{\hat{A}}_j⟩-2⟨{\hat{A}}_j⟩⟨{\hat{A}}_k⟩$$

$${\dot{s}}_x$$

$${\dot{s}}_y$$

$${\dot{s}}_z$$

$$\dot{n}$$

$$P=\frac{s_x^2+s_y^2+s_z^2}{n^2}.$$

$$s_x(t)$$

$$s_y(t)$$

$$s_z(t)$$

$$n(t)$$

$$\mathbf{\alpha }=\frac{{\gamma }_{+}^2-{\gamma }_{-}^2}{4J^2-{\gamma }_{+}^2+{\gamma }_{-}^2}\left(\begin{array}{c}0\\ \frac{2J}{{\gamma }_{-}}\\ 1\\ 1+\frac{4J^2}{{\gamma }_{-}({\gamma }_{+}+{\gamma }_{-})}\end{array}\right).$$

$${\kappa }_1$$

$${\kappa }_2$$

$${\kappa }_3$$

$${\kappa }_4$$

$$k_1$$

$$k_2$$

$$i\frac{\mathrm{d}}{\mathrm{d}t}{c}_1$$

$$i\frac{\mathrm{d}}{\mathrm{d}t}{c}_2$$

$$\phi$$

$$\vartheta$$

$$s_x$$

$$s_y$$

$$s_z$$

$$\phi$$

$$\vartheta$$

$$s_x(t)$$

$$s_y(t)$$

$$s_z(t)$$

$$n(t)$$

$$s_x(0)$$

$$s_y(0)$$

Full text

I Introduction

Non-Hermitian Hamiltonians can be used as an elegant effective description of open quantum systems Moiseyev (2011). A special class of non-Hermitian Hamiltonians that gained much interest since the seminal paper by Bender and Boettcher Bender and Boettcher (1998) are $\mathcal{PT}$-symmetric Hamiltonians which fulfill $[\hat{H},\mathcal{PT}]=0$ with the parity reflection operator $\mathcal{P}$ and the time reversal operator $\mathcal{T}$. The outstanding property of these Hamiltonians is that they can support an entirely real eigenvalue spectrum despite being non-Hermitian Bender and Boettcher (1998); Bender et al. (2002).

In position space non-Hermitian Hamiltonians can be achieved by introducing complex potentials. For Bose-Einstein condensates imaginary parts of the potential have a clear physical interpretation as a localized gain and loss of particles Kagan et al. (1998). If the system is $\mathcal{PT}$ symmetric the complex potential must fulfill $V(\bm{r})=V^{*}(-\bm{r})$. This condition means that the gain and loss contributions are balanced in such a way that stationary states with real eigenvalues can exist.

Due to its simplicity a double-well potential where particles are removed from one well and injected into the other is an ideal system to study the properties of $\mathcal{PT}$-symmetric Bose-Einstein condensates. Both the controlled incoupling and outcoupling of particles have already been demonstrated experimentally. The localized particle loss can be realized by a focused electron beam which ionizes the atoms whereupon they escape the trapping potential Gericke et al. (2008); Würtz et al. (2009); Barontini et al. (2013). Incoupling of particles can be realized by feeding from a second condensate Robins et al. (2008) in a Raman superradiance-like process Döring et al. (2009); Schneble et al. (2004); Yoshikawa et al. (2004).

This system has been extensively studied in the mean-field limit where it is described by the $\mathcal{PT}$-symmetric Gross-Pitaevskii equation using discrete matrix models Graefe (2012); Haag et al. (2015), double-delta potentials Cartarius and Wunner (2012); Cartarius et al. (2012) and spatially extended double wells Dast et al. (2013a, b); Haag et al. (2014); Dizdarevic et al. (2015). In these works it was shown that the system supports stationary $\mathcal{PT}$-symmetric states, exhibits intriguing dynamical properties, and possesses a rich eigenvalue structure, including various exceptional points. Furthermore, proposals to simulate an effective $\mathcal{PT}$-symmetric double well by embedding it into a larger Hermitian structure have been formulated in the mean-field approximation Kreibich et al. (2013, 2014).

However, if the exchange of particles with the environment plays such a crucial role, one has to expect effects that cannot be captured by the mean-field approximation in which a pure condensate described by a product of identical single-particle states is assumed Ruostekoski and Walls (1998); Syassen et al. (2008); Witthaut et al. (2008, 2009, 2011); Barmettler and Kollath (2011); Labouvie et al. (2016). In fact it was shown that in case of a two-mode Bose-Einstein condensate with balanced gain and loss the purity of the condensate oscillates, i.e., it periodically drops to small values but then is almost completely restored Dast et al. (2016a, b). These purity oscillations are experimentally accessible by measuring the average contrast in an interference experiment Dast et al. (2016a). This shows that to fully understand Bose-Einstein condensates with balanced gain and loss a description that goes beyond the mean field is required.

The most characteristic property of $\mathcal{PT}$-symmetric systems is the possibility of the existence of an entirely real eigenvalue spectrum and, consequently, the existence of stationary states despite the Hamiltonian being non-Hermitian. Hence, it is natural to ask whether the many-particle system with balanced gain and loss also supports stationary states and whether these stationary states are connected to the $\mathcal{PT}$-symmetric stationary states of the mean-field limit.

In this work we study a two-mode Bose-Einstein condensate with balanced gain and loss described by a master equation in Lindblad form which is presented in Sec. II. It is shown that the stationary solutions of the mean-field limit are not the stationary states of the master equation but instead they are located close to non-oscillatory trajectories, whose properties are discussed in Sec. III. In Sec. IV the steady state of the many-particle system, defined as $\frac{\mathrm{d}{}}{\mathrm{d}t}\hat{\rho}=0$, is calculated and its influence on the dynamics of the system is studied. It is a non-equilibrium steady state since the system is subject to particle gain and loss, and thus, a compensating particle current has to be present.

II Two-mode system with balanced gain and loss

The many-particle description of a Bose-Einstein condensate with balanced gain and loss in a double-well potential was introduced in Dast et al. (2014) using a two-mode approximation, and its dynamical properties were discussed in Dast et al. (2016a) and Dast et al. (2016b). In this section the properties required for this work are briefly recapped.

The master equation in Lindblad form Breuer and Petruccione (2002); Anglin (1997) describing this system reads

[TABLE]

where the coherent dynamics in the double well is given by the Bose-Hubbard Hamiltonian Fisher et al. (1989); Jaksch et al. (1998) for two lattice sites

[TABLE]

and the localized loss at site 1 and gain at site 2 are introduced by the Lindblad terms

[TABLE]

The parameter $J$ is the tunneling strength between the two sites. The on-site interaction strength between the particles is given by $U$. By demanding that the time derivative of the total particle number vanishes when the expectation values of the particle number are identical at both sites, we obtain a relation for the gain and loss rates $\gamma_{\mathrm{gain}}$ and $\gamma_{\mathrm{loss}}$ Dast et al. (2014),

[TABLE]

with $N_{0}$ the initial number of atoms. In the following $\gamma_{\mathrm{loss}}$ and $\gamma_{\mathrm{gain}}$ are always chosen such that they fulfill this relation, and the abbreviation $\gamma$ will be used to characterize the strength of the gain and loss.

As shown in Dast et al. (2016b) the Bogoliubov backreaction method Anglin and Vardi (2001); Vardi and Anglin (2001) can describe the dynamics of this open quantum system accurately for a limited time span. The basic idea is to formulate equations of motion for the first-order moments $\langle{\hat{a}_{j}^{\dagger}\hat{a}_{k}}\rangle$ which couple to the second-order moments $\langle{\hat{a}_{j}^{\dagger}\hat{a}_{k}\hat{a}_{l}^{\dagger}\hat{a}_{m}}\rangle$ via the interaction term. The coupling of the second-order moments to the third-order moments is eliminated by the approximation

[TABLE]

such that a closed set of equations of motion is obtained for the first- and second-order moments.

For the two-mode system this method is most conveniently formulated using the Bloch representation, where the expectation values $s_{j}=2\langle{\hat{L}_{j}}\rangle$ and $n=\langle{\hat{n}}\rangle$ of the four Hermitian operators

[TABLE]

and the covariances

[TABLE]

with $\hat{A}_{j}\in\{\hat{L}_{x},\hat{L}_{y},\hat{L}_{z},\hat{n}\}$ are used. The resulting set of equations of motions of the first-order moments reads

[TABLE]

with $\gamma_{-}=(\gamma_{\mathrm{loss}}-\gamma_{\mathrm{gain}})/2$ and $\gamma_{+}=(\gamma_{\mathrm{loss}}+\gamma_{\mathrm{gain}})/2$. The ten equations of motion of the covariances can be found in Dast et al. (2016b). Their specific form is not required for the understanding of this paper. The Bloch representation is well suited to discuss the purity of a condensate since it takes the particularly simple form

[TABLE]

This quantity signals the existence of off-diagonal long-range order Penrose and Onsager (1956); Yang (1962), and thus measures how close the condensate is to a pure condensate, described by a product of identical single-particle states.

In the non-interacting limit the equations of motion for the first- and second-order moments decouple and an analytic solution can be obtained Dast et al. (2016b). Note that the results obtained in this limit are exact, however, they are not equivalent to the $\mathcal{PT}$-symmetric Gross-Pitaevskii equation since Eq. (9) can predict fragmentation of the condensate. In the oscillatory regime $\gamma_{+}^{2}<4J^{2}$ the solution reads

[TABLE]

with $\omega=\sqrt{4J^{2}-\gamma_{+}^{2}}$ and the constant inhomogeneous solution, i.e. the steady state, being given by

[TABLE]

The parameters $\kappa_{i}$ in Eq. (11) are defined by the initial state $(s_{x}(0),\ s_{y}(0),\ s_{z}(0),\ n(0))^{T}$,

[TABLE]

The mean-field approximation of the two-mode master equation with balanced gain and loss (1) is obtained in the limit $N_{0}\to\infty$ but with constant macroscopic interaction strength $g=U(N_{0}-1)\approx UN_{0}$. To obtain the Gross-Pitaevskii equation an initially pure condensate is assumed although mean-field theories exist where the bosons reside in several different single-particle states Cederbaum and Streltsov (2003). In the Gross-Pitaevskii-like mean-field limit the condensate is described by a product of identical single-particle states, thus, it is defined by two complex mean-field coefficients $\psi=(c_{1},c_{2})^{T}$. It was shown in Dast et al. (2014) that the mean-field limit is the discrete $\mathcal{PT}$-symmetric Gross-Pitaevskii equation Graefe et al. (2008); Graefe (2012)

[TABLE]

After normalizing the state and choosing an arbitrary global phase only two degrees of freedom remain for a mean-field state which can be expressed using the angles $\varphi$ and $\vartheta$,

[TABLE]

chosen in such a way that the Bloch representation of the first-order moments of a mean-field state with $N_{0}$ particles is given by the spherical coordinates

[TABLE]

In this representation the two stationary $\mathcal{PT}$-symmetric states of Eq. (14) read

[TABLE]

where the upper (lower) sign is for the ground (excited) state of the system.

III Non-oscillatory states

In this section the non-oscillatory states of the many-particle dynamics are discussed in the non-interacting limit $U=0$ and it is shown that they can be interpreted as the equivalent of the $\mathcal{PT}$-symmetric stationary states of the Gross-Pitaevskii equation.

We start from the solutions in the oscillatory regime $\gamma_{+}^{2}<4J^{2}$ given by Eq. (11). As can be directly seen all oscillatory terms vanish for $\kappa_{3}=0$ and a non-oscillating trajectory is given by

[TABLE]

where the only remaining time dependence is the exponential decay towards the steady state $\bm{\alpha}$.

The condition $\kappa_{3}=0$ can be fulfilled by pure initial states which are expressed by the two angles $\varphi$ and $\vartheta$ as defined in Eq. (16). In the following it will be shown that for a specific initial particle number $N_{0}$ there are two pure states which fulfill this condition. They are obtained by solving the set of equations

[TABLE]

The angle $\vartheta$ is directly obtained from the third equation,

[TABLE]

where the steady state (12) and the relation for balanced gain and loss (5) were used.

Before calculating the remaining angle $\varphi$, the two parameters $\kappa_{1}$ and $\kappa_{2}$ must be determined. Equation (19d) yields the value of the parameter $\kappa_{2}$,

[TABLE]

whereas two values for $\kappa_{1}$ follow from the requirement that the initial state must be pure,

[TABLE]

Since a short calculation shows that $s_{y}(0)>0$ in the parameter regime considered, and the two values for $s_{x}(0)$ only differ in their sign, the angle $\varphi$ takes the following two values,

[TABLE]

To interpret the non-oscillatory states as the many-particle equivalent of the stationary $\mathcal{PT}$-symmetric states, it is necessary that they become equal in the limit $N_{0}\to\infty$. In this limit the argument of the arccosine function in Eq. (20) vanishes,

[TABLE]

and in Eq. (23) only the first fraction remains,

[TABLE]

This shows that for $N_{0}\to\infty$ indeed the non-oscillatory states become equal to the $\mathcal{PT}$-symmetric states given by Eq. (17). To be more precise, the state $(\varphi_{-},\vartheta)$ becomes the ground state and $(\varphi_{+},\vartheta)$ the excited state of the mean-field system.

Plotting the spherical coordinates as a function of $\gamma$ for a constant particle number $N_{0}=100$ as done in Fig. 1 yields the characteristic structure of an exceptional point of order $2$. Note that in Fig. 1 as well as in the following figures the value of $J$ is set to 1, without loss of generality. The two non-oscillatory states (solid lines) coalesce slightly before the value of $\gamma$ at which the $\mathcal{PT}$-symmetric ground and excited state (dotted lines) become equal.

Note that both the non-oscillatory and the stationary $\mathcal{PT}$-symmetric states coalesce at values of $\gamma$ that are slightly smaller than the critical point $\gamma_{+}^{2}=4J^{2}$ at which the oscillatory regime ends, as indicated by the vertical black line in Fig. 1. Only in the mean-field limit $N_{0}\to\infty$, the coalescence of all states and the critical point coincide.

The special role of the non-oscillatory states becomes clearly evident in Fig. 2, where the trajectories of the reduced components,

[TABLE]

are shown. In this representation the squared norm of the vector $(s_{x}^{\prime},\ s_{y}^{\prime},\ s_{z}^{\prime})^{T}$ is equal to the purity of the state $P=s_{x}^{\prime 2}+s_{y}^{\prime 2}+s_{z}^{\prime 2}$.

The left panel shows six trajectories for $\gamma=1.5$ and $N_{0}=100$ in red lines which start on the surface of the sphere (which means they are initially pure) at $\vartheta=\pi/2$ and $\varphi$ equally distributed between [math] and $2\pi$. Furthermore, the two non-oscillating trajectories are plotted as blue lines. One immediately recognizes that the dynamics is symmetric with respect to the $s_{y}^{\prime}$-$s_{z}^{\prime}$ plane. This can be understood by replacing $\kappa_{1}$ with $-\kappa_{1}$ in the general solution given by Eq. (11). Then one obtains exactly the same trajectory with opposite sign of the component $s_{x}(t)$, which leads to this symmetry.

Every trajectory has a structure similar to a cone and finally reaches the steady state $\bm{\alpha}$ given by Eq. (12), which lies approximately on the $s_{y}^{\prime}$ axis ($s_{y}^{\prime}\approx 0.75$, corresponding to $P\approx 0.56$). Within the cones lie the trajectories of the non-oscillatory states, which are encircled by all oscillating trajectories. Since $\gamma_{-}$ is the exponential decay rate towards the steady state, the decay is faster for bigger values of $\gamma$ resulting in less narrow windings. If $\gamma$ is increased, the non-oscillatory states on the surface of the sphere are moved towards the $s_{y}^{\prime}$ axis, i.e., approach $\varphi=\pi/2$ (cf. Fig. 1).

The right panel of Fig. 2 shows the inner cone and the non-oscillatory state with $s_{x}^{\prime}>0$ of the left panel in more detail. Furthermore the trajectory of the $\mathcal{PT}$-symmetric ground state is added (green line). This trajectory encircles the non-oscillatory state at a very small distance.

In the mean-field limit initially pure states will stay pure for all times. Thus, they will not leave the surface of the sphere and the two $\mathcal{PT}$-symmetric states are elliptic fixed points, which are encircled by all trajectories. In the many-particle system such fixed points do not exist, instead a non-oscillating trajectory exists which is encircled by all states. In this sense, the non-oscillatory states can be interpreted as the many-particle equivalent of the $\mathcal{PT}$-symmetric stationary states.

If interaction between the particles is taken into account, the $\mathcal{PT}$-symmetric states still show only weakly pronounced oscillations. However, no clearly distinguished non-oscillating trajectories exist in their vicinity. Therefore, it is not possible to define a many-particle equivalent of the $\mathcal{PT}$-symmetric states as done in the non-interacting limit.

IV Non-equilibrium steady state

In this section the non-equilibrium steady state of the two-mode system with balanced gain and loss is studied. It is defined as a constant solution of the master equation, $\frac{\mathrm{d}{}}{\mathrm{d}t}\hat{\rho}=0$. Before investigating the steady state numerically in the presence of interaction, the analytically solvable non-interacting limit is discussed.

IV.1 Non-interacting limit

In the non-interacting limit a constant solution, i.e. a steady state, exists for $\gamma_{+}^{2}\neq 4J^{2}+\gamma_{-}^{2}$ and is given by Eq. (12). Furthermore this solution is an attractor in the parameter regime $\gamma_{+}^{2}<4J^{2}+\gamma_{-}^{2}$ that every trajectory will finally reach. Due to the prefactor the components $\alpha_{2}$, $\alpha_{3}$, and $\alpha_{4}$ diverge for $\gamma_{+}^{2}\to 4J^{2}+\gamma_{-}^{2}$.

The purity of the steady state can be calculated using Eq. (10),

[TABLE]

For every physical state the purity must not exceed one, therefore the steady state is only a physical state in the parameter regime in which $P\leq 1$ holds, which is fulfilled for $\gamma_{+}^{2}\leq 4J^{2}+\gamma_{-}^{2}$. This result is consistent with the dynamical behavior: The steady state has a purity smaller than one, thus being physical, in the parameter regime $\gamma_{+}^{2}<4J^{2}+\gamma_{-}^{2}$, in which it acts as an attractor. At the critical point $\gamma_{+}^{2}=4J^{2}+\gamma_{-}^{2}$ no constant solution exists. Although for $\gamma_{+}^{2}>4J^{2}+\gamma_{-}^{2}$ the state $\bm{\alpha}$ is a constant solution of the equations of motion, it is no longer physical since its purity exceeds one. Such a state cannot be reached dynamically which is consistent with the fact that $\bm{\alpha}$ is no longer an attractor.

Using the relation for balanced gain and loss (5), the purity of the steady state can be reformulated,

[TABLE]

Since $N_{0}\gg\gamma^{2}/4J^{2}$ the purity increases approximately quadratically with $\gamma$

[TABLE]

and for $\gamma=2J$ the steady state is a perfectly pure condensate with $P=1$.

Due to $\gamma_{-}=\gamma/(N_{0}+2)$ the components $\alpha_{2}$ and $\alpha_{4}$ are much larger than the component $\alpha_{3}$. This can be clearly seen by looking at the components divided by the particle number of the steady state $\alpha_{4}$,

[TABLE]

Consequently the purity is almost exclusively determined by $\alpha_{2}$. This component is proportional to the tunneling current from site 2 to site 1, $j_{2\to 1}=J\alpha_{2}$. Thus, $\alpha_{2}/\alpha_{4}$, i.e. the tunneling current of the steady state relative to its total particle number, is positive describing a flux from site 2, where particles are injected, to site 1, where particles are removed. It increases with the strength of the in- and outcoupling $\gamma$. For large particle numbers $N_{0}$ this increase is approximately linear. The quantity $\alpha_{3}/\alpha_{4}$ is the imbalance of the particles in the steady state relative to the total particle number. For $N_{0}\gg\gamma^{2}/4J^{2}$ it is negligible, i.e., the expectation value of the particle number at the two sites is approximately the same.

To discuss the properties of the steady state in more detail the eigenvalues and eigenvectors of the reduced single-particle density matrix $\sigma_{jk}=\langle{\hat{a}_{k}^{\dagger}\hat{a}_{j}}\rangle$ are calculated Yang (1962). The two eigenvalues are determined by the purity, $\lambda_{1/2}=\frac{1}{2}(1\pm\sqrt{P})$, and the corresponding eigenvectors read

[TABLE]

The elements of the two eigenvectors can be interpreted as coefficients of a single-particle state. Of course, the similarity of the two eigenvectors stems from the fact that they are orthogonal since the single-particle density matrix is Hermitian. For $\gamma_{-}\to 0$ the two components of both eigenvectors are equal up to a phase, and for increasing values they diverge approximately linearly since $\gamma_{-}\ll 1$. In the limit of large particle numbers $\gamma_{-}$ can be neglected and the eigenvectors are given by

[TABLE]

Since the tunneling current between the two sites is given by $j_{2\to 1}=2Jr_{1}r_{2}\sin(\beta_{1}-\beta_{2})$ (with $r_{j}\exp(i\beta_{j})$ being the $j$th component of $\bm{u}_{1/2}$) the approximated expression for $\bm{u}_{1}$ describes a single-particle state with maximum tunneling current from site 2 to site 1 while $\bm{u}_{2}$ has a maximum current in the opposite direction.

Due to the incoupling of particles at site 2 and the outcoupling at site 1, a compensating tunneling current from 2 to 1 is required for the steady state. Because the eigenvalue to the eigenstate $\bm{u}_{1}$, which has a tunneling current from 2 to 1, is larger than the other eigenvalue an effective current from site 2 to 1 is achieved. For increasing values of $\gamma$ the eigenvectors remain almost unchanged but the eigenvalue $\lambda_{1}$ increases from $1/2$ to $1$ since the purity $P$ of the steady state increases from [math] to $1$ as can be seen in Eq. (29). Consequently $\lambda_{2}$ decreases towards [math] since $\lambda_{2}=1-\lambda_{1}$. This means a stronger compensating flux is produced by the change in the purity and not by a change of the eigenvectors $\bm{u}_{1/2}$.

IV.2 Numerical results

To obtain the non-equilibrium steady state of the two-mode system with balanced gain and loss numerically, one can solve the equation $\frac{\mathrm{d}{}}{\mathrm{d}t}\hat{\rho}=0$ using an iterative method such as the loose generalized minimal residual method (LGMRES) Baker et al. (2005). The crucial parameter for the numerical cost is the dimension of the Fock basis, which is determined by the maximum amount of particles at each site. To obtain results in a reasonable time span the dimension of the Fock basis of a state vector at a single site is chosen to be $25$, limiting the maximum amount of particles at each site to $24$. It is necessary to ensure that the contributions of states close to this limit are very small since only then it is a good approximation to truncate the basis. For time evolutions this can be achieved by choosing the initial amount of particles $N_{0}$ small enough. However, when searching for the steady state, the particle number is not known in advance but instead is obtained as the result of the calculation for each set of the system’s parameters such as the strength of the gain and loss or the particle interaction.

Therefore, one has to keep in mind that in the following the parameters $N_{0}$ and $g$ are not the particle number and the macroscopic interaction strength of the steady state. Instead, the steady state is reachable by an initial state with $N_{0}$ particles and interaction strength $g=U(N_{0}-1)$. The initial particle number $N_{0}$ only determines the ratio of $\gamma_{\mathrm{gain}}$ and $\gamma_{\mathrm{loss}}$. Since the particle number of the steady state $n$ differs, so does the value of $g=U(n-1)$, whereas the interaction strength between two particles $U$ is held constant.

To analyze the statistical properties of the steady state the probability of finding a certain amount of particles at site 1 and at site 2 is shown in Fig. 3 for $N_{0}=5$, $g=0.5$ and $\gamma=0.5$.

The probabilities at the two sites are almost identical and can be compared with the behavior found in the case of gain and loss at a single site described by

[TABLE]

A study of this master equation in the context of quantum optics can be found, e.g., in Walls and Milburn (2008). Using the ansatz

[TABLE]

shows that the steady state of this system is given by

[TABLE]

By inserting the condition for balanced gain and loss (5) into the probability of finding $j$ particles at the single site (35), one obtains the solid line in Fig. 3, which lies almost perfectly on top of the probabilities found in the two-mode system.

Furthermore, Fig. 3 shows the probability of finding a certain amount of particles at both sites. This quantity cannot be obtained from the probabilities at the single sites since they might be correlated. To check if the correlations are important for the probabilities of finding particles at both sites, a separable state consisting of the steady state of the single-mode calculation is constructed,

[TABLE]

with $p_{j}$ given by Eq. (35). The probability of finding $j$ particles in the two-mode system is the expectation value of the operator $\hat{X}_{j}=\sum_{k=0}^{j}|{k,j-k}\rangle\langle{k,j-k}|$, which yields

[TABLE]

Comparing this probability distribution with the numerically obtained probabilities shows a very good agreement. Only for small particle numbers less than $5$, one obtains a visible discrepancy between the two approaches, and the probabilities given by Eq. (37) are slightly bigger. However, this discrepancy does not arise due to correlations but instead is mainly a consequence of the slightly larger probabilities at each site.

This shows that for the probability distribution of finding a certain amount of particles in the system, the steady state of the two-mode system with balanced gain and loss is very similar to a product of the steady states of the single-mode calculation. However, when looking at other observables there are fundamental differences. In particular the separable state (36) does not have any correlations between the two sites and as a result both $s_{x}$ and $s_{y}$ vanish which is not true for the actual steady state of the two-mode system.

For larger particle numbers the iterative approach to find the steady state is no longer feasible since it becomes numerically too costly. To overcome this limitation the Bogoliubov backreaction method is used. As shown in Dast et al. (2016b) this method is a very good approximation of the actual dynamics for a limited time span. An approximation of the first- and second-order moments of the steady state is obtained by a root search of the equations of motion

[TABLE]

which is solved numerically. Note that in general the Bogoliubov backreaction method is not able to produce accurate results for the dynamics of an arbitrary initial state up to the time at which it reaches the steady state since this exceeds the time span in which this method is reliable. However, it can be expected that the steady state obtained via the root search is nevertheless a good approximation because the short-time behavior of the steady state itself should be well captured. By comparing the results of the Bogoliubov backreaction method with the results obtained by the iterative approach, this assumption was confirmed.

The first-order moments in Bloch representation of the steady state for an initial particle number $N_{0}=100$ are shown in Fig. 4 for different values of the interaction strength $g$.

For comparison the analytically obtained solutions in the non-interacting limit are plotted. One immediately recognizes that the critical value of $\gamma$ at which the steady state ceases to exist becomes considerably smaller for stronger interactions $g$.

Another difference is that the component $s_{x}$, which is equal to zero for the steady state in the non-interacting limit, becomes negative with absolute values that are almost as large as those of $s_{y}$. By contrast, the values of $s_{y}$ are nearly independent of the interaction strength $g$ and are exclusively determined by the strength of the in- and outcoupling $\gamma$. Since $s_{y}$ is the particle flux from site 2 to site 1, which compensates the outcoupling at site 1 and the incoupling at site 2, this behavior is intuitively understandable. The component $s_{z}$ of the steady state behaves similarly to the non-interacting limit. It is close to zero, which means that the steady state is almost balanced, i.e., the particles are equally distributed at both sites. With interaction the total particle number of the steady state does not reach values as large as in the non-interacting limit. This behavior is in agreement with the previous finding that the steady state vanishes earlier for stronger interactions since larger particle numbers imply a stronger macroscopic interaction $g=U(n-1)$. To emphasize that the steady state does not have a definite particle number, but instead there is a broad probability distribution, the uncertainty of the expectation value of the total particle number,

[TABLE]

is plotted alongside the particle number in Fig. 4(d). This shows that the uncertainty of the particle number is almost as large as the expectation value itself.

Note that at $\gamma=0$ the steady state is not well-defined since in this case an infinite amount of stationary solutions exist. Thus, the results shown in Fig. 4 and in the following are strictly only valid for $\gamma>0$.

The purity of the steady state is shown in Fig. 5(a) for different interaction strengths.

In the non-interacting limit the steady state is perfectly pure, i.e. $P=1$, at the critical value of $\gamma$ at which it vanishes. With interaction, however, the steady state disappears at substantially smaller purities.

Calculating the eigenvectors of the reduced single-particle density matrix shows that similar to the non-interacting limit all absolute values of the components stay close to $1/\sqrt{2}$. A slightly different behavior is found for the relative phase. In the case $g=0$ the relative phase is a constant, whereas with interaction the phase increases with the gain-loss parameter $\gamma$. This means that the tunneling current from site 2 to 1 of the single-particle state to the larger eigenvalue decreases. Equally the tunneling current of the second eigenvector in the other direction decreases. The complete tunneling current $s_{y}$, which must compensate the in- and outcoupling of particles, however, increases with $\gamma$ since the steady state becomes more pure and the eigenvalues of the single-particle density matrix are given by $\lambda_{1,2}=\frac{1}{2}(1\pm\sqrt{P})$. Thus, as in the non-interacting case, the increasing compensating current is generated by the steady state becoming more pure and not by a change in the eigenvectors of the single-particle density matrix.

To obtain an overview of the parameter values in which the steady state exists, the root search is repeated for various values of $g$ and $\gamma$ and the result is plotted in Fig. 5(b). The steady state exists for the parameter values in the colored area below the solid black line and the color indicates its purity. It can be seen that the critical value of $\gamma$ at which the steady state vanishes decreases rapidly for stronger particle interactions $g$ and so does the purity.

Up to now the steady state that is reachable in the two-mode system with balanced gain and loss by an initial state with $N_{0}$ particles and macroscopic interaction strength $g=U(N_{0}-1)$ was calculated. However, the macroscopic interaction of the steady state is then given by $g_{\mathrm{steady}}=g(n_{\mathrm{steady}}-1)/(N_{0}-1)$ with the particle number of the steady state $n_{\mathrm{steady}}$. Since the particle number $n_{\mathrm{steady}}$ can become very large close to the critical value of $\gamma$ where the steady state vanishes, this has a marked effect on the system.

A different approach is to ask whether a steady state exists for a specific macroscopic interaction $g$. To achieve this, the parameter $U$ is replaced by

[TABLE]

in the equations of motion of the Bogoliubov backreaction method, with $n(t)$ being the time-dependent particle number. The steady state is then again found via the root search defined in Eq. (38). Note that for this root search the ratio of the gain and loss contributions is still given by Eq. (5), i.e., they are held constant.

The purity obtained in this manner is shown in Fig. 5(c) for different values of $g$. The crucial difference to the results obtained previously is that the steady state now exists for larger values of $\gamma$ resulting in a purer steady state. In particular, a perfectly pure steady state is also possible in presence of interaction as can be seen for $g=0.5$.

Calculating the steady state for various values of $\gamma$ and $g$ yields the results shown in Fig. 5(d). The main difference to Fig. 5(b) is that the critical value of $\gamma$ at which the steady state vanishes decreases much slower for small values of $g$. Consequently, purities equal or close to one are achieved in the presence of interaction. In the parameter regime $\gamma\in[1.5,2]$ and $g\in[0,0.75]$ the steady state is almost pure, thus, there exists a pure condensate which is actually stationary. For stronger interaction strengths $g$ the results obtained by the two approaches become similar.

V Conclusion

In this work the role of stationary states in the many-particle description of a Bose-Einstein condensate with balanced gain and loss was studied.

It was shown that in the many-particle description the outstanding property of the $\mathcal{PT}$-symmetric states of the mean-field description is not that they are stationary but instead that they show only weak oscillations. Without interaction two distinguished trajectories exist in the vicinity of both $\mathcal{PT}$-symmetric states, which show no oscillations at all. For large particle numbers these non-oscillatory states become equal to the $\mathcal{PT}$-symmetric mean-field states. By looking at their time evolution using the Bloch sphere formalism, one immediately recognizes their special role since all other trajectories encircle one of the two non-oscillating trajectories. The fact that in the mean-field limit the $\mathcal{PT}$-symmetric stationary states appear as elliptic fixed points which are encircled by oscillating trajectories motivates the interpretation of the non-oscillatory states as the many-particle equivalent of the $\mathcal{PT}$-symmetric states.

Furthermore, the non-equilibrium steady state of the master equation with balanced gain and loss was calculated. In the non-interacting limit the steady state was obtained analytically. The calculation showed that without coupling to the environment the purity of the steady state is minimum but it increases for stronger values of the gain-loss parameter, and eventually the state becomes completely pure at a critical value. In this parameter regime the steady state is an attractor that every trajectory finally reaches in the limit $t\to\infty$.

With interaction the steady state was calculated numerically showing a similar behavior. However, the steady state now only exists for much smaller gain-loss contributions and it only reaches substantially smaller purities. The reason for this behavior is that the particle number of the steady state increases with the gain-loss parameter, and thus the influence of the interaction becomes very strong. By fixing the macroscopic interaction strength of the steady state, it was shown that a pure steady state can also be reached in the presence of interaction. In fact there is a whole parameter regime in which a pure stationary condensate exists. Finally, calculating the eigenvectors of the single-particle density matrix of the non-equilibrium steady state showed that the increase of the effective current is not generated by a change in the eigenvectors but instead by the steady state becoming more pure, thus, one eigenvector becoming more important than the other one.

Compared with previous studies Witthaut et al. (2009, 2011) where a Bose-Einstein condensate subject to localized particle loss was investigated, we were able to find a true steady state in the presence of particle interaction, whereas in those systems only an exponentially decaying quasi-steady state exists in the non-interacting limit and in an approximation for weak interactions. Furthermore the purification process driven by particle gain and loss Dast et al. (2016a, b) as applied in this study is fundamentally different from that of a system with only loss Witthaut et al. (2009, 2011).