Numerical observation of Hawking radiation from acoustic black holes in atomic Bose-Einstein condensates

We report numerical evidence of Hawking emission of Bogoliubov phonons from a sonic horizon in a flowing one-dimensional atomic Bose-Einstein condensate. The presence of Hawking radiation is revealed from peculiar long-range patterns in the density-density correlation function of the gas. Quantitative agreement between our fully microscopic calculations and the prediction of analog models is obtained in the hydrodynamic limit. New features are predicted and the robustness of the Hawking signal against a finite temperature discussed.


Introduction
Back in 1974, Hawking [1,2] showed that black holes are not completely 'black' objects, but rather emit thermal radiation at a temperature inversely proportional to their mass. This amazing prediction, crucial to establishing the connection between black holes and thermodynamics [3], represents a genuine quantum effect in a gravitational context and is widely considered as a milestone of modern theoretical physics. Despite its conceptual importance, the weak intensity of Hawking radiation has so far prevented any direct experimental observation.
On the basis of a formal analogy between the propagation of waves in inhomogeneous and moving media and the propagation of fields on a curved space-time background, Unruh predicted in 1981 the occurrence of Hawking radiation in any system developing a horizon for some wavy perturbation [4]. In the following years, many systems have been proposed as candidates for actual experimental detection of this class of effects [5], e.g. superfluid liquid helium [6], atomic Bose-Einstein condensates (BECs) [7]- [9], degenerate Fermi gases [10], slow light in moving media [11], traveling refractive index interfaces in nonlinear optical media [12] and surface waves in water tanks [13].
In the present paper, we shall focus our attention on the specific case of Bogoliubov phonons propagating on top of a moving atomic BEC. Atomic BECs are among the cleanest systems where quantum physics can be investigated [14]- [16]: the temperature can in fact be made so low that the behavior of matter is dominated by the dual particle-wave nature of its constituents and the quantum dynamics is not masked by spurious thermal effects. Furthermore, in contrast to other quantum coherent condensed-matter systems such as superfluid liquid helium, quantitative theories able to describe the collective dynamics from a microscopic standpoint are available.
The possibility of creating in these systems black hole-like configurations that allow the study of the analog of Hawking radiation has been discussed in the last few years by several authors [7]. In a recent work [17] based on the gravitational analogy, we predicted that in this case a very characteristic pattern should appear in the correlation function for the density fluctuations of a condensate as a consequence of the Hawking effect. As demonstrated by a number of recent experiments [18]- [23], the measurement of density correlations is a powerful tool to extract information on the microscopic physics of atomic gases. Since the Hawking effect 3 consists of pairs of correlated phonons being emitted in opposite directions from the horizon, the quantum correlations within a pair propagate across the condensate and result in long-range density correlations between distant points on opposite sides of the horizon.
The present paper reports numerical simulations that nicely confirm our predictions of [17] and further support the promise of experimental detection of Hawking radiation from the density correlations rather than from the phonon flux. We also point out qualitative features that can be used to univocally distinguish the Hawking signal from fluctuations of different nature. Differently from most previous works on analog models, the present calculations are based on the application of microscopic many-body techniques to an experimentally realistic system. Since our numerical calculation never relies on the gravitational analogy, it represents independent evidence of the existence of the Hawking effect in a realistic condensed-matter system. The quantitative intensity of the Hawking signal is in agreement with the gravitational analogy.
The structure of the paper is as follows. In section 2, we introduce the physical system under investigation and in section 3, we summarize the theoretical method used for the calculations. The numerical evidence of Hawking radiation is presented in section 4 and quantitatively analyzed in section 5, where an extensive comparison with analytical results based on the gravitational analogy is made. The effect of a finite initial temperature on the Hawking effect is discussed in the following section 6. A brief discussion of the actual observability of the predicted effect with state-of-the-art systems and detection schemes is given in section 7. Conclusions are finally drawn in section 8.

An acoustic black hole in a flowing condensate
The physical system that we consider is sketched in figure 1(a): an elongated atomic BEC which is steadily flowing at a speed v 0 along an atomic waveguide. The transverse confinement is assumed to be tight enough for the transverse degrees of freedom to be frozen [16] and the system dynamics to be accurately described by a one-dimensional model based on the following second-quantized Hamiltonian, Here,ˆ (x) andˆ † (x) are atomic field operators satisfying Bose commutation rules , m is the atomic mass, V (x) the external potential and g(x) is the atom-atom interaction constant [16]. If the characteristic size ⊥ of the transverse wavefunction φ(r ⊥ ) is much longer than the atom-atom scattering length a 0 , the interaction constant has the form [24,25]: which simplifies to g = 2hω ⊥ a 0 in the most relevant case of a cylindrically symmetric and harmonic transverse trapping potential of frequency ω ⊥ if interactions are weak enough not to distort its Gaussian ground state wavefunction. In this limit, also the contribution to the external potential V (x) due to the zero-point energy of the transverse ground state has a simple form 2 ×hω ⊥ /2.  Initially, the condensate has a spatially uniform density n. The external potential and the (repulsive) atom-atom interaction constant are also uniform and equal to, respectively, V (x) = V 1 and g(x) = g 1 > 0. Around t = t 0 , a step-like spatial modulation of characteristic thickness σ x is applied to both the potential and the interaction constant by suitably modifying the transverse confinement potential and/or the atom-atom scattering length via an external magnetic field tuned in the vicinity of the so-called Feshbach resonance [16]: within a short time σ t , V and g in the downstream x > 0 region are brought to V 2 and g 2 , whereas their values in the x < 0 region are kept equal to the initial ones V 1 and g 1 .
In order to suppress competing processes such as back-scattering of condensate atoms and soliton shedding from the potential step [25], the external potential V is chosen to exactly compensate the spatial jump in the Hartree interaction energy µ 1,2 = g 1,2 n, i.e. V 2 + µ 2 = V 1 + µ 1 . 6 In this way, the plane wave that describes the condensate evolution at the mean-field level [16].
In what follows, we will focus our attention on the black hole (or dumb hole) configuration, where the speed of sound c 1,2 = µ 1,2 /m in the different regions satisfies the chain inequality c 1 > v 0 > c 2 : a black hole-type sonic horizon separates a region of subsonic c 1 > v 0 flow outside the black hole from a supersonic v 0 > c 2 one inside the black hole. As one can see from the Bogoliubov dispersion plotted in figure 1(c), long-wavelength phonon excitations propagating in the supersonic region are dragged away by the moving condensate and are then unable to propagate back to the horizon. Only higher k, single-particle excitations outside the hydrodynamical window can emerge.
In agreement with previous work [26], dynamical stability of this black hole configuration has been numerically verified. A more detailed discussion on stability issues is postponed to a forthcoming publication [27].

The Wigner method
The dynamics of fluctuations around the plane-wave mean field solution (3) during and after the formation of the horizon can be numerically studied by means of the so-called truncated Wigner method for the interacting Bose field [28,29]. Application of this technique to calculate the time evolution of generic observables of an interacting Bose gas has been extensively developed and validated in [30]. For dilute gases such that nξ d 1 (the healing length is defined as usual as ξ = h 2 /mµ and d is the dimensionality of the system, d = 1 in our case), this technique is equivalent to the time-dependent Bogoliubov approach and has the advantage of being able to follow the system for longer times when the backaction of quantum fluctuations on the condensate starts to be important. A first application of this method to atomic BEC-based analog models is reported in [9].
At t = 0 well before the formation of the black-hole horizon, the condensate is assumed to be uniform and at thermal equilibrium in the moving frame at v 0 at a temperature T 0 . Within the Wigner framework, this corresponds to taking a random initial wavefunction ψ 0 (x) of the form: As both the external potential and the interaction constant are spatially uniform V (x) = V 1 and g(x) = g 1 , the Bogoliubov modes have the standard plane-wave form. The Bogoliubov coefficients u k , v k are defined in terms of the kinetic and Bogoliubov energies E k =h 2 k 2 /2m The mode amplitudes α k are independent, zero-mean Gaussian random variables such that α k = α 2 k = 0. The variance |α k | 2 = [2 tanh( k /2 k B T 0 )] −1 tends to a finite value 1/2 in the T 0 → 0 limit, so to include zero-point fluctuations. For each realization of the random amplitudes α k s, the condensate density n 0 has to be suitably renormalized so as to account for the condensate depletion [30]. 6 The random classical wavefunction ψ(x, t) is then propagated in time according to the standard GPE (4) starting from its initial value ψ(x, t = 0) = ψ 0 (x) and including the full timeand space-dependence of V (x, t) and g(x, t). Periodic boundary conditions are assumed for the spatial variable x, but an absorbing region far from the horizon is required to prevent the onset of spurious dynamical instabilities due to excitations circulating around the integration box [26,31]. Expectation values of symmetrically ordered field observables at any later time t are finally obtained as the corresponding averages of the classical field ψ(x, t) over the random initial condition ψ 0 (x). As usual, trivial one-time commutators have to be subtracted out if normally ordered quantities are required.

Numerical evidence of Hawking radiation
Our numerical study of the Hawking radiation from the acoustic black hole is based on a vast campaign of Wigner simulations of the condensate evolution after the formation of the blackhole horizon. Inspired by our recent work [17], specific attention has been paid to the correlation function of density fluctuations, i.e. the normalized, normal-ordered density-density correlation function: Three successive snapshots of G (2)  Apart from feature (i), which is the usual antibunching due to the repulsive atom-atom interactions [32], these observations illustrate a variety of effects of quantum field theory in a spatially and temporally varying background [33].
Feature (ii) is a transient effect that originates in the bulk of the internal x > 0 region and has no relation to the presence of a horizon. An identical fringe pattern (yet extending to the whole system, rather than to the x > 0 region only) is in fact obtained if a spatially uniform system is considered whose interaction constant is varied in time from g 1 to g 2 with the  same functional law ( figure 4(b)). Its physical interpretation is the following: as a consequence of the time modulation of the interaction constant g, correlated pairs of Bogoliubov phonons are generated during the short modulation time by a phonon analog of the dynamical Casimir effect [34]- [37], a quantum process which shares many analogies with the parametric emission of Faraday waves in classical fluids [38]. As the emission process takes place in a simultaneous and coherent way at all spatial positions, the two phonons are emitted with opposite wavevectors and then propagate in opposite directions. Their quantum correlations reflect into a pattern that only depends on the relative coordinate x − x and consists of fringes parallel to the x = x line, which move away in time at approximately twice the speed of sound. As usual, the longer the ramp time σ t , the weaker this dynamical Casimir signal which eventually disappears in the limit of a very long σ t .
On the other hand, features (iii) and (iv) do not depend on σ t , but only on the eventual presence of a horizon at long times. In particular, they completely disappear if the flow remains everywhere sub-sonic v 0 < c 1,2 (figure 4(a)) or, a fortiori, if a spatially uniform system is considered ( figure 4(b)). This fact, together with their shape in the (x, x ) plane and their persistence for indefinite times after the horizon formation suggests a strict link with the Hawking effect.
As anticipated in [17], the quantum correlations within a pair of Hawking phonons emitted in, respectively, the inward and outward directions translate into a correlation between the density fluctuations at distant points located on opposite sides from the horizon. Once the horizon is formed, correlated pairs of Hawking phonons are continuously emitted at all times t > t 0 on the α and β phonon branches of figures 1(b) and (c). These phonons then propagate from the horizon in, respectively, the outward and inward directions at speeds v 0 − c 1 < 0 and v 0 − c 2 > 0. A generic time τ after their emission they are located at x = (v 0 − c 1 ) τ < 0 and Feature (iv) originates from the (partial) elastic backscattering of the α Hawking particle onto branch γ : both γ and β phonons eventually propagate in the inward direction at speeds v 0 ± c 2 . Again, the analytically predicted slope (v 0 + c 2 )/(v 0 − c 2 ) (green dashed lines) agrees well with the axis of the numerically observed tongue (iv).

Quantitative analysis
The identification of features (iii) and (iv) as signatures of Hawking radiation is further confirmed by a quantitative comparison of our numerical data with the predictions [17] of the gravitational analogy. The peak value G (2) (peak) and the transverse full-width at half-maximum (FWHM) x FWHM of the tongue (iii) have been extracted from the cut of G (2) (x, x ) taken along a straight line x = x + x cut well outside the antibunching dip (i) (indicated as a black line in figures 2(b) and (c)). Examples of such cuts are shown in figure 3(b). Once the tip of the tongue has crossed the cut line, neither the peak value nor the FWHM depend any longer on the position x cut of the cut line or on the time t of the observation: as expected, the Hawking emission is in fact stationary in both space and time.
In figure 5(a), the peak intensity of the correlation signal (iii) is plotted as a function of the inverse of the thickness σ x of the horizon region, a quantity that in the gravitational analogy is proportional to the so-called surface gravity. As expected, the agreement with the gravitational prediction [17]: is quantitatively excellent in the hydrodynamic limit σ x /ξ 1,2 1, where the physics of the many-body system is dominated by the hydrodynamic modes that are included in the gravitational analogy. As usual, the parameter κ proportional to the surface gravity of the horizon is defined as In the plots, we have shown the renormalized quantity (nξ 1 ) × [G (2) (x, x ) − 1] that is universal in the dilute gas limit nξ 1 1, where our Wigner approach is exact. The actual intensity of the Hawking signal is then inversely proportional to the dilution parameter nξ 1 . Although a more sophisticated theoretical approach may be required to obtain quantitative predictions in the strongly interacting n ξ 1 1 case, the qualitative features of Hawking physics are not expected to change if a relatively strongly interacting system is used in order to maximize the intensity of the Hawking signal. The inverse FWHM of the tongue (iii) in the transverse direction is plotted in figure 5(b) as a function of the inverse thickness of the horizon region. Again, the agreement of the numerical result with the gravitational prediction of [17] is very good as long as σ x /ξ 1 1. This observation appears even more significant if one remembers that in this framework the FWHM is a good measure of the surface gravity κ, which then fully determines the Hawking temperature T H =hκ/(2πk B ) [33].

Effect of a finite initial temperature
To verify the robustness of our observations with respect to thermal effects, we have performed a series of numerical calculations starting from a finite initial temperature T 0 > 0. The results are summarized in figure 6.
As happens for thermal light at semi-reflecting mirrors [28], the partial reflection of thermal phonons on the horizon provides a further mechanism for the appearance of non-trivial density correlation between the inner x > 0 and the outer x < 0 regions which may mask the Hawking signal. This concern is ruled out by the numerical results shown in figure 6(a). Even for an initial temperature k B T /µ 1 = 0.1 higher than the Hawking temperature k B T H /µ 1 0.05 expected from the gravitational analogy [5], the (iii) and (iv) Hawking tongues remain perfectly visible and are even a bit strengthened by stimulation effects [39]. The new feature (v) that originates from thermal effects is located between the tongues (iii) and (iv) and is well distinct from them. Its attribution to partial reflection of thermal phonons is confirmed by the value of its slope, which is in very good agreement with the analytical prediction (v 0 + c 2 )/(v 0 − c 1 ) indicated as a dashed black line in figures 6(a) and (b).
As the slope of the thermal tongue is completely different from the one (v 0 − c 2 )/(v 0 − c 1 ) of the Hawking tongue, the Hawking signal can be easily isolated in a correlation image from spurious thermal and atom loss [40] effects also at temperatures higher than T H . This is even more remarkable as in this regime it appears difficult to separate the Hawking phonons from the thermal ones by looking just at the total phonon flux.
As a final check, we have performed a numerical calculation for the case of a flow which remains everywhere sub-sonic. The result is shown in figure 6(b): as expected, the Hawking tongues disappear, whereas the dynamical Casimir fringes (ii) and the thermal tongue (v) persist without being dramatically affected.

Discussion of some experimental issues
We conclude with a brief discussion on the main issues that may arise in an actual experiment that aims to observe the Hawking radiation from the density correlations. As the magnitude of the Hawking signal is quite low, of the order of (nξ 1 ) × G (2) (peak) ≈ 5 × 10 −3 , a critical point in an actual experiment will be the signal-to-noise ratio in extracting the correlation signal due to the Hawking effect from experimental noise. While many sources of noise can be avoided by carefully engineering the experimental setup, shot noise is an intrinsic consequence of the discrete (i.e. quantum) nature of atoms and therefore cannot be eliminated.
A precise measurement then requires large statistics to be collected, either by repeating the experiment many times or by creating a large number of (almost) identical BEC samples 13 on which to simultaneously perform the measurement. Most probably, both these strategies will face the difficulty of keeping the parameters of the experimental setup constant in space and/or time.
Another possibility is to use a continuous-wave atom laser beam propagating along an atomic fiber [41] and then perform a time average of the correlation signal over long times. Several different outcoupling mechanisms have been demonstrated, which allow for a wide tuning of both the density and the flow speed of the atom laser beam. Furthermore, active stabilization techniques [42] can be used to keep the parameters of the output beam stable in time during the experiment.
In this perspective, the key experimental problem that remains so far unsolved is the one of continuously reloading the mother BEC while emitting the atom laser beam. Existing devices [41,43] are in fact limited by the finite number of atoms present in the mother BEC, presently of the order of 10 5 , but attempts to continuously reload the condensate during the atom laser operation have been recently reported [44].
As the Hawking signal is inversely proportional to the diluteness parameter nξ 1 , a relatively strongly interacting system may be favorable to maximize the signal-to-noise ratio. As our predictions are basically a consequence of quantum hydrodynamics, we expect that their validity extends well beyond the case of weakly interacting BECs that has been considered by our microscopic model. Although no complete theory is available for a strongly interacting gas, the theory of Luttinger liquids [45] appears to indicate that larger correlation signals can be observed in stronger interacting systems, where the diluteness parameter is smaller.
Another crucial issue is to use an atomic detection scheme able to efficiently measure the density correlations from which to extract the Hawking signal: several techniques have been developed during the last few years to experimentally characterize local density fluctuations in ultracold atomic gases based on noise in absorption images [18,19], on microchannel plate detectors [20], or even on high-finesse optical cavities [23]. Thanks to its non-destructive nature, this last technique to detect single atoms in an atom laser beam appears to be very promising for the detection of Hawking radiation from density correlations.

Conclusions
In conclusion, we have provided clear numerical evidence of the presence of Hawking radiation in a flowing atomic BEC in an acoustic black-hole configuration. The signature of this effect was identified in the density-density correlation function for distant points located on opposite sides of the horizon. The intensity of the observed signal is in quantitative agreement with the predictions of the gravitational analogy in the hydrodynamic regime.
A crucial novelty of this work consists of the fact that all numerical results have been obtained starting from a microscopic description of the system dynamics using standard methods of many-body theory, and the gravitational analogy has provided a physical insight to interpret the numerical observations. In this way, our numerical observations can be considered as an independent proof of the existence of Hawking radiation and rule out some frequent concerns [46] on the role of short wavelength, 'trans-Planckian' physics on the Hawking emission.
We have shown that the stationary emission of phonons from the acoustic horizon persists even outside the hydrodynamical regime when the gravitational analogy can no longer be applied, and shares many of the typical properties of Hawking radiation.
The characteristic pattern of the Hawking signal in the density correlation function allows an easy identification from spurious effects of, e.g., thermal origin even at temperatures significantly higher than the Hawking temperature. This supports our proposal of using density correlations as a powerful experimental tool to detect Hawking radiation from acoustic black holes in atomic BECs in the near future.