Analog reheating of the early universe in the laboratory
Aleksandr Chatrchyan, Kevin T. Geier, Markus K. Oberthaler, Jürgen Berges, Philipp Hauke
AAnalog reheating of the early universe in the laboratory
Aleksandr Chatrchyan, ∗ Kevin Geier,
1, 2, 3, ∗ Markus K. Oberthaler, Jürgen Berges, and Philipp Hauke
2, 1, 3 Institute for Theoretical Physics, Ruprecht-Karls-Universität Heidelberg, Philosophenweg 16, 69120 Heidelberg, Germany INO-CNR BEC Center and Department of Physics, University of Trento, Via Sommarive 14, 38123 Povo (TN), Italy Kirchhoff Institute for Physics, Ruprecht-Karls-Universität Heidelberg, Im Neuenheimer Feld 227, 69120 Heidelberg, Germany (Dated: August 7, 2020)Cosmic reheating describes the transition of the post-inflationary universe to a hot and thermalstate. In order to shed light on the nature of this process, we propose a quantum simulation ofcosmic reheating in an ultracold Bose gas. In our model, the inflaton field dynamics is mappedonto that of an atomic Bose–Einstein condensate whose excitations are identified with the particlesproduced by the decaying inflaton field. The expansion of the universe as well as the oscillationsof the inflaton field are encoded in the time-dependence of the atomic interactions, which can betuned experimentally via Feshbach resonances. As we illustrate by means of classical-statisticalsimulations for the case of two spatial dimensions, the dynamics of the atomic system exhibits thecharacteristic stages of far-from-equilibrium reheating, including the amplification of fluctuationsvia parametric instabilities and the subsequent turbulent transport of energy towards higher mo-menta. The transport is governed by a non-thermal fixed point showing universal self-similar timeevolution as well as a transient regime of prescaling with time-dependent scaling exponents. Whilethe classical-statistical simulations can only capture the earlier stages of the dynamics for weak cou-plings, the proposed experimental implementation provides a protocol for the quantum simulationof the entire evolution even beyond the weak coupling regime.
I. INTRODUCTION
Cosmic inflation is a well-established paradigm thatsolves the puzzles of flatness, homogeneity, and isotropyof the universe [1, 2], and explains the generation of den-sity perturbations for structure formation [3]. Accordingto this paradigm, the early universe underwent a stageof accelerated expansion which, in a typical scenario, isdriven by a scalar field, the inflaton. Inflation is followedby a reheating phase [4], during which the inflaton de-cays, e.g., into degrees of freedom of the standard modelof particle physics. The heating process can start witha preheating stage of rapid particle production far fromequilibrium, and finally completes after the particles, dueto their interactions, approach thermal equilibrium at thereheating temperature [5, 6]. Over the recent years, it hasbeen realized that many salient features of the dynam-ics become universal in the far-from-equilibrium regime,across physical systems at vastly different scales [7–12].This universality makes it particularly appealing to gaininsight into the non-equilibrium processes during reheat-ing by use of table-top experiments, such as ultracoldatomic gases. Thanks to their high degree of control-lability, these systems are ideally suited to realize andstudy processes relevant for the early universe that areotherwise challenging to access [13–21].Here, we propose a quantum simulation of reheatingdynamics in a single-component Bose gas, leveraging theability to modulate the atomic interaction in time bymeans of Feshbach resonances [22]. To this end, we con- ∗ These two authors contributed equally.E-mail A.C. at: [email protected] K.G. at: [email protected]. sider a generic model where the inflaton is described by aminimally-coupled, relativistic scalar field in an expand-ing Friedmann–Robertson–Walker (FRW) universe [23].By taking the non-relativistic limit, we establish a closerelation between the post-inflationary universe and an ul-tracold Bose gas on the quantum level. At the heart ofour mapping lies a time-dependent rescaling of both theBose field operator and the spacetime coordinates, allow-ing us to map the dynamics of an expanding Bose gas intothat of a non-expanding one with a time-dependent inter-
Potential S i m u l a t i n g s y s t e m S i m u l a t e d s y s t e m In fl aton fi eld Parametricinstabilities
Scattering length Time
BEC
Turbulentthermalization T h e r m a l e q u ili b r i u m (P)reheating of the universe S e l f - s i m i l a r e v o l u t i o n I n fl a t i o n Figure 1. Schematic illustration of reheating dynamics in theearly universe and its simulation in the laboratory. Aftercosmic inflation, we consider the inflaton oscillating aroundthe minimum of its potential, producing particles via para-metric instabilities (“preheating”). This process drives thesystem into a turbulent state where energy is transported to-wards higher momenta in a self-similar way as the universeapproaches thermal equilibrium (“reheating”). In the simula-tion, the oscillating inflaton is mimicked by modulating thescattering length of a Bose–Einstein condensate (BEC). a r X i v : . [ c ond - m a t . qu a n t - g a s ] A ug action. This simplifies the simulation compared to pre-vious proposals by lifting the restrictions associated withthe expansion of the trap [24]. Remarkably, for the par-ticular case of two spatial dimensions, where the systemexhibits scale invariance [25–27], the dynamics becomesindependent of the expansion, such that arbitrary ex-panding backgrounds may be studied by post-processingthe data of a single experiment. As we demonstrate bymeans of classical-statistical (or truncated Wigner) sim-ulations [28–31], the proposed setup exhibits the char-acteristic stages of reheating dynamics [6], including thepreheating stage of parametric amplification of quantumfluctuations as well as the self-similar transport of energytowards higher momenta driving turbulent thermaliza-tion [7]. A schematic illustration of the reheating processand its proposed simulation in the laboratory is shownin fig. 1.Analog models of gravity build on emergent curvedspacetimes for low-energy excitations on top ofBose–Einstein condensates (BECs) or other systems [32].These can be generated, for instance, by modulat-ing trapping potentials or interactions, allowing one tosimulate cosmological particle production during infla-tion [13, 14, 33]. However, non-linear effects crucial forpost-inflationary reheating dynamics cannot easily be in-corporated into this framework. Nonetheless, severalexperimental setups have been proposed for simulatingpost-inflationary dynamics in Bose gases [16–20]. Besidesthese theoretical proposals, preheating-like dynamics hasrecently been observed in experiments with a ring-shapedcondensate after rapidly expanding the trap [24]. In oursetup, parametric resonance, arising from the inflatonfield oscillating around the minimum of its potential [5],is mimicked in a Bose gas by modulating the atomic in-teraction. This modulation induces a preheating-like res-onant amplification of sound waves on top of the BEC.Parametric instabilities in Bose gases have been studiedextensively in the literature [34–39] and the concept hasevolved into a promising tool for state preparation [40].Going beyond the well-understood linear regime of para-metric resonance, our main focus is geared towards thenon-linear stages of the dynamics, involving a secondaryamplification of excitations [5, 19, 41] that herald thedevelopment of turbulence [7].The preheating-like instabilities drive the system faraway from equilibrium into the vicinity of a non-thermalfixed point, where the dynamics becomes universal andself-similar [42]. A non-thermal fixed point character-izes the far-from-equilibrium behavior of an entire uni-versality class of physical systems as diverse as the in-flaton in the early universe, ultracold atomic gases inthe laboratory or even the quark-gluon plasma exploredin heavy-ion collisions [7, 8, 10, 43]. Such universal dy-namics far from equilibrium has been observed exper-imentally in spinor and tunnel-coupled Bose gases inform of inverse cascades transporting conserved quanti-ties towards lower momenta [11, 12]. Bidirectional cas-cades involving additionally a self-similar transport of energy towards higher momenta have been studied nu-merically in one-dimensional ( D) spinor Bose gases [44]and have recently been observed experimentally in an iso-lated three-dimensional ( D) Bose gas following a cool-ing quench [45]. Our protocol of parametrically excitinga pure single-component BEC constitutes a complemen-tary path of entering the regime of self-similar time evo-lution far from equilibrium and allows one to study thedirect energy cascade in both driven and isolated sys-tems.We distinguish a regime of driven turbulence in thepresence of continuous parametric driving and a regimeof free turbulence after the modulation of the interactionis switched off, corresponding to the decay of the inflaton.Recently, a new type of prescaling phenomenon [46, 47]has been established, which is closely related to the on-set of a hydrodynamic behavior far-from-equilibrium [46].In this regime, the dynamics of the system is governedby the fixed-point scaling function and time-dependentscaling exponents, which slowly relax to their univer-sal values. Remarkably, our numerical results indicatethat such a behavior emerges during the transition fromdriven to free turbulence, which opens the door for an ex-perimental observation of this phenomenon in cold Bosegases.At late times, when the occupancies at characteristicmomenta become of order unity, the system is expectedto thermalize at the reheating temperature [7]. This fi-nal stage of the dynamics is dominated by quantum fluc-tuations and therefore cannot be captured by means ofclassical-statistical simulations [29, 48]. In contrast, itcan be accessed by experimental studies in cold atomicsystems, which are quantum-mechanical by nature.This paper is organized as follows. In section II, wepresent our mapping of the inflaton to an ultracold Bosegas and illustrate that an expanding spacetime can beencoded in the time-dependence of the interaction. Fo-cusing on the case of a two-dimensional ( D) Bose gas,we proceed in section III to discussing an analog im-plementation of preheating, where particle production ismimicked by parametrically exciting Bogoliubov quasi-particles via a modulation of the scattering length. Weexplain and interpret the momentum spectrum result-ing from these instabilities, emphasizing the importantrole of non-linear effects leading to the onset of turbu-lent dynamics. In section IV, we analyze the universalself-similar time evolution for both the driven and thefree direct cascade, comparing the numerically extractedscaling exponents with predictions from kinetic theory.We also show that the system exhibits prescaling duringthe transition from driven to free turbulence. Moreover,we illustrate how expansion affects the dynamics and can,under certain conditions, even prevent the system fromthermalizing. In section section V, we discuss detailsrelevant for an experimental implementation, before weconclude in section VI.
II. EXPANDING SPACETIME IN BOSE GASES
The inflaton is commonly described by a relativisticreal scalar field φ ( x ) in curved spacetime with the ac-tion [49] S = (cid:90) d d +1 x √− g (cid:20) g µν ∂ µ φ∂ ν φ − V ( φ ) (cid:21) , (1)where V is the potential, g µν is the metric tensor, g itsdeterminant, and d the number of spatial dimensions.The expansion of the universe is well described with thehelp of the flat FRW metric [49]. It corresponds to ahomogeneous and isotropic universe, in accordance withthe cosmological principle, which is Minkowskian at eachtime slice. The FRW metric has the form d s = g µν d x µ d x ν = c d τ − a ( τ ) d x , (2)where c is the speed of light and τ denotes the cosmictime. The cosmic scale factor a ( τ ) grows with time in anexpanding universe and relates the comoving distance x to the proper distance r ( τ ) = a ( τ ) x .Our goal is to simulate the dynamics of the inflaton inan expanding universe using an ultracold Bose gas, whichis described by a non-relativistic complex field operator ˆΨ subject to the Hamiltonian [50] ˆ H = (cid:90) d d x (cid:20) ˆΨ † (cid:18) − (cid:126) ∇ m + U (cid:19) ˆΨ + g † ˆΨ † ˆΨ ˆΨ (cid:21) , (3)where (cid:126) is the reduced Planck constant, m the atomicmass, U an external trapping potential, and g the quar-tic coupling which determines the strength of the atomicinteractions via the s -wave scattering length.The most straightforward approach to incorporate ex-pansion in a trapped Bose gas is to physically expand thetrap geometry [15, 24, 26]. While this approach is sim-ple and direct in principle, there are practical limitationssuch as the restriction to short times or small expansionvelocities due to a finite optical imaging system.Alternatively, one can keep the trap geometry fixed,associating physical distances with the co-moving dis-tances of an expanding system. In this way, however,the redshift of momenta as well as dilution of the sys-tem due to the expansion have to be explicitly accountedfor in form of a modified kinetic term and a non-unitaryevolution, as discussed in the subsequent subsection (seealso ref. [51]). Besides possible technical difficulties inengineering such non-Hermitian Hamiltonians, this for-mulation has the additional drawback that it implies adecreasing atomic density, which diminishes the experi-mental signal with progressing lab time.We follow here a third approach, wherein an appropri-ate time-dependent rescaling of both the field operatorand the time coordinate allows us to encode the expan-sion solely in the time-dependence of the atomic interac-tion. The transformation back to the original variables isperformed in a post-processing step after the simulation. For experiments with ultracold atoms where a broad Fes-hbach resonance is available to tune the interactions, thisthird approach may be the preferred and most elegantway to incorporate a broad range of expansion scenarios.The possibility of modulating the interaction in orderto realize an effective expanding spacetime for linear ex-citations on top of the condensate has been considered,e.g., in ref. [33]. Here, we take a related point of view,which, however, differs in an essential aspect: since weare interested in simulating reheating dynamics, whichis a non-linear process, we consider the expansion of thenon-relativistic system as whole instead of generating anemergent expanding spacetime for the linear excitations.In what follows, we will introduce and discuss the corre-sponding mapping in detail. In sections III and IV, wewill present the physical phenomena of preheating andreheating that become accessible thanks to this mapping. A. Inflaton dynamics in expanding spacetime
We model the inflaton as a real scalar field, whose dy-namics is described by the action (1) with a simple poten-tial V ( φ ) = m φ / λφ / , where m is the mass and λ a quartic coupling.The classical equations of motion for the inflaton fieldcan be obtained from eq. (1) using the principle of leastaction. For the metric defined in eq. (2) they read c ¨ φ + 1 c dH ˙ φ − ∇ φa + m c (cid:126) φ + λ φ = 0 , (4)where H = ˙ a/a is the Hubble parameter and the dot de-notes the derivative with respect to the cosmic time τ .This equation is of a Klein–Gordon type with a modifiedspatial derivative term and an additional friction term,which express the redshift of momenta in the co-movingframe and the dilution of the field due to the expansion,respectively [49].Canonical quantization can be performed as usuallyin QFT on curved spacetime by transforming to con-formal time, d η = d τ /a , and conformal field variables, Φ = a ( d − / φ , [52]. The equations of motion then takethe same form as in Minkowski spacetime, c Φ (cid:48)(cid:48) − ∇ Φ + m c (cid:126) Φ + λa − d = 0 , (5)with a time-dependent coupling and effective mass m ( a ) = m a − (cid:126) ( d − c (cid:104)(cid:16) a (cid:48) a (cid:17) ( d − a (cid:48)(cid:48) a (cid:105) . (6)Here, ( · ) (cid:48) denotes the derivate with respect to confor-mal time η . The field Φ and its conjugate momentum Π = Φ (cid:48) /c are then promoted to operators that satisfy thecanonical commutation relations (cid:2) ˆΦ( η, x ) , ˆΠ( η, y ) (cid:3) = i (cid:126) δ ( x − y ) , (cid:2) ˆΦ( η, x ) , ˆΦ( η, y ) (cid:3) = (cid:2) ˆΠ( η, x ) , ˆΠ( η, y ) (cid:3) = 0 . (7) B. Mapping the inflaton to an ultracold Bose gas
Having introduced the transformation to conformalvariables in the previous section in order to quantize thesystem, we now switch back to the original variables τ and ˆ φ = a − ( d − / ˆΦ , since the non-relativistic limit ismost conveniently studied in terms of these.The mapping to the non-relativistic description of aBose gas is achieved by introducing a complex field op-erator ˆ ψ ( τ, x ) such that ˆ φ = (cid:126) √ mc (cid:104) ˆ ψ e − i (cid:126) mc τ + h . c . (cid:105) , (8)where h . c . denotes the Hermitian conjugate.In the mapping (8), we have factored out the rapid os-cillations of the inflation field on the frequency scale de-termined by its rest mass. For taking the non-relativisticlimit, we require the field to change much slower thanthese oscillations, i.e., we demand the mass term tobe the dominant contribution on the right-hand side ofthe equations of motion (4). This requires the follow-ing conditions to hold. First, the typical field valuesshould be not too large, λ (cid:104) ˆ φ (cid:105) (cid:28) m c / (cid:126) , which impliesthat the self-interactions are relatively weak. Similarly,the typical physical momenta should be non-relativistic, | p | (cid:28) mc . Finally, the expansion should not be too rapid, H = ˙ a/a (cid:28) mc / (cid:126) , since otherwise the field would be inthe over-damped regime and not perform oscillations atall, as during inflation. Note that the first two assump-tions are not necessarily satisfied by the inflaton field,especially at the early stages of reheating. We will re-turn to this issue in section III.The outlined assumptions imply the hierarchy (cid:12)(cid:12)(cid:12) (cid:104) ˆ ψ · · ·(cid:105) (cid:12)(cid:12)(cid:12) (cid:29) (cid:12)(cid:12)(cid:12) (cid:126) mc (cid:104) ˙ˆ ψ · · ·(cid:105) (cid:12)(cid:12)(cid:12) (cid:29) (cid:12)(cid:12)(cid:12) (cid:126) m c (cid:104) ¨ˆ ψ · · ·(cid:105) (cid:12)(cid:12)(cid:12) , (9)where (cid:104)·(cid:105) denotes the expectation value. These inequal-ities are to be understood to hold for any correlationfunction involving the field or its derivatives.Inserting the ansatz (8) into eq. (4) without any ap-proximation yields (cid:34) c (cid:18) ¨ˆ ψ − imc (cid:126) ˙ˆ ψ + dH (cid:18) ˙ˆ ψ − imc (cid:126) ˆ ψ (cid:19)(cid:19) − ∇ a ˆ ψ (10) + λ (cid:126) mc (cid:16) ˆ ψ e − i (cid:126) mc τ + 3 ˆ ψ ˆ ψ † ˆ ψ (cid:17)(cid:35) e − i (cid:126) mc τ + h . c . = 0 . Employing the inequalities (9) allows us to neglect twoof the first four terms in eq. (10). Furthermore, we canuse relations (9) to neglect the rapidly oscillating termsproportional to ˆ ψ e ± imc τ/ (cid:126) in the last bracket. The lat-ter describe number-changing processes in the relativistictheory and neglecting them leads to an emergent U (1) symmetry implying particle number conservation. As aresult, eq. (10) reduces to [51] i (cid:126) ˙ˆ ψ = (cid:18) − (cid:126) m ∇ a − i (cid:126) d H + g ˆ ψ † ˆ ψ (cid:19) ˆ ψ, (11) where we defined the coupling g = λ (cid:126) / m c and normal-ordered the operators, dropping an irrelevant energyshift .Equation (11) has the familiar form of the Heisenbergequations of motion generated by the Hamiltonian of anultracold Bose gas (3) in absence of an external trappingpotential. The kinetic term is proportional to a − , whichdescribes redshift of momenta in the co-moving frame,and the dilution of the system due to the expansion isexpressed by a non-Hermitian term causing the norm ofthe field to decay.In analogy to the transformation to conformal vari-ables in the relativistic system, we now introduce a newtime variable t via the relation d t = d τ /a , which we re-fer to as the lab time, and the rescaled field operator ˆΨ = ˆ ψa d/ . By virtue of relations (9), it is straightfor-ward to show that within the outlined approximationsthe canonical commutation relations (7) imply bosonicequal-time commutation relations for the field ˆΨ , (cid:2) ˆΨ( t, x ) , ˆΨ † ( t, y ) (cid:3) = δ ( x − y ) , (cid:2) ˆΨ( t, x ) , ˆΨ( t, y ) (cid:3) = (cid:2) ˆΨ † ( t, x ) , ˆΨ † ( t, y ) (cid:3) = 0 . (12)Furthermore, the equations of motion (11) take the stan-dard form for a non-relativistic bosonic field, i (cid:126) ∂ ˆΨ ∂t = (cid:18) − (cid:126) ∇ m + g eff ( t ) ˆΨ † ˆΨ (cid:19) ˆΨ , (13)with the time-dependent effective coupling g eff ( t ) = ga − d ( t ) . (14)It is this form of the equations of motion that we will usein the remainder of this work. C. The special case of expansion in 2D
Remarkably, for d = 2 , the coupling in eq. (13) be-comes independent of the scale factor a ( t ) . This is aconsequence of a dynamical symmetry in D Bose gasesknown as scale invariance [25]. While quantum anoma-lies in strongly interacting systems can lead to violationsof scale invariance [53, 54], it has been well tested ex-perimentally in the weakly interacting regime [27, 55].Thus, if scale invariance holds, the equations of motionare the same as in the case of a static spacetime. Thenature of the expansion, encoded in the scale factor a ( t ) ,then only enters the back transformation to the origi-nal cosmic temporal and spatial coordinates as well asfield variables. This makes the simulation in two spatial Indeed, normal-ordering g ˆ ψ ˆ ψ † ˆ ψ produces the commutator g [ ˆ ψ, ˆ ψ † ] ˆ ψ = ga − d δ (0) ˆ ψ , as follows from eqs. (7) and (8). Thisterm can then be absorbed into an irrelevant dynamical phasevia the transformation ˆ ψ → ˆ ψ exp { i (cid:126) δ (0) (cid:82) ga − d dτ } . dimensions particularly efficient, since the evolution ofa single experiment can conveniently be mapped to ar-bitrary expanding spacetimes in a post-processing step.For this reason, we will focus on weakly interacting DBose gases in the remainder of this work. Although ouruniverse is clearly not two-dimensional, this geometrycaptures most of the essential physics of reheating dy-namics, as we demonstrate in the following sections.
III. ANALOG PREHEATING
Having established the connection between a scalarfield in an expanding spacetime and an ultracold Bosegas, we are now in the position to discuss how the pre-heating stage of post-inflationary dynamics can be simu-lated in a Bose gas.
A. Preheating in the early universe
In this subsection, we briefly review the basic mecha-nism of preheating in the early universe. After inflation,the energy budget of the universe is stored predominantlyin the homogeneous inflaton field, which oscillates aroundthe minimum of its potential. The stage of preheatingconsists of an explosive, non-perturbative particle cre-ation from the decaying inflaton.A common mechanism for preheating is via the para-metric amplification of quantum fluctuations in the pres-ence of the effective potential induced by the inflaton [5,41]. This effect can most conveniently be understood bylinearizing the fluctuations of the inflaton field aroundits homogeneous background, φ ( τ, x ) = φ ( τ ) + δφ ( τ, x ) ,and inserting it in eq. (4). The resulting equation de-scribes damped oscillations of the background field, c ¨ φ + 1 c dH ˙ φ + m c (cid:126) φ + λ φ = 0 . (15)The equations of motion for the fluctuations, dropping allterms of quadratic or higher order in δφ , read in Fourierspace δ ¨ φ p + dHδ ˙ φ p + ω p ( φ, τ ) δφ p = 0 (16)with ω p ( φ , τ ) /c = p / (cid:126) a ( τ ) + m c / (cid:126) + λφ / andthe Fourier transform δφ p ( τ ) = (cid:82) d d x δφ ( τ, x )e − i px / (cid:126) .Equation (16) describes a collection of parametric oscil-lators for each momentum mode p , driven by the back-ground oscillations of φ ( τ ) . Modes within certain insta-bility bands satisfy resonance conditions, which leads toan exponential growth of occupancies corresponding toparticle production [41, 56]. The growth of fluctuationseventually invalidates the linearized approach of eq. (16).In models where the inflaton is coupled to other bosonicfields, particle production via parametric resonances oc-curs in a similar way [5]. For simplicity, we will focushere on the case where the inflaton decays efficiently intoits own quanta of excitation. B. Analog preheating in Bose gases
Our aim is to observe an analog of preheating dynamicsin a cold-atom Bose gas. Motivated by the discussions inthe previous sections, we mimic the state of the universeafter inflation by a homogeneous BEC. (Possible impli-cations of the presence of an external trapping potentialare discussed in section V.) It is tempting to search fora direct analog of the inflaton’s oscillations, which couldtrigger parametric instabilities in the condensate. How-ever, this effect is not present for a single-component Bosegas in its ground state.In fact, the absence of parametric resonance is aconsequence of the non-relativistic limit considered insection II. To better understand this, we go back toeq. (10), which is an exact rewriting of the origi-nal relativistic equation of motion (4). Linearizingthe field operator around the homogeneous condensate, ˆ ψ ( τ, x ) = ψ ( τ ) + δ ˆ ψ ( τ, x ) , the interaction term takesthe form ∼ g (cid:16) ψ δ ˆ ψ e − i (cid:126) mc τ + 2 | ψ | δ ˆ ψ + ψ δ ˆ ψ † (cid:17) , (17)which corresponds to the term proportional to λφ δφ ineq. (16). When taking the non-relativistic limit to deriveeq. (13), we have neglected the rapidly oscillating U (1) -violating term proportional to e − imc τ/ (cid:126) , keeping onlythe U (1) -conserving, slowly-varying part. However, sincethe former reflects the oscillations of the inflaton respon-sible for parametric resonance, this effect is no longerpresent in the non-relativistic model.To re-introduce parametric instabilities in the non-relativistic model, we add a periodic modulation of thecoupling with frequency ω , such that the linearized in-teraction term in eq. (13) takes the form ∼ g (1 + sin( ωτ )) (cid:16) | ψ | δ ˆ ψ + ψ δ ˆ ψ † (cid:17) . (18)On a formal level, this is similar to eq. (17) in thesense that we have an oscillating prefactor multiplyinga term cubic in the field, although it does not have ex-actly the same form. In addition, the modulation fre-quency ω needs to correspond to typical non-relativisticenergies, which is another difference to the relativisticcase, where the inflaton oscillates at relativistic frequen-cies ω rel (cid:39) mc / (cid:126) on the scale given by its rest mass.Despite these formal differences, modulating the inter-action as described above suffices to generate the desiredparametric resonance phenomena in the Bose gas [34–40],analogous to the ones taking place in the preheating stageof most inflationary models [5, 41]. This is demonstratedin detail in the subsequent subsections. C. Numerical study of preheating dynamics
We illustrate our implementation of preheating dynam-ics by means of numerical simulations based on eq. (13),which corresponds to a static ultracold Bose gas witha time-dependent interaction encoding the expansion ofthe system. We consider here a Bose gas in D, wherethe particular model of expansion does not explicitly en-ter the equations of motion, as discussed in section II C.(Nonetheless, we keep the number of spatial dimensions d in all formulas general.) Equation (13) is then formallyequivalent to the equation of motion generated by theHamiltonian (3), and the back transformation to the orig-inal cosmic time and field variables can be performed inpost-processing. The latter is illustrated at the exampleof a matter-dominated expansion in section IV E.Following the argumentation in section III B, we induceparametric instabilities by harmonically modulating thecoupling as g ( t ) = g (1 + r sin( ωt )) , (19)where g is a positive offset value, r is the amplitude ofthe modulation and ω is its frequency. Note that herewe have used the lab time t in the argument of the sinefunction, whereas in eq. (18) the cosmic time τ appears,which, when expressed in terms of the lab time t , re-sults in oscillations with increasing frequency. This sub-stitution is based on the simplifying assumption that theexpansion of the universe is insignificant during the du-ration of the modulation, in which case the relation be-tween both time variables becomes linear and the cor-responding proportionality constant can be absorbed inthe modulation frequency.In contrast to the post-inflationary setting, where theproduced particles back-react on the inflaton, the mod-ulation of the coupling is imposed externally in the sim-ulation. Switching off the modulation thus correspondsto the decay of the inflaton, which constitutes a free pa-rameter in the model.A particularly useful observable for the study of para-metric resonance and turbulent thermalization is thesingle-particle momentum distribution f ( t, p ) = 1 V (cid:104) ˆΨ † p ( t ) ˆΨ p ( t ) (cid:105) , (20)where V is the volume and ˆΨ p ( t ) = (cid:82) d d x ˆΨ( t, x )e − i px / (cid:126) denotes the Fourier transform of the Bose field ˆΨ( t, x ) .Importantly, this quantity can be experimentally ac-cessed in time-of-flight measurements (see section V forfurther discussions).We evaluate the quantum expectation value in eq. (20)by means of classical-statistical (or truncated Wigner)simulations [28–31]. This method takes into accountquantum fluctuations by stochastically sampling classicalfield configurations from the Wigner distribution of theinitial state. Each realization is propagated deterministi-cally according to the Gross–Pitaevskii equation (GPE)and quantum mechanical observables are obtained as sta-tistical averages over multiple realizations. Here, the ini-tial state is taken to be a spatially homogeneous BEC,and for each realization, all non-zero momentum modes p/p L M o m e n t u m d i s t r i bu t i o n f ( t , p ) + (cid:15) p = ¯ hω/ (cid:15) p = 2 × ¯ hω/ t = 0 t = 30 × t t = 46 × t t = 58 × t t = 90 × t Figure 2. Radially averaged momentum distribution f ( t, p ) as a function of the radial momentum p = | p | at differenttimes t , demonstrating preheating dynamics. The coupling ismodulated with relative strength r = 0 . at a frequency ω ,chosen such that the resonance condition (cid:15) p res = (cid:126) ω/ for themomentum p res = 12 × p L with p L = 2 π (cid:126) /L is fulfilled (seediscussion below eq. (24)). At early times t (cid:46) × t , a sin-gle narrow resonance can be observed around p res . At latertimes, a broad resonance band emerges with peaks at higherharmonics of the modulation frequency. These secondary res-onances are due to non-linear interactions, as discussed insection III E. are populated with vacuum noise corresponding to an av-erage occupancy of half a particle per mode. This mimicsquantum fluctuations acting as a seed for parametric in-stabilities. It is important to note that this approachgoes beyond a mean-field description, which fails to cap-ture parametric resonance since the occupancies of allexcited modes are exactly zero. Further details aboutthe simulation method can be found in appendix A.In what follows, we express length and time inunits of the characteristic scales x = (cid:126) / √ mn g and t = (cid:126) /n g , respectively, where n = N/V is the ho-mogeneous particle density in a system of N particlesin a volume V . Momenta are given either in units ofthe lowest non-zero momentum p L = 2 π (cid:126) /L or in unitsof the characteristic momentum p ξ = 2 π (cid:126) /ξ correspond-ing to the healing length ξ = (cid:126) / √ mn g . In a quasi- DBose gas, the interaction strength g = ˜ g (cid:126) /m is charac-terized by the dimensionless parameter ˜ g = √ πa s /a HO ,where a s is the s -wave scattering length and a HO is theoscillator length of the harmonic potential in the stronglyconfined direction [57]. If not stated otherwise, we con-sider a uniform quasi- D system of N = 10 weakly inter-acting particles with coupling ˜ g = 2 . × − in a squarebox with periodic boundary conditions. This choice ofparameters fixes the box length as L/x = 50 . Moreover,in this isotropic setting, the momentum distribution onlydepends on the magnitude p = | p | of the momentum, f ( t, p ) ≡ f ( t, p ) .Figure 2 shows the radially averaged momentum dis-tribution of a parametrically excited system describedby eq. (13). The interaction is modulated according toeq. (19) with r = 0 . , and the modulation frequency ω is chosen as the resonance frequency of the momentum | p res | = 12 × p L , as discussed below eq. (24) in the nextsubsection. In this simulation, the noise cutoff has beenchosen as p Λ = p ξ ≈ . × p L .At early times, t (cid:46) × t , we observe a singlenarrow resonance around the resonance momentum p res satisfying the resonance condition (cid:15) p res = (cid:126) ω/ (cf. sec-tion III D). Due to particle number conservation, thegrowing occupancy of the resonant momentum causesthe condensate to decay, mimicking particle productionfrom the decaying inflaton in the early universe. At latertimes, secondary resonances at higher harmonics of themodulation frequency appear as a result of non-linear in-teractions among the produced quasi-particles.In addition to the narrow resonance peaks, a transientgrowth of fluctuations at low momenta occurs at earlytimes before the primary peak becomes visible. Thisgrowth can be interpreted as a consequence of the factthat the sampled initial state corresponds to the groundstate of an ideal Bose gas. At t = 0 , the system is effec-tively quenched to a finite interaction, producing a powerlaw in the momentum distribution proportional to p − atlow momenta (cf. appendix A 1). This early-time behav-ior has, however, no influence on the preheating dynamicswe are interested in here.In the following subsections, we provide analytical in-sights into both the linear regime of parametric resonanceand the non-linear regime of secondary resonances. D. The linear regime of parametric resonance
Parametric instabilities play a crucial role in manymodern experiments with BECs [34–36, 38–40] and canconveniently be understood by adopting a semi-classicalpoint of view. To this end, we consider the GPE for thecondensate wave function Ψ( t, x ) , i (cid:126) ∂ t Ψ = (cid:18) − (cid:126) ∇ m + g ( t ) | Ψ | (cid:19) Ψ , (21)which is formally obtained as the classical equation ofmotion after replacing the Bose field operator in eq. (13)by its expectation value, ˆΨ( t, x ) → Ψ( t, x ) = (cid:104) ˆΨ( t, x ) (cid:105) .It is instructive to work in the Madelung representation, Ψ( t, x ) = (cid:112) n ( t, x ) exp ( iθ ( t, x )) , which allows us to ex-press the GPE (21) in form of the hydrodynamic equa-tions [50] ∂ t n = − (cid:126) m ∇ ( n ∇ θ ) , (22a) (cid:126) ∂ t θ = − g ( t ) n + (cid:126) m (cid:20) ∇ √ n √ n − ( ∇ θ ) (cid:21) . (22b) To obtain some intuition about the early stages ofthe evolution, we express both the density and thephase in terms of a homogeneous background withfluctuations on top of it, n ( t, x ) = n ( t ) + n ( t, x ) and θ ( t, x ) = θ ( t ) + θ ( t, x ) . Linearizing eq. (22) with re-spect to the fluctuations, yields the equations ∂ t n = 0 and ∂ t θ = − g ( t ) n / (cid:126) for the background condensate.For the fluctuations, we obtain ∂ t n = − (cid:126) n m ∇ θ , (23a) (cid:126) ∂ t θ = − g ( t ) n + (cid:126) mn ∇ n . (23b)Taking the time derivative and inserting the resultingequations into each other, fluctuations of the density andthe phase decouple to linear order. Transforming to mo-mentum space, n p ( t ) = (cid:82) d d x n ( t, x )e − i px / (cid:126) and sim-ilarly for the phase, the linearized equations can be ex-pressed as ∂ t n p + ω p ( t ) n p = 0 , (24a) ∂ t θ p + 2 n ∂ t g ( t )2 n g ( t ) + (cid:15) p , ∂ t θ p + ω p ( t ) θ p = 0 , (24b)where (cid:126) ω p ( t ) = (cid:15) p , ( (cid:15) p , + 2 n g ( t )) is a time-dependentform of the famous Bogoliubov dispersion relation [58],and (cid:15) p , = p / m denotes the dispersion relation of a freeparticle.Equation (24) describes a collection of parametric os-cillators for each momentum mode p , which are un-damped for the density and damped for the phase. Asdiscussed in appendix B, these are special cases of Math-ieu’s equation [59], which admits oscillatory solutionswith exponentially growing amplitudes ∼ e ζ p t , describingparametric resonance [41, 56]. The resonance conditionfor the momentum mode p res reads (cid:15) p res = (cid:126) ω/ , where (cid:15) p = (cid:112) (cid:15) p , ( (cid:15) p , + 2 n g ) denotes the Bogoliubov disper-sion relation and ω is the modulation frequency defined ineq. (19). That is, resonance occurs for those momentummodes p res whose energy equals half a quantum of energy (cid:126) ω/ injected in the system through the modulation.In fact there is an entire range of modes around p res which experience a positive growth rate, and the width ofthis instability band increases with the modulation am-plitude r . To leading order in perturbation theory [60](see appendix B), the growth rate of the resonant mo-mentum mode is given by ζ p res = rω (cid:16) n g (cid:126) ω (cid:17) (cid:115) (cid:18) (cid:126) ω n g (cid:19) − . (25)For (cid:126) ω (cid:28) n g , this rate simplifies to ζ p res ≈ rω/ . Inthis regime, the Bogoliubov dispersion becomes linear, (cid:15) p ≈ c s | p | / (cid:126) with the speed of sound c s = (cid:112) n g /m , andthe produced quasi-particles have the character of soundwaves. In the opposite limit, (cid:126) ω (cid:29) n g , particles withquadratic dispersion (cid:15) p ≈ n g + p / m are produced.In the simulation presented in fig. 2, parametric reso-nance is clearly visible as a pronounced peak at the mo-mentum satisfying the resonance condition. It is worth-while emphasizing that parametric instabilities can onlybe triggered if the initial occupancy is non-zero. Thisseed is not contained in the mean field analysis presentedin this subsection, but is added in the simulation in formof vacuum noise according to the truncated Wigner pre-scription.The linearized equations 23 are helpful to get an in-tuitive analytical understanding for the early stages ofthe dynamics and describe the emergence of the primaryresonant peak in fig. 2. However, as a result of the ex-ponential growth of occupancies, this approach fails todescribe the later stages where non-linear effects play afundamental role. These non-linearities are taken intoaccount by our numerical simulations, which are basedon the full GPE (21), and include secondary excitationsoutside the resonance band, as shown in fig. 2. Thesewill be discussed further in the subsequent section. E. Secondary instabilities
In this section we discuss the first non-linear cor-rections, which lead to secondary instabilities [19, 41].Within our semi-classical picture, these can be under-stood by considering the hydrodynamic equations of mo-tion 22 for the fluctuations to quadratic order, i.e., in-cluding terms O ( θ ) , O ( θ ( n /n )) and O (( n /n ) ) .This leads to ∂ t n = − (cid:126) n m ∇ θ − (cid:126) m (cid:0) n ∇ θ + ∇ n ∇ θ (cid:1) , (26a) (cid:126) ∂ t θ = − gn + (cid:126) mn ∇ n − (cid:126) m (cid:34) ( ∇ θ ) + ( ∇ n ) n + n ∇ n n (cid:35) . (26b)As can be seen, density and phase fluctuations no longerdecouple to this order.Taking the time derivative of the above equation forthe density fluctuations and inserting the expressions forthe first-order time derivatives of the fluctuations, theresult can be written in momentum space as ∂ t n p + ω p ( t ) n p = n m (cid:126) (cid:90) q θ p − q θ q u ( p , q )+ (cid:90) q n p − q n q (cid:20) v ( p , q )8 m (cid:126) n − g ( t ) m (cid:126) pq (cid:21) (27)with u ( p , q ) = p ( pq − q ) + 2( pq )( p − q ) , v ( p , q ) = p (2 p − pq + q ) − pq ) q and (cid:82) q = (cid:82) d d q/ (2 π (cid:126) ) d .In deriving eq. (27), we have neglected all cubic terms inthe fluctuations.The above momentum integrals are dominated by thecontribution from the exponentially growing unstable modes q , for which | p − q | ≈ | q | ≈ p res . For modes thatare stable on the linear level, the integrals therefore actas source terms and eq. (27) describes forced harmonicoscillators with an exponentially growing force. This ef-fect is analyzed in more detail in appendix C. There,we show that the forcing leads to an exponential growthof the momentum modes with p (cid:46) p res proportional to e ζ p res t , with a growth rate ζ p res given by eq. (25). Modes,for which in addition (cid:15) p ≈ (cid:15) p res holds, experience a res-onant amplification and are strongest enhanced.These features of secondary instabilities are capturedby our numerical simulations. In particular, one can ob-serve both the narrow peak and the broad band in thedistribution function in fig. 2 between t = 46 × t and t = 58 × t . Similar peaks at higher multiples of the reso-nance frequency appear at later times due to higher-ordercorrections to eq. (22).The perturbative analysis breaks down when n /n ≈ , i.e., when the number of excited atomsbecomes comparable to the number of condensateatoms. At this point, the exponential growth stops andturbulent dynamics sets in, which is discussed in thenext section.In appendix D, we discuss the theory of parametric in-stabilities and secondary excitations from the perspectiveof quantum equations of motion for correlation functions,which validates the use of classical-statistical approxima-tion. F. Analysis of growth rates
Figure 3 depicts the time evolution of the occupanciescorresponding to the primary and secondary resonancesannotated in fig. 2. The growth of the secondary res-onance marks the onset of the non-linear regime wherethe quasi-particles produced in the primary resonancestart to interact, and the linear picture presented in sec-tion III D breaks down.Before comparing the numerically extracted growthrates to the analytical predictions presented in the pre-vious sections, we first relate the momentum distri-bution f ( t, p ) to the hydrodynamic density and phasevariables. On the mean-field level, linear fluctua-tions of the condensate wave function, expressed as Ψ( t, x ) = Ψ ( t ) + Ψ ( t, x ) , are related to linear den-sity and phase fluctuations via Ψ ( t, x ) / Ψ ( t ) = n ( t, x ) / n + iθ ( t, x ) . The momentum distribution onthe linear level then corresponds to | Ψ p | = n (cid:34) n p (2 n ) + θ p (cid:35) . (28)Therefore, a parametric resonance where the density andphase fluctuations grow as n p ∼ θ p ∼ e ζ p t results in agrowth of the momentum distribution as f ( t, p ) ∼ e ζ p t .The oscillations of the occupancies in fig. 3 can be un-derstood from the linearized parametric oscillator equa- t/t M o m e n t u m d i s t r i bu t i o n f ( t , p ) + (cid:15) ( p ) = ¯ hω/ ζ (1) = 0 . × t − (cid:15) ( p ) = 2 × ¯ hω/ ζ (2) = 0 . × t −
32 38 ∝ e ζ (1) t
46 49 52246 ∝ e ζ (2) t Figure 3. Radially averaged momentum distribution f ( t, p ) as a function of time t showing the exponential growth ofthe primary and secondary resonances corresponding to theannotated peaks in fig. 2. We have extracted the growthrates by fitting a straight line to the quantity log f ( t, p ) , asshown in the insets. The resulting growth rate of the primaryresonance ζ (1)num = 0 . × t − agrees well with the analyticalprediction ζ (1)pert = 0 . × t − obtained from eq. (25). Thegrowth of the secondary instability starts later, but its rate ζ (2)num = 0 . × t − is approximately twice as large as the oneof the primary resonance, as expected from our discussion insection III E. The exponential growth stops when the num-ber of excited atoms becomes comparable to the number ofcondensate atoms. tions 24. Recall that the latter admit oscillatory solu-tions with exponentially growing amplitudes. However,being conjugate variables, the oscillations of n and θ are shifted in phase by approximately π/ . According toeq. (28), the momentum distribution thus corresponds tothe sum of two phase-shifted oscillating functions withslightly different initial amplitudes. This results in theresidual oscillations on top of the exponential growth ob-served in fig. 3. We have checked that the oscillation fre-quency of the primary resonance agrees with the modula-tion frequency, while the oscillations of the secondary res-onance additionally contain frequency components cor-responding to twice the modulation frequency, reflectingthe interactions between the resonantly produced quasi-particles.We have extracted the growth rate of the primary res-onance ζ (1)num = 0 . × t − by fitting an exponential func-tion to the numerical data, as shown in the insets offig. 3. The result is close to the analytical prediction (25)obtained from perturbation theory, ζ (1)pert = 0 . × t − .The secondary resonance at (cid:15) p res in fig. 3 grows at therate ζ (2)num = 0 . × t − , which is indeed approximatelytwice the growth rate of the primary resonance. The sec-ondary resonance quickly saturates, indicating the break-down of the perturbative analysis presented in the pre- vious section. This marks the transition into a highlynon-linear regime leading to the formation of turbulence,which is the topic of the next section. IV. ANALOG REHEATING
Having studied the main features of preheating dynam-ics and its analog implementation in an ultracold Bosegas, we now turn to the reheating process and the ques-tion of how the system finally approaches thermal equi-librium.Particle spectra formed during preheating are highlynon-thermal with large occupation numbers at low mo-menta. As a consequence, the system enters a turbu-lent state, characterized by a local transport of conservedquantities in momentum space. Typically, an inverse cas-cade transports particles towards lower momenta, whilea direct cascade transports energy towards higher mo-menta, constituting a key process in the context of tur-bulent thermalization [7]. At early times, this transportis driven, i.e., the oscillating inflaton acts as a sourceinjecting energy into the system at resonant momenta.Eventually, the inflaton decays, marking the transitionfrom driven to free turbulence. The system remains inthe turbulent state for a long time, until the occupancy ofcharacteristic momenta eventually becomes comparableto the vacuum expectation value given by the “quantumhalf”. In this final stage, which is dominated by quantumfluctuations, the system relaxes to thermal equilibrium,completing the reheating process [6, 7].Turbulent dynamics is accompanied by the emergenceof self-similarity and universality. This is reflected by apower-law behavior of the momentum distribution withina certain inertial range of momenta. One well-known ex-ample within the theory of weak wave turbulence is theprediction of a stationary direct cascade with a univer-sal power-law distribution f ( p ) ∝ p − d [61]. More gener-ally, self-similarity in far-from-equilibrium systems canbecome manifest in their time evolution. In such a sce-nario the evolution of the distribution function can beexpressed as f ( t, p ) = s α f S ( s β p ) , (29)where s = t/t ref and t ref is an arbitrary reference time.The scaling hypothesis (29) constitutes a significant re-duction of complexity as it allows to describe a rele-vant part of the system’s dynamics by simply rescaling asingle-variable scaling function f S as determined by thescaling exponents α and β . Remarkably, in many far-from-equilibrium scenarios, both the exponents as well asthe scaling form of the distribution are universal, whichmeans they are insensitive to microscopic details as wellas initial conditions, and depend only on a few macro-scopic system parameters like dimensionality, symmetryor the number of field components [7, 8, 10, 43, 62].As a result, if universality holds, physical systems withvastly differing energy scales can behave quantitatively0the same. This makes ultracold Bose gases a particu-larly promising target for simulating universal aspects offar-from-equilibrium dynamics like that of the inflaton inthe early universe. Universal self-similar time evolutionreflects the system being in the vicinity of a non-thermalfixed point, which acts as an attractor on the way towardsthermal equilibrium [42, 63, 64].In this section, we demonstrate by means of classical-statistical simulations that the salient features of thesephenomena can be observed in our implementation of re-heating dynamics. We first study the regimes of drivenand free turbulence separately, before considering a tran-sient prescaling regime, where the universal shape of thescaling function is maintained, but the scaling exponentschange over time. From our numerical simulations, weextract the scaling exponents as well as the scaling formof the distribution and compare the results with analyt-ical predictions from kinetic theory. Finally, we discussthe relaxation to thermal equilibrium, which is, however,not captured by our classical-statistical simulations as itis dominated by quantum fluctuations. We conclude thissection with a discussion of an example of how expan-sion may prevent the system from thermalizing, buildinga bridge to section II. A. Driven versus free turbulence
In order to drive the system into a turbulent state,we follow the protocol presented in section III C of para-metrically exciting a homogeneous BEC. Here, our fo-cus lies on the later stages of the non-linear dynamicsafter the proliferation of secondary instabilities when asmooth distribution in form of a power law has formed.We distinguish the regime of driven turbulence, realizedby continuously modulating the interaction according toeq. (19), and free turbulence, developed if the modulationis switched off shortly after the primary resonance hassaturated. In our analogy to reheating in the early uni-verse, the former case corresponds to the situation wherethe inflaton possesses enough energy to drive turbulencefor a long time, while in the latter case, the inflaton runsout of energy rather quickly at around the same timewhen turbulence sets in.To maximize the inertial range where self-similar scal-ing can be observed, it is desirable to inject energy at mo-mentum scales close to the lowest momenta supported bythe system where occupancies can become large. To thisend, the interaction is modulated with a relative ampli-tude r = 1 at frequency ω chosen such that the resonantmomentum becomes p res = 3 × π (cid:126) /L (cf. section III D).The left and right panels of fig. 4 show a compari-son of the direct cascades emerging in the regimes ofdriven and free turbulence, respectively. All parame-ters for both simulations are identical, with the exceptionthat in fig. 4a, the interaction is modulated continuously, while in fig. 4b, the modulation is switched off smoothly within two modulation periods π/ω ≈ . × t at time t = 80 × t , roughly corresponding to the time when theprimary resonance saturates. In the latter case, energyis conserved already at the onset of turbulence and thedriven regime is skipped. In both scenarios, the momen-tum distribution takes a scaling form corresponding toa power law close to f ( t, p ) ∝ p − , which is indicated bya dotted line as a guide to the eye. A power law pro-portional to p − d is characteristic for weak wave turbu-lence [61, 65, 66] and has been observed experimentallyin ref. [67].Moreover, as it can be seen in fig. 4, the distributions inthe two regimes exhibit self-similar time evolution in dif-ferent ways. In the case of driven turbulence, fig. 4a, thefront of the cascade evolves self-similarly, leaving behinda stationary distribution. Stationary turbulence arisesin the theory of weak wave turbulence as a stationarysolution of the scattering integral. However, such a con-figuration necessarily requires the existence of at least asource and a sink for energy to be injected and dissipated,respectively [61]. In the present case, the energy source isprovided by the modulation of the interaction. We haveverified that the rate of energy injection into the systemis approximately constant in this stationary regime, asexpected for driven turbulence [7]. Since energy is trans-ported locally in momentum space, unoccupied highermomentum modes play the role of an energy sink [7].This allows for the build-up of a stationary distribution,although the model lacks a mechanism of dissipation. Bycontrast, the distribution in the case of free turbulence,fig. 4b, is not stationary, but decreases as a function oftime for a given momentum, reflecting energy conserva-tion.To quantify the self-similar time evolution, we have ex-tracted the scaling exponents α and β defined in eq. (29)from our numerical data at different times using the max-imum likelihood technique described in appendix E. Inthe turbulent stage of the dynamics, the exponents areapproximately constant at late times, as shown in thelower part of fig. 4. The observed slow relaxation of theexponents at early times in fig. 4b can be interpreted asprescaling, as discussed in the subsequent subsection. Atthe latest simulated times, the exponents take the values α driven = − . ± . , β driven = − . ± . (30)for driven turbulence and α free = − . ± . , β free = − . ± . (31)for free turbulence, respectively. Rescaling the distribu-tion according to eq. (29) with the extracted exponents, We have checked that the differences between switching off themodulation suddenly or smoothly within a few modulation peri-ods are insignificant. − Radial momentum p/p ξ M o m e n t u m d i s t r i bu t i o n f ( t , p ) + ∝ p − t = 1024 × t t = 544 × t t = 304 × t t = 184 × t t = 124 × t t = 63 × t Rescaleddistribution123 α / β
200 400 600 800 1000Time t/t . . . − β (a) Direct cascade with continuously modulated interaction. − Radial momentum p/p ξ M o m e n t u m d i s t r i bu t i o n f ( t , p ) + ∝ p − t = 1024 × t t = 544 × t t = 304 × t t = 184 × t t = 124 × t t = 63 × t Rescaleddistribution24 α / β
200 400 600 800 1000Time t/t . . . − β (b) Direct cascade after switching off the modulation. Figure 4. Self-similar time evolution of the momentum distribution in form of a direct energy cascade for driven (a) and freeturbulence (b). Energy is injected at low momenta by modulating the scattering length according to eq. (19) with a relativeamplitude r = 1 at a frequency ω chosen such that (cid:126) ω/ (cid:15) p res with p res = 3 × π (cid:126) /L . In the case of continuous modulation(a), a stationary distribution with a power law close to p − develops, whose front is evolving self-similarly. If the driving isswitched off once the primary resonance has saturated, corresponding here to t = 80 × t (b), energy is propagated in a self-similar way to higher momenta, but the distribution at a given momentum decreases with time, reflecting energy conservation.A power law proportional to p − is shown in form of a dotted line as a guide to the eye. The insets show the distributionsrescaled according to eq. (29) using the numerically extracted scaling exponents displayed below the respective distributions. as described in appendix E, all data points collapse to asingle universal scaling function f S ( p ) , as shown in theinsets of fig. 4. The exponent β describes the speed ofenergy propagation towards higher momenta, which ishigher for the driven cascade than for the free cascade.Analytical predictions for this exponent from kinetic the-ory as well as the relations between α and β followingfrom conservation laws will be discussed in section IV C.The above values of the exponents are insensitive tothe details of how the far-from-equilibrium state is ap-proached. In particular, we have observed the same ex-ponents starting from an initial state with a highly occu-pied narrow window of momenta on top of a condensatebackground. Such an initial state is similar to the stateof the system at the end of the preheating stage, whena certain momentum mode is overpopulated as a conse-quence of parametric resonance.The slow power-law dynamics of the direct cascade can be challenging to capture with classical-statistical simu-lations for experimentally realistic configurations, sincethis method is known to be prone to instabilities causedby the vacuum noise [31]. These instabilities manifestin a decay of the “quantum half” and the generation ofspurious quantum pressure, resulting in an unphysicaldependence of the results on the ultraviolet (UV) cut-off [48]. Classical statistical simulations are therefore re-stricted to large occupancies and weak couplings, wherethese inevitable deficiencies are mitigated via a separa-tion of scales. To simulate reheating dynamics in theturbulent regime, we have therefore increased the particlenumber and reduced the coupling such that the validityof the classical-statistical approximation can be ensured.This is discussed in detail in appendix A 3, where we alsoassess the range of accessible coupling strengths for oursetup. In section V, we discuss how our numerical resultsrelate to realistic experimental conditions.2 B. Prescaling
So far, we have investigated the two limiting caseswhere turbulence is either driven or free. Now, we ad-dress the transient regime corresponding to the somewhatmore realistic situation where, in the beginning, turbu-lence is driven by the inflaton oscillations, but at somepoint goes over to free turbulence when the inflaton hasdecayed.This scenario is illustrated in fig. 5. Up to the time t = 256 × t , the direct cascade is driven (blue curves),as in fig. 4a. At this time, the modulation is switchedoff, mimicking the decay of the inflaton. We then ob-serve a slowing down of the speed of energy propagationand the distribution decreases in time for a given mo-mentum (red curves), reminiscent of the direct cascadeof free turbulence shown in fig. 4b.On the right-hand side of fig. 5, the numerically ex-tracted scaling exponents are shown as a function of time.After switching off the modulation, the exponent β slowlyevolves from a value close to the one reported in eq. (30)for driven turbulence to a value close to the one reportedin eq. (31) for free turbulence. By contrast, the ratio α/β changes rather quickly between these two regimes.This behavior is expected since this ratio is fixed by en-ergy conservation, which is enforced instantaneously af-ter switching off the modulation, cf. section IV C. Sur-prisingly, the self-similar scaling form of the distributionis approximately preserved during the transition, whichis clearly visible in the inset of fig. 5, where the distri-butions, rescaled according to eq. (29) with the time-dependent scaling exponents, fall on top of each other.Recently, such a situation, where the system’s dynam-ics is governed by a universal scaling function much be-fore the corresponding exponents have attained their uni-versal values, has been studied in the context of heavy-ioncollisions [46]. This phenomenon, termed prescaling, isclosely related to the emergence of a far-from-equilibriumhydrodynamic behavior, as it allows to describe the dy-namics in terms of few slowly-changing parameters [46].A different type of prescaling, where certain correlationfunctions already scale with their universal exponents atearly times while others do so only at much later times,has been studied numerically in three-component Bosegases [47]. Our results indicate that prescaling, as de-fined in ref. [46], can be observed during the transitionfrom driven to free turbulence , opening up new paths tostudying this phenomenon experimentally.In appendix F, we corroborate our analysis of prescal-ing using an alternative method based on the momentsof the distribution, which is particularly suitable for ex-tracting time-dependent scaling exponents [46]. A much shorter and less pronounced stage of prescaling occursalso in the scenario shown in fig. 4b, where the modulation isswitched off much earlier, at t = 80 × t . C. Kinetic description
We now discuss how the values of the scaling exponentsfor the direct cascade presented in eqs. (30) and (31) canbe understood from an analytical point of view.One relation between the exponents α and β followsfrom the scaling of the total energy. To see this, weassume a self-similar time evolution according to eq. (29)as well as a power-law scaling of the dispersion relation, (cid:15) p ∝ | p | z , where z ≈ in the particle regime of theBogoliubov dispersion. Then, the total kinetic energyscales as E ( t ) = V (cid:90) d d p (2 π (cid:126) ) d (cid:15) p f ( t, p ) = (cid:18) tt ref (cid:19) α − ( d + z ) β E ( t ref ) . (32)In the case of free turbulence, energy is conserved acrossthe cascade, E = const , while in the case of driven tur-bulence, it grows linearly with time, E ∝ t . This impliesthe relation α = γ + ( d + z ) β , (33)where γ = 1 in the driven regime and γ = 0 in theabsence of driving.Another relation between the scaling exponents can beobtained with the help of perturbative kinetic theory forBose gases [10, 65]. The time-evolution of the distribu-tion function in this framework is written in the form ofa Boltzmann equation (see, e.g., ref. [68]), ∂ t f ( t, p ) = C ↔ [ f ]( t, p ) + C ↔ [ f ]( t, p ) , (34)where C ↔ [ f ] = 4 πg (cid:90) p , p , p δ ( p + p − p − p ) × δ ( (cid:15) p + (cid:15) p − (cid:15) p − (cid:15) p ) × [( f p + 1)( f p + 1) f p f p − f p f p ( f p + 1)( f p + 1)] (35)describes ↔ collisions between non-condensate atomsand C ↔ [ f ] = 2(2 π ) (cid:126) n g (cid:90) p , p , p δ ( p − p − p ) × δ ( (cid:15) + (cid:15) p − (cid:15) p − (cid:15) p ) × [ δ ( p − p ) − δ ( p − p ) − δ ( p − p )] × [( f p + 1) f p f p − f p ( f p + 1)( f p + 1)] (36) Large occupation numbers usually require going beyond pertur-bative kinetic theory for the description of the dynamics neara non-thermal fixed point. An example is the dynamics of theinverse cascade leading to Bose condensation [10]. For the directcascade, however, occupation numbers of characteristic momentaare typically lower (although still much larger than one), suchthat a perturbative analysis is applicable (see also [7]). − Radial momentum p/p ξ M o m e n t u m d i s t r i bu t i o n f ( t , p ) + t = 1280 × t t = 932 × t t = 682 × t t = 505 × t t = 379 × t t = 288 × t t = 223 × t t = 177 × t t = 143 × t t = 63 × t Rescaleddistribution 200 400 600 800 1000 1200Time t/t . . . − β
200 400 600 800 1000 1200Time t/t α / β Figure 5. Prescaling at the transition from driven to free turbulence. The simulation parameters are identical to those infig. 4, but the modulation is switched off suddenly at a later time t = 256 × t . Before this time (blue curves), turbulenceis driven and both the momentum distribution as well as the scaling exponents α and β are the same as in fig. 4a. Afterswitching off the modulation (red curves), the ratio α/β quickly changes to the one expected for free turbulence, reflectingenergy conservation. The exponent β gradually changes towards the value obtained for free turbulence in fig. 4b, reducing thespeed of energy transport in the cascade. Although the scaling exponents still change in time, the distribution has alreadyattained its universal scaling form. This important hallmark of prescaling is indicated in the inset, where all data points collapseto a single curve after rescaling according to eq. (29) with the extracted time-dependent scaling exponents α ( t ) and β ( t ) , asdescribed in appendix E. arises in the presence of the condensate and describeseffective ↔ collisions between the non-condensateatoms. Here, (cid:15) and (cid:15) p are the energies of condensate andnon-condensate atoms, respectively, and we have usedthe abbreviations f p = f ( t, p ) and (cid:82) p = (cid:82) d d p/ (2 π (cid:126) ) d . Inour case, (cid:15) can be neglected since we are interested inthe direct cascade, which predominantly involves atomsthat are much more energetic than condensate atoms.In order to derive the desired relation between the scal-ing exponents, let us assume that the right-hand side ofthe Boltzmann equation (34) is dominated by one of thetwo scattering integrals, and denote the number of par-ticles involved in the corresponding scattering processesby l , i.e., we have l = 4 for ↔ scattering and l = 3 for ↔ scattering. Comparing how the left- and right-hand sides of eq. (34) scale with time, and assuming largeoccupation numbers, f p (cid:29) , one arrives at [10] ( l − α = [( l − d − β − . (37)Combining this result with eq. (33) yields β = − − − ( l − γ ( l − z + 2 . (38)For driven turbulence, γ = 1 , and a quadratic disper-sion, z = 2 , both collision terms lead to α = − and β = − / . These values are close to those in eq. (30),extracted from our simulation in the presence of a mod-ulated coupling. For free turbulence, γ = 0 , one recoversthe ratio α/β = 4 reflecting energy conservation, which is close to the numerical result in eq. (31). The valueof β in this case is β = − / for ↔ scattering and β = − / for ↔ scattering. The latter is closer toeq. (31), which is observed at late times in our simula-tions. By contrast, in the experiment reported in ref. [45],where, instead of parametrically driving a pure conden-sate out of equilibrium, a cooling quench has been appliedto an initially uncondensed Bose gas, a scaling exponentcloser to the ↔ scattering solution, β = − / , hasbeen observed. Deviations from the predicted scaling canin general be induced by a non-quadratic scaling of thedispersion relation (cid:15) p or by the time-dependence of thenumber of atoms in the condensate. D. Thermalization
As we have seen, on its way to thermal equilibriumthe system takes a detour via a non-thermal fixed point,in the vicinity of which the dynamics is dominated by aturbulent transport of energy towards higher momenta.The self-similar dynamics of the direct cascade stops oncethe occupancy of the characteristic momentum dominat-ing the kinetic energy budget becomes comparable to theexpectation value of the vacuum noise given by the “quan-tum half”. At this point, quantum fluctuations becomedominant over statistical fluctuations [7] and the systemis expected to relax to a Bose–Einstein distribution f BE ( p ) = 1exp( (cid:15) p /k B T f ) − , (39)4where (cid:15) p is the dispersion relation, k B is the Boltzmannconstant and T f is the reheating temperature.The dominance of quantum fluctuations makes thisfinal stage of the dynamics inaccessible to classical-statistical simulations. In contrast to the expected re-laxation to a Bose–Einstein distribution, at sufficientlylate times in the numerical simulation, the cascade stopsbeing self-similar, slows down, and approaches a classi-cal equilibrium distribution [8] with a temperature (cid:101) T (Λ) that depends on the UV cutoff Λ and is determined bythe equipartition theorem, f cl , th ( p ) + 1 / ∝ k B (cid:101) T (Λ) /(cid:15) p .In particular, mode occupancies drop to unphysical val-ues below the vacuum noise of / .The late-time regime dominated by quantum fluctua-tions is beyond the regime of classical-statistical simula-tion methods. Here, we resort to analytical estimates forthe reheating time and the reheating temperature basedon the system’s self-similar time evolution. Followingref. [7], we neglect the final stage of quantum relaxationto a Bose–Einstein distribution and consider thermaliza-tion as complete once the occupancy of the character-istic momentum dominating the kinetic energy budgetbecomes on the order of unity. More precisely, we de-fine the characteristic momentum ¯ p ( t ) as the momentumthat maximizes the integrand in the expression (32) forthe total kinetic energy. For an isotropic system with dis-persion relation (cid:15) p ∝ p z , where p = | p | , the characteristicmomentum is given by ¯ p ( t ) = arg max p p d + z − f ( t, p ) .If the scaling exponents as well as the momentum dis-tribution at some reference time are known, the assump-tion of self-similar time evolution according to eq. (29)is sufficient to predict the time t f when the occupancyof the final characteristic momentum p f will reach unity.The time t f can be regarded as the best possible approx-imation to the reheating time obtainable from classical-statistical simulations. However, eq. (29) is not directlyapplicable since the scaling exponents are not constantthroughout the entire time evolution. In fact, they takedifferent universal values in the regimes of driven and freeturbulence, respectively, which are interpolated during atransient regime of prescaling. Importantly, as shownin section IV B, the universal scaling form of the distri-bution is preserved during prescaling, allowing us to de-scribe the full evolution solely in terms of time-dependentscaling exponents.To properly account for time-dependent exponents α ( t ) and β ( t ) , we generalize the power-law expressions ineq. (29) as ( t/t ref ) α → exp (cid:82) tt ref α ( t (cid:48) ) d t (cid:48) /t (cid:48) , and similarlyfor β [46]. Note that the original expressions are recov-ered for constant α and β . This way, the evolution of theinitial characteristic momentum scale for turbulence p i isgiven by ¯ p ( t ) = exp (cid:20) − (cid:90) tt i β ( t (cid:48) ) d t (cid:48) t (cid:48) (cid:21) p i , (40)where t i denotes the time when turbulence sets in. Whilethe exponent β determines the scaling of the characteris- tic momentum, the evolution of its occupancy is governedby the exponent α , f ( t, ¯ p ( t )) = exp (cid:20)(cid:90) tt i α ( t (cid:48) ) d t (cid:48) t (cid:48) (cid:21) f ( t i , p i ) . (41)By virtue of eq. (41), it is in principle possible to com-pute the time t f by solving f ( t f , p f ) = 1 , given full knowl-edge of the time-dependence of the exponent α . It isinstructive, however, to rewrite this formula using somesimplifying approximations. To this end, we assume that α is approximately constant during the driven regime,taking its universal value α driven , and that, after switch-ing off the modulation at time t d , the relaxation to itsuniversal value for free turbulence α free occurs fast com-pared to the overall reheating time scale. We then recoverpower-law scaling for the occupancy of the final charac-teristic momentum, f ( t f , p f ) ≈ (cid:18) t d t i (cid:19) α driven (cid:18) t f t d (cid:19) α free f ( t i , p i ) , (42)and it reaches unity at the time t f = t d (cid:18) t d t i (cid:19) − α driven /α free [ f ( t i , p i )] − /α free . (43)Recall that a direct energy cascade is described bynegative exponents α and β with large initial occu-pancies f ( t i , p i ) , which decrease in time according toeq. (42). Furthermore, the condition t i ≤ t d ≤ t f requires t d ≤ t i [ f ( t i , p i )] − /α driven . If the latter inequality is sat-isfied as an equality, the occupancy of the characteris-tic momentum reaches unity already during driven tur-bulence before switching off the modulation, such that t f = t d . Here, we focus on the situation where t d < t f such that the system spends a dominant part of its evo-lution in the regime of free turbulence.Using an analogous line of arguments, we can estimatethe final characteristic momentum scale from eq. (40) as p f ≈ (cid:18) t d t i (cid:19) − β driven (cid:18) t f t d (cid:19) − β free p i = (cid:18) t d t i (cid:19) − β driven + α driven / ( d + z ) [ f ( t i , p i )] / ( d + z ) p i , (44)where the second equality follows after inserting our es-timate for the reheating time (43) and substituting therelation (33) between α free and β free imposed by energyconservation. The use of the latter identity makes theestimate (44) of the final characteristic momentum in-dependent of the scaling exponents in the regime of freeturbulence, reflecting the fact that p f is closely connectedto the total energy in the system [7], which is conservedafter switching off the modulation. By contrast, the re-heating time, according to eq. (43), is sensitive to thevalues of the scaling exponents in both regimes, and, inparticular, can be influenced by the non-universal behav-ior during prescaling.5Finally, we can obtain an order-of-magnitude estimatefor the reheating temperature T f by identifying the lat-ter with the typical kinetic energy of a particle with mo-mentum p f . For a dispersion relation of non-relativisticparticle, this reads k B T f ∼ p m . (45)When interpreting the estimates (43) to (45) in the cos-mological context, it is important to remember that ourtime variable t denotes the lab time and should be trans-formed back to the original cosmic time τ by integratingthe relation d t = d τ /a . Similarly, we need to take intoaccount the redshift of momenta, p → p/a , as expressedin the kinetic term of eq. (11). For the estimate of thereheating temperature (45), this means T f → T f /a ( t f ) .This back transformation of variables is illustrated fur-ther in the subsequent subsection at a specific example.Although eqs. (43) to (45) are useful to estimateasymptotic quantities without requiring to simulate thedynamics up to the point where the system thermalizes,the outlined argumentation is based on the strong as-sumption that the neglected final stage of the dynamics,which is dominated by quantum fluctuations, does nothave a significant impact on these estimates. The latterremains to be checked against physical reality, which canonly be provided by comparison to an experiment. E. Thermalization versus expansion
In an expanding spacetime, interaction rates decreasewith the expansion. If the expansion is too rapid, it istherefore possible that the system is unable to ever reachthermal equilibrium. This ability is usually quantified bycomparing the typical interaction rate Γ with the Hub-ble parameter H = ˙ a/a . If Γ (cid:29) H , thermalization ispossible, whereas for Γ (cid:28) H , the occupation numbersare expected to freeze without thermalizing as the meanfree path of the particles exceeds the horizon size of theuniverse.The full interpolation between both limiting cases, andthe associated thermalization behavior, can be studied inour proposed setup. According to our mapping eq. (13),the dynamics of the Bose gas in an expanding space-time for d = 2 is mapped into that of a static spacetimewith constant interaction g . Consequently, the simulat-ing experimental system will always thermalize eventu-ally. However, this does not necessarily imply the samefor the simulated system, whose description is recoveredonly after a back transformation to the original cosmicspatial and temporal coordinates as well as field variables.A necessary condition for thermalization to occur in thesimulated system is that the lab time t at which the cos-mic time τ becomes infinite is larger than the reheatingtime, t f < t ( τ = ∞ ) .An interesting feature of our mapping is that in thecase of a sufficiently rapid expansion, infinite cosmic − Radial momentum p/p ξ M o m e n t u m d i s t r i bu t i o n f ( τ , p ) + τ = 250 × t τ = 500 × t τ = 1000 × t τ = 2000 × t τ = 4000 × t τ = 8000 × t τ = 86 × t Lab time t/t C o s m i c t i m e τ / t Figure 6. Freeze-out of the momentum distribution of a sys-tem with matter-dominated expansion according to eq. (46)with H = 0 . × t − and ν = 2 / . The numerical data isthe same as in fig. 5, but the lab time t has been transformedback to the cosmic time τ by virtue of eq. (48). As shownin the inset, the cosmic time diverges at the finite lab time t ( τ = ∞ ) ≈ × t , which is, in particular, shorter than thereheating time of the associated simulating system. As a re-sult, the evolution of the simulated expanding system slowsdown and eventually freezes before thermalization is com-plete. time corresponds to a finite value of the lab time, i.e., t ( τ = ∞ ) < ∞ . To demonstrate this, we consider apower-law expansion [69], a ( τ ) = (cid:18) H τν (cid:19) ν , (46)with a positive exponent ν . This expression fixes a (0) = 1 and the initial Hubble parameter H (0) = H . The labtime is then given by t ( τ ) = (cid:90) τ d t (cid:48) a ( t (cid:48) ) = 1 H ν − ν (cid:104) a ( τ ) (1 − ν ) /ν − (cid:105) . (47)If ν > / , such as in the case of expansion in a matter-dominated universe ( ν = 2 / ) [49], the limit a = ∞ and τ = ∞ is reached within a finite lab time, given by t ( τ = ∞ ) = ν (2 ν − H , (48)in line with the intuition that a larger Hubble parameterleads to a faster freeze-out.Remarkably, eq. (48) associates the state of the sim-ulating system at each instant of the lab time t withthe asymptotic infinite-time state of a simulated system6characterized by a certain value of the initial expansionrate H . This is illustrated in fig. 6, where we used H = 0 . × t − and ν = 2 / , leading to the freeze-outtime t ( τ = ∞ ) ≈ × t .The ability to associate the dynamics of a single simu-lation with arbitrary expansion schemes is a special caseof the D geometry, where, owed to scale invariance, theequations of motion become independent of the scale fac-tor a . The situation is different in three spatial dimen-sions, where, according to eq. (14), the coupling becomestime-dependent after the mapping as g ∝ a − . In thepresence of a decreasing coupling, freeze-out may alreadyoccur during the simulation before reaching t ( τ = ∞ ) , inwhich case the expanding system will never thermalize. V. EXPERIMENTAL PERSPECTIVES
In this section, we summarize the experimental require-ments for observing the salient features of reheating dy-namics discussed in this paper.The simulating system is a single-component BECthat is parametrically excited by modulating the inter-action according to eq. (19) around some positive off-set value. This can be realized experimentally with thehelp of a Feshbach resonance, which allows one to tunethe s -wave scattering length using an external magneticfield [22, 34, 36, 39, 40]. The interaction strength g ofa D Bose gas is related to the s -wave scattering length a s via g = 4 π (cid:126) a s /m . Our general approach is inde-pendent of dimensionality (but see our remarks towardsthe end of this section about the absence of a turbulentcascade in D). For concreteness, we have focused hereon D geometries, which may be realized through a tightconfinement of the atomic gas in the vertical direction bya harmonic potential. In the quasi- D regime, the effec-tive interaction strength g = ˜ g (cid:126) /m is characterizedby the dimensionless parameter ˜ g = √ πa s /a HO , where a HO = (cid:112) (cid:126) /mω HO is the oscillator length of the confiningharmonic potential with frequency ω HO [57]. Besides viaFeshbach resonances, D Bose gases with tunable effec-tive interactions can therefore alternatively be realizedby changing the frequency of the harmonic trap [70].The numerical simulations presented in this paper havebeen conducted for the fixed value N ˜ g = 2 . × of theproduct of the particle number N and the interactionstrength ˜ g . This value is readily achievable in present-day experiments with uniform D Bose gases [27]. Onthe classical level, the parameters N and ˜ g enter theequations of motion only via the product N ˜ g (cf. ap-pendix A 2), but this no longer holds on the quantumlevel. In fact, our classical-statistical simulations are sen-sitive to the individual values of N and ˜ g since they con-trol the relative magnitude of mode occupancies withrespect to quantum fluctuations mimicked in form ofvacuum noise. As discussed in appendix A 3, classical-statistical simulations are restricted to large occupanciesand weak interactions, which is why our simulations of turbulent reheating dynamics have been conducted fora higher particle number than currently realizable in ex-periments. To estimate characteristic quantities for morerealistic experimental setups, we rescale our numericalresults to the reference parameters N ref = 1 × and ˜ g ref = 2 . × − . Pure BECs with almost atomscan be reached in state-of-the-art experimental setups,e.g., for K [71–73], where a wide range of interactionsis accessible through the broad Feshbach resonance atthe magnetic flux B = 560 . [74]. We emphasize thatextrapolating weak coupling results to stronger couplingsis a priori not justified. Yet, this procedure is often theonly way to be consistent with experimental aspects andcommonly used, e.g., in the context of quantum chromo-dynamics (QCD) [75]. In the field of ultracold atoms, anumber of positive examples exist, where weak couplingexpectations of non-equilibrium quantum dynamics havebeen found experimentally at strong couplings [11, 76].For the above choice N ˜ g = 2 . × and a D systemof K atoms in a square box of length L = 50 µ m , oneobtains the typical values x = 1 µ m and t ≈ . forthe characteristic length and time scales x = (cid:126) / √ mn g and t = (cid:126) /n g , respectively. Here, g is the offsetvalue around which the interaction is modulated and n = N/L is the atomic density.As shown in the previous sections, the reheating dy-namics studied in this work can be subdivided into twodistinct stages, namely, the early stage of preheating,where occupancies of resonant momenta grow exponen-tially as a result of parametric instabilities, and the sub-sequent stage of turbulent thermalization, characterizedby a self-similar transport of energy towards higher mo-menta. The time scales at which these phenomena canbe observed depend on the choice of the modulation fre-quency ω and the relative modulation amplitude r ineq. (19). The choice r = 1 , i.e., modulating the scatter-ing length with an amplitude close to its offset value, isefficient for rapidly driving the system out of equilibrium.(Note that the perturbative expression for the growthrate eq. (25) is no longer applicable in this case.) Thedirect cascade is observed best if occupancies of the ini-tial characteristic momentum, where energy is injectedin the system, are as high as possible. To this end,the modulation frequency ω should be chosen close tothe frequency of the lowest non-zero momentum modesupported by the system, which is typically deep in thephonon regime where the Bogoliubov dispersion is lin-ear. In our simulations of reheating dynamics presentedin figs. 4 and 5, we have chosen the modulation frequencysuch that the primary resonance occurs at the momen-tum p res = 3 × π (cid:126) /L , i.e., ω = 2 (cid:15) p res / (cid:126) ≈ π ×
199 Hz for the example parameters mentioned above. The onsetof turbulence, marking the end of the preheating stage,occurs at around t ≈ × t , corresponding to roughly oscillations of the scattering length. The subsequentturbulent scaling regime extends up to times of about t f ≈ × t , which can be deduced from figs. 4 and 5with the help of eq. (43) after rescaling the occupancies7to the reference particle number N = 1 × (cf. ap-pendix A 3). The estimated reheating time t f is wellwithin the reach of modern experiments, where typicallifetimes of BECs are on the order of a few seconds [72].In our analysis, we have considered the momentum dis-tribution f ( t, p ) = (cid:104) ˆΨ † p ( t ) ˆΨ p ( t ) (cid:105) as the main observable.Experimentally, this quantity can be accessed in time-of-flight (TOF) measurements [57]. In the case of driventurbulence (cf. fig. 4a), by virtue of eq. (44), we find thatthe characteristic momentum at the reheating time t f isroughly given by p f ≈ . × p ξ , which is on the order ofthe momentum corresponding to the healing length ξ . Ingeneral, this value depends on the duration of the regimeof driven turbulence and tends to be smaller in the case offree turbulence (cf. fig. 4b) or in the prescaling scenario(cf. fig. 5), where less energy is injected into the sys-tem. It is a fortunate circumstance that higher momenta,which are increasingly populated as the direct cascadeprogresses, are typically easier to resolve in TOF mea-surements since they require a shorter expansion time.As an alternative means of detection, in situ images ofthe density profile may be used. We have checked that allrelevant signatures of both the parametric resonance andthe turbulent cascade can be extracted from the quantity (cid:104) ˆ n † p ( t )ˆ n p ( t ) (cid:105) , where ˆ n p ( t ) is the Fourier transform of theoperator ˆ n ( t, x ) = ˆΨ † ( t, x ) ˆΨ( t, x ) describing the spatialdensity profile of the Bose gas.Although in this paper we have put our main focus onuniform systems, nowadays available in many laborato-ries [70, 77], we expect the characteristic features of thedynamics to be prevalent even in the presence of a har-monic trap. In general, such an external trapping poten-tial couples different momenta and thereby changes thenature of Bogoliubov quasi-particles. However, withina local density approximation, a harmonically trappedBEC in the Thomas–Fermi regime can be regarded aslocally uniform around the center [50]. Thus, while in-homogeneities of the trap may distort long-range corre-lations in the system, the dynamics on small scales or atlarge momenta, relevant for the direct cascade, remainsunaffected. We have confirmed this explicitly by meansof GPE simulations of a harmonically trapped system(not shown).We conclude this section with some remarks on di-mensionality. The observation of a turbulent cascade ina scalar BEC requires at least two spatial dimensions.This is owed to the fact that in a strictly D system, dueto kinematic constraints, there can be no ↔ scatter-ing processes redistributing momenta. If, however, thesystem is elongated, but not as strongly confined in thetransversal direction to be considered in the quasi- Dregime, turbulent energy transport sets in once transver-sal modes are excited. We have confirmed this by numer-ically simulating a D system in an elongated periodicbox.Compared to the D case, the study of reheating dy-namics in D is somewhat simpler for at least two rea-sons. First, as discussed in section II C, the effective in- teraction in D is independent of the cosmic scale factorand there is no need to continuously adjust the scatteringlength according to the specific expansion model chosen.Second, an absorption image of an atomic cloud, takenafter a TOF expansion or in situ , will always be a D pro-jection on the plane transversal to the optical axis. In D,this means that the distribution is integrated along theoptical axis, thereby mixing momenta of different mag-nitudes. We have mimicked this in a simulation of a Dsystem in a periodic box, where we checked in particularthat the scaling is robust if we integrate the momentumdistribution along one spatial dimension before perform-ing the radial average. Thus, our scheme can be readilyapplied to experimentally investigate reheating dynamicsin D, which is of fundamental interest from a cosmolog-ical point of view.
VI. CONCLUSION
In this work, we have demonstrated by means ofclassical-statistical simulations how single-component ul-tracold Bose gases with tunable self-interactions canbe employed for simulating the dynamics of post-inflationary reheating. We have shown how an expandingspacetime can be encoded in the time-dependence of theinteraction, and that in the special case of D, arbitraryexpanding backgrounds can be associated with a singleexperiment in post-processing. The resonant productionof particles in the preheating stage of the early universe,driven by the oscillating inflaton, is mimicked by a mod-ulation of the scattering length inducing parametric in-stabilities in the BEC.Following the proliferation of secondary instabilitiesdue to non-linear effects, the system enters a turbulentstate, characterized by a self-similar transport of energytowards higher momenta. We have analyzed the uni-versal scaling exponents describing this direct cascade,which in general undergoes a transition from a driven toa free regime after the inflaton has decayed. This tran-sition is characterized by a prescaling stage, where themomentum distribution has already attained its universalscaling form, but the scaling exponents are non-universaland time-dependent, yet fully describe the dynamics ofthe Bose gas. Our set-up opens the door to the experi-mental observation of both the universal direct cascadeand the prescaling regime.The approach to a non-thermal fixed point as well asprescaling represent a dramatic reduction of sensitivityto the initial state as well as a gradual loss of complexityin the system’s dynamics, leading to the emergence ofuniversal physics. Thus, from a more general point ofview, our set-up provides a distinct platform for exploringthe fundamental mechanisms of thermalization dynamicsin quantum many-body systems far-from-equilibrium.The final stage of the dynamics, when the system re-laxes to thermal equilibrium, is dominated by quantumfluctuations and therefore not captured by the classical-8statistical approximation. In contrast, it can be accessedby experimental studies in cold atomic systems, whichare quantum-mechanical by nature and not restricted tothe weak-coupling regime.
ACKNOWLEDGMENTS
The authors thank Uwe R. Fischer, Zoran Hadz-ibabic, Maurus Hans, Aleksas Mazeliauskas, Alek-sandr N. Mikheev, Julius Mildenberger, Christian-MarcelSchmied, Marius Sparn, Daniel Spitz, Helmut Strobel,Malo Tarpin, Gerard Valentí-Rojas, Celia Viermann, andTorsten V. Zache for discussions and collaborations onrelated topics. This work is funded by the DFG (Ger-man Research Foundation) – Project-ID – SFB
ISOQUANT, DFG BE / – , the ERC StartingGrant “StrEnQTh”, the ERC Advanced Grant “Entan-gleGen”, and the Province of Trento. The authors ac-knowledge support by the state of Baden-Württembergthrough bwHPC.A.C. and K.G. contributed equally to this work. Appendix A: Numerical methods1. Classical-statistical simulations
We simulate the dynamics of an ultracold Bose gas,governed by the Heisenberg equations of motion (13), bymeans of classical-statistical simulations [28, 29], knownin the literature also under the name truncated Wignersimulations [30, 31]. This method takes into accountquantum fluctuations in form of stochastic initial condi-tions, but relies on a deterministic time evolution gov-erned by semi-classical equations of motion. Quan-tum expectation values are obtained as statistical av-erages over multiple realizations. The following sum-mary of the truncated Wigner method is mainly basedon refs. [30, 31].For each realization, the initial field configuration issampled from the Wigner distribution of the correspond-ing quantum state, which is commonly taken as the Bo-goliubov state in equilibrium. Here, we consider a homo-geneous scalar BEC of N atoms at zero temperature ina box of volume V with periodic boundary conditions.The initial wave function is sampled as Ψ cl (0 , x ) = √ n e iθ + (cid:88) p (cid:54) =0 (cid:2) α p u p ( x ) − α ∗ p v ∗ p ( x ) (cid:3) . (A1)Here, the first term represents the condensate with par-ticle density n = N/V and phase θ . Due to the largeoccupancy of the condensate mode, the finite widthof its Wigner function can be neglected, and thus thesame value of the density can be used in each realiza-tion. To preserve the U (1) symmetry, the phase is ran-domly sampled from the uniform distribution over the interval [0 , π ) for each realization. The mode func-tions u p ( x ) = u p e i px / (cid:126) / √ V and v p ( x ) = v p e i px / (cid:126) / √ V are solutions of the Bogoliubov–de Gennes equations fora uniform system in a periodic box, with real coeffi-cients u p = (cid:112) ( (cid:15) p , + n g ) / (cid:15) p + 1 / and v p determinedby the normalization u p − v p = 1 (cf. section III D fordetails of the notation) [78]. In order to mimic quantumfluctuations, the quasi-particle amplitudes α p are sam-pled as complex Gaussian random variables with zeromean, satisfying α ∗ p α q = δ p , q / at zero temperature.Here, ( . . . ) denotes the ensemble average over all realiza-tions. This vacuum noise corresponds to an average occu-pancy of half a particle per mode, which is also known asthe “quantum half”. Unless stated otherwise, the vacuumnoise is cut off at the highest lattice momentum.Here, we use a simplified approach, setting u p = 1 and v p = 0 , such that the mode functions reduce to ordinaryplane waves. Effectively, this corresponds to preparingthe Bogoliubov ground state of an ideal gas, which isappropriate for our application since the precise natureof the quasi-particles seeding the parametric resonanceis unimportant. As a side effect, we observe a transientgrowth of population at low momenta in our simulationsat early times (cf. fig. 2). The resulting p − behaviorof the momentum distribution at low momenta matchesthe behavior of u p and v p for | p | → . Therefore, theobserved growth can be interpreted as an artifact of theeffective quench from an ideal Bose gas to a system withfinite interaction at time t = 0 .In the truncated Wigner method, quantum expecta-tion values are replaced by statistical averages over theWigner distribution. The latter correspond to expec-tation values of symmetrically ordered quantum opera-tors. Thus, the momentum distribution obtained fromclassical-statistical simulations corresponds to the onedefined in eq. (20) plus an extra contribution in formof the “quantum half” stemming from the commutationrelations, V | Ψ cl ( t, p ) | = f ( t, p ) + 12 . (A2)Here, Ψ cl ( t, p ) = (cid:82) d d x Ψ cl ( t, x )e − i px / (cid:126) denotes theFourier transform of the classical field Ψ cl ( t, x ) .The classical-statistical simulations presented in thiswork have typically been averaged over at least runs.The statistics is even further enhanced through radialaverages due to the isotropy of the system. The statisticalerror bars are thus typically smaller than the line widthin the plots.
2. Dimensionless Gross–Pitaevskii equation
Each realization is propagated in time according tothe GPE (21). To cast this equation into a dimen-sionless form suitable for numerical simulations, we ex-press time and length in units of the characteristic scales9 t = (cid:126) /n g and x = (cid:126) / √ mn g , respectively. In termsof the dimensionless time ˜ t = t/t , position ˜ x = x /x ,and field ˜Ψ cl = x d/ Ψ cl , the GPE takes the form i∂ ˜ t ˜Ψ cl = (cid:18) −
12 ˜ ∇ + ˜ g (cid:0) ˜ t (cid:1) | ˜Ψ cl | (cid:19) ˜Ψ cl , (A3)with the dimensionless coupling ˜ g = g × t / (cid:126) x d .In a quasi- D system (and in the absence of time-dependent modulations), the dimensionless coupling is,up to logarithmic corrections, given by ˜ g = √ πa s /a HO ,where a s is the s -wave scattering length and a HO is theoscillator length of the confining harmonic potential invertical direction [57]. If the wave function is normalizedto unity, the coupling changes as ˜ g → N ˜ g , where N isthe particle number. This shows that the parameters N and ˜ g enter the classical equations of motion (A3) onlythrough the product N ˜ g as the single relevant model pa-rameter. Moreover, also the dimensionless box length L/x = √ N ˜ g is fixed by this quantity as a consequenceof scale invariance. By contrast, on the quantum level,the relative magnitudes of the parameters N and ˜ g playa role (cf. appendix A 3).The GPE (A3) is discretized in space by means of aFourier pseudospectral discretization and propagated intime using the well-known split-step method [79]. The ac-curacy of the time evolution is further enhanced with thehelp of the optimized fourth-order time splitting schemegiven in table 2 of ref. [80].In our D simulations of turbulent reheating dynamics(cf. section IV), we use a spatial N g × N g grid with atleast N g = 512 grid points per dimension, while N g = 256 turns out to be sufficient for simulating early-time pre-heating dynamics (cf. section III). For a system of size L × L with L/x = 50 , the corresponding grid spacing ∆ x = L/N g ensures that the healing length ξ = x / √ ,which is the smallest physical length scale in the sys-tem, is well-resolved. Numerical stability for late-timedynamics is achieved by choosing the numerical time step ∆˜ t such that ∆˜ t ˜ k (cid:46) π , where ˜ k max = π/ ∆˜ x with ∆˜ x = ∆ x/x is the maximum numerical wave numbersupported by the grid [81].
3. Validity of classical-statistical simulations
Quantum dynamics is well described by classical-statistical simulations in the regime of large occupanciesand weak couplings [29, 48]. In particular, the truncatedWigner method is known to produce unphysical results ifthe number of virtual particles added to mimic quantumfluctuations becomes comparable to the number of realparticles in the system [30, 31]. At zero temperature, thefailure of the classical-statistical approximation is indi-cated by the instability of the vacuum and a resulting un-physical dependence of observables on the UV cutoff [48].This decay of the “quantum half”, which inevitably occurs − Radial momentum p/p ξ − − R e s c a l e dd i s t r i bu t i o n (cid:2) f ( p ) + (cid:3) × ˜ g / ˜ g r e f ˜ g/ ˜ g ref = 1˜ g/ ˜ g ref = 0 . g/ ˜ g ref = 0 . g/ ˜ g ref = 0 . g/ ˜ g ref = 0 . g/ ˜ g ref = 0 . g/ ˜ g ref = 0 . Figure 7. Snapshot of the radially averaged, rescaledmomentum distribution for the driven turbulent cascade(cf. fig. 4a) at time t = 640 × t for different valuesof the coupling ˜ g . The particle number is chosen ac-cording to N (˜ g ) = N ref × (˜ g/ ˜ g ref ) − , such that the product N ˜ g = N ref ˜ g ref remains constant. The vertical dashed linemarks the characteristic momentum (cf. section IV D), calcu-lated for the distribution corresponding to the smallest valueof ˜ g . The horizontal dotted lines mark the location of the“quantum half” after rescaling. Within the range of validityof classical-statistical simulations, the rescaled distributionsare expected to collapse to a single curve. Deviations canbe observed for ˜ g/ ˜ g ref (cid:38) . , indicating a breakdown of themethod for larger values of the coupling. at sufficiently late times, naturally restricts the classical-statistical method to the weak coupling regime, where theinstability is mitigated via a separation of scales. More-over, if the coupling is too strong, physical observablesrisk attaining a dependence on the UV cutoff already atrelatively early times through the spurious quantum pres-sure generated by the decaying vacuum. This is becausethe coupling controls the relative magnitude of mode oc-cupancies with respect to the vacuum noise, which, inturn, is regulated by the UV cutoff [48].To ensure that our numerical simulations of reheatingdynamics, in particular, in the turbulent scaling regimeat late times, lie within the range of validity of theclassical-statistical method, we follow a similar procedureas presented in ref. [48]. Formally, the classical statisticallimit is given by N → ∞ and g → , where N is the par-ticle number and g is the coupling, such that N g = const . To approach this limit in our D simulation, we fix thevalue N ˜ g = ( L/x ) by the choice L/x = 50 and repeatthe simulation for different values of the coupling ˜ g , usingthe parametrization N (˜ g ) = N ref × (˜ g/ ˜ g ref ) − with ref-erence parameters N ref = 1 × and ˜ g ref = 2 . × − .Figure 7 shows the radially averaged momentum distri-bution for the driven turbulence scenario (cf. fig. 4a) attime t = 640 × t for different values of the coupling ˜ g . Allconfigurations have been averaged over at least single0runs. To assess up to which value of the coupling thesimulations are valid, we rescale the curves by ˜ g/ ˜ g ref , inwhich case they are expected to lie on top of each other inthe region f ( t, p ) × ˜ g/ ˜ g ref (cid:38) . (Note that regions wherethe occupancy drops below unity are outside of the rangeof validity per se .) The vertical dashed line marks the po-sition of the characteristic momentum (cf. section IV D),calculated for the distribution with the lowest value of ˜ g ,around which the UV scaling properties of the distribu-tion are most sensitive. It can be seen that deviationsoccur for ˜ g/ ˜ g ref (cid:38) . . For larger values of the coupling,the shape of the power law changes and the distribu-tion hits the UV cutoff, such that cutoff-independenceof the results is no longer guaranteed. Furthermore, thedistribution at large momenta drops below the “quan-tum half”, as indicated by the horizontal dotted lines,which formally corresponds to the unphysical situationof negative occupancies. We have also checked, using asimilar approach as in ref. [48], that the contribution ofthe vacuum noise to the total energy budget of the sys-tem becomes significant for ˜ g/ ˜ g ref (cid:38) . , which coincideswith the values of the coupling for which deviations inthe momentum distribution are observed.While the early-time preheating dynamics presented insection III is only insignificantly affected by these defi-ciencies for the reference parameters N ref = 1 × and ˜ g ref = 2 . × − , this is no longer true for the late-timeturbulent reheating dynamics, as illustrated in fig. 7. Theresults presented in section IV have therefore been com-puted for ˜ g/ ˜ g ref = 0 . , which safely lies within the rangeof validity of classical-statistical simulations. Appendix B: Parametric resonance and Mathieu’sequation
To gain analytical insight in the mechanism behindparametric resonance, we rewrite the equation for thedensity fluctuations 24a in the form ∂ n p ( s ) ∂s + [ A p − q p cos(2 s )] n p ( s ) = 0 (B1)with the dimensionless time s = ωt/ π/ , and pa-rameters A p = (cid:15) p / ( (cid:126) ω/ and q p = r(cid:15) p , n g / ( (cid:126) ω/ ,where (cid:15) p = (cid:112) (cid:15) p , ( (cid:15) p , + 2 n g ) denotes the Bogoli-ubov dispersion relation. Equation (B1) is the stan-dard form of Mathieu’s equation [59]. Importantly, thisequation admits solutions which can be expressed asthe product of an oscillatory function and an exponen-tially growing prefactor ∼ e ζ p t , describing parametric res-onance [41, 56]. The condition for exact resonance isgiven by A p = 1 , or equivalently (cid:15) p = (cid:126) ω/ .The growth rate ζ p of unstable modes can be estimatedfor small modulation amplitudes q p (cid:28) using perturba-tion theory [60]. To leading order, it is given by ζ p = ω (cid:113) q p − ( A p − , (B2) which reduces to eq. (25) for the resonant momentum p res .Moreover, as it can be seen in eq. (B2), there is anentire range of modes around p res that experience a pos-itive growth rate, and which thus undergo parametricresonance. The width of this instability band is delim-ited by the modes satisfying A p = 1 ± q p , which increaseswith the amplitude r of the modulation.The growth rate predicted from eq. (B2) is comparedto our numerical simulations in fig. 3. Appendix C: Secondary instabilities
In this appendix, we discuss the secondary amplifica-tion of fluctuations due to the non-linear corrections ineq. (27). For modes that are stable on the linear level, theintegral on the right-hand side of eq. (27) acts as a sourceterm. It has a complicated momentum structure, but itstime evolution is straightforward since it is dominated bythe contribution from the exponentially growing unstablemodes q , for which | p − q | ≈ | q | ≈ p res . This conditioncan be satisfied if p (cid:46) p res , as can be seen from thetriangle inequality, p = | p − q + q | ≤ | p − q | + | q | ≈ p res . For p (cid:29) p res , on the other hand, the growth of theintegral is exponentially suppressed. The narrower theprimary resonance, the more rapid the integral drops at p ≈ p res .As in the previous appendix, we assume here that theresonance is not too strong, such that ζ p res (cid:28) ω . Thetime evolution of resonant modes can then be written asa product of an exponential growth and an oscillatingfunction [59, 60], n p ∝ e ζ p res t e ± iω res t , θ p ∝ e ζ p res t e ± iω res t , with the growth rate ζ p res given by eq. (B2) and ω res being the frequency of the resonant modes. As a result,eq. (27) takes the form of a forced harmonic oscillator, ¨ n p + ω p n p = e ζ p res t (cid:0) S p e iω res t + Q p e − iω res t + P p (cid:1) . (C1)Here, Q p and S p contain the momentum-dependent con-tributions of the integrals in eq. (27) corresponding topositive and negative frequencies, respectively, and P p contains the non-oscillating part. The contribution tothe amplitude of n p from the first term on the right-hand side of eq. (C1) proportional to S p is given by | n p | = e ζ p res t |S p | (cid:113) (8 ζ p res ω res ) + ( ω p + (2 ζ p res ) − ω ) , (C2)and analogous expressions can be obtained for Q p and P p , respectively. In other words, the modes in the above-mentioned momentum range are amplified as exp(2 ζ p res t ) ω p ≈ ω res . Asmoother spectrum is expected in the case of strong driv-ing, i.e., large ζ p res .The pattern described above is indeed observed infig. 2, where the secondary and higher-order peaks atmultiple frequencies are clearly visible. Appendix D: Analysis of instabilities in terms ofquantum equations of motion for correlationfunctions
In this appendix, we study the instabilities arising inthe presence of a periodically modulated coupling usingquantum equations of motion for the two-point corre-lation functions. We validate the results of section III,where quantum fluctuations have been treated within theclassical-statistical approximation.The classical action for a (non-relativistic) bosonicfield Ψ with a quartic self-interaction is given by [82] S [Ψ] = (cid:90) x,C (cid:104) i (cid:126) ∗ ˙Ψ − ˙Ψ ∗ Ψ) − (cid:126) m ∇ Ψ ∗ ∇ Ψ − g ∗ Ψ) (cid:105) . (D1)The integration is performed over the Schwinger–Keldyshclosed-time contour [83]. In the following discussion, wetreat ϕ = √ and ϕ = √ as the two inde-pendent field components.In quantum theory, the central objects are field cor-relation functions. The one-point correlation function isdefined as φ i ( x ) = (cid:104) ˆ ϕ i ( x ) (cid:105) . Throughout this appendix,we use the notation x = ( x , x ) ≡ ( t, x ) . For homoge-neous systems, the one-point function depends only ontime and we can identify it with the wave function ofthe condensate. The time-ordered connected two-pointfunction is defined as [82, 84] G ij ( x, y ) = (cid:104)T ˆ ϕ i ( x ) ˆ ϕ j ( y ) (cid:105) − φ i ( x ) φ j ( y )= F ij ( x, y ) − i C ( x − y ) ρ ij ( x, y ) , (D2)where T is the time-ordering operator along the closed-time contour and sgn C ( x − y ) takes the values ± depending on whether x is after or before y on theSchwinger–Keldysh contour. In the last step, the propa-gator has been decomposed into the spectral function ρ and the statistical propagator F , defined as F ij ( x, y ) = 12 (cid:104){ ˆ ϕ i ( x ) , ˆ ϕ j ( y ) }(cid:105) − φ i ( x ) φ j ( y ) , (D3) ρ ij ( x, y ) = i (cid:104) [ ˆ ϕ i ( x ) , ˆ ϕ j ( y )] (cid:105) . (D4)Here, [ · , · ] and {· , ·} denote the commutator and theanti-commutator, respectively. The spectral function de-scribes the spectrum of the excitations and, in particular,encodes the equal-time commutation relations, while thestatistical propagator characterizes its occupations.The equations of motion for the one- and two-pointfunctions can be obtained by a variational principle from the two-particle irreducible (2PI) effective action. Forthe two-point functions, they have the form [84] D ik ( x ) F kj ( x, y ) = (cid:90) x t dz Σ ρik ( x, z ) F kj ( z, y ) − (cid:90) y t dz Σ Fik ( x, z ) ρ kj ( z, y ) , D ik ( x ) ρ kj ( x, y ) = (cid:90) x y t dz Σ ρik ( x, z ) ρ kj ( z, y ) , (D5)where D ik ( x ) = [ − i (cid:126) σ ,ik ∂ x − M ik ( x )] , M ij ( x ) = δ ij (cid:104) − (cid:126) ∆ x m + g (cid:16) φ k ( x ) φ k ( x ) + F kk ( x, x ) (cid:17)(cid:105) + g (cid:16) φ i ( x ) φ j ( x ) + F ij ( x, x ) (cid:17) , (D6)and Σ F and Σ ρ are the components of proper self-energy Σ ab ( x, y ) = Σ Fab ( x, y ) − ( i/
2) sgn C ( x − y )Σ ρab ( x, y ) ,which sums all non-local two-point one-particle-irreducible diagrams constructed from the interactionvertices in the action (D1) and with lines associatedto full time-ordered propagators G . Summation overrepeated indices is implied in all expressions in thisappendix. The presence of a non-vanishing one-pointfunction generates effective cubic interactions in addi-tion to the quartic one, which follows from the shift ϕ i → φ i + ϕ i . The resulting interaction part of theaction, in terms of the fields ϕ i , reads S int [ ϕ, φ ] = − g (cid:90) x,C (cid:104) ϕ i ϕ i ϕ j ϕ j + 4 φ i ϕ i ϕ j ϕ j (cid:105) . (D7)Equations (D5) are supplemented with analogousequations of motion for the one-point function φ . For thediscussion of instabilities, the back-reaction of producedexcitations on φ is negligible, and, as a result, its evolu-tion is approximately governed by the classical GPE.
1. Primary instabilities
In the linear regime, the self-energies in eq. (D5) as wellthe F -dependent corrections to the mass in (D6) can beneglected. We consider the evolution of the statisticalpropagator F , which can be written as (cid:104) − i (cid:126) σ ,ik ∂ x + δ ik (cid:16) (cid:126) ∆ x m − g φ l ( x ) φ l ( x ) (cid:17) − gφ i ( x ) φ k ( x ) (cid:105) F kj ( x, y ) = 0 . (D8)2 Figure 8. Diagrammatic illustration of the local (left) andnon-local (right) one-loop contributions to the self-energy.The second diagram accounts for secondary instabilities.
It is convenient to introduce the “normal” and the“anomalous” correlation functions, defined as [82] (cid:101) n ( x, y ) = 12 (cid:104){ ˆΨ( x ) , ˆΨ † ( y ) }(cid:105) − (cid:104) ˆΨ( x ) (cid:105)(cid:104) ˆΨ † ( y ) (cid:105) = 12 F ij ( x, y )( δ ij + σ ,ij ) = (cid:101) n ∗ ( y, x ) , (D9) (cid:101) m ( x, y ) = 12 (cid:104){ ˆΨ( x ) , ˆΨ( y ) }(cid:105) − (cid:104) ˆΨ( x ) (cid:105)(cid:104) ˆΨ( y ) (cid:105) = 12 F ij ( x, y )( σ ,ij + iσ ,ij ) = (cid:101) m ( y, x ) , (D10)which are complex-valued functions. Assuming homo-geneity, these functions depend only on the relative spa-tial coordinate, i.e., (cid:101) n = (cid:101) n ( t, t (cid:48) , x − y ) . The evolutionequations (D8) for the components of the statistical prop-agator, Fourier transformed with respect to this relativecoordinate, are given by (cid:104) i (cid:126) ∂ t − p m − g ( t ) | Ψ( t ) | (cid:105)(cid:101) n p ( t, t (cid:48) ) − g ( t )Ψ ( t ) (cid:101) m ∗ p ( t, t (cid:48) ) = 0 , (cid:104) i (cid:126) ∂ t − p m − g ( t ) | Ψ( t ) | (cid:105) (cid:101) m p ( t, t (cid:48) ) − g ( t )Ψ ( t ) (cid:101) n ∗ p ( t, t (cid:48) ) = 0 . These equations have the same form as the ones for thelinearized (classical-statistical) fluctuation field δ Ψ p ( t ) .In particular, the linear combinations ψ ∗ ( t ) (cid:101) n p ( t, t (cid:48) ) + ψ ( t ) (cid:101) m ∗ p ( t, t (cid:48) ) and ( i/ (cid:101) m ∗ p ( t, t (cid:48) ) /ψ ( t ) − (cid:101) n p ( t, t (cid:48) ) /ψ ∗ ( t )) satisfy the wave equations (24). Therefore, ignoring thetime-dependence of the coupling, one has (cid:101) n p ( t, t (cid:48) ) ∝ e i ( ± ω p − g n )( t − t (cid:48) ) , (cid:101) m p ( t, t (cid:48) ) ∝ e i ( ± ω p − g n )( t + t (cid:48) ) , where ω p corresponds to the Bogoliubov dispersion re-lation. The dependence on t (cid:48) in the above expression isdetermined from the symmetry properties of the func-tions (cid:101) n p ( t, t (cid:48) ) and (cid:101) m p ( t, t (cid:48) ) , given in (D9).The unstable modes exhibit an exponential amplifica-tion, (cid:101) n p ( t, t (cid:48) ) ∝ (cid:101) m p ( t, t (cid:48) ) ∝ exp( ζ p ( t + t (cid:48) )) , as explainedin section (III D). Here we again employed the symmetryproperties of (cid:101) n and (cid:101) m . The same holds for the compo-nents of the spectral function. However, at t = t (cid:48) thespectral function is determined solely by the commuta-tion relations, as follows from its definition (D4), ρ ij ( x, y ) | x = y = − i σ ,ij δ ( x − y ) . (D11)In the next subsection we include the first non-linear cor-rections to describe the secondary amplification of fluc-tuations.
2. Secondary instabilities
As the fluctuations become sufficiently strong, the lin-ear approximation for the dynamics becomes inapplica-ble. The non-local one-loop correction, which is diagram-matically illustrated in the right part of fig. 8, accountsfor secondary instabilities [19, 41]. The correspondingself-energy, which can be calculated by means of stan-dard Feynman diagram techniques, reads Σ ij ( x, y ; φ, G ) = − g ( x ) g ( y )2 (cid:126) (cid:104) φ i ( x ) φ j ( y ) G ac ( x, y ) G ac ( x, y ) + 2 φ b ( x ) φ d ( y ) G ij ( x, y ) G bd ( x, y )+2 φ b ( x ) φ d ( y ) G id ( x, y ) G bj ( x, y ) + 2 φ b ( x ) φ j ( y ) G ic ( x, y ) G bc ( x, y ) + 2 φ i ( x ) φ d ( y ) G aj ( x, y ) G ad ( x, y ) (cid:105) . (D12)The other one-loop correction is the tadpole diagram, shown in the left part of fig. 8, which is explicitly included inthe definition (D6) of the effective mass . As can be seen from that definition, this correction becomes relevant onlylater, when the components of F ij become of the order of the condensate density | Ψ | (see also [41]).Including the one-loop self-energy in eq. (D5) leads to integrals over the past evolution, i.e., memory integrals, onthe right-hand side. Due to the exponential growth of fluctuations, these integrals are, however, dominated by thelatest times [19, 84], which allows to simplify them. Note that Σ contains terms proportional to G αβ ( x, y ) G γδ ( x, y ) = F αβ ( x, y ) F γδ ( x, y ) − ρ αβ ( x, y ) ρ γδ ( x, y ) (cid:124) (cid:123)(cid:122) (cid:125) Σ F − i C ( x − y ) F αβ ( x, y ) ρ γδ ( x, y ) + ρ αβ ( x, y ) F γδ ( x, y ) (cid:124) (cid:123)(cid:122) (cid:125) Σ ρ . The (cid:82) Σ ρ F term in the limit of z → x represents a “tadpole”-like contribution, which becomes important only at3later times [41]. For the (cid:82) Σ F ρ term one obtains − (cid:90) y t dz Σ Fac ( x, z ) ρ cb ( z, y ) ≈ − (cid:90) y y − c dz Σ Fac ( x, z ) ρ cb ( z, y ) ≈ i (cid:126) c Σ Fac ( x, y ) σ ,cb , where c − = O ( g n ) is some parameter of the order of the oscillation frequency of the condensate. In Σ F , there willbe terms ∝ F and ∝ ρ . We first consider the ∝ F part in Σ F , which is equivalent to Σ | G → F . Inserting eq. (D12)into the equations of motion (D5) and using the transformation (D9) and (D10), one arrives at the following equationsfor the momentum modes of the normal and the anomalous two-point functions, (cid:104) i (cid:126) ∂ t − p m − g ( t ) | Ψ( t ) | (cid:105)(cid:101) n p − g ( t )Ψ ( t ) (cid:101) m ∗ p = icg ( t ) g ( t (cid:48) ) (cid:90) q (cid:104) t )Ψ ∗ ( t (cid:48) ) (cid:16)(cid:101) n q (cid:101) n ∗ p − q + (cid:101) m q (cid:101) m ∗ p − q (cid:17) + 4Ψ ∗ ( t )Ψ ∗ ( t (cid:48) ) (cid:101) n q (cid:101) m p − q + 4Ψ( t )Ψ( t (cid:48) ) (cid:101) n q (cid:101) m ∗ p − q + 2Ψ ∗ ( t )Ψ( t (cid:48) ) (cid:101) n q (cid:101) n p − q (cid:105) , (D13) (cid:104) i (cid:126) ∂ t − p m − g ( t ) | Ψ( t ) | (cid:105) (cid:101) m p − g ( t )Ψ ( t ) (cid:101) n ∗ p = − icg ( t ) g ( t (cid:48) ) (cid:90) q (cid:104) t )Ψ( t (cid:48) ) (cid:16) n q n ∗ p − q + m q m ∗ p − q (cid:17) + 4Ψ ∗ ( t )Ψ( t (cid:48) ) (cid:101) m q (cid:101) n p − q + 4Ψ( t )Ψ ∗ ( t (cid:48) ) (cid:101) m q (cid:101) n ∗ p − q + 2Ψ ∗ ( t )Ψ ∗ ( t (cid:48) ) (cid:101) m q (cid:101) m p − q (cid:105) , (D14)where we have omitted the arguments t and t (cid:48) from thetwo-point functions.The ∝ ρ part in Σ F , which is equivalent to Σ | G → iρ/ ,leads to an analogous expression on the right-hand side.However, the components ρ are of order one in the t (cid:48) → t limit (in momentum space) due to eq. (D11). As aresult, that contribution is negligible compared to theexponentially growing one from the ∝ F part. Thisalso justifies the use of classical-statistical simulations forstudying the evolution of equal-time observables.The momentum integrals on the right-hand side areanalogous to the ones in eq. (27) and therefore lead toa secondary exponential amplification, i.e., (cid:101) n p ( t, t (cid:48) ) ∝ exp(2 ζ p ( t + t (cid:48) )) , in the momentum range p (cid:46) p res . Onecan check that the terms on the right-hand side os-cillate as functions of t either as exp( − ig n t ) or as exp( − i ( g n ± ω res ) t ) . Hence, one expects enhancedamplification near ω p = 0 and ω p = 2 ω res . The resonantfrequencies correspond to the conservation of energy (inaddition to the conservation of momentum) in the ver-tices of the Feynman diagrams.To conclude, in this appendix we have analyzed theprimary and secondary instabilities within the non-equilibrium 2PI framework and explained the validity ofthe classical-statistical approximation, justifying its usein section III to describe the regime when typical occu-pation numbers of excitations become sufficiently large. Here, we neglect the time-dependence of g , which is suppressedfor small values of r . Appendix E: Maximum likelihood technique forextracting scaling exponents
In this appendix, we describe the maximum likelihoodtechnique used in section IV to extract the scaling ex-ponents α and β describing self-similar time evolutionaccording to eq. (29) from our numerical data.To extract the exponents, we use numerical fits, sim-ilar to the procedure in ref. [85]. More specifically, atdifferent times t , the distribution functions at t and t/ are compared. Assuming the behavior of eq. (29), wetest different values of the exponents. For each pair ofexponents, we define the (normalized) likelihood W ( α, β ) ∝ exp (cid:104) − χ ( α, β ) / (2 ¯ χ ) (cid:105) , (E1)where χ ( α, β ) = (cid:82) [( f ( t (cid:48) , p ) − α f ( t (cid:48) / , β p )) p l ] d p (cid:82) [ f ( t (cid:48) , p ) p l ] d p . (E2)Here ¯ χ is the minimal observed value of χ ( α, β ) , whichcorresponds to the best fit. We use the third momentof the radial distribution function because for d = 2 and ω p = p / m , it is proportional to the energy density permomentum interval. The integrals are performed overlarge momenta for which f ( p ) (cid:38) . We find the best-fitting values of the exponents and identify the fittingerrors with the Gaussian widths of the marginal likeli-hoods (cid:82) W dβ and (cid:82)
W dα .The variable t always denotes the time elapsed from thebeginning of the simulation. This includes the stage ofparametric resonance, which is unrelated to turbulence.While irrelevant for late times, a more reasonable choiceof the initial time, closer to the beginning of turbulentdynamics, improves the scaling analysis at early times,4by allowing to avoid the extraction of large exponentsdecreasing with time. For the considered set of parame-ters in section IV, we set t = 63 × t in the fitting routinedescribed above. Appendix F: Analysis of prescaling using themethod of moments
To gain a deeper understanding of the phenomenonof prescaling, we provide here a complementary scalinganalysis using the method of moments [46]. As explainedin appendix E, the maximum likelihood technique usedto extract the scaling exponents in figs. 4 and 5 locallycompares the distributions at two reference times t and t > t . By iterating over all times t , time-dependentscaling exponents are obtained that best collapse thepairs of distributions f ( t , p ) and f ( t , p ) on top of eachother. By contrast, the method introduced in ref. [46]relies on the moments as global properties of the distri-bution and allows one to extract instantaneous scalingexponents α ( t ) and β ( t ) that do not depend on a refer-ence time. In what follows, we briefly outline the methodof moments and apply it to the prescaling scenario dis-cussed in section IV B.The n -th moment of the distribution f ( t, p ) is definedas M ( n ) ( t ) = V (cid:90) d d p (2 π (cid:126) ) d (cid:18) pp (cid:19) n f ( t, p ) , (F1)where V is the volume, p = | p | , and p is an arbitrary mo-mentum scale to make the moment dimensionless. Notethat in an isotropic system, as considered here, the dis-tribution in fact only depends on the magnitude of themomentum, f ( t, p ) = f ( t, p ) . For each moment M ( n ) ,the integrand is peaked around a certain characteristicmomentum whose scaling properties are probed. In par-ticular, the zeroth moment is the total particle numberand the second moment is proportional to the total ki-netic energy for a system with quadratic dispersion.The most straightforward way of introducing time-dependent scaling exponents α ( t ) and β ( t ) is to makethe replacements α → α ( t ) and β → β ( t ) in the scal-ing ansatz eq. (29), such that the latter reads f ( t, p ) = s α ( t ) f S (cid:0) s β ( t ) p (cid:1) with s = t/t ref . However, the exponentsdefined in this way depend on the reference time t ref . Tolift this dependency, it is advantageous to define α ( t ) and β ( t ) instead in terms of the scaling functions s α ( t ) = exp (cid:20)(cid:90) tt ref α ( t (cid:48) ) d t (cid:48) t (cid:48) (cid:21) (F2)and s β ( t ) defined in an analogous way. The scaling ansatzin eq. (29) then generalizes to f ( t, p ) = s α ( t ) f S [ s β ( t ) p ] . (F3)For constant exponents α and β , the power-law scalingof eq. (29) is recovered. Time t/t − − S c a li n g f un c t i o n s α ( t ) ∝ t − ∝ t − n = 1 , n = 2 n = 1 , n = 3 n = 1 , n = 4 n = 2 , n = 3 n = 2 , n = 4 n = 3 , n = 4 Time t/t − S c a li n g f un c t i o n s β ( t ) ∝ t − / ∝ t − /
360 380 4000 . . Figure 9. Scaling functions s α ( t ) (upper panel) and s β ( t ) (lower panel) with respect to the reference time t ref = 85 × t extracted from the moments of orders ≤ n < n ≤ . Thevertical dotted line marks the time t ≈ . × t whenthe modulation is switched off instantaneously. Before thispoint, the oscillatory behavior of the moments is directly re-flected in the evolution of the scaling functions. As driventurbulence develops, their time averages approach power lawswith exponents close to the predictions from kinetic theory, α driven = − and β driven = − / (dashed lines). After themodulation is switched off, the scaling functions exhibit akink and the oscillations vanish. In the subsequent evolu-tion, the scaling functions extracted from different momentsevolve asynchronously, until they adopt a power-law behav-ior again for times t (cid:38) × t with exponents close to thepredictions from kinetic theory in the regime of free turbu-lence, α free = − and β free = − / (dashed-dotted lines).The system evolves self-similarly where all curves have thesame slope, as analyzed in fig. 10. Inserting the above scaling ansatz into eq. (F1), it isstraightforward to derive that the moments scale withtime as M ( n ) ( t ) = s α ( t ) s d + nβ ( t ) M ( n ) ( t ref ) . (F4)Given a pair of moments M ( n ) ( t ) and M ( n ) ( t ) with n (cid:54) = n , it is thus possible to express the scaling func-5tions s α ( t ) and s β ( t ) in terms of these moments, s α ( t ) = (cid:34) M d + n ( n ) ( t ) /M d + n ( n ) ( t ref ) M d + n ( n ) ( t ) /M d + n ( n ) ( t ref ) (cid:35) / ( n − n ) , (F5a) s β ( t ) = (cid:20) M ( n ) ( t ) /M ( n ) ( t ref ) M ( n ) ( t ) /M ( n ) ( t ref ) (cid:21) / ( n − n ) , (F5b)which allows one to access the scaling exponents as α ( t ) = 1 n − n dd log t log M d + n ( n ) ( t ) M d + n ( n ) ( t ) , (F6a) β ( t ) = 1 n − n dd log t log M ( n ) ( t ) M ( n ) ( t ) . (F6b)Note that the dependency on the reference time t ref dropsout as a consequence of the derivatives.Since each moment condenses the information aboutthe distribution into a single number, it is necessary toexamine many different moments, each sensitive to a dif-ferent characteristic momentum, to certify that a givendistribution scales as a whole. Figures 9 and 10 showthe scaling functions and time-dependent scaling expo-nents, respectively, extracted using the method of mo-ments outlined above for all possible combinations of mo-ments ≤ n < n ≤ . We exclude the zeroth momentfrom our analysis since the associated particle numberconservation is incompatible with the energy-conservingdirect cascade. All simulation parameters are the sameas in fig. 5, with the exception that the modulation isswitched off at a later time, t = 56 × π/ω ≈ . × t ,allowing one to better distinguish the regimes of drivenand free turbulence. The classical-statistical simulationpresented in figs. 9 and 10 has been conducted on a spa-tial grid of × grid points to increase the UV reso-lution and the data has been averaged over runs. Toreduce numerical errors of the calculated moments due toinstabilities near the UV cutoff, we restrict the integral ineq. (F1) to absolute momenta | p | ≤ × p ξ ( | p | ≤ π (cid:126) / ∆ x ,where ∆ x is the lattice spacing) before (after) switchingoff the modulation. As explained in appendix E, we shiftthe origin of time as t → t − × t to improve the scalinganalysis at early times.Figure 9 shows the scaling functions s α ( t ) and s β ( t ) with respect to the reference time t ref = 85 × t . As tur-bulence develops, the scaling functions approach powerlaws s α ( t ) ∝ ( t/t ref ) α and s β ( t ) ∝ ( t/t ref ) β . In the drivenregime, before switching off the modulation, the corre-sponding exponents are close to the analytical predictionsfrom kinetic theory, α driven = − and β driven = − / (see section IV C). At the time when the modulation isswitched off instantaneously, the scaling functions ex-hibit a kink and continue evolving asynchronously forsome time. This indicates that the shape of the distribu-tion is not preserved exactly, but slightly readjusts dur-ing the transition from driven to free turbulence. Thechanges in the shape of the distribution are more pro-nounced at lower momenta, which is reflected by the fact
500 1000 1500 2000 2500
Time t/t − . − . − . − . − . − . S c a li n g e x p o n e n t ¯ β ( t )
400 420 440 − β ( t )
500 1000 1500 2000 2500
Time t/t R a t i oo f e x p o n e n t s ¯ α ( t ) / ¯ β ( t ) n = 1 , n = 2 n = 1 , n = 3 n = 1 , n = 4 n = 2 , n = 3 n = 2 , n = 4 n = 3 , n = 4 Figure 10. Time-dependent scaling exponents extracted fromthe moments of orders ≤ n < n ≤ . The inset in theupper panel shows the instantaneous exponent β ( t ) , whichstrongly oscillates due to the modulation. The time-averagedquantities ¯ β ( t ) and ¯ α ( t ) / ¯ β ( t ) are displayed in the main plotsin the upper and lower panel, respectively. The data hasbeen smoothened using simple moving means and the shadedregions show the corresponding moving standard deviations.The vertical dotted line represents the time when the mod-ulation is switched off. Before this point, in the regime ofdriven turbulence, the exponents are approximately constantand close to the analytical predictions from kinetic theory, β driven = − / and α driven /β driven = 2 (dashed lines). Afterthe modulation is stopped, the exponents jump discontinu-ously and the exponents extracted from different combina-tions of moments exhibit discrepancies. For t (cid:38) × t ,they converge and continue evolving as a single curve, certi-fying the existence of a prescaling regime of self-similar timeevolution. The exponent β ( t ) gradually approaches the uni-versal value β free = − / (dashed-dotted line in the upperpanel) predicted from kinetic theory in the regime of free tur-bulence, while the ratio of the exponents quickly adjusts tothe prediction α free /β free = 4 (dashed-dotted line in the lowerpanel), reflecting energy conservation. that those curves in fig. 9 extracted from lower momentsreact stronger at the transition. Despite these local ad-justments, it is remarkable that the global scaling formof the distribution is approximately preserved during thetransition from driven to free turbulence, as can be seenin fig. 5. For t (cid:38) × t , the scaling functions again6approach a power law with exponents close to the ana-lytical predictions for free turbulence, α free = − and β free = − / . In the region where the scaling func-tions extracted from different combinations of momentsall evolve in parallel, the distribution scales self-similarlyas a whole. To certify that the distribution also ex-hibits prescaling, i.e., a self-similar evolution with time-dependent scaling exponents, we analyze the change ofthe scaling functions as quantified by the exponents α ( t ) and β ( t ) in what follows.In the driven regime, the scaling functions oscillate intime, as indicated in the inset in the lower panel of fig. 9.These oscillations originate from the modulated inter-action, which causes the distribution and therefore themoments to oscillate in time (cf. fig. 3). As a result, theinstantaneous scaling exponents α ( t ) and β ( t ) , which arederivatives of the scaling functions according to eq. (F6),exhibit the strong oscillatory behavior shown for β ( t ) inthe inset in the upper panel of fig. 10. To connect withour previous results obtained using the maximum like-lihood technique (cf fig. 5), it is convenient to considerinstead the time-averaged exponents ¯ α ( t ) and ¯ β ( t ) , de-fined as ¯ β ( t ) = 1log(1 + T /t ) (cid:90) t + Tt β ( t (cid:48) ) d t (cid:48) t (cid:48) , (F7)and analogously for ¯ α ( t ) , where T = 2 π/ω is the modu-lation period. After switching off the modulation, the os-cillations of the moments vanish and the time-averagingof the exponents is not required any more.Figure 10 shows the quantities ¯ β ( t ) (upper panel) and ¯ α ( t ) / ¯ β ( t ) (lower panel) as a function of time. The deriva-tives in eq. (F6) are sensitive to small fluctuations in thedata due to statistical uncertainties, resulting in a largespread of the data points. This is especially true wheremoments of lower orders are involved since they probe the distribution at smaller momenta with a lower density ofstates. To better visualize the trend of the exponents, thedata has been smoothened by calculating simple movingmeans (separately in the regimes of driven and free tur-bulence) using a window size of ( ) periods involving ( ) data points for n > ( n = 1 ). The shadedregions represent the corresponding moving standard de-viations.Before the modulation is switched off, the system is inthe state of driven turbulence. Indeed, for t (cid:38) × t ,both ¯ β ( t ) and ¯ α ( t ) / ¯ β ( t ) are approximately constant andfor all considered moments close to the analytical pre-dictions β driven = − / and α driven /β driven = 2 (cf. sec-tion IV C). After suddenly switching off the modulation,the exponents jump discontinuously, reflecting the kinkin the scaling functions (cf. fig. 9). Furthermore, thereare discrepancies in the values of the exponents extractedfrom different combinations of moments after the mod-ulation is stopped. As discussed above, this indicatesa slight readjustment of the shape of the distribution,which is consistent with the behavior in fig. 5, where theexponents exhibit large error bars in this regime. Fortimes t (cid:38) × t , the exponents extracted from differ-ent moments converge to a single curve, certifying a self-similar evolution of the distribution as a whole. Duringthis prescaling stage, the exponent ¯ β ( t ) evolves gradu-ally towards the analytical prediction for free turbulence, β free = − / . By contrast, the ratio ¯ α ( t ) / ¯ β ( t ) adjustsrather quickly to the value α free /β free = 4 , expressingenergy conservation.In conclusion, our analysis certifies the existence of aprescaling regime during the transition from driven tofree turbulence, where the momentum distribution scalesself-similarly with time-dependent scaling exponents ap-proaching gradually their universal values. [1] A. H. Guth, Inflationary universe: A possible solution tothe horizon and flatness problems, Phys. Rev. D , 347(1981).[2] A. Starobinsky, A new type of isotropic cosmologicalmodels without singularity, Phys. Lett. B , 99 (1980).[3] A. A. Starobinsky, Dynamics of phase transition in thenew inflationary universe scenario and generation of per-turbations, Phys. Lett. , 175 (1982).[4] L. Kofman, A. D. Linde, and A. A. Starobinsky, Re-heating after inflation, Phys. Rev. Lett. , 3195 (1994),arXiv:hep-th/9405187 [hep-th].[5] L. Kofman, A. D. Linde, and A. A. Starobinsky, Towardsthe theory of reheating after inflation, Phys. Rev. D56 ,3258 (1997), arXiv:hep-ph/9704452 [hep-ph].[6] M. A. Amin, M. P. Hertzberg, D. I. Kaiser, andJ. Karouby, Nonperturbative dynamics of reheating af-ter inflation: A review, Int. J. Mod. Phys.
D24 , 1530003(2014), arXiv:1410.3808 [hep-ph].[7] R. Micha and I. I. Tkachev, Turbulent thermalization, Phys. Rev.
D70 , 043538 (2004), arXiv:hep-ph/0403101[hep-ph].[8] J. Berges, K. Boguslavski, S. Schlichting, and R. Venu-gopalan, Universality far from equilibrium: From super-fluid Bose gases to heavy-ion collisions, Phys. Rev. Lett. , 061601 (2015), arXiv:1408.1670 [hep-ph].[9] B. Nowak and T. Gasenzer, Universal dynamics on theway to thermalization, New J. Phys. , 093052 (2014),arXiv:1206.3181 [cond-mat.quant-gas].[10] A. P. Orioli, K. Boguslavski, and J. Berges, Universal self-similar dynamics of relativistic and nonrelativistic fieldtheories near nonthermal fixed points, Phys. Rev. D92 ,025041 (2015), arXiv:1503.02498 [hep-ph].[11] M. Prüfer, P. Kunkel, H. Strobel, S. Lannig, D. Lin-nemann, C.-M. Schmied, J. Berges, T. Gasenzer, andM. K. Oberthaler, Observation of universal dynamics ina spinor Bose gas far from equilibrium, Nature (London) , 217 (2018), arXiv:1805.11881 [cond-mat.quant-gas].[12] S. Erne, R. Bücker, T. Gasenzer, J. Berges, and J. Schmiedmayer, Universal dynamics in an isolatedone-dimensional Bose gas far from equilibrium, Na-ture (London) , 225 (2018), arXiv:1805.12310 [cond-mat.quant-gas].[13] P. O. Fedichev and U. R. Fischer, ’cosmological’ quasi-particle production in harmonically trapped superfluidgases, Phys. Rev.
A69 , 033602 (2004), arXiv:cond-mat/0303063 [cond-mat].[14] U. R. Fischer and R. Schutzhold, Quantum simulationof cosmic inflation in two-component Bose–Einstein con-densates, Phys. Rev.
A70 , 063615 (2004), arXiv:cond-mat/0406470 [cond-mat].[15] M. Uhlmann, Y. Xu, and R. Schutzhold, Aspects of cos-mic inflation in expanding Bose–Einstein condensates,New J. Phys. , 248 (2005), arXiv:quant-ph/0509063.[16] C. Neuenhahn and F. Marquardt, Quantum simulation ofexpanding space–time with tunnel-coupled condensates,New J. Phys. , 125007 (2015), arXiv:1208.2255 [cond-mat.quant-gas].[17] C. Neuenhahn, A. Polkovnikov, and F. Marquardt, Lo-calized phase structures growing out of quantum fluctua-tions in a quench of tunnel-coupled atomic condensates,Phys. Rev. Lett. , 085304 (2012), arXiv:1112.5982[cond-mat.quant-gas].[18] A. Posazhennikova, M. Trujillo-Martinez, and J. Kroha,Inflationary quasiparticle creation and thermalizationdynamics in coupled Bose–Einstein condensates, Phys.Rev. Lett. , 225304 (2016), arXiv:1603.04898 [cond-mat.quant-gas].[19] T. V. Zache, V. Kasper, and J. Berges, Inflation-ary preheating dynamics with two-species condensates,Phys. Rev. A95 , 063629 (2017), arXiv:1704.02271 [cond-mat.quant-gas].[20] S. Robertson, F. Michel, and R. Parentani, Nonlin-earities induced by parametric resonance in effectively1D atomic Bose condensates, Phys. Rev.
D98 , 056003(2018), arXiv:1802.00739 [cond-mat.quant-gas].[21] M. Wittemer, F. Hakelberg, P. Kiefer, J.-P. Schröder,C. Fey, R. Schützhold, U. Warring, and T. Schaetz,Phonon pair creation by inflating quantum fluctuationsin an ion trap, Phys. Rev. Lett. , 180502 (2019),arXiv:1903.05523 [quant-ph].[22] C. Chin, R. Grimm, P. Julienne, and E. Tiesinga, Fesh-bach resonances in ultracold gases, Rev. Mod. Phys. ,1225 (2010), arXiv:0812.1496 [cond-mat.other].[23] A. D. Linde, Chaotic inflation, Phys. Lett. , 177(1983).[24] S. Eckel, A. Kumar, T. Jacobson, I. B. Spielman, andG. K. Campbell, A rapidly expanding Bose-Einstein con-densate: an expanding universe in the lab, Phys. Rev. X8 , 021021 (2018), arXiv:1710.05800 [cond-mat.quant-gas].[25] V. Gritsev, P. Barmettler, and E. Demler, Scalingapproach to quantum non-equilibrium dynamics ofmany-body systems, New J. Phys. , 113005 (2010),arXiv:0912.2744 [cond-mat.quant-gas].[26] S.-Y. Chä and U. R. Fischer, Probing the scale invarianceof the inflationary power spectrum in expanding quasi-two-dimensional dipolar condensates, Phys. Rev. Lett. , 130404 (2017), arXiv:1609.06155 [cond-mat.quant-gas].[27] R. Saint-Jalm, P. C. M. Castilho, É. Le Cerf, B. Bakkali-Hassani, J. L. Ville, S. Nascimbene, J. Beugnon, andJ. Dalibard, Dynamical symmetry and breathers in a two- dimensional Bose gas, Phys. Rev. X , 021035 (2019),arXiv:1903.04528 [cond-mat.quant-gas].[28] S. Yu. Khlebnikov and I. I. Tkachev, Classical decayof inflaton, Phys. Rev. Lett. , 219 (1996), arXiv:hep-ph/9603378 [hep-ph].[29] J. Berges and T. Gasenzer, Quantum versus classical sta-tistical dynamics of an ultracold Bose gas, Phys. Rev. A76 , 033604 (2007), arXiv:cond-mat/0703163 [cond-mat.other].[30] A. Sinatra, C. Lobo, and Y. Castin, The truncatedWigner method for Bose-condensed gases: limits of va-lidity and applications, J. Phys. B , 3599 (2002),arXiv:cond-mat/0201217 [cond-mat.stat-mech].[31] P. Blakie, A. Bradley, M. Davis, R. Ballagh, and C. Gar-diner, Dynamics and statistical mechanics of ultra-coldBose gases using c-field techniques, Adv. Phys. , 363(2008), arXiv:0809.1487 [cond-mat.stat-mech].[32] C. Barceló, S. Liberati, and M. Visser, Analogue gravity,Living Rev. Relativ. , 3 (2011), arXiv:gr-qc/0505065[gr-qc].[33] P. Jain, S. Weinfurtner, M. Visser, and C. W. Gardiner,Analogue model of a FRW universe in Bose–Einsteincondensates: Application of the classical field method,Phys. Rev. A76 , 033616 (2007), arXiv:0705.2077 [cond-mat.other].[34] K. Staliunas, S. Longhi, and G. J. de Valcárcel, Fara-day patterns in Bose–Einstein condensates, Phys. Rev.Lett. , 210406 (2002), arXiv:cond-mat/0204517 [cond-mat.stat-mech].[35] P. Engels, C. Atherton, and M. A. Hoefer, Observationof Faraday waves in a Bose–Einstein condensate, Phys.Rev. Lett. , 095301 (2007), arXiv:cond-mat/0701028[cond-mat.other].[36] S. E. Pollack, D. Dries, R. G. Hulet, K. M. F. Mag-alhães, E. A. L. Henn, E. R. F. Ramos, M. A. Cara-canhas, and V. S. Bagnato, Collective excitation of aBose–Einstein condensate by modulation of the atomicscattering length, Phys. Rev. A , 053627 (2010),arXiv:1004.2887 [cond-mat.quant-gas].[37] I. Vidanović, A. Balaž, H. Al-Jibbouri, and A. Pel-ster, Nonlinear Bose–Einstein-condensate dynamics in-duced by a harmonic modulation of the s-wave scatteringlength, Phys. Rev. A , 013618 (2011), arXiv:1106.4686[cond-mat.quant-gas].[38] J. C. Jaskula, G. B. Partridge, M. Bonneau, R. Lopes,J. Ruaudel, D. Boiron, and C. I. Westbrook, An acous-tic analog to the dynamical Casimir effect in a Bose-Einstein condensate, Phys. Rev. Lett. , 220401(2012), arXiv:1207.1338 [cond-mat.quant-gas].[39] J. H. V. Nguyen, M. C. Tsatsos, D. Luo, A. U. J. Lode,G. D. Telles, V. S. Bagnato, and R. G. Hulet, Parametricexcitation of a Bose–Einstein condensate: From Faradaywaves to granulation, Phys. Rev. X , 011052 (2019),arXiv:1707.04055 [cond-mat.quant-gas].[40] Z. Zhang, K.-X. Yao, L. Feng, J. Hu, and C. Chin,Pattern formation in a driven Bose–Einstein conden-sate, Nat. Phys. , 652 (2020), arXiv:1909.05536 [cond-mat.quant-gas].[41] J. Berges and J. Serreau, Parametric resonance in quan-tum field theory, Phys. Rev. Lett. , 111601 (2003),arXiv:hep-ph/0208070 [hep-ph].[42] J. Berges and G. Hoffmeister, Nonthermal fixed pointsand the functional renormalization group, Nucl. Phys. B813 , 383 (2009), arXiv:0809.5208 [hep-th]. [43] C.-M. Schmied, A. N. Mikheev, and T. Gasenzer, Non-thermal fixed points: Universal dynamics far from equi-librium, Int. J. Mod. Phys. A , 1941006 (2019),arXiv:1810.08143 [cond-mat.quant-gas].[44] C.-M. Schmied, M. Prüfer, M. K. Oberthaler, andT. Gasenzer, Bidirectional universal dynamics in a spinorBose gas close to a nonthermal fixed point, Phys. Rev. A99 , 033611 (2019), arXiv:1812.08571 [cond-mat.quant-gas].[45] J. A. Glidden, C. Eigen, L. H. Dogra, T. A. Hilker,R. P. Smith, and Z. Hadzibabic, Bidirectional dynamicscaling in an isolated Bose gas far from equilibrium,arXiv:2006.01118 [cond-mat.quant-gas] (2020).[46] A. Mazeliauskas and J. Berges, Prescaling and far-from-equilibrium hydrodynamics in the quark-gluon plasma,Phys. Rev. Lett. , 122301 (2019), arXiv:1810.10554[hep-ph].[47] C.-M. Schmied, A. N. Mikheev, and T. Gasenzer, Prescal-ing in a far-from-equilibrium Bose gas, Phys. Rev. Lett. , 170404 (2019), arXiv:1807.07514 [cond-mat.quant-gas].[48] J. Berges, K. Boguslavski, S. Schlichting, and R. Venu-gopalan, Basin of attraction for turbulent thermaliza-tion and the range of validity of classical-statisticalsimulations, J. High Energy Phys. (5), 54,arXiv:1312.5216 [hep-ph].[49] E. Kolb and M. Turner,
The Early Universe , Frontiers inphysics, Vol. 6 (Addison-Wesley, 1990).[50] L. P. Pitaevskii and S. Stringari,
Bose–Einstein Conden-sation and Superfluidity , 1st ed., International Series ofMonographs on Physics (Oxford University Press, 2016).[51] A. Suárez and P.-H. Chavanis, Hydrodynamic representa-tion of the Klein-Gordon-Einstein equations in the weakfield limit: General formalism and perturbations anal-ysis, Phys. Rev.
D92 , 023510 (2015), arXiv:1503.07437[gr-qc].[52] V. Mukhanov and S. Winitzki,
Introduction to QuantumEffects in Gravity (Cambridge University Press, 2007).[53] M. Olshanii, H. Perrin, and V. Lorent, Example of aquantum anomaly in the physics of ultracold gases, Phys.Rev. Lett. , 095302 (2010), arXiv:1006.1072 [cond-mat.quant-gas].[54] M. Holten, L. Bayha, A. C. Klein, P. A. Murthy, P. M.Preiss, and S. Jochim, Anomalous breaking of scale in-variance in a two-dimensional Fermi gas, Phys. Rev. Lett. , 120401 (2018), arXiv:1803.08879 [cond-mat.quant-gas].[55] C.-L. Hung, X. Zhang, N. Gemelke, and C. Chin, Ob-servation of scale invariance and universality in two-dimensional Bose gases, Nature (London) , 236(2011), arXiv:1009.0016 [cond-mat.quant-gas].[56] P. B. Greene, L. Kofman, and A. A. Starobinsky, Sine–Gordon parametric resonance, Nucl. Phys.
B543 , 423(1999), arXiv:hep-ph/9808477 [hep-ph].[57] Z. Hadzibabic and J. Dalibard, Two-dimensional Bosefluids: An atomic physics perspective, Riv. del NuovoCimento , 389 (2011), arXiv:0912.1490 [cond-mat.quant-gas].[58] N. Bogolyubov, On the theory of superfluidity, J. Phys.(USSR) , 23 (1947).[59] N. W. McLachlan, Theory and application of Mathieufunctions , repr. ed. (Clarendon Pr., 1951).[60] A. H. Nayfeh,
Perturbation methods (Wiley, 1973).[61] V. E. Zakharov, V. S. L’vov, and G. Falkovich,
Kol- mogorov Spectra of Turbulence I: Wave Turbulence ,Springer Series in Nonlinear Dynamics (Springer-VerlagBerlin Heidelberg, 1992).[62] J. Berges, K. Boguslavski, A. Chatrchyan, and J. Jaeckel,Attractive versus repulsive interactions in the Bose–Einstein condensation dynamics of relativistic field the-ories, Phys. Rev.
D96 , 076020 (2017), arXiv:1707.07696[hep-ph].[63] J. Berges, A. Rothkopf, and J. Schmidt, Non-thermalfixed points: Effective weak-coupling for strongly corre-lated systems far from equilibrium, Phys. Rev. Lett. ,041603 (2008), arXiv:0803.0131 [hep-ph].[64] C. Scheppach, J. Berges, and T. Gasenzer, Matter waveturbulence: Beyond kinetic scaling, Phys. Rev.