Birth of a quasi-stationary black hole in an outcoupled Bose-Einstein condensate

We study the evolution of an initially confined atom condensate which is progressively outcoupled by gradually lowering the confining barrier on one side. The goal is to identify protocols that best lead to a quasi-stationary sonic black hole separating regions of subsonic and supersonic flow. An optical lattice is found to be more efficient than a single barrier in yielding a long-time stationary flow. This is best achieved if the final conduction band is broad and its minimum not much lower than the initial chemical potential. An optical lattice with a realistic Gaussian envelope yields similar results. We analytically prove and numerically check that, within a spatially coarse-grained description, the sonic horizon is bound to lie right at the envelope maximum. We derive an analytical formula for the Hawking temperature in that setup.


I. INTRODUCTION
An attractive feature of Bose-Einstein condensates is that of providing a convenient way of investigating analog black-hole physics in the laboratory [1,2]. It was already noted by Unruh [3,4] that Hawking radiation in a cosmic black hole [5,6] is an essentially kinematic effect that could be simulated in a quantum fluid. More specifically, it has been predicted that, for a quantum fluid passing through a sonic horizon (i.e., a subsonic-supersonic interface), phonons will be emitted into the subsonic region even at zero temperature [7][8][9][10][11][12].
A sonic black hole has been recently realized in an accelerated Bose-Einstein condensate [13]. An alternative route towards the detection of Hawking radiation may be provided by a quasi-stationary horizon, which in principle can be achieved by allowing a large confined condensate to emit in such a way that the coherent outgoing beam is dilute and fast enough to be supersonic [11,[14][15][16]. The hope is that the so far elusive Hawking radiation will be less difficult to detect unambiguously in such quasi-stationary transport scenarios. In particular, some recent works have addressed the possibility of detecting the spontaneous Hawking radiation above the stimulated signal [17][18][19].
The main goal of the present paper is to explore the actual attainability of the steady-state regime. Within a mean-field approximation, we investigate the dynamics of an initially confined condensate that begins to leak as the height of one of the confining barriers is driven from an essentially infinite to a finite value that permits a gentle yet appreciable flow of coherently outcoupled atoms. More specifically, we will focus on the case where the increasingly transparent potential is formed not by a single [16] or double [15] barrier, but by an extended optical lattice, the main reason being that the latter scenario seems more suitable for the achievement of quasistationary flow. We will see that close-to-ideal stationary flow within the permitted energy bands is achievable un-der realistic opening protocols.
The present mean-field study aims at identifying transport scenarios that offer hopes for a future detection of Hawking radiation. More conclusive predictions about the detectability of spontaneous radiation will require a study of the time-dependent Bogoliubov -de Gennes (BdG) equations, which should inform us on the actual intensity of the expected spontaneous radiation. This task is left for a future study.
Besides the motivation of realizing gravitational analogs, the achievement of stationary transport scenarios is of general interest for the investigation of atom quantum transport, in the case of both bosons [20][21][22][23][24][25] and fermions [26,27], within the emergent field of atomtronics.
This paper is arranged as follows. Section II presents the model for the gradual reduction of the optical lattice amplitude which we will be investigating. After some preliminary remarks in Section III, the main numerical results together with some theoretical arguments (that help to understand the observed trends) are presented in Section IV. The second part of that section describes the achieved quasi-stationary regime. Section V addresses the more realistic case of an optical lattice having a Gaussian envelope. Interestingly, we find that the horizon lies at the maximum of the envelope and give a theoretical explanation of that fact. The main conclusions are summarized in section VI. Appendix A provides a detailed description of the initial state of the condensate as it exists before the deconfinement procedure begins. Appendix B discusses some properties of Bloch waves in the presence of nonlinearities accounting for the interaction. Appendix C presents a perturbative treatment of the interaction. Finally, Appendix D describes the numerical method of integration and the use of absorbing boundary conditions.

II. THE MODEL
In this work we study the outcoupling of a onedimensional (1D) Bose-Einstein condensate through a finite-size repulsive optical lattice, whose intensity is gradually lowered in such a way that, within a finite time, the periodic barrier shifts from a regime of practical confinement to one of full transparency within certain atom energy bands. We focus on the mean-field dynamics, i.e., we only consider the evolution of the condensate wave function, leaving the dynamics of quasi-particles for a future study. We restrict our present study to a quasione dimensional model. The time-dependent condensate wave function Ψ(x, t) obeys the Gross-Pitaevskii (GP) equation [28,29]: where m is the atom mass and V (x, t) the time-dependent optical lattice potential. The effective one-dimensional coupling constant g = 2 ω tr a s (with a s the s-wave scattering length) is the relevant interaction strength in a setup where only the ground state of a confining transverse harmonic oscillator of frequency ω tr /2π is populated. This is the 1D mean-field regime, characterized by the condition ρa s ≪ 1, where ρ(x, t) = |Ψ(x, t)| 2 [30,31]. At the same time, ρa 2 tr /a s ≫ 1 (with a tr the transverse oscillator length) must be satisfied to stay away from the Tonks-Girardeau regime [30,32,33]. Taking the initial bulk density n 0 as a typical value for the density, we can realistically set (see Sec. IV) n 0 a s ∼ 10 −1 and n 0 a 2 tr /a s ∼ 10 3 , from which we conclude that we are safely in the 1D mean-field regime.
Equation (1) conserves the total particle number N as given by the normalization condition The condensate density is nonzero only for x > 0 because at all times we assume a sufficiently high barrier at x = 0, which is simply implemented via the hard-wall boundary condition Ψ(0, t) = 0. Initially (at times t < 0), we consider an equilibrium condensate made of N atoms occupying the region 0 < x L. Thus n 0 ≃ N/L is the initial atom density, which is defined below more precisely. We also introduce an optical lattice that spans the region L x L + L lat and whose initial amplitude V 0 is large enough for particle tunneling through the lattice to be practically forbidden. The initial wave function is stationary, Ψ(x, t) = e −iµ0t/ Ψ(x), with Ψ(x) satisfying the time-independent GP equation The initial chemical potential µ 0 is determined by the normalization condition (2). The initial healing length is defined as ξ 0 ≡ 2 /mgn 0 , where n 0 ≡ µ 0 /g. Further details on the initial condensate are given in Appendix A. At time t = 0, the optical lattice intensity starts to decrease and atoms begin to escape towards the region x L + L lat , where the potential is assumed to be negligible. On quite general grounds [15,31,34], the flow beyond the optical lattice can be expected to be supersonic.
We assume that the optical lattice is made of two fixed phase lasers of wavelength λ and whose wave vectors form an angle θ [35,36]. The time-dependent optical lattice potential is chosen so that in the lattice region (defined by where k L = π/d and d = λ/ [2 sin(θ/2)] is the lattice period, while V (x, t) = 0 everywhere else. The potential profile in Eq. (4) is somewhat idealized. A more realistic choice should include a Gaussian envelope. For simplicity, we choose to start by considering a flat-envelope optical lattice, where Bloch's theorem can be invoked with reasonable confidence. We will see that, remarkably, essentially the same results are obtained when a more realistic Gaussian envelope is used. A sketch of the time-dependent, flat-envelope optical potential and the resulting condensate flow is presented in Fig. 1.

III. PRELIMINARY REMARKS.
We note that the time-dependent amplitude V (t) evolves from V 0 ≫ µ 0 at t ≤ 0 to V ∞ µ 0 for t ≫ τ . The asymptotic behavior is determined by the initial parameters of the condensate (N, g, L), the specific form of the final potential (V ∞ , d), and the barrier lowering time scale (τ ). The initial potential amplitude, V 0 , plays almost no role provided that it is sufficiently large. More precisely, the condition of initial confinement requires µ 0 to lie well below the lowest conduction band of the initial optical lattice potential. On the other hand, the properties of the final steady state are insensitive to τ unless τ is very small (see Subsection IV A 2).
The main goal of the present work is to identify the barrier-lowering protocol that best leads to a regime of quasi-stationary outcoupled flow, by which we mean a flow regime characterized by parameters that vary slowly in time in a sense that will be specified later. As to the space dependence in that regime, we require the density to be as uniform as possible in the region 0 < x L.
In the supersonic region (x L + L lat ) we also want a uniform flow profile, even though, due to the low density, this will be more difficult to achieve. However, disturbances in the supersonic region should not affect the FIG. 1: Schematic representation of the emitting condensate setup studied in this article. Within the ideal lattice scenario, hard-wall boundary conditions are assumed at x = 0 and an optical lattice lies in the region L < x < L lat with a time-dependent amplitude such that the potential V (x, t) (represented by the semi-transparent yellow surface over the x − t plane) evolves from strongly to moderately confining. The resulting time-dependent density profile n(x, t) is represented by the grey-blue surface. The vertical axis gives the density in some (here unimportant) units. The surface V (x, t) is uplifted to provide a better vision of n(x, t). Some parameters defined in the main text are indicated. The trend towards a long-time quasi-stationary flow regime can be qualitatively appreciated.
spectrum of the Hawking radiation, which is emitted into the subsonic region. On the other hand, in the optical lattice region, the flow should be as close as possible to that of a propagating Bloch wave. At the boundary between these two regions, large gradients of the flow speed and density are likely to occur, but the current density should remain essentially uniform.
In what follows, when we refer to bands, we will be meaning the Schrödinger (non-interacting) bands, unless specified otherwise. Bands in a nonlinear context are discussed in Appendices B and C. We will see that, in virtually all cases, the relevant properties of the optical lattice are determined by its lowest band, provided the lattice is sufficiently long. An important result is that, due to the finite size of the subsonic reservoir, the quasi-stationary flow is formed only when the space-averaged chemical potential lands at a value slightly above the bottom of the lowest lattice band. Because the local chemical potential is approximately uniform along the structure, this implies that the density is small almost everywhere in the optical lattice region. As a consequence, the interaction term can be neglected, since g|Ψ(x, t)| 2 ≪ ( k L ) 2 /m; see Appendix B for details. In the non-interacting regime, Eq.
(1) becomes the usual Schrödinger equation, which for a sinusoidal potential can be transformed into a Mathieu's equation [37].
The structure of the final bands can be characterized by the dimensionless parameter v ≡ mV ∞ /8 2 k 2 L .
The nearly-free atom regime occurs when v ≪ 1. Then bands are wide and gaps are narrow. By contrast, in the tight-binding regime (v ≫ 1), bands are narrow and widely spaced. Since v ∝ V ∞ d 2 , the band structure can be modified by changing the lattice amplitude or its spacing.

IV. IDEAL OPTICAL LATTICE
In this work, the unit length is the initial bulk healing length ξ 0 defined in Section II. Accordingly, velocities are measured in units of the sound speed, c 0 ≡ gn 0 /m, and times in units of t 0 ≡ ξ 0 /c 0 . Energies are expressed in units of the initial chemical potential µ 0 = mc 2 0 . Quasione-dimensional condensates of 87 Rb are typically made of N ∼ 10 4 − 10 7 atoms and have a transverse trapping frequency ω tr ∼ 2π × 10 3 Hz and a confinement length L ∼ 10 − 400µm. The optical lattice periodicity is bounded from below, d > λ/2nm, for geometrical reasons. The value of lambda is chosen to be sufficiently far from the resonance to avoid any spontaneous emission and on the blue-side of the resonance to produce a repulsive potential (lambda < 780nm for rubidium atoms). The simulations are run up to times t ∼ 10 4 − 10 5 t 0 . As t 0 = (2n 0 a s ω tr ) −1 ∼ 10 −4 s, then t ∼ 1 − 10 s for our simulations, which is on the order of the mean lifetime of this type of condensates.
In the simulations we use a numerical scheme based on the Crank-Nicolson method to integrate the timedependent GP equation (1). Hard wall boundary conditions are assumed at x = 0. At the other end of the computational grid (located at x = L t ), we use absorbing boundary conditions. Further details of the numerical method are given in Appendix D.

A. Analysis of the simulations
To characterize the quasi-stationary regime, we use the local chemical potential defined as which can be complex. For a stationary solution, Ψ(x, t) = e −iµt/ Ψ(x), one has µ(x, t) = µ, real and independent of (x, t). The current, is also independent of (x, t) for a stationary solution, as dictated by the continuity equation ∂ t ρ + ∂ x j = 0. In the quasi-stationary regime the uniformity of j(x, t) is impossible to fulfill strictly, because the current is zero at x = 0 while the emitted atoms carry a nonzero current. Thus, there must be a current gradient and, by the continuity equation, the density has to be time dependent. The hope is however that, in the quasi-stationary regime, this time dependence is weak because the condensate leak is slow.
On the other hand, we can expect that, in the quasistationary regime, µ(x, t) is a sufficiently uniform function, with small spatial variations around its spaceaveraged mean value. To check this expectation in a quantitative manner, we introduce the space-averaged chemical potentialμ(t) together with an appropriate measure of its relative spatial fluctuation spread σ(t): We recall that µ(x, t) andμ(t) can be complex. A nonzero imaginary part of µ(x, t) reflects a leaking condensate, as revealed by the local relation Accordingly, we can define and compute the emission rate per particle as where the continuity equation has been used. The spatial average of the time-dependent chemical potential [see Eq. (8)] is mostly determined by the subsonic region, where µ(x, t) ≃ gρ(x, t).
A further rescaling of the condensate wave function Ψ(x, t) → √ n 0 Ψ(x, t) reveals more neatly the intrinsic parameters governing the system. Once the healing length ξ 0 is given, only L/ξ 0 , d/ξ 0 , τ /t 0 , V ∞ /mc 2 0 and n osc ≡ L lat /d (number of oscillations in the optical lattice) are relevant for the problem. We have already noted that V 0 plays almost no role in the limit V 0 ≫ µ 0 . We find that, for the pertinent experimental ranges, namely, L ∼ 10 − 400 µm, variations of L/ξ 0 have little effect on the properties of the quasi-stationary regime. We have noted that they have a small influence on the time needed to achieve the desired quasi-stationarity, which grows weakly with the initial size of the condensate. A similar point can be made about n osc , which becomes unimportant when it lies in the range n osc ∼ 15 − 100. The (relatively small) effect of increasing n osc even further is that there are more spatial fluctuations in the chemical potential, for two reasons: (a) the optical lattice tends to host larger µ(x, t) spatial fluctuations than the subsonic zone because of atomic reflections across the wells, an effect that is enhanced for larger optical lattices; (b) the larger the lattice, the bigger its contribution to the average chemical potential and its fluctuations.
In summary, except for the above remarks, only the parameters d/ξ 0 , τ /t 0 , V ∞ /mc 2 0 , have a noticeable effect on the transport properties of the system under study.

Role of the final band structure
As noted before, the combination of d and V ∞ fixes the final band structure. Figure 2 shows the various scenarios which one may find depending on the long-time width and position of the lowest band with respect to the initial chemical potential, µ 0 . The band structure is computed numerically. The desired steadiness of the long time behavior improves with the width of the band, as the first row in Fig. 2 reveals. In Fig. 2a, a favorable case (wide band) is presented and compared with a less favorable case in Fig. 2b, which has the same conduction band minimum but a narrower band. After a short transient, a comparison of the relative chemical-potential standard deviation σ(t), as defined in Eq. (8) and plotted in this graph, shows a clear advantage in the use of wider bands. For instance, in Fig. 2a, σ(t) ∼ 10 −4 in the stationary (long time) regime, about 10 times smaller than in Fig.  2b.
Besides, after a transition time of order ∼ 5000t 0 , all the characteristic magnitudes of the system shown in Fig.  2a vary slowly enough in time to properly view the resulting flow regime as quasi-stationary. In fact, the leak is so slow that other processes which limit the lifetime of the condensate (such as condensate decay due to inelastic collisions) operate on a shorter time scale.
If the chemical potential reaches and goes below the bottom of the conducting band in a relatively short time, then a transition occurs to an essentially confined situation where the leaking is exponentially small, corresponding to an atom transmission probability T (L lat ) ∝ exp(−κL lat ), where κ ∝ √ E min − µ, with E min the bottom of the conducting band. This situation whereby one soon reaches the regime µ < E min is not interesting for our purposes because we need some appreciable flux in order to form a useful black-hole configuration. In particular, we are typically interested in considering condensates so large that the time needed to reach the bottom of the conduction band is longer than the typical lifetime of the condensate.
A further argument can be invoked in favor of wide bands. In virtually all the cases we have addressed, Reμ(t) drops until it almost reaches the bottom of the conduction band, where leaking is slow. In the resulting regime, the density in the lattice is very small, gn r ≪ 2 k 2 L /m, wheren r is the mean density in the optical lattice, as defined precisely in subsection IV B. It is shown in Appendix C that, in this regime of low interactions, the width of the conduction band for the linear perturbations (whose evolution is governed by the BdG equations) is very close to that obtained for the linear Schrödinger equation, and the corrections are there computed. This means that the optical lattice acts like a low-pass filter, the band width being the equivalent of the cutoff frequency. The higher the cutoff, the wider is the transmission band of the lattice. As a consequence, fluctuations on the subsonic side are transmitted away through the optical lattice, which reduces the space fluctuations in the chemical potential. Another trend can be appreciated in the second row of Fig. 2. When placing µ 0 slightly below the top of the conduction band or in the first gap, as Fig. 2d illustrates, the leaking occurs faster but the quasi-stationary regime [small σ(t)] is not reached within the time span of the simulation. For the purpose of keeping σ(t) ≪ 1, a more favorable situation for the chemical potential is shown in Fig. 2c. There, for the same length d as in Fig. 2d but a higher barrier amplitude V ∞ , the chemical potential is initially placed in the final conduction band. This case clearly yields smaller fluctuations, even though the width of the conduction band is smaller. This shows that, besides having wide bands, one also needs that µ 0 be placed close to the bottom of the final conduction band in order to obtain a more favorable quasi-stationary regime. This point is further discussed in the next subsection (IV B).
Within the nearly-free particle approximation (v ≪ 1), the bottom and top of the first conduction band are given by the relations: where E R ≡ 2 k 2 L /2m is the recoil energy of the optical lattice. Given that k L = π/d, the condition that the initial chemical potential lies within the final conduction band, i.e., is guaranteed to be satisfied if The left inequality is just while the right inequality can be rewritten as: Equations (14), (15) express a sufficient condition to satisfy Eq. (12).

Non-adiabatic effects
In the favorable situation shown in Fig. 2a, the dependence on τ is not important as long as τ ≫ t 0 . A simulation is presented in Fig. 3 which shows that, in the fast regime (τ ∼ t 0 ), and due to the high-frequency excitations induced by the short lowering time scale, σ(t) remains higher than in the adiabatic case (see Fig. 2a).

B. Quasi-stationary regime
From the discussion in the previous subsection (IV A), and particularly through the situation shown in Fig. 2a, we have learned that, for wide enough bands, initial chemical potential close to the bottom of the final conduction band, adiabatic evolution (τ ≫ t 0 ), and a typical setup, the system evolves towards a quasi-stationary regime in times shorter than the lifetime of a condensate.
The quasi-stationary regime can be defined as that in which Reµ(x, t) is essentially uniform [σ(t) ≪ 1] and its global time variations take place on a time scale of the order of or greater than the condensate lifetime. The achievement of this regime is of general interest as a scenario for the study of atom quantum transport. In particular, one may expect spontaneous Hawking radiation to be detectable above a background of sufficiently small spatial fluctuations. In the present section we discuss several features of the condensate wave function for the quasi-stationary regime. For illustration purposes, all the graphs considered in this subsection have been obtained for a system with the parameters of Fig. 2a, which are sufficiently representative.
By writing the condensate wave function as Ψ( v(x, t) being the local condensate flow velocity and c(x, t) the local speed of sound. The spatial variations of both velocities are small in the subsonic and supersonic regions, but not in the lattice. We note that, in that region, c(x, t) must not be regarded as the lattice sound speed; see Appendices B and C. The profile of both quantities computed at a time, t = 4 × 10 4 t 0 , after a barrier removal time of τ = 500t 0 , is shown in Fig. 4. The subsonic zone shows an essentially flat (uniform) density and flow speed profile in the sense that the spatial fluctuations are on the order of ∼ 10 −4 n 0 for the density and ∼ 10 −3 c 0 for the flow speed, too small to be appreciated in Fig. 4. In Appendix A, an approximate analytical formula [Eq. (A21)] for the wave function of the confined condensate which fits the numerical results within this level of accuracy. This good agreement reflects the low value of the flow velocity in the condensate region.
On the other hand, in the deep central region of the optical lattice, Bloch's theorem is satisfied. We introduce the space-averaged densityn r (t) by averaging ρ(x, t) over the optical lattice after excluding 10 lattice sites at each end of the lattice. That average density, combined with numerically computed quantities that depend on the optical lattice potential, yields an effective sound velocity [see Eq. (C15)] that is plotted as a horizontal green segment spanning the averaged region in Fig. 4. It can be clearly seen that, within the optical lattice, the flow is subsonic, the horizon lying on its right edge. In the quasi-stationary regime,n r (t) decreases at a rate comparable to the inverse lifetime of the condensate, as the inset in Fig. 5 shows.
At the edges of the lattice there are strong variations of the density due to the the matching between the vastly different densities found on both the subsonic and the supersonic side.
A check on the approximate validity of Bloch's theorem in the presence of non-linear corrections, for the central part of the optical lattice and in the quasi-stationary regime, is also shown in Fig. 5. In this graph, the real part ofμ(t) is compared with the time-dependent chemical potential computed usingn r (t) and assuming zero Bloch momentum, as explained in the first paragraphs of Appendix B. The good agreement between the two curves suggests that the condensate is flowing with a very small Bloch momentum.
In the supersonic zone, once in the quasi-stationary regime, both density and flow speed profiles are almost uniform. This is hinted at in Fig. 4 but not shown explicitly. Part of the non-flat behavior is due to small spurious reflections; see Appendix D for details. Conservation of the chemical potential and the near absence of interaction effects in this zone make the flow speed almost uniform, because the chemical potential is almost fully transformed into kinetic energy. On the other hand, the supersonic density decays with time as the reservoir is depleted, but the process is such that, at each instant, the density profile remains essentially uniform (not shown). An inhomogeneous density profile would show up only on a space scale much larger than that used in this simulation. A comparison of the decaying supersonic density n d (t) and the nearly constant supersonic side speed, v d (t), is shown in Fig. 6 (subindex d stands for 'downstream' region). The emission rate per particle (not shown), as computed from Eq. (10), is practically identical to the product n d (t)v d (t), except for a numerical factor corresponding to the instantaneous total number of particles, which in the quasi-stationary regime is practically constant. This emission rate gives us the typical time scale for the variation of the number of particles of the system, which is approximately the time scale for the variation of µ. In the quasi-stationary regime considered here, it is ∼ 10 −6 − 10 −7 t −1 0 , from which we infer that the typical variation time of the chemical potential is ∼ 10 6 − 10 7 t 0 , much longer than the lifetime of a condensate ∼ 10 4 t 0 .
We conclude that the realization of this quasistationary regime needs two fundamental ingredients: the existence of a band structure and the presence of interactions. Without a band structure as that provided by the optical lattice, the condensate would continue leaking through the barrier at a fast rate. On the other hand, the presence of interactions (as reflected in the fact that ∂µ/∂n = 0) allows the condensate to stabilize its flow near the bottom of the conduction band. If the interactions in the subsonic region were negligible, the condensate would empty quickly (if µ 0 lied in the final conducting band) or it would remain confined (if µ 0 lied in the final gap).
Such trends can be appreciated in the accompanying videos that represent simulations for the same parameters as in Fig. 2a except for L = 10 µm and N = 250. The qualitative conclusions are similar. In Video 1, we see the time evolution of the density of a condensate confined by an ideal optical lattice of 30 barriers. We see that the system achieves the desired quasi-stationary regime.
On the other hand, in Video 2, we introduce a similar potential but with just a single barrier. We observe that the fluid leaks faster through the single barrier because there is no band structure providing an uplifted bottom of the conduction band to efficiently slow down the density decrease. The interaction plays the additional role of providing relaxation channels whereby the condensate lowers its energy while some collective modes are excited. The existence of Landau instabilities (see Appendix B) when µ 0 lies well above E min can be clearly appreciated in the upper right corner of Fig. 5 of Ref. 38, whose chosen parameters are similar to those of the present work. The low value of the critical velocity helps to understand the small value of the condensate Bloch momentum which we infer from the numerical results shown Fig. 5 of our present work. The appearance of instabilities can also be viewed as responsible for the fast lowering of the chemical potential after being initially prepared above the final conduction band, as shown in Fig. 2d. This interpretation is consistent with the relatively large values found for σ(t) when µ 0 is considerably above E min .

V. GAUSSIAN-SHAPED OPTICAL LATTICE
Here we perform the same analysis as in the previous section but using a more realistic optical lattice which includes a Gaussian envelope [35,39,40]: wherew = w/ cos(θ/2) (with w the laser beam width and θ the angle between the laser beams) plays the role of an effective lattice length. The time dependence of V (t) is the same as in Eq. (4). Usually,w varies in a range 10 − 200 µm. Here, L is the position of the maximum of the lattice Gaussian envelope. For consistency, we replace the hard wall at x = 0 by a Gaussian barrier of the type V L (x) = U exp(−2x 2 /w 2 L ) with w L = 2 µm and U ≫ µ 0 in order to simulate a more realistic confinement on the left side. This time-independent potential must be added to the time-dependent potential (17) which provides confinement on the right; see Appendix A for a detailed description of the initial confinement.
In the "adiabatic" regime (w ≫ d), the solutions of the linear Schrödinger equation for this type of potentials can show features similar to those found for an ideal optical lattice with the same instantaneous amplitude V (t), as can be appreciated in Fig. 7, where the transmission bands are compared. If we focus on the long-time limit (V (t) = V ∞ ), the realistic potential acquires the form where is a slowly varying function. Then, we have a locally ideal optical lattice at each point of the space with amplitude V ∞ (x). Bloch's theorem can also be applied locally and a local band structure results which is plotted as a function of space in Fig. 8. The left panel presents the setup whose single atom transmission is plotted in Fig. 7. Since the bottom of the lowest lattice conduction band is an increasing function of the periodic potential amplitude, the bottleneck for transmission across the realistic lattice occurs at the center of its Gaussian envelope. This fact explains the accurate coincidence between the bottom of both conduction bands shown in Fig. 7. We also see that, for E > E R [defined after Eq. (11)], the particle encounters a gap somewhere along the Gaussian lattice, and this explains why in Fig. 7 the transmission begins to decay for E > E R . For E min (v) < E < E R , the setup shows a plateau of essentially perfect atom transmission. The absence of interference oscillations in this region is due to the adiabatic variation of the lattice envelope. The right panel presents the different case of E R < E min (v). From the foregoing arguments, we expect not to find a conduction band, as can be numerically confirmed. We conclude that, in order to have a well defined conduction band for the realistic lattice, the condition E R > E min (v) is required, which implies V ∞ < 2.33 E R . Combining all these considerations, the conclusion is reached that a necessary condition for achieving a quasi-stationary regime is E min (v) < µ 0 . We also require µ 0 < E R to avoid having µ 0 lying too high above E min (v), which, as found for the ideal lattice, tends to generate relatively high values of σ(t). This inequality is equivalent to Eq. (15).
Here space can also be divided into three zones. In the quasi-stationary regime, both the subsonic and supersonic zones are located where the Gaussian envelope amplitude is negligible compared to the chemical potential, i.e., where [see Eqs. (8) and (19)]. In order for the subsonic side to be well differentiated, we set L ≫w. The optical lattice region is the complementary of the subsonic and supersonic zones, i.e., the region where (20) does not apply. The requirements of quasi-stationarity are similar to those formulated for the ideal optical lattice. Specifically, the quasi-stationary regime requires broad conduction bands, an initial chemical potential close to the bottom of the final conduction band, and a barrier amplitude that evolves not very fast. We also find that the condensate leaks relatively fast until Reμ(t) [Eq. (8)] approaches the bottom of the conduction band. All these features can be appreciated in Fig. 9, which is the Gaussian-envelope equivalent of Figs. 2-3. We reach a quasi-stationary state in which σ(t) ∼ 10 −4 . The bands in Fig. 9 are computed like in the ideal case, assuming a uniform barrier amplitude V (t). As noted when discussing Fig. 7, the positions of the bottom of the ideal and the realistic conduction (or transmission) bands are very similar, so the lower threshold of the transmission band can still be a good reference value to discuss the evolution ofμ(t). We notice that in the Gaussian case, the condensate leaks from the beginning of the simulation, as expected for an optical lattice with a smooth envelope. This contrasts with Fig. 2, where the condensate only begins to leak when the chemical potential is placed within the conduction band. We also study the corresponding quasi-stationary state. For that purpose, we take a snapshot of the configuration at t = 4 × 10 4 t 0 for the parameters in Fig. 9. We compare the profiles of c(x, t) and v(x, t) in Fig. 10, which is the realistic equivalent of Fig. 4. We find again essentially flat profiles for the density and flow velocity in the supersonic region, with their time evolution shown in Fig. 11. The general features of this quasi-stationary configuration are similar to those of the ideal case, but some interesting new features appear.

A. Location of the sonic horizon and related properties
We notice in Fig. 10 that the horizon seems to be placed at the maximum of the Gaussian envelope. Actually, this can be explained on quite general grounds by invoking the properties of the quasi-stationary regime and the adiabaticity conditionw ≫ d, which allows us to think in terms of a local band structure stemming from a periodic potential of local amplitude V ∞ (x). Then we can use an adiabatically space-dependent version of Eqs. (C10)-(C18) (where the sound speed, atom current and chemical potential are obtained for an infinite optical lattice) by making every parameter slowly dependent on x.
In particular, we take: where we neglect the time dependence because the system is assumed to be already in the quasi-stationary regime. The local averages for n r ,v are taken over several lattice periods. The quantities E min (x), m * (x), α 0 (x) depend on x through the amplitude of the Gaussian envelope, V ∞ (x), and therefore go through a maximum or minimum when V ∞ (x) goes through a maximum. In the quasi-stationary regime, the chemical potential is already close to the bottom of the conduction band, so the perturbative study used in Appendix C is valid. Taking spatial derivatives at the envelope maximum (x = L), while noting that the chemical potential is almost uniform, µ(x) ≃μ, and that ∂ x j(x) can be neglected (as implied by the continuity equation and quasi-stationarity), we arrive at the numerically observed result: and, as a corollary, s ′ (L) = −v ′ (L)/2. We can exploit this fact to obtain more results. For example, we can obtain the value of the density and the current at x = L as a function ofμ: which are very good approximations to the actual numerical values. Here the dependence on L of the various parameters is understood. Using (23) we can arrive at a differential equation for the time evolution ofμ. First we note that, from the continuity equation, we can write: where N L is the number of particles contained between x = 0 and x = L. As the subsonic region is in the Thomas-Fermi regime (see Appendix A), we can takeμ ≃ gN sb /L sb , where N sb is the number of particles in the subsonic region and L sb its size. As the density in the optical lattice is small, we can assume N L ≃ N sb (which impliesμ ∝ N L ) and write Eq. (24): where C is a positive constant independent ofμ. The solution of this equation is (with t 1 an integration constant), which fits the numerical data of Fig. (9) reasonably well. Finally, we can estimate the value of the Hawking temperature, which is given by: If we note that we operate in the nearly-free atom approximation (v ≪ 1) and in the weak interaction regime (gn r ≪ E R ), and derive twice the third equation (21), we obtain which gives a good estimate of the numerical value of the Hawking temperature. Noting that V ∞ ∼ µ 0 , m * ∼ m, we obtain k B T H ∼ ξ 0 µ 0 /w ∼ 10 −2 µ 0 ≪ µ 0 . The temperature of the condensate is typically of the order of µ 0 /k B , so we conclude T H ∼ 10 −2 T ≪ T .
To observe the birth of the black hole and to check that the horizon position naturally evolves towards the maximum of the optical lattice envelope, we have created a movie (Video 3) that shows the time evolution of the coarse-grained velocities (c,v) of the emitting condensate using the setup parameters of Fig. 9. At long times the predicted coincidence between the sonic horizon and the maximum of the Gaussian envelope in the stationary regime can be clearly appreciated.

VI. CONCLUSIONS
Within a mean-field description, we have investigated the process whereby an initially confined atom condensate is coherently outcoupled as the barrier on one side is gradually lowered. The goal has been to identify the barrier-lowering protocol which best leads to a quasistationary sonic black hole located at the interface between subsonic and supersonic flow. We find that the use of an optical lattice for the lowered barrier is essential to achieve a regime of quasi-stationary flow with minimal value of the fluctuation spread. First we have focused on an optical lattice of finite length and uniform amplitude. We find that the long-time band structure of the optical lattice greatly influences the asymptotic behavior of the emitted atom flow. The best quasi-stationary flow is achieved when the lowest conduction band is broad and the initial chemical potential lies not much high above the bottom of the final conduction band. When we replace the uniform amplitude of the optical lattice by a more realistic Gaussian envelope, we find that the results are similar to those of a uniform lattice with the amplitude of the envelope maximum. Quite interestingly, we argue analytically and check numerically that, in the quasistationary regime, the horizon separating the regions of subsonic and supersonic flow is pinned down right at the Gaussian maximum. We also find that the Gaussian envelope is quite efficient in guaranteeing a small deviation from the ideal stationary flow.
Whether the quasi-stationary regimes here identified can become scenarios for the detection of Hawking radiation, is something that will have to be confirmed by a future study of the quasiparticle dynamics operating against the background of seemingly favorable mean-field configurations.
In this appendix, we compute the initial profile of the condensate, which at early times (t < 0) experiences a confining time-independent potential. We require a hardwall boundary condition at x = 0, which implies, via continuity equation, that the phase of the condensate is constant in space. The amplitude A(x) ≡ |Ψ(x)| of the solution to the time-independent Gross-Pitaevskii (GP) equation Eq. (3) in the region where there is no potential, satisfies the equation As is well known, this equation can be interpreted as the equation of motion for a particle with "coordinate" A and "time" x in a certain potential Invoking "energy" conservation, the equation can be integrated as where E A is the total energy of this effective motion. Following Ref. 15, Eq. (A3) can be rewritten in terms of ρ(x) = A 2 (x) as where and e 2,3 are the zeros of The solution can be expressed in terms of elliptic functions [37]. Imposing the boundary condition ρ(0) = 0, one obtains a solution of the form In order to determine e 2,3 , another boundary condition is needed, together with the particle number normaliza-tion´dx ρ(x) = N . From (A2) and (A6), the chemical potential can be written as (A8)

Ideal confinement
The ideal confinement boundary condition is defined as A(L) = 0, and the condensate is confined between 0 and L. Using Eq. (A7) we find where K(ν) is the complete elliptic integral of the first kind [41].
Hereafter we work with the ground state (n = 1). The particle number normalization is By performing the integral in Eq. (A10) and using (A9) we have where E(ν) is the complete elliptic integral of the second kind [41]. Equations (A9) and (A11) lead to where the healing length ξ ≡ 2 L/mgN is not identical to ξ 0 defined in section II.
After eventually solving for ν, e 2 and e 3 , Eq. (A7) can be rewritten as For the chemical potential, we obtain, using (A8)-(A9) Taking into account that ν is a function of L/ξ, as given by Eq. (A13), we plot Eq. (A15) in Fig. 12. In order to find ν, Eq. (A13) must be solved numerically. However, good approximate solutions can be found. We can clearly distinguish two different regimes: L ≪ ξ and L ≫ ξ. The physical interpretation is straightforward because the ratio between the kinetic energy and the interaction energy is Then, L ≪ ξ is the Schrödinger limit in which we have ν ≃ 0. In that limit we arrive at the well-known result √g e 3 L = π and ρ(x) = e 2 sin 2 (πx/L). In all the cases considered in this work, L ≫ ξ, so we work in the limit in which interactions represent the main contribution to the chemical potential (Thomas-Fermi regime). K(ν) diverges when ν → 1 while E(ν) remains finite. From (A13) this means that ν ≃ 1. In fact, there are cases in which 1 − ν is so small that it falls below computers floating-point relative accuracy. In those cases, the only way to obtain the solution is through asymptotic expansion. One can prove that, in that limit [41], Thus, Eq. (A13) is rewritten as and from its solution we get 2K = 1 + √ 1 + r 2 . Therefore, 1 − ν = 16e −2K ≪ 1 which implies both e 2 ≃ e 3 and µ 0 ≃ ge 2 . Equation (A9) implies and then (A20) Collecting all these results, the wave function can be effectively approximated by . The function tanh quickly reaches the asymptotic value 1, which means that in the central zone of the confinement region, the solution is essentially flat.

Ideal optical lattice potential
For computational purposes, the initial rightmost boundary condition is also taken A(L) = 0 but now half a period of the lattice potential lies inside the confinement region, as explained in the main text. This artificial boundary condition does not create a problem because we take V 0 ≫ gN/L so the function inside the lattice potential is exponentially small. In order to compute the stationary solution in the region where the potential is present, a numerical solution of the GP equation has to be performed. In the situations considered in the present work, L ≫ d, so the wave function is very similar to that of the ideal confinement case. For numerical convenience, instead of fixing N and then obtaining the chemical potential, we first set n 0 to a typical experimental value of the density. Then, using µ 0 = gn 0 , we compute the number of particles N by integrating the resultant GP wave function. The computed number of particles satis-

Realistic optical lattice potential
In order to simulate a more realistic scenario, we introduce the following two potentials: on the left side, a Gaussian barrier centered at x = 0 of the form V L (x) = U exp(−2x 2 /w 2 L ), with w L = 2 µm, and on the right side, a realistic optical lattice centered at x = L, which has the We take the amplitudes of the confining potentials much larger than the chemical potential. We also take L ≫w ≫ w L , so that they are well separated in space and there is a large region where the potential is negligible and where we expect some kind of flat wave function. The width w L does not play a significant role in our simulations; we choose w L = 2 µm. In this way, we can set as boundary conditions for the numerical computation A(0) = 0 and A(L bc ) = 0, with L bc sufficiently deep in the region where V 0 exp[−2(x − L) 2 /w 2 ] ≥ µ 0 . Once we have fixed the potential and the boundary conditions for the numerical calculation, we repeat the same process of the previous subsection by fixing n 0 to a typical experimental value and using the resulting value of µ 0 to compute the number of particles N . The results of this section are partially based on Ref. 38. The time-independent GP equation in an ideal infinite optical lattice whose potential has the same form of the long-time potential of Eq. (4), V (x) = V ∞ cos 2 (k L x), reads, after rescaling the wave function and the coordi- where n r is the average atomic density, µ is the chemical potential, and E L = 4 2 k 2 L /m = 8E R .
We look for solutions of the Bloch form Ψ 0 (z) = e iqz y q (z), (B2) with y q (z + 2π) = y q (z) periodic, because the non-linear term is periodic for a Bloch-wave type solution. The normalization condition reads The Brillouin zone is placed in the region −1/2 < q < 1/2. The equation for y q is: y q +v cos(z)y q +c 2 |y q | 2 y q = α q y q , (B4) where we have allowed for a q-dependence of α defined in ( where M is a numerically enforced cut-off. After substitution of this solution in (B4) and in (B3), we get 2M + 2 equations for 2M + 2 variables (the 2M + 1 values of the Fourier coefficients c n plus the eigenvalue α). Instead of directly solving these non-linear equations, it is more efficient to minimize the quadratic sum of the 2M + 2 equations, where f j (c n , α) = 0 (with j = 1, 2 . . . 2M + 2) are the equations to be solved. It is easy to see that all these equations are real, so the coefficients c n can be chosen as real numbers and there is no need to use complex conjugates in (B6). A given Bloch solution can be unstable, either dynamically or in the sense of Landau, as explained below. The GP wave function is an extreme of the grand canonical Hamiltonian A superflow through the lattice is obtained when the actual solution Ψ 0 (z) minimizes the functional K [Ψ(z), α]. When this is not the case, the system can minimize its energy by the emission of excitations (phonons). The mean field solutions of this last type are said to exhibit Landau instabilities. These instabilities can be sought by expansion of K [Ψ(z), α] around Ψ 0 (z), i.e., Ψ(z) = Ψ 0 (z) + δΨ(z). The first-order term is automatically zero because Ψ 0 (z) solves the GP equation. The quadratic terms reads Landau instabilities correspond to negative eigenvalues of the Hermitian operator Λ. The corresponding eigenvalue equation is By absorbing the exponential plane-wave factor of the Bloch-type GP solution, u(z) = e iqz u q (z) and v(z) = e −iqz v q (z), we arrive at a new matrix operator Λ q which is periodic. Applying Bloch's theorem in the form u q (z) = e ikz u q,k (z) and v q (z) = e ikz v q,k (z) with u q,k (z) and v q,k (z) periodic in [0, 2π], the final eigenvalue equation reads Dynamical instabilities correspond to modes that grow exponentially with time. They are computed by looking for non-real eigenvalues of the BdG equations, which are formally similar to Eq. (B9): with M = σ z Λ (here σ z = diag(1, −1) is the usual Pauli matrix). Bloch's theorem also applies here and after a computation similar to that which has led to Eq. (B10), the eigenvalue equation reads with M q,k = σ z Λ q,k . In addition to the dynamical stability analysis, the real eigenvalues of this operator can be used to compute the speed of sound in the optical lattice. When q = 0, it can be proven that the the small wavevector k eigenvalues goes like ǫ = ±s|k|, with s the speed of sound (here, in units of 2 k L /m). On the other hand, the form of the Bloch-type solution of the GP equation as a function of both q and n r can be directly used to compute the sound speed without the need to solve for the BdG equations [28]. Restoring dimensions by introducing Q = 2k L q, it can be proven that where ∂ n denotes derivative with respect the mean density, n r , ∂ Q is the derivative with respect the pseudomomentum Q, with both derivatives evaluated at Q = 0, and E is an average energy density given by We will make use of this expression in Appendix C.
Appendix C: Perturbative treatment of the nonlinearity in the optical lattice In this Appendix, some of the results of Appendix B are perturbatively explored further. In Eq. (B1) there are two dimensionless parameters, v and c 2 . The former is essentially the amplitude of the potential (v ≫ 1 is the tight-binding regime while v ≪ 1 corresponds the nearly-free-particle regime) and the latter is a measure of the strength of the interaction or nonlinearity. In the cases studied in this paper, v ≪ 1. However, the forthcoming discussion applies to arbitrary values of v. The quantity δ ≡ c 2 ≪ 1 is the small parameter of our perturbation theory. We expand in powers of δ both the Bloch wave y q (z) and the displaced and dimensionless chemical potential α q , which solves Eqs. (B3)-(B4) and is defined in Eq. (B1). We obtain which transforms Eq. (B4) into [hereafter we omit the explicit z-dependence in y q (z), y H (0) q y q + δ|y q | 2 y q = α q y q In what follows, we focus on the lowest Bloch band and keep terms up to O(δ 2 ). The lowest-order term solves the linear Schrödinger equation, q . Therefore, y (0) q = φ q,0 and α (0) q = ε q,0 , where φ q,0 (z) is the corresponding eigenfunction for the lowest band, which involves Mathieu functions, and ε q,0 its eigenvalue (note the use of the index 0 for two different purposes: perturbative order is indicated in the superindex, while the band index comes in the subindex). We normalize φ q,0 according to (B3).
The first-order corrections must satisfy which, using standard perturbation techniques, leads to where {φ q,n } ∞ n=1 are the Schrödinger eigenvectors of the rest of bands and ε q,n its corresponding eigenvalues, for a given value of q.
The second-order equation reads We note that y is not needed to compute α (2) q . Specifically, we find which is always negative. Instead of invoking Mathieu functions, the numerical computation of the formulae presented in this Appendix [Eqs. (C4), (C6)] can be easily performed in a Fourier representation. As y q , φ q,n are periodic functions in [0, 2π], their Fourier expansion read y q (z) = ∞ m=−∞ c m e imz and φ q,n (z) = ∞ m=−∞ a n,m e imz . Both the c m and the a n,m coefficients can be chosen real (see previous results, Eqs. (C4), (C6) can be written as: β n a n β n = a ⊺ n ra 0 ε q,0 − ε q,n α (2) q = 3a ⊺ 0 rc (1) = −3 ∞ n=1 |β n | 2 (ε q,n − ε q,0 ). (C8) Now we can use the perturbative results (C8) to give approximate closed expressions for some parameters of the optical lattice. To first order in δ, the energy density is given by Eq. (B14) We can write where E L is defined after Eq. (B1). On the other hand, for the non-interacting chemical potential we have (assuming q ≪ 1/2): where E min is the bottom of the conduction band as defined in the main text, m * is the effective mass, and we recall Q = 2k L q. Thus we can rewrite Eq. (C10) as: Using (C12) we can rewrite (C9) as E ≃ n r E min + n r 2 Q 2 2m * + gn 2 Computing the derivatives to lowest order in δ, we arrive at: Similar results appear in Ref. 42 and references therein. The first square root is the speed of sound in the absence of the optical lattice. The second factor in the right takes into account the presence of the optical lattice and is practically unity for v ≪ 1. Specifically, we can write: and then mα (1) 0 /m * = 1 + O(v 4 ), so s ≃ gn r /m, which is the usual expression for the speed of sound.
The current is also conserved for a stationary solution of the GP equation and is given, to lowest order in δ and Q, by: In the nearly-free atom approximation, where the relative oscillations of the density around the mean value are small, we can write j ≃ n rv werev is a locally averaged flow velocity (not to be confused with the dimensionless parameter v). Then, we havev ≃ Q/m * and we can rewrite Eq. (C12) in a more appealing form: (D7) Because we ignore the value of Ψ j k+ 1 2 in the non-linear term, we use a corrector-predictor method, which consists in performing an additional step in every time iteration. In the first iteration, we use Ψ j k instead of Ψ j k+ 1 2 in order to obtain a valueΨ j k+1 . Next, we perform a new iteration taking Ψ j k+ 1 2 = Ψ j k+1 + Ψ j k /2 to obtain the final value Ψ j k+1 . The main advantage of this integration scheme is that the obtention of Ψ j k+1 only requires the resolution of a tridiagonal system of equations, which is computationally very efficient (the number of operations grows like N ).
The hard-wall boundary conditions reads Ψ l k = 0, l = 0, N + 1, and this can be easily implemented by suppressing the first and the last columns of the M matrices in (D6) On the other hand, any boundary condition imposed at the final point of the grid (x = L t ) will induce reflections which are unwanted because our goal is to simulate a semi-infinite supersonic region. To minimize those spurious reflections, one can use complex absorbing potential (CAP) at the grid boundaries [43]. Instead of that, we make use of the alternative, so-called ABC (Absorbing Boundary Conditions) [44,45]. This method is based on the linearization of the dispersion relation in the boundary in order to achieve the relation corresponding to an outgoing plane wave. Both ABC and CAP are very useful because they not only prevent the artificial reflection of the waves, but also permit to reduce the size of the supersonic zone.
The point x = L t is placed in the supersonic zone, where there is no potential. In addition, we expect that the non-linear term in (D2) can be neglected. This means that the effective Hamiltonian H in this region is the usual free (Schrödinger) Hamiltonian and Eq. (D2) can be written as which implies the dispersion relation On the supersonic side and in the quasi-stationary regime, the wave function is well peaked in momentum space around a value k 0 ≃ √ 2E min , where E min is the energy of the bottom of the first conduction band. By linearizing the dispersion relation around k 0 and expressing this relation in terms of derivatives, one can get Using the discrete version of this equation in j = N , we obtain a new equation that we add as a new row to the matrices M, which now are of size (N + 1) × (N + 1) and they will be of the form Due to the finite size of the grid and to the nonzero width of the momentum distribution in the supersonic region, the absorption is not perfect. We have found however that, in practice, the small spurious reflections have no effect on the final results.