Black hole lasers in Bose-Einstein condensates

We consider elongated condensates that cross twice the speed of sound. In the absence of periodic boundary conditions, the phonon spectrum possesses a discrete and finite set of complex frequency modes that induce a laser effect. This effect constitutes a dynamical instability and is due to the fact that the supersonic region acts as a resonant cavity. We numerically compute the complex frequencies and density-density correlation function. We obtain patterns with very specific signatures. In terms of the gravitational analogy, the flows we consider correspond to a pair of black hole and white hole horizons, and the laser effect can be conceived as a self-amplified Hawking radiation. This is verified by comparing the outgoing flux at early time with the standard black hole radiation.


Introduction
In 1981, Unruh suggested [1] that the analogue of black hole radiation could be observed in a quantum fluid that crosses the speed of sound, thereby forming a sonic horizon. Since then, several types of fluids and several types of flows were considered [2]. Very recently, a near-stationary supersonic flow engendering a pair of horizons (a black hole one followed by a white hole one) was realized in a Bose-Einstein condensate [3]. In such a background, because of the anomalous dispersion of Bogoliubov phonons, one expects to get a kind of laser effect [4]- [6] due to a self-amplification of the Hawking radiation.
The first aim of this paper is to provide the theoretical basis of this effect starting from the Bogoliubov-de Gennes equation, i.e. without making use of the gravitational analogy. In this we apply to flows containing two horizons the treatment of [7] that was only applied to a single (black hole or white hole) sonic horizon. Following [6] we then establish that the laser effect is governed by a discrete set of complex frequency modes that correspond to the resonant modes of the cavity formed by the region bordered by the two horizons. In a classical description, this discrete set governs the dynamical instability of the flow. The second aim is to compute the spectrum of the complex frequency modes by theoretical and numerical methods. The excellent agreement of the results validates both the concepts and the semi-classical methods used in the theoretical approach. Finally we consider two observables: the mean number of emitted phonons, and the two-point function of the phonon field that governs the density-density correlation pattern [8,9]. At sufficiently late time, both observables are governed by a single mode, the most unstable one. In this late time regime, their behaviours are identical to those one would obtain using a classical description of the density perturbations. At early time instead, when starting from vacuum configurations, correspondence with the quantum phonon flux emitted by a black hole horizon [7] is established.
The propagation of phonons in flows containing two horizons has already been considered. However, toroidal configurations with periodic boundary conditions were generally used [10,11]. In that case, the spectrum is very complicated because it results from a combination of two effects: the discreteness of the wave vectors defined on the torus interferes with that associated with modes that are trapped in the supersonic region. As a result, not only the analysis is difficult, but the relationships with the Hawking effect and the black hole laser effect [4]- [6] are hard to draw. On the contrary, when dealing with continuous wave vectors, the analysis of the complex frequency mode is simpler, and the relationship with Hawking radiation is easily made.
We have organized this paper as follows. In section 2 we present the Bogoliubov-de Gennes equation in a way that is suitable for our aims. In section 3 we determine the spectrum of asymptotically bound modes, and explain why complex frequency modes appear when the speed of sound is twice crossed. We then compute the eigen-frequencies using semi-classical methods. In section 4 we numerically solve the Bogoliubov-de Gennes equation and compute the spectrum. We then evaluate the mean occupation number and compare it with the spectrum one would obtain if only the black hole horizon were present. We finally evaluate the correlation pattern of the density-density two-point function.

Black-hole-white-hole geometries
Using the analogy [1,2] between the wave equation governing sound propagation in a moving fluid and that governing light propagation in a curved space-time, one can associate an acoustic geometry to the fluid flow. In 1+1 dimensions, these geometries are of the form [4] where v is the flow velocity and c the speed of sound. In this language, a black-holewhite-hole geometry is obtained when v crosses twice c. Indeed, a sonic horizon is present whenever v crosses c. Assuming that the fluid flows from right to left (v < 0), the location x H of a horizon is given by c(x H ) + v(x H ) = 0. In our case, we set the location of the white hole horizon and that of the black hole respectively to x W = −L and x B = L. These divide the x-axis in three regions, which we call I on the left of the white horizon (x < −L), II between the two horizons (−L < x < L) and III on the right of the black horizon (x > L). The flow is supersonic in the internal region II and subsonic otherwise. The opposite case where the velocity is subsonic inside does not give rise to a laser effect, and will be considered elsewhere.
In the present work, our analysis is restricted to the following stationary profiles: which generalize to two horizons those used in [7]. c H is the sound speed at the horizons, 2L is the distance between them, D determines the size of the near-horizon regions where the metric is not flat, n controls the sharpness of the transition to the flat regions, and κ W and κ B control the surface gravities. Indeed, the surface gravity is defined by [2] which yields −c H κ W for the white hole and c H κ B for the black hole. The metric (1) is then completely fixed by introducing an extra parameter q: which specifies how c + v is shared between c and v. When restricting to the cases where the range of D and q is limited. The allowed values are graphically illustrated by the shaded area in figure 1, left panel. Three velocity profiles corresponding to different couples (q, D) are plotted in figure 1, right panel.

Density perturbations in Bose Einstein condensates
We present the main steps leading to the equations for linear density perturbations in Bose-Einstein condensates [12], having in mind cases where the condensate crosses the speed of sound. We follow [7], where more details can be found. At low temperature, the quantum properties of a gas of weakly interacting atoms are efficiently described by a second quantized field satisfying the commutation relation and by its Hamiltonian where m is the atom mass, V the external potential and g the effective coupling. The last two quantities can depend on both t and x. When a significant fraction of the atoms condense, it is meaningful to expandΨ in a c-number function Ψ 0 , describing the condensed part, plus a field operatorφ, describing (relative) density perturbations over the condensateΨ In what follows, we consider elongated condensates, which means that the transverse excitations have sufficiently high energies that they are not excited. In this case, Ψ 0 andφ are effectively one-dimensional fields. For simplicity, we also assume that the condensate is infinitely long, in order to avoid discussing the discreteness of the longitudinal wave number k. The discrete case is indeed more complicated, and treated in details in [11]. Finally, we assume that the condensate is stationary, which means that Ψ 0 is of the form where µ is the chemical potential, ρ 0 (x) the mean density of condensed atoms and is their mean velocity. In this case, the Gross-Pitaevskii equation reduces to where the second equation is the continuity equation for a stationary flow.
At the linear level,φ obeys the Bogoliubov-de Gennes equation. In the present case, using (8) and (12), this equation reads where we introduced the x-dependent speed of sound and where reduces to the usual kinetic operator when the background condensate is homogeneous. The second expression is only valid in stationary flows. Equation (13) tells us that the perturbationφ is coupled to the condensate only through the functions v and c. The disappearance, from (13), of the potential V , the quantum potential − 2 2m ∂ 2 x √ ρ 0 ρ 0 , and the coupling g constitutes an essential step towards the notion of the acoustic metric of (1); see [7] for more details.
Using the canonical commutation relation which follows from equations (6) and (8), (13) is the Heisenberg equation engendered by the Hermitian Hamiltonian which is obtained by expanding (7) to second order inφ,φ † ; see appendix A.

Spectrum of bound modes and quantization
In stationary condensates, the spectrum of (18) can be characterized using very general properties. These are first, the fact that the eigenmodes are asymptotically bound (to guarantee that they have a well-defined norm since the spatial domain ofφ is infinite), second, the Hermiticity of (18), i.e. its self-adjointness with respect to the scalar product defining the eigen-modes norm, and third, that this scalar product is not positive definite, as it is the case here, see (B.15), and for bosonic fields, see e.g. [13]. When these conditions are met, the spectrum of (18) contains a continuous set of real frequency modes labelled by ω and a discrete index α, plus, possibly, a discrete and finite set of complex frequency modes that appear in pairs with complex conjugated frequencies. (If the scalar product were positive definite, as is the case for fermions, the spectrum of (18) would be purely real.) Throughout this paper, complex frequencies shall be written as λ a = ω a + iΓ a , where a is a positive integer which labels the discrete set of pairs, and where ω a and Γ a are both real and positive. The other cases are reached by complex conjugation and/or multiplication by −1. The discrete index α describes the subset of modes with the same real frequency. As explained below, in one dimension, it contains two or four modes depending on whether the flow is subsonic or supersonic.
As a result, when propagating in a stationary condensate,φ, obeying (16) and (17), can always be expanded (see appendix B) aŝ where the modes satisfy the normalization relations Correspondingly, the operatorsâ,b,ĉ satisfy the commutation relations All the other scalar products, and all other commutators, vanish. In particular, one has [b a ,b † a ′ ] = [ĉ a ,ĉ † a ′ ] = 0. Henceb a andĉ a are not destruction operators. In fact, the operatorsb a andĉ a form pairs [6] that characterize complex, i.e. with two degrees of freedom, harmonic oscillators that are unstable, and for which therefore there is no notion of quanta. (Concomitantly, one verifies that the norm dx ρ 0 [|ξ a | 2 − |η a | 2 ] vanishes.) It should be noticed that the eigenfrequencies associated withb a andĉ a are complex conjugated from each other (respectively λ a and λ * a ), something guaranteed by the Hermiticity of (18).
In terms of the above eigenmodes and operators, (18) readŝ The first line contains the usual sum over harmonic oscillators of real frequency and the c-number term accounting for the depletion, while the second line is due to the complex frequency modes. It should be noticed that the unusual form of the second c-number term follows from the unusual scalar product of (21). This Hamiltonian is not bounded from below due to the b † a c a terms. Hence the vacuum can no longer be defined as the ground state of H, as one might have expected since one is dealing with unstable oscillators. This gives rise to some ambiguity when choosing the initial 'vacuum' state (see section 3.4).
Eqs. (19) and (24) can be rewritten in a more familiar form by decomposing each couple (b a ,ĉ a ) into two couples of destruction/creation operatorsd a+ ,d † a+ andd a− ,d † a− , as shown in appendix B. The field becomeŝ and the Hamiltonian (24) We conclude with two remarks. Firstly, if one can describe the discrete set using thed a operators, there is a price to pay as the associated modes φ a+ , ϕ a+ , φ a− , ϕ a− are not frequency eigenmodes. Hence their time dependence cannot be factorized. Secondly, in usual circumstances, the discrete set of complex frequency modes is empty, as is the case in homogeneous and in near-homogeneous condensates. Nevertheless, when the condensate crosses twice the speed of sound, the discrete set is (generally) nonempty. To explain why, we first need to study the spectrum in homogeneous subsonic and supersonic flows, and then understand how to paste the corresponding modes when the flow crosses the speed of sound. When solving this, the dispersive properties ofφ become essential, as we now recall, and as first understood in [4].

Mode analysis
Inserting the field expansion (19) in (13), one obtains a system of two c-number equations The couple (φ λ , ϕ λ ) can be formed by either the real frequency modes (φ α ω (x), ϕ α ω (x)), or the complex frequency modes (ξ a (x), η a (x)) or (ψ a (x), ζ a (x)) and, accordingly, λ can be real or complex. Eliminating ϕ λ from the above system, one obtains [7] [ When the background quantities are independent of x, Fourier modes φ λ ∝ exp(ik λ x) are solutions of (28), provided k λ is a (possibly complex) root of the dispersion relation For ω < ω max , one clearly sees that, in supersonic flows, two extra (real) roots exist in the left lower quadrant. The most negative one is called, in the text, k (1) ω and the other one k (2) ω , so that k where Ω is the frequency in the comoving frame, and where c H is a typical value of the sound speed, we take to be that at the horizons, see section 2.1. We also introduced which gives the characteristic dispersive scale. When sending Λ → ∞, (29) becomes the relativistic equation Ω 2 = c 2 k 2 since the quartic term drops out. Similarly, (28) becomes the Euler equation governing sound waves. Instead, when keeping this term, the dispersion relation (29) possesses four roots, some of which can be complex. As we shall progressively see, the two extra roots are at the origin of the laser effect.
3.2.1. Modes in homogeneous condensates. In subsonic flows, |v| < c, for real λ = ω, two roots are real, and describe the right and left moving solutions. The other two are complex, conjugated to each other, and correspond to modes that are asymptotically growing or decaying, say to left. Only the first two should be used in (19), as the last two are not asymptotically bound. In supersonic flows, the situation is quite different. For real λ = ω, there exists a critical frequency ω max such that the four roots of (29) are [7] • ω < ω max : all real; • ω > ω max : two real and two complex ones, as in subsonic flows.
See figure 2 for a graphical solution of the dispersion relation (29).
When ω is real, the normalization of the (bound) modes can be easily worked out, as it is, up to a trivial factor, independent of the constant velocity v (because of Galilean invariance). Let us write where α spans over four values if |v| > c and ω < ω max and over two values elsewhere.
Using the mode normalization (20) from the mode equation (27) and the dispersion relation, one has where is independent of v. Equations (32) and (33) fix u α ω and v α ω except for a common phase. Taking both u α ω and v α ω real, one obtains In supersonic flows, when 0 < ω < ω max , the two modes associated with the extra real roots of (29) have negative norm. However, for −ω max < ω < 0, the corresponding modes have a positive norm. So, the fieldφ is expanded as ‡ where θ(z) is the Heaviside function.
In subsonic flows one obtains the same expression, but without the last term in each bracket since the extra real roots are absent.

3.2.2.
Modes in black-hole-white-hole geometries. To obtain the spectrum in condensates that cross twice the speed of sound, i.e. in flows defining the geometry of (2), we need to take into account the nontrivial propagation in the non-homogeneous regions localized near the two horizons, and the fact that each subsonic region is now ‡ In terms of the doubletsŴ defined in appendix B one haŝ defined on the half line and no longer on the full real axis. Then, we need to identify the number of globally defined independent modes. It should be noted that this set does not depend on the particular form of the mode equation, but mainly resides on the behaviour in the internal region II and in the asymptotic flat external regions I and III described in section 2.1. Therefore, the result is the same when studying a phonon field in a BEC, as in the present paper, or a generic scalar field with superluminal dispersion relation as in [6]. As explained in this reference, in geometries described by (2), there is a continuous double set (right-going and left-going waves) of real frequency modes, and the discrete set of complex frequency modes is generally not empty. Here we give a more mathematical argument yielding the same result.
In flows (2), two asymptotically flat and subsonic regions play the most important role in defining the set of modes, for both λ real and λ complex. For λ = ω real, there are three bounded modes: the two usual propagating modes associated with the real roots of (29) plus the mode that asymptotically decays. This mode should now be taken into account as it is bound in the subsonic region where it is defined. We will call the two sets of three modes φ i,I ω and φ i,III ω , respectively for the left (I) and right (III) asymptotic region, where the superscript i takes three values: u, v, CJ to characterize the asymptotic right-going u-mode, the left-going v-mode and the decaying CJ-mode, named after the paper of Corley and Jacobson [14]. In the internal region II, we have four modes that we call φ j,II ω with j = 1, 2, 3, 4. Since this region is bound, the four roots of (29) should be taken into account.
Therefore, any globally defined solution of frequency ω can be expanded in two forms by referring to its left, or right, asymptotic behaviour: Moreover, each partial wave φ i,I ω or φ i,III ω can be written in term of the φ j,II ω : Putting together the above equations, one obtains This system of four equations in six unknowns (the coefficients R i ω and L i ω ) has a twodimensional set of solutions, corresponding to two linearly independent modes. For instance, one can choose as independent solutions, φ u,in ω and φ v,in ω , the two right-going and left-going in-modes, defined as the solution of (40) with respectively L u = 1, R v = 0 and L u = 0, R v = 1. Similarly, one can chose the out-modes φ u,out ω and φ v,out ω defined respectively by L v = 0, R u = 1 and L v = 1, R u = 0. Using the scalar product of (20) to normalize these modes, the two sets are related by [6] φ u,in (41) Let us now move to the complex frequency case. Because all expressions are analytic in λ, the above analysis applies as such when replacing ω by λ = ω + iΓ. The novel aspects only come through the condition that the modes be asymptotically bound. Indeed, Γ > 0 implies that k u λ , the asymptotic wave number of the u real root of (29) (more precisely of its analytical continuation in Γ), acquires a positive imaginary part. Thus the mode diverges at x → −∞, unless one puts L u = 0. Similarly, on the right side, the v mode diverges at x → +∞, and imposing that it is bound requires R v = 0. However, these two conditions imply that the system (40) has only the trivial solution R i = L i = 0, except when its determinant vanishes. This condition defines an equation for the complex frequency λ that has a finite set of solutions: {λ a , a = 1, 2, ..., N}.

The semi-classical treatment of [6]
Having established the condition that defines the complex frequencies λ a , we now turn to the question of calculating them. Since the above algebraic method does not depend on the specific form of the mode equation, the semi-classical treatment of [6] applies to the present case. Therefore, we just summarize the results without going into the details. One should nevertheless pay attention to the fact that one is dealing with modes doublets (φ, ϕ) and to the normalization of these modes.
When the background is not homogeneous but v varies slowly with respect to the wavelength of the perturbation, the exact solutions are well approximated by their WKB approximation, which can be directly inferred from (35) The x-dependent wave numbers k α ω (x) are real roots of the dispersion relation (29) in a stationary inhomogeneous flow characterized by v(x) and c(x): Ignoring the quartic term, one obtains the dispersion relation of a massless field propagating in the acoustic metric (1). Equations are valid in the three regions I, II, III but not close to the horizons where the gradients of v and c cannot be neglected. The theoretical analysis can then be greatly simplified when taking into account the fact that the mixing between u and v modes is usually negligible [7,15]. (However the numerical analysis of the next section does not rely on this simplifying assumption.) In this approximation, the laser effect is completely due to the u-modes (for flows to the left v < 0). Considering the u in-mode in the internal region II, it can be written as a superposition of three WKB waves: where ϕ (i) −ω are the two negative frequency ϕ u -modes, corresponding to the two extra roots k ω being the most negative one, see figure 2. In [16] it was shown that the propagation from one horizon to the other is efficiently described by a unitary scattering of a two-component vector, whose first component represents the u-wave φ u ω and the other represents the trapped wave, described either by (ϕ representing (from 1 to 4) the scattering at the white hole horizon, the propagation between the white and the black hole horizon of φ u ω and (ϕ −ω ) * , the scattering at the black hole horizon and the propagation back to the white hole horizon of (ϕ (2) −ω ) * . These matrices are are the Hamilton-Jacobi actions governing the phase of the WKB modes, and L ω and R ω are the two turning points. Moreover, by unitarity, the parameters in U 1 and U 3 must satisfy For real frequency modes, the condition that the trapped mode be single valued imposes [6] where θ ω is real by unitarity. From this matricial equation, one obtains and therefore the coefficients of (44) are For complex frequency with Γ > 0, imposing that the incoming u branch vanishes, one obtains which implies Equations (49), (50) and (52) are the central equations we shall use to analyse the results of the numerical integration. Even though these equations have been obtained through a semi-classical reasoning, their validity goes beyond that of the WKB approximation, as shall be established by the remarkable agreement between the relations one can derive from them and the numerical results.
The main prediction concerns the relation between λ a , solutions of (52), and the behaviour of the coefficients of (50). Indeed, λ a solution of S 22 = 1, is a pole of B . Therefore, when Γ a is small enough (this is always true in our numerical situations), the pole is close to the real axis and the behaviour of B (2) ω for real frequencies ω is dominated by the contribution of the pole: from which Thus, |B ω | 2 should be well described by a sum of Lorentzians characterized by the complex frequencies λ a . In figure 3 we plot |B (2) ω | 2 obtained from numerical simulations and the corresponding fitted sum of Lorentzians. The agreement is excellent.
We conclude this analysis by showing that a semi-classical treatment furnishes approximate expressions for the complex frequencies λ a . Assuming |z ω | 2 , |w 2 ω | and |z ω w ω | are much smaller than 1, one can expand the condition S 22 = 1 in powers of these quantities. To zeroth order, one has λ a = ω a real, with ω a solution of the Bohr-Sommerfeld condition [6] where n BS ≡ a = 1, 2, . . . , N. A semi-classical treatment gives arg(α ωγω ) = −π. A more refined estimate [17] gives where Γ E is Euler's Gamma function.
To first order in |z ω | 2 , |w 2 ω |, |z ω w ω | and δλ a = δω a + iΓ a , S 22 = 1 gives −ωa + arg −ω − S −ω + argα ωγω ω=ωa , where T b ωa is the time for the trapped mode to make a full bounce. The phase ϑ a modulates the frequency imaginary part Γ a and the first order correction to its real part δω a , as shall be seen in section 4.2.2.

Observables: Phonon fluxes and density-density correlation patterns
Having identified the set of eigenmodes, we now consider observable quantities. We shall study both the phonon flux emitted by the black-hole-white-hole system, and the non-local density-density correlation pattern, as they illustrate very different aspects of the laser effect.
In what follows, we work in quantum settings because we shall assume that the state at the formation of the supersonic region is vacuum. In classical settings, the dynamical instability (the laser effect) would be present only if the initial densityperturbation possesses a non-vanishing overlap with some complex frequency mode, see section V.A. in [6]. For a perfectly adiabatic formation of the supersonic flow, the amplitudes associated with these modes would be zero. However, any small deviation from adiabaticity will trigger the instability, thereby engendering a behaviour very similar to that derived in quantum settings. In fact, as we shall see, the main difference between the quantum and the classical treatment is at early times because in the quantum vacuum all complex frequency modes contribute to the observables through the spontaneous excitation.
To compute observables, we first need to choose the quantum state. Because of the instability, the energy (24) is unbounded from below, and there is no clear definition of the vacuum. Nevertheless, if the formation of the supersonic region at time t 0 = 0 is adiabatic, and if the temperature condensate is low enough, i.e. much smaller than κ, the Heisenberg state of the system can be approximated by the state annihilated, at that initial time, by the destruction operatorsâ α ω ,d a+ ,d a− of (25). § Since the operatorŝ d a+ ,d a− are not stationary, the expectation values will not be stationary either, and will instead depend on the lapse of time since the formation at t 0 = 0. In what follows all expectation values will be computed in this state.
The two observables we shall study are related to the density fluctuationρ 1 =ρ−ρ 0 . Given the definition ofφ in (8),ρ 1 is given bŷ whereχ ≡φ +φ † is Hermitian. Its expansion in terms of operatorsâ,d iŝ where In the vacuum defined at t 0 = 0, the two-point function is Because of the complex frequency modes, it is not a function of t − t ′ only. This differs from the cases of a single black hole (or white hole) horizon, where the two-point function is stationary in vacuum [7].
To extract more physical information, it is convenient to return to the frequency eigenmodes appearing in (19). Defining (65) § We applied the standard reasoning to both theâ α ω and thed operators because Γ a , the imaginary part of the complex frequency λ a , is smaller than a tenth of ω a = Reλ a . the two-point function becomes The first line gives the (vacuum) contribution of the real frequency modes. Since these are only elastically scattered, see (41), they stay in their ground state and contribute neither to the emitted fluxes nor to the pattern of density-density correlations. The third line contains a mixed contribution of the growing σ a (x) and decaying modes ν a (x). It does not contribute to the fluxes at any time because σ a (x) vanishes in region I, and ν a (x) does it in III, see the paragraph after (41). Moreover, since it only depends on t−t ′ , it does not significantly contribute to the correlation pattern at late time. Therefore, both fluxes and correlation patterns are governed by the growing mode contribution, the first term of the second line in σ a (x)(σ a (x ′ )) * . It is clear that at late times, i.e. t − t 0 > 1/max(Γ a ), all observables are dominated by the mode with the largest Γ a . For instance, the equal time density-density correlation function is asymptotically given by The real part of the frequency ω a drops out and, as a result, the locus of the maxima does not propagate with time. This function is studied for various modes in section 4.4.
In case the initial state is not vacuum but a thermal state, the above two-point function would be multiplied by 2n a +1, where n a is the mean occupation number, which depends in a nontrivial way both on the temperature and on the blue shift effect due to the horizons (see equation (46) in [7]). Similar considerations apply to the quantity that is considered below. Had we worked in classical settings, we would have obtained a coherent wave whose behaviour in x and t is the same as that of 0|ρ 1 (t, x)ρ 1 (t ′ , x ′ )|0 at fixed x ′ , t ′ , see equation (54) in [6]. The main difference arises from the fact that the classical wave amplitude is fixed by initial conditions. It is also interesting to study the onset of the instability for earlier time, i.e. t − t 0 < 1/max(Γ a ). To this end, we consider the following quantity: It governs the probability that a density fluctuation of frequency ω be observed during the interval [0, T ] at some fixed location x in the subsonic region III, see [6]. When Γ a → 0, when the discrete set labeled a becomes continuous, and when T is sufficiently large, P bh−wh (ω, T ) behaves as in the Golden Rule: P bh−wh (ω, T ) ∝ Tn ω , i.e. the lapse of time T timesn ω , the mean occupation number of frequency ω. There can be pre-factors that depend on the strength of the coupling, the amplitude of the mode φ ω at x, if one deals with a derivative coupling.
It is more interesting to study P bh−wh (ω, T ) as a function of T for a given blackhole-white-hole geometry, and to compare it with P bh (ω, T ), the corresponding quantity computed in the case where only the black hole horizon would be present. When T × max(Γ a ) ≪ 1, the two observables behave in a very similar manner, even when only a few (say 10) complex frequency modes exist in the black-hole-white-hole case. This shall be seen by a numerical analysis in section 4.3, but can be also understood from the properties of the growing modes σ a (x), see the central paragraph of section V.B.5 in [6].
For classical settings, the laser effect is present only when the initial densityperturbation profile is characterized by a non-vanishing amplitude of some complex frequency mode. Since for adiabatic formation of the horizons all the amplitudes associated with those modes are zero, no laser effect appears in this case. However, even a small deviation from perfect adiabaticity causes the presence of instabilities. The main difference in the quantum framework is therefore the possibility to make the system unstable through the spontaneous excitation of complex frequency modes.

The method
In a few words, we describe how we proceeded to solve numerically (27) in the flows described by (2) and (4). As in the case of a single black hole or white hole horizon, when working with a fixed frequency, the main task is to avoid the growing mode contaminating the oscillatory modes. The way to get rid of this difficulty is by integrating from the subsonic region towards the horizon into the supersonic region where there is no growing mode [14].
Basically we solved separately the mode equation from region I to region II, and from region III to region II, using a code adapted from [7]. In each case, as in that reference, we integrated the equations for initial conditions describing the three acceptable modes in regions I and III: the right-moving u-mode, the left-moving v-mode and the decaying mode. Then for each of them we computed the amplitudes of the four oscillatory modes in the supersonic region II. The two globally defined modes are built using the procedure described in section 3.2.2. To obtain the scattering from I to III, we eliminated the amplitudes in region II.
To obtain the correlation patterns of figures 11-13 we solved the mode equation with the corresponding complex frequency λ a that we had formerly computed.

The discrete set of complex frequencies
Two strategies can be used to determine the complex frequencies λ a = ω a +iΓ a , solutions of (52). They can be obtained either by solving directly the linear system (40), or by determining the center and the width of the Lorentzians in |B (2) ω | 2 of (54). We use the latter to localize them and the former to refine the results and study how they depend on the various parameters of the system.
To give a first idea, in figure 3 we represent |B (2) ω | 2 as a function of ω/ω max (green points) for two different values of L, the distance between the horizons. A sum of Lorentzians functions (red line) is fitted to the data obtained from the numerical analysis (see (54)). The quality of the fit confirms the correctness of the theoretical analysis of section 3.3 and of [6]. In the following sections, we study the dependence of λ a on both the parameters describing the geometry (see the velocity profile (2)) and the dispersive scale Λ of (29). The dependence on n is not reported here, since it induces no significant change. Before proceeding we make some comment about the properties and the validity of the Bohr-Sommerfeld condition (55). (55) contains the difference of two wave vectors: k

The validity of the semi-classical approximation. The action appearing in
ω . These have the same sign (negative for positive ω), as they belong to the same u-branch of the dispersion relation (see figure 2), but they have opposite group velocity: k (1) ω describes a right-going mode with respect to the lab, whereas k (2) ω a left-moving one. This peculiarity give rise to an unusual phenomenon. In the usual case, eigenmodes with small (large) frequency correspond to small (large) Bohr-Sommerfeld numbers n BS , i.e to modes with few (many) nodes. In the present case instead, a large n BS corresponds to low frequency modes and vice versa. This can be understood from (55). Because the wave vectors appeared subtracted, a small n BS implies k (1) ω close to k (2) ω and this happens for ω close to ω max (see figure 2). On the contrary, a large n BS implies a large difference between k (1) ω and k (2) ω , that is, ω close to 0. This has an unusual consequence. The Bohr-Sommerfeld condition is more reliable when the action is large, that is for n BS ≫ 1. In the present case this happens for small ω/ω max . On the other hand the WKB approximation is expected to fail for small ω, as can be verified by the fact that f of (56) cannot be neglected when ω/κ ≤ 1. Both these expectations are confirmed in figure 4 where we compare the numerical results with the predictions obtained with the standard WKB approximation (green lines) and with the improved method (red lines), that is when using (56). Firstly, the agreement between the Bohr-Sommerfeld condition and the numerical results is worse for ω close to ω max and gets better when n BS increases. Secondly, at low frequency, while the quality of the The role of n of (2) is to govern the smoothness of the transition from the near-horizon region to the flat asymptotic ones. A smaller n corresponds to smoother transition, while a larger n corresponds to steeper transition and, consequently, to a larger almost-flat region between the horizons. The usefulness of n is essentially technical: it allows one to control the numerical analysis when L, the distance between the horizons, is comparable to D, the width of the transition between regions I-III of section 2.1.  standard WKB prediction (green lines) becomes worse, the improved method continues to work very well, thereby establishing its validity. Further studies about the agreement of the improved method and numerical results are in preparation [17].

4.2.2.
The birth of modes as a function of the distance between the horizons. In figure 5, left panel, we plot ω a = Re(λ a ) for the entire spectrum of complex frequency modes as a function of the half-distance between the horizons L, all other quantities being fixed. On the right panel, Γ a = Im(λ a ) is plotted as a function of L, for n BS =4 and for the same values of the other parameters. As conjectured in [6], when L grows, all ω a (L) values increase, and when there is enough 'room', a new eigenmode appears with ω a ∼ 0. As a result, the eigenfrequencies become denser close to ω max , because they neither cross it nor disappear.
It is interesting to consider more closely the birth of new modes near ω = 0. A rather difficult question to answer is the following: when extra mode appear, does the imaginary part of the frequency Γ a vanish, as found in figure 1 of [18], or not? This question is physically relevant in that it governs the continuous character of the stability of the system: if Γ a vanishes, it implies that the birth of new modes leads to a continuous behaviour in L (as one can expect on the general ground that the properties of the eigenmodes, solutions of (13), continuously depend on the parameters appearing in that equation). From figure 5, which is based on the roots of the determinant in (40), no definite conclusion can be drawn since there is a loss of numerical precision for ω → 0. A more reliable method consists in studying the behaviour of |B (2) ω | 2 on the real frequency axis. We created a gif-animation (eigenfrequencies.gif (available from stacks.iop.org/NJP/12/095015/mmedia)) of |B (2) ω | 2 as a function of ω for increasing values of L to illustrate the appearance of the new eigenfrequencies. From this it is quite clear that when a new eigenmode appears, both the real part and the imaginary part of λ a vanish, as was also found in the appendix of [19].
Another important feature confirming the theoretical analysis of section 3.3 is the presence of superposed oscillations on the overall trend of ω a . In fact, the latter is given by the solution of the Bohr-Sommerfeld equation (55), whereas the oscillations are due to the deviations governed by the sine in δω a of (58). For the imaginary part instead, there is no zeroth-order contribution, and the dominant contribution to the right plot of figure 5 is given in (57).

Asymmetric velocity profiles.
Using asymmetric velocity profiles, namely different surface gravities κ W and κ B , does not lead to substantial modifications of the spectrum. The only observable change concerns Γ a . When κ W = κ B , Γ a vanishes for some particular combination of parameters, as can be seen in figure 5, right panel. Equation (57) can indeed be rewritten as which shows that when |z ωa | = |w ωa |, Γ a vanishes when the cosine does, and this is because the scatterings at the black and the white horizons destructively interfere with each other. When the two surface gravities are different, (|z ωa | − |w ωa |) 2 is non-zero and Γ a can no longer vanish (see figure 6). The non-zero offset is significant only for sufficiently large values of ω/κ, which requires that ω max /κ be large enough, since all  ω a < ω max . Using the standard expressions z ω = e −πω/κ W and w ω = e −πω/κ B , which furnish reliable estimates [7], the ratio between the offset and the amplitude of the oscillation is Even though this function increases in ω when κ W = κ B , the consequences of this are never important because when the ratio is large, the corresponding mode will not significantly contribute to the instability since Γ a ≪ 1, as can be seen from (69).

4.2.4.
The u-v mixing and the parameter q. The mixing between the u-modes, which are right-going in the lab in subsonic flows, and the v-modes is controlled by the transmission coefficient T ω of (41). In figure 7, |T ω | 2 is plotted as a function of q for a fixed frequency ω/κ = 0.100. The transmission coefficient is almost 1 between q = 0.25 and q = 0.75. Thus, for values of q in this range the u-v mixing is negligible, as was noticed in [7] for a single black (or white) hole. In this regime the approximation discussed in section 3.3 is valid. There is, however, an important difference with respect to the single black hole case. In black-hole-white-hole geometries, the transmission coefficient deviates from zero when approaching a resonance. When working with a fixed ω, we found two sharp spikes figure 7 for q = 0.1248 and q = 0.56023 which exactly correspond to two unstable modes. Moreover, the width of the spikes is almost equal to the corresponding Γ a . We also looked for complex eigenfrequency associated with the other two local and broad minima of |T | 2 , at q = 0.3363 and q = 0.6285, and we found two eigenmodes with respectively ω a /κ = 0.090 and ω a /κ = 0.103, i.e. in the neighbourhood of the chosen frequency ω/κ = 0.100. These considerations show that close to resonances one cannot completely neglect the u-v mixing. Nevertheless, as shown in figure 4, the improved treatment of the Bohr-Sommerfeld equation produces very good estimates for the values of the real part ω a of the complex frequencies.

4.2.5.
The role of the maximal frequency ω max . In [7,15], it was shown that the deviations, due to dispersion and w.r.t. the standard Planckian distribution, of the spectrum emitted by a single BH, or WH, are mainly governed by ω max . However, the latter depends both on the UV scale Λ and the velocity profile (see section 3.2.2). The relationship takes the form Hence the same value of ω max can be reached from very different cases, and yet it was found that the fluxes are hardly sensitive to this. In the present case however, this insensitivity is lost because the number of resonances directly depends on Λ. This can be understood by studying (55), and is manifest in figure 8, where |B (2) ω | 2 is plotted as a function of ω/κ for two different values of Λ.

Growth of the asymptotic phonon fluxes
In figure 9, we represent the relevant quantities that govern the properties of the 'probability' function introduced in (68). On the left upper plot, we present the growth rates Γ a and the corresponding times T b ωa of (60) for the 13 complex frequency modes that exist in the condensate flow of figure 3, lower plot. On the right upper plot, we show the Golden Rule transition rate associated with (68) in that black-hole-white-hole geometry (dots), and the rate one would obtain if the white hole were not present. In that case, the quantity plotted is proportional to ω ×n ω , wheren ω is the mean occupation number in the black hole geometry as computed in [7]. It is given byn ω = |w ω | 2 /(1 − |w ω | 2 ), see U 3 of (46), where |w ω | 2 is well approximated by ωa of (60) as a function of ω/ω max in unit of κ −1 (points) and corresponding values of Γ a (boxes), when Lκ/c H = 25. The condensate parameters are as in figure 3, lower plot. Upper right panel, solid line: transition rate associated with (68) due to the radiation emitted by the sole black hole as a function of ω/ω max (solid line); corresponding quantity due to the discrete set of complex frequency modes in the black-hole-white-hole geometry. The oscillations are due to the cosine in (57). Lower panels: P (ω, T ) of (68) as a function of ω/ω max for an isolated black hole (solid line) and for a black-hole-white-hole pair (dashed line), respectively at time T = 30κ −1 (left plot) and T = 200κ −1 (right plot).
The pre-factor of ω is due to the normalization of the χ modes of (63). In a black-holewhite-hole, the quantity that corresponds to ω ×n ω is ω a × 2Γ a T b ωa . In both cases, these quantities vanish for ω > ω max .
In the lower plots, we show (68), at two different times, and for both the blackhole-white-hole, and the isolated black hole flows. As expected, the exponential growth of the laser effect does not show up for times T smaller than the inverse of the maximal Γ a (here T ∼ 110/κ). Moreover, the discreteness of the spectrum is not visible either at early times. To be resolved, it requires times larger than 2π∆ω a , where ∆ω a is the frequency gap between neighboring Bohr-Sommerfeld frequencies ω a . It is here of the order of ∆ω a ∼ ω max /10 ∼ 0.02κ.
On the contrary, for times of the order of 1/max(Γ a ) or greater, both the exponential growth of the laser effect and the discreteness of the spectrum show up. These allow one to distinguish the phonon flux emitted by the black-hole-white-hole pair from that emitted by the sole black hole that grows linearly in T for all values of ω. This linear growth can be observed by comparing the continuous lines of the right upper and lower plots, and constitutes a numerical validation of the Golden Rule!

Spatial properties of complex frequency modes and correlation patterns
The density-density correlation function can be calculated following the procedure described in section 3.4. As outlined there, only the largest-Γ a mode significantly contributes to the pattern at late time. Nevertheless, to appreciate the variety of cases, it is worth studying the modes and the corresponding correlation patterns also for lower Γ.
As a typical example of a mode with a large Γ, we consider the fourth BS mode in a configuration with fivecomplex eigenfrequencies. In figure 10 (upper panel) the real part of this mode is represented, while the supplementary movie laser 1.gif (available from stacks.iop.org/NJP/12/095015/mmedia) gives its evolution in time. As a second example, we choose a very different situation with q = 0.7 and a dispersion relation with higher Λ. In this case, the spectrum of complex frequencies is large. In figure 10 two modes are plotted: one with n BS = 14 and a moderate value of Γ (central panel), and one with n BS = 2 and a very small value of Γ (lower panel). When comparing the modes and their correlation patterns, some features remain the same, whereas others significantly differ.
In the first example, since Γ is high, the laser effect is strong as can be seen from the movie: in the right asymptotic region there is an exponentially growing right-going flux. This growth in time can also be inferred from the spatial decrease of the mode amplitude on the right of the black hole, because these are linearly related, see equation (54) of [6]. Considering the correlation pattern of this mode in figure 11, one sees that the dominant signals come from the inside region (the central square of size ∼ 22) where the amplitude of the mode is highest, and from the horizontal and vertical bands which describe the correlations between the escaping flux and the trapped mode.   From the mode itself and from its correlation pattern in figure 11, one can see that the u-v mixing is very low, as expected from the choice of q = 0.5. The leftgoing mode coming out of the white hole is about 2.5 orders of magnitudes smaller than the right-going one that comes out of the black hole, as can be seen from the ratio of amplitudes between the top left part and top right part in figure 11. To understand how the modes propagate, it is appropriate to consider the top left part. This region represents the correlations between u (right-going) and v (left-going) modes. Figure 13. Density-density correlation function (67) divided by e 2Γat for the mode with a low Γ represented in the lowest panel of figure 10. The bottom legend corresponds to the central region of the plot. One clearly sees the two nodes in the trapped region. One also finds them in the inside-outside correlations, both on the left and on the right of the trapped region. Finally, one notices that in the present case the u − u correlations on the right upper square are weaker than the v − v correlations on the lower left panel. This follows from the high u-v coupling between the trapped mode, which is a u-mode, and the outside v-mode.
The slope of the highest/lowest correlation lines gives the ratio between the velocities of the two modes. When the modes form a continuous set, the two velocities are the group velocities. However, in the present case, when considering a single mode, the slope of these lines corresponds instead to the ratio between the phase velocities of uand v-modes. Similarly, the pattern in the top right part are simply due to the nodes of the right moving mode evaluated at a given time, see equation (61) in [6].
We now consider the correlation pattern of figure 12, corresponding to a mode with a moderate value of Γ and n BS very large. The overall properties of the pattern are basically the same as those of figure 11. The main differences come from the high wave-vector content of the trapped mode, which produces short distance features in the central square and in the bands. Outside these regions the patterns are very similar.
In the last pattern of figure 13, corresponding to a mode with a very small value of Γ and n BS = 2, we get a different picture. The pattern is basically concentrated in the central square since the amplitude of the mode outside is related to its amplitude inside by the square root of Γ, as can be seen from (57). In the present case, the u-v mixing is so high that the pattern on the left of the white hole is higher (by a factor of about 5) than on the right of the black hole. One also notices that in the central panel (n BS = 2) there is only one node. Moreover, the beat-like shape shows that this wave is the result of the superposition of two waves with very similar and relatively high wave vectors.

The Technion experiment
Recently a black-hole-white-hole flow has been experimentally realized [3]. Our numerical code can describe this configuration, up to some limitation. In fact, the velocity profile (figure 3(b) of [3]) is rather irregular while our code works accurately only for velocity profiles as in figure 1, right panel, with an almost-flat internal region. Nevertheless, we shall proceed in order to obtain at least an estimation on the number of the unstable modes and their growing rate. As far as we know, this is the first estimate of this quantity, which is very relevant from an experimental standpoint.
The main parameters of the experimental setting are To reproduce the above setting we perform two simulations with fixed D = 1. In figure 14, |B (2) ω | 2 is plotted as a function of ω/ω max for q = 0.9 (left panel) and q = 0.48 (right panel). These two values are obtained by using, respectively, c(x = 0) and c(x → ∞) equal to their experimental values in (73) and (74). ¶ The small distance between the two horizon implies that few trapped modes are present in this situation. The values of Γ a and ω a of the complex eigenfrequencies are reported in both cases in table 1. The dominant contribution comes from the eigenfrequency with the largest ¶ We used these two values because, in our analysis, it is not possible to fix independently c(x = 0) and c(x → ∞). We expect that the actual value of the complex frequencies would lie in between the two sets we obtain. ω | 2 as a function of ω/ω max for q = 0.9, ω max /κ ≈ 3.4 (left panel) and q = 0.48, ω max /κ ≈ 2.7 (right panel). The values of the other parameters are given in the text. Green points: numerical simulation; red lines: fitted series of Lorentzians. ω a and Γ a from the fit are reported in table 1. Table 1. ω a and Γ a for the complex eigenfrequencies in the systems described in the caption of figure 14. Γ a . In the two cases this corresponds to an instability time scale, respectively τ ≈ 0.06 s and τ ≈ 0.02 s. Even if the two simulations are generated with very different values of q, the two results agree in order of magnitude. Note that the estimated instability time scale is more than twice the time scale for which the horizons are maintained in the experimental configuration (≈ 0.008 s). In order to see the laser effect, it would be needed to maintain the systems for a longer time, or to increase κ by a factor of 10, without changing the adimensional parameters of the system such as the ratio Lκ/c H .

Conclusions
In this paper, we analysed the spectrum of phonons in a BEC when the stationary flow of the condensate crosses twice the speed of sound. Our analysis is based on the Bogoliubov-de Gennes equation (13). Hence, even though the analogy with light propagation in a pair of black and white horizons is manifest, none of our results rely on this gravitational analogy.
In the limit where the condensate can be considered as infinite, i.e. when no periodic boundary condition is introduced, the spectrum of bound modes contains real frequency modes that are only elastically scattered, see (41), plus a discrete and finite set of pairs of complex frequency modes (19). These modes can be seen as the resonances of the cavity bordered by the two sonic horizons. They lead to dynamical instabilities because the scattering through each sonic horizon is anomalous, in that it mixes positive and negative norm modes, see (44). Then, using semi-classical techniques, we compute the real and the imaginary part of the complex frequencies. The real part ω a obeys a Bohr-Sommerfeld condition that introduces the discreteness of the spectrum, whereas the imaginary part Γ a is related to the norm of the scattering coefficients across the horizons, see (57). In the case of toroidal geometry, the analysis would instead be complicated by the discreteness of the wave vectors. The instability appears only when one ω a approximately matches one of the frequencies associated with that discrete set of wave vectors. This explains the presence of the instability bands found in [10,11].
In section 4, we numerically solved the Bogoliubov-de Gennes equation for real and complex frequencies. By comparing the numerical properties with those derived using the Bohr-Sommerfeld and (44), we validated the use of semi-classical methods, see figures 3-5. In section 4.3 we numerically studied the growth of number of phonons emitted by the black-hole-white-hole system. In spite of the discrete character of the trapped modes, at early time, this quantity behaves very similarly to the flux that the black hole horizon would emit in the absence of the white hole horizon. Instead, for larger time both the instability and the discreteness show up. Finally, we studied the equal time correlation pattern of density-density fluctuations associated with three different complex frequency modes, respectively with high, moderate and low value of the imaginary part Γ a . In all cases, the instability shows up most clearly in the supersonic region where the amplitude of the trapped mode is the highest one. The correlations between the trapped mode and the emitted modes display very specific properties, see figure 11 and then next two ones.
Finally, we studied (within some limitations) the experimental situation realized in June 2009 in the Technion [3]. The results are summarized in figure 14 and show that few unstable modes are found, and that the instability time scale is about ten times larger than the life time of the condensate. This means that an increase by a factor of ten of the surface gravity (without changing the adimensional parameters of the system such as the ratio Lκ/c H ) would lead to comparable time scales. The laser effect could then be observable. the structure ofŴ must bê W = n (W nân +W nâ † n ), (B.14) where W n are doublets of C-functions, and n denotes the summation over a (possibly continuous) complete sets of modes. Using ∂ T x = −∂ x , T T ρ = T ρ and the properties of Pauli matrices, one verifies that the scalar product is conserved under time evolution when W i are solution of (B.10), since For later convenience, we state here some properties of the scalar product that can be verified using the anticommutation relation of Pauli matrices and (16): Assuming that the W n have a non-zero norm, we define an orthonormal basis: which shows thatâ n andâ † n are in fact destruction and creation operators. When the condensate is stationary, one can work with eigenmodes W α λ of frequency λ, where the index α describes the set of modes with the same frequency. By definition, one has BW α λ = λW α λ . (B.24) The conservation of the scalar product gives From this it is clear that only real frequency modes can be normalized as in (B.21). However, in the presence of dynamical instabilities, complex frequency eigenmodes are present. We shall assume that these modes form a discrete and finite set of pairs of modes with conjugated frequencies that we call {λ a = ω a + iΓ a , a = 1, 2, ...N}. In place of (B.21), we choose the following pseudo-normalization for each pair, W λa |W λ * a ′ = iδ aa ′ , (B.26) and the other scalar product must vanish because of (B.25). In fact, the Hermiticity of H implies that if λ a is an eigenfrequency, then λ * a is an eigenfrequency too [16]. We call V λa and Z λa the doublets corresponding to λ a and λ * a , respectively. Summarizing, the fieldŴ can be expanded aŝ