Long time universality of black-hole lasers

For flowing quantum gases, it has been found that at long times an initial black-hole laser (BHL) configuration exhibits only two possible states: the ground state or a periodic self-oscillating state of continuous emission of solitons. So far, all the works on this subject are based on a highly idealized model, quite difficult to implement experimentally. Here we study the instability spectrum and the time evolution of a recently proposed realistic model of a BHL, thus providing a useful theoretical tool for the clear identification of black-hole lasing in future experiments. We further confirm the existence of a well-defined phase diagram at long times, which bespeaks universality in the long-time behavior of a BHL. Additionally, we develop a complementary model in which the same potential profile is applied to a subsonic homogeneous flowing condensate which, despite not forming a BHL, evolves towards the same phase diagram as the associated BHL model. This result reveals an even stronger form of robustness in the long-time behavior with respect to the initial condition which goes beyond what has been described in the previous literature.

For flowing quantum gases, it has been found that at long times an initial black-hole laser (BHL) configuration exhibits only two possible states: the ground state or a periodic self-oscillating state of continuous emission of solitons. So far, all the works on this subject are based on a highly idealized model, quite difficult to implement experimentally. Here we study the instability spectrum and the time evolution of a recently proposed realistic model of a BHL, thus providing a useful theoretical tool for the clear identification of black-hole lasing in future experiments. We further confirm the existence of a well-defined phase diagram at long times, which bespeaks universality in the long-time behavior of a BHL. Additionally, we develop a complementary model in which the same potential profile is applied to a subsonic homogeneous flowing condensate which, despite not forming a BHL, evolves towards the same phase diagram as the associated BHL model. This result reveals an even stronger form of robustness in the long-time behavior with respect to the initial condition which goes beyond what has been described in the previous literature.

I. INTRODUCTION
The spontaneous emission of radiation by the event horizon of a black hole (BH) [1], known as Hawking radiation, is one of the most celebrated predictions of modern theoretical physics, combining thermodynamics, quantum mechanics, and general relativity. However, its detection in an astrophysical context is unlikely in the foreseeable future due to its extremely low characteristic temperature, much lower than that of the cosmic microwave background. An alternative approach to study such gravitational phenomena is provided by a flowing fluid [2], in which a subsonic-supersonic interface becomes the acoustic analog of an event horizon. Based on this observation, a number of analog setups emerged in systems as different as water waves [3,4], nonlinear optical fibers [5,6] or quantum fluids of light [7,8].
Another interesting analog feature that can be observed in a Bose-Einstein condensate due to its superluminal dispersion relation is the so-called black-hole laser (BHL) effect [26], in which a configuration displaying a pair of horizons can give rise to the self-amplification of Hawking radiation due to successive reflections between them, like in a laser cavity, translated into the appearance of dynamical instabilities in the spectrum of excitations [27][28][29][30][31]. A complete experimental characterization of the black-hole laser effect is still one of the major present challenges in the gravitational analog field, since its first reported observation [32] has sparked intense discussions within the community [33][34][35][36]. There are also attempts to observe the black-hole laser effect in other analog systems with superluminal dispersion relation [37,38].
One of the most interesting features of a black-hole laser is its long-time behavior [39,40], once the initial instability has grown up to saturation and the system enters the full nonlinear regime. This regime provides a unique setup to test backreaction [41], in contrast to the case of the usual spontaneous Hawking effect, whose backreation is too small to be observed in an actual experiment. In particular, it has been seen that the longtime behavior of the system displays a well-defined phase diagram, characterized by solely the parameters of the setup and independent of the particular initial condition [40]. Specifically, for sufficiently long times, the system only exhibits two possible states: either it reaches the true ground state (GS) by evaporating away the horizons or it enters a regime of continuous emission of solitons (CES). Due to its periodic and radiating nature, the CES regime has been argued to be the closest counterpart of an actual optical laser, and the true black-hole laser [40].
Nevertheless, all the works studying the nonlinear behavior of black-hole lasers [39,40,42] have so far considered the so-called flat-profile model [11,12] which, albeit very useful for analytical calculations due to its simplicity, is quite difficult to implement experimentally. A more realistic approach was provided in Ref. [43], where it was shown that from every single black-hole solution a black-hole laser solution can be constructed. That result was used to propose an adequate theoretical description of the black-hole laser configuration in the actual experiment [32]. The advantage of this theoretical model is that it allows to isolate the genuine contribution of black-hole lasing in a setup similar to that used for the real experiment, where this contribution has so far been overshadowed by the Bogoliubov-Cherenkov radiation background [25].
In this work we present a comprehensive study of the long-time behavior of such realistic model of a black-hole laser, computing also the linear spectrum that describes its unstable behavior at short times. We confirm and extend previous results in the literature for both the linear spectrum and the long-time regime [39,40,42]. In particular, at long times we find a well-defined phase diagram between GS and CES, which hints at a possible form of long-time universality in the behavior of blackhole lasers. Furthermore, we develop a novel approach to the study of the long-time behavior of a BHL configuration that yields a rather complete understanding of the problem. Specifically, we analyze the time evolution in a complementary scenario where the same potential profile as before is suddenly applied now to a subsonic homogeneous flowing condensate which, not being in a stationary state of the new potential, evolves deterministically. We find that its long-time phase diagram is identical to that for the BHL solutions in the same setup. This reveals a stronger form of robustness with respect to the transient than previously reported. Namely, it shows that GS and CES are intrinsic states of a flowing condensate that are eventually reached independently of the early-time dynamics.
This paper is arranged as follows. In Sec. II we present the configurations analyzed throughout this work. Section III is devoted to a preliminary theoretical analysis of the BHL solutions considered in this work, analyzing both their linear and nonlinear spectra. In Sec. IV we study numerically the time evolution of the initial unstable BHL solutions, verifying what has been anticipated in the previous section. Section V deals with the time evolution of an alternative configuration in which the initial condition is a subsonic homogeneous flowing condensate, and includes the computation of its long-time phase diagram. In Sec. VI the results of the previous sections are connected and discussed. Conclusions and outlook are presented in Sec. VII.

II. ANALOG CONFIGURATIONS IN FLOWING CONDENSATES
We devote this section to present the main theoretical models used along the work; the interested reader is referred to the literature for the complete technical details [39,40,43]. We consider a one-dimensional Bose-Einstein condensate close to T = 0, described by the Gross-Pitaevskii (GP) wave function Ψ(x, t), whose associated sound and flow velocities are where g is the coupling constant, m the mass of the atoms, and we are using a density-phase decomposition of the wave function Ψ(x, t) = n(x, t)e iα(x,t) . Following the usual convention, hereafter we set units such that = m = c u = 1 and rescale the GP wave function as Ψ(x, t) → √ n u Ψ(x, t) so it becomes dimensionless, with n u , c u the asymptotic subsonic density and speed of sound, respectively.
We first revisit the BH waterfall configuration [15], depicted in Fig. 1a, created by introducing a semi-infinite attractive potential where v is the asymptotic subsonic Mach number and x H the point where the step is placed.
The corresponding stationary GP wave function describing the BH configuration is with φ 0 some phase to make the wave function continuous at x = x H . This configuration provides a realistic theoretical model of the experimental setup of Refs. [22][23][24], in which a step potential is swept with velocity v along a condensate at rest, an equivalent situation to that described here via Galilean invariance. We note that this family of solutions is uniparametric since it is determined by solely the value of v.
As explained in Ref. [43], a stationary BHL solution of arbitrary size can be associated to every stationary BH solution. The argument goes as follows: if a BH configuration is formed by a potential V BH (x) and described by a stationary GP wave function of the form where Ψ d (x) = √ n d e iv d x is here the homogeneous plane wave characterizing the supersonic region and we set x H = 0 without loss of generality, then the configuration formed by the symmetric potential admits as a stationary solution the GP wave function We note that the relations V BHL (x) = V BHL (−x) and Ψ 0 (x) = Ψ * 0 (−x) are satisfied, with Ψ 0 (x) describing a stationary BHL solution as it presents a pair of horizons, namely, a black hole for x < 0 and a white hole (WH) for x > 0. We note that the length X can be chosen arbitrarily, becoming a free parameter.
For the waterfall configuration, such a BHL solution is created from the BH solution of Eq. (3) by means of an attractive square well potential of size X, with the potential depth V 0 fine-tuned to the value V 0 = V 0 (v) given by Eq. (2). The resulting BHL configuration is described by the following stationary GP wave function which can be seen in Fig. 1b. This solution is determined by two parameters, v and X. Indeed, we can regard this BHL solution as an extension for X > 0 of the stationary soliton solution when no potential well is present [44,45]. The particular relevance of this BHL model is that, while being simple and allowing for an analytical study of a number of its properties due to its homogeneity in the supersonic region, it is similar to the actual configuration observed in the laboratory [32], created by sweeping a finite attractive well along a trapped condensate. Moreover, since it is a stationary solution, the pure black-hole lasing contribution can be isolated from other possible background phenomena such as Bogoliubov-Cherenkov radiation [46], so a study of the genuine black-hole laser properties can be carried out with the perspective of a clear identification in future experiments.
We will refer to Eq. (8) as the fine-tuned BHL solution since the value of V 0 is fixed by the value of v, and the solution inside the well is exactly a plane wave. Other types of BHL solutions, displaying a pair of BH-WH horizons, can be found for the same potential and asymptotic boundary conditions [43]. Two main families of BHL solutions can be distinguished: symmetric, Ψ S n (x) = Ψ S * n (−x), and asymmetric, Ψ A n (x) = Ψ A * n (−x), with n = 0, 1, 2 . . . The drawback is that they are not homogeneous within the well but instead are described in terms of complicated elliptic functions. Another disadvantage is that each type of solution only exists within a restricted range of values of X for a given value of v. In a general experimental situation, however, V 0 is not locked to the value V 0 (v). When exploring the complete parameter space (v, X, V 0 ), allowing now the well amplitude V 0 in Eq. (7) to be an extra degree of freedom, the fine-tuned BHL solution cease to exist. Still, one can find stationary BHL solutions, given by extensions for arbitrary V 0 of the symmetric and asymmetric families of solutions described in the previous paragraph. We will denote these solutions as generic BHL solutions. An example of a generic symmetric BHL solution is depicted in Fig. 1c. A more detailed discussion about these families of solutions can be found in Sec. III B.
Finally, the last configuration considered in this work (see Sec. V) is that in which an attractive square well is suddenly introduced in an initially homogeneous flowing condensate (IHFC). In contrast to the previous cases, this setup does not contain horizons and thus is not a BHL setup. An important feature of this model is that the described configuration is not stationary and so it evolves in time deterministically with the homogeneous plane wave Ψ 0 (x) = e ivx acting as an initial condition. This configuration is then characterized by the same set of parameters (v, X, V 0 ) as the generic BHL solutions, providing a useful tool to complement the study of the time evolution. A scheme of this configuration is depicted in Fig. 1d.

III. LINEAR AND NONLINEAR SPECTRA
Before focusing on the time evolution of the BHL solutions, we first study their linear spectrum of excitations, described by the usual Bogoliubov-de Gennes (BdG) equations, and the full nonlinear spectrum of stationary GP solutions coexisting for the same background parameters and asymptotic boundary conditions. This analysis is expected to help identifying the main features of the time evolution and elaborate some predictions a priori [39,40,42].

A. Linear spectrum
We compute in this subsection the linear spectrum of dynamical instabilities as given by the BdG equations. We restrict the analysis to the case of the the fine-tuned BHL solution of Eq. (8) since analytical results can be obtained in this background, and it exists for arbitrary cavity length X, thus allowing for the use of X as a control parameter to study the behavior of the instabilities. The results for the linear spectrum of this subsection then complement the work of Ref. [42], which is restricted to the idealized flat-profile configuration, and the work of Ref. [43], which only addressed the nonlinear spectrum of stationary GP solutions (see next section). The procedure for computing complex-frequency modes of the BdG equations in a BHL is explained in detail in Refs. [39,42].
The spectrum of dynamical instabilities is represented in the left panel of Fig. 2 as a function of the cavity length X. We observe that, for the critical lengths k 0 being the wave vector of the Bogoliubov-Cherenkov mode [46], a new dynamically unstable mode appears with a purely imaginary frequency (representing a degenerate mode since it is its own BdG conjugate, ω = −ω * ), but acquiring a non-vanishing real part of the frequency (hence becoming non-degenerate as ω = −ω * ) in the middle of the interval 2π/k 0 between the onsets of instabilities, i.e., at lengths The dominant mode (that with the largest growth rate) is typically that with the largest n, except close to the critical lengths X S n , X A n . We note the strong similarities with the results for the flat-profile configuration [42].
There is however a remarkable difference here: the system is dynamically unstable for any cavity length since already at X = 0 a mode with complex frequency and non-vanishing real part emerges. We will refer to this novel mode as the short-length (SL) mode. In contrast, the flat-profile BHL solution lacks the SL mode and is only dynamically unstable for cavity lengths X > X 0 , with X 0 the finite critical length at which the first dynamical instability appears [42].

B. Nonlinear spectrum
The results for the linear BdG spectrum can be understood in terms of the nonlinear spectrum of stationary GP solutions. A detailed study of the full nonlinear GP spectrum in a fine-tuned square well with V 0 = V 0 (v) was provided in Ref. [43]. We briefly revisit here the results of that work and connect them with the computation of the linear BdG spectrum of instabilities performed in the previous subsection. We also extend those results to arbitrary (non fine-tuned) values of the well amplitude V 0 .
The most relevant families of solutions for the present matter are three: symmetric BHL, asymmetric BHL (both already introduced in Sec. II) and the groundstate family. The first two families of solutions emerge precisely at critical lengths X S n , X A n , described in this limit by small-amplitude oscillations on top of the finetuned BHL solution. Thus, these families of solutions are smoothly connected to the fine-tuned BHL solution as a function of X, which explains important features of the BdG spectrum: the appearance of the symmetric Ψ S n (x) solution at X = X S n describes the appearance of the degenerate n-th lasing mode while the appearance of the asymmetric Ψ A n (x) solution at X = X A n describes the transformation of the n-th mode into a non-degenerate BdG mode, in the same fashion as for the flat-profile case [39,42]. In turn, the symmetric family can be divided into complete and incomplete solutions, depending on whether the soliton outside the well has reached its minimum or not, respectively. The incomplete solutions are more energetically favorable so in the following we will generally refer to the symmetric incomplete BHL solutions as simply symmetric BHL solutions, unless otherwise stated. The density profile for a symmetric and asymmetric BHL solution is depicted in the right panel of Fig. 2, along with that for the fine-tuned BHL solution.
The remaining relevant family of solutions is labeled as Ψ SH n (x) because they are described in terms of socalled shadow solitons [42] outside the cavity, instead of regular solitons as the previous BHL families. The wave functions Ψ SH n (x) are symmetric under parity, with n = 0, 1, 2 . . . representing the number of periods inside the well. The importance of this family of solutions is that the n = 0 solution Ψ SH 0 (x) describes the ground state of the system, which accumulates atoms in the central region in order to evaporate the horizons and become fully subsonic. Its density profile is represented in solid black line in the right panel of Fig. 2.
The ground state Ψ SH 0 (x) already emerges at X = 0, departing smoothly from the homogeneous subsonic plane wave, which is the ground state when no potential well is present. Thus, as a function of X, the ground state is smoothly connected to the homogeneous plane wave solution, while the fine-tuned BHL solution is smoothly connected to the soliton solution [see Eq. (8) and ensuing discussion]. This smooth departure is observed in the upper row of Fig. 3, where both the ground state and the fine-tuned BHL solution are represented for increasing cavity length. Similar scenarios in which nonlinear stationary solutions smoothly depart from the homogeneous and soliton solutions as an external potential is introduced were already studied in the context of the crossover from the Josephson effect to bulk superconducting flow within a Ginzburg-Landau description [45]. Therefore, the fine-tuned BHL solution is never the true ground state of the system. This contrasts with the case of the flat-profile configuration, represented in the lower row of Fig. 3, in which the homogeneous BHL solution is still the ground state for cavity lengths below the critical length X = X 0 . Above that length, the true ground state begins to smoothly depart as a function of X from the BHL homogeneous solution, which becomes dynamically unstable [42]. Returning to the fine-tuned BHL solution, its excited-state character, along with the conjecture that only the ground state is dynamically stable [39,42], explains the appearance of the dynamical instability at arbitrarily small nonzero X found in the previous subsection.
The stationary nonlinear solutions also play a key role in the time evolution of the system since they represent metastable states intercepted during the transient dynamics in which the system stays for some time before collapsing due to their unstable character [39,40]. Eventually, only the ground state is expected to be dynamically stable, and hence it is one of the natural ending points of the dynamics.
In the case of arbitrary V 0 , the situation is similar because the same families characterize the relevant set of stationary nonlinear GP solutions. The main difference with respect to the fine-tuned scenario is that the finetuned BHL solution no longer exists and we only have the symmetric and asymmetric configurations available as BHL solutions. An important limitation is that each BHL solution Ψ S n (x), Ψ A n (x) only exist in a restricted region of the parameter space (v, X, V 0 ). As a result, there are sets of values of (v, X, V 0 ) for which no initial BHL solutions can be found. The family of solutions Ψ SH n (x) is also present in the generic case. In particular, the ground state solution Ψ SH 0 (x) also appears for any nonzero length and, as a function of X, is smoothly connected to the homogeneous subsonic solution of X = 0. A complete study of the full spectrum of stationary GP solutions in a square well potential can be found in Ref. [47].

IV. TIME EVOLUTION OF REALISTIC BLACK-HOLE LASER SETUPS
In order to numerically study the time evolution of the BHL solutions described in Sec. II, we follow the procedure of Ref. [40]: we add some random noise to the initial stationary BHL solutions in order to trigger the dynamical instabilities, and let that initial condition evolve numerically with the time-dependent GP equation.

A. Short-time behavior
We first analyze the short times of the simulations, in which the instability has not yet grown enough to enter the full nonlinear regime and thus the dynamics can be still described within the BdG framework. We numerically find that all the BHL solutions are dynamically unstable for any nonzero cavity length. In particular, for the fine-tuned BHL configuration, we compare the numerical results with the predictions from the linear BdG spectrum computed in Sec. III A, observing a perfect agreement.
More specifically, and as expected from Fig. 2a, the evolution of the system is typically dominated by the largest n-th mode, appearing at the predicted critical lengths X S n with purely imaginary frequency, so the observed growth of the resulting density pattern is monotonous in time. An example of this behavior is dis- The ground state is depicted in dashed black line. First and second row: initial fine-tuned BHL solution with v = 0.5 and X = 1, so that X S 0 < X < X A 0 . Third and fourth row: same but for X = 2, so that X A 0 < X < X S 1 .
played in the first row of Fig. 4, where the monotonous growth of the n = 0 mode can be observed. However, at the critical lengths X A n = X S n+ 1 2 , the frequency of the n-th mode acquires a nonzero real part of the frequency and the unstable density pattern grows while oscillating. An example of this behavior is displayed in the third row of Fig. 4, where the oscillating growth of the n = 0 mode can be observed.

B. Long-time behavior
Once the initial instability has grown enough, the system enters the nonlinear regime and tends to approach the different possible stationary GP solutions described in Sec. III B while emitting sound waves and solitons in order to conserve energy and particle number. Nevertheless, we numerically find that all the stationary GP solutions are dynamically unstable except the ground state, so they eventually collapse and the dynamics is restarted again. This implies that there is a perfect correspondence between energetic and dynamical instability, confirming the conjecture of Refs. [39,42]. The above described dynamics is continued until sufficiently long times are explored, when only two possible behaviors are observed: the system either reaches GS by evaporating the horizons or evolves to the CES regime, where it self-oscillates in a time-periodic way. The results of this work then confirm the trends described in Ref. [40].
Representative examples of this dynamics can be seen in Fig. 4, where we examine the time evolution of finetuned BHL solutions. In the first two rows, we observe the monotonous growth of an initially unstable fine-tuned BHL solution towards the Ψ S 0 (x) stationary solution (see t = 100) while emitting a soliton downstream. However, such solution is also dynamically unstable and eventually collapses, finally reaching the stable ground state while emitting solitons and waves. In the last two rows, The red color in the lower row helps to identify the soliton whose evolution gives rise to the CES mechanism.
the results from another simulation are shown, in which the system grows while oscillating towards the stationary solution Ψ SH 1 (x) (see t = 200), which in turn is also unstable, collapsing for sufficiently long times and reaching once more the stable ground state after radiating. For longer cavities, the evolution becomes much more complicated since more nonlinear solutions are available to be intercepted during the transient. Nevertheless, for sufficiently long times, we have numerically observed that all fine-tuned BHL solutions always evolve towards GS.
In contrast, generic BHL solutions can indeed reach the CES regime, as shown in Fig. 5. In the first row, we observe an evolution along the same lines of Fig. 4: an initially unstable symmetric complete n = 0 BHL solution grows while oscillating towards Ψ S 0 (x) (see t = 400), emitting a soliton downstream in the process. Such stationary solution is again unstable and the system eventually collapses approaching the ground state. However, as seen in the lower row of Fig. 5, the soliton emitted upstream (marked in red; see also the latest times in Fig. 4) cannot travel now against the flow and is dragged back towards the supersonic region. The passage of the soliton through the square well leaves the system again in the same initial configuration that emits another soliton upstream, giving rise in this way to the periodic self-oscillating process leading to the CES state. The described mechanism is quite similar to that of the flat-profile configuration, characterized in detail in Ref. [40].
Even though the initial seed is stochastic, leading to very different possible transients between simulations, the specific long-time fate of the system (GS or CES) is solely determined by the background parameters (v, X, V 0 ) and not by the initial noise or the type (asymmetric-symmetric or the integer n) of initial BHL solution. Only very close to the phase transition from GS to CES some sensitivity to the initial condition can be sometimes observed. We attribute this sensitivity to the critical effect of fluctuations on the amplitude of the emitted soliton that originates the CES state.
This independence of the final state with respect to the transient implies that the equivalent of a phase diagram for the long-time behavior, characterized by the GS and CES states, can be drawn in the parameter space (v, X, V 0 ): at long times, for fixed values of (V 0 , X), there is a critical velocity v c (V 0 , X) above which the system enters the CES regime. Specifically, if the flow velocity v is below the critical velocity v c (V 0 , X), the soliton emitted upstream is able to escape and the system evolves towards the GS. Above that critical velocity, the soliton is trapped by the flow, giving rise to the CES state.
However, the numerical computation of the long-time phase diagram for generic BHL solutions presents several problems. First, there are regions in parameter space in which BHL solutions do not exist. Second, for large cavities, the high number of intermediate metastable states hit by the dynamics, along with the associated decay process, make the simulations computationally expensive due to the long times required. We develop in the next section an alternative approach to address the problem.

V. LONG-TIME PHASE DIAGRAM OF THE IHFC MODEL
We study here a complementary model with no horizons, represented in Fig. 1d, in which an attractive square well is suddenly introduced at t = 0 in a subsonic homogeneous flowing condensate, characterized by a wave function Ψ(x, 0) = e ivx , v < 1. In this way, a time-dependent dynamics is induced since, in the new potential profile, such initially homogeneous flowing condensate is not a stationary GP solution. We will refer to this setup as the IHFC model. The idea is that the resulting Hamiltonian governing the time evolution of the IHFC model is exactly the same as for the generic BHL solutions, determined by the set of parameters (v, X, V 0 ). The only difference is the initial condition at t = 0, which now is a non-stationary GP solution that will follow a deterministic evolution. This contrasts with the case of a generic BHL solution, where the dynamics is stochastic, driven by the initial random noise on top of the stationary BHL solution. Due to the strong insensitivity to the details of the transient reported in the previous section, one can expect the longtime phase diagram for the IHFC model to be identical to that for BHL solutions. Indeed, based on previous studies with localized defects [48], scenarios of emission of trains of solitons are already expected in this kind of setup.
Several advantages motivate the introduction of this model. First, the values of (v, X, V 0 ) can be chosen arbitrarily, without restriction. Second, since the ground state is smoothly connected to the initial homogeneous solution, the time required to explore the associated long-time dynamics is considerably shorter than for initial BHL solutions, which lowers the computational cost. Third, the evolution here is purely deterministic and there is no need to double check any sensitivity to the initial condition. Finally, we note that this model can be easily studied experimentally, and actually a similar setup was implemented in Ref. [49], in which a localized attractive potential was swept along a trapped condensate, observing emission of solitons.
As expected, in the numerical simulations it is observed that for long times the system only displays the same two possible behaviors. The transition between GS and CES for the IHFC model is shown in Fig. 6, where in the first row we represent a case in which the system reaches GS at long times and in the second and third rows we represent the same system but with a slightly higher flow velocity that puts the system in the CES regime. Interestingly, the mechanism giving rise to the CES state in the IHFC model is the same as that described in Fig.  5. We also remark that the times required to reach the long-time regime are considerably shorter than for initial BHL solutions; see Figs. 4, 5.
In Fig. 7 we represent the numerical results for the long-time phase diagram of the IHFC model. The critical velocities v c (V 0 , X) marking the phase boundaries are computed as follows: for fixed X and V 0 , the value of v is changed. The critical value v c is computed as the average between a point in parameter space representing a simulation showing CES, and one evolving towards GS, also originating in this way the associated error bar, chosen of the order of the numerical precision. The critical value v c (V 0 , 0) = 1 is extracted from the Landau criterion of superfluidity [50]. Indeed, this phase diagram may be viewed as showing a generalized nonlinear version of the Landau criterion for this potential profile.
In that phase diagram we also observe that, for increasing values of the well width X, the critical velocity saturates, v c (V 0 , X) ≃ v c (V 0 ). In the left panel of Fig. 8 we further check this feature by representing the values of the critical velocity v c (V 0 , X) as a function of V 0 for X = 2 (solid blue), X = 4 (dashed red), and X = 10 (dashed-doted green). We see that the function v c depends on X very mildly since the three curves almost coincide within numerical precision.

VI. ROBUSTNESS AND UNIVERSALITY OF THE LONG-TIME PHASE DIAGRAM
We now connect the results for the long-time behavior of the IHFC model of Sec. V with those of the BHL solutions of Sec. IV. As a first step, we represent in the left panel of Fig. 8 the curve V 0 (v) characterizing the fine-tuned BHL configuration as a dashed black line, observing that it is well below the critical velocities predicted by the IHFC model. This observation, along with the found independence of the long-time behavior from the transient, provides an explanation of the previously reported lack of CES for initial fine-tuned BHL solutions.
Moreover, when available, we compare the critical velocities for initial generic BHL solutions with those for the IHFC model. An example of this comparison is pro- vided in the right panel of Fig. 8, where we plot the critical curve for X = 2 for the IHFC model along with the critical velocities obtained using generic BHL solutions with the same background parameters as initial condition, finding a perfect agreement within the numerical uncertainties. The scarcity of points is due to the fact that, as explained, there are regions in the parameter space (v, X, V 0 ) where there are no BHL solutions. Specifically, the comparison can only be carried out if a BHL solution is present close to the critical velocity v c (V 0 , X). Since v c (V 0 , X) ≃ v c (V 0 ), the availability of a BHL solution at (v c (V 0 ), V 0 , X) for fixed X is controlled by V 0 , exhibiting a behavior in terms of allowed bands and forbidden gaps as a function of V 0 . The black dots in the right panel of Fig. 8 then correspond to parameters (v c (V 0 ), V 0 , X) within allowed bands, while the rest of the points of the blue curve are placed within gaps where no BHL solutions are available. For completeness, we have extended the IHFC-BHL comparison for other compatible values of V 0 and X such that (v c (V 0 ), V 0 , X) is within an allowed band, and in all cases a similar agreement for the critical velocities is observed (not shown).
These results imply a very strong form of robustness for the long-time regime since two completely different types of evolution, one stochastic, driven by random noise added to unstable stationary black-hole laser solutions, and another one purely deterministic, arising from the time-evolution of a homogeneous initial condition, converge to the same state for long-times. The particular realization depends solely on the background parameters (v, X, V 0 ) characterizing the global flow. This form of robustness goes even beyond that discussed in Ref. [40], which only involved linear random noise on top of the initial BHL configuration.
The observed robustness suggests that the CES regime is an intrinsic state of a flowing condensate and not a peculiar feature associated to some particular transient or phenomenon. Moreover, the IHFC configuration, which can be easily reproduced in the laboratory, could be used to study the long-time phase diagram of BHL configurations without the need to actually start from true BHL configurations, whose evolution in actual experiments is dominated by the background Bogoliubov-Cherenkov pattern originated at the inner horizon [25]. Nonetheless, due to the strong insensitivity to the transient details, even in the actual experiment we can expect the resulting long-time nonlinear dynamics to be dominated by the same kind of nonlinear solutions and phase diagram here described.
In addition, the fact that the long-time behavior of the realistic BHL model here analyzed is the same as that of the idealized flat-profile BHL model, suggests some form of universality for the final fate of a black-hole laser, with only two possible choices: either the system reaches the ground state by evaporating away the horizons or it enters the CES regime, self-oscillating periodically while radiating solitons.

VII. CONCLUSIONS AND OUTLOOK
In this work, we have studied the time evolution of a realistic model of a black-hole laser. We have confirmed and extended results from previous works dealing with the idealized flat-profile configuration. We have computed the spectrum of linear instabilities and found a behavior similar to that described in the literature except for appearance here of a novel mode that makes the system unstable for any nonzero cavity length. In the nonlinear regime, there is also a phase diagram for the long-time behavior in which the system either goes to the ground state by evaporating the horizons or to the CES regime, suggesting a form of universality for the long-time behavior of a black-hole laser. Actually, the CES state provides the closest counterpart of an actual optical laser device, as discussed in Ref. [40].
In order to reach a deeper understanding of the longtime dynamics, we have developed a complementary deterministic model, based on the sudden introduction of an attractive square well within a homogeneous subsonic flowing condensate, with no horizons present. We have studied the associated long-time phase diagram and observed a perfect agreement with that of the corresponding black-hole laser configurations for the same background flow. This implies a stronger form of robustness with respect to the initial condition, which extends beyond the BHL scenarios described in the literature.
From the gravitational analog perspective, the unambiguous observation of the black-hole laser effect still represents an ongoing challenge in the field. The present work provides a description of the expected properties of the black-hole laser dynamics as we are able to isolate it in our model from the Bogoliubov-Cherenkov background. Thus, the results here provided can serve as a basis for future theoretical studies trying to provide a distinctive signature of black-hole lasing in an actual experimental setup. Furthermore, since we have observed that the nonlinear regime is quite independent of the transient details, we also expect that the results of this work can be useful for the study of the nonlinear behavior of the stimulated Bogoliubov-Cherenkov wave.
Once the spontaneous Hawking effect has been finally observed, an important next step is the observation of the backreaction of the emitted Hawking radiation, which in a real black hole causes its evaporation, a process studied within a semiclassical framework. The nonlinear regime of a black-hole laser provides an ideal scenario to observe the analog of this effect: even though the full nonlin-ear dynamics of a condensate, described by the Gross-Pitaevskii equation, is very different to that of a real black hole, described by Einstein equations, still we can expect that the analog system can inspire some ideas on the final fate of a real astrophysical system. For instance, one may wonder what is the black-hole equivalent of the long-time phase diagram.
From a more general perspective, the current work provides a new scenario of emission of solitons, of interest not just within the condensed matter community [48,51,52] but also for other fields such as nonlinear optics [53,54] or quantum fluids of light [55][56][57], where similar nonlinear equations of motion are satisfied. The CES state represents an interesting scenario for the fields of quantum transport and atomtronics [58,59].
Several possible lines of research follow the results of this work. For instance, the strong robustness of the CES regime here described suggests that it is an intrinsic state of the background flow. Therefore, a detailed study of its properties is in order [60]. Future work should address the effect of quantum fluctuations on the evolution of the black-hole laser here described, devoting special attention to the stability of the CES state [60]. Finally, the present results can serve as the basis for a thorough analysis of the signatures of black-hole lasing, thus providing key elements for a future unambiguous identification.