Quantum quenches, sonic horizons and the Hawking radiation in a class of exactly solvable models

TL;DR

This paper uses an exactly solvable model to analyze the emergence of Hawking radiation analogues in cold atom systems, revealing conditions for thermal and non-thermal radiation in stationary states.

Contribution

It provides a detailed microscopic analysis of Hawking radiation in a solvable many-body model, clarifying when thermal behavior occurs in analogue gravity experiments.

Findings

Stationary states can be thermal under specific conditions.

Most experimental setups produce non-thermal Hawking-like radiation.

The work offers a microscopic interpretation of the Hawking mechanism in cold atoms.

Abstract

Taking advantage of the known exact mapping of the one-dimensional Hard Core Bose (HCB) fluid onto a non-interacting spinless fermion gas, we examine in full detail a thought experiment on cold atoms confined in a quasi-one-dimensional trap, in order to investigate the emergence of the analogue Hawking radiation. The dynamics of a gas of interacting bosons impinging on an external potential is exactly tracked up to the reach of a stationary state. Under few strict conditions on the experimental parameters, the stationary state is shown to be described asymptotically by a thermal distribution, precisely at the expected (analogue) Hawking temperature. However, we find that in most experimental conditions the emerging `Hawking-like radiation' is not thermal. This analysis provides a novel many-body microscopic interpretation of the Hawking mechanism, together with useful limits and conditions for the design of future experiments in Bose-Einstein condensates.

Figures

Equations

$$\hat{H}=\int dx$$

$$⟨O_1⟩={\int }_{-\infty }^{\infty }dk⟨k|O_1(x,-i\hslash {\partial }_x)|k⟩f(k).$$

$$⟨O_2⟩$$

$${\psi }_k^0(x)=\frac{e^{ikx}}{\sqrt{2\pi }}.$$

$${\psi }_k(x,t)$$

$$⟨{\varphi }_p|{\psi }_k^0⟩$$

$${\psi }_k(x,t)\to \{\begin{array}{ll}{\varphi }_k(x)e^{-i\frac{\hslash k^2}{2m}t}& k>0\\ \sqrt{\frac{|k|}{p}}{\varphi }_{-p}(x)e^{-i\frac{\hslash p^2}{2m}t}& k<0\end{array},$$

$$\rho (x,t)={\int }_{-k_F-k_0}^{k_F-k_0}dk|{\psi }_k(x,t){|}^2,$$

$$j(x,t)=\Re \left[i\hslash {\int }_{-k_F-k_0}^{k_F-k_0}dk{\psi }_k(x,t){\partial }_x{\psi }_k^{\ast }(x,t)\right]$$

$$V(x)=\frac{V_0}{\cosh (\alpha x{)}^2},$$

$${\varphi }_k(x)$$

$$\rho (x)={\int }_{-k_F-k_0}^{k_F-k_0}dk|{\varphi }_k(x){|}^2.$$

$$\xi$$

$${\xi }_{\pm }^2$$

$$\rho (x)=\{\begin{array}{ll}\frac{Q}{\pi }\sqrt{{\xi }^2-{\xi }_{-}^2}& \text{mbox}for\xi <-{\xi }_{-}\\ \frac{Q}{\pi }\sqrt{{\xi }^2-{\xi }_{+}^2}& \text{mbox}for\xi >{\xi }_{+}\end{array},$$

$$\rho (x)=\Theta (-{\xi }_{-}-\xi )\frac{Q}{\pi }\sqrt{{\xi }^2-{\xi }_{-}^2}+\frac{Q}{2\pi }\sqrt{{\xi }^2-{\xi }_{+}^2}+\frac{Q}{2\pi }\xi ,$$

$$\rho (x)=\frac{Q}{2\pi }\left(\sqrt{{\xi }^2-{\xi }_{-}^2}+\sqrt{{\xi }^2-{\xi }_{+}^2}\right).$$

$${\rho }_{+}$$

$$v_{+}$$

$$k_F^{eff}$$

$$k_0^{eff}$$

$$-k_F^{eff}-k_0^{eff}$$

$$k_F^{eff}-k_0^{eff}$$

$$E_{ph}(T)={\int }_0^{\infty }\frac{dk}{2\pi }\frac{\hslash ck}{e^{\beta \hslash ck}-1}=\frac{\pi }{12}\frac{(k_{B}T{)}^2}{\hslash c}$$

$$f_{\text{FD}}(k)=\frac{1}{1+e^{\beta ({ϵ}_k-\mu )}},$$

$$E(T)$$

$$h(x,x^{\prime })$$

$$\{\begin{array}{ll}\frac{1}{\sqrt{2\pi }}\frac{2k}{k+\sqrt{k^2-Q^2}}{e}^{ix\sqrt{k^2-Q^2}}& k>Q\\ \frac{1}{\sqrt{2\pi }}\left[e^{ikx}-{\left(\sqrt{1+\frac{k^2}{Q^2}}+\frac{k}{Q}\right)}^2{e}^{-ikx}\right]& k<0\end{array}$$

$$\{\begin{array}{ll}1& \text{mbox}for-k_F-k_0<k<K\\ {\left(\frac{p-k}{Q}\right)}^4& \text{mbox}forK<k<k_F+k_0\end{array}$$

$$\{\begin{array}{ll}1& \text{mbox}for-k_F-k_0<k<K\\ {\left[\frac{\sinh \left(\frac{\pi }{2\alpha }(p-k)\right)}{\sinh \left(\frac{\pi }{2\alpha }(p+k)\right)}\right]}^2& \text{mbox}forK<k<k_F+k_0\end{array}$$

$${\psi }_k(x)=\{\begin{array}{ll}\frac{1}{\sqrt{2\pi }}{T}_k{e}^{ikx}& \text{mbox}fork>0\\ \frac{1}{\sqrt{2\pi }}\left[e^{ikx}+R_k{e}^{-ikx}\right]& \text{mbox}fork<0\end{array},$$

$$T_k$$

$$R_k$$

Full text

1. Introduction

General relativity is particularly distinguished for its high internal coherence and, at present days, for its high precision in predictivity, having passed a huge number of tests, which exclude several alternative theories. One of the most recent confirmation is the direct observation of gravitational waves. However, after a century, it still resists to any attempt of conciliation with the principles of quantum physics. An interesting prediction of the combination of gravitational and quantum effects is the result derived by Hawking in 1974, i.e., the fact that black holes are expected to emit a stationary flow of particles at a black-body temperature known as the Hawking temperature, thus evaporating [1]. It is a well-known fact, however, that we have no hope of detecting such a phenomenon directly from astrophysical black holes [2]. Fortunately, Hawking radiation is not related to the dynamics of the gravitational field, but it only depends on the kinematical properties of the spacetime background. In fact, it only relies upon the presence of a horizon and not on whether the black hole space-time metric satisfies Einstein’s equations. This fact led Unruh to propose, in 1981, the realization of analogue gravity models [3], in order to reproduce the same kinematical conditions of a black (and/or white) hole in a non-gravitational system. In the years after a plethora of different analogues has been investigated: flowing water, Bose-Einstein condensates (BEC), dilute gases, fiber optics, non-linear dielectrics and others [4]. Moreover, several experimental tests have been realized, even though, up to now, a definitive confirmation of the phenomenon is still lacking [5]. Among these, experiments with cold atoms BECs look among the most promising ones because of their stability properties and the good level of precision one can reach in the experiments. For both theoretical and technical reasons, these experiments are generally performed in a one-dimensional spatial geometry, corresponding to a two-dimensional spacetime configuration. In particular, the most celebrated experiments of this kind are those recently performed by Steinhauer [6, 7], claiming to have detected the Hawking radiation, even if their interpretation is still under debate [8, 9, 10, 11].

In these systems, one of the most important difficulties both from a conceptual and a theoretical point of view is due to the fact that each particular analogue model is characterized by non-linear dispersion relations governing the spectrum of the produced quantum modes, differently from the gravitational case. This implies that, in general, it is not possible to compute exactly the spectrum of the emitted quanta (which, of course, is expected to be far from a thermal one) and usually one is led to consider approximate methods in reduced regions of the spectrum. In a Letter [12], we recently proposed and summarily discussed an exactly solvable model, devised in order to get exact information on the analogue system and, at the same time, test the limits of precision of the semiclassical methods. Here we present a detailed analysis of the model, proving the results announced in the Letter.

As previously mentioned, the experiments in BECs have always been performed in a quasi-one-dimensional geometry and, therefore, we conform to this limit. We proceed as follows: initially, we consider the analogue of the vacuum of the gravitational case, that is, a stationary flux of interacting bosons. The (repulsive) interaction between particles is essential in order to describe the excitation spectrum in terms of phonons. Accordingly, we model our system as a Tonks-Girardeau gas [13]. Then, we switch on an external potential quenching the system and we let the gas evolve in time. After a while, the system will reach a new stationary state, which, under suitable conditions, may present a sonic horizon, i.e., a point where the flow passes from a subsonic regime to a supersonic one. This evolution stands as the hydrodynamical counterpart of a gravitational collapse that leads to the formation of a horizon. The usual arguments of analogue gravity then predict that, in the final stationary state, a flux of phonons escaping the horizon should be present. Such a flux is expected to be thermal at least in a certain region of the spectrum, with a temperature proportional to the gradient of the differential velocity (the difference between the flow velocity and the sound velocity) at the sonic horizon. These predictions are usually obtained by use of a semiclassical approximation, preventing the complete control of all the physical details playing a role in the phenomenon. Instead, we investigate a microscopic model constituted by interacting quantum particles which is exactly solvable, so that the full dynamics of the model can be followed analytically by use of standard many-body techniques. This is really important in order to understand the limits of validity of the semiclassical treatment on which analogue gravity arguments are based. To his purpose, we reproduce analytically an ideal experiment step by step, so to make evident the emergence of the Hawking radiation. We also compare the exact results with the ones obtained by applying the semiclassical approximation to our model and we discuss possible experimental set-ups able to provide unambiguous evidence of the presence of the analogue Hawing radiation in Bose-Einstein condensates.

A closely related model was proposed some years ago [14] in the analogue gravity framework and, by use of semiclassical arguments, the occurrence of thermal Hawking radiation in the asymptotic state was correctly predicted. Here, through the exact solution of the model, we can follow the full dynamics of the system beyond the semiclassical approximation, providing quantitative estimates for the relevant physical observables.

2. Hard Core Bose gas

2.1. Bose-Fermi mapping

Exactly solvable models of interacting bosons in one dimension have been known for many years [15]. The most celebrated off-lattice system in this class is the $\delta$-interacting Bose gas [16] defined by the Hamiltonian

[TABLE]

where $\hat{\psi}(x)$ is the usual bosonic annihilation field operator [17]. Although the full spectrum can be calculated via Bethe-Ansatz for arbitrary interaction strength $g$ [16, 18], the structure of the wavefunction is very cumbersome and few analytical results can be effectively extracted. Studying the dynamics of this system in an external potential is even more challenging and only recently some progress has been made in this direction [19, 20]. However, in the particular (but physically relevant) limit of hard core repulsion between bosons (i.e., $g\to\infty$ in Eq. (1)) the problem becomes remarkably simple and the interacting Hamiltonian can be exactly mapped into that of a free Fermi gas [13], making the analysis considerably easier. The Bose-Fermi correspondence, valid only in one dimension, is based on two observations: $i)$ choosing a basis set labeled by all possible particle configurations in real space, the matrix elements of the kinetic term of the Hamiltonian (1) do not depend whether the field operators $\hat{\psi}(x)$ commute or anti-commute; $ii)$ the hard core constraint is automatically satisfied by fermions. Due to this mapping, which holds also in the presence of arbitrary external potentials, the full energy spectrum and all the static and dynamic correlation functions involving only density operators of the Bose gas coincide with those of the free fermion model [21]. Therefore, we are led to conclude that the Hard Core Bose fluid (HCB), also known as Tonks-Girardeau gas, represents the minimal model including the essential features of the interacting Bose gas and allowing for an exact analytical and manageable solution. It is also worth noting that, although the hard core interaction is certainly a limiting case of the more general Hamiltonian (1), the physical properties of the interacting Bose gas change smoothly up to $g\to\infty$, which does not correspond to a singular limit [16]. Indeed, recent studies on the quench dynamics of the interacting Bose gas showed that the HCB limit faithfully represents the generic behavior of the model [19, 20]. Most importantly, the HCB limit can be reached in suitably designed experiments with cold atoms (namely 87Rb) as proved in Refs. [22, 23]. The properties of this model are here investigated analytically and numerically at vanishing physical temperature, in order to unambiguously identify the thermal nature of the emerging Hawking radiation.

2.2. Quantum averages

Before exploring the physical content of the model, we first recall the explicit expressions useful to evaluate the most relevant properties we are going to investigate. Each eigenstate of a single free fermion in an external potential (which does not support bound states) is uniquely labeled by a “wavevector” $k$ characterizing the asymptotic scattering state [24]. The exact eigenstates of a collection of non-interacting fermions are then written as Slater determinants of single particle wavefunctions and are labeled by the momentum distribution $f(k)$, defining the set of occupied states. In particular, the ground state is obtained by filling all the energy levels up to the Fermi wavevector $|k|<k_{F}$. The quantum averages of one- and two-body operators [17] for a given arbitrary momentum distribution function $f(k)$ are defined in the following way:

$\bullet\,$ the average of any one-body operator $O_{1}=\sum_{n=1}^{N}O_{1}(x_{n},-i\hbar\partial_{x_{n}})$ is written in terms of the single particle eigenstates $|k\rangle$ as:

[TABLE]

$\bullet\,$ For a two-body operator $O_{2}$ expressed in terms of density operators:

[TABLE]

3. Quench dynamics

Following the line of thought originally devised by Hawking in the astrophysical context [1], we discuss a “Gedankenexperiment” in order to detect the analogue Hawking radiation in Bose-Einstein condensates. First, let us consider a single free quantum particle initially set in a plane-wave stationary state of wavevector $k$

[TABLE]

Then, at $t=0$, an external repulsive potential $V(x)$ is suddenly switched on. The ensuing time evolution is formally written as

[TABLE]

where $\phi_{p}(x)$ are the exact eigenfunctions of a particle in the external potential $V(x)$, $\epsilon_{p}=\frac{\hbar^{2}p^{2}}{2m}$ are the associated eigenvalues and $\eta\to 0^{+}$ is the usual convergence factor. Now, let us look at the long time evolution at fixed position $x$. A careful analysis of these expressions shows that, for each $k$, $\psi_{k}(x,t)$ is given by the sum of an exact scattering eigenfunction of the system in the presence of the external potential multiplied by a time dependent phase factor, plus a contribution which vanishes as $t\to+\infty$. These two terms represent, respectively, the asymptotic stationary state and a traveling wave originated during the quench (see also Appendix B).

If we now take a HCB gas instead of a single particle, according to the Bose-Fermi mapping we can equivalently consider a free spinless Fermi gas, whose dynamics is described by a Slater determinant of single-particle wavefunctions. Therefore, to describe the quench dynamics in a HCB gas, we start from a Fermi gas filling all the single-particle states up to the Fermi momentum $k_{F}$. Then we set the gas into motion to the left by rigidly shifting the momentum distribution which now includes the states $-k_{F}-k_{0}<k<k_{F}-k_{0}$. In the following we will always consider $k_{0}<k_{F}$, which corresponds to a subsonic flow. Note that uniform density $\rho_{0}$ and the initial velocity $v_{0}$ of the fluid are simply related to $k_{F}$ and $k_{0}$ by $\rho_{0}=\frac{k_{F}}{\pi}$ and $v_{0}=-\frac{\hbar k_{0}}{m}$, while the corresponding sound velocity is $c=\frac{\hbar k_{F}}{m}$. Lastly, we turn on the external potential and wait, at fixed position $x$, until a stationary state is reached. During the unitary time evolution, the many-body wavefunction preserves its Slater determinant structure with time-dependent orbitals given by (5). To quantitatively discuss the long time properties of the system we have to specify the form of the external potential $V(x)$. In the following, we consider two representative model potentials often investigated both theoretically and experimentally.

3.1. Step potential

We first consider an external potential of the form $V(x)=V_{0}\Theta(x)$ where $\Theta$ is the Heaviside function and we write for convenience $V_{0}=\frac{\hbar^{2}Q^{2}}{2m}$. This “waterfall” potential [25] coincides with the original proposal by Unruh and has also been adopted in the experiments performed by Steinhauer [6, 7].

A set of single-particle eigenfunctions are explicitly reported in Appendix A. With this choice, by evaluating the long time limit of Eq. (5), the stationary state is defined by the single-particle orbitals:

[TABLE]

with $p=\sqrt{k^{2}+Q^{2}}$. A sketch of the derivation of this asymptotic expression is reported in Appendix B. The momentum distribution defining the Slater determinant remains unchanged during the evolution and, therefore, the stationary state is defined by the orbitals belonging to the interval $-k_{F}-k_{0}<k<k_{F}-k_{0}$. As a consequence, after an initial transient, an observer at a given position $x$ will perceive the moving fluid in a stationary state described by the asymptotic single-particle wavefunctions (7). Note that, while the full evolution (5) must conserve both the energy and the total number of particles, the stationary state orbitals do not have to, since the waves formed during the quench and traveling in both directions may carry energy (and particles) to infinity at long times. Moreover, in the non-uniform stationary state, the local energy and momentum density may differ even in the asymptotic regions $x\to\pm\infty$ from their initial value before the quench.

The stationary state properties can now be analytically evaluated by use of Eq. (7). In Fig. 1 we show a few snapshots of the time evolution of the density profile

[TABLE]

as obtained by numerical integration of the Schrödinger equation for the Fermi gas flowing in a waterfall potential.

Lengths are expressed in units of $Q^{-1}$ and times in units of $\tau=\frac{m}{\hbar Q^{2}}$. Here the initial density is set to $\rho_{0}=0.5\,Q$ and the initial velocity is $0.12$ times the sound speed of the fluid $c_{0}=\frac{\hbar}{m}\pi\rho_{0}$. The quench dynamics is clearly seen in the figure: two waves are generated at $t=0$ in $x=0$ and propagate at different velocities $c_{0}\pm|v_{0}|=\frac{\hbar}{m}(k_{F}\pm k_{0})$ in the downstream ($+$) and upstream ($-$) direction. In the central region the density profile shows a rapid variation near the waterfall and several undulations develop, the density being lower upstream than downstream. These waves cannot be interpreted in terms of the usual dispersive shock waves (DSW) or simple waves, as they disappear if the same dynamics is studied by means of the mean-field Gross-Pitaevskii equation. Furthermore, as Fig. 2 shows, they persist in the stationary state, differently from what we would expect in the case of DSW. They thus reflect the presence of the hard core constraint. In Fig. 2, a blow-up of the density profile in the inner region at long times $t=50\,\tau$ is compared to the analytical solution based on Eq. (7). As it can be seen, the agreement is remarkable.

The velocity profile of the Fermi gas is defined in terms of the local mass flux

[TABLE]

as $v(x,t)=\frac{j(x,t)}{m\rho(x,t)}$. In Fig. 3, the fluid velocity is compared to the local sound speed, which, for a HCB gas, is given by $c(x,t)=\frac{\hbar}{m}\pi\rho(x,t)$.

For this parameter choice the flow is always subsonic, the ratio being less than unity. The excellent agreement between the numerical results and the analytical expressions confirms the correctness of the theoretical analysis.

At fixed initial density $\rho_{0}$, for sufficiently small values of the initial velocity $v_{0}$ the stationary state is subsonic while, for larger values of $v_{0}$, a supersonic region appears near the potential step, although the flow remains subsonic in the asymptotic region $x\to-\infty$. Nevertheless, a sonic horizon is present whenever $k_{F}-k_{0}<Q$. Further increasing $v_{0}$ the flow becomes fully supersonic beyond the sonic horizon located close to $x=0$.

Fig. 4 shows the ratio between the fluid velocity and the local sound speed in the stationary state for few values of the initial parameters $(k_{F},k_{0})$. The behavior is generally non monotonic and characterized by undulations in the downstream region, except for $k_{0}=k_{F}$ (i.e., for $v_{0}=-\pi\frac{\hbar\rho_{0}}{m}$) when the oscillations disappear and the velocity becomes constant beyond the horizon. A few snapshots of the HCB dynamics are shown in Fig. 5 for a set of parameters triggering a supersonic transition.

The analytical stationary state and numerical results long after the quench are compared in Fig. 6 and the agreement is remarkable. Note that, in this case, the development of the stationary state requires considerable longer times and the numerical integration was carried out up to times as large as $t=250\,\tau$ in order to obtain a stationary solution in the range $|x|\lesssim 10\,Q^{-1}$.

The Mach ratio $|v|/c$ long after the quench is shown in Fig. 7 displaying a supersonic transition. For this choice of parameters at $t=0$ the fluid velocity coincides with the sound speed of the uniform gas.

3.2. Repulsive barrier

We now turn to a different form of the external potential: a repulsive barrier given by

[TABLE]

where the parameter $\alpha$ governs the width and smoothness of the barrier, while $V_{0}$ is conveniently parametrized as $V_{0}=\frac{\hbar^{2}Q^{2}}{2m}+\frac{\hbar^{2}\alpha^{2}}{8m}$ [26]. Taking advantage of the exact solution of the eigenvalue problem [24] a set of exact, properly normalized eigenfunctions can be chosen as:

[TABLE]

where $\zeta=\frac{1-\tanh(\alpha x)}{2}$ and $\Gamma(a)$, $F(a,b;c;\zeta)$ are the usual Gamma and Hypergeometric functions, respectively [27]. The above expression, valid for $k>0$, represents a right-moving scattering solution with energy $\epsilon_{k}=\frac{\hbar^{2}k^{2}}{2m}$. The degenerate eigenfunction with $k<0$ is obtained by replacing $k\to-k$ and $x\to-x$. We can proceed as previously discussed by setting the free Fermi gas in a Slater determinant of plane waves with momenta in the interval $-k_{F}-k_{0}<k<k_{F}-k_{0}$ (with $0<k_{0}<k_{F}$). An analysis similar to the one carried out for the step potential (see also Appendix B) shows that, with this choice of eigenfunctions, by switching on the barrier and waiting for equilibration, the system relaxes to a state defined by a Slater determinant of the eigenstates (10) with momenta belonging to the same interval. This set of eigenstates, therefore, describes the asymptotic stationary state of the system long after the quench. The analytical form of the asymptotic densities at $x\to\pm\infty$ are reported in Appendix C.

A comparison between the numerical integration of the Schrödinger equation for long times and the asymptotic expressions obtained from Eq. (10) is shown in Figs. 8 and 9 for the density and the velocity profiles, respectively.

Here the parameter $\alpha$ has been chosen to represent a rather sharp barrier ($\alpha=Q$), while the same comparison for a smoothly varying potential with $\alpha=0.1\,Q$ is presented in Figs. 10 and 11.

As for the case of the step, we conclude that the stationary state predicted on the basis of the eigenfunctions (10) is indeed reached at long times, irrespective of the value of the parameter $\alpha$. An interesting peculiarity of the time evolution in the case of the potential barrier can be readily noticed: in the step case two density modulations originate from the defect (i.e., the waterfall placed at $x=0$) propagating in the upstream and downstream directions; for a barrier, instead, as the supersonic transition sets in, a further soliton-like wave propagating downstream is clearly visible. The soliton velocity is indeed very close to the theoretical expectation $v_{sol}\sim c+|v|=\pi\frac{\hbar Q}{m}$. This observation confirms the analysis, based on the Gross-Pitaevskii dynamics, carried out in Ref. [28] where the effect was attributed to the peculiar form of the Bogoliubov excitation spectrum in the supersonic case. Note also that the stationary state properties of the system are visibly different from those of the step potential, while they barely change with $\alpha$ in the range we have examined (modulo a trivial re-scaling of the length unit), showing that, in the case of a potential barrier, the “smooth limit” is easily achieved also for moderate values of $\alpha/Q$.

Remarkably, in the $\alpha\to 0$ limit, several analytical expressions can be obtained from the asymptotic behavior of the Hypergeometric functions (see Appendix D). The density profile in stationary state is readily evaluated by inserting the results (141, 142) into

[TABLE]

It is convenient to define a dimensionless coordinate as

[TABLE]

The density profile has different analytical expressions in three regimes:

• •

$Q>k_{F}+k_{0}$. Here all the fermions have kinetic energy lower than the height of the barrier. Total reflection occurs and the stationary state density profile is formally given by the expression

[TABLE]

while $\rho(x)=0$ elsewhere. The mass flux in this regime vanishes because particles cannot tunnel through the barrier.

• •

The most interesting regime is for $k_{F}-k_{0}<Q<k_{F}+k_{0}$. Here quantum tunneling occurs and the constant mass flux is $j=-\frac{\hbar}{4\pi}\left[(k_{F}+k_{0})^{2}-Q^{2}\right]$. The density profile is now given by

[TABLE]

where the Heaviside function $\Theta(x)=1$ for $x>0$ and vanishes for $x<0$. Note that the density profile is generally non monotonic. A sonic horizon is always present in this regime.

• •

For $Q<k_{F}-k_{0}$ the flow is fully subsonic in the stationary state. The mass flux is $j=-\frac{\hbar}{\pi}k_{F}k_{0}$ and the full density profile preserves the symmetry of the potential:

[TABLE]

A comparison between the exact solution and the asymptotic analytical results for $\alpha\to 0$ is shown in Fig. 12 where the local Mach number $\frac{v(x)}{c(x)}=\frac{j}{\pi\hbar\rho(x)^{2}}$ is displayed for $\frac{\alpha}{Q}=0.1$ and different choices of the parameters $(k_{F},k_{0})$.

The analytical expressions are able to capture the essential features of the exact solution, showing that the Mach number develops a maximum in the supersonic region, while downstream the flow may even become subsonic again (blue and red curves), in close analogy with the case of the step potential (see Fig. 4). After the peak, the density profile (and therefore also the Mach number) shows oscillations on a scale $x\sim\alpha^{-1}$ which are washed out in the $\alpha\to 0$ limit (dashed lines). The oscillations disappear both in the fully subsonic regime (cyan curve) and for $k_{0}=k_{F}$, when the Mach number is always monotonic (black curve). From the analytical expressions we can also evaluate the asymptotic uniform density and velocity in the far upstream region $x\to+\infty$: while in the subsonic case $Q<k_{F}-k_{0}$ the asymptotic density and velocity remain unchanged during the quench ($\rho_{+}=\frac{k_{F}}{\pi},v_{+}=-\frac{\hbar k_{0}}{m}$), in the interesting regime $k_{F}-k_{0}<Q<k_{F}+k_{0}$ they become

[TABLE]

Therefore, a hypothetical observer in the far upstream region, unaware of the presence of the potential barrier, would assign an effective value to the key parameters $(k_{F},k_{0})$ given by

[TABLE]

and then an effective Fermi distribution limited by the two Fermi momenta

[TABLE]

At this point, before studying the emergence of the analogue Hawking radiation in the cases presented above, it is important to note how this exact model already displays important differences from the usual semiclassical approach, for both the cases of a waterfall and that of a barrier potential in the $\alpha\to 0$ limit. Indeed, in the presence of a sonic transition, our model shows that the density and the velocities profiles can acquire several different forms in the stationary state; this, though, is profoundly different from the semiclassical approach usually adopted, which, in these cases, only admits monotonic stationary profiles, thus neglecting a whole class of solutions. Here, instead, several different dynamical features which are usually hidden can be investigated.

4. Analogue Hawking radiation

At this point, we are ready to investigate the occurrence of thermal phonons emerging from the sonic horizon in the stationary state of the model. We stress that, while in the gravitational framework the background curved metric originates from the presence of a black hole which is unrelated to the quantum field giving rise to the Hawking radiation, in the condensed matter analogy, both the effective metric and the scalar field are expressions of the same quantum fluid: the flowing Bose gas. The phonon field, which plays the role of the quantum scalar field in this framework, can indeed be considered decoupled from the fluid flow only if we limit our analysis to low energy, where the elementary excitations are known to behave as free quantum quasi-particles [29].

The energy density of a one-dimensional uniform Bose fluid in the presence of a thermal phonon branch is known to exhibit an additive contribution of the form

[TABLE]

where the usual phonon dispersion relation $\omega=c\,k$ has been used. Switching to the representation of the HCB in terms of spinless fermions, the same fluid can be described in terms of an effective Fermi-Dirac distribution

[TABLE]

leading to a thermal contribution to the energy density

[TABLE]

The first term is the ground state energy of the Fermi gas, while the second equals the thermal energy density of the two phonon branches (one for each Fermi point), as can be easily checked recalling that the sound velocity in a one-dimensional Fermi gas (or a HCB fluid) coincides with the Fermi velocity, $c=v_{F}=\frac{\hbar k_{F}}{m}$. This simple observation shows that, in one dimension, there is a one-to-one correspondence between the quantitative description of the one-dimensional HCB system in terms of a phonon gas and that of a fluid of fermionic particles, if only low energy excitations are present.

Moreover, if phonons originate from the sonic horizon and we study the flow in the far upstream region $x\to+\infty$, only quasiparticles with positive wavevectors $k>0$ will appear while left moving particles (i.e. with $k<0$) would be unaffected by the Hawking mechanism and remain at zero temperature. Therefore, we are led to a very specific expectation for the effective momentum distribution of the fermions in the region $x\to+\infty$ if analogue Hawking emission occurs: all local physical observables, in fact, should appear as if the Fermi gas was characterized by a momentum distribution of the form (24) with $\beta=(k_{B}T_{H})^{-1}$ for $k>0$ and $\beta=\infty$ for $k<0$.

Finally, another signature of the particle emission at the horizon due to the Hawking process is related to the existence of quantum correlations between the upstream and the downstream regions [30]. Physically, these correlations originate from the creation of entangled phonon pairs at the horizon propagating in opposite directions. In our microscopic model, we can check also this prediction by evaluating the density-density correlations in the stationary state:

[TABLE]

where the local density $\rho(x)=\langle\hat{\rho}(x,t)\rangle=\langle\hat{\psi}^{\dagger}(x,t)\hat{\psi}(x,t)\rangle$ is defined in Eq. (8). We now turn to the evaluation of these properties for the representative choices of potentials previously introduced.

4.1. Step potential

As already mentioned, in the far upstream region $x\to+\infty$ the asymptotic wavefunction $\psi_{k}(x,t)$ (7) simplifies and it becomes the superposition of an incident and a reflected plane wave (this is true regardless of the form of the potential). Omitting the overall time dependent phase factor, the asymptotic expression of $\psi_{k}(x,t)$ for a step potential is:

[TABLE]

while for $0<k\leq Q$ the wavefunction is exponentially small at $x\to+\infty$. When the quantum average of a physical quantity is evaluated starting from this expression, the asymptotic result at large $x$ can be formally written according to Eqs. (2,3). Disregarding the interference term between the two counter-propagating waves whose contribution vanishes for $x\to+\infty$, the limiting form coincides with that of a uniform free Fermi gas characterized by an effective momentum distribution $f(k)$ given by:

[TABLE]

where we have defined $p=\sqrt{k^{2}+Q^{2}}$, $K=\sqrt{(k_{F}-k_{0})^{2}-Q^{2}}$ for $k_{F}-k_{0}>Q$ and $K=0$ elsewhere. This momentum distribution displays a tail at positive wavevectors, denoting the presence of quasiparticles (phonons) traveling upstream in the stationary state. The analytic form of the tail coincides with the reflection coefficient of the external potential $V(x)$. According to this expression, for $k_{F}-k_{0}>Q$ the flow is purely subsonic, $K>0$ and the momentum distribution preserves the sharp discontinuity at $k=K$. Conversely, when a supersonic transition is present, $K=0$ and the two branches of $f(k)$ join smoothly at $k=0$. Although the qualitative behavior of the system conforms to the expectations based on the gravitational analogy, the quantitative details do not: a phonon flux is still present even in the absence of a sonic horizon and, most importantly, the “Hawking-like” radiation is never thermal, because the effective distribution differs from the Fermi-Dirac form (24). This result denotes a failure of the gravitational analogy, which, being based on semiclassical arguments, is not expected to faithfully represent the actual behavior of the model when a rapidly varying external potential, like a step, is switched on.

This analysis can be generalized for the case of a smooth step of the form $V(x)=\frac{V_{0}}{2}\left[1+\tanh\alpha x\right]$ where the parameter $\alpha$ controls the sharpness of the potential and the discontinuous step is recovered for $\alpha\to\infty$. By defining $V_{0}=\frac{\hbar^{2}Q^{2}}{2m}$ and calculating the reflection coefficient for this case [24], the effective momentum distribution turns out to be

[TABLE]

where $p$ and $K$ have been previously defined. We see that the general features of this effective momentum distribution remain unchanged for a smooth step, even in the limit $\alpha\to 0$, when the phonon tail reduces to $e^{-\frac{2\pi k}{\alpha}}$. Thus, in the case of an external potential which has a step form (smooth or sharp) we have emission of phonons in the $x\to+\infty$ region from the horizon but the spectrum is never thermal. Furthermore, the emission persists even if a supersonic transition is absent and disappears only for $\alpha\to 0$.

The stationary state density correlations (4) can be easily computed starting from the asymptotic form of the wavefunctions (7). It is interesting to compare the results in the fully subsonic regime ($k_{F}-k_{0}>Q$) and the ones in the supersonic regime ($0<k_{F}-k_{0}<Q$). In Figs. 13 and 14 the results are shown in the case of a sharp step for two representative choices of the parameters.

The presence of weak, but clearly-visible, correlations across the sonic horizon in Fig. 14 confirms the commonly accepted picture: phonon pairs created at the horizon and propagating upstream and downstream give rise to density correlations between the inner and the outer regions. What is more interesting, though, is that this pattern appears only for a supersonic transition (in contrast with the particle emission which is always present) but it emerges even if the radiation emitted is not thermal.

4.2. Repulsive barrier

For the repulsive barrier (9) the asymptotic wave function $\psi_{k}$ in the upstream region $x\to+\infty$ takes the form

[TABLE]

where, again, we omitted an overall time dependent phase factor, and

[TABLE]

are the transmission and reflection amplitudes. As before, we can regard this asymptotic form as an equivalent Fermi gas with effective momentum distribution $f(k)$ given by

[TABLE]

with

[TABLE]

Note that $f(k)$ is continuous in $k=k_{F}-k_{0}$ only for $k_{F}=k_{0}$. For a generic value of $\frac{\alpha}{Q}$, the effective momentum distribution (30) differs from the expected Fermi-Dirac form (24), reflecting the coupling between the quasiparticles (phonons) and the underlying metric (flowing Bose gas). Only for $\frac{\alpha}{Q}\to 0$, i.e. for very smooth barriers, phonons are excited at extremely low energies and quasiparticles behave as a free scalar field. In this limit, the effective momentum distribution acquires the suggestive form

[TABLE]

for $k_{F}-k_{0}<k<k_{F}+k_{0}$, while $f(k)=1$ for $-k_{F}-k_{0}<k<k_{F}-k_{0}$. This distribution describes a Fermi gas with a sharp jump at the left Fermi point $-k_{F}-k_{0}$, while a tail appears for $k>k_{F}-k_{0}>0$, indicating the presence of excited phonons in the HCB fluid. First, we note that both in the fully subsonic ($k_{F}-k_{0}>Q$) and in the fully supersonic regime ($Q>k_{F}+k_{0}$) the effective Fermi point $k=Q$ lies outside the interval where Eq. (32) holds. Then, for $\alpha\to 0$, $f(k)$ either vanishes (for $k_{F}-k_{0}>Q$) or is identically equal to one (for $Q>k_{F}+k_{0}$). In both cases, for $\alpha\to 0$ the effective distribution coincides with the standard zero temperature result and no phonon flux is present in the HCB fluid at $x\to\infty$. Instead, in the interesting regime where the sonic horizon forms ($k_{F}-k_{0}<Q<k_{F}+k_{0}$), we can quantitatively match the result (32) with the expected Fermi-Dirac distribution (24) at finite temperature. Indeed, let us start from Eq. (22) and linearize the energy spectrum near the effective Fermi point $k=Q$, corresponding to the chemical potential $\mu=\epsilon_{Q}$, giving the quasiparticle dispersion at low energy:

[TABLE]

with quasiparticle velocity $v_{qp}=\hbar Q/m$. Inserting such a form into the Fermi distribution (24), we obtain the expression Eq. (32) if the effective temperature $T_{H}$ is given by

[TABLE]

Therefore, as previously discussed, an observer at $x\to+\infty$ will detect a phonon field at the temperature $T_{H}$ on top of the flowing HCB fluid. As we will see, expression (34) coincides with the Hawking temperature predicted by the analogue gravity description.

We can also evaluate the density correlations (4) on the basis of the analytical form of the scattering states (10). A representative result is shown in Fig. 15 for $\alpha=0.1\,Q$ and $k_{F}=k_{0}=\frac{\pi}{5}\,Q$. The fluid velocity for the stationary state is also reported (black line) together with the local sound velocity (red line). In the bottom-right corner a blow-up of the off-diagonal correlations is displayed in order to appreciate the presence of the weak signal related to the emergence of the Hawking radiation.

An analytical understanding of the density correlations is possible in the asymptotic region $(x,x^{\prime})\to\pm\infty$ in the limit $\alpha\to 0$. Focusing on the interesting range of parameters $k_{F}-k_{0}<Q<k_{F}+k_{0}$ where the horizon is present and taking advantage of the explicit expressions of the transmission and reflection coefficients (28, 29), the asymptotic correlations can be written in terms of $k^{\pm}=k_{F}\pm k_{0}$ as

• •

$x\to+\infty,\,x^{\prime}\to+\infty$. Here the fluid is subsonic and the correlations become translationally invariant asymptotically. Setting $s=x-x^{\prime}$ we obtain, to leading order in $\alpha$,

[TABLE]

• •

$x\to-\infty,\,x^{\prime}\to-\infty$. In the downstream region density correlations show a more complex structure as a function of $s=x-x^{\prime}$

[TABLE]

• •

$x\to\pm\infty,\,x^{\prime}\to\mp\infty$. To leading order in $\alpha$, the off-diagonal correlations identically vanish far from the horizon. However, a non trivial result appears to second order, proving that density correlations between specular points in the sub- and super-sonic regions persist even at large distances from the barrier. Setting $u=x+x^{\prime}$, the asymptotic form of the correlation function is

[TABLE]

Remarkably, this exact result for the Tonks-Girardeau gas reproduces the semiclassical prediction of Ref. [31] obtained in a mean field model where the effective coupling $g$ of the Gross-Pitaevskii equation acquires different values in the upstream and downstream regions.

5. Semiclassical analysis

A semiclassical description of the Tonks-Girardeau gas can be obtained considering a Bose fluid in the second quantization formalism, described by the Hamiltonian

[TABLE]

where $W=W_{c}+W_{e}$ is the sum of a contact term

[TABLE]

and the interaction with an external potential

[TABLE]

For future reference, we consider a generic $\nu$-body interaction so that the usual pair potential corresponds to the choice $\nu=2$ in Eq. (39). The bosonic field operator $\hat{\psi}(x)$ satisfies the canonical commutation relations

[TABLE]

We can decompose the field operator as a sum of a background configuration, described by a complex function $\psi(x)$, and a quantum perturbation so that

[TABLE]

provided that $\delta\hat{\psi}(x)$ satisfies the relations (41). Inserting in the Hamiltonian and keeping terms up to order two in the perturbation, we get $\hat{H}=E+\hat{H}_{1}+\hat{H}_{2}$, with

[TABLE]

where $E$ is the reference “classical” energy and

[TABLE]

are the first quantization effective Hamiltonians. The stationary background configuration $\psi(x)$ is chosen so that the first order contribution $\hat{H}_{1}$ identically vanishes:

[TABLE]

This non-linear differential equation defines the stationary solution in the semiclassical approximation and corresponds, for $\nu=2$, to the known Gross-Pitaevskii equation for the condensate wavefunction. While the obvious choice $\nu=2$ corresponds to the physical many body Hamiltonian (1), it has been shown [32, 33, 34] that the semiclassical approximation reproduces the exact spectrum of the strongly interacting Tonks-Girardeau gas in one dimension for the alternate choice:

[TABLE]

These are the values we will consider from now on.

Note that this change in the dynamical equations does not affect the validity of the gravitational analogy. It is possible, in fact, to demonstrate the existence of an effective acoustic metric also starting from this generalized Gross-Pitaevskii equation for any $\nu$. The procedure implies the usual derivation of the acoustic metric [35] which turns out to have the canonical form, apart from the conformal factor.

5.1. The stationary configuration

We will now assume an external smooth potential of the form

[TABLE]

in the limit $\alpha\to 0$. It is convenient to introduce the dimensionless variable $z=\alpha x$. In order to solve the generalized Gross-Pitaevskii equation (48) with $\nu$ and $g$ just introduced, we make the Ansatz

[TABLE]

Inserting this form in Eq. (48), at first order in $\alpha$ we get the equation

[TABLE]

where a prime indicate derivative w.r.t. $z$. This expresses mass conservation, which, for a stationary flow in one dimension, leads to the constancy of the mass current

[TABLE]

Solving for $\phi^{\prime}(z)$ and substituting in the zeroth order equation obtained from (48) and (51), we get

[TABLE]

with

[TABLE]

The chemical potential must therefore satisfy the constraint

[TABLE]

which allows for two solutions for each compatible choice of $\mu$ and $j$. Recalling that the local velocity $v(z)$ and the sound velocity $c(z)$ for a Tonks gas are

[TABLE]

with equilibrium density $\rho(z)=A^{2}(z)$, we get for the Mach number $\beta=v/c$:

[TABLE]

The upper and lower signs correspond to a subsonic and a supersonic velocity profile respectively, and no transition (horizon) appears in general, for any potential. The only wavefunction describing a sonic transition is obtained by matching the two solutions at a point $z_{0}$ where $\mu=U(z_{0})+|\sigma|$. In order to join the solutions keeping $\beta(z)$ real, $z_{0}$ must be a maximum for the potential: $U(z_{0})=U_{max}$.

We can now compare the resulting stationary state of the semiclassical solution with the exact one in the two cases described in the previous Sections.

5.1.1. Step potential

Let us consider first a smooth step potential of the form

[TABLE]

in the limit $\alpha\to 0$. In this case there are no local maxima and, therefore, it is impossible to match the two branches; thus, we are left with either a fully supersonic or a fully subsonic regime.

When $Q>k_{F}-k_{0}$, in the exact solution the Mach number approaches $1$ when $x\to+\infty$. Therefore, from (58), we have to set $\mu=U_{0}+\sigma$. Parametrizing the mass flux in terms of $k_{F}$ and $k_{0}$ as $j=-\frac{\hbar}{4\pi}(k_{F}+k_{0})^{2}$, we obtain the asymptotic values of the density profile at $z\to\pm\infty$:

[TABLE]

where $\pm$ in the second formula correspond to the upper or the lower sign in (54). The semiclassical solution reproduces correctly the exact asymptotic density at $+\infty$, while at $-\infty$ it gives the correct profile only for $k_{F}-k_{0}=0$ (choosing the lower sign) or $k_{F}-k_{0}=Q$ (choosing the upper sign).

If $Q<k_{F}-k_{0}$ the exact solution is subsonic and we have to select the upper sign in (54) while $\mu$ can always be chosen so that the semiclassical solution fully reproduces the exact density profile.

It is important to note that the absence of a sonic horizon for a waterfall potential in the $\alpha\to 0$ limit is not a peculiarity caused by the slight modifications introduced in the Gross-Pitaevskii equation. Indeed, it can be demonstrated that also the canonical Gross-Pitaevskii equation (i.e., setting $\nu=2$ and $g$ the usual interaction parameter) does not admit a sonic transition for any choice of the parameters.

5.1.2. Repulsive barrier

Consider a barrier displaying a unique maximum in $z_{0}$:

[TABLE]

and choose $\mu=\frac{\hbar^{2}Q^{2}}{2m}+|\sigma|$. We can then take the lower sign for $z<z_{0}$ and the upper sign for $z>z_{0}$ to get a solution that passes from a supersonic regime on the left of $z_{0}$ to a subsonic regime on the right. Let us compare this solution with the exact expressions for the potential (9) previously discussed. The exact solution shows that a subsonic/supersonic transition is present for $k_{F}-k_{0}<Q<k_{F}+k_{0}$. Choosing $Q$ in this range, we parametrize the mass current in terms of $k_{F}$ and $k_{0}$ as $j=\frac{\hbar}{4\pi}[(k_{F}+k_{0})^{2}-Q^{2}]$. The density profile in semiclassical approximation $\rho(x)=A(z)^{2}$ is monotonic and its asymptotic values at $z\to\pm\infty$ are

[TABLE]

which for the given values of $\mu$ and $j$ give

[TABLE]

For a barrier of the form (9), $\rho_{+}$ coincides with the exact value, while $\rho_{-}$ is correct only for $k_{F}=k_{0}$, when the exact density profile is monotonic. Note that, in this particular case, not only the asymptotic values but also the full density profile coincides with the exact one displaying a sonic horizon at $z=0$. For future reference we note that the slope of the relative velocity at the horizon $\kappa=\frac{d}{dz}(c(z)+v(z))|_{z=0}$ is $\kappa=\frac{\hbar Q}{m}$. It is interesting to consider also the fully subsonic case $Q<k_{F}-k_{0}$, where no matching procedure is required. Correspondingly, we choose the upper sign in (54) and set $|j|=\frac{\hbar k_{F}k_{0}}{\pi}$. Now $\mu$ is unconstrained and if we parametrize it as $\mu=\frac{\hbar^{2}}{2m}(k_{F}^{2}+k_{0}^{2})$, we get that the exact density profile is reproduced by the semiclassical solution for any $k_{0}<k_{F}$.

Remarkably, in the fully subsonic case, for both the barrier and the step potential, the exact density profile in the $\alpha\to 0$ limit is always reproduced by the semiclassical one. Instead, if a horizon is present, the semiclassical approximation is correct only for strictly monotonic velocity profiles, i.e. for a very specific (fine tuned) choice of the initial velocity corresponding to a sonic flow.

Finally, it is worth to mention that we have just looked for the stationary solutions of the semiclassical equations but we have not checked whether the semiclassical dynamics does indeed drive the system towards such a solution. In the Heisenberg formalism the field operator $\hat{\psi}(x,t)$ evolves in time according to

[TABLE]

which, after separating the background configuration from the quantum fluctuations, gives

[TABLE]

where $h_{1}(x,t)$ and $\hat{H}_{2}(t)$ are formally given by Eqs. (45-47) evaluated in terms of the evolving condensate wavefunction $\psi(x,t)$. Assuming a uniform condensate $\psi(x)=\sqrt{\frac{k_{F}}{\pi}}e^{-ik_{0}x}$ as initial condition, in the $\alpha\to 0$ limit we can look for a solution of the form

[TABLE]

expressed in terms of the dimensionless variables

[TABLE]

The equations of motion for the background field in presence of the barrier (9) now become

[TABLE]

These equations have been numerically solved verifying the approach to the stationary configuration as shown in Figs. 16 and 17 in the subsonic and supersonic case respectively. The semiclassical dynamics shown in Fig. 17 can be compared with the exact one (for a barrier of width $\alpha=0.1\,Q$) of Fig. 10. While in the exact dynamics the density profile is always smooth, in the semiclassical approximation (appropriate for $\alpha\to 0$) two shocks traveling in opposite directions are clearly visible. In both cases, long after the quench, the solution approaches the analytical, stationary state profile. Note that in the semiclassical dynamics the time scale for obtaining a stationary state in a given region of space is roughly a factor two smaller than in the exact solution.

5.2. The excitation spectrum

According to the semiclassical approach, phonons are expected to be described by the quantum perturbation around the stationary background configuration, governed by the Hamiltonian $\hat{H}_{2}$. The strategy is therefore to diagonalize $\hat{H}_{2}$ via a Bogoliubov transformation

[TABLE]

where the annihilation and creation operators $\hat{a}_{\lambda}$, $\hat{a}^{\dagger}_{\lambda}$ satisfy the usual commutation relations

[TABLE]

whereas the transformation coefficients must be normalized as

[TABLE]

By imposing the diagonalization of $\hat{H}_{2}$ in the form

[TABLE]

where $\epsilon$ is the fluctuation energy, we are led to the eigenfunction equations

[TABLE]

where we omitted the dependence on $x$ for brevity. Of course, after solving for $D_{\lambda}$, $E_{\lambda}$, and $\omega_{\lambda}$, one must impose the normalization conditions (75). Notice that if $(\omega,D,E)$ is a solution satisfying the normalization condition, then $(-\omega,E^{*},D^{*})$ is another formal solution of the same equations, but with negative frequency and negative norm. Moreover, a stable equilibrium solution of (48) must allow only for positive eigenfrequencies [36].

The fluctuation energy is expressed in terms of the solution of the eigenvalue problem as:

[TABLE]

We can apply this general procedure to two physically relevant cases: the homogeneous state before the quench and the stationary state solution long after the quench.

Before the quench, in the absence of any external potential, the background configuration is

[TABLE]

corresponding to a uniform density $\rho=\frac{k_{F}}{\pi}$ flowing with negative velocity $v=-\frac{\hbar}{m}k_{0}$. Substituting this wavefunction in the eigenvalue equations, an exact solution for the coefficients $D_{p}(x)$ and $E_{p}(x)$ is given by plane waves. Correspondingly, the fluctuation operator is

[TABLE]

with healing length $\xi=k_{F}^{-1}$. The dispersion relation for the excitation spectrum is

[TABLE]

where $c=\frac{\hbar}{m}k_{F}$. This form coincides with the known Bogoliubov phonon spectrum in a moving fluid $\omega(p)$ [4, 36].

After the quench the background solution evolves and at later times is defined by the non-homogeneous, asymptotic, stationary solution of the generalized Gross-Pitaevskii equation (51, 53, 54). In this state, phonons are defined as the normal modes of the fluctuation operator on top of the new background solution. The analysis can be carried out along the same lines followed in the uniform case. In order to solve the eigenfunction equations we make the Ansatz

[TABLE]

where we omitted the label $\lambda$, and we set, as usual, $z=\alpha x$. Inserting this form in the linear system (78), we obtain an algebraic equation for the spatial derivative of the phase $\chi^{\prime}(z)$:

[TABLE]

where

[TABLE]

are the velocity profile, the sound velocity and the local healing length, respectively. The solution of the secular equation (87) gives two branches

[TABLE]

where we omitted the explicit $z$ dependence at right hand side. The corresponding eigenmodes $(D_{\lambda}^{0}(z),E_{\lambda}^{0}(z))$ are, at fixed $\omega_{\lambda}$,

[TABLE]

where the label $\lambda$ uniquely identifies each solution of the eigenvalue problem and $D^{0}_{\lambda}$, as determined by the normalization condition (75), is

[TABLE]

The normalization condition leads to the condition $\Gamma^{2}_{\lambda}(z)<1$, which, due to Eq. (90), forces the choice of the upper sign in Eq. (88). At fixed $\omega>0$, the analysis of the algebraic equation (88) for $\chi^{\prime}$ shows that two solutions are always present both in the subsonic and in the supersonic case. However, in the latter case, there are two solutions also for any $\omega>-\omega_{max}(z)$, where

[TABLE]

and $\beta(z)=\left|\frac{v(z)}{c(z)}\right|$ is the local Mach number. This shows that the excitation spectrum and the characteristic frequency $\omega_{max}$ remain finite in the limit $\alpha\to 0$; therefore, the quantity $\alpha\kappa/\omega_{max}\to 0$ (the factor $\alpha$ in the ratio is due to our definition of $\kappa$, which is taken as the derivative of the differential velocity w.r.t. the rescaled length), satisfying the criterion established in Ref. [37] for the occurrence of the analogue Hawking radiation. As previously discussed, in the presence of a sonic horizon $\beta(z)$ is monotonic. For negative frequencies $\omega<0$ the mode amplitude diverges at $z^{*}>0$ such that $\omega=-\omega_{\max}(z^{*})$ and vanishes for $z>z^{*}$. This discontinuous behavior is likely to be an artifact of the $\alpha\to 0$ limit: we expect that for small but non vanishing $\alpha$ the divergence disappears and the mode decays exponentially for $z>z^{*}$.

The analysis of the Bogoliubov spectrum then shows that the fluctuation operator, long after the quench, is explicitly given in terms of the phonon creation and annihilation operators by

[TABLE]

where $\Gamma_{\lambda}(z)$ is defined in Eq. (90). Here, $\hat{b}_{\lambda}$ and $\hat{b}^{\dagger}_{\lambda}$ are related to the modes $\hat{a}_{\lambda}$, $\hat{a}_{\lambda}^{\dagger}$ before the quench by the non-trivial time evolution of the system after switching on the external barrier.

5.3. Analogue Hawking radiation

To relate the phonon operator before and long after the quench we have to consider the evolution from the initial uniform flow, at negative times, to the final stationary state long after the sudden switching on of the smooth barrier (9). In the Heisenberg picture, the state $|\Phi\rangle$ does not evolve and remains the vacuum state of the Bogoliubov modes $\hat{a}_{p}$. The fluctuation operator at any time $t$ can be expressed as a linear combination of the bare modes with time dependent coefficients:

[TABLE]

Substituting this expansion into the equation of motion (67) we get the equations for the coefficients

[TABLE]

with initial condition at $t=0$ given by (81). Looking for a solution that remains regular in the limit $\alpha\to 0$, we are led to the Ansatz

[TABLE]

which indeed satisfies the evolution equations (97) provided that the local momentum $\chi^{\prime}_{p}(z,\tau)$ is given by

[TABLE]

where

[TABLE]

are the time dependent velocity, sound velocity and healing length of the evolving background configuration. $\rho(z,\tau)$ and $\phi^{\prime}(z,\tau)$ are solutions of the equations (70, 71), while the initial condition at the quench $\tau=0$ is $\chi^{\prime}_{p}(z,0)=p$.

For a parameter choice corresponding to a fully-subsonic, stationary state configuration ($Q<k_{F}+k_{0}$) the numerical study of these evolution equations shows that Eq. (100) has a smooth, regular solution for both positive and negative quasiparticle momenta $p$ approaching, at long times, the expected positive norm solution (88) with the same frequency $\omega_{p}$ (84). Two representative examples are shown in Fig. 18.

In this case we conclude that $\hat{b}_{\lambda}=\hat{a}_{p}$, with $\omega_{\lambda}=\omega_{p}$, so the bare operators coincide with the proper normal modes also in the final stationary solution.

The behavior is qualitatively unchanged also when a sonic horizon forms in the background solution (i.e. for $k_{F}-k_{0}<Q<k_{F}+k_{0}$), provided the quasiparticle momentum $p$ is negative. Instead, for positive $p$ the numerical solution of the differential equation (100) shows the emergence of singularities, implying that the above Ansatz does not correctly describe the actual evolution of the phonon modes.

A careful numerical study shows that, just after the quench, while the background configuration approaches the stationary solution in a neighborhood of $z=0$, the excitations $D_{p}(z,\tau)$, $E_{p}(z,\tau)$ are regular and $\chi^{\prime}_{p}(z,\tau)$ preserves the sign of its initial condition $\chi^{\prime}_{p}(z,0)=p$, as illustrated in Fig. 19. Singularities develop later in time, suggesting that the bare phonon operators $\hat{a}_{p}$ differ from the normal modes defined in the stationary background long after the quench. Thus, in this case, we have to identify the unitary transformation relating the initial and the final phonon operators.

To this end we follow the same argument put forward by Hawking in his seminal paper [1]: We compare the forward evolution just discussed with the solution of the same equations backward in time, imposing the initial condition at a time $\tau\gg 0$, where the stationary state after the quench has been already reached in a wide portion of space around $z=0$. The fluctuation field operator can then be expressed at all times either in the form (95) or as

[TABLE]

Matching these two expressions allows to express the bare modes in terms of the phonon operators long after the quench:

[TABLE]

with

[TABLE]

In these formulas $D$ and $E$ are functions of both $x$ and $t$, but the time dependence must disappear in the final expressions. Therefore, we can choose a convenient time to evaluate the transformation matrices. We fix an optimal time $\tau_{0}$ when the forward evolution after the quench has led the background function to approximate the asymptotic one in a given neighbourhood of $z\sim 0$, let us say $|z|<0.5$ (see Fig. 19).

To evaluate the unitary transformation (106) we have to solve the evolution equations (97) both forward and backward in time. It is convenient to consider wave packets of the form

[TABLE]

where $f(p)$ and $g(\lambda)$ are weight functions chosen in such a way that, at the initial condition (i.e. $\tau\to-\infty$ for the forward evolution and $\tau\to+\infty$ for the backward one) the wave packet is centered around $z\sim 0$ (i.e. close to the sonic horizon) with momenta $p$ and $\chi^{\prime}_{\lambda}(z=0)$, respectively. These wave packets are expected to move away from $z=0$ with a group velocity given by the derivative of the frequency with respect to the wavevector. Now we have to evaluate expressions like

[TABLE]

entering the transformation matrices. Note that in order to give a significant contribution to the integral, both wave packets, at time $\tau_{0}$, must be localized the same region. However, during the time evolution from the initial condition up to the matching time $\tau_{0}$, the wave packet $\tilde{D}_{\lambda}$, that moves backward in time, will proceed to the right because waves cannot propagate upstream after the sonic horizon, at least at long wavelengths. Analogously, $\tilde{D}^{*}_{p}$, moving forward in time, will be dragged inside the black hole. The only possibility to find both $\tilde{D}_{\lambda}$ and $\tilde{D}^{*}_{p}$ centered in the same region is to require that both group velocities are vanishingly small, i.e., that the momenta $p$ and $\chi^{\prime}_{\lambda}(z=0)$ are positive and small. For this choice of parameters we can assume that $D^{0}_{\lambda}(z,\tau_{0})$ is still well-approximated by its value long after the quench, i.e., by the expression (91) appropriate for the normal mode corresponding to the stationary state solution (51, 53, 54). For $\omega_{\lambda}\sim 0$, $z\sim 0$ and $\chi^{\prime}_{\lambda}(z)>0$ we get

[TABLE]

where

[TABLE]

is the Unruh analogue surface gravity (recall that in our model $v(z)$ is negative). The divergence of $\chi^{\prime}_{\lambda}(z)$ when $z\to 0$ implies that the leading contribution to the integrals defining $U_{p\lambda}$ and $V_{p\lambda}$ comes from such a region, where we can approximate $D^{0}_{p}(z,\tau_{0})\sim D^{0}_{p}(0,\tau_{0})$:

[TABLE]

which, after integration, give the transformation matrices

[TABLE]

with

[TABLE]

Recall that, since we are working in the Heisenberg picture, the ground state $|\Phi\rangle$ is unchanged during the time evolution. However, it is non-trivially related to the vacuum $|0\rangle$ of the phonon operators long after the quench, defined by $\hat{b}_{\lambda}|0\rangle=0$:

[TABLE]

where the symmetric matrix $F_{\lambda\lambda^{\prime}}$ is defined as the solution of the linear problem

[TABLE]

By substituting the previous expressions we finally get

[TABLE]

and we see that it mixes positive and negative frequencies. The ground state, in terms of the dressed phonon operators, has thus the form of a two-mode squeezed state. Notice that the negative frequencies are only defined in the supersonic region and cannot reach the region $z>0$. By tracing out them [38, 39] we then get a thermal density matrix with the Hawking temperature

[TABLE]

which is the same obtained in the exact solution of the model (34).

Let us notice that $|\Phi\rangle$ is written as a linear superposition of excited states including an infinite number of phonons of positive and negative frequencies. Since each pair has zero energy, this combination remains an eigenstate of the fluctuation Hamiltonian $\hat{H}_{2}$, which is essential to prove that the asymptotic state is a stationary solution of the evolution equations.

6. Experimental configurations

The exhaustive analysis we have presented refers to the specific case of Hard Core Bosons in one dimension, because only in this limit an analytic solution is available. Although this special system has been already experimentally reproduced several years ago [22, 23], the Tonks-Girardeau limit is rather difficult to obtain due to the tight one-dimensional confinement required. We believe that the strong coupling restriction may be relaxed and the Hawking emission may be seen also under less stringent conditions, as recent studies suggest. To this aim, we briefly examine the case of a gas of 87Rb atoms confined in a cylindrical trap of radius $a_{\perp}=0.25\,\mu$m, as an illustrative example. The trap is assumed to be flat in the longitudinal direction with total length $L\gtrsim 10\mu$m. Let us fix the initial number density of the condensate to $\rho_{0}=3.8\,10^{3}\,\mu$m*-1* and its initial velocity to $v_{0}\sim 18$ mm/s. Now we turn on an external potential of the form $V(x)=V_{0}\,e^{-(\alpha\,x)^{2}}$ with $V_{0}=\frac{\hbar^{2}Q^{2}}{2m}$ and $\alpha=0.1\,Q$. The potential amplitude corresponding to $Q=38\,\mu$m*-1* is $V_{0}\sim 3.6\,\mu$K. Then, we follow the semiclassical dynamics by integrating the usual (cubic) Gross-Pitaevskii equation in the $\alpha\to 0$ limit. As illustrated in Fig. 20 a sonic horizon forms in $x=0$ and the flow reaches a stationary state in a region of $10\,\mu$m around the horizon after a fraction of a millisecond.

According to Eq. (123) we thus expect the emission of Hawking radiation at the relatively high temperature $T_{H}\sim 100$ nK. Measuring the momentum distribution in the far upstream region or the density correlations across the horizon, the characteristic signatures of thermal emission should be clearly visible. In an experimental set-up it might be more convenient to maintain the condensate still and move the external barrier, as in the experiments of Ref. [6]. This choice is clearly equivalent, via Galileo transformation, to the case considered in the previous Sections.

7. Conclusions

The exact solution of a one-dimensional model of interacting bosons allows to unveil the microscopic origin of the phonon emission representing the condensed matter analogue of the Hawking radiation. Without referring to the gravitational analogy, we investigate the phonon flux which originates in the stationary state reached by a flowing BEC in the presence of an external potential. The fully analytical solution of the model provides several advantages over previous studies in this field: in fact, we unambiguously show that thermal emission can be expected only for a condensate flowing against a barrier potential which accelerates the particles to supersonic speeds and only in the limit of extremely smooth obstacle, i.e., when the semiclassical (WKB) limit is reached. Otherwise phonon emission does occur (even in the absence of a sonic horizon) but the thermal character is lost. This is the case of an external potential which has a step-like form. Moreover, we are able to follow the quench dynamics mimicking the experimental procedure, confirming that the analytical stationary solution is indeed reached at long times. Since density correlations across the horizon represent a widely adopted probe for the occurrence of the Hawking mechanism, we evaluate the correlation pattern; this shows the features expected on the basis of semiclassical arguments.

Although the exact solution contains all the required information about the model, we deem instructive to also develop a semiclassical study of the same set-up, in order to understand whether the currently adopted semiclassical paradigm, based on the use of the Gross-Pitaevskii equation and of the Bogoliubov analysis of the excitation spectrum, is able to reproduce the exact results. Surprisingly we find that, even in the presence of extremely smooth potentials, the semiclassical analysis lacks some important features: indeed, an entire class of stationary solutions, characterized by non-monotonic velocity profiles, cannot be obtained by means of semiclassical arguments.

Finally, we provide an illustrative example of a possible experimental set-up apt to show thermal phonon emission in a BEC at a reasonably high temperature.

The formalism introduced here is suited to further generalizations in order to address other phenomena often studied in the field, such as the black-hole laser effect or the dynamical Casimir effect. To this extent, innovative results could be developed in these directions also. Unfortunately, the extension of this model to higher dimensions is not trivial.

Appendix A

Here we report the spectrum of the non-interacting single-particle Hamiltonian

[TABLE]

with $V(x)=\frac{\hbar^{2}Q^{2}}{2m}\Theta(x)$. The eigenvalues are written in the free particle form $\epsilon_{k}=\frac{\hbar^{2}k^{2}}{2m}$, while the corresponding wave-functions acquire different forms for $0<k\leq Q$ or $|k|>Q$:

$i)$ $0<k\leq Q$. In this regime the spectrum is non-degenerate and the eigenfunctions read:

[TABLE]

where $\lambda_{k}=\sqrt{Q^{2}-k^{2}}$ and the reflection and transmission coefficients are

[TABLE]

These states are exponentially trapped in the region $x<0$. Note that the eigenfunctions $\phi_{k}(x)$ for $-Q<k\leq 0$ are not defined.

$ii)$ $|k|>Q$. Here the spectrum is doubly degenerate. The right moving wave for $k>Q$ is:

[TABLE]

and the left moving wave for $k<-Q$ is:

[TABLE]

where $q=\sqrt{k^{2}-Q^{2}}$ and the reflection and transmission coefficients are:

[TABLE]

for the right moving and

[TABLE]

for the left moving. The normalization condition between any pair of wavefunctions reads, as usual,

[TABLE]

Appendix B

We give some detail on the evaluation of the asymptotic form (7) starting from the exact time evolution (5). We first evaluate the internal products $\langle\phi_{p}|\psi^{0}_{k}\rangle$ (6) which, using the explicit expressions of Appendix A, are written as a sum of simple poles like

[TABLE]

where the Plemelj identity has been used. Inserting this sum into Eq. (5) we get:

[TABLE]

where, besides the term coming from the $\delta$-function, the integral gives also a non-vanishing contribution at long times. In fact, by defining $u=p-k$ and expanding $\alpha_{k+u}$ and $\epsilon_{k+u}$ to first order in $u$, it is straightforward to show that the second term has a finite limit as $t\to+\infty$:

[TABLE]

whose sign depends on the slope of the fermionic dispersion $\epsilon_{k}$ so that, for $\pm\epsilon_{k}^{\prime}>0$, it precisely cancels the contribution coming from the $\delta$-function. In conclusion, we can formally write:

[TABLE]

This analysis is easily extended to all external potentials which do not admit bound states, because the pole contribution in the internal product (6) comes from the asymptotic regions $x\to\pm\infty$, where every eigenfunction $\phi_{p}(x)$ acquires the form of a scattering state, i.e., can be written as a linear combination of an incident (or transmitted) and a reflected wave. Therefore, the structure of the eigenfunctions at $x\to\pm\infty$ has the same form shown in Appendix A, irrespective of the details of the external potential.

Appendix C

Exact results for the barrier (9) can be obtained starting from (10) and using that

[TABLE]

as well as the asymptotic expression

[TABLE]

when $z\to 1$ (since in our case Re$(c-a-b)=0$). These allow to find exact expressions for the asymptotic density $\rho(x)$ when $x\to\pm\infty$:

[TABLE]

Appendix D

According to the exact form (10) of the scattering states for a barrier of the form (9), the hypergeometric function $F(a,b;c;\zeta)$ has to be evaluated for arguments $(a,b,c)$ whose imaginary part grows large as $\alpha\to 0$. Here we derive the appropriate asymptotic limit of the hypergeometric function. First we use the identity (see Ref. [27], Eq. 15.3.3):

[TABLE]

Then we note that, in the case of interest (10), the variable $z$ is real and the parameters are related by $2c-a-b=1$. Therfore, we can express the hypergeometric function in terms of Legendre functions by (see Ref. [27] Eq. 15.4.17):

[TABLE]

Finally, recalling that $\zeta=\frac{1-\tanh(\alpha x)}{2}$ belongs to the interval $(0,1)$, we use the integral representation (see Ref. [40] Eq. 8.714-2) to obtain:

[TABLE]

Inserting the expressions for $a=\frac{1}{2}-i\frac{k+Q}{\alpha}$ and $b=\frac{1}{2}-i\frac{k-Q}{\alpha}$, the integral can be written as

[TABLE]

with

[TABLE]

These expressions hold for any value of $\alpha$. However, in the $\alpha\to 0$ limit they considerably simplify because the integral can be explicitly evaluated by use of the stationary phase method. Looking for the extrema of $\varphi(t)$ with $t>0$, we have to consider two distinct regimes:

• •

If $k>Q$, the unique extremum is given by

[TABLE]

• •

Instead, if $k<Q$, two solutions exist for $Q\,\tanh(\alpha x)<-\sqrt{Q^{2}-k^{2}}$:

[TABLE]

while no extremum is present for $Q\,\tanh(\alpha x)>-\sqrt{Q^{2}-k^{2}}$ ,

where, in both cases,

[TABLE]

The second derivative of the phase, evaluated at the extremum, is given by

[TABLE]

Therefore, each extremum (if any) contributes to the integral with a term

[TABLE]

The scattering states (10) then acquire an analytic form in the $\alpha\to 0$ limit. By use of the Stirling formula for the asymptotic behavior of the $\Gamma$ function, we finally obtain:

[TABLE]

where we used the shorthand notation $\varphi(t_{\pm})=\varphi_{\pm}$. The second term in curly brackets is present only for $k<Q$ and $Q\,\tanh(\alpha x)<-\sqrt{Q^{2}-k^{2}}$, while for $k<Q$ and $Q\,\tanh(\alpha x)>-\sqrt{Q^{2}-k^{2}}$ the wavefuncion vanishes to leading order in $\alpha$. This result has been derived for $k>0$. The analogous expression for $k<0$ is simply obtained by changing $k\to|k|$ and $x\to-x$ (or $\zeta\to 1-\zeta$).

As $\alpha\to 0$, the square modulus of the scattering wavefunction just reduces to

[TABLE]

for $k>Q$, while for $k<Q$ and $Q\,\tanh(\alpha x)<-\sqrt{Q^{2}-k^{2}}$ an oscillatory term still survives:

[TABLE]

These rapid oscillations in the single-particle scattering wavefunction are however washed out when we perform th integration over the wavevectors $k$ required to evaluate the averages in the Fermi gas (2).