Low frequency analogue Hawking radiation: The Bogoliubov-de Gennes model

We analytically study the low-frequency properties of the analogue Hawking effect in Bose-Einstein condensates. We show that in one-dimensional flows displaying an analogue horizon, the Hawking effect is dominant in the low-frequency regime. This happens despite non vanishing greybody factors, that is, the coupling of the Hawking mode and its partner to the mode propagating with the flow. To show this, we obtained analytical expressions for the scattering coefficients, in general flows and taking into account the full Bogoliubov dispersion relation. We discuss the obtained expressions for the greybody factors. In particular, we show that they can be significantly decreased if the flow obeys a conformal coupling condition. We argue that in the presence of a small but non-zero temperature, reducing greybody factors greatly facilitates the observation of entanglement, that is, establishing that the state of the Hawking mode and its partner is non-separable.


I. INTRODUCTION
In 1981, Unruh discovered an analogy between the propagation of sound in a flowing fluid, and that of waves around a black hole [1]. He immediately realized that this analogy could be used to build an experiment to detect, and study the mechanism of the Hawking radiation [2], which predicts that black holes spontaneously emit a flux of thermal radiation. In the recent years, this analogy has received an increasing interest from various experimental groups, and several signatures of the effect has been observed in various experimental setups [3][4][5][6][7].
In this work, we investigate the properties of the analogue Hawking effect in Bose-Einstein condensates. The differences between the spontaneous emission of phonons in a Bose-Einstein condensate and the prediction of Hawking have two origins. First, when the wavelength decreases, the dispersion relation is no longer linear [8,9], an assumption the gravitational analogy is based upon. Second, modes that propagate along with the flow are not responsible for the Hawking effect, but they affect the observables by coupling to the Hawking mode and its partner. In black hole physics, they give rise to the so-called greybody factors [10], which can significantly effect the emitted spectrum. A particular consequence of these greybody factors is that the low-frequency emission spectrum of a black hole is suppressed, because low-frequency modes do not have enough energy to escape from the gravitational potential. On the contrary, in one-dimensional acoustic flows, it was discovered that the spectrum of the Hawking effect is dominant at low frequencies [9,[11][12][13][14][15].
In this work, we study the low-frequency properties of the Hawking effect, and demonstrate that the spectrum increases at low frequencies in the most general case. We consider both greybody factors and the effect of dispersion (not considered in [13][14][15]). For this, we use a generalization of the matched asymptotic expansion method [16,17], that allows us to fully characterize the scattering at low frequencies. We discuss this scattering in the most general class of flows that can be obtained in condensates. This allows us to characterize greybody factors, but also to point out how to minimize them, and hence recover predictions very close to the ideal Hawking case.
In a companion work [18], we used the same method to characterize the low-frequency effects of dispersion. We worked in a simpler model, the linear Korteweg-de Vries model, which neglects the presence of the downstream mode. In this work on the contrary, we consider the full Bogoliubov-de Gennes equation, which describes the propagation of sound waves in a Bose-Einstein condensate. We shall not only discuss the emitted spectrum, but also how the entanglement between the Hawking mode and its partner is affected by the coupling to the downstream mode.
In section II, we briefly review the Bogoliubov approximation of the excitations of a condensate, and the definition of the scattering matrix on a varying flow. In section III, we present the method and obtain the general behavior of the coefficients of the scattering matrix. In section IV, we discuss the various consequences of greybody factors, and how to control them. The conclusion are complemented by the study of two exactly solvable examples. We work in units where = k B = 1.

A. Wave equation
In this section, we briefly present the problem of scattering of density waves in transonic flows. We refer to the literature for a more detailed account and the precise relation to black hole physics [9,[19][20][21]. We consider a gas made of N identical spinless bosons of mass m with a point-like interaction. In the one-dimensional regime [22,23] (when transverse trapping frequencies are higher than typical energy scales along the preferred direction) 1 , the bosonic field operatorΨ evolves according to V ext is the external potential, and g the (effective one-dimensional) coupling constant of inter-particle interaction. To consider the most general case, we allow both to depend on the position. However, we restrict ourselves to stationary backgrounds; hence, V ext and g are independent of time t. We assume the temperature to be low enough for the system to be in the quasi-condensate regime [23]. In this regime, the condensate wave function Ψ = ψ 0 e −iµt obeys the stationary Gross-Pitaevski equation, i.e.
where µ is the chemical potential. The condensate wave function defines the density ρ and velocity v of the flow Perturbations around this condensate are described by a phononic field operatorφ defined byΨ = ψ 0 e −iµt 1 +φ .
The field operatorφ can be decomposed into a superposition of stationary modeŝ where the sum runs over the different modes of the same frequency ω. Linearizing equation (1) inφ, we obtain the field equations for the modes φ ω and ϕ ω as To relate this equation to physical observables, we will use phase fluctuations θ ω and relative density fluctuations n ω = δρ/(2ρ) instead of the field modes φ ω , and ϕ ω . It turns out that the mode equation (6) is also easier to solve using θ ω and n ω . The two are directly related of field modes by They obey the set of equations We start by considering a homogeneous condensate flowing to the left (v < 0). Since the background quantities are independent of x, solutions of (8) are given by plane waves The frequency ω and the wave number k are related by the dispersion relation In this equation, we have defined the speed of sound c 2 = gρ/m, which gives the propagation speed of long wavelength waves. Using this velocity c, one can build a characteristic length ξ = 1 2mc . (11) ξ is the healing length. In equation (10), it determines the length below which dispersive effects become significant. In the next section, we will solve explicitly this equation for low frequencies. We first classify the various roots of equation (10), using its graphical resolution in Fig. 1. When the flow is subsonic |v| < c, there are always two solutions, one moving upstream (noted k u ), one moving downstream (noted k d ). When the flow is supersonic, and below a certain frequency 0 < ω < ω max , there are 4 solutions. Two of them correspond to the usual right and left movers, but the one going against the flow propagates too slowly, and hence is dragged by the flow. The peculiarity is that this root corresponds to a mode of negative energy (or norm, see below). In addition, we have two extra solutions. These solutions are allowed by the non-zero flow and by dispersion, i.e. they disappear in the limit ξ → 0. For this reason, we shall refer to them as the dispersive roots. One of them has a positive energy, while the other have a negative energy. We will denote them k + and k − , where the subscript refers to the sign of the energy. In addition to the dispersion relation, Figure 1: ω as a function of k from the dispersion relation (10), in units where c = 1 and ξ = 1. On the left side, in a supersonic flow (v = −1.2), on the right side, in a subsonic one (v = −0.7). The dashed line shows the number of solutions for ω = 0.1. The value of ω such that k − = k −u defines ω max . To compare with the work of [18], the dotted-dash line represents the Korteweg-de Vries approximation, that is, ω = (c + v)k + k 3 /(8m 2 c).
the mode equation (8) gives us a linear relation between the amplitudes U k and V k defined in equation (9). Using (8b), we obtain This means that the phase operatorθ = (φ −φ † )/2i decomposes into eigen-modes in the following manner:θ where the index j varies over the solution of the dispersion relation at fixed ω. The density fluctuations operatorn = (φ +φ † )/2 has a the same decomposition, with the amplitudes V k instead of U k . The canonical commutation relation [φ,φ † ] = ρδ(x − x ) now gives us the commutation relations between the operatorsb j 's. Then, to identify them as creation and annihilation operators, satisfying the canonical commutation relation we normalize the amplitudes such that (see e.g. [9] for more details). The sign in equation (15), is the sign of the norm of the corresponding mode. It indicates if that mode is associated with a creation or annihilation operator [9,25]. In equation (15), we have also defined v g , the group velocity associated with the mode of wave number k, i.e. v g = (∂ k ω) −1 . For simplicity, we can choose the phase reference such that U k is real. Then the normalization is given by . (16b)

B. Varying background and S-matrix
We now assume that the flow accelerates, or decelerates, over a region of length scale L, and centered around x = 0 for commodity. When x −L or x L, the flow is constant, that is where the subscript l (resp. r) is for x → −∞ (resp. +∞). Notice also that the condensate flows to the left, hence v r/l > 0. The 5 background functions in (17a) are not independent, since the Gross-Pitaevski equation (2) must be satisfied. The first relation is the continuity equation, which, for stationary unidimensional flows, gives ρv = const.
The second equation obtained from the Gross-Pitaevski equation (2) gives a relation between V ext (x) and the density ρ. In addition, the definition of the speed of sound gives us c 2 (x) = g(x)ρ(x)/m. Using these three relations, the background is fully characterized by two independent background functions out of equation (17a). In other words, by choosing the external potential V ext (x) and the coupling constant g(x) in (2), one can impose a certain profile for ρ, v, and c. In an experiment, these functions might be delicate to control with the necessary precision. Possibilities to have a varying coupling constant g are to adiabatically change the transverse trapping frequency [26], or to exploit a Feschbach resonance [23,27]. In this work, instead of V ext and g, we shall use v(x) and c(x) to fully determine the propagation equation of waves in the condensate through equation (8), and assume generic form for their profiles. On both asymptotic sides |x| L, the solutions of (8) are given by superpositions of plane waves. From the earlier discussion of the dispersion relation (10), there are 2 asymptotic modes on the subsonic side and 4 on the supersonic one (since we are interested in low frequencies, we always assume ω < ω max ). As read from their group velocity, 3 of them propagate towards the transition region, while 3 of them propagate away from it. The linear relation between modes going in and modes coming out defines the S-matrix: This is illustrated in Fig. 2. Because we work with normalized modes, the conservation of the norm (or energy) implies that S ∈ U (2, 1). However, as we shall see, at low frequencies, the dispersive modes θ ± tend to have the same amplitude, due to the fact that |α| ∼ |β| and |α| ∼ |β| (see equations (48) and (49) below). Therefore it will be difficult to build the scattering modes θ in ± . To circumvent the problem, we shall construct out modes instead. We then obtain the scattering coefficients by inverting the S-matrix. Since S ∈ U (2, 1), this inversion is fairly simple, and we have The minus signs follows from the fact that S ∈ U (2, 1). Moreover, this property of S also gives several relations between the coefficients, such as Both in and out basis allow for a decomposition of the field operators, i.e.
(Since we focus in this work on low frequencies, we have omitted frequencies ω > ω max . In that case the scattering becomes elastic and 2 × 2 [9].) The two possible decompositions of this operator leads to a linear relation between theâ out 's andâ in 's, inherited from the S-matrix, which reads When incoming modes are in their ground states, because the scattering mixes positive and negative norm modes, there is a spontaneous emission of phonons. This is directly obtained from equation (23). For instance, the flux of u-modes is given by And similarly for the emission of the other modes n −u , n d . When neglecting dispersion and the mode θ d , the emission spectrum of out-going modes follows a Planck distribution: where the Hawking temperature is given by the gradient of the flow at the point H where |v| = c (equivalently, by the surface gravity of the analogue horizon), that is It has been shown in several works [28][29][30][31][32][33][34][35] that (25) is maintained when L ξ, and when θ d is neglected. In this work, we relax these two assumptions, and study in depth the lowfrequency behavior of the coefficients. In the Hawking regime, the low-frequency behavior is simple, as it follows from (25) that This result is in general affected by the coupling to the downstream mode θ d . The dressing of the flux n u by the downstream mode is referred to as greybody factors in the black hole literature [10]. In the rest of the paper, we shall focus particularly on the coefficients R, and B, which are the generalization of the greybody factors 2 for an acoustic black hole in a Bose-Einstein condensate. As we shall see in section IV, these coefficients not only dress the emitted spectrum, but they also alter the amount of entanglement between the Hawking mode θ u and its partner θ −u . Lastly, we point out that our conclusions will equally apply to white hole flows (time-reverse of black holes), since the corresponding S-matrix is simply S −1 (see [33] or the appendix D of [9] for more details). Figure 2: Schematic representation of the scattering problem. We have represented the amplitudes the asymptotic behavior of θ in − on a black hole flow (see equation (19)).

A. Method and mode basis
In this paper, we will analyze the properties of the S-matrix in the limit of small frequencies. For this, we divide the problem into two regions. In the far region, for |x| L, the background is homogeneous, and solutions are a superposition of plane waves. In the near region, for |x| c/ω (where c is here the typical value of c(x) along the flow), ω can be neglected in the mode equation (8), which becomes easier to solve. Then, we identify the large x behavior of the solutions obtained for ω → 0 to the low-frequency expansion of the mode basis in a matching region such that L |x| c/ω. Therefore, this procedure is valid under the condition that this matching region is large enough. This gives us the validity regime of the low-frequency scattering coefficients we will obtain, namely This method is is known as matched asymptotic expansion. It is wildly used for second-order differential equations, in quantum mechanics [16], or black hole physics [10,38], and has also been used to analyse greybody factors in acoustic metrics [13,14]. What we present here is a generalization of this method to equation (8) (which can be recast as a fourth-order differential equation [9]), in order to take into account dispersive effects. The first step is to build the mode basis at low frequencies. To enlighten the notations along the paper, we first introduce the characteristic wave numbers on the left and right sides as We now solve the dispersion relation (10) for low ω, on both asymptotic sides. For x L, they read while for x −L, we have Note that we have kept the subleading term in O(ω) for k ± in order to obtain their group velocity. We can now build the mode basis on each side by using the above roots of the dispersion relation on plane wave solutions given by (9), and (16). We then evaluate them for |x| c/ω. On the right side we obtain and on the left we have To simplify, we called U Z the common normalization of the dispersive modes, which reads In the next section, we will solve the wave equation at low frequencies, and then obtain their asymptotic behavior for |x| L. As we shall see, the large x asymptotics of these solutions are a superposition of constant and linear terms in x, and oscillating exponentials. The existence of a matching region will allow us to identify these asymptotics to the different plane wave modes obtained in equations (32) and (33). Notice that it is necessary to keep the terms linear in x in the low-frequency expansion of θ u , θ −u and θ d , as they allow us to distinguish the two different long wavelength modes.

B. General results
We start by setting ω = 0 in equation (8). Then, using ρv = const, we integrate (8b) into By changing the constant of integration a, one obtains different possible modes. We start by setting a = 0, and will later on consider modes for a = 0. We then plug the above relation in (8a), and get We define an auxiliary field χ as The field χ then obeys the second-order equation This is already a notable result. Namely, the low-frequency behavior of phonon scattering is entirely determined by a Schrödinger-like equation in a potential, for which many tools and solvable examples are know. To start, we will show how the general asymptotic properties of this equation gives the frequency dependence of the scattering coefficients. Since we consider a transonic flow, we have c l < v l , and v r < c r . Using the definitions of (29), the solutions of (38) are asymptotically exponentials with rate q l/r . On the left side, the exponentials oscillate (c 2 l < v 2 l ), while on the right side (c 2 r > v 2 r ) they grow or decay. Explicitly, the asymptotic expansion of χ has the form We underline that the coefficients A's are independent of ω, since they are obtained from equation (38), which do not contain ω. Since all scattering states are spatially bounded, we first set A ↑ = 0. Then, we integrate χ and obtain the asymptotic behavior of θ: Here we have fixed the integration constant to vanish on the left side. This implies that A 1 is given by the integral 3 θ 1 is one particular solution of the mode equation (8) at ω = 0. The first mode we want to obtain is the long wavelength mode coming out to the right, that is θ out u . Using the inverse S-matrix in (20), and the form of the mode basis for L x 1/ω, we have We see that θ 1 is a good candidate, since it is purely oscillating on the left side (no constant or linear term in x on the left side). As previously mentioned, since we have two different long wavelength modes (θ u and θ d ), we must keep the term in O(x) in the limit ω → 0 to be able to distinguish them. As we see from equation (42), this term is higher order in ω.
Hence, one must take into account the first-order correction in ω from the wave equation, and obtain a mode with the asymptotic behavior where A 4 is independent of ω. To obtain A 4 , we first invert the relation between θ ω and n ω in (8b) at first order in ω. Using the fact that ρv is constant, we obtain The integration constant is chosen so that there is no constant or linear term in x near −∞.
We now inject this in equation (8a), and obtain We now extract A 4 by evaluating this for x → +∞. Using the asymptotic (43) and the relation (44), we show that ∂ x n ω is O(ω 2 ). Hence, identifying all terms of order O(ω), we obtain To rewrite this equation in a simpler form, we define an "effective velocity" as It is easy to see that if v(x) is constant, we simply have v eff = |v| (v eff is also independent of a choice of normalization and phase reference of θ ω ). Hence, v eff can be interpreted as some averaging of the background flow v(x) by the mode. Its interest is that it allows for a compact writing of the grebody factors. The last step to obtain the scattering coefficients is to identify N θ 1 to θ in u , where N is an overall normalizing constant. This gives us 4 equations (2 oscillating components on the left and a O(1) and O(x) on the right), and 4 unknown coefficients: α, β, R, and the normalization constant N . Solving this linear system, we obtain the scattering coefficients: This is our first main result. It gives the low-frequency behavior of the coefficient. We see that at low ω, the greybody factor R goes to a constant, i.e. R = O(ω 0 ), while α and β grow like O(1/ √ ω), in agreement with [9]. This means that conversion from short to long wavelength mode, and in particular spontaneous emission (governed by β), becomes very large at small frequencies, while conversion between long wavelength modes is bounded. We underline that this holds for arbitrary backgrounds. The only necessary point is the asymptotic behavior: the flow must make a transition from subsonic to supersonic.
The second mode we will build is the negative energy long wavelength mode coming out to the left, i.e. θ out −u . In the region L |x| c/ω, this mode reads To identify this mode to a low-frequency solution of equation (8), we need to build a solution that is linearly independent from θ 1 . Then we will identify (49) to the linear combination To do so, we first notice that for ω = 0, the couple n = 0 and θ =const. is a solution of equation (8). Since θ 1 has no constant term on the left, we see that in order to have the same asymptotic as in (49a), we can choose θ 2 = mc l /(2ωρ l ), which fixes N 2 to 1. This is, however, not enough. Indeed, as previously, we must obtain the first-order correction in ω. To do so, we write θ 2 under the form Now, by a calculation very similar to the one to obtain equation (45), we show that where a is a constant of integration. We now fix the value of a by taking x → −∞, and identifying ∂ x to the O(x) coefficient in equation (49a). Taking the limit x → +∞ with that value of a gives us the asymptotic behavior of θ 2 . This leads to with What is left to do is now to identify the linear combination (50) to the low-frequency mode (49). This gives us again a 4 × 4 linear system, which we invert to obtain To obtain the last mode, we need to identify a different linear combination (50) to We can proceed in exactly the same way as above. However, it is quicker to notice that θ out d can be obtain from θ out −u by the formal replacement c l → −c l . Either way, we obtain the rest of the scattering coefficients We see here that the coefficientsR andB, although they also describe coupling to the mode θ d , behaves quite differently from R and B, and scale like O(1/ √ ω). The reason is that R and B describes the coupling between the long wavelength modes, whileR andB encode the production of the long wavelength mode θ d when sending a dispersive mode θ ± . In addition, although the expressions ofR andB are quite similar toα andβ, they differ in an important respect. When θ d decouples, R and B becomes small, whileT becomes close to 1. In fact, the combination 1 −T v r c r /v l c l becomes small. This means that when the coupling to θ d is small,R andB are significantly smaller thatα andβ. This is possible becauseT also scales as O(ω 0 ), as R and B, since it also encodes a transition from a long wavelength mode to itself. As a last remark, we point out that the three coefficients R, B andT are related by norm conservation. Indeed, normalization of the third line of S imposes which can be checked explicitly to hold.

A. Behavior of greybody factors
There are two possibilities to minimize the effects of the greybody factors R and B. The first one is to approach the near-critical regime. When the flow speed is close to the sound speed on both sides, i.e. v r ∼ c r ∼ v l ∼ c l , the mode θ d universally decouples. More precisely, all greybody factors decrease as O(c − v). The second possibility is to tune the variations of v and c (which can be achieved by controlling V ext (x) and g(x) as mentioned in section II B) so that R and B vanish.
We saw in the preceding section that using the effective velocity v eff , the greybody factors have a simple and universal form. In particular R and B only depend on dispersion and the details of the profile such as its length L, through v eff . It is therefore useful to start by characterizing this effective velocity. The first point to notice is that in the smooth limit L ξ, v eff becomes equal to the absolute value of the velocity at the This is because, in the limit L ξ, the mode θ ω is accurately given by WKB waves, except near the horizon where v 2 (x) = c 2 (x), where there is a turning point. Hence the integral defining v eff is mainly governed by the vicinity of that point (see appendix B for more details). In the step limit L ξ, v eff can also be computed explicitly (this is done in appendix B, see equation (B7)), and we see that it is between v r and v l , depending on how symmetric the flow is.
The first important result from the expressions (48a) for R, and (55a) for B, is that the spectator mode θ d always decouple in the near-critical regime. Indeed, we see that if v r ∼ c r ∼ v l ∼ c l , then v eff ∼ v r , and hence R and B vanish. In fact, in this regime, the scattering is accurately described by the Korteweg-de Vries model, that we studied in [18]. This will be further confirmed in section IV C, where we show (in explicit examples) that the expressions for the other coefficients reduce to the one of the Korteweg-de Vries model. It is also possible to reduce significantly the mixing with θ d by carefully choosing the relative variations of v and c, even outside the near-critical regime. Indeed, we see that R = 0 (resp. B = 0) if we have v 2 eff = v r c r (resp. v 2 eff = v l c l ). To see this in more details, we first consider the smooth regime L ξ. We introduce a new set of (mutually independent) parameters M l/r is the Mach number on the left or right side, while p l/r encodes how symmetrically v and c varies off their horizon value v 0 . When p = 1, inhomogeneities are symmetrically shared between v and c. This parametrization is quite convenient, since R and B only depend on p r/l , and not on the Mach numbers themselves. Using equations (48a), and (55a), we obtain We first see from these equations that when p r = p l = 1, the two coefficients vanish, which means that θ d completely decouples (see Fig. 3). Since our treatment only provides the coefficients up to higher powers of ω, when these expressions vanish, we can only conclude that R = O(ωL/c) (or B). When relaxing the assumption L ξ, the coefficients R and B are still given by equation (61), but with p r/l defined with v eff instead of v 0 in (60). Hence, one can still obtain a vanishing R (resp. B) if v eff = v r c r (resp. v l c l ).
In the smooth limit L ξ, the fact that θ d decouples for p r = p l = 1 comes from a hidden symmetry of the phononic mode equation (8). It is known from early works in analogue gravity [1,19] that the phononic wave equation of a three-dimensional fluid reduces to the wave equation in a Lorentzian geometry. This is, however, not the case in a one-dimensional fluid [39]. A key difference is that the wave equation in a 1+1 dimensional space-time is conformally invariant while the phononic wave equation is not. However, it turns out that the two coincide under the condition that v(x)c(x) = const. (62) This implies that the phononic equation becomes conformally invariant when this condition is met. In particular, θ d exactly decouples, and hence the greybody factors R and B become trivial for all ω. This property of the phononic wave equation is shown in details in appendix C.

B. Effect of temperature and consequences on the observation of entanglement
In realistic conditions, the condensate has a finite temperature. The incoming dispersive modes have typically a short wavelength, shorter that the thermal wavelength, and hence it is reasonable to assume that they are in their ground state (but see e.g. [40] for more general states). On the other hand, the incoming long wavelength mode is in general thermally excited. Hence we assume that the incoming flux n in d of this mode follows a Planck distribution at the temperature T d : In principle, T d is given by the ambient temperature T ext multiplied by a Doppler factor due to the flow, T d = T ext (1 + |v|/c). In this state (vacuum for θ in ± and thermal for θ in d ), the emitted fluxes of phonons are given by Moreover, an initial temperature has also a significant effect on the correlation properties of the various modes. In the Hawking effect, the emitted mode θ u and its partner θ −u are quantum mechanically entangled 4 . A convenient criterion to show this entanglement is to compare autocorrelations with cross-correlations [40,42,43]. Correlations between classical excitations must obey the following bound, known as the Cauchy-Schwarz inequality: whereρ is the density matrix of the phononic state. In [6], this criterion was used to assess the entanglement in an experimental realization of the Hawking effet in an Bose-Einstein condensate. At zero temperature T ext = 0, this inequality is violated, showing that the pairs are entangled. This changes when one takes into account a temperature for the incoming mode θ d [9,40]. In that case, using equation (23) the difference between auto-correlations and cross-correlations is given by 5 To analyze the influence of the downstream mode on the Cauchy-Schwarz inequality (65), we introduce a parameter ζ uu defined as We see from (66) that depending on the sign of ζ uu , the coupling with the mode θ d increases (ζ uu > 0), or decreases (ζ uu < 0) the possibility of observing entanglement. Using the general expressions we obtained in equations (48) and (49), we see that this parameter is essentially governed by the greybody factors, since it can be expressed as Using the notations of equation (60), this simplifies into This expression has several consequences. First, ζ uu is manifestly positive. This means that the coupling to the mode θ d always reduces the possibility of observing entanglement through the Cauchy-Schwarz inequality. Second, ζ uu is independent of v eff , and in particular of the dispersive scale ξ. As a result, ζ uu does vanish if R and B do, e.g. in a smooth flow with the conformal condition (62), but it also vanishes under the simpler condition that In other words, one can minimize, and even cancel the effect of the downstream mode θ d on the Cauchy-Schwarz inequality even if the greybody factors are non-zero, if their contributions compensate in equation (68). This can be achieved by enforcing the asymptotic values of v and c to satisfy (70). Since this cancelling of ζ uu is shown for small frequencies, the condition (70) implies in general ζ uu = O(ωL/c). Moreover, under the conformal condition (62), more constraining than (70), higher powers of ω are suppressed by the dispersive scale due to conformal invariance (see appendix C), that is by a factor O(ξ/L). This means that under this latter condition, ζ uu is highly suppressed, scaling as O(ωξ/c). The variations of ζ uu with p r/l are shown in Fig. 3.

Constant flow velocity
As a first example, we consider the case where the density ρ = ρ 0 and velocity v = −v 0 are constant, and all the variation is contained in c(x) (and hence in g(x)). Although this is presumably difficult to set experimentally, this case is interesting for its mathematical simplicity. For this reason, it has been previously studied in details in the phononic regime [13][14][15] and in the step regime [11,44]. Our method allows us to generalize these results for any value of the transition size L compared to the healing length ξ. As mentioned above, since v is constant, it immediately follows that v eff = v 0 . Hence, the greybody factors are fixed irrespectively of the profile c(x). Explicitly, we find As we saw in the previous section, we observe here one of the disadvantage of having all inhomogeneities in c. Indeed, R, B, and ζ uu vanish only in the near critical limit. As discussed above, fine-tuning R, B, or ζ uu to be 0 requires to vary both v and c in a specific manner. We now turn to the computation of the other scattering coefficients. When v is constant, the equation (38) for the auxiliary field χ reduces to We now consider the profile With it, equation (38) is exactly solvable in terms of hypergeometric functions. As we saw in section III B, to obtain the scattering coefficients we need the decaying mode solution of equation (38). It is given by Using known identities of hypergeometric function (see appendix A for details), we obtain the asymptotic coefficients A 1,2,3 defined in equations (39), and (41). We shall now focus on the explicit expression of β. As we saw in equation (68),β (and hence α,α) are then readily obtained using the expressions for R and B, which we obtained above. Using the results of appendix A and equation (48c), it follows To compare this with the Hawking result of equations (25) and (26), we evaluate the surface gravity of our background. We first determine using (73) where c 2 = v 2 0 . We then evaluate the gradient at this location to deduce .
We use this to rearrange the expression for β above. We finally obtain This formula is the product of 3 terms with a clear physical interpretation. The first term is the Hawking formula at low frequencies (see (27)). The second term encodes how dispersion alter the Hawking temperature. In the limit L ξ, we have tanh(πq l L) ∼ 1, and the result become independent of the dispersive scale ξ. The last term is the modification of the emitted flux due to the coupling to the mode θ d . It is equal to 1−|R| 2 , and corresponds to the fraction of long wavelength modes that are transmitted from the horizon to infinity. In the smooth regime L ξ, equation (77) is the product of the low-frequency Planck spectrum with the (relativistic) transmission coefficient from the horizon to infinity [13,14]. This is the general result for a black hole in general relativity [10]. To summarize this behavior, it is tempting to define an effective temperature T eff , such that In the smooth limit, T eff ∼ T H , and one recovers the factorization of the relativistic case. However, since the coefficient β always grows like 1/ √ ω at low frequencies, this equation can be seen as a definition of an effective temperature T eff that takes into account greybody factors. In the near critical limit, T eff has the value found in the Korteweg-de Vries model, as expected.

Varying flow velocity and speed of sound
To discuss the generality of the conclusions just drawn, it is interesting to look for another exactly solvable example where both v and c vary. This is slightly harder due to the presence of the operator ∂ x ρ∂ x in equation (38). However, it is possible to construct a profile with similar solution as for (73) by tuning how v and c vary together, but still having their asymptotic values unspecified. For this, we first introduce a new spatial coordinate y is normalized using v 0 , the value of the velocity at the horizon, so it still has a dimension of length (this is a convenient choice, since gradients with respect to x or y are then equal on the horizon). Since |v| = 0 all along the flow, this change of variable is perfectly regular. Moreover, on both asymptotic sides, it simply amounts to a linear rescaling of x. Using this new coordinate, the equation (38) for the field χ reads We see that what matters to obtain low-frequency solutions is the evolution of the Mach number with the new coordinate y. We assume the following profile This profile is quite similar to (73), but the asymptotic values of v and c are now independent. Notice that v 0 the value of the velocity at the horizon is also an independent parameter, since the profile of v itself is not specified. To fix the ideas, we shall assume a similar shape as (81), given by The parameter ∆ gives a shift between the center of the transition of the Mach number M and that of v. Changing the value of ∆ allows us to change the value of v 0 . To illustrate what profiles are described by equations (81) and (82), we represented v and c in terms of the original spatial coordinate x in Fig. 4. Similarly to the flow (73), the solutions of equation (80) are given by hypergeometric functions. In particular, the decaying mode simply reads where we introduced the rescaled momenta Using these rescaled momenta, the expressions of A 1..4 are identical to the one in appendix A. It follows a similar expression for β than (75). To again write it in a convenient way to compare with the Hawking result, we first evaluate the surface gravity of the flow, that is . To compare with other works, we also plotted a waterfall solution (dashed curves), with M l = 5 [21,43].

This leads to
This formula has the same structure as (77), with 3 factors with the same interpretation. As previously, in the smooth limit L ξ, |β| 2 reduces to the product of a Planck spectrum at temperature T H with the relativistic transmission coefficient from the horizon to infinity.
As a last remark, we point out that the second factor in (86), which governs the dispersive corrections to the Hawking temperature, becomes 1 when πQ l L 1. In particular, the dispersive corrections are sensitive to the properties of the flow on the supersonic side. This is physically reasonable, since the dispersive modes live on that side [33]. In addition, Q l L = M 2 l − 1 × L/ξ 0 , and therefore, if the supersonic side has a large Mach number, the temperature will be close to the one predicted by Hawking, even if the condition ξ/L 1 is not well satisfied. For instance, in the experiment of [6], the length of the transition L is essentially ξ 0 , but because M l is larger than about 3, the temperature is expected to be quite close to the one predicted by Hawking, as was noticed in [43].

V. CONCLUSION
In this work we studied the analogue Hawking effect in Bose-Einstein condensates. As we show, the production of pairs of phonons via the Hawking effect is dominant in the lowfrequency regime. We developed a method to obtain analytical results for low frequencies in general flows. In a companion work [18], we analyzed the dispersive corrections to the effective temperature of emitted phonons. In this work instead, we focused on the influence of the mode propagating with the flow. This mode is not the one responsible for the Hawking effect, and is for this reason sometimes called "spectator," but it affects the observables through its coupling to the Hawking mode and its partner. In black hole physics, this mode gives rise to the so-called greybody factors. We studied here the generalization of these greybody factors in the Bogoliubov-de Gennes model.
Our first result is the general low-frequency dependence of all scattering coefficients. In particular, the spectrum of emitted phonons increases like the inverse power of the frequency (see equation (48c)), irrespectively of the coupling to the downstream mode. This is due to the fact that greybody factors tends to a constant in this regime, contrary to what happens in higher than 1 + 1 dimensional black holes [10]. This conclusion was reached in [13][14][15] using the wave equation in an acoustic metric, and is generalized here by taking into account dispersive effects. When the size of the transition from a subsonic to a supersonic flow is large compared to the healing length, the spectrum is given by the curved space-time prediction (see equations (77), (86)). More precisely, it is given by the product of the thermal spectrum with the greybody factors obtained from the relativistic equation.
Our method allows us to quantify the coupling to the downstream mode in general flows. In particular, we identified two distinct regimes where this coupling is small. First, it is always small for near-critical flows, that is, when v ∼ c on both side of the transition. In this regime, the Hawking effect is well described by a simpler model, namely the linear Korteweg-de Vries equation [18]. Second, when the product of the velocity flow with the speed of waves is constant (equation (62)), this coupling vanishes for the dispersionless equation. This means that in smooth flows, it will remain small, scaling like the ratio between the healing length and the size of the transition. This property is due to the conformal invariance of the dispersionless equation when the condition (62) is met (this is detailed in appendix C). We also show that the parameter ξ uu controlling the influence of the downstream mode on entanglement vanishes under a weaker condition (see equation (70)).
In experimental setups, if one can control with enough precision the external potential V ext and the effective 1D coupling constant g, one can reduce the coupling to the downstream mode by ensuring that the flow velocity v and the speed of sound c are related in the appropriate way. Reducing this coupling has several advantages. First, the flux of emitted phonons becomes very close to the Planck law predicted by Hawking. Second, and perhaps more importantly, it ensures that the Hawking mode and its partner are entangled for a large range of frequencies. Indeed, the downstream mode tends to reduce their entanglement via thermal effects or other sources of noise. Ensuring its decoupling significantly eases the detection of entanglement and the observation of violations of classical inequalities. superposition of oscillatory waves, while it decay exponentially on the other side, similarly to equation (39). In the WKB limit, the amplitudes on each side can be obtained using a connection formula [47]. Fortunately, this is not necessary here. Indeed, for a simple turning point (which is the case considered here) the singularity of the WKB wave is integrable, since |χ| = O((x − x H ) −1/2 ). This is enough to evaluate both A 1 and A 4 , both given by integrals of χ (or ∂ x θ ω , see (41) and (46)). We must then evaluate an integral with a slowly varying amplitude and a rapidly varying phase. In this regime, the integral is dominated by its boundary, i.e. the turning point x H . Therefore Using the definition (47) of v eff , we deduce that in the smooth limit, v eff ∼ v 0 , i.e. equation (59). In the other limit, L ξ, we can approximate the background with step function, e.g. v(x) = −v l Θ(−x) − v r Θ(x). The equation (38) for χ can be directly solved. We see that The amplitudes can be determined by imposing that both χ and ρ∂ x χ are continuous across x = 0. It is, however, enough to use the continuity of ρ∂ x χ to obtain