Analog black holes and energy extraction by super-radiance from Bose Einstein condensates (BEC) with constant density.

This paper investigates the acoustic superradiance of the density and phase fluctuations from the single vortex state of a Bose-Einstein condensate, by employing full time-domain and asymptotic frequency domain numerical calculations. The draining bathtub model of an incompressible barotropic fluid is adopted to describe the vortex. The propagation of the axisymmetric density and phase fluctuations in the condensate are governed by the massless scalar Klein-Gordon wave equation, which establishes the rotating black-hole analogy. Hence, the amplified scattering of these fluctuations from the vortex comprise the superradiance effect. A particular coordinate transformation is applied to reveal the event horizon and the ergosphere termwise in the metric and the respective asymptotic spectral solutions. A comparative analysis of the time domain and asymptotic frequency domain results are given for a range of rotational speed of the vortex and the frequency of the impinging fluctuations. The agreement at low rotational speeds of the vortex is shown to be very good, which starts to deteriorate at higher rotational speeds due to increasing constraint violations of the time-domain calculations. We further demonstrate an asymptotic upper bound for the superradiance as a function of vortex rotational speed, provided that the vortex remains stable.

Introduction

Analogies in physics enable us to observe a particular phenomenon with the same characteristic features in different systems pertaining to disparate mechanisms and space-time-energy scales. A particular example is the analogy between a cosmic black hole and the microscopic vortex state of a Bose-Einstein condensate, which casts the superradiance phenomenon of scalar waves from the black hole to the superradiance of acoustic waves from the liquid vortex. Superradiance allows for energy and angular momentum to be extracted from the vacuum, which coincides with the initial wave to have a reflection coefficient greater than one. Superradiance phenomena in rotating black holes first discovered by Zel'dovich in 1971 [1]. He showed that the amplifications occurs at the level of event horizon when a certain condition is met. In addition, the phenomena rely on the rotating black hole, described by the Kerr metric, which exhibits two key features; event horizon, described as null surface which acts as a one-way membrane and the ergoregion, stationary limit surface [2]. Since observing a cosmic scale superradiance is not a feasible option [3], we currently rely on analogous systems that can be realized at the laboratory scale. In this case a space-time geometry Kerr black hole [4], [5], [6].

The analogy was initiated by Unruh's calculations [7], who showed the equivalence between the background solution of velocity perturbations on a perfect barotropic, irrotational Newtonian fluid and the Klein-Gordon field propagating in a 4-dimensional pseudo-Riemannian manifold, in which the speed of sound plays the role of speed of light. Because superradiance phenomenon occurs in the space-time background of rotating black holes, the analogy could be set for a rotating acoustic black-hole in a liquid [8], [9], [10], [11], [12]. Thus far, theoretical and computational investigation of the analogue superradiance have been reported for liquid systems [13], [14], relativistic fluids [15], shallow water systems [16] and optic systems [17], [18], [19], [20]. Within the context of BEC, superradiance stemming from the coherent light matter interaction [21], [22], as an analogue to the Dicke effect [23] and the amplification of matter waves as a manifestation of Raman superradiant scattering were also reported theoretically [24]. In addition, various physical phenomena such as black hole lasers have been adapted to analogue black hole systems [25], [26].

The experimental studies for that purpose emerged only within the last few years. Experimental realizations of horizons were reported in water channels [27], atomic Bose-Einstein condensates (BECs) [28]. Recently, rotational superradiant scattering in a water vortex flow is reported [29], [30]. Acoustic black hole in a needle-shaped BEC of 87Rb is realized and recently spontaneous Hawking radiation, stimulated by quantum vacuum fluctuations, emanating from an analogue black hole in an atomic Bose-Einstein condensate is reported [31], [32], [33].

In the perspective of the existing literature, different studies provide theoretical models to study the superradiance in the time domain and in the frequency domain. The methodology for the time domain solution is described for the Kerr black holes [34] and charged black holes [35] and modified for the draining-tub model of an acoustic black hole in irrotational, barotropic and incompressible fluid systems [36], [37]. While in the frequency domain, the superradiance is mainly analyzed at the asymptotic with an appropriate coordinate transformation [13], [14], [38]. The present work aims to contribute in two aspects: Providing a consolidating study of the temporal and spatial features of the scattering from a BEC vortex with constant background density and characterization of the superradiance as a function of the rotational speed of the vortex. We primarily adopt the draining bathtub model (DBT) introduced by Visser [11], to describe the acoustic black-hole in BEC vortex flow. The time domain solutions are obtained by solving the Klein-Gordon equation for the propagation of acoustic waves, whereas the spectral analysis of the superradiance is conducted by asymptotic solutions of the waves at the event horizon and the spatial infinity. The time-domain solutions are obtained by implementing the numerical techniques described mainly in Refs. [34], [37], [39]. Complementary to the cited works above, our study demonstrates very good agreement between the full time-domain and asymptotic frequency domain solutions and reveal novel spectral features within the DBT model. In particular, an asymptotic upper bound for the maximum reflection coefficient with respect to rotational speed of the vortex is obtained. The paper is organized as follows: Section 2 describes briefly the BEC system and gives a theoretical formulation leading to the main (Klein-Gordon) equation. Section 2.2 and 3 are devoted to the implementation and computation of the time-domain solutions. Section 3.1 presents the asymptotic solutions in the frequency domain. The last section discusses the main results and concludes the paper.

Model

We begin by a brief description of the Bose-Einstein condensate as the physical system of interest. A quantum system of N interacting bosons in which most of the bosons occupy the same single particle quantum state, the system can be described by a Hamiltonian of the form; Here is an external potential, is the interatomic two-body potential, m is the mass of the bosons and is the boson field operator which includes the classical contribution plus excitations , where known as the wave function of the Bose-Einstein condensate.

In the non relativistic limit, most of the atoms occupy the ground state and the classical wave function, describes the system. The interatomic interaction is taken as , , where the constant a is the scattering length constant. Closed-form equation for weakly interacting bosons, with the potential defined above leads to the time dependent Gross-Pitaevskii (GP) equation: Here in hydrodynamic form the wave function can be written in terms of its magnitude and phase: Then, the density of particles is given by and the background fluid velocity is defined as . A general review on BEC formulation can be found in [40], [41].

The Draining Bathtub Model introduced by Visser et al. [11] describes the single vortex state of the BEC, where the fluid velocity has tangential and radial components, with A and B are constants to be determined.

The density and phase fluctuations are introduced respectively as and , where and define the background. For small fluctuations, the GP equation can be linearized perturbatively, leading to the propagation equations of the fluctuations as follows: where is given by In the equation above Eq. (6), the pressure term is of the order while the quantum pressure term is of the order , where R is the spatial scale [42]. This implies that for the last term on the right-hand side of Eq. (6), is negligibly small. The right side of equation 8 thus introduces a characteristic length scale, known as the healing length. In the hydrodynamic approximation, equations (5) and (6) are combined to yield where the speed of sound is defined by . For a constant background density profile, the speed of sound, c, is constant. Equation (9) then becomes equivalent to the massless Klein Gordon equation, describing the propagation of linear density and phase fluctuations in the axially symmetric space-time. The metric associated with this space-time will be

Coordinate transformations

The coordinate transformation given below is particularly useful to minimize the number of off-diagonal components of the metric, leaving only one. where and . We drop the *-superscript in the following part of the formulation. The line equation takes the form The event horizon and the ergosphere can be identified termwise in the metric now. For a stationary and axisymmetric spacetime metric, the radius of the ergosphere is given by the vanishing of and the coordinate singularity of the metric signifies the event horizon. From Eq. (12), they read as

We assign the event horizon to a making, and and scale radial coordinate by and time coordinate by .

Numerical solution in the time domain

In order to solve the Eq. (9), first we write the line element in the form; where , and . From this point, we change the notation for to Ψ to simplify the equations below, i.e. . Equation (6) with the hydrodynamic approximation applied gives the density fluctuation as . Substituting velocity from Eq. (4) and scaling by the background density and the healing length Eq. (8) gives For numerical calculations two conjugate fields are introduced as follows where , and and are the axial and azimuthal wave numbers, respectively [34], [43]. In this work, in accordance with the BEC vortex stability conditions, we consider the azimuthal wave numbers of and 1 only [44], [45]. Using the conjugate fields, we obtain a hyperbolic system of three coupled first order partial differential equations.

Results & Discussion

We first describe briefly the implementation numerical black-hole excision adopted from Refs. [46] and [37], to solve Eqs. (17)-(19). In this technique, a numerical boundary is placed beyond the event horizon (an “apparent horizon”) which excises its interior from the computational domain, thereby removing the singularity. Since no information can leave the interior of the black hole, the excision should have no effect on the physics outside. This translates to the numerical implementation as a constrained evolution in which the constraint violations occurring beyond the event horizon has to be monitored and kept below a tolerance. Consequently, the time-domain calculations are limited by a time scale until the constraint violations start to interfere with the physical computation domain beyond a tolerance. This time scale decreases with increasing rotational speed of the vortex and thus affects adversely the global settings of the time-domain calculations.

This main section is organized as follows: We first calculate the time evolution of the perturbations of the velocity potential by solving the equation set (17)-(19) and check the consistency of the numerical result at event horizon and the outer boundary. Second, the energy of the perturbations is calculated showing the superradiant and non-superradiant cases based on the ranges of the model parameters Ω and ω.

From the fields and ψ, we construct outgoing and incoming fields along the null ray At large distances purely outgoing wave is implemented such that , . We terminate the numerical simulation before the fluctuations reach to the outer boundary so that it has no effect in the obtained results. For calculating the scattering amplitude the outgoing wave is recorded at a specified location before the outer boundary. For the numerical boundary beyond the event horizon, free-end boundary conditions are applied, such that waves are free to propagate without any constraints. We remind that the relevant condition to be monitored is not this boundary but the constraint conditions evaluated at the excision horizon, in this case the event horizon. The positioning of the numerical boundary beyond the event horizon is subject to two criteria: It must be far from the singularity to reduce the amplification of the constraint violations beyond the event horizon and also far from the actual event horizon so that the constraint violations cannot propagate to the physical domain in a short time to interfere with the actual wave.

The computational spatial (radial) and time domain are set as , with discretization steps of , , respectively. The incident wave is chosen as an axisymmetric imploding Gaussian pulse modulated by a monochromatic wave, centered at with a width of and azimuthal wavenumber [43]. Here, c is the propagation speed of sound in the condensate and a is the location of the event horizon. Both parameters are scaled to unity and the amplitude N is calculated by the normalization of the incident wave on the radial domain.

We note that the location of the incident wave should be chosen numerically far enough so that the scattering outcome is independent from the location of the incident wave. The angular speed of the vortex is Ω. In the present calculations, we consider values of up to . The frequency of the incident wave is .

Fig. 1 shows the snapshots of the initial cylindrical Gaussian wavepacket at and for both superradiant () and non-superradiant () cases. The perturbation for the non-superradiant case is close to zero while in the superradiant case it gets amplified. We monitor the constraint value C, from the definition of Φ, Eq. (16) To analyze the reflection, we put a monitor at and record the outgoing wave at this location. Fig. 2 shows the constraint violations at the event horizon (left panels) and at the recording monitor's position () (right panels) for the superradiant case (). Evidently, the violations are large when the wave passes from the respective location. For note how the constraint violations at the outer boundary are also enhanced while the wave interacts with the event horizon. However, we observe that the overall magnitude of the constraint violations are about an order of magnitude higher at the event horizon compared to the outgoing wave monitor location ().

Figure 1

Incident wave at t = 0 (left) and the reflected wave at t = 100a/c for the superradiance case, amplified (dashed blue line) and non superradiant case (red line). The parameters used are r0 = 50a and b = 10a with ω = 0.7c/a, Ω = 1.4c/a.

Figure 2

Constraint violations at event horizon (a) and outer boundary (r = 70a) (b) for superradiance (m = 1). The parameters used are in Fig. 1. Dotted lines signify the time frames (tc/a = 40,110), when the wave reaches the horizon and the outer boundary.

Similarly, Fig. 3 shows the constraint violations for the non-radiant case (). Fig. 3(a) exemplifies nicely the numerical problem associated with the excision technique: After the physical wave interacted with the event horizon () the instabilities induced at the numerical boundary beyond the event horizon start to propagate out and reach to the event horizon around . Note that by this time the physical wave has reached almost to the outer boundary and is not affected by these instabilities. If the simulation had been allowed to continue further, these amplified instabilities would propagate further and render the results obsolete. The constraint violations depend on the rotational speed of the vortex (Ω) and the relative center frequency of the impinging wavepacket (). In general, high values of Ω and increases the overall magnitude of the constraint violation at the event horizon. Also the energy associated with the superfluid flow, energy of wave packet, is given as which is calculated, normalized by the energy of the initial wave and plotted in Fig. 4 for the non-radiant (dashed blue curve) and superradiant (solid red curve) cases respectively. The incident wave arrives to the event horizon approximately at .

Figure 3

Constraint violations at event horizon (a) and outer boundary (r = 70a) (b) for non-superradiance (m = 0). The parameters used are in Fig. 1.

Figure 4

Time evolution of the energy gain of the wave packet, superradiant m = 1 case and non-superradiant m = 0 case. The parameters used are in Fig. 1.

The time evolution of the density fluctuations associated with the acoustic wave propagation are plotted in Fig. 5 for superradiant and non-radiant cases which shows the detailed view of the propagation of fluctuations near the event-horizon. Sudden increase in the density fluctuations for the superradiance case, stays inside the event horizon, .

Figure 5

Density fluctuations, ρ1 in r-t plane for superradiance case, m = 1 (a) and non-superradiance case m = 0 (b).

The amplification of energy as a function of the relative center frequency of the wave is plotted in Fig. 6, for values . The amplification gradually increases with a maximum typically in the range beyond which it decreases rapidly. The amplification curve shifts up with increasing Ω.

Figure 6

The energy calculated according to Eq. (23) normalized to its initial value E(t = 0) as a function of ω/Ω where 0 < ω < Ω. The parameters used in the calculations are , 2.8c/a, 4c/a, r0 = 50a, b = 10a. Energy values recorded when the outgoing wave reached the monitor boundary at r/a = 50.

Numerical model in the frequency domain

In this section we analyze the Klein Gordon equation (Eq. (9)) in the frequency domain. Using separation of variables, the formal solution of the KG is expressed as where k and m are the axial and azimuthal wave numbers, respectively. To avoid polydromy problems [37], that is to make ψ single valued, m should be taken as an integer and k a real number defined by the boundary conditions along the z axis.

By inserting (24) into (9), we obtain a second order ODE for the radial part: We substitute with a Regge-Wheeler tortoise coordinate, , which will map to : Introducing into Eq. (25) yields the final form of the equations for and where Near the event horizon and at the far field , the asymptotic solutions are given by the harmonic functions, where and and are the amplitudes of the reflected (transmitted) waves, respectively. Here, in to achieve scattering states should be positive.

The equality of the Wronskian of these solution at asymptotics gives It shows that when the superresonance condition, , is satisfied, reflection coefficient has a magnitude larger than unity [47], [48]. Eq. (32) reveals the superradiance condition clearly (i.e. ) and gives the full spectral behavior of the reflection coefficient. Thus, we can obtain the reflection coefficient through the Fourier components of the asymptotic far field solution, which is obtained through Eq. (26) and Eq. (28).

Fig. 7(a) shows the variation of the reflection coefficient in the Ω- plane and 7(b) shows expanded slices at different Ω values (the horizontal axis spans multiple ranges). These calculations employ the model parameters () with the same values used in the time-domain calculations of section 3.

Figure 7

a) Reflection coefficient values between 1 < Ωa/c < 4.2 and Ω/15 < ω < Ω. b) Reflection coefficients as a function of ω, calculated in the range 0 < ω < mΩ. Parameters are m = 1, and = 1, 1.75, 2.5, 3.25, 4.

Fig. 8 shows the maximum value of the reflection coefficient as a function of rotational speed of the vortex. Interestingly, the frequency domain calculations predict an upper bound for the maximum superradiance, which might be connected to entropy bound as stated in [49]. We stress however that at large rotational speeds, the stability of the vortex and hence this outcome becomes arguable. Fig. 9 the compares reflection coefficient calculated in the time-domain (section 3) and frequency domain (section 3.1), respectively for different values of Ω. As noted before, the time-domain calculations are susceptible to large values of Ω due to numerical instabilities, whereas the frequency domain is immune. The comparison justifies this observation: The agreement between the time-domain and frequency domain is very good at and in the intermediate range of values.

Figure 8

Maximum reflection coefficient behavior with respect to Ω between 1 < Ωa/c < 16.

Figure 9

Dependence of the reflection coefficient to the frequency of the incident acoustic perturbation, scaled by the vortex frequency. Upper (lower) pair of curves are for Ωa/c = 2.8(=1.4). Solid curves are obtained from direct time-domain calculations. Dash-dotted lines are obtained from asymptotic frequency domain calculations.

Conclusion

In this work, we investigated the amplified scattering of acoustic waves propagating in a BEC, from a vortex state with a constant background density, by obtaining both time-domain and asymptotic frequency domain solutions numerically. Time-domain study amounts for solving the Klein-Gordon equation which governs the radial propagation of sound waves in the presence of vortex, in analogy to scalar field propagation in the curved space-time of a black-hole. It is worth to note that the classical (macroscopic) wave function of the BEC represents the classical space-time of General Relativity only when probed at long-enough wavelengths such that it behaves purely hydrodynamically. The frequency domain calculations are performed at the radial part of the solution transformed to the (-infinity, + infinity) range, and by introducing Fourier components of the asymptotic incoming and outgoing fields. The amplitudes of these fields yield the transmission and reflection coefficients at each frequency component.

The major outcomes of the study are as follows: The comparison of full time-domain solutions and asymptotic solutions in the frequency domain show that their agreement is good near the characteristic rotational speed of the vortex Ω. As Ω increases, the constraint equations introduced by the excision technique in the time-domain calculations become increasingly violated, rendering the results with large error margins. The frequency domain asymptotic solutions do not suffer from numerical instabilities and the scattering coefficient can be calculated for arbitrarily large values of Ω. This formulation predicts an upper bound of against the rotational speed of the vortex, although the stability of the vortex at high rotational speeds is not taken into account.

As a final note, the methods presented in this work are suitable to employ non-constant background density profiles. Since the sound propagation speed will also vary, it is possible to explore novel features of the superradiance phenomenon from the BEC vortex beyond the constant-density approximation. This will be pursued in subsequent studies.

Declarations

Author contribution statement

Betül Demirkaya; Tekin Dereli; Kaan Güven: Conceived and designed the analysis; Analyzed and interpreted the data; Contributed analysis tools or data; Wrote the paper.

Funding statement

Betül Demirkaya was supported by TUBITAK-BIDEB 2211 National Scholarship Program for PhD Students.

Competing interest statement

The authors declare no conflict of interest.

Additional information

No additional information is available for this paper.