Scaling approach to quantum non-equilibrium dynamics of many-body systems

Understanding non-equilibrium quantum dynamics of many-body systems is one of the most challenging problems in modern theoretical physics. While numerous approximate and exact solutions exist for systems in equilibrium, examples of non-equilibrium dynamics of many-body systems, which allow reliable theoretical analysis, are few and far between. In this paper we discuss a broad class of time-dependent interacting systems subject to external linear and parabolic potentials, for which the many-body Schr\"{o}dinger equation can be solved using a scaling transformation. We demonstrate that scaling solutions exist for both local and nonlocal interactions and derive appropriate self-consistency equations. We apply this approach to several specific experimentally relevant examples of interacting bosons in one and two dimensions. As an intriguing result we find that weakly and strongly interacting Bose-gases expanding from a parabolic trap can exhibit very similar dynamics.


Introduction
Understanding the time evolution of complex quantum systems, often in the presence of strong correlations between constituent particles, is crucial for solving many fundamental problems in physics, from the expansion of the early universe, through heavy ion collisions, to pump and probe experiments in solids. New questions about dynamical evolution arise in recently realized artificial quantum many-body systems, such as ultracold atoms in optical potentials or photons in media with strong optical nonlinearities. These systems are only weakly coupled to external heat baths and have a limited lifetime; thus many experiments require interpretation in terms of coherent quantum dynamics rather than the properties of equilibrium states. On the positive side, these systems allow remarkable control of parameters and open up exciting opportunities for controlled experiments exploring non-equilibrium many-body dynamics.
In the realm of many-body physics, low-dimensional systems have a special place. They have dramatically enhanced quantum and thermal fluctuations and exhibit most surprising manifestations of strong correlations. Rigorous theorems provide strong constraints on longrange order and often such systems cannot be analyzed using mean-field approaches even at zero temperature. Nevertheless, equilibrium properties are well understood using methods specific to low dimensions, such as Coulomb-gas representation of vortices in two dimensions or effective low-energy descriptions of one-dimensional (1D) systems, including the Luttinger liquid and sine-Gordon models (see e.g. [1]). However, such an analysis cannot be straightforwardly extended to non-equilibrium dynamics. Most equilibrium theories focus on the low-energy part of the spectrum, while non-equilibrium dynamics can couple degrees of freedom at very different energy scales [2]- [15]. It would be highly valuable to have examples of many-body dynamics of low-dimensional strongly correlated systems amenable to an unbiased analytical treatment. These examples could be used not only for analyzing experimental systems but also for testing theoretical calculations utilizing effective models or approximations, and for 3 checking the validity of new numerical approaches. In this paper, we propose such a class of non-equilibrium quantum problems with time-dependent Hamiltonians that allow for a scaling ansatz of many-body wave functions.
Scaling solutions in quantum dynamics were first discussed in the context of a single harmonic oscillator with a time-dependent frequency [16]- [20]. This problem can be reduced to a time-independent one by properly rescaling space and time. Scaling transformation of variables is possible due to the existence of a dynamical symmetry generated by dynamical invariants of the system [18,19]. There are also extensions of this approach to single-particle problems with potentials of the Coulomb and inverse square type [17,21,22]. In the context of many-body problems, scaling has first been used within mean-field approaches to bosonic systems, for the classical Gross-Pitaevskii equation [23,24], [27]- [29]. Beyond these effective one-body problems, scaling solutions exist for hard-core bosons in 1D [30] and in the unitary limit of fermionic gases with infinite scattering length [25]; these are problems for which the interaction enters a constraint on the wave function of an otherwise non-interacting system analysis. Away from these specific limits, Pitaevskii and Rosch [26] introduced a scaling ansatz for a 2D many-body system of particles interacting with contact or inverse square interaction and related the existence of such a solution to a hidden SO(2, 1) symmetry. In this paper, we further extend full many-body scaling solutions to more general types of interaction and arbitrary dimensionality. This generalization can be achieved by allowing three parameters of the system-the mass, the interaction constant and the external potential-to be time dependent. A scaling solution is possible when the interdependence of these parameters is given by an Ermakov-type equation, similar to the one discussed in earlier approaches [25,26,30], and an additional self-consistency equation that depends on the dimensionality of the system and the nature of interactions.
Dynamical control over the system parameters is possible in recently developed artificial quantum systems, such as trapped ultracold atomic gases, where the effective interaction can be tuned either using Feshbach resonances or by changing the transverse confining potential, whereas the effective mass can be changed by the application of the weak optical lattice [31]. It is also possible with photons in nonlinear optical devices, where the timedependent dispersion and Kerr nonlinearity can be achieved using electromagnetically induced transparency [32]- [36]. In this paper, we propose applications of the scaling ansatz that are experimentally relevant in the context of both these systems.
We emphasize that apart from the tunability of the parameters, no specific restrictions on the system properties are imposed. Particles can obey fermionic, bosonic or mixed statistics, interact by pairwise interaction and be subject to parabolic confining potential, to a linear potential and to a complex chemical potential. The basic idea of the scaling solution presented hereafter is to map the non-equilibrium equations of motion to an equilibrium many-body Schrödinger equation. The mapping is based on scaling functions that relate correlation functions of time-dependent systems to correlation functions of systems in equilibrium. Hence, several results known for equilibrium many-body systems can be directly translated to nonequilibrium situations. The reverse conclusion is also true: from a measurement of the system that is out of equilibrium, e.g. a quantum gas after expansion, we can deduce its initial (equilibrium) properties [37].
The paper is organized as follows. In section 2, we introduce a general formalism of scaling transformation for a many-body Schrödinger equation. In section 3, as an example of the application of our approach, we compute momentum distributions for 1D and 2D bosonic 4 gases with contact interactions released from a parabolic trap. Further details are given in the appendices, where we also discuss the relationship of our work to classical integrability of timedependent bosonic systems with contact interactions.

Scaling transformation-general approach
Our starting point is the many-body Schrödinger equation for N interacting particles in D dimensions, The external parameters (chemical potential µ(t), linear potential g(t) and trapping frequency ω(t)) and the many-body interaction potential V (x; t) depend explicitly on time. The chemical potential µ(t) = [µ(t)] + i [µ(t)] can accommodate the effects of dissipation via its imaginary parts 5 . While the dependences on the linear and chemical potentials can be removed by the Gallilei transformations and phase shifts, respectively, we note that solving the quantum problem with time dependence of the remaining parameters represents a non-trivial task. For instance, unlike in the non-interacting case, the time dependence of the mass cannot be removed by a simple redefinition of the time variable.
We address the following question: under what conditions can equation (1) (the -system) be transformed into the Schrödinger equation for a time-independent ( -) system? i ∂ (y 1 , . . . , y N ; τ ) ∂τ = H 0 (y 1 , . . . , y N ; τ ), We emphasize that so far in (2), ω 0 and m 0 are unspecified parameters; in particular, the -system can have vanishing confining potential even when the -system is confined. We assume that the time dependence of the pairwise interaction potential enters through a single We further assume that the interactions have a scaling property and are characterized by the exponent α, which we take to be the same for both -and -systems, Most generic interaction potentials (or pseudo-potentials) satisfy a scaling law (3): s-wave . Other examples 5 are ultracold fermions interacting via a p-wave channel, which gives rise to the δ pseudopotential (α = D − 1). Also logarithmic potentials can be treated; scaling the logarithmic law produces a time-dependent shift to µ(t).
To express the solution of the time-dependent Schrödinger equation (1) in terms of the solution (y 1 , . . . , y N ; τ ) of the static equation (2), we introduce the scaling ansatz, with y i = (x i /L(t)) + S(t) and τ ≡ τ (t). Direct calculation shows (see appendix A) that this ansatz is valid if the scaling functions R(t), L(t), F(t), τ (t), G(t), S(t) and M(t) satisfy a set of coupled differential equations, It is not obvious a priori that equations (5)- (12) can be satisfied simultaneously for any reasonable time dependences of system parameters m(t), v(t) and ω(t). Our next goal is to show that there are a number of non-trivial cases for which equations (5)- (12) are consistent with each other. First of all, we note that equations (5) and (6) imply that In the absence of dissipation ( [µ(t)] = 0), this condition is equivalent to the conservation of the norm of the wave function under the scaling transformation. Equation (6) allows us to express F(t) via L(t), F(t) = (m(t)/2)L/L, which can then be substituted into equation (7). This leads to the differential equation for L(t), where h(t) = m 0ṁ (t)/m(t). The term with the first derivative can be removed by the change of variables L(t) = exp[B(t)]y(t) withḂ(t) = −h/2. For y(t) we obtain where 2 (t) = 1 4 h 2 − 1 2ḣ + ω 2 (t). Equation (14) is the celebrated Ermakov equation [38] first discovered in 1880 [39]. This equation has been used primarily for tracking invariants of the time-dependent harmonic oscillator. In appendix B, we show how one can use the nonlinear superposition principle to reduce equation (14) to the linear equation. Once L(t) is known, the remaining set of equations for S(t), M(t) and G(t) can be solved directly.
In summary, to find time-dependent parameters that admit a scaling solution, one can apply the following recipe: after specifying two time-dependent functions ω(t) and m(t), one obtains a solution of the Ermakov equation (14) from which one determines time-dependent interaction strength v(t) consistent with equation (12). The solutions of the functions M(t), G(t) and S(t) can then be obtained straightforwardly provided that the functions g(t) and µ(t) are explicitly specified. Note that the complexity of our method (e.g. solving the Ermakov equation) does not depend on the number of particles N .
The initial conditions for systems (1) and (2) are related to each other through equation (4) applied at time t = 0, Generally, at t = 0 the Hamiltonians controlling the dynamics of -and -systems do not coincide. For example, they can have different confining potentials, or one system can be in a trap while the other one is in free space (ω 0 = 0). In this paper, we focus on a finite initial trapping potential, ω(0) = ω 0 > 0, for which we introduce the additional assumption that at t = 0 the two systems coincide. This means that we have m( At t > 0, the parameters of the -system begin to change in time while the parameters of the -system remain constant. Since the two systems coincide for t < 0, the initial state of the -systems at t = 0 should correspond to the equilibrium state of the -system. The existence of the scaling solution in 1D in the hard-core limit v 0 → ∞ has been established previously [30]. Within our approach, this can be understood as follows: the first equation of (12) is trivially satisfied, whereas other equations do not depend on the interaction strength and remain valid. Another special case is the 2D system with contact interactions studied previously by Pitaevskii and Rosch [26] (D = −α = 2), for which equation (12) is satisfied by constant mass and interaction.

Dynamics of a Bose gas with contact interaction released from the trap
In this section, as an example, we apply the scaling approach to an ultracold Bose gas with contact interaction that is prepared in a confined, weakly interacting initial state. The non-trivial dynamics comes from a sudden switching off of the confining potential from ω(t) = ω 0 at t = 0 to ω(t) = 0 at t > 0. The solution of the scaling equation (13) for constant mass m(t) = m 0 is then given by L(t) = (1 + ω 2 0 t 2 ) and consequently F(t) = (m 0 ω 2 0 t/2)L 2 (t). In appendix B.3, we examine additional scenarios corresponding to varying mass that exhibit similar behavior 7 of the scaling functions. Here, we also assume that µ(t) and g(t) are time-independent constants.
To characterize the non-equilibrium dynamics, it is convenient to deal with correlation functions that can be easily derived within the scaling approach (appendix D). The dynamics of the momentum distribution, for example, can be related to the single-particle density matrix g 1 of the initial state, From the asymptotic behavior of the scaling functions we can extract the long-time limit of the momentum distribution using the stationary phase approximation (SPA), Hence the momentum distribution becomes fully determined by the density distribution For a quantitative description of dynamics, we need to specify the initial correlation function, which we take from earlier analyses of effective theories for weakly interacting Bose gases in harmonic traps [40]- [42]. An important characteristic of a condensed state with a sufficiently large number of particles is the Thomas-Fermi shape of the density profile, First, we analyze the 1D case in the low-temperature regime when the coherence length is of the order of the Thomas-Fermi radius (equation (E.1) of appendix D). According to the scaling equation (12), for contact interactions, V (x, t) = v(t)δ(x) (α = −D), the interaction must be tuned inversely proportional to the scaling function, v(t) = v 0 /L(t). In figure 1(a), the results of a numerical evaluation of the momentum distributions (15) for specific initial values are shown together with results from SPA. The behavior of the p = 0 component is characterized by a steep decay on a time scale ω −1 0 followed by slowly dephasing oscillations, which are due to the finite extension of the density profile and the quadratic phase factor in (15). The corresponding period of oscillations P is determined by the Thomas-Fermi radius, Oscillations as a function of |p| at constant t can be attributed to the finite Thomas-Fermi radius as well. Here, the quadratic phase factor leads to the oscillation period growing with |p|. In agreement with the SPA prediction, the momentum distribution relaxes to a semi-circle law. This is remarkable, since such a behavior has been previously associated with 1D Bose systems in the strongly interacting limit (v 0 → ∞) [30] only. In our case, the interaction strength is initially small and then even decreases with time. We note that this cannot be understood as the effect of dilution due to expansion of the system because the effective 1D interaction parameter [43] In 2D, α = −2 and equation (12) leads to interactions which are constant in time. When the initial state is weakly interacting (appendix D), we choose an effective theory that incorporates effects of quantum and thermal fluctuations [41]. Results of numerical evaluation of equation (15) are shown in figure 1(b). The momentum distribution evolves very much like in the 1D case and is essentially determined by the initial density distribution and the associated Thomas-Fermi radius. Here, the number of particles (N = 16) is set to be smaller than in - [42], see also appendix D). Dynamical evolution is obtained from numerical integration of equation (15). The SPA represents the asymptotic t → ∞ result. Numerical errors are of the order of the line thickness. In the 1D case (a), the system parameters are N = 140, √h /(m 0 ω 0 ).
the 1D system. Therefore, the asymptotic stationary phase solution is approached slowly and oscillations dominate in the analyzed time window ω 0 t 20. We checked that in both 1D and 2D the results are robust against variation of temperature and interactions as long as phase coherence is not destroyed. The analysis of these examples leads to remarkable consequences. We note that the stationary phase regime is reached rather quickly with momentum distribution determined by the initial density distribution. Therefore specially designed initial density distributions (equilibrium or not) can be used to create specific momentum distributions, such as step-like fermionic ones, on demand. It is remarkable that such behaviour, which has been obtained previously in the strongly interacting limit, persists down to arbitrarily weak strength of interaction. This is opposite to what is realized in time-of-flight experiments of ultracold atoms released from a lattice [44], where the expansion at sufficiently large times can be regarded as free and momentum distributions get mapped to density profiles. By contrast, in our case we find that the real space density profile in the trap determines momentum distribution after expansion (see equation (15)). While we do not discuss the appropriate time evolution of ω(t), m(t) and v(t) here, we point out that the time-of-flight 'far-field' limit [44] may also be captured formally by our scaling approach when the asymptotics of L(t) are linear and the contribution of the quadratic phase factor in equation (15), m(t)L(t), vanishes in the long-time limit.

Conclusions and outlook
We have used scaling ansatz to show that certain quantum non-equilibrium problems with time-dependent parameters can be related to equilibrium problems with constant parameters provided that the time-dependent parameters satisfy a system of self-consistency equations. This approach is valid for rather general types of interactions and is not linked to the integrability of the model. However, an integrable structure, when it exists, is consistent with the scaling transformation. Solvability by the scaling ansatz is a consequence of the non-relativistic dynamical symmetry, which received considerable attention recently in relation to the nonrelativistic version of AdS/CFT correspondence [45]- [49]. The appearance of this symmetry in realistic many-body systems, which we discuss in this paper, can open up intriguing connections to the concept of AdS/CFT correspondence.
We have used the scaling approach to analyze the problem of an abrupt switching off of a confining potential for bosonic systems with contact interactions in d = 1 and 2. Such experiments can be performed using either ultracold atoms or photons in a nonlinear medium. We find that the asymptotic momentum distribution is essentially given by the initial density profile-a phenomenon that has previously been discussed only in the (Tonks-Girardeau) limit of the infinitely strong repulsive 1D Bose gas [30]. Possible future applications of the scaling ansatz include interaction quenches or transport phenomena (by considering finite linear potentials). Extensions of our method to systems with dissipation are also possible.
In our analysis, we have considered the situation when the scaling ansatz is obeyed exactly. We expect, however, that our results remain qualitatively valid even for systems with small deviations from the exactly scalable Hamiltonians. For example, weak lattice potentials should not have dramatic effects as long as the effective mass approximation is applicable. Therefore one could achieve a full description of time-of-flight experiments if the lattice potential and interactions are tuned accordingly. Moreover, it is conceivable that at a phenomenological level, the ansatz can be used even when the time and space dependences of system parameters do not fully satisfy the consistency equations. The scaling solution could then be seen as a universality class of non-equilibrium systems, very much like a renormalization group fixed point at equilibrium. It would be interesting to address this conjecture in experiments.
We consider the ansatz (4) for the transformation between the many-body Schrödinger equation with time-dependent parameters (equation (1)) and equation (2) with time-independent coefficients. Calculating directlẏ where for the sake of brevity we introduced φ( and where the dot denotes the derivative with respect to t, and Substituting this into the initial Schrödinger equation (1) with time-dependent coefficients, and adding and subtracting the term A(t) i x 2 i with a yet to be determined function A(t), we regroup the different contributions in front of (y i , τ ), ∂ (y i , τ )/∂y i and y i . Each group has several contributions proportional to x 0 i , x i , x 2 i that are linearly independent and must be treated separately. This is how conditions expressed by equation (7) appear. The remaining equation has the form of a Schrödinger equation with time-dependent coefficients We note that to compensate for the terms appearing after the change x i → y i in the quadratic potential, we get terms proportional to ω 2 0 in equations (6)- (12). Now, requiring that the three unknown functions τ , L(t) and A(t) satisfẏ we obtain the remaining conditions in the set of equations (5)- (12). Under these conditions, the Schrödinger equation for the function (y, τ ) has no time-dependent coefficients. From conditions (A.6) above, we determine the function Therefore, we find that when pairwise potentials obey equation (3) and the system of equations (5)-(12) is satisfied, equation (1) is indeed mapped to equation (2).

B.1. General properties of the Ermakov and related equations
In this appendix, we briefly review some general properties of the Ermakov (sometimes spelled Yermakov) equation that plays such a fundamental role in our formalism. We also point out the relationship of this equation with the Riccati equation and with the linear differential equation with variable coefficients. The Riccati equation directly appears in our approach in some limiting cases. The Ermakov [39] equation is defined as Here, a is some t-independent constant. If there is a nontrivial solution of the second-order differential equation then the transformation puts the Ermakov equation into the form where the subscript denotes the derivative. The solution of the initial equation then follows immediately, where C 1,2 are arbitrary constants. If we take two solutions of the linear (Hill) equation to satisfy initial data x 1 (0) = x 1 ,ẋ 1 (0) =ẋ 1 while x 2 (0) = 0,ẋ 2 = 0, then a general solution of the Ermakov equation is given by a nonlinear superposition principle, where w = x 1ẋ2 − x 2ẋ1 is a constant Wronskian. Now, provided the linear equation for x(t) is satisfied, the function u(t) defined as satisfies the Riccati equation, This demonstrates that all three equations are closely related: the Ermakov equation, the linear second-order differential equation with variable coefficients and the Riccati equation. Other remarkable equations are also connected to the Ermakov equation. For example, taking a = 1 for simplicity in (B.4) and defining ξ(t) = z(t) −2 , we obtain ξξ − (3/2)(ξ ) 2 + 2ξ 4 = 0. Now, defining w(t) via ξ(t) = αẇ/w with α 2 = −1/4, we obtain a Kummer-Schwarz equatioṅ w ... w −(3/2)(ẅ) 2 = 0. In some limiting situations (e.g. ω 0 = 0, see the next appendices), the Riccati equation appears naturally in our approach, so we sketch some of its properties here. The general Riccati equation with time-dependent coefficientṡ can be transformed into the second-order differential equation by the following substitution y(t) = exp(− f (t)u(t) dt). In many cases, a particular solution of (B.10) is easier to find than the one for (B.9).
The Riccati equation has a remarkable property: if there is a known particular solution u 0 (t) of (B.9), then the general solution of (B.9) is given by where C is an arbitrary constant. The particular solution u 0 (x) corresponds to C = ∞. The property (B.11) allows the construction of many solutions of (B.9) for given functions f (t), g(t) and h(t). If, for example, f (t) = 1, g(t) is arbitrary and h(t) = −(a 2 + ag(t)), a particular solution is u 0 (t) = a, and a general solution is then for arbitrary C. For example, for f (x) = 1, g(x) = 0, h(x) = bx n , we obtain where λ is a root of λ 2 + λ + b = 0.

B.2. Relation to dynamical symmetry
The Ermakov equation has the symmetry algebra isomorphic to SL(2, R), which is isomorphic to the algebra SO(2, 1) of rotations on the surface of a one-sheet hyperboloid. The property (B.11) of the Riccati equation is related to the covariance of the Riccati equation with respect to the fractional-linear transformations generated by the action of SL(2, R) algebra: the general solution can be expressed as a combination of particular solutions. The same algebra (more explicitly, one of its form, SU(1,1)) appears as a dynamical symmetry of the quantum harmonic oscillator, where the Ermakov equation appears as well. This was first found in [16]. There, a single quantum harmonic oscillator with time-dependent frequency has been solved using the methods of (adiabatic) invariants. An adiabatic invariant in this case is a function of a solution of the Ermakov equation. This approach has led to the appearance of the Ermakov-Pinney-type equation [39] in quantum mechanics (see e.g. [20] for a recent review). In [17], the same equation appears as a certain consistency condition on the timedependent rescaling of coordinate and time in the wave function of the oscillator. It became clear that these two approaches, one based on dynamical invariants and the other on the scaling of dynamical variables, are equivalent. Indeed, the rescaling procedure can be regarded as a transformation, generated by a certain symmetry group, i.e. SL(2, R). The generators of this symmetry are operators corresponding to dynamical invariants. Therefore the successiveness of applicability of scaling transformation implies the presence of dynamical symmetry generated by the dynamical invariants [18,19]. For this symmetry to hold, one has to have a special class of potential terms in the single-particle Hamiltonian [21]. Physically interesting potentials correspond to the contact interaction, harmonic, Coulomb and inverse square laws. That is why the scaling approach has been applied to a Calogero-Sutherland model [22] and classical Gross-Pitaevski-type systems [23,24,26,27,29]. The appearance of the SU(1, 1) dynamical symmetry in our non-relativistic systems suggests a possible connection to the non-relativistic version of the AdS/CFT correspondence [45]- [49]. In fact, the Virasoro algebra of any conformal field theory contains SU(1, 1) as subalgebra.

B.3. Specific solutions for ω 0 > 0
We compare examples of decreasing trapping potential and constant, increasing and decreasing masses.
(a) Constant mass. For the case of constant mass m(t) = m 0 , we choose an exponential decrease of the potential ω(t) = ω 0 e −t/τ ω . The two independent solutions of the homogeneous equation (B.2) read x 1 (t) = J 0 (2τ ω √ ω(t)) and x 2 (t) = Y 0 (2τ ω ω 0 √ ω(t)). In figure B.1, the resulting scaling functions obeying the initial conditions L(0) = 1, F(0) = 0 are plotted. For sufficiently small τ ω , the functions are well described by the limit τ ω → 0, for which the scaling solution reduces to We choose m(t) = m 0 e t/τ m and, for the sake of simplicity, ω(t > 0) = 0. The solution then reads  We emphasize that the solutions do not depend on the dimensionality of the system; only the interaction constants, which have to fulfill the consistency equation (12), will do so.
By introducing U (t) = d dt log(c(t)m(t)), it reduces to the Riccati equatioṅ The scaling ansatz (4) implies the relationslhip between initial conditions of the two systems: The initial condition for the function U (t) is not so important for us because of the special property of the Riccati equation, related to the Bäcklund symmetry, which allows one to interrelate solutions with different initial conditions via a rational function.
We note that the same equation describes the evolution of spin in a time-dependent magnetic field. A general way to solve it is to note that under some change of variables it can be reduced to the second-order linear differential equation, Numerous explicit solutions are possible if we specify the functions ω(t), m(t).
In the 2D case, we obtain from (5)-(12) that R(t) ≡ L(t) and time-dependent parameters are connected by the constraint c(t)m(t) = c 0 . Then which is a Riccati equation for the coordinate scaling function L(t); its solution for given timedependent parameters m(t) and ω(t) then defines a solution for the time-rescaling function .
In the simplest case of m(t) = 1, ω(t) = 0, we obtain c(t) = −1/ (1 + t). This example is a many-body analogue of the solution of the Hamiltonian with potential V (x) = c(t)δ(x) found in [21] for a single-particle Schrödinger equation. Direct application of this solution can be found in the ultracold Bose gas close to the confinement-induced resonance [50].
Other examples of solutions of (B.9) can be found in the literature, see e.g. [51].

Appendix C. Classical integrability of the nonlinear Schrödinger equation (NSE) with time-dependent parameters
It is instructive to check whether the exact scaling transformation we have studied in this paper is consistent with the property of integrability of the NSE. Here, we address this question for the classical NSE.
In the zero curvature representation, the NSE is represented by the system of first-order differential equations such that the matrices U (x, t, λ) and V (x, t, λ) that depend on the spectral parameter λ satisfy the condition ∂U ∂t − which is equivalent to the compatibility condition of the system, and which is equivalent to the initial Schrödinger equation. In the case of (C.1), one can establish that Conserved quantities are constructed from the matrices U and V in a known way. This method provides a direct way to various generalizations of NSE. In particular, one can obtain some generalization where the interaction parameter c and the mass are explicitly time-dependent functions. Introducing generalization of (C.5) as one can look for generalizations of integrable NSE by appropriately choosing the functions α(x, t), β(x, t), γ (x, t), λ(x, t), µ(x, t), A, B and D. Analysis of the zero-curvature condition (C.3) in the case of inhomogeneous time-dependent functions leads to a set of equations between those functions and reveals a large class of solutions of the classical equations of motions for NSE with time-dependent coefficients. To get a consistency condition for a zerocurvature representation, we conclude that the spectral parameter should be an inhomogeneous time-dependent function.
A restricted form of this inhomogeneous time-dependentŨ -Ṽ pair has been considered in [52], where it was shown that a combination of space-time transformation together with a U (1) gauge transformation of the linear equations for theŨ -Ṽ pair and corresponding redefinition of the field variables brings the system into the form of a homogeneous timeindependent NSE system, thus showing the integrability of a time-dependent system. We note that a similar analysis has been given in [53].
Although it is more difficult to show integrability on the quantum level directly, presumably the property of integrability is not violated in that case for a specific choice of timedependent parameters that correspond to our scaling equations. A related approach based on the inhomogeneity of spectral parameters for the quantum sine-Gordon model has recently been presented in [54].

Appendix D. Scaling of correlation functions
With the scaling ansatz (4), the relationship between the single-particle correlation functions in the time-dependent and time-independent systems is derived straightforwardly, The labels in the g 1 -function refer to the time-dependent ( ) and time-independent ( ) systems. From this expression, we can readily extract the density: ρ ( ) (x, t) = g 1 (x, x, t) = (1/L(t))ρ ( ) (x/L(t); 0). The momentum distribution of a time-dependent system, defined as is then given by

(D.3)
Note that because of the quadratic term in the exponent, the integrations are nontrivial. For the two-particle density matrix, we find analogously and the two-particle correlation function reads Other useful quantities, such as non-equilibrium time-dependent correlation functions (e.g. n ( ) ( p, t, t )) or the multi-mode squeezing spectrum (S(k, k ; t, t ) = n ( ) ( p, t)n ( ) ( p , t ) ), can also be easily computed using the scaling approach.

E.1. Trapped weakly interacting Bose gases
In order to describe a condensed Bose gas in a harmonic potential, we adopt the results of previous works [40]- [42] that consider phase fluctuations on top of the mean-field solution while density fluctuations are assumed to be negligible. This is a valid approximation for a sufficiently high number of weakly interacting particles at low temperatures. The temperature range where the density fluctuations are suppressed is T d T T φ where the temperature of quantum degeneracy is T d = Nhω 0 and T φ = T dh ω 0 /µ.
Generically, the single-particle correlation can be represented as where φ(x) denotes the average over phase fluctuations. We assume the validity of the Thomas-Fermi approximation for the density where R TF = √ 2µ/m 0 /ω 0 is the Thomas-Fermi radius. In a 1D geometry, taking into account thermal fluctuations and neglecting contributions from quantum fluctuations, one obtains the phase average [40] For the 2D case, an expression similar to the 1D case can be derived. In this work, we used the complete expression obtained by Xia et al (equation (77) in [41]), which explicitly accounts for thermal and quantum fluctuations. As a result, at inter-particle distances much smaller than 2R TF , the correlations decay exponentially with a decay rate approximately given by mk B T /(2πh 2 ρ(0)). However, for the dynamics studied in this paper, we did not find significant effects from quantum corrections.

E.2. 1D and 2D uniform Bose gases
For a 1D Bose gas, it was recently shown [55] that the effective field theory (Luttinger liquid) provides an extremely accurate description for a single-body correlation function at distances beyond the inter-particle separation. If we are not interested in its large momentum behavior, it is legitimate to use this effective theory. The single-particle correlation function in timeindependent theory is then well known (see e.g. [1]). For non-zero temperatures, it is given by (we omit oscillating terms) where ξ T =hv s /T =h 2 πρ/(m 0 K T ), ρ 0 is the uniform equilibrium density, v s is the sound velocity, K is a Luttinger parameter that is related to the interaction strength c and B = (K /π) 1/2K is Popov's factor. In the 2D case, we consider a system below the Berezinskii-Kosterlitz-Thouless (BKT) transition. The correlation functions then decay algebraically with a temperature-dependent exponent, which tends to the universal value 1/4 when approaching the BKT transition from below.

F.1. Relating systems in the trap and without it
The scaling approach can be used to establish the relationship between correlation functions in the model with time-dependent parameters (the system ) and the model with time-independent parameters (the system ). As we discussed in the main text, the trapping frequency ω 0 of the time-independent system is not fixed a priori. In particular, it can be set to zero from the very beginning. The scaling transformation therefore will relate the system in the time-dependent trap and a uniform system. The set of differential equations has to be modified accordingly. The aim of this appendix is to look into the behavior of the momentum distribution in this case.
The initial conditions state that the two wave functions are equal at t = 0. This means that the density distribution of the trapped system is homogeneous, corresponding to the uniform one. This is possible if we assume the existence of a length scale l on which this condition can be satisfied. Moreover, we assume here that the Thomas-Fermi radius of a trapped system is large enough such that there is a finite region of x ∈ [ − l, l] D where the density is considered to be constant. In the absence of a trapping potential, this region is equal to the whole observation area. We assume that this region is large enough to contain a relatively large number of particles N . Using this length scale l as a sort of cut-off, we evaluate the momentum distribution in the finite window [ − l, l] for examples of 1D and 2D systems at finite temperature.

F.2. Evaluation of momentum distributions in 1D for the uniform system
In the Luttinger liquid approximation at finite temperature, we introduce ξ ± = π(x − y)/ξ T in terms of which g ( ) 1 (x, y; 0) ∼ (sinh ξ − ) −1/2K . This function decays exponentially at large distances and the limits of integration in ξ − domain can therefore be extended from [−l, l] to (−∞, ∞) to make analytic progress. One can easily realize that because of the additional structure in the exponent of equation (D.1), the expression for the momentum distribution is essentially different from the one at equilibrium. The corresponding integral is where C(t) = F(t)L 2 (t)ξ 2 T ξ + /π 2 + L(t) pξ T /π . The integration over ξ + is then performed in the finite interval [−l, l] corresponding to the size of the selected subsystem. Expression (F.2) is proportional to the equilibrium momentum distribution at t = 0 provided that we take L(0) = 1. We find whereL(t) = L(t)ξ T /π . On the basis of this expression, we have calculated a momentum distribution for various particular functions ω(t) and m(t). Solving the set of consistency equations of section 2, we obtained all the other functions v(t) L(t) and F(t). This is illustrated in figure F.1 for particular choices of time-dependent functions ω(t) and m(t) and corresponds to a particular function v(t) found from the solution of the Riccati equation. But additional simulations with various other choices of functions ω(t) and m(t) suggest that the resulting momentum distribution defined as above in equation (15) has a step-like form. The formation of an effective momenta scale is associated with the asymptotic emergence of microcanonical-type distribution.
The Luttinger liquid expression for the g 1 -correlation function is a low-energy approximation for the true behavior of the correlation function. However, in the non-equilibrium dynamics, we excite the whole spectrum and therefore the result for our time-dependent theory based on the exact equilibrium theory may appear to be different from the one based on the low-energy approximation. In what follows, we demonstrate that the long-time behavior of the momentum distribution of the time-dependent system has a bounded support in momentum space. Our arguments can be applied to any exactly solvable models.
Suppose the g 1 (x, y)-correlation function is defined as a ground-state correlator of some field operators (x), † (x): g 1 (x, y) = † (x) (y) . We also assume that the matrix elements of the operator (x) in the eigenbasis of the equilibrium problem are known. This implies that the form factors F({λ}, {µ}) = {λ}| (0)|{µ} and the norms of the eigenstates |λ and |µ are known. Here, {µ} and {λ} are the sets of numbers that characterize the eigenstates of a system of size 2l. In particular, these numbers can correspond to the solutions of the Bethe ansatz equations in the exactly solvable problems. We also assume space-and time-translation invariance. Therefore the time-dependent g 1 function can be expanded as where E λ and P λ are, respectively, the energy and momentum of the state |λ . We assume also that the set {λ} corresponds to the ground state. By introducing the coordinates ξ = x − y and η = x + y, the momentum distribution of the time-dependent system (we take for simplicity equal-time correlation function) can be written as n( p, t) = L(t) where Si(z) is a sine-integral and x σ = (P µ − P λ )l/L(t) + lp + σ F(t)l 2 . The integrand is essentially proportional to F −1 (t) sin[(P µ − P λ + pL(t))/L(t)]/[(P µ − P λ + pL(t))/L(t)] and gives the main contribution to the sum when the momentum transfer is equal to pL(t).
In figure F.2, we plot the asymptotic behavior of the momentum distributions for various values of η. Similar to the 1D case, we find a step-like distribution that is smeared off when the BKT transition is approached.