Acoustic analogue of Hawking radiation in quantized circular superflows of Bose-Einstein condensates

We propose emulation of Hawking radiation by means of acoustic excitations propagating on top of a persistent current in an atomic Bose-Einstein condensate loaded in an annular confining potential. The setting admits realization of sonic black and white event horizons. It is found that density-density correlations, representing the acoustic analogue of the Hawking radiation, are strongly affected by the perimeter of the ring-shaped configuration and number of discrete acoustic modes admitted by it. Remarkably, there is a minimum radius of the ring which admits the emulation of the Hawking radiation. We also discuss a possible similarity of properties of the matter-wave sonic black holes to the known puzzle of the stability of Planck-scale primordial black holes in quantum gravity.


I. INTRODUCTION
Black holes (BHs) are among the most fascinating objects in the Universe, the study of which may help to illuminate an intimate relation between gravity and the quantum theory. Pioneering works [1][2][3][4] of Bekenstein and Hawking had allowed to elucidate remarkable quantum properties of BHs. It was discovered that BH is not a stationary everlasting object, as it emits black-body radiation. The effective temperature of the Hawking radiation is extremely low for known astrophysical BHs, which makes it practically impossible to observe the emission effect, therefore attention has turned to settings which admit emulation of this phenomenology in other physical settings.
Acoustic BHs were first introduced by Unruh [5,6], who had demonstrated precise formal equivalence between the behavior of sound waves in a fluid flow and that of a scalar field in curved space-time. Thus, it is possible to create analogue BH in superfluids where the transition from subsonic to supersonic flow plays the role of the event horizon. Following that work, many experimental realizations were proposed for demonstrating acoustic horizons. Behaviors similar to the analog Hawking radiation have been experimentally and theoretically explored in trapped ions [7], optical fibers [8][9][10][11], electromagnetic waveguides [12], water tanks [13,14], ultracold fermions [15], exciton-polariton condensates [16] and superfluid 3 He [17]. Tremendous progress in physics of ultracold atomic gases [18] makes atomic Bose-Einstein condensates (BECs) very suitable testbeds for sonic BHs [19,20]. In this connections, recent experimental realization [21][22][23][24] of analogs of the Hawking radiation in atomic BEC is an important milestones in simulations of quantum properties of BHs.
Previous investigations [43,44] of the event horizon in BEC with toroidal geometry relied upon the use of the supersonic flow driven by a variable condensate density distribution. On the contrary, we here use the uniform condensate density, while the event horizon is created by spatiotemporal modulation of the coupling constant which is responsible for interatomic interactions in the mean-field approximation [18,19]. The coexistence of BH and white-hole event horizons in the toroidal geometry, gives rise to self-amplifying Hawking radiation and formation of background ripples [45][46][47][48][49]. Extensive investigation of related instabilities has been carried out in the 1D system with absorbing boundary conditions [50]. Here we focus on the Hawking radiation per se, rather than the self-amplification effect. As is known [6,23,51], the analog Hawking temperature for the effectively onedimensional flow is where v(x) is the local velocity of the condensate, c(x) is the local speed of sound, and the derivatives are evaluated at the position of the horizon, x = x h . To minimize the impact of the white hole horizon, we introduce a steep gradient of c(x) near the BH horizon, and a smoothed gradient near the white hole event horizon, as shown in Fig. 1.
The main objective of the present work is impact of the quantization of the superflow, imposed by the periodic boundary conditions, on the analog Hawking radiation in the toroidal BEC. The analysis reveals the following noteworthy effects: (i) Density-density correlations, representing the acoustic analog of the Hawking radiation, are strongly affected by the perimeter of the ring and the number of discrete acoustic modes available in the ringshaped geometry. (ii) Hawking radiation vanishes if the radius of the ring falls below some critical value.
FIG. 1. A schematic of the ring-shaped Bose-Einstein condensate with asymmetric acoustic event horizons. Shown are the blue isosurface of the condensate density n, and the superflowvelocity profile v (the dark blue line). Note that both n and v are spatially uniform. Azimuthally-modulated interaction strength g(x, t) is responsible for the variation of the local sound speed, c(x) (the solid black curve) along the ring. The condensate density remains constant both in the subsonic (1: v < c) and supersonic (2: v > c) regions. The arrow shows the direction of the persistent current. Acoustic black hole horizon (BHH) and white-hole event horizon (WHH) appear when the local sound speed is equal to the superflow velocity.
The rest of the paper is organized as follows. The model is introduced in Sec. II. In Sec. III we discuss initial fluctuations and the corresponding dispersion relation. In Sec. IV we present results for density-density correlations and the acoustic analog of the Hawking radiation near the BH horizon. The paper is concluded by Sec. V, where we also discuss possible similarities of properties of the analog BHs to some fundamental problems, such as stability of primordial quantum BHs on the Planck's scale.

II. THE MODEL
We numerically analyze the acoustic analog of the Hawking radiation in the framework of the onedimensional Gross-Pitaevskii equation (GPE) where ψ(x + L, t) = ψ(x, t) is the wave function of the ring-shaped condensate of length L and radius R, while V (x, t) and g(x, t) are the effective potential and 1D coupling constant, respectively. One of efficient methods for detecting the analog Hawking radiation is the analysis of density-density correlations [52][53][54]. We use here the mean-field dynamical simulations together with sampling of the quantum noise similar to the well-known truncated Wigner approximation (TWA) [55][56][57][58][59]. In our analysis of the quantum fluid dynamics, we consider the evolution of stochastic trajectories satisfying the 1D GPE (2). Similar to the standard TWA, we add extra noise to the initial wave function, thus taking it as follows: where k 0 = M v/h is determined by the flow velocity v, β k , β * k are coefficients with random values, the summation is performed over wavenumbers k, n 0 is the spatially uniform density, and u 0k , w 0k are defined as where It is relevant to briefly discuss the derivation of GPE (2) and the meaning of its parameters, as well as the meaning of the noise in Eq. (3), in more detail. As said above, we consider the system schematically displayed in Fig. 1. The azimuthally-modulated interaction strength, g(x, t), where x = Rϕ and ϕ is the azimuthal angle around the symmetry axis of the ring, is responsible for the variation of the local sound speed, c(x), along the ring. The starting point of the derivation is the 3D mean-field GPE: where M is the atomic mass, V e = V trap + V h is the external potential that consists of two terms, viz., the transverse confinement and the potential that creates the horizon, γ = 4πh 2 a s /M is the coupling constant, a s being the s-wave scattering length, which is supposed to be a function of coordinate and time.The 3D wave function Ψ(r, t) is related to the number of atoms, N , by the normalization condition, N = |Ψ| 2 d 3 r. We use N = 10 5 of 87 Rb atoms in our calculations. We consider a thin ring of radius R filled by dilute BEC, with the transverse degrees of freedom frozen by the tight confinement. Assuming the usual factorized approximation [29,[60][61][62], we set where φ(r ⊥ ) is the function of transverse coordinates (see, e.g. [63]), and the 1D density is normalized to the total number of atoms: Then averaging Eq. (5) in the transverse plane leads to the 1D GPE (2) with the 1D coupling constant defined as The state of the condensate in the ring with the uniform density is provided by the plane-wave solution [19] of Eq. (2): (9) Potential V (x, t) and nonlinearity coefficient g(x, t) must be mutually matched so as to admit the existence of the state with uniform density n 0 .
Since the wave function must be periodic, ψ(x) = ψ(x + L), possible values of k 0 , k in Eq. (3) are restricted by the length of the ring, k 0 , k = 2πm/L, m ∈ Z. This is related to the obvious quantization of the velocity circulation: where integer m is the winding number or topological charge. If we solve Eq. (2) on a grid of N p points with spacing h, then k 0 = 2πm/(N p h), where integer m coincides with the winding number.
A desirable spatial distribution of the local sound speeds in regions 1 and 2 ( Fig. 1) can be produced by making the coupling constant, g in Eq. (5), spatially and temporally inhomogeneous. This can be implemented by either tuning the scattering length a s by means of the Feshbach resonance [18,19,[64][65][66], or longitudinally varying the transverse confinement [67]. The sound speed in regions 1 and 2 is c 1,2 = g 1,2 n/M , where µ 1,2 = g 1,2 n are the effective local chemical potentials. For the plane wave (9) with uniform density n 0 to be a solution at all values of time, we need to adjust the axial potential and coupling constant in the two regions as follows: This condition is satisfied all over the ring. Finally, to cast the 1D GPE in the normalized form of Eq. (2), we apply rescaling, Boundary conditions for Eq. (2) being periodic, we used the split-step fast Fourier transform method [68] for numerical simulations. The transition from g 1 , V 1 to g 2 , V 2 in region σ BHH is smooth, provided by the potential taken, locally, as where In our simulation we used x BHH = −50, t 0 = 5, and change of V (x, t) and g(x, t) occurs in the course of time σ t . We used Heaviside function θ(t − t 0 ) to assure that V and g are uniform at t = 0. However, as we verified, our results do not change qualitatively if potential (12) is used without θ(t − t 0 ) providing σ t ≪ t 0 . Therefore, initial value of V (x, t) is V 1 . The initial condition taken as the uniform condensate with constant V, g makes it easier to add the long-wavelength noise to the entire system or some part of it, if necessary. Using this method with change of V (x, t) and g(x, t) in time and running hundreds of simulations, the initial noise added to the uniform condensate inevitably excites all possible longwavelength modes admitted by inhomogeneous g and V .
We have checked that, for σ t ∈ [0.1, 5], variation of this parameter does not change the dynamics qualitatively, therefore we here present the results for σ t = 0.5.
Handling the WHH, which always appears in the ring, is more challenging, therefore we performed simulations for two different corresponding transition functions, the first one for x > 0 chosen as that provides mirror symmetry for the potential as long One can see that this type of the transition function cannot allow an arbitrary value of the slope at the transition point, as it requires an extremely large length of the ring for large σ W HH and small |∂c/∂x| ∼ (σ WHH ) −1 . To deal with this case, we used, instead of transition function (14), a doublestep one, designed as a set of two similar potentials, see Figs. 3 and 4(b). These two functions gradually carry over into each other, allowing one to manipulate a sufficiently smooth slope near WHH. Parameters of the horizon in the dimensionless units are similar to those in the model used in Ref. [19] and close to the experimental parameters of Ref. [22] for the speeds: c 1 = 1 (1 mm/s), v = 0.74, c 2 = 0.5, σ x = 0.5, σ t = 0.5, with the corresponding dispersion relations for regions far from the horizon shown in Fig. 2. Numerical values of g 1,2 can be estimated from relation g = M c 2 /n and V 1 = 1.6µ 1 . The initial ring has the length of L = 233 µm, τ 1 = 0.73 ms, ξ 1 = 0.73 µm, c 1 = 1 mm/s, which is tantamount to L = 320, τ 1 = 1, ξ 1 = 1, and c 1 = 1 in the dimensionless units.

III. QUANTUM FLUCTUATIONS AND THE DISPERSION RELATION
In real BEC, the presence of quantum and thermal fluctuations, which are not taken into account by GPE, is inevitable. It is possible to include these fluctuations by adding Gaussian noise to the initial wave function according to TWA. The evolution of the initial fluctuations is governed by the Bogoliubov theory [18], and expectation values of symmetrically ordered observables can be obtained by taking the stochastic average over the ensemble of evolved wave functions.
The initial state of the unperturbed system is uniform, and the potential step is switched on at time t 0 > 0, therefore in our simulations we used the stochastic initial wave function (3), where the summation is performed up to some maximum wavenumber k. In the following section, we discuss an appropriate choice of this limit value in detail, and consider its impact on the analogue Hawking radiation.
It is straightforward to derive the respective dispersion relation for the acoustic waves: where wavenumbers k should satisfy the same periodic condition as wave function (9). To briefly explain the derivation of this relation and move forward, we consider the wave function in a general form, [cf. the initial condition given by Eq. (3)], where the perturbations are assumed to be small in comparison to √ n 0 , hence this expression can be rewritten as ψ( After inserting this ansatz in the GPE and linearizing with respects to ψ ′ (x, t), one obtains the set of equations: Since u k (x) = u 0k e i(kx+k0x) / √ L, w k (x) = w 0k e i(kx−k0x) / √ L and ω 0 = k 2 0 /2 + V + gn, we obtain: It is easy to see that dispersion relation (15) follows from the consistency condition of system (18), setting n = n 0 . Coefficients u 0k , w 0k (k = 0) are solutions of Eq. (18) such that the corresponding modes u k (x), w k (x) satisfy the normalization condition, L 0 |u k (x)| 2 − |w k (x)| 2 dx = 1. Therefore, expressions for the Bogoliubov coefficients become as mentioned above in (4).
To keep the number of atoms constant, we adjusted density n 0 in each simulation. The expression that defines the number of excited atoms in the uniform condensate was found in ref. [57]. At zero temperature, it can be written as Therefore, the number of atoms remaining in the ground states is N c = N − N ′ s , and n 0 = (N c + 1/2)/L. Such expressions follow from the relation between the average over ensemble and quantum average.According to Refs. [56,57], amplitudes β k , β * k of the random perturbation are distributed as follows: where σ ǫ = (4 tanh (ǫ k /2k B T 0 )) −1/2 , T 0 being the temperature of the condensate. If we split the amplitudes in real and imaginary parts, β k = β xk + iβ yk , and, accordingly, replace dβ k dβ * k → dβ xk dβ yk , then the integration of the probability density immediately gives 1, and we have two real quantities, β xk and β yk , with the same type of the distribution. To create such a distribution, we used standard randomizing functions: Clear evidence of the analogue Hawking radiation was obtained, both in the zero-temperature limit and for T 0 > 0, in Refs. [19,20]. Here, we use the distribution for T 0 = 0 and focus on this case.

IV. ANALOG HAWKING RADIATION NEAR BHH (THE BLACK HOLE HORIZON)
To identify the generation of the Hawking radiation from the acoustic BH, we used the normalized densitydensity correlation function: where averaging is performed over an ensemble of 100 GPE simulations. Increase of the number of simulations does not tangibly affect the results. All simulations started with the input in the form of the uniform condensate with quantum noise added to it, as per Eq. (3). An essential role in our investigation plays the number of modes which are used in TWA. This number was chosen to satisfy some natural restrictions. TWA was proven to be correct for dilute Bose condensates if the number of modes obeys the constraint N > N modes /2, as was shown in Ref. [56]. To provide the presence of P -modes and thermality of the outgoing flux, frequencies should obey condition ω < ω max , see Fig. 2, as was found in Ref. [69]. For a sufficiently low wavenumber the latter condition also implies that all modes satisfy regime of the linear dispersion relation, thus maintaining a clear analogy with the Hawking effect near real BH. It is worth to note that the initial noise in Eq. (3) involves summation over negative and positive wavenumbers k, but, in our investigation, only positive wavevectors are relevant, as they correspond to atoms moving towards BHH. Moreover, including modes with the negative direction of motion did not make a visible change in the obtained correlations, therefore we include only lower positive-k modes in our initial condition. First, we address the symmetric potential with initial parameters given above (L = 320, x BHH = −50, x WHH = 50, σ x = σ WHH = 0.5). Accordingly, in condition ω < ω max we have to put ω max ≈ 0.1. Thus, we are able to generate three modes with the initial frequency below ω max . The obtained results for the correlation function are illustrated in Fig. 3. Colored lines in the figure correspond to the expected correlations between different particle pairs from the dispersion relation in Fig. 2. The slope of the colored segment relative to the x-axis is determined by expressions for the sound speed of the Bogoliubov excitations: where θ y , θ g , and θ r refer to the yellow, green and red lines, respectively, in Fig. 3. The length of each segment reflects the expected length of the correlation tongues at the corresponding times. By tongues we mean the lineshaped correlation regions that have one end located near BHH and other end growing in time in some certain direction.It is seen that, in spite of the quantization of k, all the correlation tongues are visible for both potentials and they agree well with the predictions based on the dispersion relations. In Fig. 3 it is seen that the correlation pattern alters at later times [in Fig. 3(a) at t ≃ 150) because of a checkerboard pattern appears as a result of the existence of a cavity between the hori-zons [20,45]. To resolve this problem, we decreased the value of |∂c/∂x| ∼ (σ WHH ) −1 , which is responsible for mixing of modes at the horizon. Applying our doublestep potential, we could make |∂c/∂x| and thus avoid the destructive impact of such effects, by reducing the checkerboard correlations to just two crosses created by the scattering on each step. While we observe the correlation in both patterns for the current length, the presence of such decaying |∂c/∂x| is inevitable for smaller scales that are considered below. The next step is to apply the initial noise, which contains modes that are located above ω max in the dispersion relation. The corresponding correlations, produced by 10 lowest modes, are shown in Fig. 4. As above, we see the same correlations for steep and smoothed slopes at WHH, with a qualitative difference in the checkerboard correlations. Moreover, only HR − in correlation tongues (i.e., ones between excitation modes of the HR and in types) are present, and there is no evidence of P − in or HR − P tongues. Besides, obtained correlation tongues are expanding a bit faster than expected, which may be a direct consequence of being beyond the linear dipersion regime.
We have performed similar simulations for the eight times enlarged region and distance between the horizons (L = 2560 or ≈ 1.87 mm, in physical units), keeping the same parameters of both horizons as above. The larger area makes it possible to admits more initial modes that lie below ω max . We have thus performed the simulations for 23 lowest modes belonging to the linear dis-persion regime. In Fig. 5 one can see the entire range of the expected correlations, and the anticipated position of the tongues very well coinciding with the highlighted lines in the correlation pattern. These correlations feature relatively high intensity, which is about 4 × 10 −4 for the HR − P tongues. This value is still ≃ 10 times smaller than that reported in Ref. [19], which is explained by our choice of the diluteness parameter, ξ 1 n = 1/g 1 = 312.5 ≫ 1, that defines the intensity of the correlation signal. No drastic difference is seen between the correlations near BHH for both potential slopes at WHH. This fact is explained by a finite speed of excitations and relatively large distance between the two horizons. The results obtained for the initial noise containing modes over ω max in the large region (L = 2560) share main properties with two previous simulations: there is no difference between correlation patterns near BHH for FIG. 4. The same as in Fig. 3 (without the bottom panels) for L = 320 but for the inputs with ten initial modes. Note that checkerboard pattern vanishes for the smoothed white horizon as in Fig. 3.
two different slopes at WHH, at least in the course of the observation times, and the correlation pattern, as one may expect, features solely HR − in correlation tongues. The same happens with modes taken above ω max .
The inference is that the possibility to observe the analog Hawking correlations (for fixed parameters of the flow and horizon) depends on the number of modes that have their frequencies below and above ω max , as predicted by dispersion relation (15) for the uniform density. Furthermore, the effect is strongly affected by the size of the ring. The first aspect reflects the fact that, to produce the Hawking effect, one needs to have particles with the frequency from the negative-norm branch of the dispersion relation in the laboratory reference frame. This criterion was also discovered analytically for the elongated one-dimensional condensate with one horizon [69], and, remarkably, it is true for our system as well. The second aspect implies that, by enlarging the length of the ring, one increases the number of allowed modes, which, in the limit of the ring with an infinite radius leads to the continuous frequency spectrum and infinite number of modes below ω max . On the other hand, when the length of the ring decreases, the discreteness of the system's spectrum makes a great difference in the predicted observations: while some fixed number of modes stay completely within the linear dispersion region in a large-radius torus, the same set of modes can partially or even completely exceed the frequency limit for a small ring. This change is accompanied by the transition from three pairs of the correlation tongues to only HR − in surviving one, or even the absence of any correlations, for a sufficiently large number of modes above ω max ). It is expected that, in an elliptically deformed torus, the dependence of the Hawking effect on L remains the same, in the first approximation.
Another important issue concerning our system with two horizons is its stability, which is sensitive to boundary conditions and the presence of different horizons [70]. For the initial noise that does not strongly violate the frequency condition, we have found the steep horizon to be subject to eventual instability, with lifetime t decay ≃ 600 (≈ 0.438 s, in physical units) for L = 320. At larger times, simulations demonstrate that the development of the instability near WHH produces soliton-like structures, and, finally, it leads to decay of both horizons. On the other hand, the system with a smoothed horizon shows no instability (at least, in the course of long-term simulations for t ∼ 1 s.) Larger regions with a bigger distance between the horizons tend to be more stable than the smaller ones.
Due to the above-mentioned fact that the presence of the HR − P correlations tongue is determined by the number of modes below and above the critical value of ω max , the value of the lowest wavenumber is restricted by the length of the ring. Therefore, there exists a critical size making it impossible to observe the Hawking correlations, as no Bogoliubov mode may satisfy the constraint (very roughly speaking, it resembles the known property of modulational instability, which is suppressed by periodic boundary conditions if the ring's length falls below the respective critical value [71]).This size can be evaluated from the condition that the frequency, which corresponds to the lowest wavenumber pursuant to Eq. (15), is equal to ω max , see Fig. 2. For our initial parameters of BHH the critical length of the ring is L cr ≈ 73 µm, which corresponds to critical radius R cr ≈ 12 µm. It remains a challenging objective to produce evidence of the disappearance of the HR − P correlations, as the noise with the amplitude used in the above simulations suppresses all correlations at R ≃ 12 µm. Therefore, we have changed parameters of the horizon (v = 0.61) so as to make the size of the ring close to L = 320 (which is good to observe correlations), and conducted simulations for a smoothed slope at WHH. In Fig. 6 we observe the absence of HR − P correlations for the three lowest modes. On the other hand, a larger eight-fold larger region produced visible P−in correlations (and negligible HR − P ones) for the ten lowest modes.

V. CONCLUSIONS
We have investigated theoretically the possibility to generate acoustic Hawking radiation in the superfluid ring-shaped BEC. For this purpose, we have introduced the double-step potential that minimizes the emulated Hawking temperature near the WHH (white hole horizon), where instabilities may occur. The desirable region with the supersonic flow and uniform density distribution of the condensate may be designed using the spatiotemporally-modulated interaction constant. These features make the system considered here more stable and convenient for the analysis of the acoustic analog of the Hawking radiation in the rings.
We addressed basic properties of the analogue Hawking radiation in the ring-shaped matter-wave configurations with different radii. The location and shape of tongueshaped correlation patterns are well predicted by the dispersion relation for the quasi-uniform condensate in the limit of low wavenumbers. They are stable against variations of the spatiotemporal modulation of the potential and coupling constant.
An important circumstance is that the correlation pattern is sensitive to the number of modes admitted by the ring, depending on its radius, below and above the frequency limit ω max . Varying these parameters, we observed different tongue-shaped correlation patterns: from three pairs of the tongues to a single pair for a sufficiently large number of modes. The discreteness of the frequency and momentum spectra in the ring produce a dramatic effect on properties of the Hawking radiation. Below the critical radius of the ring, even the lowest mode has its frequency above the ω max , so that HR − P correlations disappear. Thus, no acoustic Hawking radiation takes place in the ring-shaped superflows with the radius falling below the critical value.
It is relevant to relate the minimum radius of sonic black hole, which can emit the acoustic analogue of the Hawking radiation, to real (astrophysical) BHs. While astrophysical BHs lose their mass through the Hawking radiation extremely slowly, the Hawking temperature dramatically increases for small BHs: T H = hc 3 /(8πGk B M ) = 6.169 × 10 −8 (M ⊙ /M ) K. In this connection, it is relevant to mention that quantum effects are believed to be crucially important for BHs with the Planck-scale mass. A well-known puzzle for the quantum theory of gravity is the final fate of such BHs. There are good grounds to assume [72] that the Hawking radiation is suppressed for sufficiently small BHs when their size, r g = 2GM/c 2 , becomes comparable to the Compton wavelength, λ C = h/(M c), associated with the BH of the Planck's mass, M ∼ m P = hc/G. Accordingly, small non-radiating primordial black holes created in great numbers in the early Universe could survive and become a key ingredient of dark matter [73][74][75][76], as it is seen today.
We hope that the minimum radius suppressing the emission of the acoustic Hawking radiation in the ring can shed new light on the unsolved puzzle of the stabilization of the primordial BHs.

VI. ACKNOWLEDGMENT
The work of B.A.M. is supported, in part, by the Israel Science Foundation through grant No. 1287/17.