Non-Linear Stationary Solutions in Realistic Models for Analog Black-Hole Lasers

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.


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) [1,2]. Within a similar scheme, Corley and Jacobson [3] showed that a bosonic field with a 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 62M /M nK, with M 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 [4], 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 as different as Fermi gases [5], ion rings [6], polaritons [7] or, in a classical context, surface waves in a water tank [8].
Of particular interest are the analogues implemented in Bose-Einstein condensates (BEC), first suggested by Garay et al. [9]. 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 [10][11][12][13][14][15][16][17][18][19][20][21][22][23], 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 [16,[24][25][26][27][28][29].
Regarding the experimental side, the first acoustic horizon in a BEC was produced by the Technion group [30] 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 [31], although there is still some discussion in the community about the interpretation of the experimental results [32][33][34]. Recently, the same group provided the first experimental evidence of the emission of Hawking radiation by measuring the entanglement of the emitted phonons [35].
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 [36] or an optical lattice [20]; the waterfall configuration described in [37] is a theoretical model of the actual experimental setups of [30,35].
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 [38,39], 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 a square attractive well can be regarded as a model for studying the experimental BHL of [31].
The scheme of the paper is the following. In Section 2, we revisit the basic theory of gravitational analogues in BEC. The general relation between BH and BHL solutions is proven in Section 3. In Sections 4 and 5, we study the different stationary states for BHL configurations with a square-well and a double delta-barrier, respectively. Conclusions are presented in Section 6. Appendix A is devoted to introducing the different elliptic functions used in this work, while Appendices B and C are devoted to the technical details of the calculations presented in the main text.

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, [14,[39][40][41].

Effective One-Dimensional Configurations
We begin by reviewing how to reach an effective one-dimensional (1D) configuration, the so-called 1D mean-field regime [38,42]. For that purpose, we consider a 3D gas of N bosons of mass m near T = 0 (more precisely, T T c , with T c the critical temperature of the condensate), described by the second-quantization Hamiltonian [43,44]: where V ext (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 3D = 4πh 2 a s /m the corresponding coupling constant and a s the s-wave scattering length [45,46]. In order to obtain an effective 1D configuration along the x-axis, we consider a total external potential of the form V ext (x, t) = V(x, t) + V tr (y, z), where V(x, t) only depends on the x coordinate, while V tr (y, z) = 1 2 mω 2 tr ρ 2 (ρ = 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: with a tr being the transverse harmonic oscillator length. After integrating over the transverse degrees of freedom, we arrive at the following 1D effective Hamiltonian: where we have absorbed the resulting zero-point energyhω tr of the harmonic oscillator, and the effective 1D constant coupling is g 1D = 2hω tr a s . More specifically, the condition for the approximation of Equation (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) hω tr , which can be simply put as n 1D (x)a s 1, with n 1D (x) the 1D-density.
In the same fashion, the 3D canonical commutation rules for the field operator: are reduced to the 1D version: As we will only deal with 1D configurations, in the following, we omit everywhere the 1D index.

Gross-Pitaevskii and Bogoliubov-De Gennes Equations
Using the Hamiltonian of Equation (3) and the corresponding canonical commutation rules of Equation (5), we write the equation of motion for the field operatorψ(x) in the Heisenberg picture: Since there is a condensate, we can perform a mean-field approximation: with ψ(x, t) the Gross-Pitaevskii (GP) wave function [45] describing the condensate andφ(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: Assuming that the depletion cloud (i.e., the cloud formed by the atoms outside the condensate) is negligible, ψ(x, t) is normalized to the total number of particles: The conservation of the norm of the GP wave function is guaranteed by the same relation as in the usual linear Schrodinger equation: 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, ψ(x, t) = A(x, t)e iφ(x,t) , is the current and: are the mean-field density and flow velocity, respectively. Interestingly, the first line of Equation (11) is the equivalent of the continuity equation for a hydrodynamical fluid. On the other hand, taking the spatial derivative in the second line gives: 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 = gn 2 2 . The only difference is the rightmost term, which is a genuine quantum feature, as it containsh, 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 [4] (see also the discussion in the next subsection).
On the other hand, to lowest order in the quantum fluctuations of the field operator,φ(x, t), one finds from Equation (6) that: 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: that are of special interest as they describe stationary configurations. In particular, the stationary wave function ψ 0 (x) obeys the time-independent GP equation: 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: 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: with the spatial derivative. The phase is simply obtained from the relation: Note that, if J = 0, neither the amplitude nor the flow velocity vanish. For a fixed value of the number of particles, there can be several different solutions for the GP Equation (16). The true ground state of the system is that minimizing the grand-canonical energy, K = E − µ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: it can be seen that Equation (16) is precisely the condition for ψ 0 (x) to be an extreme of K. In that case, K takes the simple form: Solutions of Equation (16) that 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 µ playing the role of the Lagrange multiplier.
With respect to the quantum fluctuations,φ(x, t), Equation (14) is now a stationary problem of the form: As it is a linear equation, we can expand the field operator in terms of eigenmodes: where the spinors z n satisfy the time-independent BdG equations: Due to the structure of the equations, the conjugatez n is also a mode with energy − * n . An interesting property of the eigenvalue problem of Equation (24) is that it is non-Hermitian, and thus, it can yield complex eigenvalues. In particular, eigenvalues with a 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: with σ z = 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 conjugatez n has the opposite sign to that of z n , (z n |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 n , m are orthogonal, as seen from the relation: from which it also follows that modes with complex frequency have zero norm.

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 ψ 0 (x) = √ ne iqx+φ 0 , with n the density of the condensate, q its momentum and φ 0 some phase. After removing the phase of the condensate from the field operator,φ(x, t) → e iqx+φ 0φ (x, t), it is straightforward to show that the eigenmodes of the BdG Equation (24) are plane waves with wave vector k and frequency ω, giving rise to the following dispersion relation: with c = gn/m the sound velocity, v =hq/m the constant flow velocity, ξ ≡h/mc the so-called healing length and Ω the comoving frequency. The above dispersion relation gives four different wavevectors for a given value of the frequency. In fact, Equation (27) is just the usual Bogoliubov dispersion relation for phonons in a condensate at rest, Ω(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, Ω > 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 Figure 1, where the blue (red) curves represents the sign +(−) branches of the dispersion relation, ω ± (k) = vk ± Ω(k), and also positive (negative) normalization according to the scalar product of Equation (25). Indeed, the − branch is just the dispersion relation of the conjugate modes of the + branch, ω − (k) = −ω + (−k) (see Equation (24) and the 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 a complex wave vector. On the other hand, for supersonic flows, in the window −ω max < ω < ω max , all four modes are propagating, where the threshold frequency ω max is ω max = max k ω − (k) and is marked by a horizontal dashed line in the right of Figure 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, withhk 0 = 2m √ 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 [47], characterized by the wave vector k 0 . The previous magnitudes can be extended to non-homogeneous configurations by taking c(x) ≡ gn(x)/m and v(x) as defined in Equation (12). In a 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 Equation (22) in terms of the relative quantum fluctuations,φ( with D t = ∂ t + v(x)∂ x the comoving derivative. Gathering this equation with its complex conjugate and defining the hermitic fields:χ gives a pair of equations of the form: The above fields are related to physical magnitudes asρ(x, t) = 2n(x)χ + (x, t) and φ(x, t) =χ − (x, t), withρ,φ 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 Equation (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 Equation (30), which precisely amounts to work in the hydrodynamic regime where all of the contributions arising from the quantum pressure are neglected. In this approximation, we can write the equation for the phase fluctuations as: which can be rewritten as the relativistic Klein-Gordon equation for a massless scalar fieldφ on a metric g µν , with the effective stationary metric g µν given by: Thus, the points where c(x) = v(x) are the horizons of the acoustic metric g µν (x), analogous to astrophysical event horizons. Using a simple physical picture, for acoustic phonons (long-wavelength excitations), the dispersion relation of Equation (27) has the form ω ± (k) (v ± 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 [14,40], 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 taking 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 [26]: 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.

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, Equation (18), is analogous 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: 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 Equation (18): This polynomial equation for the density only has (two) real positive roots whenever: In the rest of the work, we will assume that the condition of Equation (36) 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 ψ 0 (x) → √ n 0 ψ 0 (x), so it becomes dimensionless, and take units such thath = m = c 0 = 1, with c 0 = gn 0 /m the sound velocity associated to the density n 0 . Length, time and energy are measured in units of ξ 0 =h/mc 0 , t 0 = ξ 0 /c 0 and E 0 = mc 2 0 , respectively. Furthermore, we will refer to v 0 as 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 µ = 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 Equation (35) as: from which we immediately obtain the density of the other homogeneous solution: 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, p . By defining z ≡ 8/v 2 and observing that the function f (z) = z 2/3 − 1 − √ 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 We also study the non-homogeneous solutions of Equation (34). For that purpose, we represent the amplitude potential W(A) in Figure 2. The local minimum corresponds to the homogeneous supersonic solution A p = √ n p (which means that it is a stable fixed point) and the local maximum to the homogeneous subsonic solution A 0 = √ n 0 = 1 (which means that it is an unstable fixed point). Rewriting Equation (34) in terms of the density gives the simplified equation: 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: 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 Equation (40) 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 Figure 2. By integrating Equation (39), we find that: with x 0 some integration constant arising from the translational invariance of the problem. The solution of the previous indefinite integral is given in terms of elliptic functions. The resulting phase of the wave function is computed from Equation (19), obtaining: where: Λ(x, n 1 , n 2 , n 3 , α) ≡ n(x, n 1 , n 2 , n 3 , α)e iφ(x,n 1 ,n 2 ,n 3 ,α) (43) n(x, n 1 , with φ(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 the blue curve for A > 1 in the right plot of Figure 2. Moreover, for high values of n, n ∝ n 3/2 , and then, the solution blows up at some finite value x bu as n(x) ∼ (x − x bu ) −2 . The same reasoning holds for E A < W(A p ) (the green dashed-dotted line of Figure 2) or E A > W(1) (the red dashed line of Figure 2). These exploding solutions are not relevant for the present work, so we ignore them in the following.
Finally, we consider the degenerate cases One possible solution corresponds to the stable fixed point of the homogeneous supersonic solution n(x) = n p , described by the plane wave: 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: For n(x) = 1, we get from Equation (41): For n(x) < 1, we obtain: which is of the same form of Equation (43) after taking into account that 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: being φ 0 some constant phase. This solution represents a soliton with zero velocity [46].
On the other hand, taking n > 1 in Equation (46) gives the so-called shadow soliton solution [27]: 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 Section 4 and 5.

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 of the BH configurations usually appearing in the literature [37].
The proof is straightforward. Consider a compact BH configuration, which satisfies a time-independent GP equation of the form: 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: with A sp , q sp the supersonic amplitude and momentum and Ψ 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: which satisfies the following GP equation: where the potential V BHL (x) is given by: The wave function of Equation (52) 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 Figure 3, respectively. One of the most considered BH configurations is the flat-profile configuration (upper left of Figure 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 (x) and the external potential V (x) 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 [12][13][14] because of its analytical simplicity.  More realistic configurations are displayed in the central and right panels of the upper row of Figure 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 [36], 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 [30,35], 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 [35].
Regarding the BHL side, the associated flat-profile BHL of lower left Figure 3 has been already considered in the literature in both analytical and numerical studies [27][28][29]. Specifically, [27] 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 that the appearance of each dynamical instability is associated with the appearance of a stationary solution with lower grand-canonical energy K than the initial homogeneous plane-wave solution.
On the other hand, [28,29] 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, [29] found that the system only presents two kinds 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 with the delta-barrier and the waterfall configurations, depicted in the lower center and right panel of Figure 3, respectively. Specifically, due to energy and particle number conservation [27], we only aim at the GP stationary solutions asymptotically matching at ±∞ 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 [27]. 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 [27,29] regardless of the specific mechanism giving rise to the initial growth of the instability at short times [32][33][34].
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 [31], 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 Section 2.4, we set units in the rest of this work such thath = 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 [29].
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.

General Structure
As a first step, we present the BH solution corresponding to the waterfall configuration, where the external potential is given by: where θ(x) is the step function and x H the point where the step is placed. The corresponding GP wave function is given by: with φ 0 some phase to make the wave function continuous at x = x H . Following the procedure described in Section 3, the associated BHL configuration is created by using an attractive square well potential of size X, We note that the different stationary solutions of a condensate flowing through an attractive square well were first addressed in [38]. Here, we restrict to the specific case where V 0 is given by Equation (55), so a BHL solution as that of the lower right of Figure 3 exists, described by the GP wave function: 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 Equation (34). The asymptotic boundary conditions fix the current to J = v and the amplitude energy outside the well to E A = E 1 = 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: where W i (A), E i are the amplitude potential and the conserved amplitude energy for the i = 1, 2 regions, with µ 1 = µ = 1 + v 2 /2 and µ 2 = µ 1 + V 0 = v 2 + 1 2v 2 . Invoking the continuity of the wave function and its derivative, we find the matching condition at both edges: As E 1 is fixed, we only need to find the possible values of E 2 , which make the GP wave function satisfy both Equations (59) and (60).
The situation is schematically depicted in the left of Figure 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 = ± X 2 must be the same, the only possible solution inside the well is a cnoidal wave, as that of Equation (43), where n i , i = 1, 2, 3 are now the roots of the equation: For convention, we choose the wave function such that it is real at x = 0, φ(0) = 0. The matching of the cnoidal wave (see Equation (43)) at x = ± X 2 gives two equations: n ± X 2 , n 1 , n 2 , n 3 , α = n W (62) Since n 1 , n 2 , n 3 , n W are functions of E 2 , the above system gives two conditions for two variables, α 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 . . . representing the number of complete periods inside the well. As a result, we only need to compute the corresponding values of E 2 , α, labeled as E m 2 , α m , and the associated parameters n m W , n m i , ν m (see Equations (43) and (60)) 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.

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 type of solution only appears for some critical values of the length X = X H m , m = 0, 1, 2 . . ., and the corresponding GP wave function is: 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 with this family of solutions is depicted as a red solid line in left Figure 4, while the density profile for ψ H m (x), m = 1, 2, 3, is shown in the central panels of the same figure.

Asymptotic Shadow Solitons
The resulting GP wave function for this family of solutions reads: with φ 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 ≤ 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, The limit case n m W = 1 corresponds to E m 2 = E H 2 and Ψ SH m (x) = Ψ H m (x). In particular, the function Ψ 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 (21). Hence, the BHL solution Ψ BHL (x) of Equation (58) should be dynamically unstable for any length X > 0 of the well. However, since Ψ 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 [48], 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 ≤ E m 2 < E SH 2 , and the m-th solution only exists for lengths: The density profile of ψ SH m (x) for m = 0, 1, 2 is represented in the rightmost panel of Figure 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 ψ 0 (x) = ψ * 0 (−x), and asymmetric solutions, with no spatial symmetry.

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: For incomplete-soliton solutions, x 0 satisfies X/2 − x 0 > 0. Since the solutions outside the well are solitons, v 2 ≤ n m W ≤ 1. The lower boundary n W = v 2 gives the same solution of Equation (58), in which n m , existing for arbitrary length X. Hence, this family of solutions is continuously connected to the BHL configuration described by ψ BHL (x). A more detailed analysis of the limit n m W → v 2 (see Appendix B.3) shows that these solutions appear at the critical lengths X C m , given by Equation (A14), and can be understood as small amplitude perturbations on top of the GP wave function of Equation (58), described by the BdG equations. Precisely, in this limit, all of 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 the right of Figure 1 and the related discussion). Indeed, as the dynamically unstable modes arising from the BHL effect are expected to first show as zero-frequency BdG modes [27], the critical lengths X C m should also signal the appearance of a new dynamical instability. The upper limit, n m W = 1, gives ψ SOL m (x) = ψ H m+1 (x), merging with the ψ SH m+1 (x) solutions. Then, the m-th solution of Equation (67) only exists for: The density profile for the incomplete-soliton solutions m = 0, 1, 2 is represented in the upper row of Figure 5.

Complete-Soliton Solutions
The wave function is of the same form as that of Equation (67), but with x 0 satisfying X/2 − x 0 < 0. This family of solutions is also continuously connected to ψ 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 ψ SOL m (x) = ψ H m (x), merging with the ψ 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: while for m ≥ 1, we can only generally say that: withX m ≤ 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 Figure 5.

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. In particular, they contain an exact integer number of periods m = 1, 2, 3 inside the well (see Eq. (A15)), with m = 0 giving the trivial solution of no well, X = 0. The corresponding GP wave function reads: where δx > 0, φ L , φ R are chosen such that the wave function and its derivative are continuous and ± corresponds to the case of complete-incomplete (incomplete-complete) solitons.
In the same fashion as the symmetric families, in the limit n m W → v 2 , the asymmetric solutions are continuously connected to ψ BHL (x), appearing at the critical lengths X A m given by Equation (A18), while in the upper limit n m W → 1, they converge to the homogeneous ψ H m (x) solutions.
Hence, the m-th asymmetric solution is restricted for lengths: The density profile for the asymmetric solutions m = 1, 2, 3 is represented in the lower row of Figure 5.

General Structure
The potential corresponding to a single delta-barrier configuration is: with W(A) given by Equation (34), A p the supersonic amplitude of Equation (38) and x H the point where the barrier is placed. Such a delta potential originates a discontinuity in the derivative of ψ 0 (x) of the form The resulting GP wave function is: The constants x 0 , φ 0 are fixed by imposing the continuity of ψ 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, This configuration was studied in [39] in order to look for resonant BH configurations, which enhance the spontaneous Hawking signal [17]. 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 the lower central part of Figure 3 exists, described by the GP wave function: with x 0 , φ 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: where E i is the conserved amplitude energy for the i = 1, 2 regions, E 1 = 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: which, in terms of the amplitude, reads: Thus, we can understand the effect of the delta barriers as "instantaneously accelerating" the classical particle described by Equation (77).
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 restricts even more the possibilities, and the solitons must satisfy: Inside the cavity, the solution corresponds to a cnoidal wave, with the phase chosen such that φ(0) = 0 for convention. By joining Equations (77), (79) and (80), we find that the energy inside the cavity is related to the amplitude at the edges through: Hence, the value of E 2 is determined by the cnoidal waves arising from Equation (77) that are compatible with Equation (81). In particular, as E p ≤ E 2 < E 1 , with E p ≡ W(A p ), the amplitudes at the edges of the cavity must be in the range: where A inf < A sup < 1 are obtained from the roots of the equation: As for the square well, the possible solutions are labeled by a discrete index m = 0, 1, 2 . . . representing the number of complete periods of the cnoidal wave inside the cavity, and the wave function is determined once E m 2 , α 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 Section 4.4.

Symmetric Solutions
Symmetric solutions satisfy: and then, the matching equations at the edges of the cavity read: n W = n ± X 2 , n 1 , n 2 , n 3 , α 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 Equation (62); 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 Equation (80), we classify the solutions according to the sign of the derivatives at the internal side of the edges of the cavity, x = ±X ∓ /2. The limit values between the different solutions are obtained from: which yields a similar equation to Equation (83): 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 with these solutions are E 2 = E q ≡ 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 .

S+ Solutions
In this case: The corresponding wave function for this family of solutions is: (89) 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 ψ S+ m (x) = ψ BHL (x), which exists for an arbitrary length. Following the results of Section 4.4, we analyze the limit E m 2 → 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 Equation (A23). As for the square well, they can be understood as small perturbations on top of the GP wave function ψ 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 ψ S+ m (x) has lower grand-canonical energy than ψ BHL (x); specifically, ψ 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 ψ BHL (x), and ψ 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 of Equation (76), in the same lines of [27].
The upper limit, E m 2 = E 1 , corresponds to the soliton solution, which gives an infinite value for the cavity length X. Hence, the m-th S+ solution only exists for: The density profile of ψ S+ m (x) for m = 0, 1, 2 is represented in left Figure 6.

S− Solutions
In this case: The wave function is given by the same formal expression of Equation (89). This family of solutions is also continuously connected to ψ 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 Equation (A26); 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: while for m ≥ 1, we can only generally say that: withX m ≤ 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 Figure 6.

SD Solutions
In this case: 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 Equation (89). 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: The density profile of the SD solutions for m = 0, 1, 2 is represented in right Figure 6.

Asymmetric Solutions
For asymmetric solutions: and the matching equations at the edges of the cavity are Equation (81) and: n ± = n ± X 2 , n 1 , n 2 , n 3 , α , n ± = A 2 ± (97) Following the discussion after Equation (85), there are two possible values for A ± for E p < E 2 < E 1 . Since A + = 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 .

AC Solutions
The wave function takes the form: with x ± , φ ± chosen such that the wave function and its derivative are continuous. As the energy satisfies E p < E m 2 ≤ E q , this family of solutions is also continuously connected to the BHL solution of Equation (76). In this limit, the critical lengths of the cavity are X = X A,p m , with X A,p m given by Equation (A30), 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 Equation (A31). Therefore, the m-th AC solution only exists for: The density profile of AC solutions for m = 0, 1, 2 is represented in Figure 7.

AD Solutions
The wave function for this family of solutions is given by the same formal expression of Equation (98) and the energy satisfies E q < E m 2 < E 1 . Reasoning as for the SD solutions, the m-th solution only exists for: The density profile of AD solutions for m = 0, 1, 2 is represented in the rightmost panel of Figure 7.

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 into the stability of the system and 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 [28,29], 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 ψ BHL (x). According to the results of [29], 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 [49,50].
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 that are 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 ψ BHL (x): as discussed in [32][33][34], the time dependence of the problem is essential in the determination of the mechanism triggering the instability. However, regardless of 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 [31], much more accurate than the flat-profile configuration. In particular, following the reasoning of the above paragraph, the corresponding stationary states are expected to play a key role in the description of future extensions of the experiment [31] exploring the non-linear dynamics.
Another interesting function that appears when studying stationary solutions of the GP equation in an infinite well [20] is the incomplete elliptic integral of the second kind: with E(ν) ≡ E( π 2 , ν) being the complete elliptic integral of the second kind. Finally, the function Π(φ, m, ν) is the incomplete elliptic integral of the third kind: and it appears when computing the phase of a cnoidal wave; see Equation (43).

Complete-Soliton Solutions
The matching condition here is formally analogous to that of Equation (A8), so E m 2 is computed from Equation (A9), and α 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 ≤ n m W ≤ 1. From Equations (92) and (93), we can expect the behavior of X SH m (E 2 ) in Equation (A8) to be highly non-monotonic in the range E C 2 ≤ E 2 ≤ E H 2 for m ≥ 1. This trend can be observed in Figure A1, where X SH m (E 2 ) and X SOL m (E 2 ) are represented.  Figure A1. Plot of X SH m (E 2 ) (solid lines) and X SOL m (E 2 ) (dashed lines) as a function of ∆ ≡ (E 2 − E C 2 )/(E H 2 − E C 2 ) in the range ∆ ∈ [0, 1]. The horizontal solid and dashed-dotted lines correspond to the limit values X = X C m and X = X H m , respectively. 1 , ν The limit E m 2 → E p gives the critical lengths: with n q + the density of the A q + > A p solution to Equation (81) for E 2 = E q .