Analogue cosmological particle creation in an ultracold quantum fluid of light

The rapid expansion of the early universe resulted in the spontaneous production of cosmological particles from vacuum fluctuations, some of which are observable today in the cosmic microwave background anisotropy. The analogue of cosmological particle creation in a quantum fluid was proposed, but the quantum, spontaneous effect due to vacuum fluctuations has not yet been observed. Here we report the spontaneous creation of analogue cosmological particles in the laboratory, using a quenched 3-dimensional quantum fluid of light. We observe acoustic peaks in the density power spectrum, in close quantitative agreement with the quantum-field theoretical prediction. We find that the long-wavelength particles provide a window to early times. This work introduces the quantum fluid of light, as cold as an atomic Bose-Einstein condensate.

T he expansion of a universe stretches all length scales, including the wavelengths of the particle modes. Thus, the frequencies of the modes evolve with time 1 , which implies that the modes at early and late times are related by a Bogoliubov transformation [2][3][4][5][6] . This field theory approach avoids the microscopic details, and predicts the spontaneous production of cosmological particles, including the primordial density fluctuations which led to the acoustic peaks in the cosmic microwave background (CMB) spectrum 4,6,7 . It is particularly relevant since the acoustic peaks can be described by linear perturbation theory 8 .
The field theory approach inspired the subject of analogue cosmological particle creation, in which laboratory experiments mimic the dynamics of scalar fields in curved space times [9][10][11][12][13][14][15] . The experiments even allow for measurement over time, which is impossible in the real universe, for which there is only one time of observation. Since the model is independent of the microscopic description of the medium, various quantum fluids were proposed for the study of cosmological particle creation in analogue universes [9][10][11][12][13][14][15] . In a two-dimensional atomic Bose-Einstein condensate, a qualitative comparison with cosmological particle creation was reported 16 . In a 1-dimensional experiment not related to quantum fluids, a rapid switch in the trapping field of two ions led to phonon pair creation and formation of spatial entanglement 17 .
Analogue cosmological particle creation is a type of dynamical Casimir effect [18][19][20] , which was observed in a superconducting circuit 21 and an optical fiber 22 . The classical, stimulated version of the effect was reported in a Bose-Einstein condensate 23 , but the observation of the quantum effect in a quantum fluid has not been reported. Pairs can also be produced by a modulational instability 24 .
We simulate expanding and contracting universes in a 3-dimensional quantum fluid of light, as coherent as an atomic Bose-Einstein condensate, and we observe time-resolved analogue cosmological particle creation out of vacuum fluctuations. Our quantum fluid is a near-resonant laser pulse traversing a warm atomic vapor cell, as illustrated in Fig. 1a. Within the vapor cell, the repulsive interactions between photons are mediated by the atoms, due to Kerr nonlinearity induced by the atomic resonance 25 . The interactions are suddenly quenched to zero when the laser beam exits the vapor cell 26 . This configuration mimics an expanding universe, since a rapid reduction of the interactions causes a sudden red shift of the energy spectrum [9][10][11][12][13] . We also observe the reverse process at the cell entrance, in which the interaction suddenly appears, mimicking a contracting universe. We demonstrate that both processes produce pairs of analogue cosmological particles, which confirms the predictions of Ref. 26 .

Results
Theoretical techniques. Our approach relies on the analogy between light propagation in a Kerr nonlinear medium, and the temporal dynamic of an atomic Bose-Einstein condensate. The effective time is τ ¼ z=c, where z is the position in the direction of propagation, and c is the speed of light. This effective time is equivalent to true time for the sake of quantum mechanical quasiparticle creation 26,27 . With no approximation other than the usual paraxial and slowly-varying envelope approximations 28 , we extend the standard monochromatic limit 28 and find that our fluid is described by the 3-dimensional Gross-Pitaevskii equation where ψ is the slowly-varying envelope of the electric field, ψ 2 is the volume density of the photons, m is their effective mass, U is an external potential, and g ψ 2 is the mean-field interaction energy. The three spatial dimensions of ∇ correspond to the transverse coordinates ðx; yÞ and to z 0 ¼ γðv g t À cτÞ, which is a coordinate comoving with the laser pulse at the group velocity v g , and compressed by the factor γ ¼ ðÀv 2 g k 0 D 0 Þ À1=2 , where k 0 is the wavenumber of the laser, and D 0 is the group velocity dispersion (see 3-dimensional Gross-Pitaevskii equation in Supplementary Methods). In other words, a laser pulse viewed in the z 0 coordinate would appear stationary and compressed relative to its length in the z coordinate. We study the analogue cosmological particles using the static structure factor, in analogy with the CMB power spectrum. The static structure factor has been used to study density fluctuations in Bose-Einstein condensates 16,29 , and we apply this technique to a fluid of light. It is given by Sðk x ; k y ; k z 0 Þ ¼ hjδρðk x ; k y ; k z 0 Þj 2 i=M, where δρðk x ; k y ; k z 0 Þ is the spatial Fourier transform of the density fluctuation at time τ, and M is the total number of particles in the fluid. With this definition, a zero-temperature, non-interacting gas has SðkÞ = 1, reflecting the presence of spatial shot noise. The operatorb y k corresponds to the creation of a quasiparticle after the quench, in mode k ¼ ðk x ; k y ; k z 0 Þ oscillating at frequency ω k . In the presence of quasiparticle populations N hb y kbk i and correlations C hb kbÀk i, the static structure factor within the Bogoliubov approximation is given by (see SðkÞ including absorption in Supplementary Methods) The populations and correlations are given by where N 0 hâ y kâk i and C 0 â kâÀk are the populations and correlations before the quench, respectively,â y k corresponds to the creation of a quasiparticle before the quench, and the operators are related by the Bogoliubov transformationb k ¼ αâ k þ βâ y Àk . For our series of two quenches, Eqs. (3) and (4) are applied twice. Since each quench either starts or ends with no interactions, α and β are the same Bogoliubov coefficients which diagonalize the Hamiltonian of a weakly-interacting quantum fluid (see SðkÞ including absorption in Supplementary Methods). In the absence of quasiparticles before a given quench, the pair production is spontaneous, and Eqs. (3) and (4) become N ¼ β 2 and C ¼ αβ. On the other hand, a distribution of quasiparticles before the quench, thermal or otherwise, will stimulate additional pairs. Experimental design. To create the fluid of light, we use a laser pulse with a 4 mm Gaussian waist and a power of 100 mW, propagating in an 85 Rb vapor cell heated to 150°C. A pulse length of 100 ns is employed to avoid saturating the camera used for observation. The laser is detuned −1.5 GHz (90 natural linewidths, 6 times the Doppler broadening) from the D2 5S 1=2 ; F ¼ 3 ! 5P 3=2 transition, giving v g = 0.007c. The interaction energy and healing length ξ = 60 µm are determined by the nonlinear change in the refractive index 4n, which is computed from the experimental parameters (see Determination of 4n in Methods). By taking into account the compression factor γ, this configuration leads to a weakly interacting photon gas with a thickness of 2 mm in the z 0 coordinate, and a dimensionless interaction coefficient ρa 3 s = 7 × 10 −14 , where ρ is the average photon density, and a s is the effective scattering length.
The fluid of light is imaged on a sCMOS camera, as shown in Fig. 1c. We tune the imaging system to pick out a certain z after the cell (fixing the effective time τ after the second quench), and the camera integrates over true time (thus integrating over z'), as illustrated in Fig. 1b. According to the Fourier slice theorem, this integration in position space gives a slice in k-space 30 . Thus, an ensemble of 200 images is obtained for each τ, and the power spectrum Sðk x ; k y ; k z0 ¼ 0Þ is computed by 2-dimensional Fourier transforms within the dashed square shown in Fig. 1c. The computation partially removes the effects of any drifts such as thermal convection, and accounts for the measured quantum efficiency of the camera (see Computation of Sðk x ; k y ; k z0 ¼ 0Þ in Methods).
Observation of analogue cosmological particle creation. In Fig. 2a we observe ring patterns in Sðk x ; k y ; k z0 ¼ 0Þ, oscillating as a function of k. These oscillations are the experimental signature of analogue cosmological particle creation, in close analogy with the acoustic peaks in the angular spectrum of the CMB. Pairs of quasiparticles with momenta ± k are generated at the moment of the quench with a random overall phase, but a definite phase relationship between þk and Àk, and oscillate with various frequencies ω k . Only certain k-values interfere constructively at the observation time τ, resulting in a ring pattern. The rings shrink with τ since lower frequencies take longer to develop oscillations. The shrinking pattern of rings is described quantitatively by Eq. (2). The radius of the first minimum in Fig. 2a is seen to be in good agreement with the theoretical prediction of Eq. (2), indicated by the dashed green curve. The azimuthal averages SðkÞ of Sðk x ; k y ; k z0 ¼ 0Þ are indicated in black in Fig. 2b. The red curves are calculated from Eq. (2), taking into account the two quenches, and the variations in α, β, and ω k which result from the measured absorption (see SðkÞ including absorption in Supplementary Methods). Very good agreement between the experimental black and theoretical red curves is seen. We also determine the spatial density correlations produced by the analogue cosmological particle creation. We derive the density-density correlation function g ð2Þ ðrÞ from S k ð Þ by the 3-dimensional spherically-symmetric Fourier transform  (2), and quantitative agreement with the experimental curves is seen.
Spontaneous particle creation in the first quench. The low-k behavior of S k ð Þ provides a window into the early times before the quenches, since the frequency of these modes approaches zero, so the modes do not have sufficient time to evolve during the experiment. The first peak in SðkÞ corresponds to the frequency 1=4τ, the lowest frequency which has time to oscillate. Well below this k-value, Eq. (2) reduces to S k ð Þ ¼ 1 þ 2N 1 , where N 1 is the incoherent population before the first quench, and the unity term corresponds to the quantum shot noise, which is scale invariant (independent of k). Thus, the value of S k ð Þ gives a direct measure of N 1 . Figure 3a shows the SðkÞ curves for all τ plotted together. We observe that SðkÞ is at most 1.4 for low k, as indicated by the dashed green line, giving N 1 ≤ 0.2. This value is finite and approximately scale invariant, which implies a negligible thermal component, since a thermal population diverges like 1=k. Furthermore, it is less than unity, implying that the spontaneous contribution dominates. Thus, the analogue cosmological particle creation is spontaneous in the first quench. This is verified by the blue and green curves in Fig. 3c, which show that stimulation in the first quench by thermal noise and white noise, respectively, would produce larger values of S k ð Þ than those of the experiment, for low k.
Stimulated particle creation in the second quench. The quasiparticles spontaneously created during the first quench stimulate pair creation in the second quench. However, if the particle production in the second quench were stimulated by the firstquench quasiparticles only, S k ð Þ would oscillate about unity, as indicated by the magenta curve in Fig. 3d. Rather, SðkÞ features an upward shift relative to unity, and a downward slope for large k. The downward slope is due to the finite resolution of the imaging system, measured to be 10 µm (see Measuring the imaging resolution in Supplementary Methods), and is included in the theoretical curves. The upward shift results from absorption and spontaneous reemission of photons from the atomic medium. By the first two terms of Eq. (2), SðkÞ oscillates about the value 1 þ 2 N 1 þ N b À Á , where N b is the background population present in the fluid between the two quenches. The population spontaneously created in the first quench does not contribute to this expression, since its spectrum (given by β 2 in Eq. (3)) does not extend to large k. The upward shift in Fig. 3a suggests N b = 1.2, a value which agrees well with our estimate for spontaneous reemission (see Absorption and spontaneous reemission in Supplementary Methods). The theoretical curves in Fig. 2b include this additional stimulation. While this incoherent, flat spectrum of 1.2 quasiparticles per mode implies that the fluid is not in its ground state, like a finite-temperature Bose-Einstein condensate, it does not negate the oscillatory behavior of S k ð Þ, and it even enhances the visibility of the oscillations. We can control this population by tuning the atomic density, the pulse duration, intensity, and detuning. In Fig. 3e we verify that this population vanishes for long weak pulses as expected, due to the finite coherence time of the spontaneous reemission. The unity value of S k ð Þ confirms that the fluid is shot-noise limited, when the effect of the atomic medium is absent.
Although our fluid of light is not in thermal equilibrium between the two quenches, we can put an upper limit on the effective temperature of the thermal component before the second quench. The blue curve in Fig. 3d includes thermal stimulation with an effective temperature 2mc 2 s = 30 mK, where c s is the speed of sound for the Bogoliubov quasiparticles, which results in a greatly enhanced first peak. Since this enhanced peak is absent from the experimental curve, we estimate the effective temperature of the thermal component to be less than 2mc 2 s , as in an atomic Bose-Einstein condensate. For the second quench, the thermal component does not diverge like 1=k since the zerotemperature static structure factor in the fluid of light goes to zero for low k (Ref. 31 ).
Interference pattern and the dispersion relation. Figure 3a exhibits a beating pattern in the envelope of the various curves, resulting from interference between analogue cosmological particles created in the two quenches. The theoretical curves in Fig. 3b show a similar pattern. The envelope has nodes and antinodes at ω 12k p ¼ πp=2τ 12 , where p is an integer (see Beating pattern in Supplementary Methods). By identifying each k p as shown in Figs. 3a, 4 points on the dispersion relation are found, as indicated by blue points in Fig. 4a. These points agree well with the dispersion relation in the medium, calculated from the interactions, and indicated by the blue curve.
Individual modes. Figure 4b shows the curves of Fig. 3a, one above the other. By plotting the SðkÞ values along the dashed line, we obtain the time dependence of a given mode k, as shown in Fig. 4c. Each mode is seen to oscillate sinusoidally after the second quench, with no visible damping. The frequencies of the oscillations, indicated by the black curve in Fig. 4a, agree well with the free-particle spectrum, which is relevant after the second quench.
Comparison with the CMB power spectrum. We compare and contrast our observed spectra with the CMB power spectrum in Fig. 5. Since 1990, several successive space missions have improved the resolution of the CMB measurements [32][33][34][35] , and Fig. 5b shows the latest results 35 . In the CMB spectrum, the oscillations occur as a result of the well-defined phase between the cosmological particles 36,37 , which is also true for our spectra. In the early universe, the density fluctuations relevant for the acoustic peaks oscillated until the effective time of observation τ, when the photons decoupled from matter 38 . The mode number l shown in Fig. 5b is proportional to k, when mapped back to the density fluctuations in the early universe, and the first peak likely corresponds to ω k τ % π (Refs. [39][40][41]. In contrast, the Bogoliubov transformation predicts a first peak in the CMB spectrum at ω k τ % π=2 (Ref. 5 ). Figure 5a shows our spectra for τ>0. For visual comparison with the CMB spectrum, k is divided by k 0 , which satisfies ω k 0 τ ¼ π=2, where ω k is the magenta curve of Fig. 4a. There are features which are common to our spectra and that of the CMB, in addition to the oscillations: the decay of the oscillations for large k or l, and the approximately scale-invariant region for small k or l. The oscillations in the CMB spectrum decay for large l due to damping by photon diffusion 39 . In contrast, the oscillations in our spectra decay because the β Bogoliubov coefficient in Eq. (4) decreases for high k. The scale-invariant region of the CMB spectrum arises from quantum fluctuations [42][43][44] , assuming that the inflation model is correct [45][46][47][48] . Similarly, the scale invariant part of our spectra reflects the quantum nature of the particle production, as a result of the vacuum of incoming particles. However, our incoming vacuum is a property of our shot-noise limited laser, as opposed to red-shifting of the modes which possibly occurred during inflation 42,49,50 . Red-shifting was observed in the laboratory 51 , and it would be interesting to combine it with analogue cosmological particle production.

Discussion
This work establishes the paraxial fluid of light as a quantum fluid. The results demonstrate that quantum field theory applies to a system in which a spatial coordinate plays the role of time. The effective temperature is less than twice the interaction energy, which is comparable to many atomic Bose-Einstein condensates.
On the other hand, the apparatus is an order of magnitude simpler, smaller, and less expensive. The direct detection of the photon fluid is also an advantage.
In conclusion, we observe both spontaneous and stimulated analogue cosmological particle creation in a quantum fluid of light. The particle production in the first quench is seen to be spontaneous, while the second includes stimulation by the first quench quasiparticles, as well as by an incoherent background. We quantitatively confirm the quantum field-theoretical prediction. The long wavelength part of the spectrum provides a window into early times before the particle creation. From an alternative perspective, we observe the spontaneous and stimulated dynamical Casimir effects in a quantum fluid.

Methods
Determination of 4n. The interaction between photons is quantified by the nonlinear contribution to the index of refraction n, given by 4n ¼ n I ð Þ À nð0Þ. We would like to express 4n in terms of easily measurable quantities. The index of refraction is given by n ¼ , where χ is the atomic susceptibility. Since n % 1, one obtains 4n ¼ 4 Re χ À Á Â Ã =2. Furthermore, the absorption coefficient is given by α a ¼ k 0 Im χ À Á . Thus, Also, χ is proportional to i À 2δ=Γ , where δ is the detuning from resonance, Γ is the linewidth, and I sat is the saturation intensity. This gives The detuning is δ = −1.5 GHz relative to the 85 Rb cooling transition, and the self-broadened linewidth is Γ=2π = 16 MHz (Ref. 52 , for the vapor cell temperature of 150°C, corresponding to an atomic density of 1 × 10 20 m −3 . The intensity is given by I ¼ 2P=πw 2 , where the waist of the beam is w = 4 mm, and the laser power P decays exponentially due to absorption, from 100 mW at the entrance to the vapor cell, to 20 mW at the exit. We estimate I sat to be the far-detuned, π-polarized value, 25 Wm −2 (Ref. 53 . The absorption coefficient is given by α a ¼ À lnT ð Þ=L, where T = 0.2 is the transmission through the vapor cell of length L = 10 mm. The wavenumber is given by k 0 ¼ 2π=λ, where λ = 780 nm. Equation (7) yields 4n = −8.6 × 10 −6 and −1.7 × 10 −6 for the entrance and exits of the vapor cell, respectively. We have neglected the effect of optical pumping into the dark ground state. Via measurements of 4n and α a , we find the optical pumping time to be a few microseconds, so optical pumping is negligible during the 100 ns pulse employed in this work.
Computation of Sðk x ; k y ; k z0 ¼ 0Þ. The power spectrum (static structure factor) Sðk x ; k y ; k z0 Þ of a system of M particles (photons in our case) is given by S k x ; k y ; k z0 ¼ δρ k x ; k y ; k z0 where δρðk x ; k y ; k z 0 Þ ¼ R dx dy dz 0 δρðx; y; z 0 Þ e Àiðk x xþk y yþk z 0 z 0 Þ and the density fluctuation δρðx; y; z 0 Þ ¼ ρðx; y; z 0 Þ À hρðx; y; z 0 Þi. Setting k z0 ¼ 0, one obtains the 2-dimensional Fourier transform δρ k x ; k y ; k z0 ¼ 0 ¼ Z dx dy δe ρ x; y À Á e Ài k x xþk y y ð Þ ð9Þ where e ρ x; y À Á ¼ R dz 0 ρ x; y; z 0 À Á is the number density integrated in the z0 direction, a quantity we measure directly on the camera. The density fluctuation δe ρ ¼ e ρ À e ρ 5 is computed for each image, where e ρ 5 is the average of 5 adjacent images rather than the average e ρ over the entire ensemble. This technique reduces the effects of drifts in the experimental parameters during the 7 s required to obtain the ensemble. As mentioned in relation to Fig. 2c, the relative density fluctuation δe ρ= e ρ is on the order of 10 −3 , so small drifts can play a role. For example, thermal convection of the 85 Rb gas may induce small changes in the shape of the fluid of light from image to image. The 2-dimensional Fourier transform in Eq. (9) is computed for each image within the dashed square of Fig. 1c. The power spectrum Sðk x ; k y ; k z0 ¼ 0Þ is computed by Eq. (8). The use of e ρ 5 rather than e ρ reduces the fluctuations by a factor of 4/5. Thus, Sðk x ; k y ; k z0 ¼ 0Þ is multiplied by 5/4 to correct this effect. Furthermore, the finite quantum efficiency Q = 0.485 of the camera tends to randomize the photon density and bring Sðk x ; k y ; k z0 ¼ 0Þ closer to unity. Thus, Sðk x ; k y ; k z0 ¼ 0Þ À 1 is multiplied by the factor 1=Q.

Data availability
Source data are provided with this paper. The images generated in this study have been deposited in Zenodo (https://doi.org/10.5281/zenodo.6438403). Source data are provided with this paper.