{'fraction_non_alphanumeric': 0.03791303064406029, 'fraction_numerical': 0.011964282419792263, 'mean_word_length': 4.624520078862717, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Analog cosmological reheating in an ultracold Bose gas

Cosmological reheating describes the transition of the post-inflationary universe to a hot and thermal state. In order to shed light on the underlying dynamics of this process, we propose to quantum-simulate the reheating-like dynamics of a generic cosmological single-field model in an ultracold Bose gas. In our setup, the excitations on top of an atomic Bose-Einstein condensate play the role of the particles produced by the decaying inflaton field after inflation. Expanding spacetime as well as the background oscillating inflaton field are mimicked in the non-relativistic limit by a time dependence of the atomic interactions, which can be tuned experimentally via Feshbach resonances. As we illustrate by means of classical-statistical simulations for the case of two spatial dimensions, the dynamics of the atomic system exhibits the characteristic stages of far-from-equilibrium reheating, including the amplification of fluctuations via parametric instabilities and the subsequent turbulent transport of energy towards higher momenta. The transport is governed by a non-thermal fixed point showing universal self-similar time evolution as well as a transient regime of prescaling with time-dependent scaling exponents. While the classical-statistical simulations can capture only the earlier stages of the dynamics for weak couplings, the proposed experiment has the potential of exploring the evolution up to late times even beyond the weak coupling regime.

Cosmological inflation is a well-established paradigm that solves the puzzles of flatness, homogeneity, and isotropy of the universe [1,2], and explains the generation of density perturbations for structure formation [3]. According to this paradigm, the early universe underwent a stage of accelerated expansion which, in a typical scenario, is driven by a scalar field, the inflaton. Inflation is followed by a reheating phase [4], during which the inflaton decays, e.g., into degrees of freedom of the standard model of particle physics. The heating process can start with a preheating stage of rapid particle production far from equilibrium, and finally completes after the particles, due to their interactions, approach thermal equilibrium at the reheating temperature [5,6]. Over the recent years, it has been realized that many salient features of the dynamics become universal in the far-from-equilibrium regime, across physical systems at vastly different scales [7][8][9][10][11][12]. This universality makes it particularly appealing to gain insight into the nonequilibrium processes during reheating by use of tabletop experiments, such as ultracold atomic gases. Due to their high degree of controllability, these systems are ideally suited to realize and study processes relevant for the early universe that are otherwise challenging to access [13][14][15][16][17][18][19][20][21].
Here, we propose a quantum simulation of reheatinglike dynamics in a single-component Bose gas, leveraging the ability to modulate the atomic interaction in time by means of Feshbach resonances [22]. The description of an ultracold Bose gas formally coincides with the nonrelativistic limit of a scalar field in an expanding Friedmann-Robertson-Walker (FRW) spacetime, where the expansion is encoded in the time dependence of the interaction. This simplifies the simulation of such a system compared to approaches with explicit expansion [23] by lifting the restrictions associated with the expansion of the trap. Remarkably, for the particular case of two spatial dimensions, where the Bose gas exhibits scale invariance [24][25][26], the dynamics becomes independent of the expansion, such that arbitrary expanding backgrounds may be studied by post-processing the data of a single experiment. It is our aim to reproduce the postinflationary dynamics of generic single-field models [7] in an ultracold Bose gas. Relativistic effects, such as particle production from the coherently oscillating inflaton field during the preheating stage, are absent in the atomic system, but preheating-like dynamics can be introduced also in the non-relativistic theory by a sinusoidal modulation of the coupling, mimicking an oscillating background field. As we demonstrate by means of classical-statistical (or truncated Wigner) simulations [27][28][29][30], the proposed setup exhibits the characteristic stages of reheating dynamics [6], including the preheating stage of parametric amplification of quantum fluctuations as well as the self-similar transport of energy towards higher momenta driving turbulent thermalization [7]. In this respect, our  Figure 1. Schematic illustration of post-inflationary reheating dynamics in the early universe and the simulation of an analogous process in an ultracold Bose gas. We consider a scenario where a single-component homogeneous "inflaton" field oscillates around the minimum of its potential, producing particles via parametric instabilities ("preheating"). Later, the system enters a turbulent state where energy is transported towards higher momenta in a self-similar way as the universe approaches thermal equilibrium ("reheating"). In the simulation, expanding spacetime as well as the oscillating inflaton field are mimicked in the non-relativistic limit by modulating the scattering length of a Bose-Einstein condensate (BEC), whose excitations play the role of particles produced by the decaying inflaton.
setup provides an analog to cosmological reheating, as illustrated schematically in Fig. 1.
Analog models of gravity build on emergent curved spacetimes for low-energy excitations on top of Bose-Einstein condensates (BECs) or other systems [31]. These can be generated, for instance, by modulating trapping potentials or interactions, allowing one to simulate cosmological particle production during inflation [13,14,32]. The mapping to an emergent expanding spacetime is, however, valid only for linear excitations, and incorporating non-linear effects into this framework, which are crucial for post-inflationary reheating dynamics and responsible for thermalization, can be challenging. Nonetheless, several experimental setups have been proposed for simulating post-inflationary dynamics in Bose gases [16][17][18][19][20]. Besides these theoretical proposals, preheating-like dynamics has recently been observed in experiments with a ring-shaped condensate after rapidly expanding the trap [23]. In our setup, a preheating-like parametric amplification of sound waves on top of the BEC is induced by modulating the atomic interaction, which mimics the oscillations of the inflaton field around the minimum of its potential [5]. Parametric instabilities in Bose gases have been studied extensively in the literature [33][34][35][36][37][38] and the concept has evolved into a promising tool for state preparation [39]. Going beyond the well-understood linear regime of parametric resonance, our main focus is geared towards the non-linear stages of the dynamics, involving a secondary amplification of excitations [5,19,40] that herald the development of turbulence [7].
The preheating-like instabilities drive the system far away from equilibrium into the vicinity of a non-thermal fixed point, where the dynamics becomes universal and self-similar [41]. A non-thermal fixed point characterizes the far-from-equilibrium behavior of an entire universality class of physical systems as diverse as the inflaton in the early universe, ultracold atomic gases in the laboratory, or even the quark-gluon plasma explored in heavy-ion collisions [7,8,10,42]. Such universal dynamics far from equilibrium has been observed experimentally in spinor and tunnel-coupled Bose gases in form of inverse cascades transporting conserved quantities towards lower momenta [11,12]. Bidirectional cascades involving additionally a self-similar transport of energy towards higher momenta have been studied numerically in one-dimensional (1D) spinor Bose gases [43] and have recently been observed experimentally in an isolated three-dimensional (3D) Bose gas following a cooling quench [44]. Our protocol of parametrically exciting a pure single-component BEC constitutes a complementary path of entering the regime of self-similar dynamics far from equilibrium and allows one to study the direct energy cascade in both driven and isolated systems.
We distinguish a regime of driven turbulence in the presence of continuous parametric driving and a regime of free turbulence after the modulation of the interaction is switched off, corresponding to the decay of the inflaton. Recently, a new type of prescaling phenomenon [45,46] has been established, which is closely related to the onset of a hydrodynamic behavior far from equilibrium [45]. In this regime, the dynamics is governed by the fixed-point scaling function and time-dependent scaling exponents, which slowly relax to their universal values. Remarkably, our numerical results indicate that such a behavior emerges during the transition from driven to free turbulence, which opens the door for an experimental observation of this phenomenon in cold Bose gases.
At late times, when the occupancies at characteristic momenta become of order unity, the system is expected to thermalize [7]. This final stage of the dynamics is dominated by quantum fluctuations and therefore cannot be captured by means of classical-statistical simulations [28,47]. In contrast, it can be accessed by experimental studies in cold atomic systems, which are quantum-mechanical by nature.
For reasons of consistency with standard cosmology, the inflaton must necessarily decay into more familiar forms of matter and gauge fields, e.g., those that constitute the standard model of particle physics, before the universe thermalizes at the reheating temperature. In contrast, reheating in our setup involves the decay of the condensate into its own fluctuations only, as in the relativistic λϕ 4 model that has been studied in Ref. [7].
Despite this simplicity, these models capture many characteristic features of the reheating process, which can be directly probed in our proposed experiment. This paper is organized as follows. In Sec. II, we present our approach of simulating expansion in a Bose gas via a time-dependent interaction, establishing an analogy to the inflaton field in the early universe after inflation. We proceed in Sec. III ("analog preheating") with discussing preheating-like dynamics, where particle production is mimicked by parametrically exciting Bogoliubov quasi-particles via a modulation of the interaction. Though our approach works in principle for any dimensionality, we focus here on the case of a spatially two-dimensional (2D) Bose gas. We explain and interpret the momentum spectrum resulting from these instabilities, emphasizing the important role of non-linear effects leading to the onset of turbulent dynamics. In Sec. IV ("analog reheating"), we analyze the universal self-similar time evolution for both the driven and the free direct cascade, comparing the numerically extracted scaling exponents with predictions from kinetic theory. We also show that the system exhibits prescaling during the transition from driven to free turbulence. Moreover, we discuss how expansion affects the dynamics and can, under certain conditions, even prevent the system from thermalizing. In Sec. V, we discuss details relevant for an experimental implementation, before we conclude in Sec. VI.

II. EXPANDING SPACETIME IN BOSE GASES
The inflaton is commonly described by a relativistic real scalar field φ(x) in curved spacetime with the action [48] where V is the potential, g µν is the metric tensor, g its determinant, and d the number of spatial dimensions. The expansion of the universe is well described with the help of the flat FRW metric [48]. It corresponds to a homogeneous and isotropic universe, in accordance with the cosmological principle, which is Minkowskian at each time slice. The FRW metric has the form where c is the speed of light and τ denotes the cosmic time. The cosmic scale factor a(τ ) grows with time in an expanding universe and relates the comoving distance x to the proper distance r(τ ) = a(τ )x. Our goal is to observe reheating-like dynamics in an ultracold Bose gas in analogy to reheating in the early universe. A Bose gas is described by a non-relativistic complex field operatorΨ subject to the Hamiltonian [49] where is the reduced Planck constant, m the atomic mass, U an external trapping potential, and g the quartic coupling which determines the strength of the atomic interactions via the s-wave scattering length. In what follows, we briefly discuss possible strategies of how to realize an expanding spacetime according to Eq. (2) in an ultracold Bose gas, before we detail the approach chosen in this work. The most straightforward approach to incorporate expansion in a trapped Bose gas is to physically expand the trap geometry [15,23,25]. While this approach is simple and direct in principle, there are practical limitations such as the restriction to short times or small expansion velocities due to a finite optical imaging system.
Alternatively, one can keep the trap geometry fixed, associating physical distances with the co-moving distances of an expanding system. In this way, however, the redshift of momenta as well as dilution of the system due to the expansion have to be explicitly accounted for in form of a modified kinetic term and a non-unitary evolution, as discussed in the subsequent subsection (see also Ref. [50]). Besides possible technical difficulties in engineering such non-Hermitian Hamiltonians, this formulation has the additional drawback that it implies a decreasing atomic density, which diminishes the experimental signal with progressing laboratory time.
We follow here a third approach, wherein the expansion of spacetime is encoded solely in the time dependence of the atomic interaction [13,25,32]. For experiments with ultracold atoms where a broad Feshbach resonance is available to tune the interactions, this third approach may be the preferred way of studying a broad range of expansion scenarios. In what follows, we will discuss our approach of incorporating expansion. In Secs. III and IV, we will then illustrate the dynamics of preheating and reheating that becomes accessible in this setup.

A. Relativistic versus non-relativistic dynamics in expanding spacetime
We model the inflaton as a real scalar field, whose dynamics is described by the action (1) with a potential where m is the mass and λ a quartic coupling.
The classical equations of motion for the inflaton field can be obtained from Eq. (1) using the principle of least action. For the metric defined in Eq. (2) they read where H =ȧ/a is the Hubble parameter and the dot denotes the derivative with respect to the cosmic time τ . This equation is of a Klein-Gordon type with a modified spatial derivative term and an additional friction term, which express the redshift of momenta in the co-moving frame and the dilution of the field due to the expansion, respectively [48].
The analogous equations for the non-relativistic system have a similar structure. To see this, we factor out the fast oscillations due to the mass term in Eq. (4) from the canonically quantized inflaton fieldφ(τ, x), defining a slowly varying complex fieldψ(τ, x) via the relation where h.c. denotes the Hermitian conjugate. As shown in Appendix A, provided typical momenta, expansion velocities, and field values are small, this new field evolves according to the equation [50] i with coupling g = λ 4 /8m 2 c. Equation (6) is reminiscent of the Heisenberg equations of motion generated by the Hamiltonian of an ultracold Bose gas (3) in absence of an external trapping potential. The kinetic term is proportional to a −2 , which describes the redshift of momenta in the co-moving frame, and the dilution of the system due to the expansion is expressed by a non-Hermitian term causing the norm of the field to decay. In analogy to the transformation to conformal variables known in quantum field theory on curved spacetime [51], we now introduce a new time variable t via the relation dt = dτ /a 2 , which we refer to as the laboratory time, and the rescaled field operatorΨ =ψa d/2 . Equation (6) then takes the standard form of the equations of motion for a non-relativistic bosonic field generated by the Hamiltonian (3), with the time-dependent effective coupling The transformation back to the original variables is performed in post-processing. Equations (7) and (8) describe a Bose gas in an expanding spacetime, where the expansion enters only in the time dependence of the interaction term via the scale factor a(t). We stress that the correspondence to Eq. (4) is restricted to the non-relativistic limit and therefore not capable of capturing relativistic effects, such as the resonant amplification of fluctuations during the preheating stage. As discussed in the next section, we can mimic an analogous process in the non-relativistic system by periodically modulating the interaction. Despite these restrictions, it is worth emphasizing that the relation between Eqs. (4) and (7) holds on the non-linear level, i.e., it does not rely on any linearization. This is in contrast to the analog gravity framework, where a mapping to Eq. (4) is achieved by engineering an emergent expanding spacetime for linear excitations on top of the condensate [32]. Such a mapping, however, typically breaks down on the non-linear level. Thus, our approach based on Eq. (7) is well suited for simulating reheating-like dynamics, where non-linear effects are essential.

B. The special case of expansion in 2D
For d = 2, the effective coupling in Eq. (8) becomes independent of the scale factor a(t). This is a consequence of a dynamical symmetry in 2D Bose gases known as scale invariance [24]. While quantum anomalies in strongly interacting systems can lead to violations of scale invariance [52,53], it has been well tested experimentally in the weakly interacting regime [26,54]. Thus, if scale invariance holds, the equations of motion are the same as in the case of a static spacetime. The nature of the expansion, encoded in the scale factor a(t), then enters the transformation only back to the original temporal and spatial coordinates as well as field variables. This makes the simulation in two spatial dimensions particularly efficient, since the evolution of a single experiment can conveniently be mapped to arbitrary expanding spacetimes in a post-processing step. For this reason, we will focus on weakly interacting 2D Bose gases in the remainder of this work. Although our universe is clearly not 2D, this geometry captures most of the essential physics of reheating dynamics, as we demonstrate in the following sections.

III. ANALOG PREHEATING
Having discussed the connection between a scalar field in an expanding spacetime and an ultracold Bose gas, we are now in the position to address the question of how an analog of the preheating stage of post-inflationary dynamics can be simulated in a Bose gas.

A. Preheating in the early universe
In this subsection, we briefly review the basic mechanism of preheating in the early universe. After inflation, the energy budget of the universe is assumed to be stored predominantly in the homogeneous inflaton field, which oscillates around the minimum of its potential. The stage of preheating consists of an explosive, non-perturbative particle creation from the decaying inflaton.
A common mechanism for preheating is via the parametric amplification of quantum fluctuations in the presence of the effective potential induced by the inflaton [5,40]. This effect can most conveniently be understood by linearizing the fluctuations of the inflaton field around its homogeneous background, φ(τ, x) = φ 0 (τ ) + δφ(τ, x), and inserting it in Eq. (4). The resulting equation describes damped oscillations of the background field, The equations of motion for the fluctuations, dropping all terms of quadratic or higher order in δφ, read in Fourier space with ω 2 p (φ 0 , τ )/c 2 = p 2 / 2 a 2 (τ ) + m 2 c 2 / 2 + λφ 2 0 /2 and the Fourier transform δφ p (τ ) = d d x δφ(τ, x)e −ipx/ . Equation (10) describes a collection of parametric oscillators for each momentum mode p, driven by the background oscillations of φ 0 (τ ). Modes within certain instability bands satisfy resonance conditions, which leads to an exponential growth of occupancies corresponding to particle production [40,55]. The growth of fluctuations eventually invalidates the linearized approach of Eq. (10). Note that in more realistic models, couplings of the inflaton to other bosonic fields can lead, in a similar way, to a production of those degrees of freedom via parametric resonances as well [5]. Here, we focus on the simplest scenario corresponding to a decay of the inflaton solely into its own quanta of excitation.

B. Analog preheating in Bose gases
Our aim is to observe an analog of preheating dynamics in a cold-atom Bose gas. Motivated by the discussions in the previous sections, we mimic the state of the universe after inflation by a homogeneous BEC. (Possible implications of the presence of an external trapping potential are discussed in Sec. V.) It is tempting to search for a direct analog of the inflaton's oscillations, which could trigger parametric instabilities in the condensate. However, this effect is not present for a single-component Bose gas in its ground state.
In fact, the absence of parametric resonance is a consequence of the non-relativistic limit considered in Sec. II. To re-introduce parametric instabilities in the non-relativistic model, we add a periodic modulation of the coupling with frequency ω as g → g(1 + sin ωτ ). As discussed in Appendix B, this procedure does not map one-to-one to the parametric resonance scenario described by Eq. (10). For instance, the modulation frequency ω is non-relativistic in the simulation, while according to Eq. (9), the inflaton oscillates at relativistic frequencies ω rel mc 2 / on the scale given by its rest mass. Furthermore, there is no dynamical mechanism of back reaction of the fluctuations on the modulation frequency or amplitude, as contained in the full model (4).
Despite these differences, modulating the interaction as described above suffices to generate the desired parametric resonance phenomena in the Bose gas [33][34][35][36][37][38][39], analogous to the ones taking place in the preheating stage of most inflationary models [5,40]. This is demonstrated in detail in the subsequent subsections.

C. Numerical study of preheating dynamics
We illustrate our implementation of preheating dynamics by means of numerical simulations based on Eq. (7), which corresponds to a static ultracold Bose gas with a time-dependent interaction. We consider here a Bose gas in two dimensions, where the expansion does not explicitly enter the equations of motion, as discussed in Sec. II B. (Nonetheless, we keep the number of spatial dimensions d in all formulas general.) Equation (7) is then formally equivalent to the equation of motion generated by the Hamiltonian (3), and the back transformation to the original time and field variables can be performed in post-processing. The latter is illustrated at the example of a matter-dominated expansion in Sec. IV E.
Following the argumentation in Sec. III B, we induce parametric instabilities by harmonically modulating the coupling as where g 0 is a positive offset value, r is the amplitude of the modulation and ω is its frequency. Note that here we have used the laboratory time t in the argument of the sine function instead of the cosmic time τ , which, when expressed in terms of the laboratory time t, results in oscillations with increasing frequency. This substitution is based on the simplifying assumption that the expansion is insignificant during the duration of the modulation, in which case the relation between both time variables becomes linear and the corresponding proportionality constant can be absorbed in the modulation frequency.
In contrast to the post-inflationary setting, where the produced particles back-react on the inflaton, the modulation of the coupling is imposed externally in the simulation. Switching off the modulation thus corresponds to the decay of the inflaton, which constitutes a free parameter in the model.
A particularly useful observable for the study of parametric resonance and turbulent thermalization is the single-particle momentum distribution where V is the volume andΨ p (t) = d d xΨ(t, x)e −ipx/ denotes the Fourier transform of the Bose fieldΨ(t, x). Importantly, this quantity can be experimentally accessed in time-of-flight measurements (see Sec. V for further discussions). We evaluate the quantum expectation value in Eq. (12) by means of classical-statistical (or truncated Wigner) simulations [27][28][29][30]. This method takes into account quantum fluctuations by stochastically sampling classical field configurations from the Wigner distribution of the initial state. Each realization is propagated deterministically according to the Gross-Pitaevskii equation (GPE) and quantum mechanical observables are obtained as statistical averages over multiple realizations. Here, the initial state is taken to be a spatially homogeneous BEC, Momentum distribution f (t, p) + 1 2 p =hω/2 Radially averaged momentum distribution f (t, p) as a function of the radial momentum p = |p| at different times t, demonstrating preheating dynamics. The coupling is modulated with relative strength r = 0.25 at a frequency ω, chosen such that the resonance condition pres = ω/2 for the momentum pres = 12 × pL with pL = 2π /L is fulfilled [see discussion below Eq. (16)]. At early times t 30 × t0, a single narrow resonance can be observed around pres. At later times, a broad resonance band emerges with peaks at higher harmonics of the modulation frequency. These secondary resonances are due to non-linear interactions, as discussed in Sec. III E. See the video accompanying the arXiv version of this paper for a qualitative illustration of the dynamics of a single realization [56].
and for each realization, all non-zero momentum modes are populated with vacuum noise corresponding to an average occupancy of half a particle per mode. This mimics quantum fluctuations acting as a seed for parametric instabilities. It is important to note that this approach goes beyond a mean-field description, which fails to capture parametric resonance since the occupancies of all excited modes are exactly zero. Further details about the simulation method can be found in Appendix C.
In what follows, we express length and time in units of the characteristic scales x 0 = / √ mn 0 g 0 and t 0 = /n 0 g 0 , respectively, where n 0 = N/V is the homogeneous particle density in a system of N particles in a volume V . Momenta are given either in units of the lowest non-zero momentum p L = 2π /L or in units of the characteristic momentum p ξ = 2π /ξ corresponding to the healing length ξ = / √ 2mn 0 g 0 . In a quasi-2D Bose gas, the interaction strength g 0 =g 2 /m is characterized by the dimensionless parameterg = √ 8πa s /a HO , where a s is the s-wave scattering length and a HO is the oscillator length of the harmonic potential in the strongly confined direction [57]. If not stated otherwise, we consider a uniform quasi-2D system of N = 10 6 weakly interacting particles with couplingg = 2.5 × 10 −3 in a square box with periodic boundary conditions. This choice of parameters fixes the box length as L/x 0 = 50. Moreover, in this isotropic setting, the momentum distribution depends only on the magnitude p = |p| of the momentum, Figure 2 shows the radially averaged momentum distribution of a parametrically excited system described by Eq. (7). The interaction is modulated according to Eq. (11) with r = 0.25, and the modulation frequency ω is chosen as the resonance frequency of the momentum |p res | = 12 × p L , as discussed below Eq. (16) in the next subsection. In this simulation, the noise cutoff has been chosen as p Λ = p ξ ≈ 70.7 × p L . To gain a qualitative impression of the induced preheating dynamics, see the video accompanying the arXiv version of this paper, where we show the evolution of both the density and the phase of a single realization [56].
At early times, t 30 × t 0 , in Fig. 2, we observe a single narrow resonance around the resonance momentum p res satisfying the resonance condition pres = ω/2 (cf. Sec. III D). Due to particle number conservation, the growing occupancy of the resonant momentum causes the condensate to decay, mimicking particle production from the decaying inflaton in the early universe. At later times, secondary resonances at higher harmonics of the modulation frequency appear as a result of non-linear interactions among the produced quasi-particles.
In addition to the narrow resonance peaks, a transient growth of fluctuations at low momenta occurs at early times before the primary peak becomes visible. This growth can be interpreted as a consequence of the fact that the sampled initial state corresponds to the ground state of an ideal Bose gas. At t = 0, the system is effectively quenched to a finite interaction, producing a power law in the momentum distribution proportional to p −1 at low momenta (cf. Appendix C 1). This early-time behavior has, however, no influence on the preheating dynamics we are interested in here.
In the following subsections, we provide analytical insights into both the linear regime of parametric resonance and the non-linear regime of secondary resonances.

D. The linear regime of parametric resonance
Parametric instabilities play a crucial role in many modern experiments with BECs [33][34][35][37][38][39] and can conveniently be understood by adopting a semi-classical point of view. To this end, we consider the GPE for the condensate wave function Ψ(t, x), which is formally obtained as the classical equation of motion after replacing the Bose field operator in Eq. (7) by its expectation value, It is instructive to work in the Madelung representation, Ψ(t, x) = n(t, x) exp (iθ(t, x)), which allows us to express the GPE (13) in form of the hydrodynamic equa- To obtain some intuition about the early stages of the evolution, we express both the density and the phase in terms of a homogeneous background with fluctuations on top of it, n(t, x) = n 0 (t) + n 1 (t, x) and (14) with respect to the fluctuations, yields the equations ∂ t n 0 = 0 and ∂ t θ 0 = −g(t)n 0 / for the background condensate.
For the fluctuations, we obtain Taking the time derivative and inserting the resulting equations into each other, fluctuations of the density and the phase decouple to linear order. Transforming to momentum space, n 1 p (t) = d d x n 1 (t, x)e −ipx/ and similarly for the phase, the linearized equations can be expressed as ) is a time-dependent form of the famous Bogoliubov dispersion relation [58], and p,0 = p 2 /2m denotes the dispersion relation of a free particle.
Equation (16) describes a collection of parametric oscillators for each momentum mode p, which are undamped for the density and damped for the phase. As discussed in Appendix D, these are special cases of Mathieu's equation [59], which admits oscillatory solutions with exponentially growing amplitudes ∼ e ζpt , describing parametric resonance [40,55]. The resonance condition for the momentum mode p res reads pres = ω/2, where p = p,0 ( p,0 + 2n 0 g 0 ) denotes the Bogoliubov dispersion relation and ω is the modulation frequency defined in Eq. (11). That is, resonance occurs for those momentum modes p res whose energy equals half a quantum of energy ω/2 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 of this instability band increases with the modulation amplitude r. To leading order in perturbation theory [60] (see Appendix D), the growth rate of the resonant momentum mode is given by For ω n 0 g 0 , this rate simplifies to ζ pres ≈ rω/8. In this regime, the Bogoliubov dispersion becomes linear, p ≈ c s |p|/ with the speed of sound c s = n 0 g 0 /m, and the produced quasi-particles have the character of sound waves. In the opposite limit, ω n 0 g 0 , particles with quadratic dispersion p ≈ n 0 g 0 + p 2 /2m are produced.
In the simulation presented in Fig. 2, parametric resonance is clearly visible as a pronounced peak at the momentum satisfying the resonance condition. Likewise, the excitation of a single dominant wave length in the linear stage of the dynamics is qualitatively confirmed in the evolution of both the density and the phase of a single realization, as can be seen in the video accompanying the arXiv version of this paper [56] as well as in the snapshots 1 shown in the central panels in the lower part of Fig. 1. It is worthwhile emphasizing that parametric instabilities can be triggered only if the initial occupancy is non-zero. This seed is not contained in the mean field analysis presented in this subsection, but is added in the simulation in form of vacuum noise according to the truncated Wigner prescription.
The linearized equations (15) are helpful to get an intuitive analytical understanding for the early stages of the dynamics and describe the emergence of the primary resonant peak in Fig. 2. However, as a result of the exponential growth of occupancies, this approach fails to describe the later stages where non-linear effects play a fundamental role. These non-linearities are taken into account by our numerical simulations, which are based on the full GPE (13), and include secondary excitations outside the resonance band, as shown in Fig. 2. These will be discussed further in the subsequent section.

E. Secondary instabilities
In this section we discuss the first non-linear corrections, which lead to secondary instabilities [19,40]. Within our semi-classical picture, these can be understood by considering the hydrodynamic equations of motion (14) for the fluctuations to quadratic order, i.e., including terms O(θ 2 1 ), O(θ 1 (n 1 /n 0 )) and O((n 1 /n 0 ) 2 ). This leads to 1 Since n 1 and θ 1 are conjugate variables, they are phase shifted such that the density fluctuations reach their maximum when the phase fluctuations cross zero, and vice versa. To create a better visual impression, the central snapshot of the phase in Fig. 1 has been shifted forward in time by a quarter of an oscillation period until the phase fluctuations reach their next maximum.
As can be seen, density and phase fluctuations no longer decouple to this order. Taking the time derivative of the above equation for the density fluctuations and inserting the expressions for the first-order time derivatives of the fluctuations, the result can be written in momentum space as In deriving Eq. (19), we have neglected all cubic terms in the fluctuations. The above momentum integrals are dominated by the contribution from the exponentially growing unstable modes q, for which |p − q| ≈ |q| ≈ p res . For modes that are stable on the linear level, the integrals therefore act as source terms and Eq. (19) describes forced harmonic oscillators with an exponentially growing force. This effect is analyzed in more detail in Appendix E. There, we show that the forcing leads to an exponential growth of the momentum modes with p 2p res proportional to e 2ζp res t , with a growth rate ζ pres given by Eq. (17). Modes, for which in addition p ≈ 2 pres holds, experience a resonant amplification and are strongest enhanced.
These features of secondary instabilities are captured by our numerical simulations. In particular, one can observe both the narrow peak and the broad band in the distribution function in Fig. 2 between t = 46 × t 0 and t = 58 × t 0 . Similar peaks at higher multiples of the resonance frequency appear at later times due to higher-order corrections to Eq. (14).
The perturbative analysis breaks down when n 1 /n 0 ≈ 1, i.e., when the number of excited atoms becomes comparable to the number of condensate atoms. At this point, the exponential growth stops and turbulent dynamics sets in, which can be observed qualitatively in the video accompanying the arXiv version of this paper [56]. A typical snapshot of the density and the phase of a single realization after the onset of turbulence is shown in the right panels in the lower part of Fig. 1. In the next section, we present a quantitative analysis of turbulent dynamics in momentum space.
In Appendix F, we discuss the theory of parametric instabilities and secondary excitations from the perspective of quantum equations of motion for correlation functions, which validates the use of the classical-statistical approximation.

F. Analysis of growth rates
agrees well with the analytical prediction 2ζ (17). The growth of the secondary instability starts later, but its rate ζ is approximately twice as large as the one of the primary resonance, as expected from our discussion in Sec. III E. The exponential growth stops when the number of excited atoms becomes comparable to the number of condensate atoms. annotated in Fig. 2. The growth of the secondary resonance marks the onset of the non-linear regime where the quasi-particles produced in the primary resonance start to interact, and the linear picture presented in Sec. III D breaks down.
Before comparing the numerically extracted growth rates to the analytical predictions presented in the previous sections, we first relate the momentum distribution f (t, p) to the hydrodynamic density and phase variables.
On the mean-field level, linear fluctuations of the condensate wave function, expressed as Ψ(t, x) = Ψ 0 (t) + Ψ 1 (t, x), are related to linear density and phase fluctuations via x). The momentum distribution on the linear level then corresponds to Therefore, a parametric resonance where the density and phase fluctuations grow as n 1 p ∼ θ 1 p ∼ e ζpt results in a growth of the momentum distribution as f (t, p) ∼ e 2ζpt . The oscillations of the occupancies in Fig. 3 can be understood from the linearized parametric oscillator equations (16). Recall that the latter admit oscillatory solutions with exponentially growing amplitudes. However, being conjugate variables, the oscillations of n 1 and θ 1 are shifted in phase by approximately π/2. According to Eq. (20), the momentum distribution thus corresponds to the sum of two phase-shifted oscillating functions with slightly different initial amplitudes. This results in the residual oscillations on top of the exponential growth observed in Fig. 3. We have checked that the oscillation frequency of the primary resonance agrees with the modulation frequency, while the oscillations of the secondary resonance additionally contain frequency components corresponding to twice the modulation frequency, reflecting the interactions between the resonantly produced quasi-particles.
We have extracted the growth rate of the primary resonance ζ 0 by fitting an exponential function to the numerical data, as shown in the insets of Fig. 3. The result is close to the analytical prediction (17) obtained from perturbation theory, 2ζ 0 . The secondary resonance at 2 pres in Fig. 3 grows at the rate ζ 0 , which is indeed approximately twice the growth rate of the primary resonance. The secondary resonance quickly saturates, indicating the breakdown of the perturbative analysis presented in the previous section. This marks the transition into a highly non-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-like dynamics along with a possible implementation in an ultracold Bose gas, we now turn to the reheating process and the question of how the system finally approaches thermal equilibrium.
Particle spectra formed during preheating are highly non-thermal with large occupation numbers at low momenta. As a consequence, the system enters a turbulent state, characterized by a local transport of conserved quantities in momentum space. Typically, an inverse cascade transports particles towards lower momenta, while a direct cascade transports energy towards higher momenta, constituting a key process in the context of turbulent thermalization [7]. At early times, this transport is driven, i.e., the oscillating inflaton acts as a source injecting energy into the system at resonant momenta. Eventually, the inflaton decays, marking the transition from driven to free turbulence. The system remains in the turbulent state for a long time, until the occupancy of characteristic momenta eventually becomes comparable to the vacuum expectation value given by the "quantum half". In this final stage, which is dominated by quantum fluctuations, the system relaxes to thermal equilibrium, completing the reheating process [6,7]. Turbulent dynamics is accompanied by the emergence of self-similarity and universality. This is reflected by a power-law behavior of the momentum distribution within a certain inertial range of momenta. One well-known example within the theory of weak wave turbulence is the prediction of a stationary direct cascade with a universal power-law distribution f (p) ∝ p −d [61]. More generally, self-similarity in far-from-equilibrium systems can become manifest in their time evolution. In such a scenario the evolution of the distribution function can be expressed as where s = t/t ref and t ref is an arbitrary reference time.
The scaling hypothesis (21) constitutes a significant reduction of complexity as it allows to describe a relevant part of the system's dynamics by simply rescaling a single-variable scaling function f S as determined by the scaling exponents α and β. Remarkably, in many farfrom-equilibrium scenarios, both the exponents as well as the scaling form of the distribution are universal, which means they are insensitive to microscopic details as well as initial conditions, and depend only on a few macroscopic system parameters like dimensionality, symmetry or the number of field components [7,8,10,42,62]. As a result, if universality holds, physical systems with vastly differing energy scales can behave quantitatively the same. This makes ultracold Bose gases a particularly promising target for simulating universal aspects of far-from-equilibrium dynamics like that of the inflaton in the early universe. Universal self-similar time evolution reflects the system being in the vicinity of a non-thermal fixed point, which acts as an attractor on the way towards thermal equilibrium [41,63,64].
In this section, we demonstrate by means of classicalstatistical simulations that the salient features of these phenomena can be observed in our implementation of reheating dynamics. We first study the regimes of driven and free turbulence separately, before considering a transient prescaling regime, where the universal shape of the scaling function is maintained, but the scaling exponents change over time. From our numerical simulations, we extract the scaling exponents as well as the scaling form of the distribution and compare the results with analytical predictions from kinetic theory. Finally, we discuss the relaxation to thermal equilibrium, which is, however, not captured by our classical-statistical simulations as it is dominated by quantum fluctuations. We conclude this section with a discussion of an example of how expansion may prevent the system from thermalizing, building a bridge to Sec. II.

A. Driven versus free turbulence
In order to drive the system into a turbulent state, we follow the protocol presented in Sec. III C of parametrically exciting a homogeneous BEC. Here, our focus lies on the later stages of the non-linear dynamics after the proliferation of secondary instabilities when a smooth distribution in form of a power law has formed. We distinguish the regime of driven turbulence, realized by continuously modulating the interaction according to Eq. (11), and free turbulence, developed if the modulation is switched off shortly after the primary resonance has saturated. In our analogy to reheating in the early universe, the former case corresponds to the situation where the inflaton possesses enough energy to drive turbulence for a long time, while in the latter case, the inflaton runs out of energy rather quickly at around the same time when turbulence sets in.
To maximize the inertial range where self-similar scaling can be observed, it is desirable to inject energy at momentum scales close to the lowest momenta supported by the system where occupancies can become large. To this end, the interaction is modulated with a relative amplitude r = 1 at frequency ω chosen such that the resonant momentum becomes p res = 3 × 2π /L (cf. Sec. III D).
The left and right panels of Fig. 4 show a comparison of the direct cascades emerging in the regimes of driven and free turbulence, respectively. All parameters for both simulations are identical, with the exception that in Fig. 4(a), the interaction is modulated continuously, while in Fig. 4(b), the modulation is switched off smoothly 2 within two modulation periods 2π/ω ≈ 8.2 × t 0 at time t = 80 × t 0 , roughly corresponding to the time when the primary resonance saturates. In the latter case, energy is conserved already at the onset of turbulence and the driven regime is skipped. In both scenarios, the momentum distribution takes a scaling form corresponding to a power law close to f (t, p) ∝ p −2 , which is indicated by a dotted line as a guide to the eye. A power law proportional to p −d is characteristic for weak wave turbulence [61,65,66] and has been observed experimentally in Ref. [67].
Moreover, as it can be seen in Fig. 4, the distributions in the two regimes exhibit self-similar time evolution in different ways. In the case of driven turbulence, Fig. 4(a), the front of the cascade evolves self-similarly, leaving behind a stationary distribution. Stationary turbulence arises in the theory of weak wave turbulence as a stationary solution of the scattering integral. However, such a configuration necessarily requires the existence of at least a source and a sink for energy to be injected and dissipated, respectively [61]. In the present case, the energy source is provided by the modulation of the interaction. We have verified that the rate of energy injection into the system is approximately constant in this stationary regime, as expected for driven turbulence [7]. Since energy is transported locally in momentum space, unoccupied higher momentum 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. By contrast, the distribution in the case of Momentum distribution f (t, p)   Energy is injected at low momenta by modulating the scattering length according to Eq. (11) with a relative amplitude r = 1 at a frequency ω chosen such that ω/2 = pres with pres = 3 × 2π /L. In the case of continuous modulation (a), a stationary distribution with a power law close to p −2 develops, whose front is evolving self-similarly. If the driving is switched off once the primary resonance has saturated, corresponding here to t = 80 × t0 (b), energy is propagated in a selfsimilar way to higher momenta, but the distribution at a given momentum decreases with time, reflecting energy conservation. A power law proportional to p −2 is shown in form of a dotted line as a guide to the eye. The insets show the distributions rescaled according to Eq. (21) using the numerically extracted scaling exponents displayed below the respective distributions. free turbulence, Fig. 4(b), is not stationary, but decreases as a function of time for a given momentum, reflecting energy conservation.
To quantify the self-similar time evolution, we have extracted the scaling exponents α and β defined in Eq. (21) from our numerical data at different times using the maximum likelihood technique described in Appendix G. In the turbulent stage of the dynamics, the exponents are approximately constant at late times, as shown in the lower part of Fig. 4. The observed slow relaxation of the exponents at early times in Fig. 4(b) can be interpreted as prescaling, as discussed in the subsequent subsection. At the latest simulated times, the exponents take the values α driven = −0.95 ± 0.17, β driven = −0.46 ± 0.06, (22) for driven turbulence, and α free = −1.05 ± 0.11, β free = −0.27 ± 0.03, for free turbulence, respectively. Rescaling the distribution according to Eq. (21) with the extracted exponents, as described in Appendix G, all data points collapse to a single universal scaling function f S (p), as shown in the insets of Fig. 4. The exponent β describes the speed of energy propagation towards higher momenta, which is higher for the driven cascade than for the free cascade. Analytical predictions for this exponent from kinetic theory as well as the relations between α and β following from conservation laws will be discussed in Sec. IV C.
The above values of the exponents are insensitive to the details of how the far-from-equilibrium state is approached. In particular, we have observed the same exponents starting from an initial state with a highly occupied narrow window of momenta on top of a condensate background. Such an initial state is similar to the state of the system at the end of the preheating stage, when a 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 simulations for experimentally realistic configurations, since this method is known to be prone to instabilities caused by the vacuum noise [30]. These instabilities manifest in a decay of the "quantum half" and the generation of spurious quantum pressure, resulting in an unphysical dependence of the results on the ultraviolet (UV) cutoff [47]. Classical statistical simulations are therefore restricted to large occupancies and weak couplings, where these inevitable deficiencies are mitigated via a separation of scales. To simulate reheating dynamics in the turbulent regime, we have therefore increased the particle number and reduced the coupling such that the validity of the classical-statistical approximation can be ensured. This is discussed in detail in Appendix C 3, where we also assess the range of accessible coupling strengths for our setup. In Sec. V, we discuss how our numerical results relate to realistic experimental conditions.

B. Prescaling
So far, we have investigated the two limiting cases where turbulence is either driven or free. Now, we address the transient regime corresponding to the somewhat more realistic situation where, in the beginning, turbulence is driven by the inflaton oscillations, but at some point goes over to free turbulence when the inflaton has decayed.
This scenario is illustrated in Fig. 5. Up to the time t = 256 × t 0 , the direct cascade is driven (blue curves), as in Fig. 4(a). At this time, the modulation is switched off, mimicking the decay of the inflaton. We then observe a slowing down of the speed of energy propagation and the distribution decreases in time for a given momentum (red curves), reminiscent of the direct cascade of free turbulence shown in Fig. 4(b). On the right-hand side of Fig. 5, the numerically extracted scaling exponents are shown as a function of time. After switching off the modulation, the exponent β slowly evolves from a value close to the one reported in Eq. (22) for driven turbulence to a value close to the one reported in Eq. (23) for free turbulence. By contrast, the ratio α/β changes rather quickly between these two regimes. This behavior is expected since this ratio is fixed by energy conservation, which is enforced instantaneously after switching off the modulation, cf. Sec. IV C. Surprisingly, the self-similar scaling form of the distribution is approximately preserved during the transition, which is clearly visible in the inset of Fig. 5, where the distributions, rescaled according to Eq. (21) with the time-dependent scaling exponents, fall on top of each other.
Recently, such a situation, where the system's dynamics is governed by a universal scaling function much before the corresponding exponents have attained their universal values, has been studied in the context of heavy-ion collisions [45]. This phenomenon, termed prescaling, is closely related to the emergence of a far-from-equilibrium hydrodynamic behavior, as it allows to describe the dynamics in terms of few slowly changing parameters [45]. A different type of prescaling, where certain correlation functions already scale with their universal exponents at early times while others do so only at much later times, has been studied numerically in three-component Bose gases [46]. Our results indicate that prescaling, as defined in Ref. [45], can be observed during the transition from driven to free turbulence 3 , opening up new paths to studying this phenomenon experimentally.
In Appendix H, we corroborate our analysis of prescaling using an alternative method based on the moments of the distribution, which is particularly suitable for extracting time-dependent scaling exponents [45].

C. Kinetic description
We now discuss how the values of the scaling exponents for the direct cascade presented in Eqs. (22) and (23) can be understood from an analytical point of view.
One relation between the exponents α and β follows from the scaling of the total energy. To see this, we assume a self-similar time evolution according to Eq. (21) as well as a power-law scaling of the dispersion relation, p ∝ |p| z , where z ≈ 2 in the particle regime of the Bogoliubov dispersion. Then, the total kinetic energy scales as (24) In the case of free turbulence, energy is conserved across the cascade, E = const, while in the case of driven turbulence, it grows linearly with time, E ∝ t. This implies the relation where γ = 1 in the driven regime and γ = 0 in the absence of driving. Another relation between the scaling exponents can be obtained with the help of perturbative 4 kinetic theory for Bose gases [10,65]. The time evolution of the distribution Momentum distribution f (t, p)  Fig. 4, but the modulation is switched off suddenly at a later time t = 256 × t0. Before this time (blue curves), turbulence is driven and both the momentum distribution as well as the scaling exponents α and β are the same as in Fig. 4(a). After switching off the modulation (red curves), the ratio α/β quickly changes to the one expected for free turbulence, reflecting energy conservation. The exponent β gradually changes towards the value obtained for free turbulence in Fig. 4(b), reducing the speed of energy transport in the cascade. Although the scaling exponents still change in time, the distribution has already attained its universal scaling form. This important hallmark of prescaling is indicated in the inset, where all data points collapse to a single curve after rescaling according to Eq. (21) with the extracted time-dependent scaling exponents α(t) and β(t), as described in Appendix G.
function in this framework is written in the form of a Boltzmann equation (see, e.g., Ref. [68]), where describes 2 ↔ 2 collisions between non-condensate atoms and arises in the presence of the condensate and describes effective 1 ↔ 2 collisions between the non-condensate atoms. Here, 0 and p are the energies of condensate and non-condensate atoms, respectively, and we have used the abbreviations f p = f (t, p) and p = d d p/(2π ) d . In our case, 0 can be neglected since we are interested in the direct cascade, which predominantly involves atoms that are much more energetic than condensate atoms.
In order to derive the desired relation between the scaling exponents, let us assume that the right-hand side of the Boltzmann equation (26) is dominated by one of the two scattering integrals, and denote the number of particles involved in the corresponding scattering processes by l, i.e., we have l = 4 for 2 ↔ 2 scattering and l = 3 for 1 ↔ 2 scattering. Comparing how the left-and righthand sides of Eq. (26) scale with time, and assuming large occupation numbers, f p 1, one arrives at [10] ( Combining this result with Eq. (25) yields For driven turbulence, γ = 1, and a quadratic dispersion, z = 2, both collision terms lead to α = −1 and β = −1/2. These values are close to those in Eq. (22), extracted from our simulation in the presence of a modulated coupling. For free turbulence, γ = 0, one recovers the ratio α/β = 4 reflecting energy conservation, which is close to the numerical result in Eq. (23). The value of β in this case is β = −1/6 for 2 ↔ 2 scattering and β = −1/4 for 2 ↔ 1 scattering. The latter is closer to Eq. (23), which is observed at late times in our simulations. By contrast, in the experiment reported in Ref. [44], where, instead of parametrically driving a pure condensate out of equilibrium, a cooling quench has been applied to an initially uncondensed Bose gas, a scaling exponent closer to the 2 ↔ 2 scattering solution, β = −1/6, has been observed. Deviations from the predicted scaling can in general be induced by a non-quadratic scaling of the dispersion relation p or by the time dependence of the number of atoms in the condensate.

D. Thermalization
As we have seen, on its way to thermal equilibrium the system takes a detour via a non-thermal fixed point, in the vicinity of which the dynamics is dominated by a turbulent transport of energy towards higher momenta. The self-similar dynamics of the direct cascade stops once the occupancy of the characteristic momentum dominating the kinetic energy budget becomes comparable to the expectation value of the vacuum noise given by the "quantum half". At this point, quantum fluctuations become dominant over statistical fluctuations [7] and the system is expected to relax to a Bose-Einstein distribution where p is the dispersion relation, k B is the Boltzmann constant and T f is the reheating temperature. The dominance of quantum fluctuations makes this final stage of the dynamics inaccessible to classicalstatistical simulations. In contrast to the expected relaxation to a Bose-Einstein distribution, at sufficiently late times in the numerical simulation, the cascade stops being self-similar, slows down, and approaches a classical equilibrium distribution [8] with a temperature T (Λ) that depends on the UV cutoff Λ and is determined by the equipartition theorem, f cl,th (p) + 1/2 ∝ k B T (Λ)/ p . In particular, mode occupancies drop to unphysical values below the vacuum noise of 1/2. The late-time regime dominated by quantum fluctuations is beyond the regime of classical-statistical simulation methods. Here, we resort to analytical estimates for the reheating time and the reheating temperature in our system based on the self-similar time evolution. Following Ref. [7], we neglect the final stage of quantum relaxation to a Bose-Einstein distribution and consider thermalization as complete once the occupancy of the characteristic momentum dominating the kinetic energy budget becomes on the order of unity. More precisely, we define the characteristic momentump(t) as the momentum that maximizes the integrand in the expression (24) for the total kinetic energy. For an isotropic system with dispersion relation p ∝ p z , where p = |p|, the characteristic momentum is given byp(t) = arg max p p d+z −1 f (t, p).
If the scaling exponents as well as the momentum distribution at some reference time are known, the assumption of self-similar time evolution according to Eq. (21) is sufficient to predict the time t f when the occupancy of the final characteristic momentum p f will reach unity.
The time t f can be regarded as the best possible approximation to the reheating time obtainable from classical-statistical simulations. However, Eq. (21) is not directly applicable since the scaling exponents are not constant throughout the entire time evolution. In fact, they take different universal values in the regimes of driven and free turbulence, respectively, which are interpolated during a transient regime of prescaling. Importantly, as shown in Sec. IV B, the universal scaling form of the distribution is preserved during prescaling, allowing us to describe the full evolution solely in terms of time-dependent scaling exponents.
To properly account for time-dependent exponents α(t) and β(t), we generalize the power-law expressions in Eq. (21) as (t/t ref ) α → exp t t ref α(t ) dt /t , and similarly for β [45]. Note that the original expressions are recovered for constant α and β. This way, the evolution of the initial characteristic momentum scale for turbulence p i is given bȳ where t i denotes the time when turbulence sets in. While the exponent β determines the scaling of the characteristic momentum, the evolution of its occupancy is governed by the exponent α, By virtue of Eq. (33), it is in principle possible to compute the time t f by solving f (t f , p f ) = 1, given full knowledge of the time dependence of the exponent α. It is instructive, however, to rewrite this formula using some simplifying approximations. To this end, we assume that α is approximately constant during the driven regime, taking its universal value α driven , and that, after switching off the modulation at time t d , the relaxation to its universal value for free turbulence α free occurs fast compared to the overall reheating time scale. We then recover power-law scaling for the occupancy of the final characteristic momentum, and it reaches unity at the time Recall that a direct energy cascade is described by negative exponents α and β with large initial occupancies f (t i , p i ), which decrease in time according to Eq. (34). Furthermore, the condition t i ≤ t d ≤ t f requires such that the system spends a dominant part of its evolution in the regime of free turbulence.
Using an analogous line of arguments, we can estimate the final characteristic momentum scale from Eq. (32) as where the second equality follows after inserting our estimate for the reheating time (35) and substituting the relation (25) between α free and β free imposed by energy conservation. The use of the latter identity makes the estimate (36) of the final characteristic momentum independent of the scaling exponents in the regime of free turbulence, reflecting the fact that p f is closely connected to the total energy in the system [7], which is conserved after switching off the modulation. By contrast, the reheating time, according to Eq. (35), is sensitive to the values of the scaling exponents in both regimes, and, in particular, can be influenced by the non-universal behavior during prescaling. Finally, we can obtain an order-of-magnitude estimate for the reheating temperature T f in our system by identifying the latter with the typical kinetic energy of a particle with momentum p f . For a dispersion relation of non-relativistic particle, this reads When interpreting the estimates (35) to (37) in the cosmological context, it is important to remember that our time variable t denotes the laboratory time and should be transformed back to the cosmic time τ by integrating the relation dt = dτ /a 2 . Similarly, we need to take into account the redshift of momenta, p → p/a, as expressed in the kinetic term of Eq. (6). For the estimate of the reheating temperature (37), this means T f → T f /a 2 (t f ). This back transformation of variables is illustrated further in the subsequent subsection at a specific example.
Although Eqs. (35) to (37) are useful to estimate asymptotic quantities without requiring to simulate the dynamics up to the point where the system thermalizes, the outlined argumentation is based on the strong assumption that the neglected final stage of the dynamics, which is dominated by quantum fluctuations, does not have a significant impact on these estimates. The latter remains to be checked against physical reality, which can be provided by comparison to an experiment.

E. Thermalization versus expansion
In an expanding spacetime, interaction rates decrease with the expansion. If the expansion is too rapid, it is therefore possible that the system is unable to ever reach  and ν = 2/3. The numerical data are the same as in Fig. 5, but the laboratory time t has been transformed back to the cosmic time τ by virtue of Eq. (40). As shown in the inset, the cosmic time diverges at the finite laboratory time t(τ = ∞) ≈ 444 × t0, which is, in particular, shorter than the reheating time of the associated simulating system. As a result, the evolution of the simulated expanding system slows down and eventually freezes before thermalization is complete. thermal equilibrium. In a cosmological setting, this ability is usually quantified by comparing the typical interaction rate Γ with the Hubble parameter H =ȧ/a. If Γ H, thermalization is possible, whereas for Γ H, the occupation numbers are expected to freeze without thermalizing as the mean free path of the particles exceeds the horizon size of the universe.
The full interpolation between both limiting cases, and the associated thermalization behavior, can be studied in our proposed setup. According to Eq. (7), for d = 2, the dynamics of a Bose gas in an expanding spacetime can be described in terms of a constant interaction g. Consequently, the simulating experimental system will always thermalize eventually. However, this does not necessarily imply the same for the simulated expanding system, whose description is recovered only after a back transformation to the original spatial and temporal coordinates as well as field variables. A necessary condition for thermalization to occur in the simulated system is that the laboratory time t at which the cosmic time τ becomes infinite is larger than the reheating time, t f < t(τ = ∞).
An interesting feature of our transformation of the time coordinate is that in the case of a sufficiently rapid expansion, infinite cosmic time corresponds to a finite value of the laboratory time, i.e., t(τ = ∞) < ∞. To demonstrate this, we consider a power-law expansion [69], with a positive exponent ν. This expression fixes a(0) = 1 and the initial expansion rate H(0) = H 0 . The laboratory time is then given by If ν > 1/2, such as in the case of expansion in a matterdominated universe (ν = 2/3) [48], the limit a = ∞ and τ = ∞ is reached within a finite laboratory time, given by in line with the intuition that a larger expansion rate leads to a faster freeze-out. Remarkably, Eq. (40) associates the state of the simulating system at each instant of the laboratory time t with the asymptotic infinite-time state of a simulated system characterized by a certain value of the initial expansion rate H 0 . This is illustrated in Fig. 6, where we have used The ability to associate the dynamics of a single simulation with arbitrary expansion schemes is a special case of the 2D geometry, where, owing to scale invariance, the equations of motion become independent of the scale factor a. The situation is different in three spatial dimensions, where, according to Eq. (8), the coupling becomes time-dependent as g ∝ a −1 . In the presence of a decreasing coupling, freeze-out may already occur during the simulation before reaching t(τ = ∞), in which case the expanding system will never thermalize.

V. EXPERIMENTAL PERSPECTIVES
In this section, we summarize the experimental requirements for observing the salient features of reheating dynamics discussed in this paper.
The simulating system is a single-component BEC that is parametrically excited by modulating the interaction according to Eq. (11) around some positive offset value. This can be realized experimentally with the help of a Feshbach resonance, which allows one to tune the s-wave scattering length using an external magnetic field [22,33,35,38,39]. The interaction strength g 3D of a 3D Bose gas is related to the s-wave scattering length a s via g 3D = 4π 2 a s /m. Our general approach is independent of dimensionality (but see our remarks towards the end of this section about the absence of a turbulent cascade in one dimension). For concreteness, we have focused here on 2D geometries, which may be realized through a tight confinement of the atomic gas in the vertical direction by a harmonic potential. In the quasi-2D regime, the effective interaction strength g 2D =g 2 /m is characterized by the dimensionless parameterg = √ 8πa s /a HO , where a HO = /mω HO is the oscillator length of the confining harmonic potential with frequency ω HO [57]. Besides via Feshbach resonances, 2D Bose gases with tunable effective interactions can therefore alternatively be realized by changing the frequency of the harmonic trap [70].
The numerical simulations presented in this paper have been conducted for the fixed value Ng = 2.5 × 10 3 of the product of the particle number N and the interaction strengthg. This value is readily achievable in presentday experiments with uniform 2D Bose gases [26]. On the classical level, the parameters N andg enter the equations of motion only via the product Ng (cf. Appendix C 2), but this no longer holds on the quantum level. In fact, our classical-statistical simulations are sensitive to the individual values of N andg since they control the relative magnitude of mode occupancies with respect to quantum fluctuations mimicked in form of vacuum noise. As discussed in Appendix C 3, classicalstatistical simulations are restricted to large occupancies and weak interactions, which is why our simulations of turbulent reheating dynamics have been conducted for a higher particle number than currently realizable in experiments. To estimate characteristic quantities for more realistic experimental setups, we rescale our numerical results to the reference parameters N ref = 1 × 10 6 and g ref = 2.5 × 10 −3 . Pure BECs with almost 10 6 atoms can be reached in state-of-the-art experimental setups, e.g., for 39 K [71][72][73], where a wide range of interactions is accessible through the broad Feshbach resonance at the magnetic flux B = 560.7 G [74]. We emphasize that extrapolating weak coupling results to stronger couplings is a priori not justified. Yet, this procedure is often the only way to be consistent with experimental aspects and commonly used, e.g., in the context of quantum chromodynamics (QCD) [75]. In the field of ultracold atoms, a number of positive examples exist, where weak coupling expectations of non-equilibrium quantum dynamics have been found experimentally at strong couplings [11,76].
For the above choice Ng = 2.5 × 10 3 and a 2D system of 39 K atoms in a square box of length L = 50 µm, one obtains the typical values x 0 = 1 µm and t 0 ≈ 0.6 ms for the characteristic length and time scales x 0 = / √ mn 0 g 0 and t 0 = /n 0 g 0 , respectively. Here, g 0 is the offset value around which the interaction is modulated and n 0 = N/L 2 is the atomic density.
As shown in the previous sections, the reheating dynamics studied in this work can be subdivided into two distinct stages, namely, the early stage of preheating, where occupancies of resonant momenta grow exponentially as a result of parametric instabilities, and the subsequent stage of turbulent thermalization, characterized by a self-similar transport of energy towards higher mo-menta. The time scales at which these phenomena can be observed depend on the choice of the modulation frequency ω and the relative modulation amplitude r in Eq. (11). The choice r = 1, i.e., modulating the scattering length with an amplitude close to its offset value, is efficient for rapidly driving the system out of equilibrium. (Note that the perturbative expression for the growth rate Eq. (17) is no longer applicable in this case.) The direct cascade is observed best if occupancies of the initial characteristic momentum, where energy is injected in the system, are as high as possible. To this end, the modulation frequency ω should be chosen close to the frequency of the lowest non-zero momentum mode supported by the system, which is typically deep in the phonon regime where the Bogoliubov dispersion is linear. In our simulations of reheating dynamics presented in Figs. 4 and 5, we have chosen the modulation frequency such that the primary resonance occurs at the momentum p res = 3 × 2π /L, i.e., ω = 2 pres / ≈ 2π × 199 Hz for the example parameters mentioned above. The onset of turbulence, marking the end of the preheating stage, occurs at around t ≈ 80 × t 0 , corresponding to roughly 10 oscillations of the scattering length. The subsequent turbulent scaling regime extends up to times of about t f ≈ 1500 × t 0 , which can be deduced from Figs. 4 and 5 with the help of Eq. (35) after rescaling the occupancies to the reference particle number N = 1 × 10 6 (cf. Appendix C 3). The estimated reheating time t f is well within the reach of modern experiments, where typical lifetimes of BECs are on the order of a few seconds [72].
In our analysis, we have considered the momentum distribution f (t, p) = Ψ † p (t)Ψ p (t) as the main observable. Experimentally, this quantity can be accessed in time-of-flight (TOF) measurements [57]. In the case of driven turbulence [cf. Fig. 4(a)], by virtue of Eq. (36), we find that the characteristic momentum at the reheating time t f is roughly given by p f ≈ 1.5 × p ξ , which is on the order of the momentum corresponding to the healing length ξ. In general, this value depends on the duration of the regime of driven turbulence and tends to be smaller in the case of free turbulence [cf. Fig. 4(b)] or in the prescaling scenario [cf. Fig. 5], where less energy is injected into the system. It is a fortunate circumstance that higher momenta, which are increasingly populated as the direct cascade progresses, are typically easier to resolve in TOF measurements since they require a shorter expansion time. As an alternative means of detection, in situ images of the density profile may be used. We have checked that all relevant signatures of both the parametric resonance and the turbulent cascade can be extracted from the quantity n † p (t)n p (t) , wheren p (t) is the Fourier transform of the operatorn(t, x) =Ψ † (t, x)Ψ(t, x) describing the spatial density profile of the Bose gas.
Although in this paper we have put our main focus on uniform systems, nowadays available in many laboratories [70,77], we expect the characteristic features of the dynamics to be prevalent even in the presence of a harmonic trap. In general, such an external trapping poten-tial couples different momenta and thereby changes the nature of Bogoliubov quasi-particles. However, within a local density approximation, a harmonically trapped BEC in the Thomas-Fermi regime can be regarded as locally uniform around the center [49]. Thus, while inhomogeneities of the trap may distort long-range correlations in the system, the dynamics on small scales or at large momenta, relevant for the direct cascade, remains unaffected. We have confirmed this explicitly by means of GPE simulations of a harmonically trapped system (not shown).
We conclude this section with some remarks on dimensionality. The observation of a turbulent cascade in a scalar BEC requires at least two spatial dimensions. This is owed to the fact that in a strictly 1D system, due to kinematic constraints, there can be no 2 ↔ 2 scattering processes redistributing momenta. If, however, the system is elongated, but not as strongly confined in the transversal direction to be considered in the quasi-1D regime, turbulent energy transport sets in once transversal modes are excited. We have confirmed this by numerically simulating a 3D system in an elongated periodic box.
Compared to the 3D case, the study of reheating dynamics in two dimensions is somewhat simpler for at least two reasons. First, as discussed in Sec. II B, the effective interaction in 2D is independent of the cosmic scale factor and there is no need to continuously adjust the scattering length according to the specific expansion model chosen. Second, an absorption image of an atomic cloud, taken after a TOF expansion or in situ, will always be a 2D projection on the plane transversal to the optical axis. In three dimensions, this means that the distribution is integrated along the optical axis, thereby mixing momenta of different magnitudes. We have mimicked this in a simulation of a 3D system in a periodic box, where we checked in particular that the scaling is robust if we integrate the momentum distribution along one spatial dimension before performing the radial average. Thus, our scheme can be readily applied to experimentally investigate reheating dynamics in three dimensions, which is of fundamental interest from a cosmological point of view.

VI. CONCLUSION
In this work, we have demonstrated by means of classical-statistical simulations how single-component ultracold Bose gases with tunable self-interactions can be employed for simulating the characteristic stages of the dynamics of post-inflationary reheating. We have shown how expansion can be encoded in the time dependence of the interaction, and that in the special case of 2D, arbitrary expanding backgrounds can be associated with a single experiment in post-processing. The resonant production of particles in the preheating stage of the early universe, driven by the oscillating inflaton, is mimicked by a modulation of the scattering length inducing parametric instabilities in the BEC.
Following the proliferation of secondary instabilities due to non-linear effects, the system enters a turbulent state, characterized by a self-similar transport of energy towards higher momenta. We have analyzed the universal scaling exponents describing this direct cascade, which in general undergoes a transition from a driven to a free regime after the inflaton has decayed. This transition is characterized by a prescaling stage, where the momentum distribution has already attained its universal scaling form, but the scaling exponents are non-universal and time-dependent, yet fully describe the dynamics of the Bose gas. Our set-up opens the door to the experimental observation of both the universal direct cascade and the prescaling regime.
The approach to a non-thermal fixed point as well as prescaling represent a dramatic reduction of sensitivity to the initial state as well as a gradual loss of complexity in the system's dynamics, leading to the emergence of universal physics. Thus, from a more general point of view, our set-up provides a distinct platform for exploring the fundamental mechanisms of thermalization dynamics in quantum many-body systems far from equilibrium.
The final stage of the dynamics, when the system relaxes to thermal equilibrium, is dominated by quantum fluctuations and therefore not captured by the classicalstatistical approximation. In contrast, it can be accessed by experimental studies in cold atomic systems, which are quantum-mechanical by nature and not restricted to the weak-coupling regime. A.C. and K.T.G. contributed equally to this work.

Appendix A: Non-relativistic limit of the Klein-Gordon equation in expanding spacetime
In this appendix, we derive the non-relativistic limit of the equations of motion of the inflaton field (4) on the quantum level and discuss in detail the relevant approximations that enter the derivation.
Canonical quantization can be performed, as usually in quantum field theory on curved spacetime, by transforming to conformal time, dη = dτ /a, and conformal field variables, Φ = a (d−1)/2 φ, [51]. The equations of motion then take the same form as in Minkowski spacetime, with a time-dependent coupling and effective mass Here, (·) denotes the derivate with respect to conformal time η. The field Φ and its conjugate momentum Π = Φ /c are then promoted to operators that satisfy the canonical commutation relations Φ (η, x),Π(η, y) = i δ(x − y),
To this end, we consider the complex field operatorψ(τ, x) defined in Eq. (5), where we have factored out the rapid oscillations of the inflation field on the frequency scale determined by its rest mass. For taking the non-relativistic limit, we require the field to change much slower than these oscillations, i.e., we demand the mass term to be the dominant contribution on the right-hand side of the equations of motion (4). This requires the following conditions to hold. First, the typical field values should be not too large, λ φ2 m 2 c 2 / 2 , which implies that the self-interactions are relatively weak. Similarly, the typical physical momenta should be non-relativistic, |p| mc. Finally, the expansion should not be too rapid, H =ȧ/a mc 2 / , since otherwise the field would be in the over-damped regime and not perform oscillations at all, as during inflation. Note that the first two assumptions are not necessarily satisfied by the inflaton field, especially at the early stages of preheating (see Sec. III).
The outlined assumptions imply the hierarchy where · · · denotes the expectation value. These inequalities are to be understood to hold for any correlation function involving the field or its derivatives. Inserting the ansatz (5) into Eq. (4) without any approximation yields ψ 3 e −2imc 2 τ / + 3ψψ †ψ e −imc 2 τ / + h.c. = 0.
Employing the inequalities (A4) allows us to neglect two of the first four terms in Eq. (A5). Furthermore, we can use relations (A4) to neglect the rapidly oscillating terms proportional toψ 3 e ±2imc 2 τ / in the last bracket. The latter describe number-changing processes in the relativistic theory and neglecting them leads to an emergent U (1) symmetry implying particle number conservation. As a result, Eq. (A5) reduces to Eq. (6), where we have normal-ordered the operators, dropping an irrelevant energy shift 5 .
By virtue of the relations (A4), it is straightforward to show that within the outlined approximations the canonical commutation relations (A3) imply bosonic equal-time commutation relations for the fieldΨ =ψa d/2 , (A6)

Appendix B: Formal differences of parametric oscillations in the relativistic and non-relativistic model
To pinpoint the formal differences between the parametric oscillations described by Eq. (10) and those induced in the Bose gas by periodically modulating the interaction, we consider Eq. (A5), which is an exact rewriting of the original relativistic equation of motion (4). Linearizing the field operator around the homogeneous condensate,ψ(τ, x) = ψ 0 (τ ) + δψ(τ, x), the interaction term takes the form which corresponds to the term proportional to λφ 2 0 δφ in Eq. (10). When taking the non-relativistic limit to derive Eq. (7), we have neglected the rapidly oscillating U (1)violating term proportional to e −2imc 2 τ / , keeping only the U (1)-conserving, slowly-varying part. However, since the former reflects the oscillations of the inflaton responsible for parametric resonance, this effect is no longer present in the non-relativistic model.
To re-introduce parametric instabilities in the nonrelativistic model, we add a periodic modulation of the coupling with frequency ω, such that the linearized interaction term in Eq. (7) takes the form On a formal level, this is similar to Eq. (B1) in the sense that we have an oscillating prefactor multiplying a term cubic in the field, although it does not have exactly the same form and therefore the mapping is not exact.
complex Gaussian random variables with zero mean, satisfying α * p α q = δ p,q /2 at zero temperature. Here, (. . . ) denotes the ensemble average over all realizations. This vacuum noise corresponds to an average occupancy of half a particle per mode, which is also known as the "quantum half". Unless stated otherwise, the vacuum noise 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 ordinary plane waves. Effectively, this corresponds to preparing the Bogoliubov ground state of an ideal gas, which is appropriate for our application since the precise nature of the quasi-particles seeding the parametric resonance is unimportant. As a side effect, we observe a transient growth of population at low momenta in our simulations at early times (cf. Fig. 2). The resulting p −1 behavior of the momentum distribution at low momenta matches the behavior of u 2 p and v 2 p for |p| → 0. Therefore, the observed growth can be interpreted as an artifact of the effective quench from an ideal Bose gas to a system with finite interaction at time t = 0.
In the truncated Wigner method, quantum expectation values are replaced by statistical averages over the Wigner distribution. The latter correspond to expectation values of symmetrically ordered quantum operators. Thus, the momentum distribution obtained from classical-statistical simulations corresponds to the one defined in Eq. (12) plus an extra contribution in form of the "quantum half" stemming from the commutation relations, Here, Ψ cl (t, p) = d d x Ψ cl (t, x)e −ipx/ denotes the Fourier transform of the classical field Ψ cl (t, x).
The classical-statistical simulations presented in this work have typically been averaged over at least 64 runs. The statistics is even further enhanced through radial averages due to the isotropy of the system. The statistical error bars are thus typically smaller than the line width in the plots.

Dimensionless Gross-Pitaevskii equation
Each realization is propagated in time according to the GPE (13). To cast this equation into a dimensionless form suitable for numerical simulations, we express time and length in units of the characteristic scales t 0 = /n 0 g 0 and x 0 = / √ mn 0 g 0 , respectively. In terms of the dimensionless timet = t/t 0 , positionx = x/x 0 , and fieldΨ cl = x d/2 0 Ψ cl , the GPE takes the form with the dimensionless couplingg = g × t 0 / x d 0 .
In a quasi-2D system (and in the absence of timedependent modulations), the dimensionless coupling is, up to logarithmic corrections, given byg = where a s is the s-wave scattering length and a HO is the oscillator length of the confining harmonic potential in vertical direction [57]. If the wave function is normalized to unity, the coupling changes asg → Ng, where N is the particle number. This shows that the parameters N andg enter the classical equations of motion (C3) only through the product Ng as the single relevant model parameter. Moreover, also the dimensionless box length L/x 0 = √ Ng is fixed by this quantity as a consequence of scale invariance. By contrast, on the quantum level, the relative magnitudes of the parameters N andg play a role (cf. Appendix C 3).
The GPE (C3) is discretized in space by means of a Fourier pseudospectral discretization and propagated in time using the well-known split-step method [79]. The accuracy of the time evolution is further enhanced with the help of the optimized fourth-order time splitting scheme given in Table 2 of Ref. [80].
In our 2D simulations of turbulent reheating dynamics (cf. Sec. IV), we use a spatial N g × N g grid with at least N g = 512 grid points per dimension, while N g = 256 turns out to be sufficient for simulating early-time preheating dynamics (cf. Sec. III). For a system of size L × L with L/x 0 = 50, the corresponding grid spacing ∆x = L/N g ensures that the healing length ξ = x 0 / √ 2, which is the smallest physical length scale in the system, is well-resolved. Numerical stability for late-time dynamics is achieved by choosing the numerical time step ∆t such that ∆tk 2 max 2π, wherek max = π/∆x with ∆x = ∆x/x 0 is the maximum numerical wave number supported by the grid [81].

Validity of classical-statistical simulations
Quantum dynamics is well described by classicalstatistical simulations in the regime of large occupancies and weak couplings [28,47]. In particular, the truncated Wigner method is known to produce unphysical results if the number of virtual particles added to mimic quantum fluctuations becomes comparable to the number of real particles in the system [29,30]. At zero temperature, the failure of the classical-statistical approximation is indicated by the instability of the vacuum and a resulting unphysical dependence of observables on the UV cutoff [47]. This decay of the "quantum half", which inevitably occurs at sufficiently late times, naturally restricts the classical-statistical method to the weak coupling regime, where the instability is mitigated via a separation of scales. Moreover, if the coupling is too strong, physical observables risk attaining a dependence on the UV cutoff already at relatively early times through the spurious quantum pressure generated by the decaying vacuum. This is because the coupling controls the relative magnitude of mode occupancies with respect to the vacuum noise, which, in turn, is regulated by the UV cutoff [47].
To ensure that our numerical simulations of reheating dynamics, in particular, in the turbulent scaling regime at late times, lie within the range of validity of the classical-statistical method, we follow a similar procedure as presented in Ref. [47]. Formally, the classical statistical limit is given by N → ∞ and g → 0, where N is the particle number and g is the coupling, such that N g = const. To approach this limit in our 2D simulation, we fix the value Ng = (L/x 0 ) 2 by the choice L/x 0 = 50 and repeat the simulation for different values of the couplingg, using the parametriza- Figure 7 shows the radially averaged momentum distribution for the driven turbulence scenario [cf. Fig. 4 for the distribution with the lowest value ofg, around which the UV scaling properties of the distribution are most sensitive. It can be seen that deviations occur for g/g ref 0.1. For larger values of the coupling, the shape of the power law changes and the distribution hits the UV cutoff, such that cutoff-independence of the results is no longer guaranteed. Furthermore, the distribution at large momenta drops below the "quantum half", as indicated by the horizontal dotted lines, which formally corresponds to the unphysical situation of negative occupancies. We have also checked, using a similar approach as in Ref. [47], that the contribution of the vacuum noise to the total energy budget of the system becomes significant forg/g ref 0.1, which coincides with the values of the coupling for which deviations in the momentum distribution are observed.
While the early-time preheating dynamics presented in Sec. III is only insignificantly affected by these deficiencies for the reference parameters N ref = 1 × 10 6 and g ref = 2.5 × 10 −3 , this is no longer true for the late-time turbulent reheating dynamics, as illustrated in Fig. 7. The results presented in Sec. IV have therefore been computed forg/g ref = 0.01, which safely lies within the range of validity of classical-statistical simulations.
resonance. The width of this instability band is delimited by the modes satisfying A p = 1±q p , which increases with the amplitude r of the modulation.
The growth rate predicted from Eq. (D2) is compared to our numerical simulations in Fig. 3.

Appendix E: Secondary instabilities
In this appendix, we discuss the secondary amplification of fluctuations due to the non-linear corrections in Eq. (19). For modes that are stable on the linear level, the integral on the right-hand side of Eq. (19) acts as a source term. It has a complicated momentum structure, but its time evolution is straightforward since it is dominated by the contribution from the exponentially growing unstable modes q, for which |p − q| ≈ |q| ≈ p res . This condition can be satisfied if p 2p res , as can be seen from the triangle inequality, For p 2p res , on the other hand, the growth of the integral is exponentially suppressed. The narrower the primary resonance, the more rapid the integral drops at p ≈ 2p res .
As in the previous appendix, we assume here that the resonance is not too strong, such that ζ pres ω. The time evolution of resonant modes can then be written as a product of an exponential growth and an oscillating function [59,60], n 1 p ∝ e ζp res t e ±iωrest , θ 1 p ∝ e ζp res t e ±iωrest , with the growth rate ζ pres given by Eq. (D2) and ω res being the frequency of the resonant modes. As a result, Eq. (19) takes the form of a forced harmonic oscillator, n 1 p + ω 2 p n 1 p = e 2ζp res t S p e 2iωrest + Q p e −2iωrest + P p . (E1) Here, Q p and S p contain the momentum-dependent contributions of the integrals in Eq. (19) corresponding to positive and negative frequencies, respectively, and P p contains the non-oscillating part. The contribution to the amplitude of n 1 p from the first term on the righthand side of Eq. (E1) proportional to S p is given by and analogous expressions can be obtained for Q p and P p , respectively. In other words, the modes in the abovementioned momentum range are amplified as exp(2ζ pres t) and a resonant amplification occurs if ω p ≈ 2ω res . A smoother spectrum is expected in the case of strong driving, i.e., large ζ pres .
The pattern described above is indeed observed in Fig. 2, where the secondary and higher-order peaks at multiple frequencies are clearly visible.

Appendix F: Analysis of instabilities in terms of quantum equations of motion for correlation functions
In this appendix, we study the instabilities arising in the presence of a periodically modulated coupling using quantum equations of motion for the two-point correlation functions. We validate the results of Sec. III, where quantum fluctuations have been treated within the classical-statistical approximation.
The classical action for a (non-relativistic) bosonic field Ψ with a quartic self-interaction is given by [82] (F1) The integration is performed over the Schwinger-Keldysh closed-time contour [83]. In the following discussion, we treat ϕ 1 = √ 2 Re Ψ and ϕ 2 = √ 2 Im Ψ as the two independent field components.
In quantum theory, the central objects are field correlation functions. The one-point correlation function is defined as φ i (x) = φ i (x) . Throughout this appendix, we use the notation x = (x 0 , x) ≡ (t, x). For homogeneous systems, the one-point function depends only on time and we can identify it with the wave function of the condensate. The time-ordered connected two-point function is defined as [82,84] where T is the time-ordering operator along the closedtime contour and sgn C (x 0 − y 0 ) takes the values ±1 depending on whether x 0 is after or before y 0 on the Schwinger-Keldysh contour. In the last step, the propagator has been decomposed into the spectral function ρ and the statistical propagator F , defined as Here, [·, ·] and {·, ·} denote the commutator and the anti-commutator, respectively. The spectral function describes the spectrum of the excitations and, in particular, encodes the equal-time commutation relations, while the statistical propagator characterizes its occupations.
The equations of motion for the one-and two-point functions can be obtained by a variational principle from the two-particle irreducible (2PI) effective action. For the two-point functions, they have the form [84] and Σ F and Σ ρ are the components of proper self-energy , y), which sums all non-local two-point one-particleirreducible diagrams constructed from the interaction vertices in the action (F1) and with lines associated to full time-ordered propagators G. Summation over repeated indices is implied in all expressions in this appendix. The presence of a non-vanishing one-point function generates effective cubic interactions in addition to the quartic one, which follows from the shift The resulting interaction part of the action, in terms of the fields ϕ i , reads Equations (F5) are supplemented with analogous equations of motion for the one-point function φ. For the discussion of instabilities, the back-reaction of produced excitations on φ is negligible, and, as a result, its evolution is approximately governed by the classical GPE.

Primary instabilities
In the linear regime, the self-energies in Eq. (F5) as well the F -dependent corrections to the mass in (F6) can be neglected. We consider the evolution of the statistical propagator F , which can be written as It is convenient to introduce the "normal" and the "anomalous" correlation functions, defined as [82] n(x, y)  which are complex-valued functions. Assuming homogeneity, these functions depend only on the relative spatial coordinate, i.e., n = n(t, t , x − y). The evolution equations (F8) for the components of the statistical propagator, Fourier transformed with respect to this relative coordinate, are given by These equations have the same form as the ones for the linearized (classical-statistical) fluctuation field δΨ p (t). In particular, the linear combinations ψ * (t) n p (t, t ) + ) satisfy the wave equations (16). Therefore, ignoring the time dependence of the coupling, one has where ω p corresponds to the Bogoliubov dispersion relation. The dependence on t in the above expression is determined from the symmetry properties of the functions n p (t, t ) and m p (t, t ), given in (F9).
The unstable modes exhibit an exponential amplification, n p (t, t ) ∝ m p (t, t ) ∝ exp(ζ p (t + t )), as explained in section (III D). Here we again employed the symmetry properties of n and m. The same holds for the components of the spectral function. However, at t = t the spectral function is determined solely by the commutation relations, as follows from its definition (F4), In the next subsection we include the first non-linear corrections to describe the secondary amplification of fluctuations.

Secondary instabilities
As the fluctuations become sufficiently strong, the linear approximation for the dynamics becomes inapplicable. The non-local one-loop correction, which is diagrammatically illustrated in the right part of Fig. 8, accounts for secondary instabilities [19,40]. The corresponding self-energy, which can be calculated by means of standard Feynman diagram techniques, reads The other one-loop correction is the tadpole diagram, shown in the left part of Fig. 8, which is explicitly included in the definition (F6) of the effective mass . As can be seen from that definition, this correction becomes relevant only later, when the components of F ij become of the order of the condensate density |Ψ| 2 (see also [40]). Including the one-loop self-energy in Eq. (F5) leads to integrals over the past evolution, i.e., memory integrals, on the right-hand side. Due to the exponential growth of fluctuations, these integrals are, however, dominated by the latest times [19,84], which allows to simplify them. Note that Σ contains terms proportional to The Σ ρ F term in the limit of z 0 → x 0 represents a "tadpole"-like contribution, which becomes important only at later times [40]. For the Σ F ρ term one obtains where c −1 = O(g 0 n 0 ) is some parameter of the order of the oscillation frequency of the condensate. In Σ F , there will be terms ∝ F 2 and ∝ ρ 2 . We first consider the ∝ F 2 part in Σ F , which is equivalent to Σ| G→F . Inserting Eq. (F12) into the equations of motion (F5) and using the transformation (F9) and (F10), one arrives at the following equations for the momentum modes of the normal and the anomalous two-point functions, − 2g(t)|Ψ(t)| 2 n p − g(t)Ψ 2 (t) m * p = icg(t)g(t ) q 4Ψ(t)Ψ * (t ) n q n * p−q + m q m * p−q + 4Ψ * (t)Ψ * (t ) n q m p−q + 4Ψ(t)Ψ(t ) n q m * p−q + 2Ψ * (t)Ψ(t ) n q n p−q , i ∂ t − p 2 2m − 2g(t)|Ψ(t)| 2 m p − g(t)Ψ 2 (t) n * p = −icg(t)g(t ) q 4Ψ(t)Ψ(t ) n q n * p−q + m q m * p−q + 4Ψ * (t)Ψ(t ) m q n p−q + 4Ψ(t)Ψ * (t ) m q n * p−q + 2Ψ * (t)Ψ * (t ) m q m p−q , (F14) where we have omitted the arguments t and t from the two-point functions.
The ∝ ρ 2 part in Σ F , which is equivalent to Σ| G→iρ/2 , leads to an analogous expression on the right-hand side. However, the components ρ are of order one in the t → t limit (in momentum space) due to Eq. (F11). As a result, that contribution is negligible compared to the exponentially growing one from the ∝ F 2 part. This also justifies the use of classical-statistical simulations for studying the evolution of equal-time observables.
The momentum integrals on the right-hand side are analogous to the ones in Eq. (19) and therefore lead to a secondary exponential amplification, i.e., n p (t, t ) ∝ exp(2ζ p (t + t )), in the momentum range p 2p res . One can check that the terms on the right-hand side os-cillate 6 as functions of t either as exp(−ig 0 n 0 t) or as exp(−i(g 0 n 0 ± 2ω res )t). Hence, one expects enhanced amplification near ω p = 0 and ω p = 2ω res . The resonant frequencies correspond to the conservation of energy (in addition to the conservation of momentum) in the vertices of the Feynman diagrams.
To conclude, in this appendix we have analyzed the primary and secondary instabilities within the nonequilibrium 2PI framework and explained the validity of the classical-statistical approximation, justifying its use in Sec. III to describe the regime when typical occupation numbers of excitations become sufficiently large. extracted from the moments of orders 1 ≤ n1 < n2 ≤ 4. The vertical dotted line marks the time t ≈ 458.6 × t0 when the modulation is switched off instantaneously. Before this point, the oscillatory behavior of the moments is directly reflected in the evolution of the scaling functions. As driven turbulence develops, their time averages approach power laws with exponents close to the predictions from kinetic theory, α driven = −1 and β driven = −1/2 (dashed lines). After the modulation is switched off, the scaling functions exhibit a kink and the oscillations vanish. In the subsequent evolution, the scaling functions extracted from different moments evolve asynchronously, until they adopt a power-law behavior again for times t 1000 × t0 with exponents close to the predictions from kinetic theory in the regime of free turbulence, α free = −1 and β free = −1/4 (dashed-dotted lines).
The system evolves self-similarly where all curves have the same slope, as analyzed in Fig. 10. ments outlined above for all possible combinations of moments 1 ≤ n 1 < n 2 ≤ 4. We exclude the zeroth moment from our analysis since the associated particle number conservation is incompatible with the energy-conserving direct cascade. All simulation parameters are the same as in Fig. 5, with the exception that the modulation is switched off at a later time, t = 56×2π/ω ≈ 458.6×t 0 , allowing one to better distinguish the regimes of driven and free turbulence. The classical-statistical simulation presented in Figs. 9 and 10 has been conducted on a spatial grid of 768 × 768 grid points to increase the UV resolution and the data have been averaged over 100 runs. To reduce numerical errors of the calculated moments due to instabilities near the UV cutoff, we restrict the integral in Eq. (H1) to absolute momenta |p| ≤ 3×p ξ (|p| ≤ π /∆x, where ∆x is the lattice spacing) before (after) switching off the modulation. As explained in Appendix G, we shift the origin of time as t → t−63×t 0 to improve the scaling analysis at early times. Figure 9 shows the scaling functions s α (t) and s β (t) with respect to the reference time t ref = 85 × t 0 . As turbulence develops, the scaling functions approach power laws s α (t) ∝ (t/t ref ) α and s β (t) ∝ (t/t ref ) β . In the driven regime, before switching off the modulation, the corresponding exponents are close to the analytical predictions from kinetic theory, α driven = −1 and β driven = −1/2 (see Sec. IV C). At the time when the modulation is switched off instantaneously, the scaling functions exhibit a kink and continue evolving asynchronously for some time. This indicates that the shape of the distribution is not preserved exactly, but slightly readjusts during the transition from driven to free turbulence. The changes in the shape of the distribution are more pronounced at lower momenta, which is reflected by the fact that those curves in Fig. 9 extracted from lower moments react stronger at the transition. Despite these local adjustments, it is remarkable that the global scaling form of the distribution is approximately preserved during the transition from driven to free turbulence, as can be seen in Fig. 5. For t 1000 × t 0 , the scaling functions again approach a power law with exponents close to the analytical predictions for free turbulence, α free = −1 and β free = −1/4. In the region where the scaling functions extracted from different combinations of moments all evolve in parallel, the distribution scales self-similarly as a whole. To certify that the distribution also exhibits prescaling, i.e., a self-similar evolution with time-dependent scaling exponents, we analyze the change of the scaling functions as quantified by the exponents α(t) and β(t) in what follows.
In the driven regime, the scaling functions oscillate in time, as indicated in the inset in the lower panel of Fig. 9. These oscillations originate from the modulated interaction, which causes the distribution and therefore the moments to oscillate in time (cf. Fig. 3). As a result, the instantaneous scaling exponents α(t) and β(t), which are derivatives of the scaling functions according to Eq. (H6), exhibit the strong oscillatory behavior shown for β(t) in the inset in the upper panel of Fig. 10. To connect with our previous results obtained using the maximum likelihood technique (cf Fig. 5), it is convenient to consider instead the time-averaged exponentsᾱ(t) andβ(t), defined asβ and analogously forᾱ(t), where T = 2π/ω is the modulation period. After switching off the modulation, the oscillations of the moments vanish and the time-averaging of 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 derivatives in Eq. (H6) are sensitive to small fluctuations in the data due to statistical uncertainties, resulting in a large spread of the data points. This is especially true where moments of lower orders are involved since they probe the distribution at smaller momenta with a lower density of states. To better visualize the trend of the exponents, the data have been smoothed by calculating simple moving means (separately in the regimes of driven and free turbulence) using a window size of 8 (16) periods involving 512 (1024) data points for n 1 > 1 (n 1 = 1). The shaded regions represent the corresponding moving standard deviations. Before the modulation is switched off, the system is in the state of driven turbulence. Indeed, for t 200 × t 0 , bothβ(t) andᾱ(t)/β(t) are approximately constant and for all considered moments close to the analytical predictions β driven = −1/2 and α driven /β driven = 2 (cf. Sec. IV C). After suddenly switching off the modulation, the exponents jump discontinuously, reflecting the kink in the scaling functions (cf. Fig. 9). Furthermore, there are discrepancies in the values of the exponents extracted from different combinations of moments after the modulation is stopped. As discussed above, this indicates a slight readjustment of the shape of the distribution, which is consistent with the behavior in Fig. 5, where the exponents exhibit large error bars in this regime. For times t 1000 × t 0 , the exponents extracted from different moments converge to a single curve, certifying a self-similar evolution of the distribution as a whole. During this prescaling stage, the exponentβ(t) evolves gradually towards the analytical prediction for free turbulence, β free = −1/4. By contrast, the ratioᾱ(t)/β(t) adjusts rather quickly to the value α free /β free = 4, expressing energy conservation.
In conclusion, our analysis certifies the existence of a prescaling regime during the transition from driven to free turbulence, where the momentum distribution scales self-similarly with time-dependent scaling exponents approaching gradually their universal values.