Non-linear stationary solutions in realistic models for analog black-hole lasers

TL;DR

This paper proposes realistic models for black-hole lasers in Bose-Einstein condensates, demonstrating how to create large supersonic regions and analyzing stationary states using experimentally feasible setups.

Contribution

It introduces experimentally implementable models for black-hole lasers, extending configurations symmetrically to achieve large supersonic regions and analyzing their stationary states.

Findings

Black-hole laser configurations can be symmetrically extended to have arbitrarily large supersonic regions.

Proposed models include an attractive square well and a double delta-barrier for experimental realization.

Identified the true ground state and discussed dynamical instabilities in these setups.

Abstract

From both a theoretical and an experimental point of view, Bose-Einstein condensates are good candidates for studying gravitational analogues of black holes and black-hole lasers. In particular, a recent experiment has shown that a black-hole laser configuration can be created in the laboratory. However, the most considered theoretical models for analog black-hole lasers are quite difficult to implement experimentally. In order to fill this gap, we devote this work to present more realistic models for black-hole lasers. For that purpose, we first prove that, by symmetrically extending every black-hole configuration, one can obtain a black-hole laser configuration with an arbitrarily large supersonic region. Based on this result, we propose the use of an attractive square well and a double delta-barrier, which can be implemented using standard experimental tools, for studying black-hole lasers. We also compute the different stationary states of these setups, identifying the true ground state of the system and discussing the relation between the obtained solutions and the appearance of dynamical instabilities.

Equations

$$\begin{gathered} \hat{H}=\int {\mathrm{d}}^3\mathbf{x}{\hat{\Psi }}^{†}(\mathbf{x})\left[-\frac{{\hslash }^2}{2m}{\nabla }^2+V_{\mathrm{ext}}(\mathbf{x},t)\right]\hat{\Psi }(\mathbf{x}) \\ +\frac{g_{3D}}{2}{\hat{\Psi }}^{†}(\mathbf{x}){\hat{\Psi }}^{†}(\mathbf{x})\hat{\Psi }(\mathbf{x})\hat{\Psi }(\mathbf{x}) \end{gathered}$$

$$\hat{\Psi }(\mathbf{x})\simeq \hat{\psi }(x)\frac{e^{-\frac{{\rho }^2}{2a_{\mathrm{tr}}^2}}}{\sqrt{\pi }{a}_{\mathrm{tr}}},a_{\mathrm{tr}}=\sqrt{\frac{\hslash }{m{\omega }_{\mathrm{tr}}}}$$

$$\begin{gathered} {\hat{H}}_{1D}=\int \mathrm{d}x{\hat{\psi }}^{†}(x)\left[-\frac{{\hslash }^2}{2m}{\partial }_x^2+V(x,t)\right]\hat{\psi }(x) \\ +\frac{g_{1D}}{2}{\hat{\psi }}^{†}(x){\hat{\psi }}^{†}(x)\hat{\psi }(x)\hat{\psi }(x) \end{gathered}$$

$$[\hat{\Psi }(\mathbf{x}),{\hat{\Psi }}^{†}({\mathbf{x}}^{\prime })]=\delta (\mathbf{x}-{\mathbf{x}}^{\prime }),$$

$$[\hat{\psi }(x),{\hat{\psi }}^{†}(x^{\prime })]=\delta (x-x^{\prime }),$$

$$i\hslash {\partial }_t\hat{\psi }(x,t)=[\hat{\psi }(x,t),\hat{H}]=\left[-\frac{{\hslash }^2}{2m}{\partial }_x^2+V(x,t)\right]\hat{\psi }(x,t)+g{\hat{\psi }}^{†}(x,t)\hat{\psi }(x,t)\hat{\psi }(x,t)$$

$$\hat{\psi }(x,t)=\psi (x,t)+\hat{\phi }(x,t),$$

$$i\hslash {\partial }_t\psi (x,t)=\left[-\frac{{\hslash }^2}{2m}{\partial }_x^2+V(x,t)+g|\psi (x,t){|}^2\right]\psi (x,t)$$

$$N=\int \mathrm{d}x|\psi (x,t){|}^2$$

$${\partial }_t|\psi (x,t){|}^2+{\partial }_{x}J(x,t)=0,J(x,t)=-\frac{i\hslash }{2m}\left[{\psi }^{\ast }(x,t){\partial }_x\psi (x,t)-\psi (x,t){\partial }_x{\psi }^{\ast }(x,t)\right]$$

$${\partial }_{t}n(x,t)+{\partial }_x[n(x,t)v(x,t)]$$

$$-\hslash {\partial }_t\varphi (x,t)$$

$$n(x,t)=A^2(x,t),v(x,t)=\frac{\hslash {\partial }_x\varphi (x,t)}{m},$$

$$m{\partial }_{t}v(x,t)=-{\partial }_x\left[\frac{1}{2}m{v}^2(x,t)+V(x,t)\right]-\frac{1}{n(x,t)}{\partial }_{x}P(x,t)+{\partial }_x\left[\frac{{\hslash }^2}{2m\sqrt{n(x,t)}}{\partial }_x^2\sqrt{n(x,t)}\right],P(x,t)=\frac{g{n}^2(x,t)}{2}$$

$$i\hslash {\partial }_t\hat{\Phi }(x,t)$$

$$M(x,t)$$

$$G(x,t)$$

$$\hat{\psi }(x,t)=[{\psi }_0(x)+\hat{\phi }(x,t)]e^{-i\frac{\mu }{\hslash }t}$$

$$\mu {\psi }_0(x)=\left[-\frac{{\hslash }^2}{2m}{\partial }_x^2+V(x)+g|{\psi }_0(x){|}^2\right]{\psi }_0(x)$$

$${\partial }_{x}J(x)=0$$

$$\mu A(x)=-\frac{{\hslash }^2}{2m}{A}^{\prime \prime }(x)+\frac{m{J}^2}{2A^3(x)}+V(x)A(x)+g{A}^3(x)$$

$$\varphi (x)=\int \mathrm{d}x\frac{mv(x)}{\hslash }=\int \mathrm{d}x\frac{mJ}{\hslash A^2(x)}$$

$$K[\psi ]=\int \mathrm{d}x{\psi }^{\ast }(x)\left[-\frac{{\hslash }^2}{2m}{\partial }_x^2+V(x)-\mu \right]\psi (x)+\frac{g}{2}|\psi (x){|}^4,$$

$$K[{\psi }_0]=-\int \mathrm{d}x\frac{g}{2}|{\psi }_0(x){|}^4.$$

$$i\hslash {\partial }_t\hat{\Phi }(x,t)$$

$$M(x)$$

$$G(x)$$

$$\hat{\Phi }(x,t)$$

$$z_n(x)$$

$$M(x)z_n(x)={ϵ}_n{z}_n(x),{ϵ}_n=\hslash {\omega }_n$$

$$(z_1|z_2)\equiv ⟨z_1⟩{\sigma }_z{z}_2=\int {\mathrm{d}}^3\mathbf{x}{z}_1^{†}(\mathbf{x}){\sigma }_z{z}_2(\mathbf{x})=\int {\mathrm{d}}^3\mathbf{x}{u}_1^{\ast }(\mathbf{x})u_2(\mathbf{x})-v_1^{\ast }(\mathbf{x})v_2(\mathbf{x}),$$

$$({ϵ}_n-{ϵ}_m^{\ast })(z_m|z_n)=0,$$

$${\left[\omega -vk\right]}^2={\Omega }^2(k)=c^2{k}^2+\frac{{\hslash }^2{k}^4}{4m^2}=c^2{k}^2\left[1+\frac{(k\xi {)}^2)}{4}\right]$$

Full text

I Introduction

Hawking radiation is one of the most intriguing results of theoretical Physics; using a semiclassical model in which fields are quantized on top of a classical gravitational background, Hawking predicted the spontaneous emission of radiation by the event horizon of a black hole (BH) Hawking (1974, 1975). Within a similar scheme, Corley and Jacobson Corley and Jacobson (1999) showed that a bosonic field with superluminal dispersion relation in a metric with two horizons can give rise to a dynamical instability, the so-called black-hole laser (BHL) effect. The problem is that the observation of such phenomena seems unlikely in the near future due to the small effective temperature of emission, $T_{H}\simeq 62M_{\odot}/M~{}\text{nK}$, with $M_{\odot}$ the mass of the Sun and $M$ the mass of the black-hole. For instance, the microwave background temperature is $2.7$ K, well above the Hawking temperature $T_{H}$.

An alternative way to study these effects was suggested by Unruh Unruh (1981), who proved that a subsonic-supersonic interface in a quantum fluid is the acoustic analog of an event horizon in a BH. This pioneering work opened the door to the study of gravitational problems in the laboratory and since then, many analog setups have been proposed in systems such as different as Fermi gases Giovanazzi (2005), ion rings Horstmann et al. (2010), polaritons Nguyen et al. (2015) or, in a classical context, surface waves in a water tank Euvé et al. (2016).

Of particular interest are the analogues implemented in Bose-Einstein condensates (BEC), first suggested by Garay et al. Garay et al. (2000). The main advantages of this kind of setup are the low temperature, the relative ease of handling, and the deep understanding of the quantum excitations. The analogue of the Hawking radiation in this system is the spontaneous emission of entangled phonons by the acoustic horizon into the subsonic and supersonic regions Leonhardt et al. (2003a, b); Balbinot et al. (2008); Carusotto et al. (2008); Recati et al. (2009); Macher and Parentani (2009a); Coutant and Parentani (2010); de Nova et al. (2014a); Finazzi and Carusotto (2014); Busch and Parentani (2014); de Nova et al. (2014b, 2015a); Michel et al. (2016a, b), similar to the particle-antiparticle creation at the event horizon of a BH. In addition, a flowing condensate presenting a finite-size supersonic region (giving rise to a pair of acoustic horizons) provides the analog of a black-hole laser Barceló et al. (2006); Jain et al. (2007); Coutant and Parentani (2010); Finazzi and Parentani (2010); Michel and Parentani (2013, 2015); de Nova et al. (2016).

Regarding the experimental side, the first acoustic horizon in a BEC was produced by the Technion group Lahav et al. (2010) with the help of a sharp negative potential created by a laser that locally accelerate the atoms. Within this kind of setup, the first observation of the BHL effect was reported Steinhauer (2014), although there is still some discussion in the community about the interpretation of the experimental results Tettamanti et al. (2016); Wang et al. (2016); Steinhauer and de Nova (2017). Recently, the same group provided the first experimental evidence of the emission of Hawking radiation by measuring the entanglement of the emitted phonons Steinhauer (2016).

Most of the theoretical works present in the BEC analog literature deal with an extremely idealized model, the so-called flat-profile configuration, in which the background condensate is homogeneous and the horizons are created through a very specific spatial dependence of the coupling constant and the external potential. Although this simple model is able to capture the essential features of Hawking radiation, it is quite unrealistic from an experimental point of view. More realistic models study the formation of acoustic BHs considering the flow of a condensate through a localized obstacle, modeled by a delta barrier Kamchatnov and Pavloff (2012), or an optical lattice de Nova et al. (2014b); the waterfall configuration described in Ref. Larré et al. (2012) is a theoretical model of the actual experimental setups of Refs. Lahav et al. (2010); Steinhauer (2016).

The goal of this Article is to extend the previous results and provide more realistic theoretical models also for analog BHLs. For that purpose, we prove that each BH configuration can be symmetrically extended to provide a BHL configuration. By applying this result to the waterfall and the delta-barrier configurations described above, we obtain two new different black-hole laser configurations that are created by using an attractive square well and a double delta-barrier, respectively. For these configurations, we compute the different families of non-linear stationary states that characterize the stability of the system as well as its long-time behavior. We note that, although stationary transport scenarios in a square well or a double delta-barrier have already been studied in the literature Leboeuf and Pavloff (2001); Zapata et al. (2011), to the best of our knowledge this is the first time that they have been explicitly proposed for modeling black-hole lasers.

Apart from the intrinsic interest of finding new models from a theoretical point of view, these configurations are also expected to be very useful in practice; in particular, the case where the supersonic region is created using an square attractive well can be regarded as a model for studying the experimental BHL of Ref. Steinhauer (2014).

The scheme of the paper is the following. In Sec. II we revisit the basic theory of gravitational analogues in BEC. The general relation between BH and BHL solutions is proven in Sec. III. In Secs. IV, V we study the different stationary states for BHL configurations with a square-well and a double delta-barrier, respectively. Conclusions are presented in Sec. VI. Appendix A is devoted to introduce the different elliptic functions used in this work while Appendices B, C are devoted to the technical details of the calculations presented in the main text.

II Gravitational analogues in Bose-Einstein condensates

We first provide in this section a general introduction to Bose-Einstein condensates and gravitational analogues. For more details, see for instance Refs. Macher and Parentani (2009b); Recati et al. (2009); Zapata et al. (2011); de Nova et al. (2015b).

II.1 Effective one-dimensional configurations

We begin by reviewing how to reach an effective one-dimensional (1D) configuration, the so-called 1D mean-field regime Leboeuf and Pavloff (2001); Menotti and Stringari (2002). For that purpose, we consider a 3D gas of $N$ bosons of mass $m$ near $T=0$ (more precisely, $T\ll T_{c}$, with $T_{c}$ the critical temperature of the condensate), described by the second-quantization Hamiltonian Fetter and Walecka (2003); Dickhoff and Van Neck (2005)

[TABLE]

where $V_{\rm{ext}}(\mathbf{x},t)$ is the external potential and the interaction between atoms is taken into account for low momentum by a short-range potential, with $g_{\rm{3D}}=4\pi\hbar^{2}a_{s}/m$ the corresponding coupling constant and $a_{s}$ the $s$-wave scattering length Pitaevskii and Stringari (2016); Pethick and Smith (2008).

In order to obtain an effective 1D configuration along the $x$-axis, we consider a total external potential of the form $V_{\rm{ext}}(\mathbf{x},t)=V(x,t)+V_{\rm{tr}}(y,z)$, where $V(x,t)$ only depends on the $x$ coordinate while $V_{\rm{tr}}(y,z)=\frac{1}{2}m\omega_{\rm{tr}}^{2}\rho^{2}$ ($\rho=\sqrt{y^{2}+z^{2}}$ being the radial distance to the $x$-axis) represents a transverse harmonic trap, very usual in experimental setups. If the non-linear interacting term is sufficiently small, we can treat it perturbatively and assume that the transverse motion is frozen to the corresponding harmonic oscillator ground state, and hence use the following approximation for the field operator

[TABLE]

with $a_{\rm{tr}}$ being the transverse harmonic oscillator length. After integrating over the transverse degrees of freedom, we arrive at the following 1D effective Hamiltonian:

[TABLE]

where we have absorbed the resulting zero-point energy $\hbar\omega_{\rm{tr}}$ of the harmonic oscillator and the effective 1D constant coupling is $g_{\rm{1D}}=2\hbar\omega_{\rm{tr}}a_{s}$. More specifically, the condition for the approximation of Eq. (2) to be valid is that the non-linear interacting term is small compared to the transverse confinement energy scale, $g_{1D}n_{1D}(x)\ll\hbar\omega_{\rm{tr}}$, which can be simply put as $n_{1D}(x)a_{s}\ll 1$, with $n_{1D}(x)$ the 1D-density.

In the same fashion, the 3D canonical commutation rules for the field operator

[TABLE]

are reduced to the 1D version

[TABLE]

As we will only deal with 1D configurations, in the following we omit everywhere the 1D index.

II.2 Gross-Pitaevskii and Bogoliubov-de Gennes equations

Using the Hamiltonian (3) and the corresponding canonical commutation rules (5), we write the equation of motion for the field operator $\hat{\psi}(x)$ in the Heisenberg picture:

[TABLE]

Since there is a condensate, we can perform a mean-field approximation

[TABLE]

with $\psi(x,t)$ the Gross-Pitaevskii (GP) wave function Pitaevskii and Stringari (2016) describing the condensate and $\hat{\varphi}(x,t)$ representing the quantum fluctuations of the field operator. The time evolution of the GP wave function is described by the time-dependent GP equation, a non-linear Schrödinger equation of the form:

[TABLE]

Assuming that the depletion cloud (i.e., the cloud formed by the atoms outside the condensate) is negligible, $\psi(x,t)$ is normalized to the total number of particles:

[TABLE]

The conservation of the norm of the GP wave function is guaranteed by the same relation as in the usual linear Schrodinger equation:

[TABLE]

with $J(x,t)$ the current.

It is quite instructive to rewrite these equations in terms of the amplitude and phase of the wave function, $\psi(x,t)=A(x,t)e^{i\phi(x,t)}$,

[TABLE]

where $J(x,t)=n(x,t)v(x,t)$ is the current and

[TABLE]

are the mean-field density and flow velocity, respectively. Interestingly, the first line of Eq. (11) is the equivalent of the continuity equation for a hydrodynamical fluid. On the other hand, taking spatial derivative in the second line gives:

[TABLE]

which can be regarded as the analog of the Euler equation for the velocity of a potential flow since the pressure of a uniform condensate at equilibrium is $P=\frac{gn^{2}}{2}$. The only difference is the rightmost term, which is a genuine quantum feature as it contains $\hbar$ and it is often called the quantum pressure term. However, in the hydrodynamic regime, where the density of the condensate varies on a large scale compared to the other terms, one can neglect the contribution of the quantum pressure and recover the same equations as for an ideal potential fluid flow; this is the key point of the gravitational analogy since the original analogy was precisely established for ideal potential fluid flows Unruh (1981) (see also discussion in the next subsection).

On the other hand, to lowest order in the quantum fluctuations of the field operator, $\hat{\varphi}(x,t)$, one finds from Eq. (6) that:

[TABLE]

which are known as the Bogoliubov-de Gennes (BdG) equations.

For time-independent potentials, $V(x,t)=V(x)$, we can look for particular solutions of the form

[TABLE]

that are of special interest as they describe stationary configurations. In particular, the stationary wave function $\psi_{0}(x)$ obeys the time-independent GP equation:

[TABLE]

Note that the above equation is a non-linear eigenvalue problem. The presented amplitude-phase decomposition greatly simplifies the stationary problem as the continuity equation is reduced to

[TABLE]

so the current $J(x)=n(x)v(x)=J$ is constant. Using this fact, we can rewrite the equation for the amplitude as a purely real second-order differential equation:

[TABLE]

with ′ the spatial derivative. The phase is simply obtained from the relation:

[TABLE]

Note that, if $J\neq 0$, neither the amplitude or the flow velocity vanish.

For a fixed value of the number of particles, there can be several different solutions for the GP equation (20). The true ground state of the system is that minimizing the grand-canonical energy, $K=E-\mu N$, with $E$ the energy of the state (i.e., the expectation value of the Hamiltonian evaluated for the GP wave function) and $N$ the total number of particles. Indeed, by rewriting the expression for $K$ as a functional for the GP wave function:

[TABLE]

it can be seen that Eq. (20) is precisely the condition for $\psi_{0}(x)$ to be an extreme of $K$. In that case, $K$ takes the simple form:

[TABLE]

Solutions of Eq. (20) which are not a local minimum of $K$ are energetically unstable as any perturbation would induce the system to decay to a lower energy state. In physical terms, we can understood the minimization of $K$ as the minimization of the expectation value of Hamiltonian $H$ with the constraint of fixed total number of particles $N$, with the chemical potential $\mu$ playing the role of Lagrange multiplier.

With respect to the quantum fluctuations, $\hat{\varphi}(x,t)$, Eq. (II.2) is now a stationary problem of the form:

[TABLE]

As it is a linear equation, we can expand the field operator in terms of eigenmodes:

[TABLE]

where the spinors $z_{n}$ satisfy the time-independent BdG equations:

[TABLE]

Due to the structure of the equations, the conjugate $\bar{z}_{n}$ is also a mode with energy $-\epsilon^{*}_{n}$.

An interesting property of the eigenvalue problem (34) is that it is non-Hermitian and thus it can yield complex eigenvalues. In particular, eigenvalues with positive imaginary part correspond to dynamical instabilities, i.e., exponentially growing modes: the presence of such dynamical instabilities in a finite region of a condensate flow are the origin of the black-hole laser effect, discussed in the next section.

Moreover, there is a Klein-Gordon type scalar product associated to the BdG eigenvalue problem, given by:

[TABLE]

with $\sigma_{z}=\rm{diag}(1,-1)$ the corresponding Pauli matrix. Note that this scalar product is not positive definite so the norm of a given solution $z_{n}$, defined as $(z_{n}|z_{n})$, can be positive, negative or zero. In fact, the norm of the conjugate $\bar{z}_{n}$ has the opposite sign to that of $z_{n}$, $(\bar{z}_{n}|\bar{z}_{n})=-(z_{n}|z_{n})$.

The utility of this scalar product is that, as usual, two modes $z_{n},z_{m}$ with different eigenvalues $\epsilon_{n},\epsilon_{m}$ are orthogonal, as seen from the relation

[TABLE]

from which it also follows that modes with complex frequency have zero norm.

II.3 Analog configurations

Gravitational analogues in BEC appear when considering stationary condensate flows. We note that, although for illustrative purposes we restrict here to 1D configurations, the following discussion can be straightforwardly adapted for general 3D stationary flows. First, we analyze 1D homogeneous stationary flows, characterized by GP plane waves of the form $\psi_{0}(x)=\sqrt{n}e^{iqx+\phi_{0}}$, with $n$ the density of the condensate, $q$ its momentum and $\phi_{0}$ some phase. After removing the phase of the condensate from the field operator, $\hat{\varphi}(x,t)\rightarrow e^{iqx+\phi_{0}}\hat{\varphi}(x,t)$, it is straightforward to show that the eigenmodes of the BdG equations (34) are plane waves with wave vector $k$ and frequency $\omega$, giving rise to the following dispersion relation:

[TABLE]

with $c=\sqrt{gn/m}$ the sound velocity, $v=\hbar q/m$ the constant flow velocity, $\xi\equiv\hbar/mc$ the so-called healing length and $\Omega$ the comoving frequency. The above dispersion relation gives four different wavevectors for a given value of the frequency. In fact, Eq. (37) is just the usual Bogoliubov dispersion relation for phonons in a condensate at rest, $\Omega(k)$, shifted by the Doppler effect due to the fluid velocity $v$. For convention, we take the flow velocity and comoving frequency as positive ($v,\Omega>0$) throughout this work. In this way, the flow is supersonic when $v>c$ and subsonic when $v<c$.

The dispersion relation for subsonic (supersonic) flows is schematically represented in left (right) panel of Fig. 1, where the blue (red) curves represents the sign $+(-)$ branches of the dispersion relation, $\omega_{\pm}(k)=vk\pm\Omega(k)$, and also positive (negative) normalization according to the scalar product of Eq. (35). Indeed, the $-$ branch is just the dispersion relation of the conjugate modes of the $+$ branch, $\omega_{-}(k)=-\omega_{+}(-k)$ [see Eq. (34) and related discussion]. For subsonic flows, for a given real frequency, there are only two propagating modes (i.e., modes with purely real wavevector) and the other two solutions have complex wave vector. On the other hand, for supersonic flows, in the window $-\omega_{\rm max}<\omega<\omega_{\rm max}$ all the four modes are propagating, where the threshold frequency $\omega_{\rm max}$ is $\omega_{\rm max}=\max_{k}\omega_{-}(k)$ and is marked by a horizontal dashed line in right Fig. 1. Outside this window, we recover essentially the same structure of subsonic flows and only two modes are propagating. The presence of negative energy modes for $-k_{0}<k<0$ in the $+$ branch of the supersonic dispersion relation, with $\hbar k_{0}=2m\sqrt{v^{2}-c^{2}}$, arises due the energetic instability of supersonic flows, as first argued by Landau. As a result, the introduction of a time-independent perturbation in a supersonic flow gives rise to the emission of Bogoliubov-Čerenkov radiation Carusotto et al. (2006), characterized by the wave vector $k_{0}$.

The previous magnitudes can be extended to non-homogeneous configurations by taking $c(x)\equiv\sqrt{gn(x)/m}$ and $v(x)$ as defined in Eq. (12). In an similar way, we say that the flow is subsonic where $v(x)<c(x)$ and supersonic where $v(x)>c(x)$. It is precisely in this context where the gravitational analogy with astrophysical black holes arises. For that purpose, we rewrite Eq. (II.2) in terms of the relative quantum fluctuations, $\hat{\varphi}(x,t)\equiv\psi_{0}(x)\hat{\chi}(x,t)$,

[TABLE]

with $D_{t}=\partial_{t}+v(x)\partial_{x}$ the comoving derivative. Gathering this equation with its complex conjugate and defining the hermitic fields

[TABLE]

gives a pair of equations of the form:

[TABLE]

The above fields are related to physical magnitudes as $\hat{\rho}(x,t)=2n(x)\hat{\chi}_{+}(x,t)$ and $\hat{\phi}(x,t)=\hat{\chi}_{-}(x,t)$, with $\hat{\rho},\hat{\phi}$ the density and phase fluctuations, respectively. Note that the first line of the above equation results from linearizing the continuity equation while the second line results from linearizing the equation for the phase [see Eq. (11)].

Now, if we assume that the background condensate varies on a sufficiently large scale, in the long-wavelength limit we can neglect the contribution of $T(x)$ at the r.h.s. of the second line of Eq. (40), which precisely amounts to work in the hydrodynamic regime where all the contributions arising from the quantum pressure are neglected. In this approximation, we can write the equation for the phase fluctuations as:

[TABLE]

which can be rewritten as the relativistic Klein-Gordon equation for a massless scalar field $\hat{\phi}$ on a metric $g_{\mu\nu}$,

[TABLE]

with the effective stationary metric $g_{\mu\nu}$ given by

[TABLE]

Thus, the points where $c(x)=v(x)$ are the horizons of the acoustic metric $g_{\mu\nu}(x)$, analog to astrophysical event horizons. Using a simple physical picture, for acoustic phonons (long-wavelength excitations) the dispersion relation (37) has the form $\omega_{\pm}(k)\simeq(v\pm c)k$, so they are dragged away by a supersonic flow and hence trapped in the supersonic side of an acoustic horizon in the same way as light is trapped inside the event horizon of a black hole. Nevertheless, although the above derivation was done for phonons, the analog of the Hawking effect still holds when considering the complete superluminal dispersion relation in a black hole (BH) configuration Macher and Parentani (2009b); Recati et al. (2009), where modes with a sufficiently large wave vector in the supersonic region can travel upstream and escape unlike in gravitational black holes, where nothing escapes. A BH configuration is defined as that with two asymptotic homogeneous regions, one subsonic and one supersonic, with flow traveling from subsonic to supersonic, while if the flow goes from supersonic to subsonic, we have a white-hole (WH) configuration, the time reversal of a BH (which just amounts to take the complex conjugate of the GP wave function). Invoking the continuity of the wave function, a BH configuration implies that, at least, one acoustic horizon is formed, that is, a point where $v(x)=c(x)$.

In the same fashion, a configuration with two asymptotic homogeneous subsonic regions and displaying a pair of acoustic horizons (corresponding to a black and a white hole) is the analog of a black-hole laser (BHL). Specifically, the BHL effect in this setup is characterized by the appearance of dynamical instabilities in the BdG spectrum. As supersonic flows are energetically unstable, one can expect this instability to occur for sufficiently large supersonic regions between the two horizons. A more physical insight of the process can be given using a semiclassical picture Finazzi and Parentani (2010): negative energy radiation emitted at the BH impacts at the WH and, due to the superluminal dispersion relation, some of the reflected modes are able to travel upstream and hit the BH again, stimulating further emission. This process originates a self-amplifying emission that gives rise to a dynamical instability in the flow.

II.4 Solutions of the homogeneous Gross-Pitaevskii equation

We finally end this section by reviewing the different stationary solutions of the homogeneous GP equation as they are the building blocks of most of the theoretical analog models due to their analytical tractability, serving also as a basis for the calculations presented in this work. In homogeneous problems, the external potential $V(x)$ is constant and it can be reabsorbed into the definition of the chemical potential. In that case, the resulting equation for the amplitude, Eq. (22), is analog to the Newtonian equation of motion of a classical particle in a potential with the role of position and time played here by the amplitude of the wave function and the spatial coordinate $x$, respectively. Then, it can be integrated to obtain

[TABLE]

The quantity $E_{A}$ is the “energy” of the classical particle and we refer to it as the amplitude energy while $W(A)$ is the corresponding amplitude potential.

As a first step, we study the equilibrium points of $W(A)$ as they give the homogeneous plane wave solutions, $A(x)=A$, which can be obtained from the zeros of Eq. (22):

[TABLE]

This polynomial equation for the density only has (two) real positive roots whenever:

[TABLE]

In the rest of the work, we will assume that condition (46) is fulfilled and denote the largest root as $n=n_{0}$; the associated flow velocity is constant and equal to $v_{0}=J/n_{0}$. In order to simplify the calculations, we rescale the wave function as $\psi_{0}(x)\rightarrow\sqrt{n_{0}}\psi_{0}(x)$ so it become dimensionless and take units such that $\hbar=m=c_{0}=1$, with $c_{0}=\sqrt{gn_{0}/m}$ the sound velocity associated to the density $n_{0}$. Length, time and energy are measured in units of $\xi_{0}=\hbar/mc_{0}$, $t_{0}=\xi_{0}/c_{0}$ and $E_{0}=mc^{2}_{0}$, respectively. Also, we will refer to $v_{0}$ as simply $v$ for simplicity.

In this system of units, the amplitude of the homogeneous solution with density $n_{0}$ is just $A=A_{0}=1$ and the associated current simply reads $J=v$, while the chemical potential is $\mu=1+v^{2}/2$. Indeed, $v$ also represents now the value of the Mach number of the flow [that is, the dimensionless ratio between the flow velocity and the speed of sound, $v(x)/c(x)$].

With the help of these considerations, we rewrite Eq. (45) as

[TABLE]

from which we immediately obtain the density of the other homogeneous solution:

[TABLE]

By construction, $n_{p}<n_{0}=1$, which implies that $v<1$ and hence the homogeneous solution $n=n_{0}=1$ is necessarily subsonic (note that the limit value $v=1$ corresponds to the degenerate case $A_{p}=A_{0}=1$). The flow velocity of the solution $n=n_{p}$, $v_{p}$, is obtained from the conserved current $n_{p}v_{p}=J=v$. As in these units the sound speed is just the square root of the density, $c(x)=\sqrt{n(x)}=A(x)$, the Mach number of the solution $A=A_{p}$ is $v_{p}/c_{p}=v/n_{p}^{3/2}$. By defining $z\equiv 8/v^{2}$ and observing that the function $f(z)=z^{2/3}-1-\sqrt{1+z}$ increases monotonically for $z>0$, we conclude that the solution $A=A_{p}$ corresponds to a supersonic flow since $f(z)$ only has one zero at $z=8$ ($v=1$).

We also study the non-homogeneous solutions of Eq. (44). For that purpose, we represent the amplitude potential $W(A)$ in Fig. 2. The local minimum corresponds to the homogeneous supersonic solution $A_{p}=\sqrt{n_{p}}$ (which means that it is an stable fixed point) and the local maximum to the homogeneous subsonic solution $A_{0}=\sqrt{n_{0}}=1$ (which means that it is an unstable fixed point).

Rewriting Eq. (44) in terms of the density gives the simplified equation:

[TABLE]

with the densities $n_{i},~{}i=1,2,3$ computed from the zeros of the equation $W(A)=E_{A}$, equivalent to obtaining the roots of the following polynomial equation in terms of the density

[TABLE]

Several cases can be distinguished depending on the value of the amplitude energy $E_{A}$. First, for $W(A_{p})<E_{A}<W(1)$, the three roots of Eq. (50) are real and we order them such that $0<n_{1}<n_{2}<n_{3}$. The case $n_{1}<n(x)<n_{2}$ corresponds to the oscillating solution represented by the the closed blue line of right Fig. 2. By integrating Eq. (49), we find that:

[TABLE]

with $x_{0}$ some integration constant arising from the translational invariance of the problem. The solution of the previous indefinite integral is given in term of elliptic functions. The resulting phase of the wave function is computed from Eq. (23), obtaining:

[TABLE]

where:

[TABLE]

with $\phi(0)$ some global phase. We refer the reader to Appendix A for the precise definition of the different elliptic functions used along this work.

The case $n(x)>n_{3}$ for the same value of $E_{A}$ corresponds to a solution that grows indefinitely, see blue curve for $A>1$ in right plot of Fig. 2). Moreover, for high values of $n$, $n^{\prime}\propto n^{3/2}$ and then the solution blows up at some finite value $x_{bu}$ as $n(x)\sim(x-x_{bu})^{-2}$. The same reasoning holds for $E_{A}<W(A_{p})$ (green dashed-dotted line of Fig. 2) or $E_{A}>W(1)$ (red dashed line of Fig. 2). These exploding solutions are not relevant for the present work so we ignore them in the following.

Finally, we consider the degenerate cases $E_{A}=W(A_{p})$ and $E_{A}=W(A_{0})$. For $E_{A}=W(A_{p})$, $n_{1}=n_{2}=n_{p}$. One possible solution corresponds to the stable fixed point of the homogeneous supersonic solution $n(x)=n_{p}$, described by the plane wave:

[TABLE]

The other possible solution corresponds to $n(x)>n_{3}$, which blows up in a similar way to the exploding solutions described above.

The other degenerate case is $E_{A}=W(1)$, where the roots satisfy $n_{2}=n_{3}=1$ and $n_{1}=v^{2}$. For $n(x)=1$, the solution is the subsonic plane wave

[TABLE]

For $n(x)\neq 1$, we get from Eq. (51)

[TABLE]

For $n(x)<1$, we obtain:

[TABLE]

which is of the same form of Eq. (53) after taking into account that $\text{sn}(u,1)=\tanh(u)$. The phase of the wave function can be obtained analytically in a simple form and we can write the total wave function as:

[TABLE]

being $\phi_{0}$ some constant phase. This solution represents a soliton with zero velocity Pethick and Smith (2008).

On the other hand, taking $n>1$ in Eq. (56) gives the so-called shadow soliton solution Michel and Parentani (2013):

[TABLE]

Although this solution also blows up at a finite value of $x$, it is quite relevant for the computation of stationary states in BHL configurations; see Secs. IV, V.

III General relation between black holes and black-hole lasers

After introducing the basic concepts and techniques of gravitational analogues in BEC, we proceed to prove one of the central results of this work, which states that every compact BH solution can be used to produce a BHL configuration with an arbitrary large homogeneous supersonic region, explicitly showing the mechanism to construct such BHL configuration. We define a compact BH solution as that in which a homogeneous supersonic flow is reached at a finite point, $x=x_{H}$. Indeed, this is the situation of all the BH configurations usually appearing in the literature Larré et al. (2012).

The proof is straightforward. Consider a compact BH configuration, which satisfies a time-independent GP equation of the form:

[TABLE]

For simplicity, we consider that the BH is produced only with the help of an external potential but the generalization to situations in which the coupling constant (like the flat-profile configuration) or the mass are space-dependent is trivial. By definition of compact BH configuration, $V_{C}(x>x_{H})=V_{sp}$ is homogeneous and the GP wave function is of the form:

[TABLE]

with $A_{sp},q_{sp}$ the supersonic amplitude and momentum and $\Psi_{C}(x)$ the part of the wave function that describes the subsonic-supersonic transition. Without loss of generality, we choose the origin of coordinates such that $x_{H}=-X/2$, with $X>0$.

The idea for obtaining a BHL configuration is to replicate the same structure of the potential and the GP wave function for $x>0$. This can be done by taking the spatial and time reverse of the wave function and the potential. Explicitly, we consider the GP wave function

[TABLE]

which satisfies the following GP equation

[TABLE]

where the potential $V_{BHL}(x)$ is given by:

[TABLE]

The wave function of Eq. (62) describes a BHL configuration with a homogeneous supersonic flow in a region of size $X$. Indeed, since $X$ is not fixed, we can construct a supersonic region of arbitrary length with this solution. Therefore, we conclude that every BH solution can be extended to produce a BHL configuration.

Now, with the help of the previous result, we study the BHL configurations arising from some well-known BH configurations, shown in the lower and upper row of Fig. 3, respectively. One of the most considered BH configurations is the flat-profile configuration (upper left of Fig. 3), in which the GP wave function is a global plane wave so the condensate density $n$ and flow velocity $v$ are homogeneous in all of the space. In order to fulfill the homogeneity condition, the 1D coupling strength $g\left(x\right)$ and the external potential $V\left(x\right)$ must satisfy that $g(x)n+V(x)$ is constant. In particular, in order to construct a BH solution, $g(x)$ is chosen to be a step function with a downstream value $g(x)=g_{2}$ such that the resulting flow is supersonic. Although the experimental implementation of this configuration is extremely challenging due to the required high precision in the control of both the external potential and the local coupling constant, it is considered in many theoretical works Balbinot et al. (2008); Carusotto et al. (2008); Recati et al. (2009) because of its analytical simplicity.

More realistic configurations are displayed in central and right panels of the upper row of Fig. 3, corresponding to the delta-barrier and waterfall configurations, respectively; the point is that these configurations only require simple external potentials that are achievable with the use of standard experimental tools as blue-detuned (for repulsive potentials) or red-detuned (for attractive potentials) lasers. For instance, the delta-barrier configuration models the BH arising from the flow of a condensate through a localized obstacle Kamchatnov and Pavloff (2012), represented by a repulsive delta potential; by Galilean invariance, this configuration is similar to launching the obstacle against the condensate. On the other hand, the waterfall configuration uses an attractive step potential to accelerate the flow and create a supersonic current. In fact, this model provides a realistic description of the actual setups of Refs. Lahav et al. (2010); Steinhauer (2016), in which a negative step potential created with the help of a laser is swept along a trapped condensate, finding a good agreement with the experimental data Steinhauer (2016).

Regarding the BHL side, the associated flat-profile BHL of lower left Fig. 3 has been already considered in the literature in both analytical and numerical studies Michel and Parentani (2013, 2015); de Nova et al. (2016). Specifically, Ref. Michel and Parentani (2013) provided a detailed study not only of the associated dynamical instabilities of the flow but also of the different non-linear stationary solutions existing for fixed chemical potential and current as they describe the potential quasi-stationary states of the system for long times, once the dynamical instability has already grown up. In fact, it was shown the appearance of each dynamical instability is associated to the appearance of a stationary solution with lower grand-canonical energy $K$ than the initial homogeneous plane-wave solution.

On the other hand, Refs. Michel and Parentani (2015); de Nova et al. (2016) extended the previous analytical work with numerical simulations of the time evolution of the initially unstable homogeneous solution in order to study the non-linear saturation of the instability. In particular, Ref. de Nova et al. (2016) found that the system only presents two kind of asymptotic behaviors: it either reaches the GP ground state or a regime of continuous emission of solitons (CES) in which the system emits trains of solitons with perfect periodicity, providing in this way the closest analog of an actual optical laser.

Following these results, in order to provide more realistic models for black-hole lasing, we study in the rest of the work the BHL configurations associated to the delta-barrier and the waterfall configurations, depicted in the lower center and right panel of Fig. 3, respectively. Specifically, due to energy and particle number conservation Michel and Parentani (2013), we only aim at the GP stationary solutions asymptotically matching at $\pm\infty$ the corresponding subsonic plane wave solution. Several reasons motivate this choice: first, solutions with lower grand-canonical energy and continuously connected to the initial BHL solution are expected to also characterize the appearance of dynamically unstable modes Michel and Parentani (2013). In addition, there should not be substantial differences in the linear regime with respect to the usual flat-profile BHL. Moreover, thinking in a realistic implementation, the computed stationary solutions should still govern the late time dynamics once the system enters into the non-linear regime Michel and Parentani (2013); de Nova et al. (2016) regardless the specific mechanism giving rise to the initial growth of the instability at short times Tettamanti et al. (2016); Steinhauer and de Nova (2017); Wang et al. (2016).

Although the explained protocol to create BHL configurations applies for any arbitrary compact BH configuration, for illustrative purposes we focus on these two particular cases as they provide simple analytical models of realistic experimental scenarios and extend well-known models in the literature. Of particular interest is the BHL resulting from the waterfall configuration as it corresponds to an attractive potential well and hence it is expected to capture the essential features of the actual BHL configuration of the experiment of Ref. Steinhauer (2014), in which the laser cavity is created by sweeping along the condensate the effective potential well arising from the combination of the background trap and the negative step potential.

In order to simplify the notation and match the results of Sec. II.4, we set units in the rest of this work such that $\hbar=m=c_{0}=1$, where $c_{0}$ is the asymptotic subsonic speed of sound. Once in these units, it is easy to check that for the BHL configurations considered in this work, the problem is completely determined by only two parameters: the asymptotic subsonic flow velocity $v$ (which is also the subsonic Mach number in these units) and the size of the supersonic region $X$, since the amplitudes of the different potentials are functions of $v$. This contrasts to the case of the flat-profile BHL configuration, where there are three degrees of freedom, $v,X,c_{2}$, with $c_{2}$ the supersonic sound speed de Nova et al. (2016).

For both configurations, we first present the general structure of the problem and then describe the main features and the conditions of existence for the different families of stationary solutions; technical details of the computations are given in Appendices B, C, to which the interested reader is referred.

IV Attractive square well

IV.1 General structure

As a first step, we present the BH solution corresponding to the waterfall configuration, where the external potential is given by:

[TABLE]

where $\theta(x)$ is the step function and $x_{H}$ the point where the step is placed. The corresponding GP wave function is given by:

[TABLE]

with $\phi_{0}$ some phase to make the wave function continuous at $x=x_{H}$. Following the procedure described in Sec. III, the associated BHL configuration is created by using an attractive square well potential of size $X$,

[TABLE]

We note that the different stationary solutions of a condensate flowing through an attractive square well were first addressed in Ref. Leboeuf and Pavloff (2001). Here we restrict to the specific case where $V_{0}$ is given by Eq. (65), so a BHL solution as that of lower right Fig. 3 exists, described by the GP wave function:

[TABLE]

In order to find the remaining stationary solutions that asymptotically match the subsonic plane-wave solution $e^{ivx}$ on both sides, we use the phase-amplitude decomposition and consider Eq. (44). The asymptotic boundary conditions fix the current to $J=v$ and the amplitude energy outside the well to $E_{A}=E_{1}=\frac{1}{2}+v^{2}$. Accordingly, two different regions can be distinguished: region 1 corresponds to the exterior of the square well, $|x|>X/2$, while region 2 corresponds to its interior, $|x|<X/2$. The equation for the amplitude reads in each region as:

[TABLE]

where $W_{i}(A),E_{i}$ are the amplitude potential and the conserved amplitude energy for the $i=1,2$ regions, with $\mu_{1}=\mu=1+v^{2}/2$ and $\mu_{2}=\mu_{1}+V_{0}=v^{2}+\frac{1}{2v^{2}}$. Invoking the continuity of the wave function and its derivative, we find the matching condition at both edges

[TABLE]

As $E_{1}$ is fixed, we only need to find the possible values of $E_{2}$ which make that the GP wave function satisfy both Eq. (69) and Eq. (70).

The situation is schematically depicted in left Fig. 4: outside the well, the orbits follow the dashed black line as they must asymptotically match the subsonic plane wave on both sides, while the other curves represent possible solutions inside the well.

Thus, there are only three possible solutions outside the well: shadow solitons, with amplitude $A(x)>1$; regular solitons, with $A(x)<1$; and the homogeneous subsonic plane wave, $A(x)=1$. On the other hand, as the amplitude at $x=\pm\frac{X}{2}$ must be the same, the only possible solution inside the well is a cnoidal wave as that of Eq. (53), where $n_{i},~{}i=1,2,3$, are now the roots of the equation

[TABLE]

For convention, we choose the wave function such that it is real at $x=0$, $\phi(0)=0$.

The matching of the cnoidal wave (53) at $x=\pm\frac{X}{2}$ gives two equations:

[TABLE]

Since $n_{1},n_{2},n_{3},n_{W}$ are functions of $E_{2}$, the above system gives two conditions for two variables, $\alpha$ and $E_{2}$. Due to the periodicity of the elliptic functions, the possible solutions for a given length $X$ are discretized by an index $m=0,1,2\ldots$ representing the number of complete periods inside the well. As a result, we only need to compute the corresponding values of $E_{2},\alpha$, labeled as $E^{m}_{2},\alpha_{m}$, and the associated parameters $n^{m}_{W},n^{m}_{i},\nu_{m}$ [see Eqs. (53), (70)] to determine the wave function; the details of this calculation are given in Appendix B.

We now discuss the different families of stationary solutions depending on the three possible cases for the wave function outside the well.

IV.2 Homogeneous plane wave

In this case, the wave function outside the well is the homogeneous subsonic plane wave and it is the most simple one as it represents a limit of the other two families of solutions. This solution only appears for some critical values of the length $X=X^{H}_{m}$, $m=0,1,2\ldots$, and the corresponding GP wave function is:

[TABLE]

where the critical values $n^{H}_{1,2,3},X^{H}_{m}$ are given in Appendix B.1. Note that $m=0$ corresponds to $X=X^{H}_{0}=0$ and the trivial solution of a homogeneous subsonic plane wave in whole space as there is no attractive well.

The orbit in phase space inside the well associated to this family of solutions is depicted as a red solid line in left Fig. 4 while the density profile for $\psi^{H}_{m}(x)$, $m=1,2,3$, is shown in the central panels of the same figure.

IV.3 Asymptotic shadow solitons

The resulting GP wave function for this family of solutions reads:

[TABLE]

with $\phi_{0},x_{0}$ chosen such that the wave function and its derivative are continuous, $x_{0}$ satisfying $x_{0}-X/2<0$.

As the solutions outside the well are shadow solitons, $1\leq n^{m}_{W}<n^{SH}$, where $n^{SH}$ is the density obtained from the intersection between the dashed black line and the solid blue line in left Figure 4,

[TABLE]

The limit case $n^{m}_{W}=1$ corresponds to $E^{m}_{2}=E^{H}_{2}$ and $\Psi^{SH}_{m}(x)=\Psi^{H}_{m}(x)$. In particular, the function $\Psi^{SH}_{0}(x)$ is continuously connected to the homogeneous subsonic plane wave solution, and since it has the larger amplitude inside the well, it represents the true ground state of the system, as can be inferred from Equation (25). Hence, the BHL solution $\Psi^{BHL}(x)$ of Equation (68) should be dynamically unstable for any length $X>0$ of the well. However, since $\Psi^{BHL}(x)$ for $X=0$ corresponds to a perfect soliton, and it is not continuously connected to the ground state, the expected instability should be just the acceleration of the soliton Busch and Anglin (2000), rather than the growth of a lasing mode, at least for small cavity lengths $0<X<X^{C}_{0}$ in which there is no room for other solutions (see discussion in the next section). A detailed computation of the BdG spectrum should be carried out in order to confirm the previous hypothesis.

On the other hand, the upper limit $n^{m}_{W}=n^{SH}$ is a strict inequality, as it corresponds to $E^{m}_{2}=E^{SH}_{2}$, so the wave function is a soliton inside the well, giving rise to an infinite size $X$. Therefore, $E^{H}_{2}\leq E^{m}_{2}<E^{SH}_{2}$, and the $m$-th solution only exists for lengths:

[TABLE]

The density profile of $\psi^{SH}_{m}(x)$ for $m=0,1,2$ is represented in the rightmost panel of Fig. 4.

IV.4 Asymptotic solitons

When the solutions outside the well corresponds to solitons, things are more intriguing than in the previous cases. We mainly distinguish between symmetric solutions, where the density has even parity and the GP wave function satisfies $\psi_{0}(x)=\psi^{*}_{0}(-x)$, and asymmetric solutions, with no spatial symmetry.

IV.4.1 Symmetric solutions

For solutions with symmetric character, the solitons out of the well either both contain a minimum in the density (complete-soliton solutions) or not (incomplete-soliton solutions).

Incomplete-soliton solutions.

In this case, the GP wave functions reads:

[TABLE]

For incomplete-soliton solutions, $x_{0}$ satisfies $X/2-x_{0}>0$.

Since the solutions outside the well are solitons, $v^{2}\leq n^{m}_{W}\leq 1$. The lower boundary $n_{W}=v^{2}$ gives the same solution of Eq. (68), in which $n^{m}_{1}=n^{m}_{2}=n^{m}_{W}=v^{2}$ and $E^{m}_{2}=E^{C}_{2}\equiv E_{1}+v^{2}V_{0}=1+\frac{v^{4}}{2}$, existing for arbitrary length $X$. Hence, this family of solutions is continuously connected to the BHL configuration described by $\psi^{BHL}(x)$. A more detailed analysis of the limit $n^{m}_{W}\rightarrow v^{2}$ (see Appendix B.3) shows that these solutions appear at the critical lengths $X^{C}_{m}$, given by Eq. (125), and can be understood as small amplitude perturbations on top of the GP wave function of Eq. (68), described by the BdG equations. Precisely, in this limit, all the elliptic functions are reduced to trigonometric functions and the cnoidal wave to a regular sinusoidal wave with wavevector $k_{0}$, resulting from the corresponding BdG plane-wave solutions with zero frequency for the supersonic flow inside the well (see right Fig. 1 and related discussion). Indeed, as the dynamically unstable modes arising from the BHL effect are expected to first show as zero-frequency BdG modes Michel and Parentani (2013), the critical lengths $X^{C}_{m}$ should also signal the appearance of a new dynamical instability.

The upper limit, $n^{m}_{W}=1$ gives $\psi^{SOL}_{m}(x)=\psi^{H}_{m+1}(x)$, merging with the $\psi^{SH}_{m+1}(x)$ solutions. Then, the $m$-th solution of Eq. (81) only exists for

[TABLE]

The density profile for the incomplete-soliton solutions $m=0,1,2$ is represented in the upper row of Fig. 5.

Complete-soliton solutions.

The wave function is of the same form as that of Eq. (81) but with $x_{0}$ satisfying $X/2-x_{0}<0$. This family of solutions is also continuously connected to $\psi^{BHL}(x)$ and the limit $n^{m}_{W}=v^{2}$ gives the same length as for incomplete-soliton solutions, $X=X^{C}_{m}$. The upper limit $n^{m}_{W}=1$ gives $X=X^{H}_{m}$ and $\psi^{SOL}_{m}(x)=\psi^{H}_{m}(x)$, merging with the $\psi^{SH}_{m}(x)$ solutions. For $m=0$, $X^{C}_{0}>X^{H}_{0}=0$; however, for $m$ sufficiently large, $X^{C}_{m}<X^{H}_{m}$. Thus, the $m=0$ complete-soliton solution is limited to the range

[TABLE]

while for $m\geq 1$ we can only generally say that

[TABLE]

with $\tilde{X}_{m}\leq\min(X^{C}_{m},X^{H}_{m})$; check Appendix B.3 for the details.

The density profile for the complete-soliton solutions $m=0,1,2$ is represented in the central row of Fig. 5.

IV.4.2 Asymmetric solutions

These solutions are not symmetric with respect to the well and are characterized by one complete (incomplete) soliton at the left and one incomplete (complete) soliton at the right, one case corresponding to the spatial reverse of the density profile of the other. n particular, they contain an exact integer number of periods $m=1,2,3$ inside the well (see Eq. (126)), with $m=0$ giving the trivial solution of no well, $X=0$. The corresponding GP wave function reads:

[TABLE]

where $\delta x>0,\phi_{L},\phi_{R}$ are chosen such that the wave function and its derivative are continuous and $\pm$ corresponds to the case of complete-incomplete (incomplete-complete) solitons.

In the same fashion of the symmetric families, in the limit $n^{m}_{W}\rightarrow v^{2}$ the asymmetric solutions are continuously connected to $\psi^{BHL}(x)$, appearing at the critical lengths $X^{A}_{m}$ given by Eq. (129), while in the upper limit $n^{m}_{W}\rightarrow 1$ they converge to the homogeneous $\psi^{H}_{m}(x)$ solutions.

Hence, the $m$-th asymmetric solution is restricted for lengths

[TABLE]

The density profile for the asymmetric solutions $m=1,2,3$ is represented in the lower row of Fig. 5.

V Double delta-barrier

V.1 General structure

The potential corresponding to a single delta-barrier configuration is

[TABLE]

with $W(A)$ given by Eq. (44), $A_{p}$ the supersonic amplitude of Eq. (48) and $x_{H}$ the point where the barrier is placed. Such a delta potential originates a discontinuity in the derivative of $\psi_{0}(x)$ of the form $\psi^{\prime}_{0}(x_{H}^{+})-\psi^{\prime}_{0}(x_{H}^{-})=2Z\psi_{0}(x_{H})$. The resulting GP wave function is:

[TABLE]

The constants $x_{0},\phi_{0}$ are fixed by imposing the continuity of $\psi_{0}(x)$ and its derivative at $x=x_{H}$.

The associated BHL configuration is described by a cavity of length $X$ placed between two delta barriers,

[TABLE]

This configuration was studied in Ref. Zapata et al. (2011) in order to look for resonant BH configurations, which enhance the spontaneous Hawking signal de Nova et al. (2014a). Here, we focus only on looking for stationary BHL solutions, as in the previous section. In particular, by construction, a BHL solution as that of lower central Fig. 3 exists, described by the GP wave function:

[TABLE]

with $x_{0},\phi_{0}$ chosen such that the wave function and its derivative are continuous.

As in the computation for the square well, we distinguish two different regions: region 1 corresponds to the exterior of the cavity, $|x|>X/2$, while region 2 corresponds to its interior, $|x|<X/2$. In each region, we have that

[TABLE]

where $E_{i}$ is the conserved amplitude energy for the $i=1,2$ regions, $E_{1}=\frac{1}{2}+v^{2}$ fixed by the asymptotic subsonic behavior.

The wave function is continuous everywhere and the only effect of the two delta barriers is to introduce a discontinuity in the derivative of the wave function given by

[TABLE]

which, in terms of the amplitude, reads:

[TABLE]

Thus, we can understood the effect of the delta barriers as “instantaneously accelerating” the classical particle described by Eq. (87).

As a result of the above considerations, the only possible choice for the wave function outside the cavity is the soliton solution as the other solutions would monotonically increase. The same reasoning restrict even more the possibilities and the solitons must satisfy

[TABLE]

Inside the cavity, the solution corresponds to a cnoidal wave, with the phase chosen such that $\phi(0)=0$ for convention.

By joining Eqs. (87), (89) and (90), we find that the energy inside the cavity is related to the amplitude at the edges through:

[TABLE]

Hence, the value of $E_{2}$ is determined by the cnoidal waves arising from Eq. (87) that are compatible with Eq. (91). In particular, as $E_{p}\leq E_{2}<E_{1}$, with $E_{p}\equiv W(A_{p})$, the amplitudes at the edges of the cavity must be in the range

[TABLE]

where $A_{\inf}<A_{\sup}<1$ are obtained from the roots of the equation:

[TABLE]

As for the square well, the possible solutions are labeled by a discrete index $m=0,1,2\ldots$ representing the number of complete periods inside the cavity and the wave function is determined once $E^{m}_{2},\alpha_{m}$ are obtained; the details of this computation are provided in Appendix C.

In order to classify the different families of solutions, we first distinguish between symmetric and asymmetric solutions, in analogy to the discussion presented in Sec. IV.4.

V.2 Symmetric solutions

Symmetric solutions satisfy:

[TABLE]

and then the matching equations at the edges of the cavity read:

[TABLE]

with $n_{W}=A^{2}_{W}$. Since the first line gives the density at the edges $n_{W}$ as an implicit function of $E_{2}$, the second line is a similar matching condition as that of Eq. (72); the only difference is that now $n_{W}$ is a much more complicated function of $E_{2}$. In fact, there are two different solutions for $A_{W}<1$ for a given value of the amplitude energy $E_{2}$ in the range $E_{p}<E_{2}<E_{1}$, one satisfying $A_{p}<A_{W}<A_{\sup}$ and the other one $A_{\inf}<A_{W}<A_{p}$.

As the signs of the derivatives at the edges outside the cavity are fixed by Eq. (90), we classify the solutions according to the sign od the derivatives at the internal side of the cavity edges, $x=\pm X^{\mp}/2$. The limit values between the different solutions are obtained from

[TABLE]

which yields a similar equation to Eq. (93):

[TABLE]

This equation has two solutions for $A<1$, $A=A_{q}$ and, by construction, $A=A_{p}$, with $A_{q}<A_{p}$. The energies associated to these solutions are $E_{2}=E_{q}\equiv W(A_{q})$ and $E_{2}=E_{p}<E_{q}$.

Following the above considerations, we distinguish three families of solutions: $S+$, for $A_{p}<A_{W}<A_{\sup}$; $S-$, for $A_{q}<A_{W}<A_{p}$; and $SD$, for $A_{\inf}<A_{W}<A_{q}$.

V.2.1 $S+$ solutions

In this case:

[TABLE]

The corresponding wave function for this family of solutions is:

[TABLE]

We proceed to discuss the two limit values for the energy $E^{m}_{2}$. In analogy to the square well, $E^{m}_{2}=E_{p}$ gives $\psi^{S+}_{m}(x)=\psi^{BHL}(x)$, which exists for an arbitrary length. Following the results of Sec. IV.4, we analyze the limit $E^{m}_{2}\rightarrow E_{p}$ (see Appendix C.1.1 for the details), in which we find that these solutions appear at the critical lengths $X=X^{C}_{m}$, given by Eq. (134). As for the square well, they can be understood as small perturbations on top of the GP wave function $\psi^{BHL}(x)$, described by the zero-frequency BdG plane waves with wavevector $k_{0}$ in the supersonic region.

Reasoning in the same way, the critical lengths $X^{C}_{m}$ are expected to also describe the appearance of new dynamical instabilities. Indeed, the family of solutions $\psi^{S+}_{m}(x)$ has lower grand-canonical energy than $\psi^{BHL}(x)$; specifically, $\psi^{S+}_{0}(x)$ is the ground state of the system. Note that, in contrast to the square well, here the ground state is continuously connected to $\psi^{BHL}(x)$ and $\psi^{BHL}(x)$ should only be dynamically unstable for finite lengths $X>X^{C}_{0}>0$. Thus, we expect to find in this case a perfect correspondence between dynamical instabilities and stationary solutions with lower grand-canonical energy than the BHL solution (86), in the same lines of Ref. Michel and Parentani (2013).

The upper limit, $E^{m}_{2}=E_{1}$, corresponds to the soliton solution, that gives an infinite value for the cavity length $X$. Hence, the $m$-th $S+$ solution only exists for

[TABLE]

The density profile of $\psi^{S+}_{m}(x)$ for $m=0,1,2$ is represented in left Fig. 6.

V.2.2 $S-$ solutions

In this case:

[TABLE]

The wave function is given by the same formal expression of Eq. (99). This family of solutions is also continuously connected to $\psi^{BHL}(x)$ and the small amplitude limit near $E_{p}$ gives the same critical lengths $X=X^{C}_{m}$ as the $S+$ solutions.

The upper limit, $E^{m}_{2}=E_{q}$, only appears for discrete values of the length, $X=X^{q}_{m}$, given by Eq. (137); note that $X=X^{q}_{0}=0$ corresponds to a trivial single delta-barrier configuration with amplitude $2Z$.

The conditions for the existence of these solutions satisfy similar properties as those of complete-soliton solutions for the square well: for $m=0$, $X^{C}_{0}>X^{q}_{0}=0$ while for $m$ sufficiently large, $X^{C}_{m}<X^{q}_{m}$. Hence, reasoning in the same way, the $m=0$ $S-$ solution is limited to the range

[TABLE]

while for $m\geq 1$ we can only generally say that

[TABLE]

with $\tilde{X}_{m}\leq\min(X^{C}_{m},X^{H}_{m})$; check Appendix C.1.2 for the details.

The density profile of the $S-$ solutions for $m=0,1,2$ is represented in the central panels of Fig. 6.

V.2.3 $SD$ solutions

In this case:

[TABLE]

and this family of solutions is disconnected from the supersonic homogeneous solution with $E_{2}=E_{p}$. The wave function is also formally given by Eq. (99).

The limit solution $E^{m}_{2}=E_{q}$ corresponds to the upper limit of the $m+1$ $S-$ solutions, while $E^{m}_{2}=E_{1}$ is the soliton solution giving infinite cavity length so the condition for the existence of the $m$-th $SD$ solution is

[TABLE]

The density profile of the $SD$ solutions for $m=0,1,2$ is represented in right Fig. 6.

V.3 Asymmetric solutions

For asymmetric solutions:

[TABLE]

and the matching equations at the edges of the cavity are Eq. (91) and

[TABLE]

Following the discussion after Eq. (95), there are two possible values for $A_{\pm}$ for $E_{p}<E_{2}<E_{1}$. Since $A_{+}\neq A_{-}$, either $A_{-}<A_{+}$ or $A_{+}<A_{-}$; we fix them by imposing $A_{-}<A_{+}$ [the contrary case would just give the wave function resulting from the spatial inversion of $n(x)$]. This choice implies that $A_{p}<A_{+}<A_{\sup}$ and we distinguish two families of solutions according to the value of $A_{-}$: $AC$, for $A_{q}<A_{-}<A_{p}$; and $AD$, for $A_{\inf}<A_{-}<A_{q}$.

V.3.1 $AC$ solutions

The wave function takes the form:

[TABLE]

with $x_{\pm},\phi_{\pm}$ chosen such that the wave function and its derivative are continuous.

As the energy satisfies $E_{p}<E^{m}_{2}\leq E_{q}$, this family of solutions is also continuously connected to the BHL solution (86). In this limit, the critical lengths of the cavity are $X=X^{A,p}_{m}$, with $X^{A,p}_{m}$ given by Eq. (141), while the opposite limit, $E^{m}_{2}=E_{q}$, gives the critical lengths $X=X^{A,q}_{m}$, with $X^{A,q}_{m}$ given by Eq. (142).

Therefore, the $m$-th $AC$ solution only exists for

[TABLE]

The density profile of $AC$ solutions for $m=0,1,2$ is represented in Fig. 7.

V.3.2 $AD$ solutions

The wave function for this family of solutions is given by the same formal expression of Eq. (108) and the energy satisfies $E_{q}<E^{m}_{2}<E_{1}$. Reasoning as for the $SD$ solutions, the $m$-th solution only exists for:

[TABLE]

The density profile of $AD$ solutions for $m=0,1,2$ is represented in the rightmost panel Fig. 7.

VI Conclusions and outlook

In this work we have analyzed the use of more realistic models for black-hole lasers in Bose-Einstein condensates. First, we have proven a general result that associates a black-hole laser configuration to every compact black-hole solution. As an application, we have proposed two new black-hole laser configurations based on the waterfall and the delta-barrier configurations usually considered for studying analog black holes. In order to characterize them, we have provided a complete classification of the different families of non-linear stationary solutions as they are key to understand the stability of the system as well as its non-linear behavior.

Future works should explore in greater detail these configurations. For instance, a computation of the linear BdG spectrum would provide a further insight on the stability of the system and on the dynamics of the system at short times. Once obtained, a natural task would be to relate the appearance of dynamical instabilities with some of the families of stationary solutions presented in this work, following the ideas outlined in the main text.

On the other hand, in a similar way to Refs. Michel and Parentani (2015); de Nova et al. (2016), the non-linear black-hole lasing regime could be explored by a numerical simulation of the time-dependent Gross-Pitaevskii equation describing the evolution of the instability of the initial black-hole laser solutions $\psi^{BHL}(x)$. According to the results of Ref. de Nova et al. (2016), only two scenarios are expected for late times: either the system converges to the ground state of the system or it enters in a regime of continuous emission of solitons (CES). The characterization of the resulting phase diagram would provide more numerical data that could be useful for the elaboration of a more quantitative theory of the CES regime, which is currently lacking. We note that the production of such a soliton laser is of potential interest in quantum transport scenarios or the emergent field of atomtronics Seaman et al. (2007); Labouvie et al. (2015).

From an experimental point of view, the two black-hole laser configurations here presented describe more realistic scenarios than the typically used flat-profile configuration, as they provide simple models of external potentials easy to implement with standard experimental tools. Consistently, a more realistic numerical simulation should also take into account the complete time evolution of the configuration from the beginning, not just starting from the black-hole laser solution $\psi^{BHL}(x)$: as discussed in Refs. Tettamanti et al. (2016); Wang et al. (2016); Steinhauer and de Nova (2017), the time dependence of the problem is essential in the determination of the mechanism triggering the instability. However, regardless the specific transient of the system, the obtained stationary solutions should still be of great relevance in the non-linear dynamics occurring at long times after the onset of the instability. Among all of them, the family of $SH$ solutions for the attractive well and $S+$ solutions for the double delta-barrier are of special importance as they represent the true ground state of the system.

In addition, as a direct application of the results of the work, the black-hole laser model using an attractive square well is particularly interesting as it is expected to also provide a good description of the actual experimental configuration of Ref. Steinhauer (2014), much more accurate than the flat-profile configuration. In particular, following the reasoning of the above paragraph, the corresponding stationary states are expected play a key role in the description of future extensions of the experiment Steinhauer (2014) exploring the non-linear dynamics.

Acknowledgements.

I thank I. Carusotto, R. Parentani, F. Sols, J. Steinhauer and I. Zapata for useful discussions on the topic. This work has been supported by the Israel Science Foundation.

Appendix A Elliptic functions

We briefly present in this Appendix the definition of the different elliptic functions appearing throughout this work; see Refs. Abramowitz (1988); Byrd and Friedman (1971) for more details. The incomplete elliptic integral of the first kind is defined as

[TABLE]

and the Jacobi amplitude $\text{am}(u,\nu)$ is defined as its inverse function in the argument $\phi$ for fixed $\nu$ so $u=F\left[\text{am}(u,\nu),\nu\right]$.

The Jacobi elliptic functions $\text{sn}(u,\nu),\text{cn}(u,\nu)$ are defined as:

[TABLE]

As a consequence of the previous relations, $\text{sn}(u,\nu)$ and $\text{cn}(u,\nu)$ are periodic functions with period $4K(\nu)$, $K(\nu)\equiv F(\frac{\pi}{2},\nu)$ being the complete elliptic integral of the first kind, and satisfy $\text{sn}^{2}(u,\nu)+\text{cn}^{2}(u,\nu)=1$. In particular, in the limit of $\nu=0$, they reduce to the usual trigonometric functions, $\text{sn}(u,0)=\sin~{}u$, $\text{cn}(u,0)=\cos~{}u$ and $K(0)=\frac{\pi}{2}$ while in the upper limit, $\nu=1$, $\text{sn}(u,1)=\tanh~{}u$ and $K(1)=\infty$.

Another interesting function that appears when studying stationary solutions of the GP equation in an infinite well de Nova et al. (2014b) is the incomplete elliptic integral of the second kind

[TABLE]

with $E(\nu)\equiv E(\frac{\pi}{2},\nu)$ being the complete elliptic integral of the second kind.

Finally, the function $\Pi(\phi,m,\nu)$ is the incomplete elliptic integral of the third kind:

[TABLE]

and it appears when computing the phase of a cnoidal wave, see Eq. (53).

Appendix B Computation of the non-linear solutions of the attractive square well

We provide in this Appendix the technical details of the computation of the different families of stationary states for the attractive square well.

B.1 Homogeneous plane wave

As the wave function outside the well is the homogeneous subsonic plane wave, the amplitude at the edges of the well is fixed to $n_{W}=1$, so the amplitude energy inside the well is also fixed to the value:

[TABLE]

The roots of Eq. (71) are now $n^{H}_{2}=n_{W}=1$ and

[TABLE]

As $E_{2}$ is fixed, this solution only exists for certain critical lengths $X=X^{H}_{m}$, computed from Eq. (72):

[TABLE]

while $\alpha_{m}=(m+1)K(\nu_{m})$.

B.2 Subsonic shadow solitons

In the case where the wave function outside the well corresponds to shadow solitons, the matching condition (72) takes the form:

[TABLE]

from which we find that the amplitude energy $E^{m}_{2}$ is computed through the implicit equation

[TABLE]

and $\alpha_{m}=(m+1)K(\nu_{m})$.

B.3 Subsonic solitons

B.3.1 Symmetric solutions

Incomplete-soliton solutions.

In this case, the matching condition (72) gives:

[TABLE]

from which we obtain that

[TABLE]

and $\alpha_{m}=mK(\nu_{m})$.

The limit $n^{m}_{W}\rightarrow v^{2}$ gives $\nu_{m}\rightarrow 0$ and then all the elliptic functions are reduced to trigonometric functions. Specifically, for small values of the parameter

[TABLE]

we find from Eqs. (70), (71) that

[TABLE]

$M_{p}$ being the supersonic Mach number. Then, as $K(0)=\frac{\pi}{2}$, we find from Eq. (121) that the lengths in this limit are:

[TABLE]

Complete-soliton solutions.

The matching condition here is formally analog to that of Eq. (118) so $E^{m}_{2}$ is computed from Eq. (119) and $\alpha_{m}$ takes the same value; the difference is that now one has to take into account that the solutions outside the well are solitons and then $n^{m}_{W}$ is in the range $v^{2}\leq n^{m}_{W}\leq 1$. From Eqs. (102), (103), we can expect the behavior of $X^{SH}_{m}(E_{2})$ in Eq. (118) to be highly non-monotonic in the range $E^{C}_{2}\leq E_{2}\leq E^{H}_{2}$ for $m\geq 1$. This trend can be observed in Fig. 8, where $X^{SH}_{m}(E_{2})$ and $X^{SOL}_{m}(E_{2})$ are represented.

B.3.2 Asymmetric solutions

The matching condition (72) reads:

[TABLE]

with $m=1,2,3\ldots$ the number of complete periods inside the well. The value of $E^{m}_{2}$ is computed from

[TABLE]

and then

[TABLE]

Following a similar calculation as for the symmetric families, the small amplitude limit $n^{m}_{W}\rightarrow v^{2}$ gives the critical lengths

[TABLE]

Appendix C Computation of the non-linear solutions of the double delta-barrier

C.1 Symmetric solutions

We provide in this Appendix the technical details of the computation of the different families of stationary states for the double delta-barrier.

C.1.1 $S+$ solutions

The matching equations (95) give the same formal relations of Eq. (118) so $\alpha_{m}=(m+1)K(\nu_{m})$ and the possible energies inside the cavity are obtained from:

[TABLE]

We recall that $n_{W}(E_{2})$ is now a much more complicated function of $E_{2}$.

Following the discussion associated to Eq. (122), we compute the solution in the small amplitude limit near $\Psi^{BHL}(x)$, expanding all the magnitudes in terms of the parameter $\delta^{m}_{2}\equiv E^{m}_{2}-E_{p}$. The roots $n^{m}_{1,2}$ are computed from Eq. (123), where the supersonic Mach number is now

[TABLE]

The remaining root is $n^{m}_{3}\simeq v^{2}/n^{2}_{p}$. On the other hand, as

[TABLE]

from Eq. (87) we get:

[TABLE]

Using these results, we find that the critical lengths are

[TABLE]

where we have used the identity

[TABLE]

C.1.2 $S-$ solutions

The matching equations (95) give the same formal relations of Eq. (120) so $\alpha_{m}=mK(\nu_{m})$ and the possible energies inside the cavity are obtained from:

[TABLE]

The small amplitude limit near $E_{p}$ for this family of solutions gives the same critical lengths $X=X^{S-}_{m}(E^{m}_{2}\rightarrow E_{p})=X^{C}_{m}$ as the $S+$ solutions. The upper limit, $E^{m}_{2}=E_{q}$, only appears for discrete values of the length,

[TABLE]

$n^{q}_{i}$ being the corresponding roots for $E_{2}=E_{q}$, with $n^{q}_{1}=n_{q}$.

The behavior of $X^{S-}_{m}(E_{2})$ in Eq. (136) in the range $E_{p}\leq E_{2}\leq E_{q}$ is also highly non-monotonic for $m\geq 1$, similar to that of $X^{SH}_{m}(E_{2})$ in Fig. 8.

C.1.3 $SD$ solutions

The matching equations (95) are the same as for $S+$ solutions so $\alpha_{m}=(m+1)K(\nu_{m})$ and the possible energies $E^{m}_{2}$ inside the cavity are computed from Eq. (130).

C.2 Asymmetric solutions

C.2.1 $AC$ solutions

In this case, the asymmetric matching condition of Eq. (107) gives:

[TABLE]

from which we obtain that:

[TABLE]

and

[TABLE]

The limit $E^{m}_{2}\rightarrow E_{p}$ gives the critical lengths

[TABLE]

while, in the opposite limit, $E^{m}_{2}=E_{q}$,

[TABLE]

with $n^{q}_{+}$ the density of the $A^{q}_{+}>A_{p}$ solution to Eq. (91) for $E_{2}=E_{q}$.

C.2.2 $AD$ solutions

The asymmetric matching condition of Eq. (107) gives now:

[TABLE]

which yields

[TABLE]

and

[TABLE]

where $E_{q}<E^{m}_{2}<E_{1}$.