Nonequilibrium steady states of Bose-Einstein condensates with a local particle loss in double potential barriers

TL;DR

This paper studies the stability and bistability of non-equilibrium Bose-Einstein condensates with local particle loss and double barriers, revealing hysteresis phenomena through an exactly solvable mean-field model.

Contribution

It introduces an exactly solvable mean-field model for BECs with local loss and double barriers, demonstrating bistability and hysteresis in non-equilibrium steady states.

Findings

Existence of parameter regions with two stable steady states.

Observation of unidirectional hysteresis phenomena.

Exact solutions enabling stability analysis.

Abstract

We investigate stability of non-equilibrium steady states of Bose-Einstein condensates with a local one-body loss in the presence of double potential barriers. We construct an exactly solvable mean-field model, in which the local loss and the potential barriers take the form of a delta function. Using the exact solutions of our model, we show that there are parameter regions in which two steady-state solutions are dynamically stable, i.e., the model exhibits bistability. We also find that unidirectional hysteresis phenomena appear when the local-loss rate is varied in some parameter region.

Equations

$$i\hslash \frac{\partial }{\partial t}\psi (x,t)$$

$$U(x)$$

$$\left[-\frac{{\hslash }^2}{2M}\frac{d^2}{d{x}^2}+U(x)-\mu +g|\Psi (x){|}^2\right]\Psi (x)=0,$$

$$\Psi (x)\stackrel{x\to \pm \infty }{\to }\sqrt{n_{\infty }}{e}^{-iM{v}_{\infty }|x|/\hslash }{e}^{i{\phi }_{\pm }},$$

$$\mu$$

$$\Psi (\pm L+0)=\Psi (\pm L-0),\Psi (+0)=\Psi (-0),$$

$$\frac{{\hslash }^2}{2M}\left[{\frac{d\Psi (x)}{dx}|}_{x=\pm L+0}-{\frac{d\Psi (x)}{dx}|}_{x=\pm L-0}\right]=U_0\Psi (\pm L),$$

$$\frac{{\hslash }^2}{2M}\left[{\frac{d\Psi (x)}{dx}|}_{x=+0}-{\frac{d\Psi (x)}{dx}|}_{x=-0}\right]=-\frac{i\hslash {\gamma }_0}{2}\Psi (0).$$

$$\frac{\partial }{\partial t}n(x,t)$$

$$n(x,t)$$

$$J(x,t)$$

$$\frac{d}{dt}N(t)$$

$${\Psi }_{\mathrm{PW}}(x)$$

$$v_{\infty }$$

$${\Psi }_{\mathrm{DS}}(x)$$

$$v_{\infty }$$

$${\Psi }_{\mathrm{GS}}(x)$$

$$f(x)$$

$$v_{\infty }$$

$$\Psi (x)$$

$$\frac{n_{\mathrm{out}}(x)}{n_{\infty }}$$

$$\times {\tanh }^2\left[\sqrt{1-{\left(\frac{v_{\infty }}{v_{\mathrm{s}}}\right)}^2}\frac{x-L+x_{+}}{\xi }\right],$$

$${\phi }_{\mathrm{out}}(x)$$

$$-{\tan }^{-1}\left[\frac{G(x+x_{+})}{v_{\infty }/v_{\mathrm{s}}}\right]+{\tan }^{-1}\left[\frac{G(L+x_{+})}{v_{\infty }/v_{\mathrm{s}}}\right],$$

$$\frac{x_{+}}{\xi }$$

$$\times {\tanh }^{-1}\left[\sqrt{\frac{n_L/n_{\infty }-(v_{\infty }/v_{\mathrm{s}}{)}^2}{1-(v_{\infty }/v_{\mathrm{s}}{)}^2}}\right],$$

$$G(x)$$

$$v_{\infty }$$

$$\frac{n_{\mathrm{in}}^{(1)}(x)}{n_{\infty }}$$

$${\phi }_{\mathrm{in}}^{(1)}(x)$$

$$\times \Pi [m_{1}A/(n_0/n_{\infty });\mathrm{am}({\Delta }^{1/4}x/\xi |m_1)|m_1],$$

$$m_1$$

$$\frac{n_{\mathrm{in}}^{(2)}(x)}{n_{\infty }}$$

$${\phi }_{\mathrm{in}}^{(2)}(x)$$

$$\times \Pi \left[\frac{B-n_0/n_{\infty }}{n_0/n_{\infty }};\mathrm{am}({\Delta }^{1/4}x/\xi |m_2)|m_2\right],$$

$$m_2$$

$$A$$

$$B$$

$$\Delta$$

$$n_L$$

$${\Psi }_{\mathrm{out}}(x)$$

$${\Psi }_{\mathrm{in}}(x)$$

$$\tanh \left(\frac{x_0}{\xi }\right)e^{i{\phi }_0}$$

Full text

I Introduction

Ultracold gases are well known as coherent quantum systems with high controllability Bloch2008 . Ultracold gases are confined in a vacuum chamber by using magnetic fields or laser beams such that they are well decoupled from environments. This means that ultracold gases are regarded as isolated quantum systems Polkovnikov2011 . Many interesting phenomena have been studied in the context of isolated quantum systems, such as thermalization Rigol2008 ; Polkovnikov2011 ; Gring2012 ; Kaufman2016 ; Mori2018 and many-body localization Schreiber2015 ; Smith2016 ; Choi2016 ; Nandkishore2015 ; Altman2015 ; Abanin2018a .

Recent technological advances in ultra-cold atom experiments allow us to introduce couplings to the environment, namely, dissipation, in a well-controlled manner Syassen2008 ; Barontini2013 ; Yan2013 ; Labouvie2015 ; Patil2015 ; Labouvie2016 ; Luschen2017 ; Tomita2017 . This means that we can switch ultracold gases from isolated systems to controllable open many-body quantum systems Diehl2008 ; Verstraete2009 ; Daley2014 ; Ashida2016 . The dissipation can be regarded as continuous measurements. When the dissipation is strong compared to other energy scales of the systems, quantum Zeno effects occur Misra1977 , which suppress coherent processes such as tunneling. These effects have been observed in ultracold-gas experiments Syassen2008 ; Mark2012 ; Yan2013 ; Barontini2013 ; Labouvie2015 ; Patil2015 ; Labouvie2016 ; Tomita2017 . It is also noteworthy that the controllable dissipations provide us new possibilities for exploring novel quantum systems, such as $\mathcal{PT}$ symmetric systems Bender1998 ; Ruter2010 ; Konotop2016 ; Xiao2017 ; Li2019 and non-Hermitian quantum systems Bender2007 ; El-Ganainy2018 ; Gong2018 .

Recently, the experimental group at Technische Universitat Kaiserslauten observed bistability in a Bose-Einstein condensate (BEC) with a local particle loss confined in a one-dimensional optical lattice Labouvie2016 . The local particle loss can be realized by focusing an electron beam on the central site of the optical lattice. They prepared two different initial conditions. One is that the central site of the optical lattice is occupied by the particles and the other is that the central site is almost empty. Measuring the particle number of the central site by using scanning electron microscopy techniques, they observed two different stable states. In the small (strong) dissipation regime, the occupied (empty) state is realized regardless of the initial conditions. On the other hand, at the intermediate dissipation strength, the two different stable states are realized depending on the initial states. This means that the system exhibits bistability.

This experiment can be understood as a problem of stability of supercurrents under particle losses. Because the local particle loss induces a density difference between the central site and the others, the supercurrent flows from the surrounding sites into the central sites. The results observed in the experiment indicate that particle losses produce nontrivial effects on superfluidity. In fact, our previous work also showed that global three-body losses induce supercurrent decay in a ring trap Kunimi2019 .

In previous theoretical studies Brazhnyi2009 ; Sels2018a , it has been shown that in the absence of optical lattice potentials, which are described by a real-number external field in the Gross-Pitaevskii (GP) equation, the system does not exhibit a discontinuous jump in the density under a local one-body loss associated with the bistability when the strength of the dissipation is varied. This is contrary to an experiment Labouvie2016 , in which an optical lattice potential is present. In this work, we construct a simple model that is analytically solvable and exhibits the discontinuous jump associated with bistability. Specifically, we use a one-dimensional GP equation with a local one-body loss and double potential barriers, which are, respectively, described by pure imaginary and real delta function potentials. On the basis of semi-analytical solutions of our model, we indeed show that the inclusion of the double potential barriers leads to the emergence of bistability accompanied by the discontinuous jump. In addition, we find unidirectional hysteresis phenomena in our systems. These phenomena are called anomalous hysteresis Yamamoto2012 ; Yamamoto2013_2 ; Yamamoto2013 .

This paper is organized as follows. In Sec. II, we explain the problem that we consider and its formulation based on a dissipative GP equation, which describes a BEC with a local particle loss. In Sec. III.1, using the exact solution of the GP equation, we briefly review important properties of the BEC in the absence of double potential barriers. In Sec. III.2, we obtain exact solutions of the GP equation in the presence of double potential barriers in order to discuss the stability of nonequilibrium steady states of the BEC. In Sec. III.3, we show that our system exhibits anomalous hysteresis phenomena. In Sec. IV, we summarize our results. In the Appendixes, we explain how to perform the stability analysis of stationary solutions of the GP equation and the details of the derivations of the exact solution of the GP equation.

II Model

In this paper, we consider a one-dimensional GP equation with a local one-body loss term and double-potential-barrier terms,

[TABLE]

where $M$ is the mass of the atom, $g>0$ is the two-body interaction coefficient, and $\psi(x,t)$ is the order parameter of the BEC. The dissipation term takes the form of the delta function localized at $x=0$ and $\gamma_{0}\geq 0$ is the strength of the dissipation. The two potential barriers located at $x=\pm L$ are added to mimic the density dips near the local loss created by the optical lattice in the experiment Labouvie2016 . Their strength is denoted $U_{0}\geq 0$. This dissipative GP equation can be derived by the mean-field approximation of the Lindblad equation with the local one-body loss term (see details in the Supplemental Material of Ref. Barontini2013 ).

In Sec. III, we consider the stability of nonequilibrium steady states of a BEC, in which a stationary supercurrent flows into the location of the particle loss. Such states are represented as solutions of the time-independent GP equation, which is derived by inserting $\psi(x,t)=\Psi(x)e^{-i\mu t/\hbar}$ into Eq. (1),

[TABLE]

where $\mu$ is the chemical potential.

We set the boundary condition at $x\to\pm\infty$ as (see also Fig. 1)

[TABLE]

where $n_{\infty}\geq 0$ is the mean particle density at $|x|\rightarrow\infty$, $v_{\infty}\geq 0$ is the magnitude of the flow velocity at $|x|\rightarrow\infty$, and $\varphi_{\pm}$ is the phase. From this boundary condition, we obtain the chemical potential:

[TABLE]

The velocity $v_{\infty}$ is determined by the boundary conditions due to the delta functions, which are given by

[TABLE]

We check the stability of the obtained stationary solutions by the numerical simulations of the time-dependent GP equation. For the details see Appendix A.

At the end of this section, we remark on a crucial difference between our model and the actual experimental setup. In our setup, the particles are lost at the origin and provided at $|x|\rightarrow\infty$ [see Eq. (4)]. This fact can be easily seen by writing down the equation of continuity,

[TABLE]

where $n(x,t)$ and $J(x,t)$ are the particle density and the current density, respectively. Integrating Eq. (9) over $(-\infty,+\infty)$ yields

[TABLE]

where $N(t)\equiv\int^{+\infty}_{-\infty}dx[n(x,t)-n_{\infty}]$ is the total particle number difference at time $t$ Note1 . The first and second terms on the right-hand side of Eq. (12) represent the gain of the particles from the boundaries and the third one represents the loss of the particles at $x=0$. This equation shows that nonequilibrium steady states can be realized when the loss and gain of the particles are balanced.

In the experiment, the BEC is confined in the trap potential with the local particle loss. Because there is no particle source, in contrast to our theoretical setup, the total particle number in the trap monotonically decreases. Hence, strictly speaking, the stationary states cannot exist except in a vacuum state (no particle in the trap). However, according to the inset in Fig. 2 (a) in Ref. Labouvie2016 , we can see that the particle number at the central site is almost stationary over the time scale $40$-$60\;{\rm ms}$. In this time scale, the particle loss and the hopping from the adjacent sites to the central site are balanced. As long as we focus on the vicinity of the central site, the systems can be approximated as nonequilibrium steady states. Stationary states in our model correspond to these nonequilibrium steady states.

Another difference is the width of the local dissipation term. As described above, we assume that the local dissipation is given by the delta function. This treatment can be justified when the width of the dissipation is much smaller than the healing length. However, in the experiment, the width of the dissipation is about $O(0.1\mu{\rm m})$ Barontini2013 . Because the healing length of the experiment is $O(0.1\mu{\rm m})$, the dissipation in the experiment cannot be regarded as the delta function. We also remark on the effects of the finite width in Sec. III.2.

III Results

III.1 In the absence of double potential barriers

For the reader’s convenience we first review exact solutions in the absence of double potential barriers, which have been derived in some previous works Brazhnyi2009 ; Sels2018a before showing our results.

There are three kinds of exact solutions in the absence of potential barriers ($U_{0}=0$). One is a plane-wave (PW) solution:

[TABLE]

The second is a dark soliton (DS) solution,

[TABLE]

where $\xi\equiv\hbar/\sqrt{Mgn_{\infty}}$ is the healing length. The last one is a gray soliton (GS) solution:

[TABLE]

where $v_{\rm s}\equiv\sqrt{gn_{\infty}/M}$ is the sound velocity. We can easily check that these expressions satisfy the GP equation (3). We note that the PW and DS solutions exist for arbitrary parameters and the GS solution exists for $\gamma_{0}>2v_{\rm s}$.

Here, we define the density at the origin as $n_{0}\equiv|\Psi(x=0)|^{2}$, which corresponds to the density at the central site in the experiment Labouvie2016 . We plot $n_{0}$ as a function of $\gamma_{0}$ in Fig. 2, which we call an $n_{0}$-$\gamma_{0}$ diagram. This result shows that the system exhibits bistability in the whole $\gamma_{0}$ region. For $\gamma_{0}\leq 2v_{\rm s}$ the PW and the DS states are stable, and for $\gamma_{0}>2v_{\rm s}$ the GS and the DS states are stable. We can see the unstable PW states for $\gamma_{0}>2v_{\rm s}$. This can be understood by the velocity shown in Fig. 3. The velocity of the PW state is given by Eq. (14), which is proportional to the dissipation strength $\gamma_{0}$. When the velocity exceeds the sound velocity, which is the Landau critical velocity of uniform superfluids Landau1941 ; Wu2003 , the PW state becomes energetically unstable.

The GS state emerges at $\gamma_{0}=2v_{\rm s}\;(v_{\infty}=v_{\rm s})$. The velocity of the GS state is a monotonically decreasing function of $\gamma_{0}$ [see Eq. (19)]. We can interpret this behavior as follows. Suppose that we start with the PW state at $\gamma_{0}=0$. When we increase the dissipation strength from $\gamma_{0}=0$, the superflow velocity becomes high and then reaches the Landau critical velocity. Finally, the PW states become unstable and bifurcate into the unstable PW branch and the stable GS branch.

In the DS states, the density at the origin is always [math]. This means that the DS states do not feel the dissipation. In fact, boundary condition (8) is satisfied in the DS solution (15), for the whole $\gamma_{0}$ region. Therefore, the DS states always exist regardless of the dissipation strength.

III.2 In the presence of double potential barriers

Here, we show the results in the presence of double potential barriers. We assume that a functional form of the stationary solution is given by an even function or an odd function. Owing to this assumption, it is sufficient to consider only the $x\geq 0$ region. Because the potentials are only the delta function type, we can separately solve the GP equation in an inside region ($0\leq x\leq L$) and an outside region ($x>L$). After obtaining the solutions of each region, we connect them by using the boundary conditions (6), (7), and (8). Such techniques for solving the GP equation with delta-function potentials have been developed in the context of Josephson junction systems Baratoff1970 ; Sols1994 ; Hakim1997 ; Kovrizhin2001 ; Pham2002 ; Pavloff2002 ; Kagan2003 ; Astrakharchik2004 ; Seaman2005 ; Bilas2005 ; Danshita2006 ; Danshita2007 ; Watanabe2009 ; Sykes2009 ; Takahashi2009 ; Piazza2010 ; Kato2010 ; Watabe2013 ; Cominotti2014 . For convenience, we introduce the following variables

[TABLE]

First, we consider the even-function case. The solution of the outside region is given by

[TABLE]

where $\varphi_{L}\equiv\varphi(x=L)$ and $n_{L}\equiv n(x=L)$ are determined using the boundary conditions below. $v_{\infty}$ is given by

[TABLE]

This relation can be derived by using the assumption of an even function, the expression of the current density, and the boundary conditions (8). The details of the derivation of the outside solution and Eq. (25) are summarized in Appendix B.

In the inside region, we find four types of inside solutions. However, only two solutions appear in the parameter regions of our interest, where $0\leq\gamma_{0}/v_{\rm s}\leq 4$ and $0\leq n_{0}/n_{\infty}\leq 1$. Then we consider two types of solutions:

[TABLE]

where we have set the origin of the phase as $\varphi_{\rm in}^{(i)}(x=0)=0$ and used the Jacobi elliptic functions ${\rm sn}(x|m)$ and ${\rm nd}(x|m)\equiv 1/{\rm dn}(x|m)$, the incomplete elliptic integral of the third kind $\Pi(n;\phi|m)$, and the Jacobi amplitude function ${\rm am}(x|m)$. The notations for the Jacobi elliptic functions and the elliptic integrals follow by Abramowitz and Stegun Abramowitz_Stegun . We also used the following quantities:

[TABLE]

From the above results and boundary condition (6), $n_{L}$ and $\varphi_{L}$ are determined by

[TABLE]

Next, we consider the odd-function case. From this assumption, we obtain $\Psi(x=0)=0$. This means that the odd-function solution does not depend on $\gamma_{0}$ (see the descriptions of the DS in Sec. III.1). From the equation of continuity, the current density is independent of $x$. In this case, $J(x)=0$ because $\Psi(0)=0$. Therefore, the odd-function solution does not carry a supercurrent and we can take $\Psi(x)$ as a real function without loss of generality. The solution is given by

[TABLE]

where $\varphi_{0}=0\text{ or }\pi$, and $x_{0}$ and $m_{0}$ are constants. $\varphi_{0}$ and $x_{0}$ are determined by boundary condition (6):

[TABLE]

Although the functional forms of the exact solution have been derived, $n_{0}$ (for the even-function case) and $m_{0}$ (for the odd-function case) have not been determined yet. These variables can be determined by solving boundary condition (7). Unfortunately, we cannot solve Eq. (7) analytically. We solve Eq. (7) numerically. The details of the derivations of these solutions are reported in Appendixes B and C.

Here, we remark on the range of $L$. From Eqs. (26) and (29), we can find that the inside solutions have periodicity $2K(m_{1})\xi/\Delta^{1/4}$ and $2K(m_{2})\xi/\Delta^{1/4}$ due to the properties of the Jacobi elliptic functions, where $K(\cdot)$ is the complete elliptic integral of the first kind. If $L$ is much larger than these periods, we can expect that there are solutions that oscillate multiple times in the inside region. To avoid the complexity of the problem, we restrict the range of $L$ to $0\leq L\lesssim 3.3$, which means that the number of oscillations in the inside region is less than 1.

In the presence of double potential barriers, we find five types of $n_{0}$-$\gamma_{0}$ diagrams. The parameter region for the $n_{0}$-$\gamma_{0}$ diagrams is shown in Fig. 4.

A typical type 1diagram is shown in Fig. 5. In type 1, we have two stable branches. One is the even function (upper branch) and the other is the odd function (lower branch). The type 1 solution tends to exist in a region where $U_{0}$ is small. This means that type 1 can be interpreted as perturbed $U_{0}=0$ states. In fact, the $n_{0}$-$\gamma_{0}$ diagram in Fig. 4 is similar to that in Fig 2 except for the existence of the unstable PW branch.

Type 2 emerges in the region adjacent to type 1. A typical $n_{0}$-$\gamma_{0}$ diagram is shown in Fig. 6. In type 2, we can see the discontinuous jump between the upper branch and the lower branch. A similar discontinuous jump has been observed in experiments Labouvie2016 . In contrast, there is no discontinuous jump in the absence of potential barriers (see Fig. 2). This result means that the discontinuous jump is due to the effects of potential barriers.

We show a typical $n_{0}$-$\gamma_{0}$ diagram of type 3 in Fig. 7. In type 3, the upper and lower branches are completely separated. We can see a saddle-node bifurcation in the upper branch, in which two fixed points collide with each other and are annihilated Strogatz . This behavior is similar to that of Josephson junction systems. Theoretically, these systems have been studied using the GP equation or the Ginzburg-Landau equation with a single potential barrier Baratoff1970 ; Sols1994 ; Hakim1997 ; Kovrizhin2001 ; Pham2002 ; Pavloff2002 ; Kagan2003 ; Astrakharchik2004 ; Seaman2005 ; Bilas2005 ; Danshita2006 ; Danshita2007 ; Watanabe2009 ; Sykes2009 ; Takahashi2009 ; Piazza2010 ; Kato2010 ; Watabe2013 ; Cominotti2014 . In fact, our system can be regarded as a connection of two reverse Josephson junction systems via local loss. The upper branch is reflected by the properties of the Josephson junction, i.e., superfluidity.

Type 4 emerges in a narrow region surrounded by types 2, 3, and 5 (see inset in Fig. 4). A typical $n_{0}$-$\gamma_{0}$ diagram is shown in Fig. 8. In type 4, the upper and lower branches are similar to those of type 2 and one additional branch emerges between the upper and the lower branches.

A typical $n_{0}$-$\gamma_{0}$ diagram of type 5 is shown in Fig. 9. Type 5 is located between type 1 and type 4. Type 5 is similar to type 4 except for the upper branch. The upper branch of type 5 is similar to that of type 1.

From the above results, we can see the bistability for the whole $\gamma_{0}$ region in types 1 and 5 and partial regions in types 2, 3, and 4. The difference between the presence and the absence of potential barriers is the existence of the discontinuous jump, which can be seen in types 2, 3, and 4.

Comparing our results with the experimental ones, we find that our results are in part inconsistent with the experiment Labouvie2016 . In the small-t$\gamma_{0}$ region, while only one stable state was observed in the experiment, there are two stable states in our model, one of which is the DS state. One possible reason for this discrepancy is that the local particle loss is modeled as a delta-function form.

III.3 Anomalous hysteresis

In addition to the bistability, the present system exhibits a nontrivial hysteresis phenomenon, which is called anomalous hysteresis Yamamoto2012 ; Yamamoto2013_2 ; Yamamoto2013 . A feature of anomalous hysteresis is unidirectionality. In conventional hysteresis phenomena, if we observe a discontinuous jump from an initial phase to another phase upon changing the parameters sufficiently slowly, another jump going back to the initial phase exists along the reverse path in the parameter space. However, in anomalous hysteresis, the discontinuous jump exists only in one direction. This phenomenon has been predicted in quantum phase transitions of dipolar or multicomponent Bose gases in an optical lattice Yamamoto2012 ; Yamamoto2013 and frustrated magnets Yamamoto2013_2 , and it can be understood within the framework of the Ginzburg-Landau theory.

Here, we focus on type 3. The processes of anomalous hysteresis are shown in Fig. 10. First, we prepare the initial state at point (1) shown in Fig. 10 (a). Then we increase $\gamma_{0}$ sufficiently slowly. When the dissipation strength reaches the critical value, the discontinuous jump occurs from the upper branch to the lower branch. After the discontinuous jump, we decrease the dissipation strength and, finally, reach point (2) shown in Fig. 10 (a). Next, let us consider the inverse process; that is, the initial state is point (2) in Fig. 10 (b) and the goal is point (1) in Fig. 10 (b). However, this process is impossible because the lowest branch is stable for the whole $\gamma_{0}$ region. This means that we cannot reach point (1) starting from point (2) as long as we consider sufficiently slow changes of the parameters. This is nothing but the anomalous hysteresis phenomenon mentioned above.

Here, we discuss the time scale of changing parameters. For example, let us consider the case shown in Fig. 10. The energy difference between the upper stable branch and the lower stable branch is given by the energy of the dark soliton, which is of the order of the chemical potential. In the actual experiments, the chemical potential is typically of the order of $1\;{\rm kHz}$. The inverse of this energy scale gives us the time scale of changing parameters. Therefore, we should change the parameters within a time scale slower than $1{\rm ms}$. This condition can be easily satisfied in cold-gas experiments. We also remark on the adiabatic condition of the system. The adiabatic condition of the present system, which corresponds to the condition where no excitation is present, is given by the Bogoliubov excitation. We can roughly estimate the times cale to be $100\;{\rm ms}$ for the system size $O(10\xi)$. This means that Bogoliubov excitations are present in the experiment Labouvie2016 because the experimental timescale is shorter than $100{\rm ms}$. Nevertheless, the hysteresis loop can be clearly observed. This indicates that the adiabaticity is not a necessary condition but a sufficient condition for observing the hysteresis loop.

At the end of this section, we discuss the feasibility of observing anomalous hysteresis in experiments. Thus far, anomalous hysteresis has not been observed experimentally for the following reasons. In the case of dipolar or multicomponent Bose gases in an optical lattice Yamamoto2012 ; Yamamoto2013 , the temperature in the optical lattice has not been lowered enough to observe anomalous hysteresis. In the case of frustrated magnets Yamamoto2013_2 , it is difficult to tune the parameters to the optimal values for observing anomalous hysteresis.

In contrast to the previous works, there is no difficulty with our model in achieving sufficiently low temperatures and optimal values of the parameters. However, anomalous hysteresis has not been observed in experiments Labouvie2016 . There are a few possible reasons for this discrepancy. One is the effects of the harmonic trap. The presence of the trap potential may affect the hysteresis because it changes the boundary condition of the system. In our systems, we fix the wave function as the plane wave at $\infty$. This means that the particles are provided from the bath. This situation is different from the experimental setup, which is isolated from the environment except for the local loss. We also remark that this boundary condition produces an additional nonlinearity. The combination of the boundary conditions at $\infty$ and at the origin determines the velocity at $\infty$ [see Eq. (25)]. The velocity depends on the density at the origin. This constraint does not exist in the experimental setup. This difference may affect the existence of anomalous hysteresis. Another one is the effects of optical lattices. The optical lattice extends over the entire system. In contrast to this, in our system, the double delta potentials are localized near the center of the system. This difference may affect the hysteresis. In addition to these points, the width of the local dissipation may affect the stability as discussed in Sec. III.2.

IV Summary and future prospects

We have investigated the stability of a BEC with a local one-body loss in double potential barriers by using the mean-field approximation. We obtained the exact solutions of the GP equation in the presence of delta-function potentials with the pure imaginary and real coefficients, which are written by the Jacobi elliptic functions. We showed that there is a wide parameter region, in which two nonequilibrium steady states are dynamically stable, i.e., our model exhibits bistability. We also found the anomalous hysteresis phenomenon in our system.

As a future plan, we will investigate the effects of the width of the local dissipation and the optical lattice potentials. These effects may change the stability of the present system. By studying these effects, we may clarify the origin of the bistability observed in the experiment.

It is interesting to extend our analysis to strongly correlated regimes. Our model is based on the mean-field theory, which can be justified only in weakly correlated regimes. Strongly correlated nonequilibrium states are one of the most difficult problems in various fields. As a topic related to bistability, negative differential conductivity is theoretically studied by using anti-de Sitter space and conformal field theory correspondence Nakamura2012 .

Another extension is to consider the effects of local multi-body losses, for example, two-body and three-body losses. Particularly, controllable global two-body losses have been realized using the photo-association laser Tomita2017 . By developing this kind of experimental technique, controllable local two-body losses will be experimentally realized.

Acknowledgements.

M. K. thanks D. A. Takahashi for his lecture on the Jacobi elliptic function, S. Takada for useful discussion, and T. Tomita for useful comments. M. K. was supported by Grant-in-Aid for JSPS Research Fellow No. JP16J07240. I. D. was supported by JSPS KAKENHI Grants No. 15H05855, No. 25220711, No. 18K03492, and No. 18H05228, by a research grant from CREST, JST, and by the Q-LEAP program of MEXT, Japan.

Appendix A STABILITY ANALYSIS

Here, we explain how to perform the stability analysis of the stationary state. To do this, we investigate real-time dynamics. However, we do not use Eq. (1) because of some technical reasons described below.

The original problem is defined by an infinite-sized system. However, this system is not tractable numerically. Instead of considering the infinite systems, we consider the finite-size system $(-L_{\rm s},+L_{\rm s})$, where we take $L_{\rm s}$ to be about $100\xi$. The equation considered here is given by

[TABLE]

where we have introduced the spatially varying dissipation term $\Gamma(x)$. The reason we introduce the dissipation term is to avoid effects of the reflection of the boundary, which does not exist in the original problem. The functional form of the dissipation $\Gamma(x)$ is the same as that used in Ref. Reeves2015 . The parameters are set to $L_{\rm d}=L_{\rm s}/2$ and $W=10\xi$. We note that the choice of these parameters is insensitive to the results as long as $L_{\rm d},W\gg\xi$ are satisfied. We also introduce the time dependence of the chemical potential to converge to the stationary solution at the long time. The boundary condition at the edge of the system is given by

[TABLE]

We numerically solve Eq. (39) by using the fourth-order Runge-Kutta method. The centered difference method is used for the space discretization. We use the number of meshes $N_{x}=2001$-$64001$. In this calculation, we approximate the delta function as the Kronecker $\delta(x-x_{j})\simeq(1/\Delta x)\delta_{i,j}$, where $x_{i}\equiv\Delta x\times i\;[i=-(N_{x}-1)/2,\cdots,+(N_{x}-1)/2]$ and $\Delta x$ is the mesh size. We write the discretized wave function at mesh $i$ and time $t$ as $\psi_{i}(t)$. We have checked that the analytically obtained stationary solutions and the numerically obtained stationary solutions are in good agreement.

The procedure of the stability analysis is as follows. We use the initial conditions as the exact solution plus small random noise. That is, the initial condition is given by $\psi_{j}(0)=\psi_{\rm exact}(x_{j})+\epsilon_{j}^{\rm R}+i\epsilon_{j}^{\rm I}$, where $\psi_{\rm exact}(x_{j})$ is the exact solution at mesh $j$ and $\epsilon_{j}^{\rm R}$ and $\epsilon_{j}^{\rm I}$ are real values. We set $-10^{-4}\leq\epsilon_{j}^{\rm R},\epsilon_{j}^{\rm I}\leq 10^{-4}$. Then we numerically calculate the real-time dynamics. After long-time evolution [typically $1000\tau\sim 10000\tau$, where $\tau\equiv\hbar/(gn_{\infty})$], we compare the final state with the initial state.

A typical example of the time evolution is shown in Fig. 11 (a). We see the dynamics of $n_{0}(t)$ for type 3. We can see that the lowest branch (1) and uppermost branch (4), shown by red lines, are stable against a small perturbation in the initial states. On the other hand, branches (2) and (3), shown by blue lines, are unstable. The instability sets in at $t\sim 6000\tau$ for branch (2) and $t\sim 20\tau$ for branch (3), respectively. In order to quantify the instability, we calculate the quantity Cassidy2009 :

[TABLE]

where $\psi_{i}(t)$ is the wave function at mesh $i$ at time $t$. When $\lambda(t)$ becomes exponentially large, dynamical instability occurs. Figure. 11 (b) shows the time evolution of $\lambda(t)$ for the same parameter as in Fig. 11 (a). The results show that the values of $\lambda(t)$ for branches (1) and (4) are less than $1$ at all times, while those for branches (2) and (3) are exponentially large after the instability occurs. From these results, we can conclude that branches (1) and (4) are stable and branches (2) and (3) are unstable. In the same manner, we can judge the stability of the exact solutions with other parameters.

Appendix B DETAILS OF THE DERIVATION OF THE EXACT SOLUTIONS FOR THE EVEN-FUNCTION CASE

In this Appendix, we describe the details of the derivation of the exact solutions for the even-function case. As we described in Sec. III.2, it is sufficient to consider only the region of $x>0$.

First, we derive Eq. (25). From boundary condition (4) and the equation of continuity (9), we obtain the current density in stationary states as

[TABLE]

where ${\rm sgn}(\cdot)$ is the sign function. The boundary condition due to the local loss potential (8), can be written as

[TABLE]

where we have used the assumption of an even function. Using the second Eq. (47) and the expression of the current density

[TABLE]

we obtain Eq. (25); $v_{\infty}=(n_{0}/n_{\infty})\gamma_{0}/2$.

Then we consider solving the GP equation. We define

[TABLE]

It can be easily shown that $C(x)$ is a constant for $0\leq x\leq L$ and $x>L$. Substituting $\Psi(x)=\sqrt{n(x)}e^{i\varphi(x)}$ and $J(x)$ into Eq. (49), we obtain

[TABLE]

Here, we consider the outside region ($x>L$). In this region, we obtain $C(x)=(1/2)gn_{\infty}^{2}+Mv_{\infty}^{2}n_{\infty}$ from the boundary condition at $x\rightarrow\infty$ (4). Equation (50) in the outside region reduces to

[TABLE]

We mention that Eq. (51) is related to the problem of classical mechanics. Here, we consider a classical particle under the potential $V(x)$. In this case, the energy of the system is given by

[TABLE]

where $x(t)$ is the position of the classical particle at time $t$, we set the mass of the particle to $m=1/2$, and $E$ is the total energy. When $V(x)-E=-(x-1)^{2}(x-a)$  $(0<a<1)$, this equation has the same structure as Eq. (51). We assume that $x(t\to\infty)=1$, which corresponds to the boundary condition of the density $n_{\rm out}(x)/n_{\infty}\to 1$ at $x\to\infty$. From this, we can obtain the information on the motion under the potential $V(x)$ in an intuitive way. Figure 12 shows the potential. From Eq. (52), the motion is possible only if $V(x)-E\leq 0$. Here, we set the initial condition $x(0)=b$. When $b\geq 1$, we obtain $dx(t)/dt\leq 0$. When $b<1$, we have two cases: one is that $x(t)$ is monotonically approaching $1$, and the other is that $x(t)$ is bounced at $x(t)=a$ and goes to $1$. The difference comes from the sign of the initial condition $dx(t)/dt|_{t=0}$.

From the above discussion, we can expect that there are three types of solutions in the outside region. From Eq. (51) we obtain

[TABLE]

Here, we consider the case $n_{L}/n_{\infty}<1$. In this case, we can show $n_{L}\leq n_{\rm out}(x)\leq n_{\infty}$ from Eq. (51) and perform the integral in Eq. (53); then we obtain Eq. (21),

[TABLE]

To perform the integral, we used the integral formula

[TABLE]

where this formula is valid for $(bp-aq)p>0$. In the case of $n_{L}/n_{\infty}>1$, we can obtain a different solution, whose functional form is given by replacing $\tanh$ with $\coth$ in Eq. (21). However, we cannot find the parameter region where this solution satisfies the boundary conditions. Therefore, we do not consider the case $n_{L}/n_{\infty}>1$ in the text.

The phase of the outside region can be obtained by integrating Eq. (48). Its expression is given by

[TABLE]

To perform the integral, we used the mathematical formulas:

[TABLE]

where $F(x)$ is a smooth function.

We can obtain the constraint of the velocity $v_{\infty}$ from the above results. From Eq. (51), $n_{\rm out}(x)/n_{\infty}\geq(v_{\infty}/v_{\rm s})^{2}$ must hold. Using $n_{\rm out}(x)/n_{\infty}\leq 1$, we obtain the relation

[TABLE]

This means that the velocity of the stationary solution is always subsonic. This is consistent with the well-known results for the condition of the existence of a gray soliton in uniform systems.

Now, we consider the inside region ($0<x<L$). Using the first Eq. (47), Eq. (25), and Eq. (50), we can determine $C_{\rm in}\equiv C(x)\;(\text{for }x<L)$ in the inside region:

[TABLE]

From Eq. (61), we can rewrite (50) in the inside region as

[TABLE]

where $A$ and $B$ were defined by Eqs. (32) and (33):

[TABLE]

We can integrate Eq. (62) in a similar manner to the case of the outside region. The corresponding potential of the classical mechanics is given by

[TABLE]

In this case, the initial condition is given by $x(0)=x_{0}$, which corresponds to $n(x=0)=n_{0}$. To perform the integral, we need to know the relation between $A$ and $B$. When $\Delta\geq 0$, we obtain $A\geq B$ from Eqs. (32), (33), and (34). Therefore, we classify the solutions as four types: solution 1, $\Delta\geq 0,\;B\leq x_{0}\leq A$; solution 2, $\Delta\geq 0,\;x_{0}\leq B\leq A$; solution 3, $\Delta\geq 0,\;B\leq A\leq x_{0}$; and solution 4, $\Delta<0$. The behavior of the potential is shown in Fig. 13.

Here, we consider solution 1. From the inequalities $\Delta\geq 0$, $B\leq n_{0}/n_{\infty}\leq A$, and Eq. (60), this solution exists in the region

[TABLE]

We plot the parameter region in Fig. 14. Solution1 [Eq. (26)] can be obtained by integration of Eq. (62),

[TABLE]

where we have used formula 17.4.63 in Ref. Abramowitz_Stegun . The phase (69), is also obtained by integrating Eq. (48). To perform this, integral, we used the relations:

[TABLE]

where $F(\phi|m)$ is the incomplete elliptic integral of the first kind.

The region where solution 2 exists is derived by $\Delta\geq 0,\;n_{0}/n_{\infty}\leq B\leq A$, and Eq. (60):

[TABLE]

The expression of solution 2 is given by

[TABLE]

To obtain Eqs. (74) and (75), we used formula 17.4.62 in Ref. Abramowitz_Stegun .

The region of solution 3 is derived by $\Delta\geq 0,\;B\leq A\leq n_{0}/n_{\infty}$, and Eq. (60):

[TABLE]

The expression of solution 3 is given by

[TABLE]

where ${\rm sc}(x|m)\equiv{\rm sn}(x|m)/{\rm cn}(x|m)$ and we have used formula 17.4.64 in Ref. Abramowitz_Stegun .

The region of solution 4 is derived by $\Delta<0$ and Eq. (60):

[TABLE]

The expression of solution 4 is given by

[TABLE]

where we have defined

[TABLE]

Here, we have used formulas 16.18.4 and 17.4.71 in Ref. Abramowitz_Stegun .

What remains to do is to determine the parameters $n_{L}$, $\varphi_{L}$, and $n_{0}$ by connecting the inside and the outside solutions via the boundary conditions. $n_{L}$ and $\varphi_{L}$ are determined by the first expression of Eq. (6), that is, $n_{L}^{(i)}=n_{\rm in}^{(i)}(x=L)$ and $\varphi^{(i)}_{L}=\varphi_{\rm in}^{(i)}(x=L)$. The explicit expressions for the density are given by

[TABLE]

$n_{0}$ is determined by boundary condition (8), which reduces to

[TABLE]

Equation (95) is automatically satisfied due to the expression of the current density (48). Equation (94) reduces to

[TABLE]

Because $n_{L}^{(i)}$ is a function of $n_{0}$, Eq. (96) is a one-variable equation of $n_{0}$ for fixed $\gamma_{0}$ and $U_{0}$. Therefore, the problem of solving the GP equation (nonlinear differential equation) reduces to solving the one-variable equation (96). Because we cannot obtain the analytical solutions of Eq. (96), we solve this equation numerically.

As mentioned in the text, we cannot find the parameter region where solutions 3 and 4 satisfy the boundary conditions. This means that $n_{0}/n_{\infty}$ moves only $0\leq n_{0}/n_{\infty}\leq 1$.

Appendix C DETAILS OF THE DERIVATION OF THE EXACT SOLUTIONS FOR THE ODD-FUNCTION CASE

In this Appendix, we show the expression of boundary condition (7) for the odd-function case. Substituting Eqs. (36), (37), and (38) into Eq. (7), we obtain

[TABLE]

where $\varphi_{0}$ has been determined by Eq. (38).