Momentum correlations as signature of sonic Hawking radiation in Bose-Einstein condensates

We study the two-body momentum correlation signal in a quasi one dimensional Bose-Einstein condensate in the presence of a sonic horizon. We identify the relevant correlation lines in momentum space and compute the intensity of the corresponding signal. We consider a set of different experimental procedures and identify the specific issues of each measuring process. We show that some inter-channel correlations, in particular the Hawking quantum-partner one, are particularly well adapted for witnessing quantum non-separability, being resilient to the effects of temperature and/or quantum quenches.


Introduction
While the possibility of quantizing gravitation remains elusive, some noticeable progresses have been made in the description of the interaction between the space-time metric and a quantum field. In particular, the dynamical Casimir effect [1,2] and Hawking radiation from black holes [3] both correspond to a quantum creation of entangled pairs of particles induced by (strong) space-time inhomogeneities, and have both been predicted in the framework of quantum field theory in curved spacetime. In the case of Hawking radiation, the prospect of an experimental study in the genuine astrophysical context seems hopeless, because the radiation has a thermal spectrum at a temperature T H = κ/2πc , where c is the speed of light and κ the black hole horizon's surface gravity (we use units such that k B = 1) and in the standard situation of a black hole formed after a gravitational collapse, T H is much lower than the temperature of the microwave background radiation [4]. However, the phenomenon of Hawking radiation has a robust kinematic origin, and elaborating on the close analogy of a transonic flow structure with the gravitational metric near a black hole event horizon, Unruh proposed to observe Hawking radiation in a condensed matter context [5]: this idea is often considered as the birth of the field of analogue gravity [6,7].
Many physical realizations have been proposed for observing analog Hawking radiation (see, e.g., [7][8][9][10][11][12]), among which the implementation of a sonic horizon in the flow of a Bose-Einstein condensate (BEC) rapidly appeared as quite promising [13]: the low temperature of the system and its paradigmatic quantum nature makes it an ideal playground for studying this phenomenon. However, a direct observation of the analogous sonic radiation in this system is still hindered by thermal effects and difficult to identify unambiguously. The recognition of this difficulty motivated the authors of Refs. [14,15] to propose the detection of density correlation as an evidence for Hawking emission of correlated pairs of particles from the horizon: Indeed, in an analogous system, contrarily to the gravitational case, the experimentalist is a super-observer who can make measurements from both sides of the horizon. The correlation between the Hawking particle and its partner were shown in Refs. [14,15] to induce a characteristic peak in the correlator of density fluctuations which could be used to demonstrate the existence of analogous Hawking radiation in a BEC system. The physical interpretation of this correlated signal is similar to the one initially given by Hawking [3,16]: Just at the event horizon, vacuum fluctuations produce pairs of virtual quasi-particles, one with negative norm and one with positive norm. The negative norm quasi-particle propagates in the region inside the black hole where it can exist as a real quasi-particle (and is often denoted as the "partner"). The other quasi-particle of the pair is denoted as the "Hawking quantum"; it can escape to infinity, where it constitutes a part of the Hawking radiation. In a BEC the quasi-particles are Bogoliubov excitations which correspond to density fluctuations: hence the emission of the correlated pairs of particles induces density correlations. An interesting aspect of these correlations is that they are resilient to finite temperature effects [17,18].
In the same line of idea, we proposed in Ref. [19] to study correlations in momentum space as evidence of Hawking radiation. The physical idea is the same as the one behind the study of density correlations in real space, but the specifics are different, with a number of advantages: first, the practical implementation of this type of experiment is well documented [20][21][22][23]. Also, it was shown in Ref. [19] that the momentum correlation signal is particularly well adapted to the study of Hawking radiation, being even less affected by the background temperature than the real-space correlation signal, and offers a clear and robust signature of the entangled nature of the Hawking pairs. In the present paper we develop and explicit the results presented in Ref. [19]. We detail the theoretical description of the quantum fluctuation of the system and precise how a local Fourier analysis can be performed. This leads us to underline some characteristics of the experimental detection scheme which are crucial for the detection of entanglement (cf. the discussions in Appendix A and at the end of Sec. 3.1). We also extend the treatment of Ref. [19] in order to include what we denote as "non adiabatic effects" and "in situ measurements" (Sec. 3.3). We show that the results presented in Ref. [19] are robust with respect to this more general treatment, and that new correlation lines appear which, at variance with the previous ones, show no signature of non-separability.
Another important motivation of our work is the recent experimental study of Steinhauer [24] who studied an acoustic BEC black hole in one of the models discussed below (the so called "waterfall model" [25]) and presented results on entanglement similar to the ones discussed below.
The paper is organized as follows: Sec. 2 presents several black hole configurations in a quasi one dimensional BEC system. In Sec. 3 we compute the corresponding theoretical momentum correlation functions and the non-separability signals in the different configurations in a variety of situations: In particular we present the adiabatic and non-adiabatic regimes and also address in both cases the effects of temperature. These results are compared in Sec. 4 with the ones obtained in the absence of sonic horizon. In section 5 we discuss the limitations of our theoretical approach and finally we present our conclusion in section 6. Some technical points are presented in the Appendices: in Appendix A we discuss a rigorous windowed Fourier analysis which induces important restrictions to the measurement process; in Appendix B we give the form of the most general correlation functions and in Appendix C we discuss the specific case of a subsonic flow in the presence of a localized external potential.

Black hole configurations and their theoretical description 2.1 Quasi one-dimensional sonic black holes
In this work we consider a system where bosons are confined in one dimension by a harmonic transverse potential of angular frequency ω ⊥ . We denote by x the longitudinal degree of freedom and assume no trapping along x. In this configuration the theoretical description of the system is conveniently worked out in the quasi one-dimensional limit where the particles are described by a one dimensional (1D) quantum fieldΨ(x). According to the Bogoliubov prescription one writes the field operator as the sum of a main contribution (a classical field Ψ 0 ) and a small quantum remnantΨ (x) = Ψ 0 (x) +ψ(x) .
(1) Ψ 0 (x) describes the condensate order parameter and verifies the stationary 1D Gross-Pitaevskii equation with g 1d = 2 ω ⊥ a [26], where a is the 3D s-wave scattering length. In (2) µ is the chemical potential and U (x) a possible longitudinal external potential. In the absence of external potential, for a static homogeneous system of constant linear density |Ψ 0 | 2 = n one gets µ = g 1d n. Useful quantities are the sound velocity in the uniform system c = µ/m and the healing length ξ = /mc. In 1D a description based on Equations (1) and (2) is not quite legitimate, both in the high and in the low density limit: at large density, transverse excitations of the condensate cannot be discarded and the quasi 1D description (2) fails; at low density, phase fluctuations destroy the long range order and the possibility of a true Bose-Einstein condensation which is at the heart of the Bogoliubov description (1). In the remaining of this section we stick to the simple approach embodied by Equations (1) and (2) and we postpone the discussion of its limitations to Sec. 5.
We denote as a black hole configuration a 1D configuration in which the asymptotic upstream flow is subsonic (with constant density n u ) and the asymptotic downstream one is supersonic (with constant density n d ). Typically n u = n d and when a region of the flow is denoted for instance as subsonic, this means that in this region the density of the condensate is constant, and its velocity V u is smaller than the asymptotic sound velocity c u = g 1d n u /m.
Several black hole configurations have been proposed in Refs. [14,15,17]. The specific form of the order parameter is always of the type: where K u,d = mV u,d / , V u being the asymptotic upstream flow velocity and V d the downstream one (V u and V d are both positive). We also introduce the healing lengths ξ α = /(mc α ) and the Mach numbers M α = V α /c α (α = u or d depending if one considers upstream or downstream quantities). The functions φ u and φ d verify |φ d (x)| = 1 and lim x→−∞ |φ u (x)| = 1. The asymptotic upstream and downstream flows are respectively subsonic and supersonic, meaning that M u < 1 < M d . A remark on the location of the event horizon is in order here. First, as in any analogous system, its actual position is energy-dependent: we will even see n(x) n u → ←n d Figure 1: Waterfall configuration. The flow is incident from the left with an asymptotic density n u and a (subsonic) velocity V u . The downstream (x > 0) velocity V d is supersonic. The downstream density n d is constant and lower than n u . The region x > 0 is shaded in order to recall that it corresponds to the interior of the black hole.
below that the horizon disappears at large energy. This effect being taken into account, the customary procedure is to do a semi-classical analysis and to define as the true horizon the large wavelength one. In this case the horizon is the point where the flow velocity is equal to the local speed of sound. However, the flow varies rapidly around x = 0 in the configurations we study below, and a quantity such as the local speed of sound is ill defined in this region. As a result, the position of the event horizon cannot be unambiguously defined. This is not a drawback of the model: what really matters is that the asymptotic upstream and downstream flows are truly respectively sub-and super-sonic.

The "waterfall" configuration
We first consider one of the realistic configurations introduced in Ref. [25] and realized experimentally in Ref. [24]. In this configuration, which we denote as "waterfall", the 1D flow of a BEC is subject to an external potential which is a step function of the form U (x) = −U 0 Θ(x), where Θ is the Heaviside function (and U 0 > 0). In this case, a stationary profile with a flow which is subsonic upstream and supersonic downstream, i.e., a black hole configuration, has been identified in Refs. [24,27] and is schematically represented in Fig. 1. The upstream profile is exactly one half of a dark soliton and the downstream one corresponds to a flow with constant density and velocity (cf. Fig. 1). In the waterfall configuration one has φ d (x) = −i and φ u (x) = cos θ tanh(x cos θ/ξ u ) − i sin θ, where sin θ = M u . One has also U 0 /g 1d n u = 1 2 (M 2 u + M −2 u ) − 1 and V u = c d < c u < V d which indeed corresponds to a black hole type of horizon (M u < 1 < M d ).

The "δ peak" configuration
In this configuration the 1D flow of a BEC is subject to an external potential which is a Dirac distribution of the form U (x) = κ δ(x), where κ > 0. In this case, a stationary profile with a flow which is subsonic upstream and supersonic downstream, i.e., a black hole configuration, has been identified in Refs. [28,29], and it has been shown in Ref. [30] how this configuration can be reached dynamically. The downstream one corresponds to a flow with constant density and velocity and the upstream profile is a fraction of a dark soliton, with φ u (x) = cos θ tanh[(x − x 0 ) cos θ/ξ u ] − i sin θ, where sin θ = M u and x 0 depends on M u and κ (see details in [25]).

The "flat profile" configuration
We finally present a model configuration first introduced in Ref. [15], which has been denoted as "flat profile" in Ref. [25]. Although this configuration is not likely to be realized experimentally, it has been demonstrated in Ref. [25] that it yields a density correlation signal which is quite similar to the one obtained in the more realistic waterfall and δ−peak configurations. We will use the flat profile configuration to present our results below (in Sec. 3) for pedagogical reasons, because it leads to a simpler phenomenology for the momentum correlation than the other configurations.
In the flat profile configuration one has n u = n d ≡ n 0 and K u = K d ≡ K 0 and the φ α functions of Eq. (3) assume a very simple value: φ u (x) = φ d (x) = 1. This means that Ψ 0 (x) is a plane wave for all x. A horizon can still be realized in this case by tuning the values of the external potential U (x) and of the non-linear constant g 1d (x) such that These values are chosen so that a flow with Ψ 0 (x) = √ n 0 exp(iK 0 x) is solution of Eq. (2) for all x. This imposes We finally note that in the flat profile configuration one has c d < V d = V u < c u . This corresponds to a sonic black hole horizon since the upstream and downstream Mach numbers verify It is important to notice that, at variance with the cases of the waterfall and of the δ peak configurations, where, once an asymptotic flow is fixed (say, the upstream one) all the characteristics of the flow are uniquely determined, in the case of the flat profile configuration the values of M u and M d can be chosen independently one from the other. As a result, one cannot directly compare the results of, say the measure of non separability, of a waterfall and a δ peak configuration, but each of them can be compared with a flat profile configuration. This will be done in Figs. 4, 5, 6 and 7.

The excitation spectrum of a homogeneous condensate
In the case of a static homogeneous condensate, the dispersion relation of longitudinal excitations is the standard Bogoliubov one: u|in → ←u|out out in out in Figure 2: Dispersion relations. The left plot corresponds to a subsonic flow. The right plot corresponds to a supersonic flow; it is shaded in order to recall that it describes the situation inside the black hole. In each plot the horizontal dashed line is fixed by the chosen value of ω. The q (ω)'s are the corresponding abscissas, with ∈ {u|in, u|out} in the left plot and ∈ {d1|in, d1|out, d2|in, d2|out} in the right plot. The direction of propagation (left or right) of the eigen-modes is represented by an arrow. The lower diagram illustrates schematically the "ingoing" and "outgoing" terminology used in the text and in the two upper plots.
In a region where the condensate flows with a velocity V this is modified to In this case ω is the energy of the elementary excitation in the frame where the obstacle is at rest, while ω B is the frequency measured in the frame of the fluid. The momentum of the excitation relative to the background flow is q, and its momentum in the frame of the obstacle is q +mV . In a black hole configuration, the upstream and the downstream channels are characterized by dispersion relations of the type (7) (with the appropriate values of V , ξ and c), they are illustrated in Fig. 2. In this figure the upper left (upper right) plot represents the asymptotic upstream subsonic (downstream supersonic) case. The part of the dispersion relation represented by a dashed line correspond to negative norm modes, as explained in the following section. In all this work we only consider the ω > 0 part of the dispersion relations, this is made possible by the ω ↔ −ω symmetry of the spectrum [31]. One sees in the figure that, upstream, the waves of one of the channels are directed towards the horizon (ingoing waves, denoted as u|in) whereas the waves in the other channel propagate away from the horizon (outgoing waves: u|out). The definition of which mode is ingoing and which is outgoing depends on the side of the horizon that is considered, this is illustrated by the lower diagram of Fig. 2. Downstream, there are two ingoing waves (d1|in and d2|in) and two outgoing waves (d1|out and d2|out). Note that the two d2 channels disappear at large energy, when ω > Ω, see Fig. 2. Ω is the energy of an elementary excitation whose group velocity in the frame where the condensate is at rest (∂ω B /∂q) is equal to the flow velocity V d (such an equality is only possible for supersonic flows). Excitations with momentum larger than the one of this excitation (which is denoted as q * in Fig. 2) move faster than the flow and are able to escape the "black hole".

The wave function in real space
The field operatorψ(x) associated in the Schrödinger representation to the particles which are out of the condensate [as defined by Eq. (1)] is expanded over the scattering modes: s create an excitation of energy ω in one of the three scattering modes (L = U , D1 or D2), they obey the following commutation relation: Each of the three scattering modes (U , D1 or D2) is initiated by one of the three entrance channels (u|in, d1|in or d2|in). Far from the horizon, the density and the velocity of the flow are position-independent and the corresponding wave functions are mere plane waves of the following form: • Deep in the upstream subsonic region, i.e., for x < 0, x −ξ u : • Deep in the downstream supersonic region, i.e., when x > 0, x ξ d : Note that theū L 's,w L 's, q 's, S i,j 's, U 's and the W 's in Eqs. (10) and (11) all depend on ω. For instance q (ω) is defined on Fig. 2 for ∈ {u|out, u|in, d1|out, d1|in, d2|out, d2|in}. We chose a normalization of the the coefficients U and W -the so called "Bogoliubov amplitudes" -such that The sign + (−) in (12) refers to positive (negative) norm modes. All the upstream modes (u|in and u|out) have a positive norm. In the downstream region, the d1|in and d1|out modes have a positive norm while the d2|in and d2|out ones have a negative norm. The normalization (12) ensures that a positive (negative) mode carries a current +1 (−1). The explicit expression of the the coefficients U (ω) and W (ω) corresponding to the normalization (12) is given in Ref. [25]. Expressions (10) are not valid close to the horizon due to (i) the modification of the density profile which is position-dependent in vicinity of the horizon 1 and (ii) to the occurrence of evanescent modes (with complex momenta solutions of Eq. (7)) which are of importance near the horizon. Let us for instance discuss the physical content of the last of Eqs. (11). It describes a d2|in wave incoming from +∞ (last term of the r.h.s., the corresponding group velocity is negative, cf. Fig. 2) which is back scattered into a d1|out and a d2|out wave, with respective reflection amplitudes S d1,d2 and S d2,d2 . Part of this wave is also transmitted in the x < 0 region as a u|out wave with transmission amplitude S u,d2 : the corresponding expression far from the horizon is displayed in the last of Eqs. (10). The scattering amplitudes form the S matrix which is 3 × 3 for energies ω lower than the threshold Ω defined on Fig. 2: Current conservation reads For ω > Ω, the last row and the last column of (13) vanish because the d2 mode disappears.
In this case the S matrix is 2 × 2 and satisfies SS † = 1.
In Eqs. (10) and (11) we did not write the contribution of the evanescent modes since they decay exponentially and are negligible far from the horizon (when |x| ξ (u,d) ), but we fully take these modes into account in the expression of theū's and thew's near the horizon (around x = 0); this is needed for an accurate computation of the S matrix. In a given configuration (say "waterfall"), the elements of the S matrix are determined for each value of ω by enforcing continuity of the wave functions (and of their spatial derivatives) of the elementary excitations at x = 0. This represents an easy numerical task which consists in solving a linear 4 × 4 system for each value of ω.

The wave function in momentum space
Because of the presence of the negative norm/negative energy d2|in mode, stationary black hole configurations such as the ones presented in section 2.1 are meta-stable 2 . Indeed, when reaching the horizon, the d2|in waves give rise to a radiation of u|out quanta in the upstream region, which constitutes the spontaneous Hawking radiation. In BEC systems however, this radiation is not easily detected. The reason is that the occupation of the Hawking radiating modes is approximately of thermal type, with an effective temperature much lower than the true temperature of the system (typically by a factor 10), and the Hawking signal is thus drowned in the thermal noise. This circumstance led the authors of Refs. [14,15] to look for density correlations as alternative evidence of the Hawking effect.
The idea, checked in Refs. [14,15] is that outgoing waves generated by the same d2|in mode are all correlated. Moreover, since the corresponding amplitudes are governed by the S matrix which describes how ingoing waves hitting the horizon generate outgoing ones, the knowledge of the S matrix makes a detailed description of the correlation signal possible. This point has been checked thoroughly in Refs. [17,18,25]. In particular, the characteristic peaks of the density correlations correspond to the Hawking quantum (u|out) -partner (d2|out) correlations (for points situated on both sides of the horizon), and also to correlations along the u|out − d1|out (again, for points situated on both sides of the horizon) and d2|out − d1|out (for both points inside the horizon) channels.
In a BEC, momentum correlations could be more precisely detected than density correlations, by following a procedure used in Ref. [22] in a similar context, in the case of the dynamical Casimir effect. This is the reason why we proposed in Ref. [19] to demonstrate the existence of sonic Hawking radiation by the means of correlations in momentum space.

A local Fourier transform
By appropriately introducing a local Fourier transform in both the subsonic (exterior of the black hole) and supersonic (black hole interior) regions, we shall explicitly construct the momentum correlator and analyze, in Sec. 3, its nontrivial qualitative features, which are in correspondence with those present in the density correlation signal. This signal concerns the occupation number in the momentum representation:N (K) =ψ † (K)ψ(K), whereψ(K) is the Fourier transform ofψ(x) 3 : From expression (8) and the mode analysis presented in Secs. 2.2 and 2.3, it is clear that the momentum distribution has a different form in the far-upstream and far-downstream regions. Hence, instead of (15), it is more appropriate to perform a specific mode analysis in each of these regions [32][33][34]. This can be done by using a window function selecting the desired region of space. The precise form of this window function is irrelevant, but for concreteness we will consider a Gaussian. In the upstream region for instance, one takes a window and the corresponding windowed Fourier transform iŝ This procedure is meant to select the momentum components which can be identified from (10). For this purpose, the parameters of the window function have to be chosen in order to work in the appropriate region of space. This is done by taking X u < 0, σ u > 0 and letting X u and σ u respectively tend to −∞ and +∞, imposing for the ratio X u /σ u → −1 which allows that X u + σ u = C st −ξ u . This procedure is illustrated in Fig. 3. It is important to take the limit of large σ u and |X u | for obtaining a sharp (δ-like) distribution in momentum increasing |X u | and σ u Figure 3: Schematic representation of the behavior of the parameters X u and σ u of the window function Π u (x) defined in Eq. (16). The shaded zone is the region where the window function notably differs from zero. The spacing between X u + σ u and the origin is large compared to ξ u , and this ensures that the Fourier analysis, which is made local thanks to the window function, is performed in the deep subsonic region.
space. In Eq. (16) the normalization parameter Λ u is introduced to effectively describe the finite efficiency of the experimental detection apparatus. More precisely, it is argued in Sec. 4 that the quantity λ u = σ u Λ 2 u / √ 2π describes the rate of detection of particles in the window [X u − σ u , X u + σ u ]: if λ u = 0 none of the particles are detected, if instead λ u = 1, they are all detected. It is interesting to stress that the explicit results given in Sec. 3 for the normalized momentum correlation function (29) do not depend on the efficiencies λ u and λ d of the detectors.
Of course, a local Fourier transform similar to (17) is performed downstream, using a different window Π d (x) with parameters Λ d , X d (> 0) and σ d , leading to the downstream Fourier transformψ d (K). We will see in Appendix A that the parameters of the upstream and of the downstream window function have to be chosen in a specific manner when one wants to consider a specific correlation signal. However, these precise conditions -which are of importance for the theoretical analysis of the experimental detection scheme -can be replaced, once met, by the following schematic, but natural rules: These rules are less rigorous than the correct mathematical procedure (17) which uses window functions for defining the local Fourier transforms, but it is shown in Appendix A that, provided some simple and natural redefinitions (in particular of the singular Dirac distributions) are considered, they yield the correct result. Hence, we shift discussion of the more rigorous results to Appendix A and present our results in the main text using these simplified rules. Performing the schematic Fourier transform in the x < 0 region we get, whereas in the x > 0 one getŝ In the two above integrals the q 's are functions of ω computed as schematically represented in Fig. 2. The ω integration yields factors |∂q /∂ω| which can be absorbed in a re-definition of the U's and of the W's: where ω (q) is the reciprocal function of q (ω). The "tilde Bogoliubov coefficients" satisfy the following normalization: Note that in the integrals defining the expressions (18) and (19) ofψ u (K) andψ d (K), all the terms (Bogoliubov coefficients and coefficients of the S-matrix) involving a d2-subscript cancel when ω > Ω. Defining and for the downstream Fourier transform [instead of (19)], All the terms in the expressions (22), (23), (24) and (25)

The case of the flat profile configuration
In this configuration -presented in Sec. 2.1.3 -the fact that K u = K d ≡ K 0 makes it possible to express (22), (23), (24) and (25) It is then possible -and useful, see sec. 3 below -to regroup the momentum operators in Eqs. (22), (23), (24) and (25) in terms of k < 0 and k > 0 contributions, and to writê Each of the above expressions (k < 0 and k > 0) contains terms coming from both the subsonic and the supersonic regions. This procedure can induce a problem of normalization. As discussed in Sec. 4 after Eq. (92) this problem is easily solved by using an appropriate overall multiplicative factor which we do not include for readability. Furthermore, the problem disappears when one considers normalized quantity such as g 2 (K, Q) defined below [Eq. (29)].

Momentum correlations in the presence of a sonic horizon
The momentum-momentum correlation signal is embodied in the function whereN (K) =ψ † (K)ψ(K). We also consider in some details the normalized quantity In the definitions (28) and (29) the normally ordered product eliminates the diagonal shot noise contribution. The computation of the two-body momentum correlation (28) in the generic case is quite tedious because the upstream and downstream Fourier transform are different, as explained in Sec. 2.4.1. Besides the formulas assume different forms depending of the values of K and Q relative to K u and K d , cf. Eqs. (22), (23), (24) and (25). As a result, one has to consider nine different cases. Although we will treat all the black-hole configurations presented in Sec. 2.1 (plots encompassing all the different cases for the waterfall, δ peak and flat profile configurations are given in Figs. 4 and 5), to simplify the presentation we give here the explicit results in the "flat profile" configuration where the background density is a uniform plane wave (cf. Sec. 2.1.3). In this case one has to consider only 4 cases: the four quadrants in the (k, q) plane (where k = K − K 0 and q = Q − K 0 ) and one can rewrite G 2 in terms of k and q and then perform explicit computations using the expressions (26) and (27). The theoretical evaluation of the momentum correlation function is performed in order to match with an experimental detection scheme which consists of opening the trap and letting the elementary excitations be converted into particles expelled from both ends of the condensate, according to a process known as "phonon evaporation" [35]. If this process is assumed to be gentle and adiabatic, each elementary excitation is converted into a single particle. More precisely, one has, in this case, for the positive norm u and d1 modes: while for the negative norm d2 modes: With the normalization (12), this corresponds, for instance for a positive norm mode, to a perfect transmutation of an elementary excitation into a particle state which carries a unit current. For the "tilde Bogoliubov coefficients" (20), the prescription (31) yields for the positive norm states ( ∈ {u|in, u|out, d1|in, d1|out}) whereas one gets for the negative norm modes ( ∈ {d2|in, d2|out}).
As demonstrated in Ref. [22], after an adiabatic opening of the trapping potential, a measure of the velocity distribution of the emitted particles gives access to the momentum distribution within the condensate and to the correlators defined in Eqs. (28) and (29). The process can be sudden, in which case the adiabatic hypothesis breaks down. More precisely, if t open is the characteristic time during which the trap is opened and an elementary excitation is converted into a real particle, the adiabatic approximation fails for excitations of energy ω verifying: ω −1 t open : for low lying excitations (those for which ω → 0) the opening of the trap always seems abrupt [35]. A quantitative study of this phenomenon will be presented in a future publication [36].
An outline of the complications introduced by non adiabatic effects is postponed to Section (3.3), but until then we give the results after an adiabatic expansion, in which case, owing to (33) and (34), the expression ofψ given in (26), (27) andψ In expressions (35) and (36) the "tilde Bogoliubov coefficients" which are non zero assume the limiting values (33) and (34).

Zero temperature
Let us first consider the case where the system is initially in its vacuum state, i.e., the vacuum of excitations. This is the zero temperature case. We find for the one-body momentum distribution: and These relations can be cast under the form where Note that within the adiabatic approximation, the T = 0 momentum signal (40) would cancel in the absence of horizon, since the d2 mode and all the corresponding elements of the S-matrix would disappear in this case. This is confirmed by the study of Sec. 4. In expression (39), the δ-peak contribution is singular: one has a δ(0) [as in (37) and (38)]. This is due to the schematic nature of the rules R1 and R2 defined in section 2.4.1 for the Fourier transform. The more rigorous local Fourier transform in terms of a window function -explained in the first part of Sec. 2.4.1-yields a finite contribution, as shown in Appendix A: see, e.g., Eqs. (108) and (109) which are the rigorous versions of Eqs. (37) and (38).
Our main interest in this work is the study of the correlation signal in momentum space, that is of the two-points correlation functions G 2 and g 2 defined in Eqs. (28) and (29). The most robust signals are those present even at T = 0. These are As already noted for the one-body signal, a first obvious outcome of this computation is that the T = 0 correlations in momentum space disappear in the absence of sonic horizon, since in this case the d2 mode does not exist and all the corresponding elements of the S matrix cancel in (41), (42) and (43). Looking more in detail into the results, one sees that in the (k < 0, q < 0) and in the (k > 0, q > 0) sectors one has terms with a δ 2 (k − q) correlation which we henceforth denote as diagonal. These terms are simply of the form N 2 (k)δ 2 (k − q). Again, the occurrence of the highly singular squared δ distribution in the expressions (41), (42) and (43) is an artifact of our schematic Fourier transform rules R1 and R2 which disappears when one considers windowed Fourier transforms, see Appendix A.
Besides the diagonal term, one has correlation lines along the u|out − d2|out (Hawking quanta-partner), d2|out − d1|out and u|out − d1|out channels. All these correlations are also present mutatis mutandis in the density fluctuation sector [14,15,18,25]. As in the interpretation of black-hole radiation first given by Hawking [3], the existence of these correlations is an indication of the fact that the vacuum of the ingoing modes is not the same as the vacuum of outgoing ones, and this results in spontaneous emission of pairs of outgoing quasi-particles even in the absence of incoming ones [37,38]. This is possible even in our stationary setting because the number of created quasiparticles of positive energy equals the number negative energy ones, i.e. energy is conserved [25,39].  The results corresponding to the zero temperature correlation signals (41), (42) and (43)  The graphical rule for drawing the correlation lines in these figures is simple. Let's consider the u|out − d2|out correlation line for instance. From Eq. (41) one sees that the relative wave vectors k and q are correlated when ω u|out (k) = ω d2|out (−q), and also along the curve obtained by exchanging the roles of k and q. This corresponds to the two curves q = −q d2|out (ω u|out (k)) and q = q u|out (ω d2|out (−k)). Going to the absolute wave vectors K and Q, these correlation lines correspond to the curves Q − K d = −q d2|out (ω u|out (K − K u )) and Q − K u = q u|out (ω d2|out (−K + K d )). These two curves are symmetrical with respect to the diagonal. In the case of the flat profile configuration they meet on the diagonal at the point of coordinate (K u , K d = K u ). Figures 4 and 5 indicate the location of the relevant correlation signal, but not its amplitude. The theoretical evaluation of this amplitude can be done either using the schematic rules R1 and R2 [but then yields to singular expressions as in (41), (42) and (43)] either using windowed Fourier transforms (but then depends on the specific choice of the windows, cf. Appendix A). One can circumvent these difficulties and obtain a non-ambiguous result by working with the rescaled g 2 function (29), as explained now.
Let us first consider "diagonal correlation terms" which are intra-channel correlation in the u|out, d1|out and d2|out channels corresponding to the diagonal lines in Figs. 4 and 5. These diagonal terms are conveniently isolated by studying correlation functions of the type where N (K) u|out and G 2 (K, K) u|out are the u|out contributions to N (K) and to the diagonal part of G 2 (K, K). One obtains straightforwardly a result which follows from Wick's theorem and which is also valid for the other channels, even at finite temperature, for non-adiabatic opening of the trap, and also when considering Fourier transforms less schematic than in Eqs. (18) and (19) (cf. Appendix A). At zero temperature, in the adiabatic regime we consider in the present subsection, the only non-diagonal contributions to G 2 are of the type u|out-d2|out, u|out-d1|out and d1|out-d2|out. One considers here for instance correlation functions of the type (with obvious notations) From expressions (37) and (41) one gets For obtaining the r.h.s. of Eq. (47) we used the pseudo-unitarity condition (14) and also the fact that While this relation is easily verified in the schematic framework used in the main text (which originates from the rules R1 and R2 of Sec. 2.4.1), it appears less straightforward when studying the more rigorously defined local Fourier transform (see Appendix A). In this case, the equivalent of Eq. (48) is Eq. (124) and is equal to unity only if Eqs. (126) are verified, i.e., if the window functions fulfill specific conditions, the physical content of which is discussed in Appendix A. We will assume that these conditions are met in the following (or equivalently we keep on using the schematic rules R1 and R2 of Sec. 2.4.1). For the other inter-channel correlators the normalized two-body functions read: and The study of g 2 is of interest because the occurrence of entanglement and the quantum nature of the Hawking process can be tested through the violation of the Cauchy-Schwarz inequality, as recently studied in a similar context in Refs. [32,33,[40][41][42][43][44][45][46]. For instance, the Cauchy-Schwarz inequality is violated along the characteristic Hawking quanta-partner correlation lines u|out-d2|out of Figs. 4 or 5 if (see, e.g., [47]) As already noticed after Eq. (45), the right-hand side of inequality (51) is equal to 2 for all temperature. From expression (47) one sees that the Cauchy-Schwarz inequality is violated at T = 0 along the u|out-d2|out correlation line for those values of K and Q such that, at energy ω = ω u|out (k) = ω d2|out (−q), one has |S d2,d2 (ω)| > 1. S d2,d2 is the scattering amplitude from the d2|in mode towards the d2|out mode; its modulus can be larger than unity without violating the pseudo-unitarity condition (14). |S d2,d2 (ω)| diverges as ω −1/2 when ω → 0 [25] and the Cauchy-Schwarz inequality is thus always violated for ω > 0. The same holds true for the violation of the Cauchy-Schwarz inequality along the d1|out-d2|out channels at T = 0 since the formula (50) yields for g 2 (K, Q) d1|out−d2|out the same result than (47) does for g 2 (K, Q) u|out−d2|out . From Eq. (49) it is clear that the Cauchy-Schwarz inequality is not violated at T = 0 along the u|out-d1|out correlation channel. A remark should be made here concerning the experimental detection process. Our theoretical analysis corresponds to windowed upstream and downstream momentum detection. It is thus theoretically possible to distinguish an upstream and a downstream component in the signal in momentum space. This is the reason why it is legitimate to identify, for instance a u|out − d2|out component in the total G 2 (K, Q), or a u|out component in the total N (K) , as done in Eq. (46). However, some apparatuses may mix the upstream and downstream signals in their detection scheme. In this case, the overlap in momentum space of the domain of existence of the u|out and d2|out signals forbids the use of a definition as precise as (46). In this case, a clear separation can only be done for the momenta located outside of the overlap region. In the waterfall configuration represented in the left plot of Fig. 4 for instance, with such a detection apparatus, one cannot study the u|out − d2|out correlation signal in regions where the d2|out and u|out momenta overlap, i.e, when K or Q lie in a segment for which the blue and red diagonal curves of the figure overlap. In this respect, such an apparatus would not have access to the u|out − d2|out correlation signal at all for the flat profile configuration, since in this case the region of momenta issued from the u|out channel is completely contained in the region of momenta issued from the d2|out channel. To recall this possible experimental issue, when we plot below quantities such as g 2 (K, Q) shade the region of momenta where an ambiguity is possible (cf. Figs. 6 and 7). We stress however that not all experimental detection schemes have to be subject to this drawback.

Finite temperature
The analog stationary black hole configuration we consider is thermodynamically unstable, and cannot support a thermal state. However, a thermal-like occupation of the states can be defined by considering a time dependent process of formation of the horizon starting from a uniform thermal subsonic uniform configuration, following the procedure already considered in Ref. [18] (see also [17]): The condensate has initially a uniform density n u , a flow velocity V u and a sound velocity c u and is at thermal equilibrium at temperature T in the moving frame. The horizon is then adiabatically switched on, either by ramping down the scattering length in the downstream region -leading to a flat profile configuration -, or by ramping up an external potential -leading to a waterfall or a delta peak configuration.
Note that one could choose other prescriptions for defining the occupation of the scattering modes. For instance the initial state (uniform density n u with uniform velocity V u ) could be in thermal equilibrium in the frame of the obstacle (and not in the frame where the condensate is at rest, as considered above); this would modify the explicit expression of the n L 's in (52). Precisely one would have in this case: n U (ω) = n th (ω), n D1 (ω) = n th (ω u|out (q d1|in (ω))) and n D1 (ω) = n th (ω u|out (−q d2|in (ω))). In the following we just use the contraction rules without specifying the expressions of the n L 's, and the formulas written below are thus generally valid.
In the adiabatic regime where Eqs. (33) and (34) hold, we find, for negative k: and, for positive k: As in the zero temperature case, formulas (53) and (54) can be cast under the form (39), with here and Then, one obtains for the momentum correlation: In the first term of the last line of this formula there is no ambiguity: |S d1,u | 2 and n U have to be evaluated at ω d1|out (k) = ω u|in (q), | U u|in | 2 has to be evaluated at q and | U d1|out | 2 has to be evaluated at k.
We remark that at T = 0 there appear in-out correlators, which were absent at T = 0. We shall focus here on the most robust correlators, those already present at T = 0, and evaluate the intensity of the correlation signal along these lines by using the normalized correlation function g 2 which is window independent. The Hawking quanta-partner correlator (47) gets modified to where [from (55)] and The u|out − d1|out correlator (49) becomes and the d1|out − d2|out correlator (50) takes the form where [from (56)] Eqs. (60), (63) and (64) are particularly important because they allow us to study how an initial nonzero temperature affects the violation of the Cauchy-Schwarz inequality, g 2 > 2, see e.g. Eq. (51). The results in the waterfall, δ-peak and flat profile configurations are presented in Figs. 6 and 7. We saw in the previous subsection that at T = 0 the Cauchy-Schwarz inequality is always violated along the u|out − d2|out and d1|out − d2|out channels for all ω > 0, while it is not violated in the u|out − d1|out one. In a more realistic case in which there is an initial nonzero temperature, the amount of entanglement is, as expected, reduced with respect to the T = 0 case. In particular, there is always a region, for small enough and large enough momenta (when the corresponding frequency is close to 0 or to Ω), where g 2 < 2. In both the u|out − d2|out and d1|out − d2|out channels, however, there is always an intermediate ω region in which the Cauchy-Schwarz inequality is violated provided the temperature is not too large.
We clearly see that the most favorable case for demonstrating quantum entanglementi.e. spontaneous Hawking radiation -corresponds to the violation of the Cauchy-Schwarz inequality along the Hawking quanta -partner (u|out − d2|out) channel in the waterfall configuration, which is exactly the situation which has been recently studied experimentally by Steinhauer [24]. We remark that this configuration corresponds to the case where the region of overlap of the signals in momentum space is the smallest.
We also note from Fig. 7 that, for a δ-peak potential, the Cauchy-Schwarz inequality is violated up to a temperature slightly larger than the chemical potential. This has been already noticed in Ref. [32]. The signature of quantum correlation is more robust for a waterfall configuration: one sees on Fig. 6 that the Cauchy-Schwarz inequality is violated up to a temperature of order of twice the chemical potential.

Non adiabatic effects and in situ measurements.
Since the adiabatic limit described by Eqs. (33) and (34) is rather idealized, especially for small frequencies (as discussed earlier, for low lying excitations the assumption of adiabaticity is always violated), it is useful to describe how the results presented in the previous two subsections are modified by non-adiabatic effects. Concretely, instead of Eqs. (33) and (34), at late times after the opening of the trap, the Bogoliubov coefficients will behave as where the α and β coefficients satisfy the same normalization condition (21) as the corresponding U, W "tilde Bogoliubov coefficients" (20). Note that our approach encompasses also the situation which we denote as in situ below. In this situation the measurement process gives a direct access to the actual value of the tilde Bogoliubov coefficients without any evaporative process which would affect these coefficients as described in Eqs. (33) and (34) for the adiabatic limit, or possibly differently in a non adiabatic regime. In the following, we use the generic terminology "non adiabatic" for describing both the case of in situ measurements and of non-adiabatic modification of the situation (33), (34). To ease the presentation we have -only temporarily -considered in (66) a single set of coefficients for the upstream u region and also a single one for the downstream d region, but it is clear that in a rigorous analysis (presented below) each mode will have its own α and β coefficients. Hence the expressions (26) and (27) remain valid, but the U and W's are replaced by α and β coefficients, following the rules (66).
The reason for the change of notation (66) is threefold. First, the expression of the nonadiabatic coefficients maybe quite non-trivial and different from the ones of the U's and W's. This will occur for instance after a step of dynamical Casimir amplification where the system is artificially submitted to a rapid quench. Also, this notation allows for the possibility that, during the opening of the trap, the initial "tilde Bogoliubov coefficients" get modified in a manner less trivial than (33) and (34). Second, this (momentarily) simplified notations where the mode indices are omitted permits a simple presentation of the main features of the G 2 function in the non-adiabatic case (see Appendix B). Finally, in the weakly non adiabatic case all the β's are small compared with the α's, whereas keeping the previous U's and W's notations would make it difficult to identify the small terms in the expression for G 2 .
The results for the correlator are much more complex than those presented in the previous two subsections in the adiabatic limit. We give here the general structure of the correlators, the explicit results -with also account of finite temperature -are given in Appendix B: We recall that in the adiabatic limit, corresponding to Eqs. (33) and (34) with the substitution (66), all the β's are zero. In the more general case considered here they are not; the α's and β's, with the substitution (66) -satisfy the normalization (21) -implying that α is bigger than its adiabatic value. We see from the results presented in Appendix B that the terms already present in the adiabatic regime are now multiplied by a factor α 4 . In particular, the finite temperature adiabatic terms of (57), (58) and (59)  . If we consider the weakly non adiabatic regime, this means that they are now larger than the corresponding adiabatic value. New sub-leading terms appear, of order α 2 β 2 terms (among which an antidiagonal term) and also higher order β 4 contribution. This results in a very complicated pattern.

Violation of the Cauchy-Schwarz inequality along the Hawking quantumpartner correlation lines
Having derived the structure of the correlation lines, we now turn to the intensity of the correlation signal, and to the possible violation of the Cauchy-Schwarz inequality in the non adiabatic regime. We shall restrict our attention to the study of the incidence of non-adiabaticity on the Hawking quantum -partner u|out-d2|out correlator; this will bring pieces of information valid for all the other correlation lines. The results will be given first at T = 0 4 , then at finite temperature for completeness.
• From the results in Appendix B, in the case where both k and q are negative, the zero temperature contribution of the u|out-d2|out modes to G 2 reads [see Eq. (128)] The arrow in this equation indicates that its right hand side is not the sole contribution to the u|out-d2|out correlation signal: not only should it be supplemented by a contribution in which the roles of k and q are exchanged, but also new non-adiabatic terms arise (see below). The contribution (70) corresponds to the signal which already exists in the adiabatic case, whose intensity is here modified by the α coefficients. Note that, contrarily to the schematic presentation of Appendix B, we consider here the most general case, and have explicitly written the mode-dependence of the α coefficients. This will remain the case in the rest of the section.
In the negative k sector, the zero temperature contributions of the u|out and d2|out modes to N (k) read In the calculation of g 2 (K, Q) u|out−d2|out the α coefficients of Eqs. (70) and (71) factorize out and we get for the normalized correlator the same result as in zero temperature adiabatic case, Eq. (47), namely The same factorization of the α coefficients occurs also at finite temperature, and the adiabatic expression (60) is thus also valid in the non-adiabatic regime. We stress that this important signal is thus quite robust, not being affected by possible non adiabatic effects, and the violation of Cauchy-Schwarz inequality is also observable from in situ measurements. From the results presented in Appendix B one sees that in the non-adiabatic regime three more couples of u|out-d2|out correlation lines appear with respect to the adiabatic situation. They correspond to the conditions ω d2|out (k) = ω u|out (−q) for k > 0 and q > 0 , ω d2|out (q) = ω u|out (k) for k < 0 and q > 0 , and to the same conditions in which the roles of k and q are exchanged. All these correlation lines are displayed in Fig. 8, together with the ones which already exist in the adiabatic case 5 . We now evaluate the intensity of the signal corresponding to each of these lines.
• In the k and q > 0 sector, we get, at T = 0 [cf. Eq. (134)]: and leading to the following contribution to the normalized correlator: As was previously the case for the α coefficients, the β's here also factorize out. A similar factorization will occur in all the subsequent cases considered in this section. At finite temperature, expression (76) becomes where here The important information here is that, even at zero temperature, the new u|out-d2|out correlation line which appears in the k and q > 0 sector due to non-adiabatic effects is not associated to a non separable signal [cf. Eq. (76)]. As expected -and can be verified from expression (77)-thermal effects do not modify this situation. As we will now see, this conclusion remains also valid for all the other correlation lines which where not present in the adiabatic regime.
• For k < 0 and q > 0, one of the contributions to G 2 (k, q) u|out−d2|out reads [cf. Eq. (137)] and at T = 0 this corresponds to a normalized correlation signal: At finite temperature this contribution modifies to: • For k > 0 and q < 0 another contribution to G 2 (k, q) u|out−d2|out reads [cf. Eq. (137)] giving at T = 0 This form of writing the result has been obtained by using the pseudo-unitarity of the S matrix. It makes clear that no violation of the Cauchy-Schwarz inequality occurs along the correlation line considered here. At finite temperature one obtains where here To sum up, we recall that the above study is a partial focus on a subpart of the whole correlation pattern, concerning only the most important Hawking quantum-partner signal (in our terminology, the u|out − d2|out signal). We used it to demonstrate that the most interesting correlation is the one already present in the T = 0 adiabatic case: the Cauchy-Schwarz inequality can only be violated along this line. The same is true for the d1|out−d2|out signal.

Momentum correlations in the absence of sonic horizon
In this section we consider a configuration where the upstream and the downstream regions are both subsonic. In this case there is no horizon, but one can still be in a configuration where the upstream and downstream non linear coefficients are different (as in the flat profile configuration), also the system can be affected by the presence of an external potential (as in the waterfall and delta peak configuration). If this external potential is localized (i.e., tends rapidly enough to zero at infinity), and if the nonlinear coefficient keeps the same value in all the system, then the type of flow considered is rather simple. More precisely, as demonstrated in Appendix C, the upstream flow velocity and density at −∞ are the same as the flow velocity and density at +∞: V u = V d and n u = n d . As a result, the general formulas given in the present section simplify, this is explained in Appendix C.
In the situation we consider in the present section, the d2 negative-norm mode disappears since the downstream region is subsonic: there is now a single downstream mode which we simply denote by "d". The obstacle is thus characterized by a S-matrix which is 2 × 2 and unitary: In this case, instead of Eq. (8) one haŝ where K α = mV α / as in (8), V α being the value of the upstream (α = u) or downstream (α = d) asymptotic flow velocity. We perform the usual (fake) Fourier transform onψ, using the rules R1 and R2 of section (2.4.1). Collecting separately the k < 0 and k > 0 contributions we get ψ(k < 0) = U u|out S u,ubU + S u,dbD ω u|out (k) This expression corresponds to Eqs. (26) and (27) in which the negative norm d2 mode has been suppressed. From it one can compute the density distribution in momentum space and also the momentum correlation function. In particular, in the adiabatic limit one gets and Since one is in a situation where the flow is everywhere subsonic, one can define a bona fide temperature state where, for all the modes, the occupation number of a state of energy ω is n U = n D = n th (ω). In the idealized case of perfect transmission one has S u,u = 0 = S d,d and S u,d = 1 = S d,u , then Eqs. (90) and (91) reduce to N (k) = 2 n th (ω(k)) .
The factor 2 is spurious. It comes from the fact that one does two Fourier transforms, one upstream and one downstream, and that there is a kind of built-in double counting in this approach. The problem is suppressed if one considers upstream and downstream windowed Fourier transforms, as done in Appendix A. In this case, in the absence of the obstacle, instead of Eqs. (90) and (91) one gets How the correct treatment of the windowed upstream and downstream Fourier transforms leads to the precise form of the two terms in the big parenthesis of the r.h.s. of (93) is explained in Appendix A.
From the present analysis one is thus led to define finite efficiencies of the upstream and downstream particle detectors. A finite efficiency corresponds to a measurement process in which a fraction of particles are missed in the detection of the momentum signal. This is described theoretically by the normalization of the windowed Fourier transforms (16): defining the efficiencies as λ α = σ α Λ 2 α / √ 2π (α = u or d) with λ α ∈ [0, 1] one casts formula (93) under the form N (k) = 1 2 (λ u +λ d )n th (ω(k)). So, when both efficiencies are zero, there is no particle detected and no momentum signal, and when, on the contrary, the detection efficiencies are both unity one gets in the absence of the obstacle N (k) = n th (ω(k)), as it should be. Note that the double counting which has been easily identified in formula (92) is present in all the previous formulas of main text (Eqs. (53) to (56) and also Eqs. (57) to (59)). It can be cured in the same simple way. For instance, considering perfect detectors, one should have added a factor 1/ √ 2 in the definition of rules R1 and R2 of section (2.4.1) for the schematic Fourier transforms. We did not do that for avoiding an overall multiplicative factor in all the formulas, and also because, as demonstrated in Appendix A this double counting (and also the finite efficiency of the detector) disappears in the formulas for the normalized correlation signal g 2 which is our main interest in the present work.
Finally, we give the expression for the momentum correlator, for simplicity in the adiabatic regime These expressions can be further simplified by considering, for this configuration which is everywhere subsonic, a common occupation number n th (ω) of a state of energy ω. As expected the corresponding expressions in the black hole case, Eqs. (57)(58)(59), reduce to the above in the absence of the negative norm d2 modes. Also, unlike in the black hole case all contributions disappear in the case where the initial state is the vacuum.

Limitations of the theoretical description
In this section we properly define the domain of applicability of our approach. A first limitation concerns the low density regime. As well known, in one dimension phase fluctuations prevent a true Bose-Einstein condensation. As a result, a description based on the simple separation (1) between a classical field and a small quantum correction does not allow to properly estimate the large x behavior of the one body density matrix Ψ † (x)Ψ(0) of a homogeneous 1D system (see, e.g., Ref. [50]) and, at low density, phase fluctuations blur the sharp correlations of Figs. 4 and 5, cf. Refs. [51] and [52]. We nonetheless argue that (1) it is still useful for understanding the qualitative behavior of some observables important for analyzing correlations in the system: for instance (1) yields the correct two-body correlation Ψ † (0)Ψ † (x)Ψ(x)Ψ(0) of a homogeneous system [53]. The relevance of (1) depends on the characteristic length involved in the spatial correlation considered. If the correlation characteristic length is smaller than the phase coherence length L φ = ξ exp π n 2maω ⊥ then (1) may be used. This is what happens for the two-body correlation : this quantity is non trivial (i.e., different from the square n 2 of the linear density n = Ψ † (0)Ψ(0) ) only in a range of distances x < ξ, typically much smaller than L φ . More precisely, L φ is exponentially larger than ξ, and the separation (1) is thus valid, when where a ⊥ = /mω ⊥ is the transverse harmonic oscillator length. In the large density limit the 1D description also fails, not because of lack of BEC as just discussed, but because the transverse degrees of freedom of the system are not completely frozen. A realistic 3D black-hole configuration has been first considered in Ref. [54], with also account for 3 body losses. In the present discussion we focus on the treatment of transverse excitations. Assuming that Bose condensation is total, but not disregarding the transverse degrees of freedom, one considers the dynamics of the system as described by a Gross-Pitaevskii equation for the classical field Ψ 0 ( r, t): where r = x e x + r ⊥ , r ⊥ denoting the transverse coordinate, V ⊥ ( r ⊥ ) = 1 2 mω 2 ⊥ r 2 ⊥ being the transverse potential and g = 4π 2 a/m. In (98) the normalization is chosen so that ρ 0 ( r, t) = |Ψ 0 ( r, t)| 2 is the density of particles. At equilibrium, in the so-called Thomas-Fermi limit [55], the Laplacian term in (98) can be omitted and the density has a cylindrical symmetry with Here µ is the chemical potential fixed by the normalization: d 2 r ⊥ ρ 0 (r ⊥ ) = n ; µ = 2 ω ⊥ √ a n [56,57]. The Thomas-Fermi approximation holds in the large density limit a n 1 [28,58]. In this limit, which has been denoted as "3D cigar" in Ref. [59], the classical field description is accurate, that is, the quantum fluctuations around Ψ 0 are small.
In the cylindrical geometry we consider here, the excitation spectrum has several branches corresponding to density fluctuations of the form δρ 0 ( r, t) = δρ(r ⊥ )e imθ e i(qx−ωt) , where θ is a polar angle in the transverse plane. For each branch the lowest state is obtained for m = 0 and q = 0 and its energy reads ω n = ω ⊥ 2n(n + 1), with n ∈ N. Taking into account possible longitudinal excitations one gets in the long wave-length limit [56,57] where c TF = µ/2 m is the sound velocity in the Thomas-Fermi limit (which has been measured by the MIT group [60]) and R ⊥ = 2 c TF /ω ⊥ is the transverse extension of the condensate (in the same limit). Eqs. (100) and (101) describe a lower mode with sonic-like dispersion relation and gaped transverse excited states which behave quadratically at low q. Note that the low q expansion displayed in Eq. (100) does not correspond to what is expected from the usual Bogoliubov dispersion relation (6). This is a hint that the longitudinal dynamics of the system is modified in the Thomas-Fermi limit. Of course, the hydrodynamic result (100) is limited to the region q R −1 ⊥ and cannot provide a reliable description of the whole excitation spectrum. But the departure from the usual Bogoliubov dispersion is confirmed by numerical solutions of Bogoliubov-de Gennes equations [61,62] which are valid for the whole range of wave vectors and for a range of densities larger than those based on the Thomas-Fermi approximation. These computations show that when increasing the linear density starting from a value n ∼ a −1 (i.e., when one goes deeper in the Thomas-Fermi regime) the dispersion relation develops a plateau in the region q ∼ 1/R ⊥ . This is interpreted as a tendency of the excitations to explore the radial parts of the condensate where the density is lower and where the local sound velocity accordingly decreases. This effect is not taken into account in the theoretical analysis presented in the main text where we work in a regime which has been denoted as "1D mean field" in Ref. [59]. At zero temperature this corresponds to the regime where the condition (97) is supplemented by n a 1 .
In this regime µ = 2 ω ⊥ an [26] which is much smaller that the energy 2 ω ⊥ of the first transverse excited state 6 , one can thus safely neglect transverse excitations and the transverse density profile is not of the type (99), but has rather a Gaussian shape. It is interesting to evaluate the actual range of parameters corresponding to the fulfillment of conditions (97) and (102), which is the regime of validity of our approach. For a transverse trap of frequency of 1 kHz, one gets for 23 Na (a ⊥ /a) 2 = 1.7 × 10 −5 , for 87 Rb (a ⊥ /a) 2 = 2.6 × 10 −4 and for He * (a ⊥ /a) 2 = 2.2 × 10 −5 [64]. Hence the domain of validity of the 1D mean field approximation used in the present work ranges over four orders of magnitudes in density.
In present time experiments, when the 1D mean field regime fails, this is mostly due to the fact that the the linear density is large, and in this case (102) may be violated. Then, the transverse density profile has the Thomas-Fermi shape (99). It is thus of interest to briefly and qualitatively discuss the features appearing in momentum correlators in acoustic black holes due to the transverse modes (101) 7 .
A first remark is in order here: the new transverse modes are typically not coupled to the modes studied in the present work. The reason for this is that the potential V (x) used to implement the sonic horizon does not couple modes with different transverse quantum numbers. Only small imperfections and nonlinear effect would induce such a coupling, and the results presented in the present work would remain almost unaffected. If the dispersion relation were limited to expression (101), i.e., were of the Klein-Gordon type, new outgoing modes would appear which would be populated by the time dependent formation process of the horizon, as in the gravitational context. In the present case however, the transverse dispersion relation (101) encompasses terms of higher order in q, and, as a result, new transverse incoming modes appear in the supersonic region, of the d2|in type, as illustrated in Fig. 9. As a result of the existence of these new incoming modes, the Hawking radiation process would occur also in the transverse sector even in the stationary context, and consequently new correlations lines appear which should add to the ones studied in the main text. However, in the regime (102) they should correspond to a very weak signal.

Conclusions
In this work we have investigated in detail the two-body momentum correlation in a quasi 1D BEC in the presence of a sonic horizon. Our modeling of the measurement process shows that the measurements have to be performed with some care: (1) the spatial windows selected for the Fourier analysis must be chosen carefully in order not to damp the expected signal (see Appendix A and also Ref. [34]) and (2) a separation of the upstream and downstream signals in the detection scheme favors the highlighting of quantum non separability (end of Sec. 3.1). Once the appropriate requirements are met, the normalized correlator (29) appears to be a robust quantity, making it possible to test quantum entanglement in a rich variety of situations, namely, in situ, or after artificially inducing a quench in the system or at the end of an adiabatic expansion after opening of the trap. Among the possible implementations of a sonic horizon we have studied, the largest quantum correlation signal, observed between the Hawking quantum and its partner, is realized in the so-called "waterfall configuration". In this configuration the Cauchy-Schwarz inequality is violated up to a temperature larger than the chemical potential and therefore should be experimentally testable in a finite temperature setting [the violation is still present at T = 1.5 × mc 2 u , see Fig. 6, upper left panel].

A Rigorous local Fourier transform
In this appendix we present the precise form of the Fourier transforms performed by using the window functions which, in the upstream region, are of the form (16). We first note that the Fourier transform of the window function is: where δ (1) is an approximation of the Dirac distribution (tending towards the δ function when σ u → ∞). One also has where δ (2) is another approximation of the Dirac δ-distribution defined by (104) and verifying δ 2π. So, instead of the approximate formula (18) one getŝ Formula (19) is modified in a similar way: In the following we present the results for the flat profile configuration. In this case K u = K d ≡ K 0 and we note k = K − K 0 , q = Q − K 0 . We first evaluate the one-body term N (K) which has contributions coming from both ψ † u (K)ψ u (K) and ψ † d (K)ψ d (K) . The double integral over ω and ω defining these terms is reduced to a single integral by means of the contractions (52). In this integral one can safely discard overlap terms such as Π u (k − q u|out (ω))Π u (k − q u|in (ω)) when σ u → ∞ since Π u (K) is proportional to δ (1) u (K). One thus gets terms generically of the form of the one resulting from the contraction of the first term of the integral in the r.h.s. of (105) with its hermitian conjugate, which reads: In the integral of the second term of (107) one has made the change of variable p = q u|out (ω) and all the ω-dependent terms have to be evaluated at ω u|out (p). In the last term of (107) one has used the fact that |Π u (k − p)| 2 is proportional to δ (2) u (k − p) and all the ω-dependent terms have to be evaluated at ω u|out (k). As can be checked for instance by comparison with the similar contribution to N (k) in Sec. 3.2, using the correct windowing for the Fourier transform, instead of the singular term δ(k − k) obtained with the schematic rules R1 and R2, one gets now a factor σ u Λ 2 u / √ 8π for the terms issued from the upstream window and a factor σ d Λ 2 d / √ 8π for the terms issued from the downstream window. As an illustration, the formulas equivalent to (37) and (38) (which, we recall, correspond to the adiabatic and zero temperature situation) read here and It now remains to evaluate the two-body function G 2 (K, Q). This involves four integrations over ω, two of which disappear when using the contraction rules (52). In the contributions to G 2 (K, Q) one has to distinguish the diagonal terms -i.e., intra-channel correlations -and the crossed ones -inter-channel. The evaluation of the diagonal terms is simpler, and we only state the results: Instead of the singular term δ 2 (k − q) (such as obtained for instance in the diagonal terms of (57) and (58)) one gets a term Λ 4 u (σ u /4 √ π)δ (1) u (k − q) for the contributions from the upstream windowing and a term Λ 4 for the contributions from the downstream windowing.
One now has all the tools for determining the effect of the windowing on the evaluation of the intra-channel correlation signals, of the type g 2 (K, K) u|out for instance: where N (K) u|out and G 2 (K, K) u|out are the u|out contributions to N (K) and to G 2 (K, K).
With the correct rules presented above, one gets g 2 (K, K) u|out = 2, as in the main text. The same result holds true for all the intra-channel correlation terms. We present the evaluation of the inter-channel terms in more detail, because it is less straightforward than the one of the diagonal terms, and also because it involves considerations relevant to the experimental detection scheme. Let us focus on the u|out-d2|out contribution for instance. As done above [Eq. (107)] we illustrate the general case by studying one of the many contributions to G 2 | u|out−d2|out . In the four field quantity (30), one has a product of four integrals of the type (105) and (106). For instance, one of the double contractions of terms issued from these integrals is The contractions are evaluated using the finite temperature rules (52), and the contribution of the term corresponding to (111) can be written as the products of two independent integrals, I and J: The two integrals have similar forms. Each appears with a prefactor which disappears when the product I × J is performed: we thus drop this prefactor in the following and denoteĨ and J the integrals where this prefactor is removed. We now focus on the evaluation ofĨ; after a change of variable p = q u|out (ω) it reads where A(ω) = (1 + n D2 )S d2,d2 S * u,d2 U * u|out W d2|out |∂ω u|out /∂p|. For presenting the results it is easier to work in a simple regime where the dispersion relations are dispersionless; this will be assumed in the remaining of this appendix, the general result being given in the final formula (126). In this case ∂ω u|out /∂k = V u − c u ≡ V u|out and ∂ω d2|out /∂k = V d + c d ≡ V d2|out and 8 one can write q d2|out (ω u|out (p)) = −γ p, where γ ≡ −V u|out /V d2|out > 0, cf. Fig. 10 (the notation γ is temporarily introduced to make the computations easier to follow). Then (114) reads I(k, q) = Λ u Λ d σ u σ d 2 ∞ 0 dp 2π A(ω u|out (p)) exp{T (p, k, q)} where T (p, k, q) = − σ 2 with P (k, q) = σ 2 u k + σ 2 d γq − 2i(X u + γX d ) ω k p q d2|out (ω u|out (p)) = −γ p ω d2|out (k) = V d2|out · k ω u|out (k) = V u|out · k .
It suffices to evaluate the integral (115) for σ u and σ d → ∞, which is the relevant limit as explained in the main text (cf. Sec. 2.4.1). In this case A(ω u|out (p)) is a weakly dependent function of p compared with the rapidly varying exponent, and the integral in (115) can be computed by means of the steepest descent method. This amounts to evaluate A(ω u|out (p)) for p = P (k, q) and to compute the remaining Gaussian integral. The result reads I(k, q) = Λ u Λ d A(ω u|out (P ))V d2|out δ u.d (V u|out k + V d2|out q) exp{−Z(k, q)} , where is again an approximation of the Dirac δ-function. The term exp{−Z} in (119) induces a damping which is not present in the schematic approach presented in the main text. This term can be removed if one imposes X u = −γX d , i.e., This relation has a simple physical interpretation: the time taken by an elementary excitation pertaining to the u|out channel to go from the horizon to the center (X u < 0) of the upstream detection zone has to be the same as the time taken by its partner (d2|out channel) to go from the horizon to the center (X d > 0) of the downstream detection zone. Note that this relation depends on the signal one is interested in (here the u|out − d2|out channel): For other channels (say the u|out − d1|out channel) the condition (121) will be modified and the centers of the window functions have to be shifted accordingly.
If the condition (121) is not fulfilled, the measured correlation will be damped compared to the perfect result [presented in the main text]. Note also that when this condition is fulfilled, since q = γ k from the δ u.d contribution in (119) one has P (k, q) = k and, in this equation, A(ω) is evaluated at ω u|out (k) as expected. Once condition (121) is realized, the product I × J =Ĩ ×J is found to be equal to The same contribution evaluated with the less rigorous approach presented in the main text yields to a very similar expression, where the Λ 2 u Λ 2 d prefactor is missing and where the term δ (3) u.d (V u|out k + V d2|out q) 2 is replaced by δ 2 (ω u|out (k) − ω d2|out (−q)). Finally, we consider the evaluation of the normalized inter-channel correlator When evaluating the fraction appearing in the r.h.s. of (123) along the line ω u|out (k) = ω d2|out (−q) one obtains a ratio identical to the one obtained in Eq. (60) of the main text, multiplied by a factor This term is equal to unity, as it should, only if i.e., if the width of the window functions Π u (x) and Π d (x) are in the same ratio (121) as their center.
We recall that we have used a simplified linear dispersion relation for deriving the relations (121) and (125). However, there is dispersion in the system; this means that these relations have to be adapted for each k and q along a specific correlation line: they should read in the u|out − d2|out case considered here X u ∂ω u|out /∂k = X d ∂ω d2|out /∂q , and σ u |∂ω u|out /∂k| = σ d ∂ω d2|out /∂q .
The same condition has been already derived by de Nova, Sols and Zapata in Ref. [34].

C Subsonic flow in the presence of a localized obstacle
In this appendix we consider the scattering of a stationary subsonic flow onto a localized external potential and we assume that the downstream flow is also subsonic. This is a special case of the situation considered in section 4. The configuration is illustrated in Fig. 11. We show in this case that -when the non-linearity coefficient g is x-independent -the farupstream velocity and density are equal to the far-downstream velocity and density: V u = V d and n u = n d .
n u n(x) U ext (x) x V u V d n d Figure 11: Sketch of the situation considered in the present appendix. The far upstream and far downstream asymptotic flows are both subsonic. The obstacle is represented by a localized external potential U ext (x).
Let's initially assume that the far upstream Mach number (M u = V u /c u ) and the far downstream one (M d = V d /c d ) are both less than unity, possibly different, and also that V u (n u ) may be different from V d (n d ). In our stationary setting, from the conservation of current and the definition of the speed of sound one gets where the last equality defines the quantity X.
The equality of the upstream and downstream chemical potentials reads 1 2 mV 2 u + gn u = 1 2 mV 2 d + gn d .