Resonant Hawking radiation in Bose-Einstein condensates

We study double-barrier interfaces separating regions of asymptotically subsonic and supersonic flow of Bose condensed atoms. These setups contain at least one black hole sonic horizon from which the analog of Hawking radiation should be generated and emitted against the flow in the subsonic region. Multiple coherent scattering by the double-barrier structure strongly modulates the transmission probability of phonons, rendering it very sensitive to their frequency. As a result, resonant tunneling occurs with high probability within a few narrow frequency intervals. This gives rise to highly non-thermal spectra with sharp peaks. We find that these peaks are mostly associated to decaying resonances and only occasionally to dynamical instabilities. Even at achievable nonzero temperatures, the radiation peaks can be dominated by the spontaneous emission, i.e. enhanced zero-point fluctuations, and not, as often in analog models, by stimulated emission.

3 extremely difficult. An alternative proposal to indirectly measure Hawking radiation relies on the squeezed character of the state [6,7], which could be observed at currently attainable temperatures by counting atoms on both sides of the horizon at coinciding times. However, it should be pointed out that the main contributions to these correlations are due to the stimulated amplification of pre-existing phonons [7,8], and not to the spontaneous amplification of vacuum fluctuations.
Here we study a new method that specifically aims at detecting the spontaneous contribution to Hawking radiation. This approach relies on the strong frequency dependence of resonant tunneling through a double-barrier structure. Such a sonic BH analogue behaves as a Fabry-Perot resonator for quasi-particles, with the peculiar feature that quasi-particles propagate linearly against a condensate background, which is itself governed by a nonlinear equation. The Hawking emitted phonon spectrum shows peaks at frequencies different from zero. We find that at currently achievable temperatures, thermal noise could be weak enough not to blur the characteristics of this resonant radiation, and a time-of-flight (TOF) experiment could allow for its detection. Our proposal is quite similar to the BH laser setup [9,10,12,13] in that sharp peaks are found in both cases. However, although dynamical instabilities may appear occasionally, they are not a necessary feature of these types of setups (see section 5). We propose here to focus on situations that are dynamically stable, i.e. where all peaks are due to resonances. The dynamical stability of the flow is likely to be a valuable asset in actual experiments. A systematic study of resonances and instabilities is currently under investigation.
This paper is organized as follows. In section 2, we present a mean-field study of the considered setup, discuss the main results and some general features. Section 3 is devoted to formulating the scattering problem of Bogoliubov quasi-particles propagating against the condensate background and to identifying the essential features of Hawking radiation. In section 4, we present and discuss the numerical results obtained for the Hawking radiation spectrum emitted from a double delta-barrier interface separating a subsonic and a supersonic region. Section 5 deals with the general distinction between quasi-normal modes (QNMs; or resonances) and dynamical instabilities, both of which are candidates to explain the sharp peaks in the radiation spectrum. Finally, section 6 is devoted to a summary and discussion of the main results. The main text is complemented by two appendices. Appendix A presents the analytical calculation of the mean-field model. Appendix B deals with the analytical resolution of the quasi-particle eigenvalue problem in the inhomogeneous region on the subsonic and supersonic sides near the double-barrier interface.

Formulation of the model: a condensate wave function
We study an atom transport setup that is schematically depicted in figure 1. A quasi-onedimensional (1D) bosonic condensate, which occupies the left region x < 0, is allowed to leak to the right through two identical delta potential barriers. The leftmost barrier is conventionally placed at x = 0 separated by a distance d from the second barrier. We will see that this doublebarrier setup behaves as a resonant structure. For convenience, we neglect quantum fluctuations of the condensate stemming from its 1D character [18,19].
We may decompose the Heisenberg second-quantized field operator (x, t) = e −iµt/h 0 (x) + δˆ (x, t) A typical BH structure discussed in this paper. The two vertical arrows represent the delta barriers. The profiles of the speed of sound (proportional to √ ρ(x), where ρ(x) is the local condensate density) and the speed of flow (proportional to 1/ρ(x)) are shown. The condensate flows to the right, while Hawking radiation is emitted into the left (subsonic) region. In this particular case, three event horizons exist since the two curves intersect at three points. The flow regime is such that the condensate density exhibits a depression between the barriers. into a stationary condensate wave function e −iµt/h 0 (x), with µ being the chemical potential and δˆ (x, t) its fluctuations. In this section, we focus on the condensate behavior. At low temperatures and densities, the mean-field equation that governs the stationary flow of a Bose-Einstein condensate is the time-independent Gross-Pitaevskii (GP) equation for the condensate wave function: We wish to study the effect of a potential consisting of two delta barriers of equal strength, although for comparison we will occasionally consider the single-barrier case (see [20]). Thus, we will assume an external potential V ext (x) that takes one of the following two forms: where c u = √ g u n u /m, with n u ≡ lim x→−∞ | 0 (x)| 2 , is the speed of sound on the subsonic, upstream (x < 0) side, and z is the dimensionless strength of each barrier. The healing length on the asymptotic upstream side is ξ u =h/mc u .
In both the single-and double-barrier cases, solutions can be found in which the condensate velocity is supersonic on the right of the second barrier and subsonic from −∞ to some point in the vicinity of the leftmost (x = 0) barrier (see appendix A for a general discussion of these solutions). Those flow profiles must have one or more horizons, defined as points where the local condensate velocity equals the local speed of sound. This subsonic-supersonic scenario is the same that has been shown to display the bosonic analogue of Andreev reflection, where the supersonic side plays the role of the normal fluid [11]. Thus, one may generally speak of Andreev-Hawking processes when dealing with scattering events undergone by elementary excitations of condensate flow through a subsonic-supersonic interface.
Such sonic analogues relying on condensate flow have a superluminal dispersion relation at higher frequencies and their horizons do not imply a strict causal disconnection among regions, but they are sufficient to produce Hawking radiation analogues (see, however, [21] for a discussion of a scenario where horizons would not be needed). For the case of a single delta barrier the only horizon lies on the near left of the barrier. The double-barrier case is richer: there appears one horizon in the vicinity of the x = 0 barrier (which may lie on either side) and possibly one or more horizon pairs. The leftmost horizon, and the only one in the single-delta barrier case, can be viewed as a BH analogue, because there the flow goes from the subsonic to the supersonic side. The possible additional pairs of horizons appearing in a double-barrier structure are the analogues of white hole (WH)/BH pairs, where a WH is the time-reversed version of a BH. WH analogues seem to be extremely difficult to generate experimentally. There is a debate in the literature on whether WHs are stable at all, a question to which linear stability analysis has not yet provided an answer [8,22].
To clearly identify what constitutes Hawking radiation, the experiment should ideally be done in such a way that the condensate wave function is stationary and asymptotically flat: [23]. In appendix A, we show that with those homogenous boundary conditions, there is one or more mean-field solutions for two delta potential barriers, and just one for a single delta barrier. A sketch of a typical density for two delta barriers is shown in figure 1.
In figure 2, we plot the region of the (z, q) plane for which stationary solutions with homogeneous boundary conditions exist to the double-barrier problem. We recall that z is the dimensionless parameter characterizing the strength of the two identical delta barriers, while q is the condensate momentum on the subsonic (upstream) side. We note that, for a given value of q, there is a minimum and a maximum barrier strength between which a solution is guaranteed to exist for that particular value of q. This contrasts with the behavior of the single-barrier case, for which only one solution exists for a given value of q. Interestingly, the resulting line z(q) line for the single-barrier case coincides with the upper boundary of the shaded region in figure 2. Figure 3 shows a representative sample of solutions in the (d, z) plane, where d is the distance between barriers. Each curve represents a value of the condensate momentum q. Inspection of figure 3 reveals that, for a given q, there is a finite interval of allowed z values. That range has been shown in figure 2. Figure 3 reveals that, for given z and q, a multiplicity of d solutions exist. The two lowest d values can be characterized by two different values of ρ(0), which we call ρ max and ρ min . The upper d solutions come also in pairs and are regularly separated by a distance difference equal to the period of the nonlinear oscillations between the two barriers, which is the same for both ρ max and ρ min .
An interesting feature is that, at small barrier separations, a single q solution exists for a given d and z. For higher d (slightly above 2 for the sample of curves shown in figure 3),  two q solutions exist for given d and z. For still higher d, three q solutions exist for given d and z, and so forth. Since the barrier separation d and the barrier strength z are expected to be experimentally adjustable parameters, the clear trend is that, the larger the distance, the higher the number of allowed stationary solutions, a fact that could translate itself into instabilities.

Bogoliubov analysis
A general and detailed introduction to the Hawking radiation physics in Bose-Einstein condensates, within a Bogoliubov-de Gennes (BdG) description, can be found in [8,14,24,25]. In this section we mainly introduce the notation; see those works for a more complete presentation.
The quantum fluctuation part introduced in equation (1), δˆ (x, t), can be subject to a canonical transformation resulting in the expansion where u i (x) and v * i (x) are the components of the wave function of the bosonic quasi-particle created byγ † i . These u, v components satisfy the bosonic BdG equations: while theγ ,γ † operators satisfy, for real ω, where the normalization ν i can be set to ±1.
In equation (5), g 1D = 2h 2 a 0 /(ma 2 ⊥ ) is the effective 1D coupling constant, where a 0 is the s-wave scattering length and a ⊥ the transverse confinement harmonic oscillator length, and V ext (x) is an externally imposed potential. We note that the sum over i in the first equation may be interpreted as including both ω i ≷ 0, with ν i > 0, or both ν i ≷ 0, with ω i > 0.
The role of the horizon is appreciated when doing a WKB-type (Wentzel, Brillouin and Kramers) approximation for the excitations (see [24,25]). Close to the horizon are the turning points of the low-energy classical trajectories where the WKB solution is not appropriate. By matching the solutions on both sides and assuming some general scale separation (see also [26]), one can show that in the case of just one horizon the profile of emitted Hawking radiation approaches 1/ω at low frequencies. This important result guarantees the (approximately) thermal radiation profile. Thus, at low frequencies, and in the absence of other scattering obstacles, the zero-point radiation has a 1/ω behavior that makes it in principle difficult to distinguish it from truly thermal behavior.
Due to the presence of delta barriers in our chosen configuration, WKB-type approximations cannot be used uncritically. However, thanks to a theorem on dark-soliton perturbation theory (see [27]), we know that we have at our disposal a complete set of solutions on the left of the x = 0 barrier. In the flat supersonic region (x > d), the solutions are even simpler to work out. See appendix B for both cases, x < 0 and x > d. The only non-analytically solvable problem lies between barriers (0 < x < d), but this is a finite region where numerical integration of the BdG equations (5) requires only a moderate computational effort. The next step is to identify the relevant scattering states. This has been done for Bose-Einstein condensates in studies on BH analogues [8,14] and on bosonic Andreev reflection [11]. In this paper, we focus on the scattering states and their connection to QNMs (or  (7)) on the subsonic (left) and supersonic (right) sides. The blue/red branches correspond to positive/negative normalization as defined in equation (6). resonances) and dynamical instabilities. A systematic study of complex-energy eigenmodes is left for a future work [28]. In the asymptotic regions, propagating modes obey the Bogoliubov dispersion relation. Following [14], we label the asymptotic regions with indices u for upstream (x → −∞) and d for downstream (x → ∞). The upstream dispersion relation is where v u is the upstream flow velocity, and similarly for downstream ω d (k). A graph of this relation is shown in figure 4 for both the subsonic and supersonic sides. The branches shown in blue/red correspond to the +/− of equation (7) and can be shown to lead to positive/negative normalization. Modes are named after the sign of their group velocity (in/out according to whether they approach/leave the scattering structure), location (u/d) and, in the supersonic case, 1-2 stands for modes with normal/anomalous (i.e. positive/negative) normalization. Importantly, the anomalous d2 modes exist only for frequencies ω < ω max .
Here, we are mainly interested in the scattering state with frequency ω characterized by the incoming channel d2-in. Its wave function reads where where q u and q d are the upstream and downstream condensate momenta, w = dω/dk denotes the relevant group velocities, and the S-matrix elements are shown. In equations (8)-(10), u α (k) and v α (k) are the spinor components for a quasi-particle of the uniform Bose gas, in channel α at momentum k. A similar decomposition can be made for the other in-modes. We note that the ks of the previous equation are solutions of equation (7) for a given ω lying in the interval 0 < ω < ω max , each solution defining a particular mode. Not shown in the previous equation but necessary for the matching, an evanescent solution exists on the subsonic side. The same holds on the supersonic side for ω > ω max . To avoid double counting one can either choose all the positive normalization modes (for both in and out), as is done in [11], and then one has to deal with negative frequency modes, or as is usual in this Hawking context, choose only positive frequency modes and then interchangeγ i ↔ −γ † i for the negative normalization ones 5 . In this work, we adopt the convention of ω i > 0 and write The normalization chosen in equations (8)-(10) guarantees that modes are normalized to unit quasi-particle current and so [γ I,ω ,γ † I ,ω ] = δ I I δ(ω − ω ). An identical expression may be written changing in → out. Standard scattering theory arguments show that the S(ω)-matrix coefficients connect the in-modes to the out-modes: Due to the pseudo-Hermitian character of the BdG equations (5), this S(ω)-matrix obeys, for frequencies ω < ω max , a pseudo-unitary relation S † (ω)ηS(ω) = η with η = diag(1, 1, −1). This enforces non-standard relations among the transmission and reflection coefficients, [8,11,14]).
The case ω > ω max is simpler because the S(ω) matrix is unitary and there are no anomalous reflections or transmissions. Thus all asymptotic states can be chosen with ν > 0. If we rewrite R I (ω) := S I I (ω), T I J (ω) := S I J (ω) with I, J different and taking values I, J = u, d1, then the usual unitarity relation |R I (ω)| 2 + |T I J (ω)| 2 = 1 applies. The upstream phonon flux spectrum can be computed from these considerations and for ω < ω max shows the remarkable form [8,14] dI uout dω Assuming that we can populate the ingoing fluxes with comoving thermal populations, namely where k α-in (ω) is the solution to equation (7) for given ω, with v α = v u , v d the flow velocities. Equation (14) reveals that, even at T = 0, a non-zero upstream flow of energy (phononic Hawking radiation) must be expected. The spectrum for ω > ω max is of the same form as equation (14) but with the last term removed, i.e. without any zero-point energy flux. Figure 6 shows some frequency spectra of the upstream phononic flow (see equation (14)) for the setups and condensate solutions depicted previously in figure 5, with a one-to-one correspondence between the graphs in the two figures. The top-left graph of figure 6 shows a structureless profile with a peak at ω = 0 followed by a 1/ω tail, that is typical of Hawking radiation profiles in the absence of barriers, and very similar to what is obtained in the singledelta barrier case. The message is that two nearby barriers behave similarly to a single barrier of double strength. We may also note that the zero-temperature contribution (thick blue line) is not easily distinguishable from the total contribution at non-zero temperature (thin red line), because at low frequencies they both follow the same thermal law 1/ω. This property seems to be common to all structures which do not permit one or more nonlinear oscillations of the condensate between the barriers. An exception to this trend occurs when the delta barrier strength z is close to its upper limit. Then a bigger separation among barriers is possible and, as a consequence, a peak at non-zero ω may develop (see the upper-right graph in figure 6). A general trend that can be clearly appreciated in the rest of the graphs is that the larger the separation the greater the number of peaks that appear. In particular, the two bottom graphs in figure 6 exhibit a double peak in the allowed frequency interval. The reader might wonder why panels (e) and (f) of figure 6 look somewhat different since the parameters hardly vary, the only difference being a relative change in the inter-barrier distance of approximately 0.01. The reason is that the net amplification factor results from a rather complicated expression involving interferences between the scattering coefficients on both sides of the scattering region. See equation (69) and figure 5, right panel, of [9], where a similar sensitive dependence has been found. We note that spontaneous phonons are generated at the event horizons, and after multiple scattering by barriers and horizons, they are partly emitted into the subsonic side. Those scattering events include normal processes (u and v components do not mix) and Andreev (or anomalous) processes (u converts into v or vice versa; see [11]). The bigger the separation between the barriers, the greater the number of quasi-bound states that can be accommodated between them and the greater the number of peaks that appear in the Hawking radiation spectrum.

Hawking radiation spectra
An important point that will be discussed in greater depth in the next section is that peaks can be due to resonances (also called QNMs) or instabilities, so more information is needed to classify them. An unequivocal experimental signature of Hawking radiation (associated with zero-point quantum fluctuations of an otherwise stationary state) would undoubtedly be favored by a transport regime where the classical flow is dynamically stable. We recall that the S(ω) matrix is pseudo-unitary in the Hawking sector 0 < ω < ω max . The general spectral theorem on these matrices [29] guarantees that eigenvalues come in pairs of s i , 1/s * i or have unit modulus, |s j | = 1. A general trend we have observed is that in the region ω ω max /3 most of the scattering matrices show only one unit-modulus eigenvalue. By contrast, in the region ω ω max /2, most of the S-matrices have three unimodular eigenvalues. Nevertheless, the theorem mentioned before guarantees that the determinant will be a pure phase. Then, as in a conventional phase shift analysis of quantum mechanical scattering, a phase jumping upwards when crossing through a peak from low to high ω can be interpreted as a resonance.  In contrast, a jump downwards will reveal an instability. Hence, the most convenient way to detect resonances while distinguishing them from instabilities is to simultaneously plot the phase of the determinant of the S(ω) matrix. This corresponds to the thin green line appearing  (14)), the so-called Hawking radiation. The parameters of the various graphs are identical to their counterparts in figure 5. Thick blue: the current spectrum at k B T = 0. Thin red: finite temperature (in units of µ), 0.1 for the two uppermost curves and 0.3 for the rest. Dotted green: phase of the determinant of the S(ω) matrix. hω max /µ = 0.9299 for all the graphs except (b), which hashω max /µ = 0.9859. We draw the reader's attention to the fact that the peak in panel (d) corresponds to a dynamical instability, as can be seen from the drop of the green line as it approaches the peak for lower values of ω, and not to a resonance, as is the case in the other panels.
in all the graphs of figure 6. All curves show a clear resonance behavior, with the exception of the middle-right phase curve, which reveals an instability, as indicated by the sudden drop in the green line when traversing the instability. Most of the cases that we have studied show 13 QNM behavior, but the occurrence of instabilities is not so rare that it can be ignored. From our analysis, one cannot rule out the existence of strong instabilities with a large imaginary part in the eigenfrequency, since these would generate not easily detectable structures in the frequency spectrum. However, from the systematic analysis of [9,10], where such instabilities were not found, one can conjecture that also in the present case those strong instabilities will not be found. An analysis of the time-dependent GP equation in real time could establish that this is indeed the case. These considerations underline the need for a systematic search of instabilities and QNMs as poles of the propagator in the complex energy plane [28].
The examples shown in figures 5 and 6 were chosen because they clearly reveal general trends. Unfortunately, in a TOF expansion, which correlates the long-time density distribution with the momentum distribution of the initial state, their peaks barely stand out above the background of the depletion cloud, and even less when thermal fluctuations are taken into account. We have included in figure 7 a set of two more favorable setups that show large signals for zero-point Hawking radiation. The two bottom graphs show the momentum distribution of the initial trapped state and are computed in the approximation where only the subsonic flow is included and boundary effects are neglected. When a TOF measurement is performed in such systems, the depletion contribution is negligible at momentum values, which however reveal clear resonant peaks. Zero-point Hawking radiation nearly exhausts those peaks at low temperatures and still gives the main contribution to the area under the peak at temperatures as high as 0.9µ. This fact could allow for the unequivocal detection of Hawking radiation.

Quasi-normal modes (QNMs; resonances) and dynamical instabilities
A discussion of instabilities in BEC BH analogues can be found in [9,24,25,30]. For a study of QNMs in a similar context, see [31]. General studies of QNM in gravitational BHs and in optical analogues can be found in [32] and [33,34], respectively. A more complete study for the setup considered in this paper will be presented [28].
We have said that peaks in the spectrum of Hawking radiation may be due to resonances or instabilities. The former are characterized by poles of the analytical continuation of the S(ω) matrix in the complex ω-plane with Im(ω) < 0, whereas the latter have Im(ω) > 0. Moreover, both types of complex modes may have |Im(ω)| so large that they are not clearly appreciated in the radiation spectrum, yet they can hide important instabilities. While a systematic search for poles is left for a future work, we discuss here some general features of the behavior of QNMs and instabilities in the complex k-plane.
When a linear stability analysis is made of a stationary solution of the GP equation (2), instabilities appear as solutions of the BdG equations (5) whose frequency ω has a positive imaginary part. Because complex frequency implies complex asymptotic wave numbers, these solutions must be localized in real space. For a given complex frequency, there are always four complex ks (counting each different k with its multiplicity), two with Im(k) > 0 and two with Im(k) < 0.
Like in the search for bound states in conventional quantum mechanics, a discrete, possibly empty set of modes is to be expected. The mode with the largest Im(ω) dominates the longtime behavior and its inverse is a good measure of the decay time of the condensate due to that instability. There is a well-defined procedure to quantize such unstable modes [35], which show up as essentially free particles (instead of oscillators). In open systems such as the one analyzed at present, there are however some subtleties in that quantization procedure which, to  (7), for the two asymptotic sides (left is subsonic and right is supersonic) of a BH configuration. Units of ξ −1 u (left) or ξ −1 d (right) are used. The paths shown result from tracking the k solutions as Im(ω) evolves from zero to a positive value, leaving the real part fixed. The labels and colors correspond to those of figure 4, but we have included the evanescent (green) and exploding (brown) solutions. The parameters are chosen so that the structure has a current q = 0.01 and Re(ω) is fixed to a value 0.8 < ω max . our knowledge, have not been properly addressed in the literature. Specifically, some spectral properties of the BdG operators are used that can only be guaranteed for finite-dimensional operators [36].
At a given real frequency between 0 and ω max we may have four complex k solutions on each side of the interface. These are shown by circles in figure 8. On the upstream side (left graph), they correspond to an incoming, outgoing, evanescent and exploding solution. Downstream (right graph) we have two incoming and two outgoing solutions, corresponding to the normal (d1) and anomalous (d2) channels. For instance, a conventional retarded scattering state is characterized by one incoming channel matching all possible outgoing and evanescent solutions, with no amplitude for exploding or other incoming solutions.
An instability is a bound state made exclusively of spatially decaying solutions that happen to match at a particular complex frequency (with positive imaginary part). Localized solutions correspond to Im(k) < 0 on the upstream side and to Im(k) > 0 on the downstream side. Figure 8 shows how the various k solutions evolve as ω varies from Im(ω) = 0 to Im(ω) > 0, whereas Re(ω) remains constant. Eventual instabilities can only be obtained by matching the evolved 'u-out' and 'evanescent' solutions upstream (Im k < 0) and the evolved 'd1-out' and 'd2-out' downstream (Im k > 0), with vanishing amplitudes for the other solutions. We emphasize that the wave function of these unstable modes involves exclusively analytical continuations of outgoing and evanescent solutions.
QNM modes, on the other hand, are obtained from the analytical continuation to the Im(ω) < 0 half-plane of the retarded Green's function of the time-dependent version of the BdG equations (5) (see [33,34] for a discussion in an optical context but precisely translatable to the present one). More specifically, these QNMs can be obtained as an analytical continuation to Im(ω) < 0 of scattering states involving only outgoing and evanescent solutions, i.e. without the intervention of incoming waves. Like for instabilities, a discrete number of QNMs is to be expected. Therefore, in order to find resonances, we are only interested in the evolution of 'u-out' and 'evanescent' solutions on the upstream side, and 'd1-out' and 'd2-out' on the downstream side, ruling out the other solutions. Interestingly, this is exactly the same set of solutions that are relevant in the search for instabilities. We conclude that instabilities and resonances are obtained from the matching of the same set of solutions albeit in different regions of the complex energy plane.
Hence, QNM/instabilities correspond to poles of the S(ω) matrix when analytically continued to the lower/upper complex plane from the real energy line. This can be readily seen if one allows for a general coefficient in the in-wave component of the scattering solution equation (8), for example. So both resonances and instabilities are candidates to explain peaks in the square of a given S(ω) matrix coefficient. The frequency of those possible underlying modes must of course have an imaginary part much smaller than their real part. Since the determinant of the S(ω) matrix is a pure phase (as noted in the previous section), it can be used to discriminate between both behaviors. Specifically, a jump upwards/downwards of the phase when traversing the peak with increasing frequency implies a QNM/instability. Figure 6, together with other not shown data, reveals that peaks in the Hawking spectrum are mostly due to resonances, instabilities being more the exception than the norm. While this sampling is encouraging, a systematic search for QNMs and instabilities in various stationary flows is left for future work.
Finally, we note that for both QNMs and instabilities, the boundary conditions described here are different from those adopted in [30,31].

Discussion and conclusion
We have studied the flow of an atomic condensate through a double-barrier interface separating regions of subsonic and supersonic flow. Such a setup provides a scenario where Hawking radiation into the subsonic side is enhanced at some frequencies due to multiple scattering of quasi-particles by the two barriers and the modulations of the condensate. The resulting highly non-thermal Hawking radiation presents peaks at frequencies that may lie well above the working temperature and thus can be unambiguously interpreted as stemming from quantum fluctuations of the quasi-particle vacuum. The non-thermal Hawking spectrum emitted by the double-barrier interface represents an important advantage over the cases of single or zero barrier, where the low-frequency zero-point radiation has a thermal character which makes it more difficult to distinguish it from genuinely thermal radiation.
Our calculation is based on a model of stationary flow of condensate and quasi-particles through a double-barrier structure with open boundary conditions, where the condensate density is asymptotically flat on both sides of the structure, while quasi-particle motion is described by scattering states characterized by incoming channels that are thermally populated. While a stationary scattering picture of transport has proved to hold predictive power in electron systems, it still represents an idealized scenario in cold atom contexts, where stationary circuits have not yet been developed and where transport of finite-sized condensates is mostly investigated within a time-dependent scheme [37]. As long as steady condensate transport is still an item for the future, it will be of interest to perform numerical simulations of time-dependent transport that may reveal those features of the Andreev-Hawking phenomena that we have explored here.
An important question is whether, in currently achievable setups, the resonance frequency can be tuned to lie well above the currently attainable temperatures that would characterize the incoming quasi-particle population. We have seen that, while in a TOF experiment the contribution from Hawking radiation peaks may be easily overshadowed by the contribution of the depletion cloud, setups can be designed where the resonant Hawking peaks are sufficiently sharp to be clearly visible even at temperatures comparable to the chemical potential.
We note that x 0 0, which means that ρ (0 − ) < 0, where ρ(x) := | 0 (x)| 2 → 1 as x → −∞; see below. The chemical potential and uniform current read where 0 q < 1 is required to obtain j > 0 and subsonic flow. On the supersonic side, the condensate wave function is a simple plane wave, which by current conservation has to be of the form 0 (x) = A min e i ( q x/A 2 min +φ a ) , A min := where φ a is a phase to be determined later. There is a simple analogy that permits us to qualitatively understand the behavior of the condensate wave function. Shifting to an amplitude and phase representation and splitting the GP equation (A.1) into its real and imaginary components in the regions where V ext (x) = 0, we may write The third equation, which comes from the imaginary part of equation (A.1), is but the continuity equation. The second equation, after multiplication by A (x) and one integration, can be written as (cf equation (A.3)) E being an integration constant. Therefore (see e.g. [42,43]), by looking at the position as a time coordinate and the amplitude as the position of a fictitious particle, this equation expresses the energy conservation for this particle under a force deriving from the potential  4) and [19]). Integration of this equation with the initial condition lim x→−∞ A(x) = 1 leads to (A.2). We note that the integration constant must be chosen as E max ≡ V q (A max ) = V q (1). We will also use the quantity E min ≡ V q (A min ). A delta barrier zδ(x − x 0 ) appears as an instantaneous kick, occurring at 'time' x 0 , which adds a positive amount to the speed of the particle, which now becomes A (x + 0 ) = A (x − 0 ) + 2z A(x 0 ). The 'energy' of the particle changes instantaneously without changing the phase φ(x 0 ) or 'position' A(x 0 ). In both the single-and double-barrier cases, the first kick, δ(x), must take place when the particle has negative 'speed', A (0 − ) < 0; otherwise the kick would increase the energy with no chance of ending at A min . The case of a single delta-barrier at the origin, V 1 (x), leads to a simple change in energy, which can be solved to obtain z(q) (or q(z), as the dependence is monotonic; see figure 2): and the parameter x 0 can be extracted from (A.2) by imposing | 0 (x 0 )| = A min . The case of two delta barriers, V 2 (x), is more involved. In this case, the energy after the first kick must lie between the two limiting values E min < E < E max , because if bigger than E max , it would mean that the particle has acquired a positive speed, and the second kick could never send it to A min . To end up in A min with zero kinetic energy, A (d − ) < 0 is required. We conclude that, for a given current q, the maximum strength z that allows for a solution is (A.7), which is the unique z value for a single delta-barrier. From the matching at x = d, it follows that the energy for the particle motion between kicks (0 < x < d) is Analysis of the x = 0 kick leads to another equation that, combined with (A.8) computed at x = 0 + , leads to the following two equations: which allows us to solve for A(0), A (0 + ) as a function of q, z. If we introduce the densities ρ(x) = A 2 (x), ρ 0 := ρ(0), ρ 0 := ρ (0 + ), then equation (A.9) can be rewritten as ρ 0 = 2zρ 0 − q 2 + 1/2 − E(q, z) z , (A.10) whose plot is shown in the right graph of figure A.1 6 . The minimum z that allows for a solution can be computed by enforcing that the two lines cross tangentially (see the blue line in figure 2). From the above discussion, a method can be outlined to compute the possible distances d between the barriers for given q, z between the relevant limits: one takes one of the two solutions of (A.10) and propagates it (clockwise in the blue line of the right graph of figure A.1) until ρ(x) = A 2 min , ρ (x) < 0. The 'time' x needed is the minimum distance between barriers. An integer number of periods can be added, where the period is the time needed to get back to the initial point. Two families of solutions result, corresponding to the two solutions that can be chosen from equation (A.10).

B.2. Supersonic region
In this case, there is no need to use the technique of squaring the Jost functions, as the solution can be trivially worked out. We only quote the results, using the notation of equation (A.4) and introducing p := q/A min , c := A min , which are, respectively, the condensate wave vector and the speed of sound in the supersonic region: where G ± (k, q) is a shorthand notation for G ± (k, q) ≡ k 2 /2 + c 2 k 2 (k 2 + 4c 2 ) ± We are assuming, without loss of generality, that ε > 0. The first solution (B.10) is always valid and is normalized to |u(x, t)| 2 − |v(x, t)| 2 = 1. The second solution is propagating only when ε < ω max and its normalization is anomalous: |u(x, t)| 2 − |v(x, t)| 2 = −1. As expected, the spectrum is that of Bogoliubov quasi-particles for the supersonic region: ε = pk ± k 2 4 2c 2 + k 2 2 . (B.12)