Low-energy effective theory of non-thermal fixed points in a multicomponent Bose gas
Aleksandr N. Mikheev Christian-Marcel Schmied Thomas Gasenzer







aleksandr.mikheev@kip.uni-heidelberg.de
Kirchhoff-Institut für Physik, 
             Ruprecht-Karls-Universität Heidelberg,
             Im Neuenheimer Feld 227, 
             69120 Heidelberg, Germany


christian-marcel.schmied@kip.uni-heidelberg.de
Kirchhoff-Institut für Physik, 
             Ruprecht-Karls-Universität Heidelberg,
             Im Neuenheimer Feld 227, 
             69120 Heidelberg, Germany
Dodd-Walls Centre for Photonic and Quantum Technologies,
             Department of Physics, 
             University of Otago, 
             Dunedin 9016, 
             New Zealand


t.gasenzer@uni-heidelberg.de
Kirchhoff-Institut für Physik, 
             Ruprecht-Karls-Universität Heidelberg,
             Im Neuenheimer Feld 227, 
             69120 Heidelberg, Germany
             
             
                      





Non-thermal fixed points in the evolution of a quantum many-body system quenched far out of equilibrium manifest themselves in a scaling evolution of correlations in space and time. We develop a low-energy effective theory of non-thermal fixed points in a bosonic quantum many-body system by integrating out long-wavelength density fluctuations. The system consists of N distinguishable spatially uniform Bose gases with U(N)-symmetric interactions. The effective theory describes interacting Goldstone modes of the total and relative-phase excitations. It is similar in character to the non-linear Luttinger-liquid description of low-energy phonons in a single dilute Bose gas, with the markable difference of a universal non-local coupling function depending, in the large-N limit, only on momentum, single-particle mass, and density of the gas. Our theory provides a perturbative description of the non-thermal fixed point, technically easy to apply to experimentally relevant cases with a small number of fields N. Numerical results for N=3 allow us to characterize the analytical form of the scaling function and confirm the analytically predicted scaling exponents. The predicted and observed exponentially suppressed coherence at short distances takes the form of that of a quasicondensate in low-dimensional equilibrium systems. The fixed point which is dominated by the relative phases is found to be Gaussian, while a non-Gaussian fixed point is anticipated to require scaling evolution with a distinctly lower power of time.





03.65.Db 	
03.75.Kk, 	

05.70.Jk, 	

47.27.E-, 	


47.27.T- 	





Low-energy effective theory of non-thermal fixed points in a multicomponent Bose gas
    Thomas Gasenzer
    December 30, 2023
====================================================================================







§ INTRODUCTION



Relaxation of quantum many-body systems after a quench far out of equilibrium has been  studied intensely during recent years. Little is known about the general structure of possible evolutions. Various scenarios have been proposed and observed, including prethermalization <cit.>, generalized Gibbs ensembles <cit.>, many-body localization <cit.>, critical and prethermal dynamics <cit.>, decoherence and revivals <cit.>,  turbulence <cit.>, and the approach to a non-thermal fixed point <cit.>.
The rich spectrum of different possible phenomena highlights the capabilities of quantum dynamics as compared to what is possible in classical statistical ensembles.

An important difference concerns the phase angle of the quantum mechanical wave function and the associated superposition principle. In the case of a quantum many-body system, the phase angle gives rise to interference effects and decoherence and encodes the collective dynamics of the fundamental field degrees of freedom. For example, long-range coherence and thus stiffness of the phase forms the basis of sound excitations on the top of a Bose-Einstein condensate of (weakly) interacting particles. This is related to local (quasi) particle number conservation reflecting a U(1) symmetry of the underlying model description.

Here we focus on universal scaling dynamics in the relaxation of a dilute Bose gas quenched far out of equilibrium. Universal dynamics depends on a few basic symmetry properties only and thus can be classified independently of the details of microscopic properties and initial conditions. Scaling dynamics has been discussed for classical systems almost as long as spatial scaling alone. From dynamical critical phenomena <cit.> the discussion extended to coarsening and phase-ordering kinetics <cit.>, glassy dynamics and ageing <cit.>, (wave-)turbulence <cit.>, and its variants in the quantum realm of superfluids <cit.>. Different types of pre­thermal and universal dynamics after quenches of quantum many-body systems far out of equilibrium have been studied recently 
<cit.>, many of them in the context of ultracold Bose gases.
Non-thermal fixed points have been proposed, without <cit.> and with <cit.> reference to order-parameters, topological defects, and ordering kinetics, paving the way to a unifying description of universal dynamics.






A major part of the analytical work on non-thermal fixed points is based on scalar model systems with quartic interactions between N-component Bose fields, O(N)-symmetric under orthogonal transformations in the space of field components. A non-perturbative large-N approximation <cit.> allows for a description of scaling at non-thermal fixed-points <cit.>.

Here we develop a low-energy effective field theory (EFT) of such systems. We describe the linear phase-angle excitations around ground states with broken U(N) symmetry. Their non-linear bare interactions are, in general, non-local and, at the fixed point, characterized by the momentum-dependent coupling which scales analogously to the resummed couplings in the non-perturbative theory. We use this to outline a complementary approach to non-thermal fixed points which is based on a leading-order coupling expansion, practically applicable for small N≥1, and discuss consequences for N=1. 

We furthermore study numerically universal scaling dynamics in an N=3 component (d=3)-dimensional dilute Bose gas to test our analytical predictions. Our simulations corroborate the predicted scaling exponents and at the same time point to modifications of the pure scaling behavior at relatively early evolution times. While the non-thermal fixed point is conjectured to be approached asymptotically in time, prescaling of the short-distance correlations which are more easily accessible in experiment, can be seen at much earlier evolution times <cit.>.


Our paper is organized as follows:
In the remainder of this section, we introduce the model system under consideration. In LEEFT, we develop the non-linear low-energy effective theory and discuss the limit of infinitely many number of components N →∞. In this limit, in Scaling, we make predictions for the scaling at a non-thermal fixed point using the kinetic equation derived in the same section. Finally, in TWA, we show the numerical results for N = 3 component Bose gas in d = 3 dimensions, partially reproducing data from <cit.>, and discuss the central properties of the fixed point. We summarize and draw our conclusions in Conclusions. The appendices contain further details.
  



 §.§ Model


We consider a system of N spatially uniform Bose gases of identical particles. The different gases are distinguished by, e.g., the hyperfine magnetic level the gas atoms are in.
They are described by a U(N)-symmetric Gross-Pitaevskii (GP) model with quartic contact interaction in the total density,


    H_U(N) =   ∫d^dx  [ -Φ_a^†∇^2/2mΦ_a  +  g/2  (Φ_a^†Φ_a)^2] .


Here and in the following, we use units implying ħ=1, space-time field arguments are suppressed, and m is the particle mass.
Summation over the a=1,...,N Bose fields which are obeying standard commutators [Φ_a(𝐱,t),Φ_b^†(𝐲,t)]=δ_abδ(𝐱-𝐲) is implied.
The identical interspecies and intraspecies contact interactions are parametrized by the coupling g. As a result, the model exhibits a full U(N) symmetry under unitary transformations 𝒰∈U(N) of the fields, Φ_a→𝒰_abΦ_b. The generalization to inhomogeneous systems in a trap is possible but disregarded here.




 §.§ Universal scaling dynamics at a non-thermal fixed point


The nearly condensed gas is assumed to exhibit long-range order in the total phase, while domain walls and topological defects of any kind are assumed to be subdominant during the time evolution. 
Relative phases, though, can be strongly excited, representing particlelike and holelike Goldstone excitations with single-particle dispersion as further discussed below. 

This can be achieved in a system with sufficiently large N by, e.g., a strong cooling quench or an instability: 
An extreme version of a cooling quench would be to first tune adiabatically to a chemical potential 0 < -μ≪ k_B T_c, with the condensation temperature T_c, and then remove all particles with energies higher than a certain energy scale ω(k_q) ∼ |μ | <cit.>. 
Such an extreme initial condition, in experiment, can alternatively be prepared by means of an instability <cit.>.
In both cases, the crucial condition is to build up strong over­occupation such that the majority of particles and energy is around the small but non-zero momentum scale k_q, ideally close to the healing-length scale set by the chemical potential μ <cit.>. 
This induces a far-from-equilibrium evolution starting from modes with the comparatively high excitation energy ω(k_q). 

During this evolution, the majority of particles is transported to lower momenta while the major part of the energy is deposited by a few particles being scattered to even higher-momentum modes, where they eventually form a thermal tail.
The highly occupied modes which take up most of the particles have excitations energies ω(k) below the scale set by μ. This allows us to use a low-energy effective theory description.

From this type of initial conditions, the system can approach a non-thermal fixed point and show universal scaling evolution <cit.>. This universal scaling in space and time at the fixed point is expected, e.g., in the occupation numbers n_a(𝐤,t)=⟨Φ_a^†(𝐤,t)Φ_a(𝐤,t)⟩ of the Bose field excitations in each component a in momentum space, in the form of


    n_a(𝐤,t) = (t/t_ref)^αn_S,a([t/t_ref]^β𝐤).


The scaling function n_S,a(𝐤)=n_a(𝐤,t_ref) depends only on a single argument, and defines the scaling form together with the exponents α, β. 
The reference time t_ref is chosen within the scaling regime.

As an example, we study, in NTFPEvolution, the evolution of a Bose gas with N=3 components in d=3 dimensions, starting from a far-from-equilibrium initial condition (grey `box' in (a)) with momentum distributions n_a(𝐤,t_0)=n_0Θ(k_q-|𝐤|), a=1,2,3, while all phases θ_a(𝐤,t_0) are random.
The data shown here and in the following are partially reproduced from <cit.>. The figure shows that, for times t≳ t_ref=31 t_Ξ, and within a limited range of low momenta, the evolution of the momentum distribution exhibits approximate scaling according to NTFPscaling.  Rescaling the momentum occupation spectra at different times they fall onto a single universal scaling function as shown in NTFPEvolutionb, in accordance with the findings of <cit.> for the case N=1. At late evolution times, t_ref=200 t_Ξ<t<350 t_Ξ, we extract α=1.62±0.37, β=0.53±0.09, see WindowFit1.

The scaling function is characterized by a plateau up to an inverse coherence-length scale k_Λ(t), which rescales in time according to k_Λ(t) ∼ t^-β. At momenta larger than this inverse coherence-length scale, |𝐤|=k≫ k_Λ(t), the scaling function takes the power-law form n_a(k)∼ k^-ζ, with ζ≃ d+1=4, confirming earlier predictions <cit.>.

Taking the Fourier transform of the momentum distribution n_a(k) one obtains the first-order spatial coherence function g^(1)_a(𝐫,t)=⟨Φ_a^†(𝐱+𝐫,t)Φ_a(𝐱,t)⟩,
which, at short distances Ξ≲ r≪ V^-1/3 (with volume V), takes an exponential form g^(1)(𝐫,t) ∼ e^-k_Λ(t) |𝐫| (see NTFPEvolutionc). Here, Ξ=  [2 m  ρ^(0) g]^-1/2 denotes the healing length corresponding to the total density. The exponential form at comparatively short distances is reminiscent of the build-up of a quasicondensate with a rescaling coherence length scale.

The scaling evolution in the vicinity of a non-thermal fixed point in a dilute Bose gas has been described in terms of kinetic equations with many-body T-matrices derived from a non-perturbative approximation of the underlying field-theoretic equations of motion <cit.>. 
As shown in <cit.>, the non-perturbative collisional properties of the bosons become relevant in the low-energy region of momenta below the healing-length scale k_Ξ = Ξ^-1, i.e., for momenta k≪ k_Ξ, where the occupation number is strongly overoccupied (cf. Figs. <ref>a and <ref>b) and exhibits momentum scaling n_a(k)∼ k^-ζ for momenta above k_Λ which defines the transition from the plateau to the power-law fall-off.

Here, we present an alternative approach in which we first reformulate the theory in terms of phase excitations only which are the relevant degrees of freedom at momentum scales k≪ k_Ξ. 
This leads to a low-energy effective theory which  takes the form of a non-linear Luttinger liquid, with density fluctuations integrated out at quadratic order, inducing cubic and quartic interactions of the phase excitations. 
The theory is used to obtain a Boltzmann-type kinetic equation with a T-matrix evaluated in leading perturbative order.
Close to the non-thermal fixed point, where the far-from-equilibrium dynamics is dominated by the over-occupied low-energy excitations at k≪ k_Ξ, this perturbative approximation remains valid because the resulting T-matrix is power-law suppressed at the low momenta transferred between these modes.
Similar to equilibrium systems in d<3 dimensions we find quasicondensate-type coherence close to the non-thermal fixed point, characterized by g^(1)(𝐫,t) ∼ e^-k_Λ(t) |𝐫| of the time evolving system in three spatial dimensions.
 









§ LOW-ENERGY EFFECTIVE THEORY



In this section, we provide details of the derivations of the low-energy effective theory for the U(N)-component dilute Bose gas (<ref>) with U(N)-symmetric contact interactions, which forms the central result of this paper.




 §.§ Collective excitations



We start our analytical discussion of the universal scaling dynamics with a brief summary of the collective low-energy linear excitations of the model ONGPH. From ONGPH, one obtains the classical action in terms of fluctuating complex fields φ_a(x), with x≡(𝐱,x^0)≡(𝐱,t), 


    S[φ] =   ∫_x{ı/2[φ_a^* (x) ∂_t φ_a (x) - φ_a (x) ∂_t φ_a^* (x) ] 
      
         
      - 1/2 m[∇φ_a^* (x) ]·[∇φ_a (x) ] 
      - g/2[φ_a^* (x) φ_a (x) ]^2} ,


where ∫_x=∫d^d+1x and sums over double indices are implied.
Using the Madelung field representation in terms of polar coordinates <cit.>,


    φ_a (𝐱,t) = √(ρ_a(𝐱,t)) exp{ıθ_a (𝐱,t)},


with densities ρ_a and phases θ_a, the Lagrangian reads as


    ℒ 
      = -∑_a {ρ_a ∂_t θ_a 
       + 1/2m[ ρ_a (∇θ_a)^2  + (∇√(ρ_a))^2 ] } 
       - g/2ρ^2,


where ρ=∑_aρ_a.
From this, a continuity equation relating the density to the particle current 𝐣_a=ρ_a∂_𝐱θ_a/m, and an equation for the phase θ_a follow, which, in the limit of small fluctuations θ_a and δρ_a=ρ_a-ρ_a^(0) about the uniform ground-state densities ρ_a^(0)= ⟨Φ_a^†(x)Φ_a(x) ⟩ [Expanding around the time-independent homogeneous solution has to be distinguished from expanding around the condensate mean field which, during the evolution in the vicinity of an NTFP, typically grows algebraically. The density corresponding to the field expectation value squared needs to be distinguished from the zero mode of the total density which is proportional to the total particle number and given by the sum of the condensate density and the integral over all non-zero momentum occupancies, ρ^(0)_a = |⟨Φ_a⟩|^2 + ∑_k⃗⟨Φ^†_a (k⃗) Φ_a(k⃗)⟩. Hence, conservation of the total particle number in a closed system ensures the constancy of ρ^(0)_a.
Furthermore, one should note that, in the presence of a condensate in the ground state of, e.g., a 3D system, one would have ρ_a^(0)= |⟨Φ_a(x) ⟩|^2, while the approach here also applies to lower-dimensional systems where, in equilibrium, the ground state bears a quasi-condensate.],
reduce to the linearized equations of motion


    ∂_t θ_a    = 1/4 m ρ_a^(0)∇^2 δρ_a - g ∑_b δρ_b ,   ∂_t δρ_a = -ρ_a^(0)/m∇^2 θ_a.
  

In Fourier space, taking a further time derivative, they can be combined to the Bogoliubov-type wave equation for the θ_a,


    ∂_t^2 θ_a(𝐤,t) + 𝐤^2/2m(𝐤^2/2mδ^ab + 2 g ρ^(0)_b ) θ_b (𝐤,t) = 0 ,


where Einstein's sum convention is implied. 
While, for N=1, we recover the Bogoliubov dispersion, for general N, diagonalization of the coefficient matrix yields the eigenfrequencies of N-1 Goldstone and one Bogoliubov mode [Obviously, all of these modes are gapless Goldstone excitations. Nevertheless, in order to distinguish them, hereinafter we will only refer to the quadratic modes as Goldstone ones, whereas the linear one will be addressed as Bogoliubov.], 


    ω_c(𝐤) 
         ≡ω_G(𝐤) 
      = 𝐤^2/2m,    c=1,…,N-1 ,
      
    ω_N (𝐤) 
         ≡ω_B(𝐤)
      = √(𝐤^2/2m(𝐤^2/2m + 2 g ρ^(0))) ,


where ρ^(0)=∑_aρ_a^(0) is the total condensate density (see Appendix <ref> for details). Note that the Goldstone theorem <cit.> predicts, due to the spontaneous breaking of the U(N) →U(N-1), 2N - 1 gapless Goldstone modes. However, only N of these modes, with frequencies GoldstoneFreq and BogHiggsFreq, are independent due to the absence of Lorentz invariance and thus particle-hole symmetry <cit.>.

The general solution can be written as a linear combination of the eigenmodes, 


    θ_a(𝐤,t)
      = ∑_b=1^N C^(a)_b e_ba cos(ω_b (𝐤)  t + ϑ^(a)_b ) ,


where e_ba=(𝐞_b)_a are the eigenvector components (<ref>) and (<ref>), and C^(a)_b and ϑ^(a)_b are real constants defined by the initial conditions. We find that the free-particle excitations are the relative phases between the different components, while the Bogoliubov excitation is in the total phase. 

We can always choose the eigenbasis such that all modes (equations) decouple due to the U(N) symmetry of the action. On the other hand, in experimentally realistic scenarios, this symmetry may be explicitly broken by the initial density matrix ρ̂_0 which fixes the mean densities ρ^(0)_a. Therefore, it is suggestive to stick with the initial “physical” basis (see also the discussion in Appendix <ref>). The solution for the density fluctuations has the same form, i.e., 


    δρ_a(𝐤,t)
      = ∑_b=1^N D^(a)_b e_ba cos(ω_b (𝐤)  t + ϕ^(a)_b ).


Taking the time derivative of deltarhoasolution, we notice that, at low energies, i.e., for k ≪ k_Ξ, where 


    k_Ξ = [2 m  ρ^(0) g]^1/2


is a momentum scale set by the inverse healing length corresponding to the total density, the Bogoliubov mode contribution to the time derivative of fluctuations dominates, i.e.,


    ∂_t δρ_a (𝐤,t) ∼ω_N (𝐤,t) δρ_a (𝐤,t).
  

Then, according to (<ref>), below k_Ξ, 


    δρ_a(𝐤)/ρ^(0)_a∼|𝐤|/k_Ξθ_a(𝐤) 
      ≪θ_a(𝐤)
         (k ≪ k_Ξ).


Here we assumed that, at low energies, 


    |D^(a)_N ω_N(𝐤) sinϕ^(a)_N | ≳|D^(a)_b ω_b(𝐤) sinϕ^(a)_b |  ,


i.e., that either both terms in ∂_t δρ_a are of the same order or that the Bogoliubov term dominates, which  depends, in general, on the initial conditions. 




 §.§ Luttinger-liquid-type effective action






  §.§.§ Derivation of the effective action


In order to demonstrate the implications of the collective excitations introduced above for the possible non-equilibrium dynamics of the system we use the density-phase representation to obtain a low-energy effective field theory of the model ONGPH with all ρ^(0)_a>0.
As for a single-component GP system, the quartic interactions proportional to the square of the local total density ρ=∑_aφ_a^*φ_a give rise to a reduction of density fluctuations at low energies. While this does not necessarily imply a suppression of the fluctuations of the single-species densities below some energy scale, such a suppression may be generated dynamically for certain initial conditions. In this section, we will derive the effective theory with a single-species density fluctuations suppression, while the more general case is considered in App. <ref>.

To begin with, we again adopt the Madelung representation and expand the Lagrangian (<ref>) up to the second order in density fluctuations. The resulting action reads:


    S[θ, δρ] 
      = S_G[θ, δρ]+ S_nG[θ, δρ] ,


where

    S_G
      = -∫_x{𝒥_a δρ_a    
      + 1/2δρ_a g^abδρ_b 
      + ρ_a^(0)[∂_t θ_a 
      + (∇θ_a )^2/2 m]  } ,
      
    𝒥_a 
          = g ρ^(0) + ∂_t θ_a + 1/2 m(∇θ_a )^2 , 
      
    
      g^ab   = g - δ^ab/4 m ρ_a^(0)∇^2 , 
      
    S_nG   = ∫_x (∇δρ_a)^2/8 mρ_a^(0)(δρ_a/ρ_a^(0))
       +h.o.t.


and we introduced ρ^(0)≡∑_a ρ_a^(0).
Note that the neglected higher-order terms (h.o.t.) involve powers of the density fluctuations only, while all phase-dependent terms have been taken into account. Hence, there is no approximation made concerning the size of the phase fluctuations. This is relevant, e.g., for Bose gases in d < 3 dimensions where even below the BEC crossover the quasicondensate exhibits large phase fluctuations.

As it was mentioned above, the derivation of the low-energy effective action consists of integrating out density fluctuations. This can be formally done adopting a Feynman-Vernon influence functional approach <cit.> explained in App. <ref>. In this paper, for the sake of simplicity, we assume that the initial state can be well-described by a product state of a (Gaussian) density matrix of phase fluctuations and a ground-state density fluctuations matrix, i.e., 


    ρ̂(t_0) ≈ρ̂_δρ (t_0) ×ρ̂_θ (t_0) ,   ρ_δρ (t_0) ≈ |Ω⟩⟨Ω |.


While this approximation might look extreme, it should be
noted that the kinetic theory which we are going to use in the following disdains details of the initial conditions regardless.

With the aforementioned approximations, the derivation of the low-energy effective action is reduced to the standard procedure within zero-temperature quantum field theory:


    Z = ∫ D θ  D δρ e^ı S[θ,δρ]≡∫ D θ  e^ı S_eff[θ].


In general, the integral (<ref>) is non-Gaussian and thus cannot be computed exactly. We can, nevertheless, compute it perturbatively,


    ∫ D δρ  e^ı S_G + ı S_nG = Z_0 ∑_n=0^∞ı^n ⟨ (S_nG)^n ⟩/n!,


where Z_0 = ∫ D δρ  e^ı S_G, and restrict ourselves to some finite number of terms in the above expression. For instance, at the lowest order:


    exp(ı S_eff[θ]) ≈∫ D δρ exp(ı S_G[θ,δρ]).


Equivalently, one can truncate the expansion (<ref>) at the Gaussian order from the very beginning, also resulting in (<ref>).

Since the kernel app:g_ab contains a spatial derivative, it is more convenient to proceed in Fourier space (see Not for details of the notation):


    S_G[θ,δρ]
      = -1/2∫_𝐤,t 
    𝐤',t'δρ_a(𝐤,t) g^ab(𝐤,t;𝐤',t') δρ_b(𝐤',t')   
      
       
      - ∫_𝐤,t{𝒥_a (𝐤,t) δρ_a(𝐤,t) 
      + ρ_a^(0)/2 m𝐤^2 θ_a(𝐤,t) θ_a(-𝐤,t) } ,


where the total-derivative term ρ^(0)∂_t θ_a is dropped, and where


    𝒥_a(𝐤,t) 
         =  gρ^(0)  (2 π)^d δ (𝐤) + ∂_t θ_a (-𝐤,t)  
      
    
      + 1/2 m   ∫_𝐤'𝐤' ( 𝐤 
      + 𝐤') θ_a(𝐤',t) θ_a(-𝐤 - 𝐤',t) ,
      
    
      g^ab(𝐤,t;𝐤',t') 
      =    g(1 +  δ^ab 𝐤^2/2k_Ξ,a^2)(2π)^d δ (𝐤 + 𝐤')  δ (t-t').


Here k_Ξ,a = [2 m  ρ_a^(0) g]^1/2 is a momentum scale taking the form of the inverse healing length of a single component.
We recall that no condensate is required such that the interpretation of 1/k_Ξ,a as a healing length is in general not applicable.

We can absorb the first term in the definition of 𝒥_a, app:Jk, by going over to a (grand-)canonical formulation, effectively shifting the energy of the zero-mode by a constant, 


    gρ^(0)  (2 π)^d δ (𝐤) + ∂_t θ_a (-𝐤,t)
      →∂_t θ_a (-𝐤,t).


In our path-integral formulation in terms of the fundamental Bose fields this can be achieved by shifting all densities by a background term, see App. <ref>. As a result of the above, the current field simplifies to


    𝒥_a(𝐤,t) 
         = ∂_t θ_a (-𝐤,t)  
      
    
      + 1/2 m   ∫_𝐤'𝐤' ( 𝐤 
      + 𝐤') θ_a(𝐤',t) θ_a(-𝐤 - 𝐤',t).



Note that the operator app:gkk is diagonal both in 𝐤-space and in t-space so that the only non-trivial step is to invert the matrix part of the kernel.

Before we perform the Gaussian integral over the δρ_a to obtain the effective action, we rescale the density fluctuations,


    δρ_a →√(N ρ_a^(0)/ρ^(0))δρ_a' ,


which multiplies the path-integral measure by an irrelevant constant. 
With this, we obtain a modified kernel in app:SGp, 


    g̃^ab(𝐤,t;𝐤',t') 
      =    N g [(ρ_a^(0)ρ_b^(0))^1/2/ρ^(0) +  δ^ab 𝐤^2/2 k^2_Ξ]
      
       × (2π)^d δ (𝐤 + 𝐤')  δ (t-t')


and also need to take into account the rescaling in the linear term, 𝒥_a(𝐤,t) δρ_a(𝐤,t) →𝒥_a(𝐤,t) (Nρ_a^(0)/ρ^(0))^1/2δρ_a'(𝐤,t).



To absorb these factors which take into account the deviation from the mean density per component ρ^(0)/N, we also rescale the remaining fluctuating phase-angle fields as 


    θ_a' = √(N ρ_a^(0)/ρ^(0)) θ_a ,


with the condition


    ∫𝒟θ 𝒪 (θ)   e^ı S_eff [θ] 
    = ∫𝒟θ'  𝒪 (θ)   e^ı S'_eff [θ'],


where the argument of the operator 𝒪 is understood to be expressed in terms of the rescaled field. Equation RescCond yields exp(ı S_eff[θ]) = const.×exp(ı S'_eff[θ']), with an insignificant multiplicative constant. This implies a corresponding rescaling of 𝒥_a →𝒥_a'.

As discussed in detail in Inversegab, the inverse of g̃^ab is 


    (g̃^-1)^ab(𝐤,t;𝐤',t') 
      = 1/N2 k_Ξ^2/g 𝐤^2(δ^ab - √(ρ_a^(0)ρ_b^(0))/ρ^(0)/1 + 𝐤^2/2 k_Ξ^2)
    (2 π)^d δ (𝐤 + 𝐤') δ (t - t').


After completing the squares and Gaussian integration over the δρ', and taking into account the rescaling of 𝒥_a, this kernel enters the contribution


    1/2∫_𝐤𝐤'tt'𝒥'_a(𝐤,t)
      (g̃^-1)^ab(𝐤,t;𝐤',t') 
      𝒥'_b(𝐤',t')


to the effective action which is quadratic in the 𝒥_a. The effective action finally becomes:


    S'_eff[θ'] 
      = ∫_𝐤,t1/2 N{1/g_1/N 
      (𝐤)(δ^ab - k_Ξ,a k_Ξ,b/k^2_Ξ/1 + 𝐤^2/2 k^2_Ξ) 
      × 𝒥'_a (𝐤,t) 𝒥'_b (-𝐤,t) 
     - ρ^(0)𝐤^2/mθ'_a(𝐤,t) θ'_a (-𝐤,t)},


where


    𝒥'_a(𝐤,t)
      =  
      ∂_t θ'_a (-𝐤,t)  
    + (ρ^(0)/N ρ^(0)_a)^1/2∫_𝐤'𝐤' ( 𝐤 + 𝐤')/2 mθ'_a(𝐤',t) θ'_a(-𝐤 - 𝐤',t) ,


with momentum-dependent “coupling function”


    g_1/N(𝐤) = g 𝐤^2/2 k^2_Ξ.


The effective action can be split into a Gaussian and a non-Gaussian interaction part.
Omitting primes in denoting the fields and the action in what follows, we get


    S_eff=S_eff,G+S_eff,nG^(3)+S_eff,nG^(4) ,


with


    S_eff,G[θ] 
      = ∫_𝐤,t1/2 N{1/g_1/N 
      (𝐤)(δ^ab - k_Ξ,a k_Ξ,b/k^2_Ξ/1 + 𝐤^2/2 k^2_Ξ) 
      ∂_tθ_a (𝐤,t) ∂_tθ_b (-𝐤,t) 
     - ρ^(0)𝐤^2/mθ_a(𝐤,t) θ_a (-𝐤,t)},


and the three- and four-wave interaction parts


    S_eff,nG^(3)[θ] 
         = ∫_𝐤𝐤',t1/N^3/21/g_1/N 
      (𝐤)(δ^ab - k_Ξ,a k_Ξ,b/k^2_Ξ/1 + 𝐤^2/2 k^2_Ξ) 
      k_Ξ/k_Ξ,b𝐤' ( 𝐤' - 𝐤)/2 m∂_tθ_a (-𝐤,t) θ_b(𝐤',t) θ_b(𝐤 - 𝐤',t),
     
    
      S_eff,nG^(4)[θ] 
         = ∫_𝐤𝐤'𝐤”,t1/2N^21/g_1/N 
      (𝐤)(
     δ^ab k_Ξ^2/k_Ξ,a^2
      - 1/1 + 𝐤^2/2 k^2_Ξ) 
      𝐤' ( 𝐤' + 𝐤)/2 m𝐤” ( 𝐤” - 𝐤)/2 m
       × θ_a(𝐤',t) θ_a(-𝐤 - 𝐤',t)  
      θ_b(𝐤”,t) θ_b(𝐤 - 𝐤”,t).




Note that the quadratic term in SeffGthetaa, for a single-component gas (N=1), reduces to the standard Luttinger-liquid  action [The Hamiltonian of the (1D) (Tomonaga-)Luttinger-liquid model <cit.> reads H_TLL=(c_s/2π)∫ dx {K^-1[πΠ(x,t)]^2+K[∂_xθ(x,t)]^2}, with speed of sound c_s=√(v_Jv_N), Luttinger parameter K=√(v_J/v_N), phase stiffness v_J=πρ_0/m, and density stiffness v_J=g/π. 
The canonically conjugate momentum is defined in terms of the density fluctuation as Π(x,t)=δρ(x,t).
This corresponds to the Lagrangian L_TLL=(2g)^-1∫ dx {[∂_tθ(x,t)]^2-[c_s∂_xθ(x,t)]^2}, and thus to the canonically conjugate field Π(x,t)=g^-1∂_tθ(x,t).]
accounting for non-interacting sound modes on the background of a (quasi) condensate, useful for describing equilibrium low-temperature Bose gases in one and two spatial dimensions <cit.>.
In the following, the non-linear terms SeffnG3thetaa and SeffnG4thetaa, however, will be crucial to describe the kinetic transport of phase excitations close to the non-thermal fixed point.

In equilibrium, the Luttinger liquid serves as a low-energy effective theory for near-zero-temperature Bose gases in less than three dimensions, where no spontaneously broken field is available as a Bogoliubov mean field to expand around  <cit.>. 
Here, we will exploit that the same low-energy effective theory, with non-linear terms taken into account, serves to describe transport in momentum space within a kinetic formulation of the time evolution.
We emphasize that the formulation of the low-energy effective theory is independent of the particular state the system is in, provided the assumption of having a constant high mean density in each field component to expand about is well satisfied.
Hence, far-from-equilibrium phenomena can be described by our theory provided the mode occupancies of the state are dominated by the low-energy range of momenta below k_Ξ.





  §.§.§ Large-N limit of the Luttinger-liquid-type effective action



To allow comparisons with previous formulations of scaling at a non-thermal fixed point, including <cit.>, it is necessary to consider the limit N →∞. Using the fact that k_Ξ,a/k_Ξ∼ 1/√(N) one easily obtains that all the “tensor structures” in SeffGthetaa, SeffnG3thetaa, and SeffnG4thetaa become diagonal to the leading order in 1/N [One has to be careful, however. Even though each off-diagonal term is of the order O(1/N), there are N - 1 of them. Therefore, the total contribution is of the order O(1). Furthermore, the very presence of the Bogoliubov mode is due to these off-diagonal elements, and hence neglecting them will result in losing the linear dispersion mode, which we will see immediately. Thus, strictly speaking, the following effective action accounts only for a part of the full set
of possible dynamical excitations.].

Summarizing the above derivations and going back to the non-rescaled phase fields, we obtain, in the large-N limit, the following low-energy effective action (see  vertices for details)


    S_eff[θ] 
      = ∫_k⃗, k⃗', 𝒞1/2 θ_a (k⃗, t) ı D_ab^-1(k⃗, t; k⃗', t') θ_b(k⃗', t')
      
       - ∫_{k⃗_i }, 𝒞k_Ξ,a^2/k_Ξ^2k⃗_1 ·k⃗_2/2m g_1/N (k⃗_3) θ_a (k⃗_1, t) θ_a (k⃗_2, t) ∂_t θ_a (k⃗_3, t)  δ(∑_i=1^3k⃗_i) 
      
       + ∫_{k⃗_i }, 𝒞k_Ξ,a^2/k_Ξ^2(k⃗_1 ·k⃗_2)   (k⃗_3 ·k⃗_4)/8m^2  g_1/N (k⃗_1 - k⃗_2) θ_a(k⃗_1,t) ⋯θ_a(k⃗_4,t) δ(∑_i=1^4k⃗_i) ,
 

with free inverse propagator, defining δ_k⃗, -k⃗'≡ (2 π)^d  δ (k⃗ + k⃗'),


    ı D_ab^-1(k⃗, t; k⃗', t') 
    = δ_k⃗, -k⃗'/Ng_1/N (k⃗) δ_abδ_𝒞 (t - t') (-∂_t^2 - (k⃗^2/2m)^2 ).


We note that the large-N limit can be taken either for a fixed total density ρ^(0), or keeping the single-component densities ρ_a^(0) constant. 
In the former case, the single-component densities vanish for N→∞, such that one eventually obtains a classical gas of distinguishable bosons. As a result, the condition δρ_a/ρ_a^(0)≪θ_a, cf. CondLEEFT, breaks down asymptotically as the particles become distinguishable and lose mutual coherence.

Hence, the low-energy effective theory, in the large-N limit, rather applies to the situation of fixed, i.e., N-independent, densities ρ_a^(0) equal in each of the N effectively decoupled copies of field components.
To avoid, in this case, that the gas loses its diluteness and becomes a dense fluid as the total density ρ^(0) = Nρ_a^(0) [Hereinafter, we assume, for simplicity, that all components are equally occupied, ρ^(0)  = N ρ^(0)_a. Since in the large-N limit there is no coupling between different modes, this does not effect the following analysis.] increases, one needs to adjust the coupling g.
A sensible choice is to keep the relevant Gross-Pitaevskii coupling, i.e., the healing-length energy scale corresponding to the total density, k_Ξ^2/(2m)=gρ^(0), fixed and thus the rescaled coupling λ=gN.
For this choice, the low-energy effective action becomes invariant under changes of N. 
Moreover, its regime of validity, with excitations on scales longer than the healing length [cf. CondLEEFT], remains fixed at large N.

We remark that rescaling the density with N and the coupling with 1/N implies that, in d>1 spatial dimensions, the diluteness parameter of the total gas scales as 


    ζ_d = [ρ^(0)]^1/da_d∼ N^1/d-1 ,


where a_d∼ g is the length scale associated with the coupling g in d dimensions.
For example, in d=3 dimensions, the diluteness of the total gas increases with N such that, in the large-N limit, the theory becomes perturbative, and quantum fluctuations can be neglected.

Taking the single-component densities fixed as discussed above, we will write all expressions, in what follows, in terms of the Goldstone coupling, 


    g_G (𝐤)≡ N g_1/N(𝐤).


Finally, we note that the interaction terms can be conveniently represented by the diagrams depicted in Fig. <ref>

 





§ KINETIC THEORY AND SCALING ANALYSIS



In the previous section, we have derived the effective action (<ref>) to describe dynamics of the U(N)-symmetric Gross-Pitaevskii model below the scale k_Ξ. To do so, we integrated out density fluctuations which, as was shown in Sect. <ref>, are suppressed by the factor k/k_Ξ as compared to phase fluctuations. Furthermore, in the limit N →∞, the effective description reduces to N independent and identical copies described by the action (<ref>). 

 Our aim now is to apply the low-energy effective theory developed in luttingerleeft to universal scaling dynamics at non-thermal fixed points. A first step in this direction is to derive a kinetic description of transport in momentum space, as has been done previously within the non-perturbative approach based on the fundamental Bose fields <cit.>.
We first set up kinetic equations and scattering integrals, in analogy to standard derivations (see, e.g., <cit.> for details).

The procedure is then to state a scaling hypothesis of the form NTFPscaling for the momentum distribution function whose time evolution is governed by the kinetic equation.
Hence, while the distribution evolves by rescaling in time and space, it is assumed to keep a stationary scaling form which is governed by a fixed-point equation corresponding to a rescaled form of the kinetic equation.
As a result, the analysis focuses, mathematically speaking, on the fixed-point configuration, providing the properties of the quasi-stationary scaling form and the scaling exponents governing the time evolution.

Most of the discussion is kept within the large-N limit, except for the case of N=1 discussed separately in SingleCompScaling. 




 §.§ Kinetic description



The kinetic description focuses on the evolution of equal-time two-point correlators, specifically on occupation-number distributions of quasiparticles in momentum space (assuming identical excitations in all modes, e.g., in a large-N approximation),


    f_a(𝐤,t) = ⟨θ_a(𝐤,t) θ_a(-𝐤,t) ⟩.


One starts with Kadanoff-Baym dynamic equations for the Schwinger-Keldysh time-ordered two-point Greens functions G_ab(x,y)=⟨𝒯_𝒞θ_a(x) θ_b(y)⟩. After decomposing into spectral and Keldysh statistical components, one introduces a quasiparticle ansatz for the spectral functions.
Finally, performing a gradient expansion, in Wigner coordinates, along the relative-time direction, and sending the initial time t_0 → -∞, one obtains, in two-loop approximation of the two-particle irreducible effective action (Luttinger-Ward functional), a kinetic equation for the time evolution of a suitably defined momentum spectrum f_𝐤≡ f(𝐤,t)≡ f_a(𝐤,t) of phase excitations,


    ∂_tf(𝐤,t)
      = I[f](𝐤,t) ,


where I[f] is a scattering integral to be specified below.

From the action Seff4, we infer (see Appendix <ref> for details) that the scattering integral has the form


    I[f](𝐤,t)=    I_3(𝐤,t)+I_4(𝐤,t) ,


with


    I_3 (𝐤,t) 
      ∼   ∫_k_0,p,q |γ(𝐤,𝐩,𝐪)|^2  
      ρ̃_k ρ̃_p ρ̃_q δ (k + p - q) 
      
       ×[ (f_k + 1) (f_p + 1) f_q  -  f_k f_p (f_r + 1) ] ,
      
    
      I_4 (𝐤,t) 
      ∼   ∫_k_0,p,q,r |λ(𝐤,𝐩,𝐪,𝐫)|^2
      ρ̃_k ρ̃_p ρ̃_q ρ̃_r δ (k + p - q - r) 
      
       × [ (f_k + 1) (f_p + 1) f_q f_r -f_k f_p (f_q + 1) (f_r + 1) ]  ,


where the (d+1)-dimensional delta distributions imply energy and momentum conservation, with k_0=ω(k)=k^2/2m denoting the quasiparticle frequencies, and where the scattering matrices are defined as


    γ(𝐤,𝐩,𝐪) 
         = (𝐤·𝐩) ω(𝐪)/m g_G (𝐪) + perm^s ,
      
    λ(𝐤,𝐩,𝐪,𝐫) 
         = (𝐤·𝐩) (𝐪·𝐫)/2m^2 g_G (𝐤 - 𝐩) 
      + perm^s ,


with momentum-dependent coupling function which, in the large-N limit, takes the form g_G(𝐤) = Ng 𝐤^2/(2 k^2_Ξ) [cf. Eqs. app:g1N and g_Gold].
Assuming that the quasiparticle (on-shell) approximation is valid, i.e., that the spectral functions ρ̃(k)=ρ̃(k_0,𝐤) describe free stable quasiparticles, we find that for our EFT the bare retarded Green's function has the form


    G_R^(0)   (𝐤, ω) 
      = g_G (𝐤)/(ω + ı 0^+)^2 - ω^2(𝐤)
       = g_G (𝐤)/2 ω (𝐤)[ 1/ω - ω(𝐤) + ı 0^+ 
      - 1/ω + ω(𝐤) + ı 0^+].


Using the definition of the spectral function in terms of G_R, ρ̃ (𝐤,ω) = -2  Im  G_R (𝐤,ω), as well as the relations


    Im1/x + ıϵ 
      = - ϵ/x^2 + ϵ^2 ,   lim_ϵ→ 0ϵ/x^2 + ϵ^2 
      = πδ (x) ,


we obtain the bare spectral function of our EFT Seff4,


    ρ̃^(0) (𝐤, ω) = π g_G (𝐤)/ω(𝐤)[ δ (ω - ω(𝐤)) - δ (ω + ω(𝐤)) ] ,


so that the T-matrices adopt the following form:


    |T_3(𝐤,𝐩,𝐪)|^2 = |γ(𝐤,𝐩,𝐪)|^2 
       g_G (𝐤)  g_G (𝐩)  g_G (𝐪)/8  ω (𝐤)  ω (𝐩)  ω (𝐪) ,


    |T_4(𝐤,𝐩,𝐪,𝐫)|^2 = |λ(𝐤,𝐩,𝐪,𝐫)|^2   
      g_G (𝐤) ⋯ g_G (𝐫)/2 ω (𝐤)  ⋯ 2 ω (𝐫).






 §.§ Scaling at a non-thermal fixed point



We have now the tools at hand to proceed to the main physical results of this paper, predictions for the scaling behavior of the considered model at a non-thermal fixed point.
The procedure is as follows.
The kinetic equation derived above is assumed to describe the transport process in momentum space of low-energy phase excitations towards infrared scales.
As illustrated by our numerical data obtained for the evolution after a strong quench (see NTFPEvolution), this transport assumes a universal scaling form after having passed a brief period of fast initial evolution.

The main task of the analytical work presented here is to formulate and evaluate fixed-point equations which allow one to predict the scaling functions and exponents characterizing this universal dynamics.
We focus on the analysis at the fixed point, i.e., we assume that the scaling evolution has already been fully established.
This allows us to insert a scaling hypothesis for the solution into the kinetic equation and apply power-counting techniques to determine the unknown scaling exponents <cit.>.
In this way we will re-derive previously known scaling exponents and complement them with results for situations which have not been considered before, in particular for the experimentally relevant case of a single field component, N=1, including in particular the three-vertex term I_3.




  §.§.§ Spatio-temporal scaling


Using the T-matrices given above, we can perform a scaling analysis of universal dynamics of the system at a non-thermal fixed point. 
To this end, we consider the scattering integrals corresponding to each vertex. 
According to the scaling hypothesis, the quasiparticle spectral distribution satisfies the scaling form


    f (𝐤,t) = s^α/β f(s 𝐤, s^-1/βt),


and so do the T-matrices,


    |T_l(𝐤_1,…,𝐤_l;t)| 
      = s^-m_l |T_l(s 𝐤_1,…,s 𝐤_l;s^-1/βt)| ,


where the index l∈{3,4} denotes the l-point vertex.
From the above, one obtains the spatio-temporal scaling of the scattering integral <cit.>,


    I_l[f](𝐤,t) = s^-μ_l I_l[f](s 𝐤,s^-1/β t).


Using I_4 and I_3 together with distr_scaling and T_scaling, one derives the expression for a scaling exponent μ_l for each collision term:


    μ_l = (l - 2)  d - z + 2m_l - (l  - 1) α/β,


where 


    2m_l =  2m̃_l + l  (γ-z),


with 


    m̃_3    = z + 2 - γ,
    m̃_4    = 4 - γ,


and γ being defined as the scaling exponent of g_G:


    g_G (𝐤) = s^-γ g_G (s 𝐤).


Note that we here assume that the scaling of the coupling g_G appearing in the retarded Green's function GR, and in the couplings, gamma3 and lambda4, is the same, i.e., the momenta as well as their differences are in the same scaling regime.
This should be guaranteed if the external momentum 𝐤 of the scattering integrals is in the infrared (IR) scaling limit and thus  the integrals are dominated by momenta in this regime.


If the distribution function f obeying the scaling form distr_scaling is to solve the kinetic equation QKinEq for a given μ=μ_l, the scaling exponents need to satisfy the relation


    α = 1 - βμ.


In addition, in the presence of global conservation laws for the integral ∫_𝐤f(𝐤,t) (quasiparticle number) or the integral ∫_𝐤ω(𝐤)f(𝐤,t) (quasiparticle energy), the following constraints apply for the scaling exponents:


    α   = β d,    number conservation; 
    α'    = β' (d + z),    energy conservation.


We note that the conservation of quasiparticle number is, in general, not expected.
In the low-energy regime of a Bose gas considered here, the quasiparticles are gapless phonon and relative-phase modes as described above whose number can change due to the cubic interaction terms appearing in the action Seff4.
It appears nevertheless as one of the key properties of non-thermal fixed points that the transport underlying the universal scaling dynamics is related to conservation laws which can include that of quasiparticle number.


Here and in the following, primed exponents α', β' refer to transport conserving energy.
Collecting the above results, we obtain


    μ_3    = d + 4 - 2 z + γ - 2 α/β ,
    μ_4    = 2 d + 8 - 5 z + 2γ - 3 α/β ,


which, together with (<ref>), gives:


    l = 3:       β· (d + 4 - 2 z + γ) = 1 + α ,
       l = 4:       β· (2 d + 8 - 5 z + 2 γ) = 1 + 2 α.


One should note that, in principle, mu3rel and mu4rel provide the closed system of equations allowing to determine α and β:


    α   = (d + 4 - 3 z + γ)/z ,
    β   = 1/z.


However, since, depending on the dimensionality d and the momentum region of interest, one type of interaction can dominate over the other one, it is more reasonable to analyze the two terms independently. 
In order to close the system in that case, an additional relation is required, which, in fact, is provided by either  quasiparticle number conservation qp_cons or energy conservation en_cons within the scaling regime. 
For particle number conservation, we obtain


    l = 3:    β =    1/4 - 2 z + γ ,
       l = 4:    β =    1/8 - 5 z + 2 γ ,


while, for energy conservation, one gets


    l = 3:     β' =    1/4 - 3 z + γ ,
       l = 4:     β' =    1/8 - 7 z + 2 γ.


Substituting, in the large-N limit, the scaling exponent z = 2 of the free Goldstone dispersion app:GoldstoneFreq and γ = 2 of the coupling function app:g1N, the T-matrices scale with


    m_3 = 2 ,    
    m_4 = 2 ,


and the resulting scaling exponents read as


    β   = 1/2 ,    (z=2)
    
    α   = d/2 ,


for both three- and four-point vertices, and


    β'    = -1/2 ,
    α'    = -(d + z)/2 ,


for the four-point vertex, while, at the same time, for the three-point vertex, no valid solution exists.
We point out that the above exponents are equivalent to the respective exponents derived in <cit.> using the large-N resummed kinetic theory for the fundamental Bose fields, for the case of a dynamical exponent z=2, and a vanishing anomalous dimension η=0 (see also ScalingFunctionBosePhase below).

The question which arises here is, whether both, three- and four-wave interactions are equally relevant.
To answer this, we need to compare the spatio-temporal  scaling properties of the scattering integrals I_3 and I_4, for given fixed-point solutions f(𝐤,t).
Focusing on the infrared transport of quasiparticles, for which α/β=d, we obtain, from mu3 and mu4, inserting also the relation γ=2(z-1) valid for both, free (z=2) and Bogoliubov (z=1) quasiparticles, that


    -μ_3    = d - 2  ,
    
    -μ_4    = d - 4 +  z  .


Hence, we find, for z=2 in the large-N limit, that μ_3=μ_4, which implies that the relative importance of the scattering integrals I_3 and I_4 is not expected to change in time.






  §.§.§ Scaling function



A further important question concerns the form of the scaling function.
This form is determined by the stationary fixed-point equation, i.e., self-consistently, by the infrared scaling properties of the scattering integrals I_3 and I_4, for given fixed-point solutions f(𝐤,t_ref) at a fixed reference time t_ref.
Inserting s=(t/t_ref)^β into the scaling hypothesis distr_scaling yields the scaling form


    f (𝐤,t) = (t/t_ref)^α f_S([t/t_ref]^β𝐤) ,


with universal scaling function


    f_S(𝐤)=f (𝐤,t_ref).


For the spatial form of the scaling function we make the scaling ansatz


    f_S(𝐤)=s^κf_S(s𝐤)  ,


which we will motivate in more detail later on.
Inserting fScalingFunctionForm into the kinetic equation QKinEq yields the stationary fixed-point equation


    (α - βκ)f_S(p)    = t_refI[f_S](p).


If both sides of FPEqkappa are non-zero, they must scale in the same way in p.
Taking into account the fixed-time scaling


    I_l[f](𝐤,t_ref) = s^-μ_κ,l I_l[f](s 𝐤,t_ref).


of the scattering integral I_l defined by the scaling exponent μ_κ,l, one finds that the momentum exponent κ of the scaling function f_S  is 


    κ=-μ_κ,l ,


provided that the scaling of I_l determines that of I.

Hence, we need to determine the exponents μ_κ,l, which, in general, can take different values for the three-wave and the four-wave collision integrals.
As we will show in the following, κ>0 such that, in the scaling limit of small momenta, all distribution functions in the scattering integrals can be assumed to be in the semi-classical regime where f_p≫1, etc., for the other momenta appearing.
Hence, only the terms quadratic (cubic) in f can be assumed to contribute to I_3 (I_4).
As a result, power counting of I_3 and I_4 in this semi-classical regime, taking into account that the scaling 


    |T_l(𝐤_1,…,𝐤_l;t_ref)| 
      = s^-m_l |T_l(s 𝐤_1,…,s 𝐤_l;t_ref)| ,


of the T-matrices T_3 and T_4, respectively, is valid also at a fixed reference time, gives


    l = 3:        μ_κ,3 =    4+d+γ-2z-2κ ,
      
       l = 4:        μ_κ,4 =    8+2d+2γ-5z-3κ.


These scaling exponents, together with kappafrommukappa, yield the momentum exponent of the quasiparticle scaling function as


    l = 3:        κ=κ_3 =    4+d+γ-2z ,
      
       l = 4:        κ=κ_4 =    4+d+γ-5z/2 ,


given that either of the scattering integrals with l=3,4 dominates in the region of momenta considered.

For a given κ_l, and with the large-N exponents z=2 and γ=2 inserted, one finds that


    μ_κ,3-μ_κ,4 = κ_l- d≥1.


As a result, the four-wave scattering integral is expected to dominate at small momenta, k→0, over the three-wave term, such that


    κ=κ_4 = d+1     (z=2)


results as the momentum scaling exponent of the quasiparticle distribution f(𝐤,t)∼ k^-κ at the non-thermal fixed point.

This  appears to contradict the previous analysis of spatio-temporal scaling [cf. Eqs. Iltimemomentumscaling and mu3QP and mu4QP], which, for z=2 showed equal importance of I_3 and I_4 while, for z=1, the integral I_3 dominated the late-time scaling in the IR.
We emphasize, however, that the spatio-temporal scaling exponents betaQ3 and betaQ4 were obtained from scaling relations which are independent of the precise form of the scaling function f(𝐤,t) but only require the scaling relation fScalingForm. 
Hence, the question as to which of the vertices dominates the transport can be independently answered from the question of which vertex is responsible for the shape of the scaling function.

In contrast, in this section, we argue that the scaling function f(𝐤,t_ref) at a fixed reference time, in the IR scaling limit, i.e., for k_Λ(t_ref)→0, is the solution of a fixed-point equation dominated by I_4(𝐤,t_ref).
We note that, at a given early time of the evolution, when k_Λ is still large and this spatial scaling limit has not been reached yet, our estimate may not apply to the evolution starting from a given initial condition.
We will discuss this in more detail when comparing with numerical results in the following section.







  §.§.§ Relation between phase and phase-angle scaling forms



The result (<ref>), as well as the temporal scaling exponents (<ref>) and (<ref>) derived above appear, remarkably, in full agreement with earlier findings from non-perturbatively resummed kinetic theory (see <cit.>). To show this rigorously, however, still requires a translation between the quasiparticle distribution f_a(𝐤,t) characterizing the phase-angle excitations considered here and the scaling of the particle number distribution n_a(𝐤,t) encoded in the fundamental Bose field,


    n_a(𝐤,t)
      =⟨Φ_a^†(𝐤,t)Φ_a(𝐤,t)⟩≃ρ^(0)_aℱ[C_a(𝐱-𝐲,t)](𝐤) ,


where the Fourier transform ℱ (cf. Appendix <ref>) is taken with respect to 𝐫=𝐱-𝐲 of the on average translation-invariant phase correlator


    C_a(𝐫,t)
         =⟨ e^-ıθ_a(𝐱+𝐫,t)e^ıθ_a(𝐱,t)⟩,


independent of 𝐱.

Defining the (single-source) Schwinger-Keldysh generating functional of correlation functions of the phase angle field θ_a,


    Z[J] =  ∫ D θ  e^ı S_eff[θ] + ı∫_x,𝒞 J_a(x)θ_a(x),
  

we can express the phase correlator as


    C_a(𝐫,t)
         = Z[ {J_a=[-δ(𝐱-𝐫)+δ(𝐱)]δ_𝒞(x_0-t),0}],


where we take into account that the generating functional, for J=0, evaluates to unity on the closed time path, Z[0]=1, and that all other J_b, b≠a, vanish. Using the standard representation of the generating functional in terms of connected Greens functions [Note that the summation goes from n=1, as opposed to the standard zero-T expansion, which reflects that W[0] ≡ 0 on the Schwinger-Keldysh contour.],


    Z[J] 
    = 
    exp{∑_n=1^∞∫_{ x_i }_i=1^n,𝒞ı^n/n!  G_θ^(n)(x_1,…,x_n) J(x_1) ⋯ J(x_n) },


with


    G_θ,a^(n)(𝐱_1,…,𝐱_n,t)
         = ⟨θ_a(𝐱_1,t)⋯θ_a(𝐱_n,t)⟩_c


being an equal-time n-point connected correlator of the phase-angle field θ_a, we can express the two-point phase correlator as


    C_a(𝐫,t)
      = exp{∑_n=1^∞∑_𝐱_i=0,𝐫
            c_nG_θ,a^(n)(𝐱_1,…,𝐱_n,t)} ,


where 


    c_n= ı^n/n!exp{∑_i=1^nıπδ_r,|𝐱_𝐢|}


contains a combinatorial factor ı^n/n! as well as a sign which depends on the coordinates of the cumulants.

The above result provides a strong constraint on the scaling of the cumulants of the phase angle given that the phase correlator C_a satisfies the scaling relation


    C_a(𝐫,t)
      = s^α/β-dC_a(s^-1𝐫,s^-1/βt) ,


which follows directly from the scaling of n_a(𝐤,t). It should be emphasized that the above scaling exponents α and β are those associated with the scaling of fundamental bosonic fields Φ, for which the particle number conservation is expected due to U(1) symmetry. Then, assuming the absence of accidental cancellations, the relation (<ref>) between the two-point correlator and the phase-angle connected correlation functions thus implies 


    exp{ G^(n)_θ,a   ({𝐱_i},t)}|_𝐱_i∈{0,𝐫}
       =
    s^-α/β + d.exp{ G^(n)_θ,a({ s^-1𝐱_i},s^-1/β t) }|_𝐱_i∈{0,𝐫},
 

which, after taking the logarithm and recalling that α = d β is satisfied, yields


    .G_θ,a^(n)({𝐱_i},t)|_𝐱_i∈{0,𝐫}
      = .G_θ,a^(n)({s^-1𝐱_i},s^-1/βt)|_𝐱_i∈{0,𝐫},


or, after taking s = t^β,


    .G_θ,a^(n)({𝐱_i},t)|_𝐱_i∈{0,𝐫}
      = G_S,θ,a^(n)([t/t_ref]^-β𝐫)


with β independent of the correlator order n.

To summarize, we have shown that if density fluctuations are suppressed and hence (<ref>) and (<ref>) are applicable, the two-point equal-time correlator of fundamental Bose fields can be expressed in terms of connected correlation functions of the phase field as (<ref>). This furthermore implies, recalling that the particle number associated with the fundamental field Φ_a is conserved, a strong constraint on the scaling of the phase field correlators, (<ref>). This proves that, under the above assumptions, the scaling exponents derived from the Luttinger-liquid-like low-energy EFT coincide with those of the fundamental bosonic particles.  

Finally, let us consider the leading-order term, n=2 [Recall that the expectation value of the phase field is zero, so that the n=1 term does not contribute.], in (<ref>) which corresponds to the Gaussian approximation, 


    C_a(𝐫,t)
      ≃exp{G^(2)_θ,a(𝐫,t) - G^(2)_θ,a(0,t)}.


If the Gaussian approximation is sufficient, which is going to be discussed below, we can infer the scaling of the phase correlator C_a from the scaling of the 2-point equal-time cumulant G^(2)_θ,a of the angle. According to (<ref>), (<ref>), and (<ref>), we can take the momentum-space correlator (<ref>) to have the form


    f_a(𝐤,t) 
      ∼ k_Λ(t) |𝐤|^-(d+1) ,


for momenta k_IR < |𝐤| < k_UV in the power-law region. Here, the IR scale is expected to be associated with k_Λ(t), while the UV scale reflects the presence of a thermal tail that is typically developed in quenched ultracold atomic systems. To see (<ref>) explicitly, we note that (<ref>) and (<ref>) require f(𝐤,t_ref) ∼ |𝐤|^-(d+1). On the other hand, (<ref>) implies f(𝐤,t) ∼ t^α (t^β |𝐤|)^-(d+1), which combined with the aforementioned constraint α = d β results in (<ref>), where k_Λ(t) ∼ t^-β.

A simple scaling analysis of the Fourier transform of (<ref>) in d dimensions gives the spatial two-point function G^(2)_θ,a(𝐫,t)=ℱ[f_a](𝐫,t) as 


    G^(2)_θ,a(𝐫,t)-C_θ,a(0,t)
      =-const.× k_Λ(t) |𝐫| ,


which is normalized to some constant at r=0, and we added a minus sign reflecting that the correlations are expected to fall off at larger distances. Combining these results, we find that the two-point correlation function of the fundamental field should have a pure exponential decay at intermediate distances, which is consistent with our numerical results (cf. Fig. <ref>c),


    C_a(𝐫,t)
      ≃exp{
            -k_Λ(t)  |𝐫|},
               (k_UV^-1≲ 
            |𝐫| ≲ k_IR^-1).


This corresponds to the angle-averaged spatial first-order coherence function g^(1)(𝐱-𝐲,t)=⟨Φ^†(𝐱,t)Φ(𝐲,t)⟩ of the Bose field Φ(𝐱,t)


    g^(1)_a(𝐫,t) = ρ_a^(0)exp{-k_Λ(t)   |𝐫|} ,   (k_UV^-1≲ 
       |𝐫| ≲ k_IR^-1),


with uniform particle density ρ^(0)_a and background phase θ^(0)_a.
Its Fourier transform yields the occupation number distribution of the bosons,


    n_a^G(𝐤,t)=C_dk_Λ(t)/[k_Λ(t)^2+|𝐤|^2]^ζ/2 ,    (k_IR≲ 
    |𝐤| ≲ k_UV)


with normalization constant C_d=ρ^(0)_aΓ([d+1]/2)/π^(d+1)/2 and momentum exponent ζ=d+1.

Precisely this scaling form has been assumed in the previous non-perturbative analysis   <cit.>, where the scaling at large momenta |𝐤|≫ k_Λ(t) was predicted to be given by ζ=d+1 for z=2 as recovered in the Gaussian approximation assumed here. Furthermore, n_a(𝐤,t) was assumed to satisfy a scaling relation of the form (<ref>), 


    n_a(𝐤,t) = s^α/β n_a(s 𝐤,s^-1/βt) ,


such that the time-dependent constant scales as k_Λ(t)∼ t^-β.
The  exponents were predicted to be β=1/z=1/2 and α=β d, assuming z=2 and taking into account the conservation of the momentum integral over n_a(𝐤,t), cf. <cit.>.

Hence, provided the validity of the Gaussian approximation (<ref>), the scaling properties derived here in the large-N limit, defined by the exponents z=2, α=d/2, β=1/2, and κ=d+1, as well as by the scaling form (<ref>), are consistent with the scaling properties derived within the non-perturbative approach of <cit.>, for the case of a dynamical exponent z=2, and a vanishing anomalous dimension η=0.

As we will discuss in more detail in the following, the Gaussian approximation, for the situations studied here, is supported by the scaling analysis of the effective action. One expects this limit to be reached, ideally, for infinite system size in position as well as momentum space, at t→∞. In this limit, the interaction part of the action vanishes as compared to the part quadratic in the phase angle fields, such that also any correlations are expected to become Gaussian. This will, however, not be the general case. 
As we will argue, for sufficiently small β, the respective fixed point is expected to have a non-Gaussian character.

In summary, we observe that the values of the exponents obtained above are in a agreement with the results of <cit.>, under the assumption of a vanishing anomalous scaling, η = 0. 
With this, we conclude that, if 
(a) the kinetic description is adequate;
(b) the T-matrix has the form (<ref>) and (<ref>) for the three- and four-point interactions, respectively;
(c) the quasiparticle number or energy of the system are conserved in the scaling regime;
then our low-energy EFT predicts the universal self-similar dynamics close to the non-thermal fixed point in agreement with the previously developed non-perturbative approach in the special limit of infinite N.





  §.§.§ The case N=1 of a single field component



One can also apply our scaling analysis to the limiting case N=1 in which the action, in the next-to-leading order of a derivative expansion, reads as


    S_eff[θ] 
    =
    1/2 g   ∫_x,𝒞{θ(x) (-∂^2_t + c_s^2 ∇^2 ) θ(x) 
       +
    1/m∂_t θ(x) (∇θ (x) )^2 + 1/4 m^2(∇θ(x))^4 } ,


with c_s^2 = ρ^(0) g/m being the speed of sound. 
It describes, in the IR limit, the scattering of modes with linear Bogoliubov dispersion, with z=1 and γ=0.
Equivalently, the tensor couplings in Eqs. SeffnG3thetaa and SeffnG4thetaa reduce to
g^-1_1/N(𝐤) [1 - (1 + 𝐤^2/2 k^2_Ξ)^-1]=g^-1(1 + 𝐤^2/2 k^2_Ξ)^-1, i.e., at low energies k≪ k_Ξ, the coupling g_1/N(𝐤) in Seff4 is replaced by the constant bare coupling g. 

Given the scaling behavior of the dispersion and the coupling, the scaling analysis of the kinetic equation proceeds, for N=1, as described above in the large-N limit. 
In particular, the scaling relations between the different exponents, for general z and γ, are equally valid in both cases.
Inserting γ=0 and z = 1 into (<ref>) and (<ref>), one finds for quasiparticle transport the exponents 


    l = 3:   β = 1/2 ,
    
       l = 4:   β = 1/3 ,


and, from (<ref>) and (<ref>), that -μ_3=d-2>-μ_4=d-3. Hence, with increasing time, the scattering integral I_4, for z<2, loses out against I_3 such that the transport is predicted to occur, in the scaling limit, with β=1/2 according to (<ref>). As a result, despite z=1, the same values α=d/2 and β=1/2 describe transport of quasiparticles towards the IR, dominated by one-to-two scattering.

Determining the spatial exponent κ by using the same arguments as for the case z=2, which imply that  the 4-wave interactions dominate, one finds κ=κ_4=d+3/2. Note, however, that according to the arguments leading to the dominance of I_3 in the long-time limit, one rather expects, in the case N=1, the three-wave interaction to also dominate the stationary fixed-point equation. Under these circumstances, we rather predict


    κ=κ_3 = d+2     (z=1)


to represent the momentum scaling exponent in the large-time scaling limit t →∞. It should be remarked that this exponent describes the scaling of a vortex-dominated state, for d≥2, as found in many simulations for the case N=1 (cf., e.g., <cit.>).
This is rather unexpected as the Luttinger-liquid-based kinetic approach chosen here does not take into account the compact property of the phase and thus the excitations of vortices. The question as to whether the above result is just a coincidence or consequence of some physics thus remains open.




  §.§.§ Upper “critical” dimension



Given the above predictions for the scaling behavior at non-thermal fixed points we can use standard power-counting arguments to obtain the respective upper “critical” dimension. With this, we refer to the spatial dimension above which non-Gaussian terms in the action scale to zero at the non-thermal fixed point.

Note, however, that there is no critical point reached in the conventional sense: the scaling limit is expected to be attained, in the thermodynamic limit, at infinite evolution times. In practice, however, scaling already shows up at intermediate times, while the approach of the precise scaling function defining the scaling limit takes asymptotically long times. This can be conjectured from the results of our simulations shown in Fig. <ref>, in particular also from the observation that the higher-order powers in k_Λ(t)r contributing to the scaling form of g^(1)(r,t) take much longer to form than the low-order ones (see <cit.>).

To obtain an estimate for the upper “critical” dimension, we analyze the relative scaling of the quadratic and higher-order contributions to the effective action (<ref>). Considering the Schwinger-Keldysh action S(t)=∫_𝒞,t_ref^tdt'L(t') integrating from a reference time t_ref to the present evolution time t and back, its scaling is defined as


    S(s^-1/β t) = s^d_S S(t),


with canonical scaling dimension d_S=[S] of the action. Using the dynamical canonical scaling dimension of the field θ_a at the non-thermal fixed point,


    [θ_a] = -α/2β=-d/2 ,


we obtain (see Appendix <ref> for details) the canonical scaling of the quadratic, cubic and quartic parts of the action (<ref>),



    [S^(2)]    = 2z-γ-1/β   =    0 ,
      
    
      2[S^(3)]    = d+4+2z-2γ-2/β   =    d+4-2z ,
      
    
      [S^(4)]    = d+4-γ-1/β   =    d+4-2z ,



where we considered twice the scaling dimension of the 3-vertex as it occurs in even multiples in any diagram contributing to the two-particle irreducible (2PI) effective action, and inserted, in the respective second equations, γ=2z-2 and 1/β=2.

For the non-thermal fixed point to be a Gaussian fixed point in the IR scaling limit, the relative scaling conditions


    2[S^(3)]>[S^(2)] ,    
      [S^(4)]>[S^(2)]


need to be fulfilled.
This is the case in any dimension d>0, as d+4-2z>0 for z<(d+4)/2.

A remark is in order. The above analysis is to be understood, in general, as giving a first estimate as, strictly speaking, time-dependent correlation functions  need to be analyzed, e.g., within a functional renormalization-group approach <cit.> by studying the flows of Γ^(n). Nonetheless, for the fixed point considered here, the above arguments are unlikely to fail for our fixed point in the large-N limit since 1/β=z and, thus, any scaling in time is governed by the same exponent. The Gaussianity of the fixed point is further corroborated by the fact that also the scattering described by the integrals I_l[f] scales to zero in the infinite-time limit. This can be inferred from their scaling exponents μ_l defined in (<ref>) and (<ref>), since, for z=2, μ_3=μ_4=2-d, such that -μ_l-α/β=-2<0 and thus I_l[f](𝐤,t)=t^-2I_l[f]([t/t_0]^1/2𝐤,t_0). One should, however, be careful since that may rather indicate a critical slow down.  
Furthermore, we emphasize that Gaussianity of the non-thermal fixed point here refers to phase quasiparticles only, while in terms of the fundamental fields the fixed point can easily appear to be non-Gaussian.

In this context, we furthermore point out that the scaling of the T-matrix with m_l=2 is crucial for the validity of the perturbative Boltzmann approach.
As laid out in detail in Ref. <cit.>, the overoccupied momentum distributions would lead to a divergence of the scattering integral in the low-momentum limit and thus to a breakdown of the perturbative kinetic Boltzmann description of the transport process. 
In the context of wave turbulence, this problem is well known to lead to what is called strong wave turbulence <cit.>, requiring a modified kinetic equation.
Within the formulation in fundamental fields, a non-perturbative approximation of the kinetic equation is required to remedy this problem for the present case of a non-thermal fixed point, implying a T-matrix which scales to zero at low-momentum transfers <cit.>.
Here, however, the T-matrix scales to zero at low-momentum transfers already within the perturbative Boltzmann approximation, implying a Gaussian non-thermal fixed point and thus validity of the approach also in the scaling limit. 

Given the applicability of the above arguments, comparing (<ref>) and (<ref>), we find that the 4-vertex is irrelevant as compared to the free action for z<(d+4)/2, independent of the value of β. Leaving β arbitrary, we find that the condition for the non-thermal fixed point to be non-Gaussian in d≤ d_up gives the upper critical dimension,   


    d_up=2z-6+β^-1 ,


which is determined solely by the 3-vertex. 
To obtain an upper critical dimension d_up≥1 requires


    β^-1≥ d_up+6-2z.


For example, considering the scaling exponent β^-1≃5 found in <cit.> in d=2 dimensions and assuming the above general arguments to be applicable there, the value of β  implies d_up=2z-1≥1 for z≥1.
Hence, for the fixed point to be a non-Gaussian one in d=2 dimensions, the relevant z is required to be equal to or exceed z=3/2.  

We finally remark that one can analyze, in an analogous manner, the stationary spatial scaling of the cubic and quartic parts of the Lagrangian as compared to the Gaussian part.
This shows as well that, for the cases considered here, the upper critical dimension vanishes or is negative.
Hence, we expect that in the scaling limit, the Gaussian approximation considered in the previous sections is valid. 
Nevertheless, the interaction part of the action is crucial for the universal scaling transport to occur at all. 
In the course of the rescaling, leading closer to the fixed point, these interactions become increasingly suppressed, vanishing asymptotically at the fixed point.




 §.§ Concluding remarks on the scaling dynamics


To close the topic, let us briefly summarize the results obtained above and compare them with already known ones. 

For the model ONGPH, the scaling exponent β has previously been predicted, within a large-N approximation, by means of a scaling analysis of non-perturbative wave-Boltzmann-type kinetic equations governing the time evolution of n_a(𝐤,t) <cit.>. As the scaling evolution is subject to a conservation law, α and β are not independent. In the large-N limit studied there, the scaling dynamics corresponds to transport in momentum space which conserves the density, ρ_a=∫_𝐤n_a(𝐤,t)≡ const. [At finite times, on the way towards the fixed point, the integral only includes the long-wave-length modes up to a scale k_Λ(t) evolving with the distribution as k_Λ(t)∼ t^-β and is expected to be approximately conserved. Furthermore, as we showed, the transport shifts relative-phase excitations to the IR which only in the N→∞ limit decouple from each other and correspond to single-particle excitations in each component.], i.e., one finds α=β d in d dimensions. As β is positive, the transport is directed toward the infrared.

Analyzing this transport to lower k by means of kinetic equations describing the four-wave scattering between free quasiparticle modes, it was found that β=1/z <cit.>. This result is independent of  d and N and related to the dynamical scaling exponent z governing the dispersion ω_qp(𝐤)∼ k^z of the quasiparticles entering the kinetic description of the transport.

Our low-energy effective theory derived above provides an alternative, perturbative way of predicting scaling evolution of the type NTFPscaling, which is complementary to the non-perturbative approach chosen in <cit.>.
According to our theory, the occupation numbers at  k≲ k_Ξ are  dominated by the phase excitations n_a(𝐤,t)≃ρ^(0)_aℱ[⟨ e^-ıθ_a(𝐱,t)e^ıθ_a(𝐲,t)⟩](𝐤), where ℱ denotes the Fourier transform with respect to 𝐱-𝐲 of the on average translation invariant phase correlator (cf. Appendix <ref>). We find analytically that the dynamical exponent z, in the large-N limit, is z=2, corresponding to the Goldstone dispersion GoldstoneFreq and confirming the numerical results of, e.g., <cit.>.



The effective action (<ref>) includes non-linear interaction terms ∼θ^l, l=3,4, between the phase modes. These interactions lead to the non-linear transport describing scaling evolution of the type NTFPscaling. This can be shown in terms of kinetic wave-Boltzmann equations governing the phase excitation spectra ⟨θ_a(𝐤,t) θ_a(-𝐤,t)⟩, in analogy to <cit.>, cf. Kinetic. As we outline in more detail in ScalingLEEFT, for the four-wave interaction term, which dominates the spatial scaling in the IR limit, the scaling analysis of our kinetic integrals in a large-N approximation yields β=1/2 and ζ=d+1.
Hence, the infrared scale is predicted to algebraically decrease in time, k_Λ(t)∼ t^-β.

We have, moreover, studied the case of N=1. We find, in particular, the same value for the exponent β=1/2 is predicted despite the different dynamical exponent z=1 of the sound-wave dispersion, linearly scaling at low momenta, see SingleCompScaling. The cases 1≤ N<∞ are less straightforward and will be the subject of a forthcoming publication. 

Our analytical results are, within errors, consistent with the truncated Wigner simulations shown in NTFPEvolution. This is seen when taking into account that in leading approximation, the phase and phase-angle correlators of the translationally invariant system are related by ⟨ e^-ıθ_a(𝐱,t)e^ıθ_a(𝐲,t)⟩≃exp{⟨θ_a(𝐱,t) θ_a(𝐲,t)⟩}.
From this relation, the phase-angle correlator is approximately given by


    ⟨θ_a(𝐫+𝐱,t) θ_a(𝐱,t)⟩≈ -k_Λ(t) |𝐫| + 𝒪(|𝐫|^2).


If this is valid in d dimensions, up to additive terms involving cutoffs, the phase correlator scales as ⟨θ_a(𝐤,t) θ_a(-𝐤,t)⟩∼ |𝐤|^-(d+1), as analytically predicted within our low-energy effective theory (ζ=d+1).

We find that the consistency between our perturbatively obtained scaling exponent and the non-perturbative results of <cit.> crucially relies on the momentum scaling of the effective coupling app:g1N which is shown here to be related to the scaling ∼𝐤^2 of the hydrodynamic kinetic energy.
In fact, g_1/N depends only on momentum, density, and particle mass, but not on the microscopic coupling g. These same properties were found for the universal effective coupling function which enters the non-perturbative kinetic scattering integrals derived previously <cit.>.

We stress that, whereas in <cit.> 
scaling evolution after quenches of the gap and interaction parameters from a thermal initial state have been discussed in the context of a large-N approximation applied to a scalar O(N) model, as well as perturbative, ε-expansion, and functional renormalization group, the resulting scaling forms do not allow for the scaling to be discussed here.
While the cited work accounts for scaling in the context of critical coarsening, initial-slip dynamics and ageing phenomena <cit.>, for the type of quenches considered there, the same (initial-slip) scaling exponent results for the correlation and response functions <cit.>.
It is an interesting question beyond the scope of this work as to how initial-slip scaling manifests in the context of a non-thermal fixed point.


We finally emphasize that the non-thermal fixed point identified within the low-energy effective theory derived here turns out to have a Gaussian character with respect to the relevant degrees of freedom. As we discussed in more detail in UpperCritDim, the contributions to the action which are cubic and quartic in θ, exhibit subdominant scaling as compared to the Gaussian part and thus rescale to zero in the scaling limit.
One thus obtains, by standard scaling arguments for the bare action,  an upper `critical' dimension which is negative for the non-thermal fixed point with β=1/2, dominated by the sound and relative-phase excitations.
We anticipate that fixed points with smaller β, as observed in <cit.>, involving the interaction of non-linear excitations and topological defects, may have a positive upper `critical' dimension and thus have a non-Gaussian character
[The scaling limit of the non-thermal fixed point is reached at late evolution times provided that the initial condition is tuned such that no competing processes can prevent this.
For a given system size this is typically ensured by a sufficiently extreme out-of-equilibrium initial state, such as the box distribution used here.]. 
























§ NUMERICAL SIMULATIONS FOR N=3 AND D=3



In the remainder of this article we complement our numerical results presented in Ref. <cit.> and briefly in Intro with further details and observables.
We have solved the coupled Gross-Pitaevskii equations resulting from the Hamiltonian ONGPH for N=3, in d=3 spatial dimensions, by means of a spectral split-step algorithm.
We compute the time evolution of the correlation functions within the semi-classical truncated Wigner approximation <cit.>.

We start with a far-from-equilibrium initial condition given by a “box” momentum distribution n_a(𝐤,t_0) = n_0 Θ(k_q - |𝐤|), where, in the thermodynamic limit one can express the occupation number as  n_0=(4π k_q^3)^-1ρ^(0), and all phases θ_a(𝐤,t_0) are chosen randomly on the circle.
As before, a = 1,2,3 enumerates the three components. 
We choose the momentum cutoff for the initial condition to be k_q = 1.4 k_Ξ leading to an initial occupation number of n_0≃ 2350.
All simulations are performed on a grid with N_g = 256^3 grid points using periodic boundary conditions. 
The corresponding physical volume of our system is V = N_gΞ^3. 
The total particle number is 𝒩 =ρ^(0)V= 6.7 × 10^9, i.e., we have 𝒩_a = 2.23 × 10^9 particles in each of the three components. 
The correlation functions are averaged over 144 trajectories.
Varying the UV cutoff of the grid does not change our numerical results.

The initial condition in our simulations is chosen to take an extreme out-of-equilibrium form.
The modes are strongly occupied, on the order of n_0∼ζ_d^-3/2. 
Here, ζ_d=aρ^(0)^1/3 is the diluteness parameter which measures the scattering length a in units of the inter-particle spacing.
The occupancy distribution is cut off at a maximum scale k_q which is on the order of the relevant (healing) length scale set by the interactions, i.e., of k_Ξ=(8π aρ^(0))^1/2.
In physical terms, this means that all particles are placed, predominantly, at a momentum scale with kinetic energy  equal to the interaction energy or chemical potential μ=gρ^(0) the gas has stored if it was fully Bose condensed.
The phases of the excited modes are random, i.e., no significant coherence prevails.

In practice, one may prepare such an initial state, e.g., by applying a strong cooling quench to an equilibrated system or by allowing the system to pass through an instability after a parameter quench <cit.>.

During the intermediate to long-time evolution, the system can approach a non-thermal fixed point and obey the scaling evolution NTFPscaling within a range of modes 0<k≲ k_Ξ.

In NTFPEvolutiona, we depict the evolution of the occupation numbers of our system, starting from the above-described  initial condition (gray box).
Figure <ref>b shows that, for times t≳ t_ref=31 t_Ξ, the momentum distribution, to a good approximation, undergoes a scaling evolution according to NTFPscaling. 
Rescaling the distributions at different times during this period they collapse to a single universal scaling function. 

We have determined the respective scaling exponents α and β resulting from least-square rescaling fits of the occupation spectra within a time window [t_ref, t_ref + Δ t] with Δ t = 146 t_Ξ. 
The resulting t_ref-dependent exponents are shown in WindowFit1, confirming, within errors, the relation α=β d in d=3 dimensions.
We find similar results (not shown) for the other components.
For all three components, the exponent β is slightly larger than the value β=1/2  analytically predicted in the large-N limit and for N=1.

Analyzing the time evolution of the zero-momentum mode occupation n_1 (k=0,t) gives direct access to the scaling exponent α as the universal dynamics according to NTFPscaling reduces to n_1(k=0,t) = (t/t_ref)^α f_S,1(k=0).
Figure <ref> shows that at late times, t ≳ 200 t_Ξ, the evolution is  governed by α≃ d/2 = 3/2 consistent with our analytic prediction in the case of a conserved quasiparticle number.

In NTFPEvolutionc, we show the first-order spatial coherence function g^(1)_1(𝐫,t)=⟨Φ_1^†(𝐱+𝐫,t)Φ_1(𝐱,t)⟩ for different times.
They take, to a good approximation, an exponential form at short distances r ≳Ξ.
Similar results are found (not shown) for the components a=2,3. 
Such an exponential form of the coherence function is reminiscent of that of a quasicondensate in an equilibrium gas in one spatial dimension <cit.>.
Hence, the coherence function signals the buildup of a non-equilibrium quasicondensate in three spatial dimensions, rescaling in time and space towards a longer-range coherence.
Further studies of the correlation properties of this state concerning its relation to an equilibrium quasicondensate are beyond the scope of this work and are to be presented elsewhere.

In the following, we want to briefly comment on the form of the scaling function of the momentum distribution n_a(k).
Therefore, we compare the momentum distribution, at different times, with two different scaling functions (see InverseOcc).
The first one, quoted in app:NTFPscalingnkG, is obtained within the Gaussian approximation PhaseCorrCumGauss of the relation between the phase-angle and the phase correlators.
It corresponds to the purely exponential first-order coherence function  app:NTFPscalingg1.
The second scaling function used takes the form 


    n_a(𝐤,t)=C_dk_Λ(t)[4k_Λ^d+1(t)+|𝐤|^d+1]^-1.


Taking the Fourier transform of NTFPscalingnk, for d=3, gives 

    g^(1)(𝐫,t) ≈ρ^(0)e^-k_Λ(t) |𝐫| sinc(k_Λ(t) |𝐫|) ,


which has an oscillatory contribution due to the cardinal-sine function.
We find that at late times, our data are approximatively described by the scaling function NTFPscalingnk.
The scaling form app:NTFPscalingnkG, however, differs from the data for all evolution times.
Hence, the corresponding first-order coherence function has an additional oscillatory contribution, which becomes visible at larger distances.
See <cit.> for more details.
Accounting for the oscillatory part within our analytical treatment requires a more refined analysis of non-linear phase excitations which is beyond the scope of this work.

We  remark that the scaling form NTFPscalingnk is only one exemplary choice which takes into account the most striking features of n_a(k).
To capture all details of the data, the scaling form might involve additional terms. 
To determine the exact scaling function is beyond the scope of this work. 













We furthermore analyze the scaling evolution of the momentum distribution C_12(k,t) = ⟨|Φ_1^†Φ_2(k,t)| ^2⟩ and the corresponding spatial coherence function g_12(r,t) of the relative phase operator which we show in RelativePhase12.
The rescaled data collapse in a similar manner as the single-component occupation numbers, with extracted exponents presented in WindowFit12.
The scaling function g^(1)_12 in position space (RelativePhase12c) shows a form reminiscent of an exponential decay comparable to the first-order coherence function.


We finally analyze, in HydroDecompEvo, the time evolution of the total momentum distribution n_tot(k,t)=∑_a=1^3n_a(k,t) in hydrodynamic decomposition.
This is obtained by decomposing the kinetic energy distribution in momentum space into quantum-pressure (q), spin (s), nematic (n), incompressible (i) and compressible (c) parts, as defined in detail in HydroDecomp.
The decomposition provides additional information on the character of the hydrodynamic flows corresponding to the phase field θ_a.

We point to the observation that the incompressible flow arising from vortical motion is subdominant as compared to the nematic and spin parts which determine the leading scaling of the occupation number as n_tot(k)∼ k^-4.
The vortex flow contributes significantly only at very small momenta, where the plateau appears in the total number distribution at very low momenta. 
In this regime, as is seen at times t/t_Ξ=31, 61, the nematic and spin parts fall off towards zero momentum but this is compensated for by the incompressible part.
Note that the distributions shown in HydroDecompEvo result from the tensor decomposition of the current which is a four-point function of the Bose fields. 
Hence,  n_kin(k,t) which is the sum of the hydrodynamic parts, deviates from n_tot(k,t) in the infrared <cit.>. 

We emphasize that the spin and nematic fluctuations are determined by the fluctuations of the relative phases between the field components a=1,2,3.
The corresponding momentum distributions take the same form as the scaling function shown in RelativePhase12b. The universal scaling of C_12(k,t) confirms our hypothesis that the scaling behavior is dominated, for non-zero N, by the relative phases of the components. 

We additionally find that the incompressible as well as the spin and nematic parts of the hydrodynamic decomposition also exhibit universal scaling according to NTFPscaling (see HydroDecompRescaled). 
The late-time value of β≃ 0.5 extracted for the universal scaling of the spin and nematic parts (see insets HydroDecompRescaledb and c) demonstrates the crucial role of the relative phases to the scaling behavior of the system.
Interestingly the scaling behavior of the subdominant vortical flow in the system seems to be characterized by a scaling exponent β≃ 0.58 nearly unaffected by variations of the reference time (see inset HydroDecompRescaleda).
For smaller N, in particular N=1, vortices are expected to play a more prominent role <cit.>  in the evolution following a quench of the type considered here
[The transport associated with a non-thermal fixed point in the N=1 case has recently been associated with a `transport peak' occurring in the spectral decomposition of the statistical correlation function, independent of the thermal Bogoliubov excitations in the system <cit.>.].


We finally mention that during the approach of the scaling limit and thus of the non-thermal fixed point, the system shows prescaling  <cit.> (see also <cit.> in a perturbative high-energy context).
Our analysis demonstrates that (cf. in particular HydroDecompEvo) effects of non-linear excitations such as vortices, contributing to the incompressible flow, can induce scaling violations of the single-component occupation numbers and coherence functions.
On the other hand, also other contributions such as 1/N corrections to the scaling analysis could explain certain systematic deviations of the exponents found at the largest simulation times from the predictions presented in this work.
A more detailed finite-N analysis is the subject of future work.







§ CONCLUSIONS



We present a low-energy effective theory for the interacting phase-angle excitations of a U(N)-symmetric Gross-Pitaevskii model of N Bose fields with local quartic self-coupling of the total particle density. The theory provides a perturbative formulation of far-from-equilibrium low-energy universal scaling dynamics at a non-thermal fixed point. This is complementary to the non-perturbative approach on the basis of fundamental Bose fields chosen previously, while being technically easier to evaluate in the scaling regime.  

Our approach provides the leading-order dynamical scaling exponent z=2 of the quasiparticle eigenfrequencies in the large-N limit.
This result for z closes a longstanding gap in the non-perturbative kinetic-theory formulation of non-thermal fixed points <cit.>.
Applying the theory in the large-N limit, we recover the universal scaling exponents at a non-thermal fixed point predicted previously within the non-perturbative approach (cf. <cit.>).

We find analytically that, in leading order, the first-order coherence function, close to a non-thermal fixed point, falls off exponentially in space, g^(1)(𝐫,t) ∼ e^-k_Λ(t) |𝐫|, with a coherence-length scale rescaling as k_Λ(t)∼ t^-β in time, with β=1/2 both, for the case N=1 and N→∞.
This is reminiscent of equilibrium systems in d<3 dimensions where such an exponentially reduced long-range phase order indicates the presence of a quasicondensate, with the there static coherence length depending on the coupling and/or temperature.
Also in these situations, without spontaneous symmetry breaking and a field expectation value singling out a zero-momentum condensate mode, the states are characterized by strong occupancies of the low-energy momentum modes of the system.

Our analytical predictions are corroborated by the results of Truncated-Wigner simulations for N=3 in d=3 dimensions.
Considering g^(1)(r,t) at short distances r=|𝐫| where it falls off exponentially, we are able to confirm, with high accuracy, the analytic predictions β=1/2 and α=d/2 for the spatio-temporal scaling exponents, effectively leaving little space for anomalous deviations.

This finding is consistent with our analytical result that the corresponding non-thermal fixed point, which is approached in the scaling limit of infinite evolution times, has a Gaussian character as we argue by standard scaling arguments applied to the bare action of the low-energy effective theory.
In contrast, to obtain a positive upper “critical” dimension a much smaller exponent β is required as, e.g., was found numerically for an anomalous fixed point dominated by vortex interactions in two spatial dimensions <cit.>.


We emphasize that the numerically determined coherence function, however, shows oscillations at larger distances which are not covered by the analytical approach which rests on a homogeneous background phase θ^(0)_a.
An extension of the theory to a non-uniform θ^(0)_a(𝐫,t) seems viable but is beyond the scope of this work.

When evaluating our low-energy effective description for the case of a single-component Bose gas, N=1, we find the same spatio-temporal scaling exponent β=1/2, despite the fact that the dynamical exponent is z=1 in this case.
This analytical result corroborates earlier numerical findings presented in <cit.>.


Our results support the conjecture, that related universal scaling in O(N)-symmetric relativistic models (cf. <cit.>), while the O(N) symmetry is broken by the evolution, is connected with the appearance of an approximately conserved charge due to a remaining U(1) symmetry not broken by the flow <cit.>.
This was also seen in numerical simulations <cit.>.
In consequence, we presume that β=1/2 would apply also for these systems, independent of N, reflecting  relative phase fluctuations between the field components and their universal transport toward low momenta. 
As the non-thermal fixed-point scaling relation α=dβ, to a good approximation, holds also for the relative-phase correlators C_ab(k,t) shown in RelativePhase12 as well as for the spin-contributions to the energy spectrum seen in HydroDecompRescaled, a further emerging symmetry is expected to play an important role.
This symmetry has been conjectured to be related to the suppression of density fluctuations in the system <cit.> and will be further discussed elsewhere.

From our analytical results in the large-N limit and for the case N=1 together with our numerical results for N=3, we propose that the universality class corresponding to the scaling exponent β is related to the dynamical breaking of a U(1) symmetry and therefore independent of N.
Our approach offers itself for a refinement using field-theoretic renormalization-group techniques and for applications to small-N spin systems available in experiment.

In the following Appendices, we provide further details of
the calculations leading to our results presented in the main
text.






§ ACKNOWLEDGEMENTS
 
The authors thank I. Aliaga Sirvent, J. Berges, K. Boguslavski, R. Bück­er, I. Chantesana, S. Diehl, S. Erne, F. Essler, S. Heupts, M. Karl, P. Kunkel, S. Lannig, D. Linnemann, A. Mazeliauskas, M. K. Oberthaler, J. M. Pawlowski, A. Piñeiro Orioli, M. Prüfer, R. F. Rosa-Medina Pimentel, J. Schmiedmayer, T. Schröder, H. Strobel, and C. Wetterich for discussions and collaboration on the topics described here. 
This work was supported by the Horizon-2020 framework programme of the European Union (FET-Proactive, AQuS, No. 640800, ERC Advanced Grant EntangleGen, Project-ID 694561),  by Deutsche Forschungsgemeinschaft (SFB 1225 ISOQUANT), by Deutscher Akademischer Austauschdienst (No. 57381316), and by Center for Quantum Dynamics, Heidelberg University.
C.-M.S. thanks the Dodd-Walls Centre for Photonic and Quantum Technologies, University of Otago, Dunedin, New Zealand, for hospitality and support. 
T.G. thanks the Erwin Schrödinger International Institute for Mathematics and Physics, Wien, Austria, for hospitality and support within their topical program Quantum Paths.






























§ NOTATION


In this paper, we adopt the (-,-,-,+) signature convention for the metric so that the scalar product of two (d + 1)-vectors x = (x_1,…,x_d,x_0) = (𝐱,x_0) and y = (y_1,…,y_d,y_0) is given by x y = x_0 y_0 - 𝐱𝐲. Define the Fourier transform as:


    ℱ[f(x)](p)    ≡ f(p) = ∫_x e^-ı p x f(x) , 
    ℱ^-1[f(p)](x)    ≡ f(x) = ∫_p e^ı p x f(p) ,


with the short-hand notations ∫_x ≡∫d^d+1 x and ∫_p ≡∫d^d+1 p/(2 π)^d+1.
We also introduce the Fourier transform over only space (time) variables:


    f(𝐤,t)    = ∫_𝐱 e^ı𝐤 𝐱 f(𝐱,t) = ∫_ω e^ıω t f(𝐤,ω) ,
    
    f(𝐱,ω)    = ∫_t e^-ıω t f(𝐱,t) = ∫_𝐤 e^-ı𝐤 𝐱 f(𝐤,ω).








§ LOW-ENERGY EFFECTIVE ACTION OUT OF EQUILIBRIUM



In many cases, the description of a system of interest can be significantly simplified if one only considers the physics above some length scale 1/Λ, or, equivalently, below a momentum and thus energy scale Λ. 
This typically is the case if there is a separation of scales and the dynamics distinguishes fields with different masses (“heavy” vs. “light” field modes) or if some degrees of freedom are suppressed compared to the others (“fast” vs. “slow” modes). 
A low-energy effective field theory is then obtained by integrating out the fast (i.e., heavy) fluctuations. 
The matching condition reads as


    ρ (t) 𝒪(t) = ρ_eff (t) 𝒪(t),


where 𝒪(t) is an operator measuring only slow (light) degrees of freedom.

In the following, we will discuss the evaluation of the above expectation value within the Schwinger-Keldysh approach and the Feynman-Vernon influence-functional approach <cit.> (for details, see also <cit.>).
Assuming that the initial density matrix ρ̂(t_0) at a distant past time t_0 splits into phase and density fluctuations, and that the density fluctuation part can be taken in Gaussian approximation, it becomes clear that the computation of the low-energy effective action can be done as in a standard zero-temperature approach.
The reason behind this is that, when moving towards a non-thermal fixed point in a closed system, the system loses its memory about the details of the initial state, similarly as driven-dissipative <cit.> and equilibrium systems. 
As a consequence the above assumptions are expected to be well-justified.

Let us consider the example where φ (x) describes the slow modes and ψ (x) the fast ones with the following action:


    S[φ,ψ] = S_φ[φ] + S_ψ[ψ] + S_int[φ,ψ] ,


where S_φ[φ] and S_ψ[ψ] are the (semi-)classical actions that do not contain any mixing terms and thus describe independent dynamics of the fields, while S_int[φ,ψ] corresponds to interaction between two fields. 

The total density matrix of the whole system is given by


    ρ[φ,ψ;t] = ⟨φ^+,ψ^+| ρ̂(t)  | φ^-,ψ^-⟩,


where ± denotes that the field is defined on the positive or negative branch of the closed time path, respectively. According to the above, the reduced (or effective) density matrix, which is the object of interest, is then defined by tracing out the fast modes at the time t considered,


    ρ_eff[φ;t] = ∫ d ψ ⟨φ^+,ψ | ρ̂(t)  | φ^-,ψ⟩ .


Quite often, the initial state can be chosen to be described by a Gaussian density matrix. Higher-order correlations in this case are built up during the evolution. Similarly, hereafter we assume that the initial density matrix is factorized, i.e.,


    ρ̂(t_0) = ρ̂_φ(t_0) ×ρ̂_ψ(t_0) ,


where ρ̂_φ(t_0) and ρ̂_ψ(t_0) are the initial density operators of the fields φ and ψ, respectively. The condition (<ref>) implies the absence of correlation between φ and ψ at t = t_0. Nevertheless, due to the presence of the interaction term S_int[φ,ψ], the latter will be generated at later times t > t_0. In our case, φ→θ and ψ→δρ so that (<ref>) means that there is no correlation between phase and density fields in the initial configuration. In the absence of topological defects, this is typically an adequate assumption.

The evolution of the reduced matrix ρ_eff[φ,t] then reads as


    ρ_eff[φ̅;t] = ∫ [d φ_0^+] [d φ_0^-]  (φ̅;t  | φ_0;t_0) ρ_eff[φ_0;t_0] ,


where we use the bar to label the value at the time t. The transition amplitude can be expressed in terms of a closed-time path functional integral


    (φ̅;t  | φ_0;t_0) = ∫_φ_0^+^φ̅^+𝒟' φ^+∫_φ_0^-^φ̅^-𝒟' φ^-  e^ı S_eff[φ] ,


with


    S_eff[φ] = S_φ[φ^+] - S_φ[φ^-] + S_IF[φ^+,φ^-]


being the full effective action while S_IF is the influence action. The latter is given by


    𝔉[φ^+,φ^-] ≡ e^ı S_IF[φ^+,φ^-] = ∫ [d ψ̅] ∫ [d ψ_0^+] [d ψ_0^-]  ρ_ψ [ψ_0;t_0] 
    ×∫_ψ_0^+^ψ̅𝒟' ψ^+∫_ψ_0^-^ψ̅𝒟' ψ^-  e^ı{ S_ψ[ψ^+] + S_int[φ^+,ψ^+] - S_ψ[ψ^-] - S_int[φ^-,ψ^-] } ,


where we have introduced the so-called Feynman-Vernon influence functional 𝔉 <cit.>, which can be conveniently rewritten as


    𝔉[φ^+,φ^-] = ⟨𝒯_𝒞 e^ı S_int[φ,ψ]⟩_ψ ,


with average being defined with respect to the ψ-field density matrix. Note that the Schwinger-Keldysh contour can be depicted as shown in Fig. <ref>. 



Now, given the effective action S_eff [Eqs. (<ref>), (<ref>)], we can define the effective generating functional Z_eff[J,R;ρ̂_eff,0]:




    Z_eff[J,R;ρ̂_eff,0] 
       = 
    ∫ [d φ_0^+] [d φ_0^-]  ρ_eff[φ_0;t_0] ∫_φ_0^+^φ_0^-𝒟' φ  e^ı{ S_eff[φ] +  ∫_x,𝒞φ (x) J(x) + 1/2∫_xy,𝒞φ(x) R(x,y) φ(y) }
       =
    ∫ [d φ_0^+] [d φ_0^-]  ρ_eff[φ_0;t_0] ∫_φ_0^+^φ_0^-𝒟' φ  e^ı{ S_φ[φ^+] - S_φ[φ^-] + S_IF[φ^+,φ^-] +  ∫_x,𝒞φ (x) J(x) + 1/2∫_xy,𝒞φ(x) R(x,y) φ(y) }
       =
    ∫ [d φ_0 d ψ_0]  ρ[φ_0,ψ_0;t_0] ∫_φ_0^+,ψ_0^+^φ_0^-,ψ_0^-𝒟' φ𝒟' ψ  e^ı{ S_full[φ,ψ] +  ∫_x,𝒞φ (x) J(x) + 1/2∫_xy,𝒞φ(x) R(x,y) φ(y) } ,


where we have introduced 


    S_full[φ,ψ] 
    ≡(S_φ[φ^+] + S_ψ[ψ^+] + S_int[φ^+,ψ^+]) -
    (S_φ[φ^-] + S_ψ[ψ^-] + S_int[φ^-,ψ^-]) ,


and the condition ψ̅^+ = ψ̅^- = ψ̅ is implied. In the last line, we also used the following property:


    ρ_eff[φ_0;t_0] 
    = 
    ∫ d ψ ⟨φ_0^+,ψ | ρ̂(t_0)  | φ_0^-,ψ⟩ 
    = 
    ρ_φ[φ_0;t_0] ∫ d ψ ⟨ψ  |  ρ̂_ψ(t_0)   |  ψ⟩ 
    =
    ρ_φ[φ_0;t_0] ρ̂_ψ(t_0) = ρ_φ[φ_0;t_0].


By comparing the first line of (<ref>) with the last one we get the final expression for the effective action:


    exp{ı S_eff[φ]}   =
    exp{ı(S_φ[φ^+] - S_φ[φ^-]) }∫ [dψ_0^+] [dψ_0^-]  ρ_ψ[ψ_0;t_0] 
       ×∫_ψ_0^+^ψ_0^-𝒟' ψ exp{ı(S_ψ[ψ^+] + S_int[φ^+,ψ^+] - S_ψ[ψ^-] - S_int[φ^-,ψ^-]) } ,




with condition ψ̅^+ = ψ̅^- = ψ̅ being, again, implied.
As a result, assuming that the initial density matrix of the heavy modes ρ_ψ[ψ_0;t_0] takes the form of a Gaussian ground states, the integral in seff_final falls into two independent factors, the positive and negative closed-time path zero-temperature contributions. 
It follows that the influence functional 𝔉 is formally equivalent to the Keldysh vacuum generating functional, and the influence action infpathint is the vacuum effective action. 
Since in this work we concentrate on dynamics at the non-thermal fixed point, which is formally reached only at t →∞, we can use the memory loss property of non-thermal fixed points, as well as the fact that density fluctuations are suppressed, to choose the reference time t_0 such that the density matrix of the corresponding density fluctuations can be taken in Gaussian approximation in the ground state, thus reducing the computation of the effective action to the zero-temperature case, as pointed out above.





§ EIGENVALUES OF THE MASS MATRIX



In this Appendix, we determine the eigenvalues of the mass matrix


    M^ab = (𝐤^2/2 m)^2 δ^ab + 2 g ρ^(0)_n 𝐤^2/2mδ^nb,
   

where no assumptions on the values of ρ^(0)_b were made.
We define a = (𝐤^2/2 m)^2, b_n = 2 g ρ^(0)_n (𝐤^2/2m) such that this matrix can be rewritten as


    M = a   I_N × N + J_N × N 𝒫 = [ a + b_1     b_2       ⋯     b_N;     b_1 a + b_2       ⋯     b_N;       ⋮       ⋮       ⋱       ⋮;     b_1     b_2       ⋯ a + b_N ].


Here, I_N × N is the unit matrix, J_N × N is the matrix of ones, and 𝒫 = (b_1,…,b_N ).

We now take a look at the second term. This matrix contains N identical rows, and hence it has the eigenvalue λ = 0 with an algebraic multiplicity N - 1. 
It is easy to check that the corresponding eigenvectors are given by 


    𝐞_1   = ( -1, 1, 0,...,0 )^T ,   𝐞_2 = ( -1, 0, 1,...,0 )^T , …,
    𝐞_N-1   = ( -1, 0, 0,...,1 )^T ,


meaning that the geometric multiplicity is also N - 1. 
The remaining eigenvalue,


    J_N × N 𝒫 𝐞_N = (∑_i b_i )𝐞_N ,
 

corresponds to the eigenvector


    𝐞_N = ( 1, 1,…, 1 )^T.


Hence, we conclude that M possesses the eigenvalue λ = a with multiplicity N - 1 and λ = a + ∑_i b_i with multiplicity 1. 
We thus have two types of excitation modes, with 
eigenfrequencies being the square roots of these eigenvalues,


    ω_1(𝐤) 
         = ω_2 (𝐤) = … = ω_N-1 (𝐤) = 𝐤^2/2m ,
      
    ω_N (𝐤) 
         = √(𝐤^2/2m(𝐤^2/2m + 2 g ρ^(0))).







§ SPONTANEOUS SYMMETRY BREAKING



In this appendix we discuss the breaking of the U(N) symmetry by the particular choice of initial conditions we apply. 
In the respective ground state, this corresponds to spontaneous symmetry breaking. 
Consider the N-component Gross-Pitaevskii Lagrangian in the absence of external potential:


    ℒ =     ı/2[φ^*_a(x) ∂_t φ_a(x) - φ_a(x) ∂_t φ^*_a(x)]
       - 1/2m[∇φ_a^*(x)] ·[∇φ_a(x)] - g/2[φ_a(x) φ_a(x)]^2 ,


On a classical level, a mean-field expectation value of φ can be derived from the classical potential:


    𝒱(|φ_a|) =  g/2(∑_b |φ_a|^2)^2,


which yields a zero-valued ground state. This is not surprising since (<ref>) corresponds to a vacuum QFT. In condensed matter physics, one is typically interested in field theories with non-zero background densities. To obtain a finite background density, one introduces Lagrange multipliers associated with the U(1) charges for each component,


    j_a^0 
      = ∂ℒ/∂(∂_0 φ_a)δφ_a/δα + c.c.
      = φ_a^*φ_a=ρ_a


where ∂_0≡∂_t, and the aforementioned U(1) symmetries read 


    φ_a → e^-ıαφ_a ,   φ^*_a → e^ıαφ^*_a.


Introducing the chemical potentials one replaces


    H → H - μ_a N_a,


with particle number N_a=ρ_aV in volume V, and thus, in non-relativistic case,


    ℒ→ℒ + μ_a φ_a^*φ_a.


We note that this, in general, explicitly breaks the U(N) symmetry into a product U(1)^× N. The classical potential then takes the form


    𝒱(|φ_a|) = -μ_a |φ_a|^2 + g/2(∑_b |φ_a|^2)^2.


Varying with respect to, e.g., φ_b^* one gets the conditions


    (μ_b -  g ∑_a |φ_a|^2 ) φ_b = 0,    b = 1,...,N.


It follows directly that, in order to have a non-vanishing amplitude φ_b of component b, the associated chemical potential must be μ_b=gρ, with total density ρ=∑_aρ_a. Hence, the chemical potentials for all non-vanishing densities need to be equal and given by the total density. Therefore, in the ground state, the densities will be non-zero in the subspace having identical chemical potentials with the largest value, while the rest are zero. These non-vanishing ground-state densities ρ_a^(0) add up to the total density ρ^(0)=μ/g set by the chemical potential μ for the respective components.

In order to have a non-zero total density without specifying densities of each component and thus remain a U(N) symmetry, a chemical-potential term is introduced as


    H → H - μ N ⟷ℒ→ℒ + μ φ_a^*φ_a ,


which is equivalent to fixing only the total density. We note that the above is equivalent to replacing, with μ = g ρ^(0), 


    -g/2(∑_b φ_b^*φ_b)^2 → 
    -g/2(∑_b φ_b^*φ_b - ρ^(0))^2.


Let us close the discussion with the following remark. Above we tried to preserve a full U(N) symmetry. On the other hand, this symmetry might be explicitly broken in experimentally relevant scenarios. Indeed, in the experiment, particle number is typically specified for each component. One can take care of this by choosing an appropriate initial density matrix ρ̂(t_0) which fixes the expectation values of the particle number operators at t = t_0. In this paper, we disregard this subtlety and only take it into account by choosing the vacuum with ρ_a^(0)≠ 0 for each component. 





§ INVERSE OF THE KERNEL G̃^AB



To derive the inverse of g̃^ab defined in app:G_inv, we consider a general N × N matrix of the form U(u,v) = u I + 𝐯⊗𝐯 , where 𝐯 = (v_1,v_2,…,v_N)^T and ⊗ denotes an outer (tensor) product. To find the inverse of U(u,v) we use the ansatz U^-1(u,v) = x I - 𝐲⊗𝐲 , which yields the condition


    x u δ_ij - u  y_i y_j + x  v_i v_j - (𝐯·𝐲)  v_i y_j = δ_ij ,


where we used (𝐯⊗𝐯) · (𝐲⊗𝐲) = (𝐯·𝐲) 𝐯⊗𝐲.
We note that, if 𝐯 = 0, then obviously also 𝐲 = 0 such that we get x = 1/u and the above condition becomes


    y_i y_j - v_i/uv_j/u + (𝐯·𝐲) v_i/u y_j = 0.
 

This is a system of non-linear equations, which, in principle, cannot be solved analytically. However, we note that the above system is almost symmetric in 𝐲↔𝐯/u, which gives a hint to use the following ansatz as a solution: 𝐲 = c 𝐯 , where c > 0 is some positive constant. Using this we obtain


    (c^2 - 1/u^2 + c^2 ∑_j v_j^2/u) v_i v_j = 0 ,


such that


    c^2 = 1/u (u +∑_j v_j^2 )→ y_i = v_i/√(u (u + ∑_j v_j^2 )).


According to app:G_mod we can use these results with u = (N g /2)(𝐤^2/ k^2_Ξ) and v_a = (N gρ_a^(0)/ρ^(0))^1/2, a=1,…,N, to obtain the inverse of g̃^ab as given in app:G_inv.







§ REARRANGEMENT OF 3- AND 4-VERTICES



The three- and four-wave coupling terms SeffnG3thetaa and SeffnG4thetaa, respectively, can be rewritten using


    ∫_𝐤,𝐤'𝐤' · (𝐤'-𝐤)/g_1/N(𝐤)∂_t θ_a(-𝐤,t) θ_a (𝐤',t) θ_a (𝐤 - 𝐤',t) 
       =-∫_𝐤_2,𝐤_3𝐤_2 · (𝐤_2 + 𝐤_3)/g_1/N (-𝐤_3)∂_t θ_a (𝐤_3,t) θ_a (𝐤_2,t) θ_a (-𝐤_2 - 𝐤_3,t) 
       =
    -∫_𝐤_1,𝐤_2,𝐤_3𝐤_1 ·𝐤_2/g_1/N(𝐤_3)δ (𝐤_1 + 𝐤_2 + 𝐤_3) θ_a (𝐤_1, t) θ_a (𝐤_2,t) ∂_t θ_a (𝐤_3,t)  ,


and


    ∫_𝐤,𝐤',𝐤”   g_1/N^-1 (𝐤) 𝐤' · (𝐤 + 𝐤') [𝐤”· (-𝐤 + 𝐤”)]θ_a(𝐤') θ_a(-𝐤 - 𝐤') θ_a(𝐤”) θ_a(𝐤 -𝐤”) 
       =- ∫_𝐤_1, 𝐤_2, 𝐤_3g_1/N^-1 (𝐤_1 - 𝐤_2)(𝐤_1 ·𝐤_2) [𝐤_3 · (𝐤_1 + 𝐤_2 + 𝐤_3))]θ_a(𝐤_1) θ_a(𝐤_2) θ_a(𝐤_3) θ_a(-𝐤_1 - 𝐤_2 - 𝐤_3) 
      
       = - ∫_𝐤_1,𝐤_2, 𝐤_3, 𝐤_4g_1/N^-1 (𝐤_1 - 𝐤_2)(𝐤_1 ·𝐤_2) (𝐤_3 ·𝐤_4)θ_a(𝐤_1) θ_a(𝐤_2) 
      θ_a(𝐤_3) θ_a(𝐤_4) δ (𝐤_1 + 𝐤_2 + 𝐤_3 + 𝐤_4) ,


where we also used that g_1/N (𝐤) = g_1/N (-𝐤).
Taking the above terms together gives the effective action Seff4.






§ NON-LINEAR SIGMA MODEL



In Sect. <ref> in the main text, the derivation of the low-energy effective action required suppression of the density fluctuations for each component. While possible, such an assumption is definitely not always satisfied. For that reason, we outline, in this appendix, following <cit.>, the derivation of the effective theory which is independent of the aforementioned assumption and only requires the suppression of the total density fluctuations.




  §.§.§ Field parametrization


To begin with, let us recall the basic statement of the (Lorentz-invariant) Goldstone theorem. Suppose the Lagrangian is left invariant by a symmetry group G with n(G) generators, whereas the vacuum is only invariant under a subgroup H with n(H) generators. Then, there are n(G) - n(H) Nambu-Goldstone gapless bosons. In other words, Goldstone excitations live in the coset space G/H. In our case, G is U(N) and H is U(N-1). Having the above in mind, we write down the following mathematical statement:


    U(N)/U(N-1) = S^2N-1.


In principle, we can now use the standard field parametrization procedure as described in, e.g., <cit.>. However, since our theory is not Lorentz-invariant, the Goldstone theorem cannot be applied directly <cit.>. Furthermore, recalling that the Bogoliubov excitation is associated with a global phase, it is tempting to set it apart directly. Since, a phase can always be seen as a point on a circle S^1, we would like to, in some sense, divide S^2N-1 into an S^1 piece and the rest. Using the mathematical language, this can be done by a Hopf-like fibration. One considers the S^2N-1 sphere as an S^1 bundle over the complex projective space ℂℙ^N-1. 

The elements of the projective space ℂℙ^N-1 are the complex planes through the origin. For a complex plane span_ℂ(w=(Z_1,… Z_N)^T)⊆ℝ^2N all points are associated with each other and represented by a representative geometry. In principle, the choice of the representative is completely arbitrary, though usually it is possible to give a construction that relates this abstract space to geometric shapes. To do so, we normalize w: w' = w/|w|. We note that w^' now lies on the S^2N-1 sphere. The intersection of the complex plane span_ℂ(w) and S^2N-1 defines a circle of points S^1. Any point on that circle can now be chosen as the representative. Here, we just take all points of the circle to be associated with one point of the complex projective space. This defines the following projection:


    π:    S^2N-1→ S^2N-1/S^1≅ℂℙ^N-1,    w' → w'/S^1.


This projection cannot be inverted globally but only locally. S^2N-1 can be understood as the set of representatives called projective space with a circle S^1 attached to each point. Mathematically, this, together with some requirements on the smoothness of the attachment, defines a fiber bundle on the ℂℙ^N-1 with fiber S^1. Compared to the product ℂℙ^N-1× S^1 the above fiber bundle has to be seen as a collection of circles (S^1_z)_z∈ℂℙ^N-1, where copies are distinguished by the base point z. This type of decomposition of a sphere is also called Hopf fibration. 

The fiber bundle should locally look like a product. Therefore, around the ground state we are allowed to write the point of the sphere as a product,


    ψ = e^ıθχ,   with  ψ∈ S^2N-1,  e^iθ∈ S^1,  χ∈ℂℙ^N-1.


Next, we parametrize the complex projective space using homogeneous coordinates. Choosing some point (1,z), z ∈ℂ^N-1, one can find a representative of the complex plane spanned by this point,


    χ(z) = 1/√(1+z ^†z)[ 1; z ] , 
      with  z ∈ℂ^N-1 ,


which lies on a unit sphere. We are now able to assemble our field parametrization:


    φ = √(ρ)ψ = √(ρ) e ^ıθ1/√(1+z ^†z)[ 1; z ] ,    z ∈ℂ^N-1 ,


with the ground state chosen φ^(0) = (√(μ/g),0,…,0)^T, for the sake of simplicity. 





  §.§.§ Classical picture


Plugging the parametrization (<ref>) into the Gross-Pitaevskii Lagrangian, one can then find that the Lagrangian now takes the following form:


    ℒ
    =
       - 
    ρ∂_t θ + ı/2ρ/1 + z^† z(z^†∂_t z - h.c.) 
       - 
    1/2m{(∇√(ρ))^2 + 
    ρ [ (∇θ +ı/2(∇ z^†)z- h.c/(1+z^† z)^2)^2 
       +
    (∇ z^†)(∇ z)/1+z^† z - (∇ z^†)z· z^† (∇ z)/(1+z^† z)^2]} - g/2ρ^2 + μρ.


We first analyze the classical picture behind (<ref>). To that end, we expand the Lagrangian up to the second order in fields and use μ = g ρ^(0). The result reads:


    ℒ^(2) 
    = 
       -δρ∂_t θ  +i/2ρ^(0)(z^†∂_t z - h.c.)- g/2δρ^2 
       -
    1/2m[1/4ρ^(0)(∇δρ)^2 +ρ^(0)(∇θ)^2 +ρ^(0)∇ z ∇ z^†].


The corresponding Euler-Lagrange equations take the following form:



    ∂_t δρ   = 
    - ρ^(0)/m∇^2 θ ,
    ∂_t θ   =
    1/4 m ρ^(0)∇^2 δρ - g δρ ,
    ı∂_t z 
       =
    -ρ^(0)/2m∇^2 z ,



and the same (complex conjugate) equation for z^†. Combining (<ref>) with (<ref>), we get:



    ∂_t^2 θ   = ∇^2/2m(-∇^2/2m  + 2 g ρ^(0)) θ ,
    ı∂_t z    = -ρ^(0)/2m∇^2 z.



With this, one obtains a global phase excitation with Bogoliubov dispersion and N-1 coupled (particle-hole) z-excitations with free-particle dispersion. 




  §.§.§ Effective action


Following the same procedure as in luttingerleeft we now want to integrate out heavy degrees of freedom. In this case, those are the longitudinal density fluctuations δρ since they correspond to positive eigenmodes of the potential's curvature. By integrating them out, one can find, to the leading order in derivative expansion, the following effective Lagrangian <cit.>:


    ℒ_eff 
    =
       ρ^(0)(-∂_tθ + ı/2z^†∂_tz - ∂_t z^† z/1 + z^† z) 
       -
    ρ^(0)/2m(∇θ + ı/2z^†∇ z - ∇ z^† z/1 + z^† z)^2 
       -
    ρ^(0)/2m( ∇ z^†∇ z/1 + z^† z - (∇ z^† z)(z^†∇ z)/(1+z^† z)^2) 
       -
    1/2g(-∂_tθ + ı/2z^†∂_t z - ∂_t z^† z/1 + z^† z)^2 + h.o.t.


We note that the resulting effective Lagrangian has a form of a (non-relativistic) non-linear sigma model (NLSM) with the Fubini-Study metric. As it was shown in <cit.>, different types of NLSM also appear in the experimentally relevant system of an ultracold Bose gas with spin-spin interactions. Furthermore, such models are known to enjoy stable topologically non-trivial configurations, and it was observed that dynamical properties of a non-thermal fixed point can change dramatically in presence of topological defects (vortices) in a Bose gas <cit.>. 
Finally, to close the discussion, let us mention that out-of-equilibrium dynamics of non-linear sigma models has been considered in the context of disordered metals and superconductors <cit.>. The discussion of far-from-equilibrium dynamics within such non-linear sigma models is beyond the scope of this work and of interest for future research.





§ KINETIC EQUATION



In this Appendix, we discuss the derivation of the quantum Boltzmann equation for the Luttinger-liquid-like effective action (<ref>) following the standard routine (see, e.g., <cit.>). We start with the Dyson equation written in the following form:


    ∫_z,𝒞    G_0^-1(x,z) G(z,y) - ∫_z,𝒞[Σ(x,z) + ı R(x,z) ] G(z,y)
    
    =     δ_𝒞(x-y).


In the absence of explicit external sources, the term R(x,y) has a support only for x^0 > t_0 reflecting the initial density matrix contribution. The initial condition for G(x,y) has to be specified at x^0,y^0 = t_0, while the subsequent evolution for times t > t_0 follows from (<ref>) with R = 0. Therefore, we can set R(x,y) to zero at t>0.
See EFT_general for a more detailed discussion of the relevance of the initial state within our low-energy effective field theory description. 



Recalling that θ(x) has a zero macroscopic value, we note the classical inverse propagator G_0(x,y) coincides with the free propagator (<ref>), which in real space reads as [We changed the notation of the mean density to n_a^(0) for this section to avoid a possible confusion with the spectral function.]


    ı G^-1_0(x,y; ⟨θ⟩ = 0) = -𝒦(𝐱 - 𝐲) ∂_x^0^2 δ_𝒞(x^0 - y^0) + n_a^(0)/m∇^2_𝐱δ_𝒞 (x - y) ,


with 


    𝒦(𝐱 - 𝐲) = ∫_𝐤e^-ı𝐤 (𝐱 - 𝐲)/g_G(𝐤).


Equation (<ref>) then takes the form


    ∂_x^0^2 ∫_𝐳𝒦(𝐱 - 𝐳) G(𝐳,x^0;y) - n_a^(0)/m∇^2_𝐱 G(x,y) + 
    ı∫_z,𝒞Σ(x,z) G(z,y) = -ıδ_𝒞(x-y).


Next, we decompose the propagator G(x,y) and the proper self-energy Σ(x,y) into spectral and statistical parts:



    G(x,y)    = F(x,y) - ı/2ρ(x,y) _𝒞(x^0 - y^0),
    Σ(x,y)    = -ıΣ^(0)(x) δ(x-y) + Σ^F(x,y) - ı/2Σ^ρ(x,y) _𝒞(x^0 - y^0),



where we also introduced the local part of the self-energy Σ^(0)(x). Let us analyze each term in (<ref>) separately. The first term takes the form


    ∫_𝐳   𝒦(𝐱 - 𝐳) ∂^2_x^0[F(𝐳,x^0;y) - ı/2ρ(𝐳,x^0;y) _𝒞(x^0 - y^0) ] 
    
    =    ∫_𝐳𝒦(𝐱 - 𝐳) [∂_x^0^2 F(𝐳,x^0;y)  
    - ı/2_𝒞(x^0 - y^0) ∂_x^0^2 ρ(𝐳,x^0;y) - 2 ıδ_𝒞(x^0 - y^0) ∂_x^0ρ(𝐳,x^0;y) - ıρ(𝐳,x^0;y) ∂_x^0δ_𝒞(x^0 - y^0) ]  
    
    =    ∫_𝐳𝒦(𝐱 - 𝐳) [∂_x^0^2 F(𝐳,x^0;y) - ı/2∂_x^0^2 ρ(𝐳,x^0;y) _𝒞(x^0 - y^0) ] - ıδ_𝒞(x - y),


where we have used that [The following statements should, of course, be understood in terms of distributional derivatives]:


    ∂_x^0_𝒞(x^0 - y^0) = 2 δ_𝒞(x^0 - y^0),   ρ(𝐳,x^0;y) ∂_x^0δ_𝒞(x^0 - y^0) = -δ_𝒞(x^0 - y^0) ρ(𝐳,x^0;y),


and


    δ_𝒞(x^0    - y^0) ∫_𝐳𝒦(𝐱 - 𝐳) ∂_x^0ρ(𝐳,x^0;y) = ıδ_𝒞(x^0 - y^0) ⟨[∫_𝐳𝒦(𝐱 - 𝐳) ∂_x^0θ(𝐳,x^0), θ(𝐲,y^0)] ⟩
       =
    - ıδ_𝒞(x^0 - y^0) ⟨.[θ(𝐲,y^0),π_θ(𝐱,x^0) ] |_x^0=y^0⟩ = δ_𝒞 (x^0 - y^0) δ(𝐱 - 𝐲) = δ_𝒞(x - y).


Since the second term and the local part of the third one in (<ref>) are trivial, it is only left to consider the non-local part of the latter:


    ı   ∫_z,𝒞(Σ^F(x,z) - ı/2Σ^ρ(x,z) _𝒞(x^0 - z^0) ) (F(z,y) - ı/2ρ(z,y) _𝒞(z^0 - y^0) ) 
       = ∫_𝐳(∫_t_0^x^0 d z^0  Σ^ρ(x,z) F(z,y)
    - ı/2_𝒞(x^0 - y^0) ∫_y^0^x^0 d z^0  Σ^ρ(x,z) ρ(z,y)  
    - ∫_t_0^y^0 d z^0  Σ^F(x,z)ρ(z,y) ),


where we have used _𝒞(x) = θ_𝒞(x) - θ_𝒞(-x) and properties of the Schwinger-Keldysh contour. 

Collecting the results obtained above and comparing the real and imaginary parts of the resulting equation, we find the following set of equations for statistical and spectral functions:




    (𝒟̂_x F)(x,y)    = -∫_t_0^x^0 dz  Σ^ρ(x,z) F(z,y) + ∫_t_0^y^0 dz  Σ^F(x,z) ρ(z,y), 
    
    (𝒟̂_x ρ)(x,y)    = -∫_y^0^x^0 dz  Σ^ρ(x,z) ρ(z,y),



where we have adopted the short-hand notation ∫_t_1^t_2 d z = ∫_t_1^t_2 d z^0 ∫ d^d z. Above, we also introduced an integro-differential operator 𝒟̂_x related to the real-space Green's function and defined as


    (𝒟̂_x h)(x,y) = ∫_𝐳𝒦(𝐱 - 𝐳) ∂_x^0^2 h(𝐳,x^0;y) - n_a^(0)/m∇_𝐱^2 h(x,y) + Σ^(0)(x)h(x,y).
  

The system of coupled equations (<ref>) describing the evolution of a two-point Green's function is commonly referred as the Kadanoff-Baym equations. Note that in the absence of approximations on the self-energy, these equations are exact. These equations, however, are not feasible to solve analytically. Therefore, some further approximations are required, leading to Boltzmann-type kinetic equations. 
Note that, in this work, we intend to perform a scaling analysis of the system at non-thermal fixed points.
This can be achieved using the Kadanoff-Baym equations directly <cit.>, while, however, the scaling analysis of the kinetic equations <cit.> is more straightforward and employs standard procedures <cit.>.

To proceed, we perform the so-called gradient expansion <cit.> of the Kadanoff-Baym equations. To do so, we first interchange x and y in the equations and subtract the result from (<ref>) and (<ref>), respectively. Using the relation between retarded/advanced and spectral components one then gets:




    {(𝒟̂_x - 𝒟̂_y ) F }(x,y)    = ∫_zθ(z^0) (F(x,z) Σ^A(z,y) + G^R (x,z) Σ^F(z,y) - Σ^R(x,z)F(z,y) - Σ^F(x,z)G^A(z,y)),
    {(𝒟̂_x - 𝒟̂_y ) ρ}(x,y)    = ∫_z(G^R(x,z) Σ^ρ(z,y) + ρ(x,z) Σ^A(z,y) - Σ^ρ(x,z)G^A(z,y) - Σ^R(x,z)ρ(z,y)).



Next, we introduce center and relative coordinates:


    X_xy^μ = x^μ + y^μ/2,      s_xy^μ = x^μ - y^μ.


The idea is to approximate a finite time description of the system by employing the above equations with t_0 → -∞. This then allows for a systematic expansion with respect to derivatives of X^μ and powers of s^μ. The analysis of the right-hand side of (<ref>) is standard and can be found in, e.g., <cit.>, and we therefore concentrate on the left-hand side. From (<ref>) it immediately follows that



    ∂_x^0^2    = ∂_s^0_xy^2 + ∂_s^0_xy X^0_xy^2 + ∂_X^0_xy^2,
    ∂_y^0^2    = ∂_s^0_xy^2 - ∂_s^0_xy X^0_xy^2 + ∂_X^0_xy^2,



which yields


    ∫_𝐳{𝒦(𝐱 - 𝐳) ∂_x^0^2 F(𝐳,x^0;y) - 𝒦(𝐳 - 𝐲) ∂_y^0^2 F(x;𝐳,y^0) } 
    =
    ∫_𝐳   {∂^2_s_xy^0 X_xy^0[𝒦(𝐱 - 𝐳) F(𝐳,x^0;y) + 𝒦(𝐳 - 𝐲) F(x;𝐳,y^0)] 
       + 
    (∂_s_xy^0^2 + ∂_X_xy^0^2 ) [𝒦(𝐱 - 𝐳) F(𝐳,x^0;y) 
    - 𝒦(𝐳 - 𝐲) F(x;𝐳,y^0)] } .


Next, following the idea of the Wigner distribution function, we perform a Fourier transform with respect to relative coordinates, i.e., we define:


    F(𝐤,ω,X)    = ∫_-2 X^0^2 X^0 d s^0   e^-ıω s^0∫_-∞^∞ d^ds   e^ı𝐤𝐬 F (X + s/2,X - s/2),
    ρ̃(𝐤,ω,X)    = -ı∫_-2 X^0^2 X^0 d s^0   e^-ıω s^0∫_-∞^∞ d^ds   e^ı𝐤𝐬ρ(X + s/2,X - s/2),


and similarly for Σ^F and Σ^ρ. Here, the imaginary unit ı factor was added to have ρ̃(X,𝐤,ω) real [Recall that F(x,y) is even, whereas ρ(x,y) is odd with respect to relative coordinates, which implies that the Fourier transform of the former is real, while the latter one is purely imaginary.]. Finally, while, for given finite X^0, relative coordinate s^0 is limited by ± X^0, it is standard to extend the integration limits in (<ref>) to ±∞, which is related to the assumption t → -∞ that has been done above. With all that in mind, we perform a Fourier transform of both sides of (<ref>) with respect to the relative coordinate s_xy. Again, we focus on (<ref>) since the rest is either similar or standard for any (scalar) theory. To that end, we first perform a spatial transformation, e.g.,


    ∫_𝐳∫_𝐬_xy    e^ı𝐤𝐬_xy𝒦(𝐬_xz) F (𝐗_zy + 𝐬_zy/2, X_xy^0 + s^0_xy/2; 𝐗_zy - 𝐬_zy/2, X_xy^0 - s^0_xy/2) 
    
    =
       ∫_𝐳∫_𝐬_xy∫_𝐩 e^ı𝐤𝐬_xy - ı𝐩𝐬_xz𝒦(𝐗_zy + 𝐬_zy/2) F (𝐩,X_xy^0 + s^0_xy/2; 𝐗_zy - 𝐬_zy/2, X_xy^0 - s^0_xy/2) 
    
    =
       ∫_𝐳∫_𝐬_xy∫_𝐩∫_𝐪 e^ı𝐤𝐬_xy - ı𝐩𝐱 + ı𝐪𝐲 + ı(𝐩 - 𝐪)𝐳𝒦(𝐩) F(𝐪,x^0,y^0) = ∫_𝐬_xy∫_𝐩 e^ı (𝐤 - 𝐩) 𝐬_xy𝒦(𝐩) F(𝐩,x^0,y^0) = 𝒦(𝐤) F(𝐤,x^0,y^0),


where 𝒦(𝐤) = (g_G(𝐤))^-1, which immediately follows from (<ref>). Here, we have already assumed no dependence on spatial central coordinates due to homogeneity of the system. We thus conclude that after the spatial Wigner transform the RHS of (<ref>) becomes:


    (<ref>)→∂_s_xy^0 X_xy^0^2 [𝒦(𝐤) F(𝐤,x^0,y^0) + 𝒦(𝐤) F(𝐤,x^0,y^0)] + (∂_s_xy^0^2 + ∂_X_xy^0^2 ) [𝒦(𝐤) F(𝐤,x^0,y^0) - 𝒦(𝐤) F(𝐤,x^0,y^0) ] 
    
    =
    2 𝒦(𝐤) ∂_s_xy^0 X_xy^0^2 F(𝐤,x^0,y^0) .


As it was already mentioned, the remaining part of the analysis is standard and can be found in <cit.>. In short, we finally end up with the set of the so-called transport equations:




    2 k^0/g_G(𝐤)∂ F(k,X^0)/∂ X^0   = Σ̃^ρ(k,X^0) F(k,X^0) - Σ^F(k,X^0) ρ̃(k,X^0),
    2 k^0/g_G(𝐤)∂ρ̃ (k,X^0)/∂ X^0   = 0 ,



which are the leading order equations with respect to ∂_X^μ and k^μ expansion. Note that (<ref>) implies that the spectral remains unchanged during the evolution.

It remains to compute the self-energies Σ̃^ρ and Σ^F in some chosen approximation. In this paper, we restrict ourselves to only first non-trivial contributions to the self-energy, see Fig. <ref> for a diagrammatic representation. We discuss the adequacy of such an approximation a posteriori when discussing the properties of the scattering matrix, see ScalingLEEFT.

Before we proceed, let us define the real-space couplings:



    γ(𝐱_1,…,𝐱_3) 
       =
    ∫_𝐩_1,…,𝐩_3 e^ı(𝐩_1𝐱_1 + … + 𝐩_3 𝐱_3 )γ(𝐩_1,…,𝐩_3) (2 π)^d  δ(∑_i=1^3 𝐩_i ),
    λ(𝐱_1,…,𝐱_4)    =
    ∫_𝐩_1,…,𝐩_4 e^ı(𝐩_1 𝐱_1 + … + 𝐩_4 𝐱_4 )λ(𝐩_1,…,𝐩_4) (2 π)^d δ(∑_i=1^4𝐩_i ).



Note that the sign might seem to be off (c.f. Appendix <ref>), but it reflects the definition


    ∫_{𝐤_i } g_l(𝐤_1, …, 𝐤_l) θ(𝐤_1,t) …θ(𝐤_l,t) ≡∫_{𝐱_i } g_l(𝐱_1, …, 𝐱_l) θ(𝐱_1,t) …θ(𝐱_l,t).




We are now in the position to analyze the diagrams. First, we note that the left, so-called tadpole, diagram is typically local and only gives a shift of a mass. However, since the interactions of our low-energy effective action are non-local, one has to be a bit more careful. Because the interaction term is still local in time, so is the tadpole. Therefore, it only contributes to Σ^F and has only a central coordinate X^0 = t. This implies that in Wigner space, i.e., after taking a Fourier transform with respect to relative coordinates, it is proportional to δ (k^0), as one expects from a time-local tadpole. Lastly, we notice that it only enters the transport equations (<ref>) together with ρ̃(X^0,k), which is zero at k^0 = 0. Hence, we conclude that the tadpole diagram does not contribute to the kinetic equation.

Next, we consider the central (sunset) diagram. Its symbolic expression reads as


    Σ_2-loop(x,y) 
    =
    -∫_𝐱_1,𝐱_2,𝐱_3 
    𝐲_1,𝐲_2,𝐲_3λ(𝐱,𝐱_1,𝐱_2,𝐱_3) λ(𝐲, 𝐲_1,𝐲_2,𝐲_3) G(𝐱_1, t;𝐲_1, t') G(𝐱_2, t; 𝐲_2, t') G(𝐱_3, t;𝐲_3, t').


Introducing the short-hand notation, G(𝐱_i, t;𝐲_i, t') ≡ G_i, and employing the definition (<ref>), it is straightforward to show that


    G_1 G_2 G_3 = [F_1 F_2 F_3 - 1/4(F_1 ρ_2 ρ_3 + ρ_1 F_2 ρ_3 + ρ_1 ρ_2 F_3 )] + [F_1 F_2 ρ_3 + F_1 ρ_2 F_3 + 
     ρ_1 F_2 F_3 - 1/4ρ_1 ρ_2 ρ_3 ] (-ı/2_𝒞(t - t') ),


which combined with (<ref>) yields



    Σ_2-loop^F(x,y) 
       = 
    -∫_𝐱_1,𝐱_2,𝐱_3 
    𝐲_1,𝐲_2,𝐲_3λ(𝐱,𝐱_1,𝐱_2,𝐱_3) λ(𝐲, 𝐲_1,𝐲_2,𝐲_3) [F_1 F_2 F_3 - 1/4(F_1 ρ_2 ρ_3 + ρ_1 F_2 ρ_3 + ρ_1 ρ_2 F_3 )],
    Σ^ρ_2-loop(x,y) 
       =
    -∫_𝐱_1,𝐱_2,𝐱_3 
    𝐲_1,𝐲_2,𝐲_3λ(𝐱,𝐱_1,𝐱_2,𝐱_3) λ(𝐲, 𝐲_1,𝐲_2,𝐲_3) [F_1 F_2 ρ_3 + F_1 ρ_2 F_3 + ρ_1 F_2 F_3 - 1/4ρ_1 ρ_2 ρ_3 ].



Now we perform a Wigner transformation. Since all the terms have a similar structure, we only consider a single of them:


    ∫_{𝐱_i }, {𝐲_i }
     s_xyλ(𝐱,𝐱_1,𝐱_2,𝐱_3) λ(𝐲, 𝐲_1,𝐲_2,𝐲_3) ρ_1 F_2 ρ_3 e^-ı k s_xy 
    =
       ∫λ(𝐩,𝐩_1,𝐩_2,𝐩_3) λ(𝐪,𝐪_1,𝐪_2,𝐪_3) ρ_1 F_2 ρ_3   e^-ı k s_xy + ı (𝐩𝐱  + … + 𝐪_3 𝐲_3)
       ×
    (2 π)^2d δ( 𝐩 + 𝐩_1 + 𝐩_2 + 𝐩_3 ) δ(𝐪 + 𝐪_1 + 𝐪_2 + 𝐪_3).


Here, we suppressed the integration variables to ease the notation. Next, we notice that the integration over {𝐱_i } and {𝐲_i } can be transformed into the integration over {𝐬_xy,i} and {𝐗_xy,i} via a unitary transformation. Furthermore, we use that


    e^ı𝐩_i𝐱_i + ı𝐩_i𝐱_i = e^ı (𝐩_i + 𝐪_i) 𝐗_xy,i + ı(𝐩_i - 𝐪_i) 𝐬_xy,i/2,


which, recalling that, for a spatially homogeneous system, F_i and ρ_i do not depend on spatial central coordinates, yields a series of delta-functions (2π)^d δ(𝐩_i + 𝐪_i). We thus obtain:


    (<ref>)   =
    ∫λ(𝐩,𝐩_1,𝐩_2,𝐩_3) λ(𝐪,-𝐩_1,-𝐩_2,-𝐩_3)   e^ı (-k s_xy + 𝐩𝐱 + 𝐪𝐲 + 𝐩_1 𝐬_xy,1 + 𝐩_2 𝐬_xy,2 + 𝐩_3 𝐬_xy,3)ρ_1 F_2 ρ_3 (2 π)^2d δ( 𝐩 + 𝐩_1 + 𝐩_2 + 𝐩_3 ) 
       =×δ(𝐪 - 𝐩_1 - 𝐩_2 - 𝐩_3) 
    =
    ∫ |λ(𝐩,𝐩_1,𝐩_2,𝐩_3)|^2 
    e^ı (- k s_xy + 𝐩𝐬_xy + 𝐩_1 𝐬_xy,1 + 𝐩_2 𝐬_xy,2 + 𝐩_3 𝐬_xy,3)ρ_1 F_2 ρ_3 (2π)^d δ (𝐩 + 𝐩_1 + 𝐩_2 + 𝐩_3) 
       =
    ı^2 ∫ |λ(-𝐤,𝐩_1,𝐩_2,𝐩_3)|^2 e^ı(-k^0 + p_1^0 + p_2^0 + p_3^0) s_xy^0ρ(𝐩_1,p_1^0,X^0) F(𝐩_2,p_2^0,X^0) 
    ρ(𝐩_3,p_3^0,X^0) (2 π)^d  δ (-𝐤 + 𝐩_1 + 𝐩_2 + 𝐩_3) 
       =
    -∫_p,q,r |λ(-𝐤,𝐩,𝐪,𝐫)|^2 
    ρ(p,X^0) F(q,X^0) ρ(r,X^0) (2 π)^d+1 δ(k - p - q - r),
 

where in the last line we have changed the notation of dummy variables. Here, we adopted the notation |λ(𝐤,𝐩,𝐪,𝐫)|^2 ≡λ(𝐤,𝐩,𝐪,𝐫) λ(-𝐤,-𝐩,-𝐪,-𝐫). Repeating the same procedure for each term in (<ref>) and recalling that for Σ^ρ an additional factor -ı has to be added, we conclude that up to 2 loops the transport equations take the form:



    2 k^0/g_G(𝐤)∂F_k/∂ X^0 
    =   ∫_p,q,r |λ(-𝐤,𝐩,𝐪,𝐫)|^2 [ρ̃_k  F_p F_q F_r  + 1/4ρ̃_k (F_p ρ̃_q ρ̃_r + ρ̃_p F_q ρ̃_r + ρ̃_p ρ̃_q F_r ) 
       -
    F_k F_p F_q ρ̃_r  - F_k F_p ρ̃_q F_r - F_k ρ̃_p F_q F_r - 1/4 F_k ρ̃_pρ̃_q ρ̃_r ] (2 π)^d+1 δ (k - p - q - r),
    2 k^0/g_G(𝐤)∂ρ̃_k/∂ X^0 
    =    0,



where F_k ≡ F(k,X^0), ρ̃_k ≡ρ̃(k,X^0). We now write, without any loss of generality, the generalized relation


    F(k,X^0) = (f(k,X^0) + 1/2) ρ̃(k,X^0) ,


which in the equilibrium case reduces to the Callen-Welton fluctuation-dissipation theorem. Changing in (<ref>) the integration variable p → -p and employing that for a real scalar theory one has 


    F(-p,X) = F(p,X),   ρ̃(-p,X) = -ρ̃(p,X),


we obtain, after some simple algebra,


    2 k^0/g_G(𝐤)ρ̃_k  ∂_t f_k = ∫_p,q,r|λ (𝐤,𝐩,-𝐪,-𝐫) |^2 (2 π)^d+1 δ(k + p - q - r) [(f_k + 1)(f_p + 1) f_q f_r - f_k f_p (f_q + 1) (f_r + 1)  ] ρ̃_pρ̃_kρ̃_qρ̃_r,


where we have relabeled the central coordinate X^0 to t, reflecting the fact that it plays the role of time in our limit, and used |λ(-𝐤,-𝐩,𝐪,𝐫)|^2 ≡ |λ(𝐤,𝐩,-𝐪,-𝐫)|^2.

At last, we adopt the quasiparticle (or on-shell) approximation, i.e., assume that the spectral function takes the form of a free one, ρ̃(p) ≈ρ̃^(0)(p) [Note that, for a non-local theory as given by the vertices depicted in Fig. <ref>, the validity of the on-shell approximation needs to be verified
after all. This question is beyond the scope of the present work and is the subject of further study.]. From (<ref>) one can derive:


    G^R,(0)(k) = g_G(𝐤)/(ω + ı 0^+)^2 - ω_𝐤^2 = g_G(𝐤)/2 ω_𝐤(1/ω - ω_𝐤 + ı 0^+ - 1/ω + ω_𝐤 + ı 0^+),


with k = (𝐤, ω) and ω_𝐤 = 𝐤^2/2m. Recalling that


    lim_ϵ→ 01/x + ıϵ = -lim_ϵ→ 0ϵ/x^2 + ϵ^2 = -πδ(x)


and using the relation between spectral and retarded Green functions, one gets


    ρ̃(𝐤,ω) = π g_G(𝐤)/ω_𝐤[δ (ω - ω_𝐤) - δ (ω + ω_𝐤)].


It should be emphasized that naively it might seem that (<ref>) violates the sum rule. We note, however, that the prefactor suggests that the sum rule has to be enhanced with a prefactor 2 k^0/g_G(𝐤) prefactor, similar to relativistic systems. This is further corroborated by (<ref>) and the commutation relation of the phase field. 

Finally, we substitute this into (<ref>) and integrate over ∫_0^∞d k^0/2 π. The left-hand side then reads as


    LHS = ∂_t ∫_0^∞d k^0/2 π2 k^0/g_G(𝐤)π g_G(𝐤)/ω_𝐤[δ (ω - ω_𝐤) - δ (ω + ω_𝐤)] f(k,t) = ∂_t f(𝐤,t),


where we have defined the occupation number distribution function as


    f(𝐤,t) ≡ f(𝐤,ω_𝐤,t),


This definition is motivated by the equilibrium limit of (<ref>), for which f(X^0,k) becomes a thermal distribution function. The right-hand side is a bit more elaborate, yet very straightforward as well. First, we notice that, since integration over p, q, and r goes from -∞ to +∞, both δ-functions in (<ref>) contribute, and we thus have to deal with 2^3 = 8 different terms. In order to analyze them, it is convenient to split the frequency integrals, e.g., ∫_-∞^∞ d p^0 … =  ∫_-∞^0 d p^0 … + ∫_0^∞ d p^0 … with p → -p variable change for the negative part and then use f(𝐩,-p^0,t) = -[f(𝐩,p^0,t) + 1 ]. For completeness, let us consider one term explicitly:


    -   ∫_𝐩,𝐪,𝐫∫_0^∞ d k^0 ∫_0^∞d p^0/2π∫_-∞^0 d q^0/2π∫_0^∞d r^0/2π |λ(𝐤,𝐩,-𝐪,-𝐫)|^2 π^4 g_G(𝐤) g_G(𝐩) g_G(𝐪) g_G(𝐫)/ω_𝐤ω_𝐩ω_𝐪ω_𝐫δ(k^0 - ω_𝐤) δ(p^0 - ω_𝐩) 
       ×δ(q^0 + ω_𝐪) δ(r^0 - ω_𝐫) 
    (2 π)^d+1 δ (k + p - q - r) [(f_k + 1)(f_p + 1) f_q f_r - f_k f_p (f_q + 1) (f_r + 1) ] 
       = 
    ∫_𝐩,𝐪,𝐫 
    |λ(𝐤,𝐩,-𝐪,-𝐫)|^2 
    g_G(𝐤) g_G(𝐩) g_G(𝐪) g_G(𝐫)/2 ω_𝐤 2ω_𝐩 2 ω_𝐪 2ω_𝐫 
    (2 π)^d+1 δ (ω_𝐤 + ω_𝐩 + ω_𝐪 - ω_𝐫) δ(𝐤 + 𝐩 - 𝐪 - 𝐫) 
       ×[(f(𝐤,t) + 1)(f(𝐩,t) + 1) (f(𝐪,t) + 1) f(𝐫,t) - f(𝐤,t) f(𝐩,t) f(𝐪,t) (f(𝐫,t) + 1)  ],


where the minus sign in the first line comes from the one in front of δ(q^0 + ω_𝐪). Finally, we assume that f(t,𝐩) has no explicit momentum dependence, i.e., it only depends on 𝐩 via ω_𝐩  (cf. <cit.>), and change the integration variable 𝐩→ -𝐩:


    (<ref>) 
    =
    ∫_𝐩,𝐪,𝐫   
    |λ(𝐤,𝐩,𝐪,-𝐫)|^2 g_G(𝐤) g_G(𝐩) g_G(𝐪) g_G(𝐫)/2 ω_𝐤 2ω_𝐩 2 ω_𝐪 2ω_𝐫 (2 π)^d+1 δ (ω_𝐤 + ω_𝐩 + ω_𝐪 - ω_𝐫) δ(𝐤 + 𝐩 + 𝐪 - 𝐫) 
       ×[(f_𝐤 + 1)(f_𝐩 + 1) (f_𝐪 + 1) f_𝐫 - f_𝐤 f_𝐩 f_𝐪 (f_𝐫 + 1)  ],


where we have used that g_G(𝐩) and ω_𝐩 are even functions and adopted the short-hand notation f(𝐤,t) ≡ f_𝐤. 

The remaining terms can be computed following exactly the same procedure. With that we conclude:


    ∂_t f_𝐤 = I[f](𝐤,t) = I_3[f](𝐤,t) + I_4[f](𝐤,t),


with


    I_4[f](𝐤,t) = 
    I^2 ↔ 2[f](𝐤,t) + 
    I^3 ↔ 1[f](𝐤,t) + 
    I^1 ↔ 3[f](𝐤,t) +
    I^4 ↔ 0[f](𝐤,t),


where




    I^2 ↔ 2[f](𝐤,t) 
    =
    ∫_𝐩,𝐪,𝐫    
    |λ(𝐤,𝐩,-𝐪,-𝐫)|^2 g_G(𝐤) g_G(𝐩) g_G(𝐪) g_G(𝐫)/2 ω_𝐤 2ω_𝐩 2 ω_𝐪 2ω_𝐫 (2 π)^d+1 δ (ω_𝐤 + ω_𝐩 - ω_𝐪 - ω_𝐫) δ(𝐤 + 𝐩 - 𝐪 - 𝐫) 
       ×[(f_𝐤 + 1)(f_𝐩 + 1) f_𝐪 f_𝐫 - f_𝐤 f_𝐩 (f_𝐪 + 1) (f_𝐫 + 1)  ] +
    (𝐩→ -𝐩, -𝐪→𝐪) + (𝐩→ -𝐩, -𝐫→𝐫),


    I^3 ↔ 1[f](𝐤,t) 
    =
    ∫_𝐩,𝐪,𝐫   
    |λ(𝐤,𝐩,𝐪,-𝐫)|^2 g_G(𝐤) g_G(𝐩) g_G(𝐪) g_G(𝐫)/2 ω_𝐤 2ω_𝐩 2 ω_𝐪 2ω_𝐫  (2 π)^d+1 δ (ω_𝐤 + ω_𝐩 + ω_𝐪 - ω_𝐫) δ(𝐤 + 𝐩 + 𝐪 - 𝐫) 
       ×[(f_𝐤 + 1)(f_𝐩 + 1) (f_𝐪 + 1) f_𝐫 - f_𝐤 f_𝐩 f_𝐪 (f_𝐫 + 1)  ] +
    (𝐩→ -𝐩, -𝐫→𝐫) + (𝐪→ -𝐪, -𝐫→𝐫),


    I^1 ↔ 3[f](𝐤,t) 
    =
    ∫_𝐩,𝐪,𝐫   
    |λ(𝐤,-𝐩,-𝐪,-𝐫)|^2 g_G(𝐤) g_G(𝐩) g_G(𝐪) g_G(𝐫)/2 ω_𝐤 2ω_𝐩 2 ω_𝐪 2ω_𝐫 (2 π)^d+1 δ (ω_𝐤 - ω_𝐩 - ω_𝐪 - ω_𝐫) δ(𝐤 - 𝐩 - 𝐪 - 𝐫) 
       ×[(f_𝐤 + 1)f_𝐩 f_𝐪 f_𝐫 - f_𝐤 (f_𝐩 + 1) (f_𝐪 + 1) (f_𝐫 + 1)  ],


    I^4 ↔ 0[f](𝐤,t)
    =
    ∫_𝐩,𝐪,𝐫    
    |λ(𝐤,𝐩,𝐪,𝐫)|^2 g_G(𝐤) g_G(𝐩) g_G(𝐪) g_G(𝐫)/2 ω_𝐤 2ω_𝐩 2 ω_𝐪 2ω_𝐫 (2 π)^d+1 δ (ω_𝐤 + ω_𝐩 + ω_𝐪 + ω_𝐫) δ(𝐤 + 𝐩 + 𝐪 + 𝐫) 
       ×[(f_𝐤 + 1)(f_𝐩 + 1) (f_𝐪 + 1) (f_𝐫 + 1) - f_𝐤 f_𝐩 f_𝐪 f_𝐫].



Above, 𝐩↔ -𝐩 means that the same term has to be taken but with the change of sign for 𝐩 and f_𝐩↔ f_𝐩 + 1 and λ(𝐤,𝐩,𝐪,𝐫) is the momentum-dependent 4-coupling constant of the action (<ref>):  


    λ(𝐤,𝐩,𝐪,𝐫) = (𝐤·𝐩) (𝐪·𝐫)/2 m^2 g_G(𝐤 - 𝐩) + perm^s.


Similarly, one can derive the 3-vertex collision integral, which in 2-loop order reads:


    I_3[f](𝐤,t) = 
    I^1 ↔ 2[f](𝐤,t) + 
    I^2 ↔ 1[f](𝐤,t) +  
    I^3 ↔ 0[f](𝐤,t),


with 




    I^1 ↔ 2[f](𝐤,t) 
       =
    ∫_𝐩,𝐪 
    |γ(𝐤,-𝐩,-𝐪)|^2 g_G(𝐤) g_G(𝐩) g_G(𝐪)/2 ω_𝐤 2ω_𝐩 2 ω_𝐪 (2 π)^d+1 δ (ω_𝐤 - ω_𝐩 - ω_𝐪) δ(𝐤 - 𝐩 - 𝐪) 
       ×[(f_𝐤 + 1)f_𝐩 f_𝐪 - f_𝐤 (f_𝐩 + 1) (f_𝐪 + 1)  ],
    
    I^2 ↔ 1[f](𝐤,t) 
       =
    ∫_𝐩,𝐪 
    |γ(𝐤,𝐩,-𝐪)|^2 g_G(𝐤) g_G(𝐩) g_G(𝐪)/2 ω_𝐤 2ω_𝐩 2 ω_𝐪 (2 π)^d+1 δ (ω_𝐤 + ω_𝐩 - ω_𝐪) δ(𝐤 + 𝐩 - 𝐪) 
       ×[(f_𝐤 + 1)(f_𝐩 + 1) f_𝐪 - f_𝐤 f_𝐩 (f_𝐪 + 1)  ] + (𝐩→ -𝐩, -𝐪→𝐪),
    
    I^3 ↔ 0[f](𝐤,t)
       =
    ∫_𝐩,𝐪 
    |γ(𝐤,𝐩,𝐪)|^2 g_G(𝐤) g_G(𝐩) g_G(𝐪)/2 ω_𝐤 2ω_𝐩 2 ω_𝐪 (2 π)^d+1 δ (ω_𝐤 + ω_𝐩 + ω_𝐪) δ(𝐤 + 𝐩 + 𝐪)  
       ×[(f_𝐤 + 1)(f_𝐩 + 1) (f_𝐪 + 1) - f_𝐤 f_𝐩 f_𝐪].



Here, |γ(𝐤,𝐩,𝐪)|^2 should be understood on-shell, i.e.,


    γ(𝐤,𝐩,𝐪) = (𝐤·𝐩) ω_𝐪/m g_G(𝐪) + perm^s.






§ CANONICAL SCALING OF THE ACTION



Let us determine the scaling dimension of each term of the effective action (<ref>). For instance, consider the quadratic part:


    S^(2)(s^-1/β t ) 
       =
    ∫_t_ref,𝒞^s^-1/βt d t' ∫d^d k/(2π)^d θ_a(𝐤,t') ı D_a^-1(𝐤,t') θ_a(-𝐤,t')  
       =
    s^-1/β + d∫_s^1/β t_ref,𝒞^t d t̃∫d^d k'/(2π)^d θ_a(s𝐤',s^-1/βt̃) i D_a^-1(s𝐤',s^-1/βt̃) θ_a(-s𝐤',s^-1/βt̃)
       ≈
    s^-1/β + d - α/β - γ + 2z∫_t_ref,𝒞^t d t̃∫d^d k/(2π)^d θ_a(𝐤,t̃) i D_a^-1(𝐤,t̃) θ_a(-𝐤,t̃) = s^d_S,2 S^(2)(t),


which yields exactly (<ref>). Here, in the last line we assumed that t_ref is sent to the remote past, such that the Schwinger-Keldysh integral is almost invariant under t_ref→ s^1/βt_ref, and used α = d β.   

Repeating the same steps for interaction terms gives


    S^(3)(s^-1/β t) ≈
    s^- 1/β + 2d - 3 α/2β + z + 2 - γ S^(3)(t) = s^d_S,3 S^(3)(t)


and

    S^(4)(s^-1/β t) ≈
    s^-1/β + 3d - 2α/β + 4 - γ S^(4)(t) = s^d_S,4 S^(4)(t),


which results in (<ref>) and (<ref>), respectively. It should be remarked that terms involving derivatives with respect to time lead to ambiguities. Indeed, on one hand, they should be rescaled with t → s^-1/β t following the definition (<ref>). On the other hand, the on-shell scaling of the frequency implies t → s^-z t. Nevertheless, since in our case z = 1/β, this gives the same result. 





§ HYDRODYNAMIC DECOMPOSITION



In this Appendix, we provide the definition of the hydrodynamic decomposition used to obtain the data shown in HydroDecompEvo.

The three-component system is described, in a hydrodynamic formulation <cit.>, by the total density ρ, the spin vector f_μ and the nematic tensor n_μν:


    ρ   =  ∑_a Φ_a^†Φ_a 
    
      f_μ   = 1/ρ∑_a, a^'Φ_a^†( f_μ)_a a^'Φ_a^' ,    
      
    
      n_μν   = 1/ρ∑_a, a^'Φ_a^†( n_μν)_a  a^'Φ_a^',


μ=x,y,z, with f_μ being the spin-1 matrices in the fundamental representation, and the nematic or quadrupole tensor representation n_μν =  ( f_μf_ν + f_νf_μ )/2.
The above densities are related, via continuity equations <cit.>, to the mass current 𝐣 defining the superfluid velocity field 𝐯, as well as to the spin and nematic currents, 𝐣^(s)_μ and 𝐣^(n)_μν, respectively: 


    𝐣   =ρ𝐯
      = 1/2 m ı∑_a[ Φ_a^†( ∇Φ_a ) 
                                                            - (∇Φ_a^†) Φ_a ],
      
    𝐣^(s)_μ   = 1/2 m ı∑_a,a'( f_μ)_a a^'[ Φ_a^†( ∇Φ_a') 
                                                - (∇Φ_a^†) Φ_a'],
      
    𝐣^(n)_μν   = 1/2 m ı∑_a,a'( n_μν)_a a^'[ Φ_a^†( ∇Φ_a') 
                                                   - (∇Φ_a^†) Φ_a'].


Note that the nematic tensor vanishes for a fully polarized (ferromagnetic) gas but becomes important in the situations considered in this work.
The spin and nematic currents can be expressed in terms of the hydrodynamic variables ρ, f_μ, n_μν, and 𝐯, as well as the gradients ∇ f_μ and ∇ n_μν. 
Furthermore, diagonalization of f_μ and n_μν shows that the six independent variables are ρ, the total phase θ which defines the fluid momentum m𝐯=∇θ, as well as three Euler angles and a polarization angle <cit.>.

For expressing the hydrodynamic energy it is useful to define the generalized velocities, corresponding to the quantum-pressure (q), the spin (s), the nematic (n), the incompressible (i), and compressible (c) parts:


    w⃗^(q)   = m^-1∇√(ρ),
       w⃗^(i,c) =    √(ρ) v⃗^(i,c),
    
    w⃗_μ^(s)   = (2m)^-1√(ρ) ∇ f_μ,
       w⃗_μν^(n) =    (2m)^-1√(2ρ) ∇ n_μν.


Here, v⃗^(i,c) are obtained by a Helmholtz decomposition of v⃗ = v⃗^(i) + v⃗^(c), with the incompressible part having a vanishing divergence,  ∇·v⃗^(i) = 0, the compressible part a vanishing curl, ∇×v⃗^(c) = 0.


Using the hydrodynamic variables we can express the energy as 


    E = E_kin
    + g/2∫_r⃗ρ^2 ,


where the kinetic part  reads as


    E_kin 
    = ∫_r⃗{1/2 m[ 
    ( ∇√(ρ))^2
    + ρ/4( ∇ f_μ)^2 
    + ρ/2( ∇ n_μν)^2 
    ] 
    + m/2ρv⃗^2 }.


Hence, in Fourier space, the kinetic-energy spectrum is given by the correlation functions of the generalized velocities


    ε_kin(k) = ε^(q) (k) + ε^(c) (k) + ε^(i) (k)+ ε^(s) (k) + ε^(n) (k),


averaged over the orientation of the momentum vector: 


    ε^(δ) (k)     = m/2∫ d Ω_𝐤⟨|w⃗^(δ) (𝐤) |^2 ⟩,   (δ=q,i,c)
    
    ε^(s) (k)     = m/2∫ d Ω_𝐤⟨w⃗_μ ^(s) (𝐤) ·w⃗_μ ^(s) (𝐤)  ⟩,
    
    ε^(n) (k)     = m/2∫ d Ω_𝐤⟨w⃗_μν ^(n) (𝐤) ·w⃗_μν ^(n) (𝐤)  ⟩,


where Einstein's sum convention is implied.
The respective total energies are obtained as E^(δ)=∫ d^dx∫ dk  k^d-1ε^(δ)(k).
The spectrum of the kinetic energy can then be used to calculate corresponding occupation numbers using the relation


    n^(δ) (k)  = 2m k^-2ε^ (δ) (k),


where δ = q, i, c, s, n.
The total occupation number is then approximately given by 

    n_tot (k)  ≈ 2m k^-2ε_kin  (k)


(see, e.g., <cit.> for a discussion of deviations).

We emphasize that the spin f_μ and quadrupole operators n_μν form SU(2) subgroups of the special unitary group SU(3), i.e., not only {f_x,f_y,f_z} form an su(2) algebra of generators of SU(2) but also, e.g., {f_x,n_yz,n_zz-n_yy} and {n_xz,n_yz,f_z} <cit.>.
On the basis of the field representations spinop and nemop, one can furthermore show that for the case of identical near-constant bosonic densities ρ_a^(0) and comparatively large phase fluctuations, the fluctuations epss and epsn of the spin and nematic operators are determined by combinations of the relative phases between the Bose fields.











apsrev4-1










98
fxundefined [1]
 ifx#1

fnum [1]
 #1firstoftwo
 secondoftwo
 
fx [1]
 #1firstoftwo
 secondoftwo
 





noop [0]secondoftwo

ref[1]@startlink#1@href
href[1]#1@endlink
anitize@url [0]`
12`$12`&12`#12`1̂2`_12`%12
startlink[1]
endlink[0]

rl [1]href #1












@bib@innerbibempty

[Gring et al.(2012)Gring,
  Kuhnert, Langen, Kitagawa,
  Rauer, Schreitl, Mazets,
  Smith, Demler, and Schmiedmayer]Gring2011a
  author author M. Gring, author M. Kuhnert,
  author T. Langen, author T. Kitagawa, author
  B. Rauer, author M. Schreitl, author I. Mazets, author D. A.  Smith, author E. Demler,  and author J. Schmiedmayer, 10.1126/science.1224953 journal journal Science volume
  337, pages 1318 (year 2012), http://arxiv.org/abs/1112.0013 arXiv:1112.0013 [cond-mat.quant-gas]
  NoStop
[Berges et al.(2004)Berges,
  Borsanyi, and Wetterich]Berges:2004ce
  author author J. Berges, author S. Borsanyi,   and author C. Wetterich, 10.1103/PhysRevLett.93.142002 journal journal Phys. Rev. Lett. volume 93, pages 142002 (year 2004)NoStop
[Langen et al.(2015)Langen,
  Erne, Geiger, Rauer,
  Schweigler, Kuhnert, Rohringer, Mazets, Gasenzer, and  Schmiedmayer]Langen2015b.Science348.207
  author author T. Langen, author S. Erne,
  author R. Geiger, author B. Rauer, author
  T. Schweigler, author
  M. Kuhnert, author W. Rohringer, author I. E.  Mazets, author T. Gasenzer,  and author J. Schmiedmayer, 10.1126/science.1257026 journal journal Science volume
  348, pages 207 (year 2015), http://arxiv.org/abs/1411.7185 arXiv:1411.7185 [cond-mat.quant-gas]
  NoStop
[Jaynes(1957a)]Jaynes1957a
  author author E. T. Jaynes, 10.1103/PhysRev.106.620 journal
  journal Phys. Rev. volume 106,  pages 620 (year 1957a)NoStop
[Jaynes(1957b)]Jaynes1957b
  author author E. T. Jaynes, 10.1103/PhysRev.108.171 journal
  journal Phys. Rev. volume 108,  pages 171 (year 1957b)NoStop
[Schreiber et al.(2015)Schreiber, Hodgman, Bordia, Lüschen, Fischer, Vosk, Altman, Schneider, and Bloch]Schreiber2015a.Science349.842
  author author M. Schreiber, author S. S. Hodgman, author P. Bordia,
  author H. P. Lüschen,
  author M. H. Fischer, author R. Vosk, author
  E. Altman, author U. Schneider,  and author I. Bloch, 10.1126/science.aaa7432
  journal journal Science volume 349, pages 842 (year 2015), http://arxiv.org/abs/1501.05661 arXiv:1501.05661 [cond-mat.quant-gas]
  NoStop
[Braun et al.(2015)Braun,
  Friesdorf, Hodgman, Schreiber, Ronzheimer, Riera, del Rey, Bloch, Eisert, and Schneider]Braun2014a.arXiv1403.7199B
  author author S. Braun, author M. Friesdorf,
  author S. S. Hodgman, author M. Schreiber, author
  J. P. Ronzheimer, author
  A. Riera, author M. del Rey, author I. Bloch, author J. Eisert,  and author U. Schneider, 10.1073/pnas.1408861112
  journal journal PNAS volume 112, pages 3641 (year 2015), http://arxiv.org/abs/1403.7199 arXiv:1403.7199 [cond-mat.quant-gas]
  NoStop
[Nicklas et al.(2015)Nicklas, Karl, Höfer, Johnson, Muessel, Strobel, Tomkovic, Gasenzer, and Oberthaler]Nicklas:2015gwa
  author author E. Nicklas, author M. Karl,
  author M. Höfer, author A. Johnson, author
  W. Muessel, author H. Strobel, author J. Tomkovic, author T. Gasenzer,  and author M. K. Oberthaler, 10.1103/PhysRevLett.115.245301 journal journal
  Phys. Rev. Lett. volume 115, pages
  245301 (year 2015), http://arxiv.org/abs/1509.02173 arXiv:1509.02173 [cond-mat.quant-gas]
  NoStop

[Navon et al.(2015)Navon, Gaunt, Smith, and  Hadzibabic]Navon2015a.Science.347.167N
  author author N. Navon, author A. L. Gaunt, author R. P. Smith,  and author Z. Hadzibabic, 10.1126/science.1258676 journal journal Science volume
  347, pages 167 (year 2015), http://arxiv.org/abs/1410.8487 arXiv:1410.8487 [cond-mat.quant-gas]
  NoStop
[Eigen et al.(2018)Eigen, Glidden, Lopes, Cornell, Smith, and Hadzibabic]Eigen2018a.arXiv180509802E
  author author C. Eigen, author J. A. P. Glidden, author R. Lopes,
  author E. A. Cornell, author R. P. Smith,  and author Z. Hadzibabic, 10.1038/s41586-018-0674-1 journal journal
  Nature volume 563, pages 221
  (year 2018), http://arxiv.org/abs/1805.09802
  arXiv:1805.09802 [cond-mat.quant-gas] NoStop
[Rauer et al.(2018)Rauer,
  Erne, Schweigler, Cataldini,
  Tajik, and Schmiedmayer]Rauer2017a.arXiv170508231R.Science360.307
  author author B. Rauer, author S. Erne,
  author T. Schweigler, author F. Cataldini, author
  M. Tajik,  and author
  J. Schmiedmayer, 10.1126/science.aan7938 journal journal
  Science volume 360, pages 307
  (year 2018), http://arxiv.org/abs/1705.08231
  arXiv:1705.08231 [cond-mat.quant-gas] NoStop
[Navon et al.(2016)Navon, Gaunt, Smith, and  Hadzibabic]Navon2016a.Nature.539.72
  author author N. Navon, author A. L. Gaunt, author R. P. Smith,  and author Z. Hadzibabic, 10.1038/nature20114 journal
  journal  volume 539, pages 72 (year 2016), http://arxiv.org/abs/1609.01271 arXiv:1609.01271 [cond-mat.quant-gas]
  NoStop
[Prüfer et al.(2018)Prüfer, Kunkel, Strobel, Lannig, Linnemann, Schmied, Berges, Gasenzer, and Oberthaler]Prufer:2018hto
  author author M. Prüfer, author P. Kunkel,
  author H. Strobel, author S. Lannig, author
  D. Linnemann, author
  C.-M. Schmied, author
  J. Berges, author T. Gasenzer,  and author M. K. Oberthaler, 10.1038/s41586-018-0659-0 journal journal
  Nature volume 563, pages 217
  (year 2018), http://arxiv.org/abs/1805.11881
  arXiv:1805.11881 [cond-mat.quant-gas] NoStop

[Erne et al.(2018)Erne,
  Bücker, Gasenzer, Berges, and Schmiedmayer]Erne:2018gmz
  author author S. Erne, author R. Bücker,
  author T. Gasenzer, author J. Berges,  and author
  J. Schmiedmayer, 10.1038/s41586-018-0667-0 journal journal
  Nature volume 563, pages 225
  (year 2018), http://arxiv.org/abs/1805.12310
  arXiv:1805.12310 [cond-mat.quant-gas] NoStop

[Hohenberg and Halperin(1977)]Hohenberg1977a
  author author P. C. Hohenberg and author B. I. Halperin, 10.1103/RevModPhys.49.435 journal journal Rev. Mod. Phys. volume
  49, pages 435 (year 1977)NoStop
[Janssen(1979)]Janssen1979a
  author author H. Janssen, in 10.1007/3-540-09523-3 booktitle Dynamical critical phenomena and related topics, Lecture Notes
  in Physics, vol. 104 (publisher Springer, Heidelberg 1979,  year 1979) p. pages 26NoStop
[Bray(1994)]Bray1994a.AdvPhys.43.357
  author author A. J. Bray, 10.1080/00018739400101505 journal
  journal Adv. Phys. volume 43,  pages 357 (year 1994), http://arxiv.org/abs/cond-mat/9501089 arXiv:cond-mat/9501089 NoStop
[Calabrese and Gambassi(2005)]Calabrese2005a.JPA38.05.R133
  author author P. Calabrese and author A. Gambassi, http://stacks.iop.org/0305-4470/38/i=18/a=R01
  journal journal J. Phys. A: Math. Gen. volume 38, pages R133 (year
  2005), http://arxiv.org/abs/cond-mat/0410357
  arXiv:cond-mat/0410357 NoStop
[Frisch(1995)]Frisch1995a
  author author U. Frisch, 10.1017/CBO9781139170666 title Turbulence: The Legacy of A. N. Kolmogorov (publisher CUP, Cambridge, UK, year 1995)NoStop
[Zakharov et al.(1992)Zakharov, L'vov, and Falkovich]Zakharov1992a
  author author V. E. Zakharov, author V. S. L'vov,  and author G. Falkovich, 10.1007/978-3-642-50052-7 title Kolmogorov Spectra of Turbulence I: Wave Turbulence (publisher Springer, Berlin, year 1992)NoStop
[Nazarenko(2011)]Nazarenko2011a
  author author S. Nazarenko, 10.1007/978-3-642-15942-8 title Wave turbulence, series Lecture Notes in Physics  No. number 825 (publisher Springer, address Heidelberg, year 2011) pp. pages XVI,
  279 S.Stop
[Svistunov(1991)]Svistunov1991a
  author author B. Svistunov, @noop  journal journal J.
  Mosc. Phys. Soc. volume 1, pages 373
  (year 1991)NoStop
[Kagan et al.(1992)Kagan,
  Svistunov, and Shlyapnikov]Kagan1992a
  author author Y. Kagan, author B. V. Svistunov,  and author G. V. Shlyapnikov, http://www.jetp.ac.ru/cgi-bin/e/index/e/74/2/p279?a=list journal journal [Zh. Eksp. Teor. Fiz. 101, 528 (1992)] Sov.
  Phys. JETP volume 74, pages 279
  (year 1992)NoStop
[Kagan and Svistunov(1994)]Kagan1994a
  author author Y. Kagan and author B. V. Svistunov, http://www.jetp.ac.ru/cgi-bin/e/index/e/78/2/p187?a=list journal journal [Zh. Eksp. Teor. Fiz. 105, 353 (1994)] Sov.
  Phys. JETP volume 78, pages 187
  (year 1994)NoStop
[Kagan(1995)]Kagan1995a
  author author Y. Kagan, @noop  title Bose-Einstein
  Condensation (publisher CUP, Cambridge, UK,  year 1995) p. pages 202NoStop
[Semikoz and Tkachev(1995)]Semikoz1995a.PhysRevLett.74.3093
  author author D. V. Semikoz and author I. I. Tkachev, 10.1103/PhysRevLett.74.3093 journal journal Phys. Rev. Lett. volume 74, pages 3093 (year
  1995)NoStop
[Semikoz and Tkachev(1997)]Semikoz1997a
  author author D. Semikoz and author I. Tkachev, 10.1103/PhysRevD.55.489 journal
  journal Phys. Rev. D volume 55,  pages 489 (year 1997), http://arxiv.org/abs/hep-ph/9507306 hep-ph/9507306 NoStop
[Berloff and Svistunov(2002)]Berloff2002a
  author author N. G. Berloff and author B. V. Svistunov, 10.1103/PhysRevA.66.013603 journal journal Phys. Rev. A volume
  66, pages 013603 (year 2002), http://arxiv.org/abs/cond-mat/0107209 cond-mat/0107209 NoStop
[Kozik and Svistunov(2004)]Kozik2004a.PhysRevLett.92.035301
  author author E. Kozik and author B. Svistunov, 10.1103/PhysRevLett.92.035301 journal journal Phys. Rev. Lett. volume 92, pages 035301 (year
  2004), http://arxiv.org/abs/cond-mat/0308193 cond-mat/0308193 NoStop
[Kozik and Svistunov(2005a)]Kozik2005a.PhysRevLett.94.025301
  author author E. Kozik and author B. Svistunov, 10.1103/PhysRevLett.94.025301 journal journal Phys. Rev. Lett. volume 94, pages 025301 (year
  2005a), http://arxiv.org/abs/cond-mat/0408241 cond-mat/0408241 NoStop
[Kozik and Svistunov(2005b)]Kozik2005a.PhysRevB.72.172505
  author author E. Kozik and author B. Svistunov, 10.1103/PhysRevB.72.172505 journal journal Phys. Rev. B volume
  72, pages 172505 (year
  2005b), http://arxiv.org/abs/cond-mat/0505020 cond-mat/0505020 NoStop
[Kozik and Svistunov(2009)]Kozik2009a
  author author E. V. Kozik and author B. V. Svistunov, 10.1007/s10909-009-9914-y journal journal J. Low Temp. Phys. volume 156, pages 215 (year 2009), http://arxiv.org/abs/0904.1379 arXiv:0904.1379 [cond-mat.stat-mech]
  NoStop
[Tsubota(2008)]Tsubota2008a
  author author M. Tsubota, 10.1143/JPSJ.77.111006 journal
  journal J. Phys. Soc. Jpn. volume
  77, pages 111006 (year 2008), http://arxiv.org/abs/0806.2737 arXiv:0806.2737 [cond-mat.other]
  NoStop
[Vinen(2006)]Vinen2006a
  author author W. Vinen, 10.1007/s10909-006-9240-6 journal
  journal J. Low Temp. Phys. volume
  145, pages 7 (year 2006)NoStop
[Aarts et al.(2000)Aarts,
  Bonini, and Wetterich]Aarts2000a.PhysRevD.63.025012
  author author G. Aarts, author G. F. Bonini,   and author C. Wetterich, 10.1103/PhysRevD.63.025012 journal journal Phys. Rev. D volume 63, pages 025012 (year 2000), http://arxiv.org/abs/hep-ph/0007357 arXiv:hep-ph/0007357 NoStop
[Lamacraft(2007)]Lamacraft2007.PhysRevLett.98.160404
  author author A. Lamacraft, 10.1103/PhysRevLett.98.160404 journal journal Phys. Rev. Lett. volume 98, pages 160404 (year 2007),  http://arxiv.org/abs/cond-mat/0611017 arXiv:cond-mat/0611017
  [cond-mat.stat-mech] NoStop
[Rossini et al.(2009)Rossini, Silva, Mussardo, and  Santoro]Rossini2009a.PhysRevLett.102.127204
  author author D. Rossini, author A. Silva,
  author G. Mussardo,  and  author G. E. Santoro, 10.1103/PhysRevLett.102.127204 journal journal Phys. Rev. Lett. volume 102, pages 127204 (year 2009), http://arxiv.org/abs/0810.5508 arXiv:0810.5508 [cond-mat.stat-mech]
  NoStop
[Dalla Torre et al.(2013)Dalla Torre, Demler, and Polkovnikov]DallaTorre2013.PhysRevLett.110.090404
  author author E. G. Dalla Torre, author E. Demler,
   and author A. Polkovnikov,  10.1103/PhysRevLett.110.090404 journal
  journal Phys. Rev. Lett. volume 110,  pages 090404 (year 2013), http://arxiv.org/abs/1211.5145 arXiv:1211.5145 [cond-mat.quant-gas]
  NoStop
[Gambassi and Calabrese(2011)]Gambassi2011a.EPL95.6
  author author A. Gambassi and author P. Calabrese, 10.1209/0295-5075/95/66007 journal journal Europhys. Lett. volume
  95, pages 66007 (year 2011), http://arxiv.org/abs/1012.5294 arXiv:1012.5294 [cond-mat.stat-mech]
  NoStop
[Sciolla and Biroli(2013)]Sciolla2013a.PhysRevB.88.201110
  author author B. Sciolla and author G. Biroli, 10.1103/PhysRevB.88.201110 journal journal Phys. Rev. B volume
  88, pages 201110 (year 2013), http://arxiv.org/abs/1211.2572 arXiv:1211.2572 [cond-mat.stat-mech]
  NoStop
[Smacchia et al.(2015)Smacchia, Knap, Demler, and Silva]Smacchia2015a.PhysRevB.91.205136
  author author P. Smacchia, author M. Knap,
  author E. Demler,  and author A. Silva, 10.1103/PhysRevB.91.205136 journal journal Phys.
  Rev. B volume 91, pages 205136
  (year 2015), http://arxiv.org/abs/1409.1883
  arXiv:1409.1883 [cond-mat.stat-mech] NoStop
[Maraga et al.(2015)Maraga,
  Chiocchetta, Mitra, and Gambassi]Maraga2015a.PhysRevE.92.042151
  author author A. Maraga, author A. Chiocchetta,
  author A. Mitra,  and author A. Gambassi, 10.1103/PhysRevE.92.042151 journal journal Phys.
  Rev. E volume 92, pages 042151
  (year 2015), http://arxiv.org/abs/1506.04528
  arXiv:1506.04528 [cond-mat.stat-mech] NoStop
[Maraga et al.(2016)Maraga,
  Smacchia, and Silva]Maraga2016b.PhysRevB.94.245122
  author author A. Maraga, author P. Smacchia,   and author A. Silva, 10.1103/PhysRevB.94.245122 journal journal Phys. Rev. B volume 94, pages 245122 (year 2016), http://arxiv.org/abs/1602.01763 arXiv:1602.01763 [cond-mat.stat-mech]
  NoStop
[Chiocchetta et al.(2015)Chiocchetta, Tavora, Gambassi, and  Mitra]Chiocchetta2015a.PhysRevB.91.220302
  author author A. Chiocchetta, author M. Tavora,
  author A. Gambassi,  and  author A. Mitra, 10.1103/PhysRevB.91.220302 journal journal Phys. Rev. B volume 91, pages 220302 (year 2015), http://arxiv.org/abs/1411.7939 arXiv:1411.7939 [cond-mat.quant-gas]
  NoStop
[Chiocchetta et al.(2016a)Chiocchetta, Tavora, Gambassi, and Mitra]Chiocchetta2016a.PhysRevB.94.134311
  author author A. Chiocchetta, author M. Tavora,
  author A. Gambassi,  and  author A. Mitra, 10.1103/PhysRevB.94.134311 journal journal Phys. Rev. B volume 94, pages 134311 (year 2016a), http://arxiv.org/abs/1604.04614 arXiv:1604.04614 [cond-mat.stat-mech]
  NoStop
[Chiocchetta et al.(2016b)Chiocchetta, Gambassi, Diehl, and Marino]Chiocchetta:2016waa.PhysRevB.94.174301
  author author A. Chiocchetta, author A. Gambassi, author S. Diehl,   and author J. Marino, 10.1103/PhysRevB.94.174301 journal journal Phys. Rev. B volume 94, pages 174301 (year 2016b), http://arxiv.org/abs/1606.06272 arXiv:1606.06272 [cond-mat.stat-mech]
  NoStop
[Chiocchetta et al.(2017)Chiocchetta, Gambassi, Diehl, and  Marino]Chiocchetta2016b.161202419C.PhysRevLett.118.135701
  author author A. Chiocchetta, author A. Gambassi, author S. Diehl,   and author J. Marino, 10.1103/PhysRevLett.118.135701 journal journal Phys. Rev. Lett. volume 118, pages 135701 (year 2017), http://arxiv.org/abs/1612.02419 arXiv:1612.02419 [cond-mat.stat-mech]
  NoStop
[Marino and Diehl(2016)]Marino2016a.PhysRevLett.116.070407
  author author J. Marino and author S. Diehl,  10.1103/PhysRevLett.116.070407 journal
  journal Phys. Rev. Lett. volume 116,  pages 070407 (year 2016), http://arxiv.org/abs/1508.02723 arXiv:1508.02723 [cond-mat.quant-gas]
  NoStop
[Marino and Diehl(2016)]Marino2016PhRvB..94h5150M
  author author J. Marino and author S. Diehl, 10.1103/PhysRevB.94.085150 journal journal Phys. Rev. B volume
  94, pages 085150 (year 2016), http://arxiv.org/abs/1606.00452 arXiv:1606.00452 [cond-mat.quant-gas]
  NoStop
[Damle et al.(1996)Damle,
  Majumdar, and Sachdev]Damle1996a.PhysRevA.54.5037
  author author K. Damle, author S. N. Majumdar,
   and author S. Sachdev, 10.1103/PhysRevA.54.5037 journal journal Phys. Rev. A volume 54, pages 5037 (year 1996), http://arxiv.org/abs/cond-mat/9511058 arXiv:cond-mat/9511058 NoStop
[Mukerjee et al.(2007)Mukerjee, Xu, and Moore]Mukerjee2007a.PhysRevB.76.104519
  author author S. Mukerjee, author C. Xu,  and  author J. E. Moore, 10.1103/PhysRevB.76.104519 journal journal Phys. Rev. B volume 76, pages 104519 (year 2007), http://arxiv.org/abs/0704.3440 arXiv:0704.3440 [cond-mat.stat-mech]
  NoStop
[Williamson and Blakie(2016a)]Williamson2016a.PhysRevLett.116.025301
  author author L. A. Williamson and author P. B. Blakie, 10.1103/PhysRevLett.116.025301 journal journal Phys. Rev. Lett. volume 116, pages 025301 (year
  2016a), http://arxiv.org/abs/1504.06404
  arXiv:1504.06404 [cond-mat.quant-gas] NoStop
[Hofmann et al.(2014)Hofmann, Natu, and Das
  Sarma]Hofmann2014PhRvL.113i5702H
  author author J. Hofmann, author S. S. Natu,  and author S. Das
  Sarma, 10.1103/PhysRevLett.113.095702 journal journal Phys. Rev. Lett. volume 113, pages 095702 (year 2014),  http://arxiv.org/abs/1403.1284 arXiv:1403.1284
  [cond-mat.quant-gas] NoStop
[Williamson and Blakie(2016b)]Williamson2016a.PhysRevA.94.023608
  author author L. A. Williamson and author P. B. Blakie, 10.1103/PhysRevA.94.023608 journal journal Phys. Rev. A volume
  94, pages 023608 (year 2016b),  http://arxiv.org/abs/1605.04016 arXiv:1605.04016
  [cond-mat.quant-gas] NoStop
[Bourges and Blakie(2017)]Bourges2016a.arXiv161108922B.PhysRevA.95.023616
  author author A. Bourges and author P. B. Blakie, 10.1103/PhysRevA.95.023616 journal journal Phys. Rev. A volume
  95, pages 023616 (year 2017), http://arxiv.org/abs/1611.08922 arXiv:1611.08922 [cond-mat.quant-gas]
  NoStop
[Berges et al.(2008)Berges,
  Rothkopf, and Schmidt]Berges:2008wm
  author author J. Berges, author A. Rothkopf,   and author J. Schmidt, 10.1103/PhysRevLett.101.041603 journal journal Phys. Rev. Lett. volume 101, pages 041603 (year 2008), http://arxiv.org/abs/0803.0131 arXiv:0803.0131 [hep-ph] NoStop

[Berges and Hoffmeister(2009)]Berges:2008sr
  author author J. Berges and author G. Hoffmeister, 10.1016/j.nuclphysb.2008.12.017 journal journal Nucl. Phys. volume
  B813, pages 383 (year 2009), http://arxiv.org/abs/0809.5208 arXiv:0809.5208 [hep-th] NoStop

[Scheppach et al.(2010)Scheppach, Berges, and Gasenzer]Scheppach:2009wu
  author author C. Scheppach, author J. Berges,   and author T. Gasenzer, 10.1103/PhysRevA.81.033611 journal journal Phys. Rev. A volume 81, pages 033611 (year 2010), http://arxiv.org/abs/0912.4183 arXiv:0912.4183 [cond-mat.quant-gas]
  NoStop
[Berges and Sexty(2011)]Berges:2010ez
  author author J. Berges and author D. Sexty,  10.1103/PhysRevD.83.085004 journal journal Phys. Rev. D volume 83, pages 085004 (year 2011), http://arxiv.org/abs/1012.5944 arXiv:1012.5944 [hep-ph] NoStop
[Piñeiro Orioli et al.(2015)Piñeiro Orioli, Boguslavski, and Berges]Orioli:2015dxa
  author author A. Piñeiro Orioli, author K. Boguslavski,  and author J. Berges, 10.1103/PhysRevD.92.025041 journal journal Phys. Rev. D volume
  92, pages 025041 (year 2015), http://arxiv.org/abs/1503.02498 arXiv:1503.02498 [hep-ph] NoStop

[Berges(2016)]Berges:2015kfa
  author author J. Berges, in 10.1093/acprof:oso/9780198768166.001.0001
  booktitle Proc. Int. School on Strongly Interacting
  Quantum Systems Out of Equilibrium, Les Houches, editor
  edited by editor T. Giamarchi et
  al. (publisher OUP, Oxford, year 2016)  http://arxiv.org/abs/1503.02907 arXiv:1503.02907 [hep-ph]
  NoStop

[Walz et al.(2018)Walz,
  Boguslavski, and Berges]Walz:2017ffj.PhysRevD.97.116011
  author author R. Walz, author K. Boguslavski,   and author J. Berges, 10.1103/PhysRevD.97.116011 journal journal Phys. Rev. D volume 97, pages 116011 (year 2018), http://arxiv.org/abs/1710.11146 arXiv:1710.11146 [hep-ph] NoStop
[Chantesana et al.(2019)Chantesana, Piñeiro Orioli, and Gasenzer]Chantesana:2018qsb.PhysRevA.99.043620
  author author I. Chantesana, author A. Piñeiro Orioli,  and author T. Gasenzer, 10.1103/PhysRevA.99.043620
  journal journal Phys. Rev. A volume 99, pages 043620 (year 2019),  http://arxiv.org/abs/1801.09490 arXiv:1801.09490
  [cond-mat.quant-gas] NoStop

[Nowak et al.(2011)Nowak,
  Sexty, and Gasenzer]Nowak:2010tm
  author author B. Nowak, author D. Sexty,  and  author T. Gasenzer, 10.1103/PhysRevB.84.020506 journal journal Phys. Rev. B volume 84, pages 020506(R) (year 2011), http://arxiv.org/abs/1012.4437v2 arXiv:1012.4437v2 [cond-mat.quant-gas]
  NoStop

[Nowak et al.(2012)Nowak,
  Schole, Sexty, and Gasenzer]Nowak:2011sk
  author author B. Nowak, author J. Schole,
  author D. Sexty,  and author T. Gasenzer, 10.1103/PhysRevA.85.043627 journal journal Phys.
  Rev. A volume 85, pages 043627
  (year 2012), http://arxiv.org/abs/1111.6127
  arXiv:1111.6127 [cond-mat.quant-gas] NoStop

[Schole et al.(2012)Schole,
  Nowak, and Gasenzer]Schole:2012kt
  author author J. Schole, author B. Nowak,   and author T. Gasenzer, 10.1103/PhysRevA.86.013624 journal journal Phys. Rev. A volume 86, pages 013624 (year 2012), http://arxiv.org/abs/1204.2487 arXiv:1204.2487 [cond-mat.quant-gas]
  NoStop

[Nowak et al.(2014)Nowak, Schole, and Gasenzer]Nowak:2012gd
  author author B. Nowak, author J. Schole,   and author T. Gasenzer,  10.1088/1367-2630/16/9/093052 journal
  journal New J. Phys. volume 16,  pages 093052 (year 2014), http://arxiv.org/abs/1206.3181v2 arXiv:1206.3181v2 [cond-mat.quant-gas]
  NoStop
[Karl et al.(2013)Karl,
  Nowak, and Gasenzer]Karl:2013mn
  author author M. Karl, author B. Nowak,   and author T. Gasenzer,  10.1038/srep02394 journal journal Sci. Rep. volume 3, eid
  2394 (year 2013), 10.1038/srep02394, http://arxiv.org/abs/1302.1122 arXiv:1302.1122 [cond-mat.quant-gas]
  NoStop
[Karl et al.(2013)Karl,
  Nowak, and Gasenzer]Karl:2013kua
  author author M. Karl, author B. Nowak,  and  author T. Gasenzer, 10.1103/PhysRevA.88.063615 journal journal Phys. Rev. A volume 88, pages 063615 (year 2013), http://arxiv.org/abs/1307.7368 arXiv:1307.7368 [cond-mat.quant-gas]
  NoStop

[Karl and Gasenzer(2017)]Karl2017b.NJP19.093014
  author author M. Karl and author T. Gasenzer, http://stacks.iop.org/1367-2630/19/i=9/a=093014
  journal journal New J. Phys. volume 19, pages 093014 (year 2017),  http://arxiv.org/abs/1611.01163 arXiv:1611.01163
  [cond-mat.quant-gas] NoStop
[Schmied et al.(2019a)Schmied, Prüfer, Oberthaler, and Gasenzer]SchmiedPhysRevA.99.033611
  author author C.-M. Schmied, author M. Prüfer,
  author M. K. Oberthaler,  and  author T. Gasenzer, 10.1103/PhysRevA.99.033611 journal journal Phys. Rev. A volume 99, pages 033611 (year 2019a)NoStop
[Schmied et al.(2019a)Schmied, Gasenzer, and Blakie]Schmied:2019abm
  author author C. M. Schmied, author T. Gasenzer,   and author P. B. Blakie,  @noop    (year 2019a), http://arxiv.org/abs/1904.13222 arXiv:1904.13222 [cond-mat.quant-gas]
  NoStop

[Berges(2002)]Berges:2001fi
  author author J. Berges, 10.1016/S0375-9474(01)01295-7 journal journal Nucl. Phys. volume
  A699, pages 847 (year 2002), http://arxiv.org/abs/hep-ph/0105311 arXiv:hep-ph/0105311 NoStop
[Berges(2005)]Berges:2004yj
  author author J. Berges, 10.1063/1.1843591 journal
  journal AIP Conf. Proc. volume 739,  pages 3 (year 2005), http://arxiv.org/abs/hep-ph/0409233 arXiv:hep-ph/0409233 NoStop
[Schmied et al.(2019)Schmied, Mikheev, and Gasenzer]Schmied:2018upn.PhysRevLett.122.170404
  author author C.-M. Schmied, author A. N. Mikheev,  and author T. Gasenzer, 10.1103/PhysRevLett.122.170404 journal journal Phys. Rev. Lett. volume 122, pages 170404 (year
  2019), https://arxiv.org/abs/1807.07514 arXiv:1807.07514 
  [cond-mat.quant-gas]NoStop
[Madelung(1926)]Madelung1926a
  author author E. Madelung, 10.1007/BF01400372 journal
  journal Z. Phys. volume 40, pages 322 (year 1926)NoStop
[Note1()]Note1
  note Expanding around the time-independent homogeneous solution
  has to be distinguished from expanding around the condensate mean field
  which, during the evolution in the vicinity of an non-thermal fixed point, typically grows
  algebraically. The density corresponding to the field expectation value
  squared needs to be distinguished from the zero mode of the total density
  which is proportional to the total particle number and given by the sum of
  the condensate density and the integral over all non-zero momentum
  occupancies, ρ ^(0)_a = | "426830A Φ
  _a "526930B |^2 + ∑@ @ _𝐤 "426830A Φ ^†_a (𝐤) Φ
  _a(𝐤) "526930B. Hence, conservation of the
  total particle number in a closed system ensures the constancy of ρ
  ^(0)_a. Furthermore, one should note that, in the presence of a condensate
  in the ground state of, e.g., a 3D system, one would have ρ
  _a^(0)= | "426830A Φ _a(x)  "526930B |^2,
  while the approach here also applies to lower-dimensional systems where, in
  equilibrium, the ground state bears a quasi-condensate.Stop
[Note2()]Note2
  note Obviously, all of these modes are gapless Goldstone
  excitations. Nevertheless, in order to distinguish them, hereinafter we will
  only refer to the quadratic modes as Goldstone ones, whereas the linear one
  will be addressed as Bogoliubov.Stop
[Goldstone(1961)]Goldstone:1961eq
  author author J. Goldstone, 10.1007/BF02812722 journal
  journal Nuovo Cim. volume 19,  pages 154 (year 1961)NoStop

[Watanabe and Murayama(2012)]Watanabe2012a.PhysRevLett.108.251602
  author author H. Watanabe and author H. Murayama, 10.1103/PhysRevLett.108.251602 journal journal Phys. Rev. Lett. volume 108, pages 251602 (year 2012),  http://arxiv.org/abs/1203.0609 arXiv:1203.0609 [hep-ph]
  NoStop
[Watanabe and Murayama(2014)]Watanabe2014a.PhysRevX.4.031057
  author author H. Watanabe and author H. Murayama, 10.1103/PhysRevX.4.031057 journal journal Phys. Rev. X volume
  4, pages 031057 (year 2014), http://arxiv.org/abs/1402.7066 arXiv:1402.7066 [hep-ph] NoStop
[Pekker and Varma(2015)]Pekker2015aARCMP.6.269P
  author author D. Pekker and author C. M. Varma, 10.1146/annurev-conmatphys-031214-014350
  journal journal Ann. Rev. Cond. Mat. Phys.  volume 6, pages 269 (year
  2015), http://arxiv.org/abs/1406.2968 arXiv:1406.2968
  [cond-mat.supr-con] NoStop
[Feynman and Vernon(1963)]Feynman:1963fq
  author author R. P. Feynman and author F. L. Vernon, Jr., 10.1016/0003-4916(63)90068-X journal journal
  Annals Phys. volume 24, pages 118
  (year 1963)NoStop

[Note3()]Note3
  note The Hamiltonian of the (1D) (Tomonaga-)Luttinger-liquid
  model <cit.> reads H_TLL=(c_s/2π
  )@ dx +.1667em{K^-1[πΠ (x,t)]^2+K[∂ _xθ (x,t)]^2},
  with speed of sound c_s=√(v_Jv_N), Luttinger parameter
  K=√(v_J/v_N), phase stiffness v_J=πρ _0/m,
  and density stiffness v_J=g/π. The canonically conjugate momentum is
  defined in terms of the density fluctuation as Π (x,t)=δρ (x,t).
  This corresponds to the Lagrangian L_TLL=(2g)^-1@ dx +.1667em{[∂ _tθ (x,t)]^2-[c_s∂ _xθ
  (x,t)]^2}, and thus to the canonically conjugate field Π
  (x,t)=g^-1∂ _tθ (x,t).Stop
[Cazalilla et al.(2011)Cazalilla, Citro, Giamarchi, Orignac, and Rigol]Cazalilla2011a
  author author M. A. Cazalilla, author R. Citro,
  author T. Giamarchi, author E. Orignac,  and author M. Rigol, 10.1103/RevModPhys.83.1405 journal journal Rev.
  Mod. Phys. volume 83, pages 1405
  (year 2011), https://arxiv.org/abs/1101.5337 	arXiv:1101.5337 [cond-mat.str-el] NoStop
[Note4()]Note4
  note One has to be careful, however. Even though each
  off-diagonal term is of the order O(1/N), there are N - 1 of them.
  Therefore, the total contribution is of the order O(1). Furthermore, the
  very presence of the Bogoliubov mode is due to these off-diagonal elements,
  and hence neglecting them will result in losing the linear dispersion mode,
  which we will see immediately. Thus, strictly speaking, the following
  effective action accounts only for a part of the full set of possible
  dynamical excitations.Stop
[Note5()]Note5
  note Hereinafter, we assume, for simplicity, that all components
  are equally occupied, ρ ^(0) = N ρ ^(0)_a. Since in the large-N
  limit there is no coupling between different modes, this does not effect the
  following analysis.Stop
[Pich(2018)]Pich:2018ltt
  author author A. Pich, in @noop  booktitle Les Houches
  summer school: EFT in Particle Physics and Cosmology Les Houches, Chamonix
  Valley, France, July 3-28, 2017 (year 2018) http://arxiv.org/abs/1804.05664 arXiv:1804.05664 [hep-ph] NoStop

[Kamenev and Levchenko(2009)]Kamenev:2009
  author author A. Kamenev and author A. Levchenko, 10.1080/00018730902850504 journal journal Advances in Physics volume 58, pages 197 (year 2009), http://arxiv.org/abs/0901.3586 arXiv:0901.3586 [cond-mat.other]
  NoStop
[Lindner and Müller(2006)]Lindner:2005kv
  author author M. Lindner and author M. M. Müller, 10.1103/PhysRevD.73.125002 journal journal Phys. Rev. D volume
  73, pages 125002 (year 2006), http://arxiv.org/abs/hep-ph/0512147 arXiv:hep-ph/0512147 NoStop
[Lindner and Müller(2008)]Lindner:2007am
  author author M. Lindner and author M. M. Müller, 10.1103/PhysRevD.77.025027 journal journal Phys. Rev. D volume
  77, pages 025027 (year 2008), http://arxiv.org/abs/0710.2917 arXiv:0710.2917 [hep-ph] NoStop

[Note6()]Note6
  note Note that the summation goes from n=1, as opposed to the
  standard zero-T expansion, which reflects that W[0] ≡ 0 on the
  Schwinger-Keldysh contour.Stop
[Note7()]Note7
  note Recall that the expectation value of the phase field is
  zero, so that the n=1 term does not contribute.Stop
[Pawlowski(2007)]Pawlowski:2005xe
  author author J. M. Pawlowski, 10.1016/j.aop.2007.01.007 journal journal Ann. Phys. volume
  322, pages 2831 (year 2007), http://arxiv.org/abs/hep-th/0512261 arXiv:hep-th/0512261 NoStop
[Gasenzer and Pawlowski(2008)]Gasenzer:2008zz
  author author T. Gasenzer and author J. M. Pawlowski, 10.1016/j.physletb.2008.10.049 journal journal Phys. Lett. volume
  B670, pages 135 (year 2008), http://arxiv.org/abs/arXiv:0710.4627 [cond-mat.other] arXiv:0710.4627
  [cond-mat.other] NoStop

[Johnstone et al.(2018)Johnstone, Groszek, Starkey,
  Billington, Simula, and Helmerson]Johnstone2018a.arXiv180106952
  author author S. P. Johnstone, author A. J. Groszek, author P. T. Starkey, author C. J. Billington, author T. P. Simula,  and author K. Helmerson, @noop  journal journal
  ArXiv e-prints  (year 2018), http://arxiv.org/abs/1801.06952v2 arXiv:1801.06952v2 [cond-mat.quant-gas]
  NoStop
[Note8()]Note8
  note At finite times, on the way towards the fixed point, the
  integral only includes the long-wavelength modes up to a scale k_Λ(t) evolving with the distribution as k_Λ(t)∼ t^-β and
  is expected to be approximately conserved. Furthermore, as we showed, the
  transport shifts relative-phase excitations to the IR which only in the N→∞ limit decouple from each other and correspond to single-particle
  excitations in each component.Stop
[Schachner et al.(2017)Schachner, Piñeiro Orioli, and Berges]Schachner:2016frd
  author author A. Schachner, author A. Piñeiro Orioli,  and author J. Berges, 10.1103/PhysRevA.95.053605
  journal journal Phys. Rev. A volume 95, pages 053605 (year 2017),  http://arxiv.org/abs/1612.03038 arXiv:1612.03038
  [cond-mat.quant-gas] NoStop

[Janssen et al.(1989)Janssen, Schaub, and Schmittmann]Janssen1989a
  author author H. K. Janssen, author B. Schaub,   and author B. Schmittmann,  10.1007/BF01319383 journal journal Z. Phys. B Cond. Mat. volume 73,  pages 539 (year 1989)NoStop
[Janssen(1992)]Janssen1992a
  author author H. Janssen, in https://doi.org/10.1142/9789814355872_0007
  booktitle From phase transitions to chaos (publisher World Scientific, Singapore, year 1992)  p. pages 68NoStop
[Calabrese and Gambassi(2002)]Calabrese2002a.PhysRevE.65.066120
  author author P. Calabrese and author A. Gambassi, 10.1103/PhysRevE.65.066120 journal journal Phys. Rev. E volume
  65, pages 066120 (year 2002), http://arxiv.org/abs/cond-mat/0203096 arXiv:cond-mat/0203096
  [cond-mat.stat-mech] NoStop
[Gambassi(2006)]Gambassi2006a.JPAConfSer.40.2006.13
  author author A. Gambassi, http://stacks.iop.org/1742-6596/40/i=1/a=002
  journal journal J. Phys. Conf. Ser. volume 40, pages 13 (year
  2006)NoStop
[Note9()]Note9
  note The scaling limit of the non-thermal fixed point is reached at late evolution
  times provided that the initial condition is tuned such that no competing
  processes can prevent this. For a given system size this is typically ensured
  by a sufficiently extreme out-of-equilibrium initial state, such as the box
  distribution used here.Stop
[Blakie et al.(2008)Blakie,
  Bradley, Davis, Ballagh,  and Gardiner]Blakie2008a
  author author P. B. Blakie, author A. S. Bradley,
  author M. J. Davis, author R. J. Ballagh,  and author C. W. Gardiner, 10.1080/00018730802564254 journal journal Adv.
  Phys. volume 57, pages 363 (year 2008), http://arxiv.org/abs/arXiv:0809.1487v2
  [cond-mat.stat-mech] arXiv:0809.1487v2 [cond-mat.stat-mech] NoStop

[Note10()]Note10
  note The transport associated with a non-thermal fixed point in
  the N=1 case has recently been associated with a `transport peak' occurring
  in the spectral decomposition of the statistical correlation function,
  independent of the thermal Bogoliubov excitations in the system <cit.>.Stop
[Piñeiro Orioli and Berges(2019)]PineiroOrioli:2018hst
  author author A. Piñeiro Orioli and author J. Berges, 10.1103/PhysRevLett.122.150401 journal journal
  Phys. Rev. Lett. volume 122, pages
  150401 (year 2019), http://arxiv.org/abs/1810.12392 arXiv:1810.12392 [cond-mat.quant-gas]
  NoStop

[Mazeliauskas and Berges(2019)]Mazeliauskas:2018yef
  author author A. Mazeliauskas and author J. Berges, 10.1103/PhysRevLett.122.122301 journal journal Phys. Rev. Lett. volume 122, pages 122301 (year 2019),  http://arxiv.org/abs/1810.10554 arXiv:1810.10554 [hep-ph]
  NoStop

[Moore(2016)]Moore:2015adu
  author author G. D. Moore, 10.1103/PhysRevD.93.065043 journal
  journal Phys. Rev. D volume 93,  pages 065043 (year 2016), http://arxiv.org/abs/1511.00697 arXiv:1511.00697 [hep-ph] NoStop

[Gasenzer et al.(2012)Gasenzer, Nowak, and Sexty]Gasenzer:2011by
  author author T. Gasenzer, author B. Nowak,   and author D. Sexty, 10.1016/j.physletb.2012.03.031 journal journal Phys. Lett. volume B710, pages 500 (year 2012), http://arxiv.org/abs/1108.0541 arXiv:1108.0541 [hep-ph] NoStop

[Calzetta and Hu(2008)]Calzetta:2008iqa
  author author E. A. Calzetta and author B.-L. B. Hu, 10.1017/CBO9780511535123 title
  Nonequilibrium Quantum Field Theory, Cambridge Monographs on
  Mathematical Physics (publisher CUP, Cambridge, UK,  year 2008)NoStop

[Sieberer et al.(2016)Sieberer, Buchhold, and Diehl]Sieberer_2016
  author author L. M. Sieberer, author M. Buchhold,
   and author S. Diehl, 10.1088/0034-4885/79/9/096001 journal journal Rept. Progr. Phys. volume
  79, pages 096001 (year 2016), http://arxiv.org/abs/1512.00637 arXiv:1512.00637 [cond-mat.quant-gas]
  NoStop
[Note11()]Note11
  note We changed the notation of the mean density to n_a^(0)
  for this section to avoid a possible confusion with the spectral
  function.Stop
[Note12()]Note12
  note The following statements should, of course, be understood in
  terms of distributional derivativesNoStop
[Branschädel and Gasenzer(2008)]Branschadel:2008sk
  author author A. Branschädel and author T. Gasenzer, 10.1088/0953-4075/41/13/135302' journal journal J. Phys. B volume
  41, pages 135302 (year 2008), http://arxiv.org/abs/0801.4466 arXiv:0801.4466 [cond-mat.other]
  NoStop

[Berges and Borsanyi(2006)]Berges:2005md
  author author J. Berges and author S. Borsanyi, 10.1103/PhysRevD.74.045022 journal journal Phys. Rev. D volume
  74, pages 045022 (year 2006), http://arxiv.org/abs/hep-ph/0512155 arXiv:hep-ph/0512155 NoStop

[Note13()]Note13
  note Recall that F(x,y) is even, whereas ρ (x,y) is odd
  with respect to relative coordinates, which implies that the Fourier
  transform of the former is real, while the latter one is purely
  imaginary.Stop
[Note14()]Note14
  note Note that, for a non-local theory as given by the vertices
  depicted in Fig. <ref>, the validity of the on-shell approximation
  needs to be verified after all. This question is beyond the scope of the
  present work and is the subject of further study.Stop
[Yukawa and Ueda(2012)]Yukawa2012a.PhysRevA.86.063614
  author author E. Yukawa and author M. Ueda,  10.1103/PhysRevA.86.063614 journal journal Phys. Rev. A volume 86, pages 063614 (year 2012), http://arxiv.org/abs/1209.5509 arXiv:1209.5509 [cond-mat.quant-gas]
  NoStop
[Di et al.(2010)Di,
  Wang, and Wei]Di2010a
  author author Y. Di, author Y. Wang,  and  author H. Wei, 10.1088/1751-8113/43/6/065303 journal journal J.
  Phys. A: Math. Theor. volume 43, pages
  065303 (year 2010)NoStop
[Hamley et al.(2012)Hamley, Gerving, Hoang, Bookjans, and Chapman]Hamley2012a.NatPh.8.305H
  author author C. D. Hamley, author C. S. Gerving, author T. M. Hoang, author E. M. Bookjans,  and author M. S. Chapman, 10.1038/nphys2245 journal
  journal Nature Physics volume 8,  pages 305 (year 2012), http://arxiv.org/abs/1111.1694 arXiv:1111.1694 [cond-mat.quant-gas]
  NoStop








