Unraveling the nature of universal dynamics in $O(N)$ theories

Many-body quantum systems far from equilibrium can exhibit universal scaling dynamics which defy standard classification schemes. Here, we disentangle the dominant excitations in the universal dynamics of highly-occupied $N$-component scalar systems using unequal-time correlators. While previous equal-time studies have conjectured the infrared properties to be universal for all $N$, we clearly identify for the first time two fundamentally different phenomena relevant at different $N$. We find all $N\geq3$ to be indeed dominated by the same Lorentzian ``large-$N$'' peak, whereas $N=1$ is characterized instead by a non-Lorentzian peak with different properties, and for $N=2$ we see a mixture of two contributions.

The study of these far-from-equilibrium universality classes in isolated systems has so far primarily focused on the properties of equal-time momentum distribution functions, f (t, p).The typical scenario is depicted in Fig. 1a.Starting with high occupation numbers, which may be obtained, e.g., from instabilities or strong cooling quenches [24,25], the system quickly approaches an attractor solution characterized by self-similar scaling, f (t, p) = t α f S (t β p), also referred to as non-thermal fixed point.During this phase, the evolution is determined by the universal exponents α, β, and the universal function f S , which are largely insensitive to system parameters and details of the initial conditions.
Scalar field theories with O(N ) symmetry have been shown to exhibit such universal dynamics.In the infrared, they are characterized by α = βd and β ≈ 1/2 in d spatial dimensions [12].The physics is linked to particle number transport towards low momenta and the growth of a zero-mode condensate.Remarkably, previous works have found α, β and the form of f S to be universal for all values of N (see Fig. 1b), including both relativistic O(N ) and non-relativistic U (N ) theories describing ultracold Bose gases [12][13][14][15][16].The origin of this univer- sality has remained so far a mystery.Both exponents and scaling function have been successfully calculated using a large-N kinetic theory, which describes elastic collisions of quasiparticles with free dispersion and a renormalized interaction [12][13][14][15]30].However, it is unclear if and why this description should apply at small N as well.At the same time, descriptions based on defects, e.g.vortices, have provided alternative explanations of related models at small N [10,11,31,32].
In this Letter, we resolve this long-standing puzzle on the universality observed in O(N ) scalar theories using instead unequal-time (two-point) correlation functions.By providing information on both occupancies and dispersion relations, these observables allow to identify the dominant far-from-equilibrium excitations [33,34], and further provide access to the universal dynamical exponent z [35][36][37][38][39][40].As our main result, we find that the previously believed N -universality actually breaks up into (at least) two clearly distinct universality classes characterized by different phenomena, which are, however, almost indistinguishable from equal-time correlators and even the dynamical critical exponent z alone.O(N ) theories and statistical function.-Weconsider an O(N )-symmetric scalar field theory for relativistic scalar fields ϕ a (t, x), a = 1, . . ., N , in d = 3 spatial dimensions with classical action (c = k B = = 1) (1) Here, t,x ≡ dt d 3 x and sum over repeated indices is implied.We consider a weak coupling λ 1 and different values of m.However, since field fluctuations generate an effective mass M > m, the exact value of m is not relevant for the infrared physics discussed in this work.
We focus first on the unequal-time statistical function, defined for a translation invariant system as the anticommutator expectation value where 'c' denotes the connected part.Introducing the center and relative time coordinates, τ ≡ (t + t )/2 and ∆t = t − t , we Fourier transform it according to F (τ, ω, p) ≡ dx d∆t e i(ω∆t−px) F (t, t , x).
The statistical function can be seen as an unequal-time generalization of the distribution function, which at low momenta p M is given by f (t, p) ≈ F (t, t, p) M [41].F contains not only information about the occupancy of excitations in the system, but also about their frequency dependence.Thus, the information contained in F can be crucial to understand which excitations dominate the dynamics of a system.
We consider far-from-equilibrium initial conditions with large (Gaussian) fluctuations up to a characteristic scale Q as given by f (t = 0, p) = n0 λ Θ(Q − p).Due to the initial 'overoccupation' of mode excitations around Q, the subsequent redistribution dynamics is dominated by transport of particles to lower momenta, and is characterized by universal scaling in the infrared as explained in the introduction.
To describe the system we employ classical-statistical simulations (TWA), which are justified in the limit of high occupancies and small couplings as considered here [42][43][44][45].We perform large-scale simulations averaging over up to 40 runs using n 0 = 100, and give all dimensionful quantities in units of Q.We use either m = 0 or m = 0.5 and extract M from our data.F is computed from a classical correlation function, and the relativetime Fourier transforms are performed using standard signal-processing methods.If not stated otherwise, we show data for τ = 1000.Further details are given in the Supplementary Material [46].Large-N peak.-Weconsider first the dynamics in the large-N limit.In general, the statistical function F exhibits several peaks at different frequencies.However, for large N , we find the signal to be clearly dominated by one single contribution, as shown for N = 8 in Fig. 1e.We refer to it as the 'large-N peak'.This peak is depicted for fixed p in the inset (red points), including error bars.It is well-described by a Lorentzian parametrized as which is shown as a black dashed line and accurately agrees with the data.The results of the fitting procedure with Eq. ( 3) are shown in the right column of Fig. 2 for N = 8.The dispersion relation (Fig. 2b) is time-independent and accurately agrees with a free-particle dispersion of the form ω largeN (p) = p 2 + M 2 .Hence, the dispersion scales quadratically at low momenta as ω largeN −M ∼ p 2 /(2M ), implying a dynamical critical exponent z = 2.
The decay rate, on the other hand, decreases with time and depends approximately linearly on momentum at low p (Fig. 2d).As shown in the inset, its time evo-lution obeys self-similar scaling given by γ largeN (τ, p) = τ −βz γ largeN,S (τ β p) with β = 1/2 and z = 2, consistent with the dispersion.Interestingly, we find that the lifetime of the large-N quasiparticles grows with time since γ largeN (τ, p) → 0, such that the form of the large-N peak approaches a δ-function.
The dominance of this peak implies F largeN (τ, ω, p) ≈ F (τ, ω, p) and A largeN (τ, p) ≈ f (τ, p)/M .Together with the scaling behavior of the dispersion and width, this leads to a self-similar evolution To understand the nature of the large-N excitations, we study the behavior of fluctuations with the classical equation of motion [46] for details).This equation has stable rotating solutions parametrized by ϕ 0 (t, x) = | ϕ 0 |e(t), where ë(t) = −M 2 e(t).They correspond to rotations in a two-dimensional hyperplane in ϕ-space, which have been observed numerically [16].By studying linear fluctuations around these solutions, we find two in-plane (phase and radial) excitations, and a set of ∼N outof-plane excitations which correspond to rotations perpendicular to the plane spanned by e(t) and ė(t).The out-of-plane excitations have a p 2 + M 2 dispersion and dominate in the large-N limit.These excitations correspond to the large-N peak discussed here.
N = 1 non-Lorentzian peak.-Weconsider now the statistical function F of a single-component theory as shown in Fig. 1d.At low momenta p 0.4 it is again dominated by a single peak [47], shown for fixed p in the inset (red points).While a Lorentzian fit (blue dashed line) with (3) fails to capture the tails of the peak, we find it to be phenomenologically well described by (black dashed curve) where the subscript 'nL' stands for non-Lorentzian.We employ this form as a fit function to extract the properties of this peak leading to the results in the left column of Fig. 2. The dispersion relation ω nL (τ, p) (Fig. 2a) is approximately linear and obeys a self-similar scaling form ω nL (τ, p) − M = τ −βz ωS (τ β p) with exponents β = 1/2 and z = 2, as shown in the inset.The decay rate γ nL (p) (Fig. 2c) is found to be instead time-independent and to scale as γ nL ∼ p z .This implies a self-similar evolution as in Eq. ( 4), notably with the same scaling exponents as for the large-N peak and Remarkably, we find the properties of this non-Lorentzian peak to be identical to the infrared peak found in Ref. [34] for a non-relativistic U (1) complex scalar theory.This can be explained by noticing that at small momenta, p M , particle number changing processes are suppressed and the O(1) theory is hence described by an emergent non-relativistic U (1) theory.The mapping can be made more rigorous by defining the non-relativistic -Which of these two distinct physical phenomena, if any, dominates for intermediate N at low momenta is a priori not obvious.
We find that the N = 2 theory is a special case.The statistical function F appears to have two distinct contributions with a similar weight at low momenta, as shown in Fig. 3a.The inset shows a fit (black dashed line) to these peaks with the sum of both functions (5) and (3).Each of these peaks is included in the inset as separate dashed lines (green and yellow), and their respective dispersions are shown in the main plot.We find that the dispersion of the left peak agrees with ω nL of the non-Lorentzian peak in O(1) theory, while the dispersion of the right peak obeys approximately ω largeN of the large-N peak.However, at low momenta both contributions overlap so strongly that they become almost indistinguishable.Therefore, their p → 0 asymptotic behavior is hard to extract.
For N ≥ 3 we find the infrared dynamics to be dominated by the large-N peak.The statistical function for N = 3 is shown in Fig. 3b and for fixed momentum in the inset.The dominant peak has a dispersion that agrees with ω largeN (p) (gray dashed line), and we confirmed that the width shows a similar behavior as in the large-N limit.We note, however, that for small N we find evidence of an additional contribution overlapping with the main peak at lower frequencies.Based on the above, this is possibly related to a non-Lorentzian contribution, which appears to quickly disappear as N or momentum increases.
Spectral peaks, we study the spectral function defined as At equal times, t = t , it is determined by the equaltime commutation relations, ρ| t=t = 0 and ∂ t ρ| t=t = δ(x − x ).For unequal times, this quantity encapsulates the linear response of the system to perturbations and thus contains information about the low-lying excitations of the system.To compute it we employ a linear response approach as described in Refs.[33,34] (see [46] for details).
In Fig. 4 we show color plots of the spectral function for N = 1 and N = 8.In general, we find that ρ shows the same peak structure as F but with different relative weights between the peaks.However, since ρ does not contain information about occupancies, the weight of the peaks does not reveal the dominant excitations.
In particular, for N = 1 we find that the non-Lorentzian peak in ρ has a very small weight [49] (upper inset of Fig. 4a), which becomes visible only at low momenta.Nevertheless, this peak has the same dispersion and width as the peak in F (blue points in Figs.2a,c).In fact, we find that the non-Lorentzian peaks in F and ρ fulfill a generalized fluctuation-dissipation relation (FDR) given by We show this in the lower inset of Fig. 4a by filtering out all peaks but the non-Lorentzian peak [46].Here, µ ≡ M is an effective chemical potential linked to the approximate conservation of particle number at low momenta.Equation ( 7) is reminiscent of the equilibrium FDR [50], except with a mode-dependent temperature T nL (τ, p).In the large-N limit, the spectral function is dominated by the large-N peak (Fig. 4b) with a dispersion and width which also match the results for the corresponding large-N peak in F (Figs. 2b,d).The two peaks can again be related through a generalized FDR as in Eq. ( 7) (lower inset of Fig. 4b).However, since the large-N peak dominates in both F and ρ, and since its width becomes narrower with time, Eq. ( 7) can be simplified in the long-time, large-N limit to This is again similar to the thermal case, except with a time-dependent non-equilibrium distribution, and is reminiscent of kinetic approximations [51].At intermediate N ≥ 3 we find small deviations from this behavior, which vanish as N → ∞.
Discussion.-Our results reveal the existence of (at least) two distinct infrared universality classes governed by two different types of phenomena, despite being characterized by very similar universal exponents.
For N ≥ 3, the dynamics is dominated by a Lorentzian large-N peak.It results from a set of excitations with relativistic dispersion p 2 + M 2 which dominate due to their number scaling with N .The fact that this physics dominates even at lower N is, however, remarkable, especially for N = 3 where only one such excitation exists.Apart from this, given that the lifetime of these quasiparticles grows with time, and the fact that they fulfill the generalized FDR of Eq. ( 8), our results validate the analytic large-N kinetic theory used in Refs.[12][13][14].
For N = 1, the dynamics is instead dominated by a peak of non-Lorentzian shape with time-dependent dispersion and time-independent quadratic width, which coincides with the findings for the non-relativistic U (1) theory [34].This hints to a common origin for these infrared excitations.We suggest that this non-Lorentzian peak in O(1) corresponds to vortex line excitations in effective non-relativistic degrees of freedom.This claim is motivated by previous works on U (1) dynamics [8,10], and O(1) dynamics in d = 2 [32], all of which saw evidence of vortex excitations in real-space snapshots of the field.
Our results for N = 2 show a mixture of two different contributions.We found evidence that one contribution has a p 2 + M 2 dispersion, analogously to the large-N peak.Analytically, however, we do not find any excitation with such a dispersion for the O(2) theory [46].The second contribution was found to share the same properties as the non-Lorentzian peak of O(1).Thus, we hypothesize that these two peaks possibly originate from vortex-type excitations and domain walls, such as those observed in Ref. [11].
Conclusion.-In this work, we have disentangled the physical origin of the infrared universal dynamics of O(N ) scalar theories.We have identified two distinct universality classes as a function of N : non-Lorentzian excitations for N = 1 and Lorentzian rotational excitations for N ≥ 3, while for N = 2 we find a mixture.In a broader context, our work shows the importance of the unequal-time statistical function F to reveal the dominant physical phenomena in far-from-equilibrium systems, and will potentially trigger future research into this direction.In particular, measuring unequal-time functions [52,53] could substantially improve our understanding of running cold-atom experiments.

Supplementary Material:
Unraveling the nature of universal dynamics in O(N ) theories

TECHNICAL IMPLEMENTATION
Here we provide more details on our simulations.As explained in the main text, we consider O(N )-symmetric relativistic scalar theories given by the action Because of the large occupation numbers at low momenta in our initial conditions quantum vacuum contributions are suppressed by a factor of λ relative to f (0, p).In the weak-coupling limit λ → 0 with λf (t, p) held fixed, the original quantum field theory is accurately mapped onto a classical-statistical field description [42][43][44].In particular, the dynamics is not governed by the lattice cutoff but by physical scales like Q instead.We ensured that our results are insensitive to descretization parameters like the lattice spacing or the volume by varying them.
In the numerical approach, classical fields ϕ a (t, x) and their conjugate momenta π = ∂ t ϕ are discretized on a lattice with lattice spacing a s and volume V = (N s a s ) 3 .Unless stated otherwise, we used N s = 256 and a s = 0.8 in units of Q.The above initial conditions (S1) are implemented by sampling the field as with independent Gaussian random numbers fulfilling c a (p) (c b (q)) * cl = V δ p,q δ a,b and similarly for c.The fields are evolved by solving the classical Hamilton equations of motion, where we use a t = 0.05a s for the time step.To reduce lattice artefacts, we use a fourth order discretization scheme for second order spatial derivatives, as described in Ref. [7].
In the framework of classical-statistical simulations, the statistical function F can be straightforwardly computed as [51] where • cl denotes average over initial conditions and over the direction of p.
To compute the spectral function we employ a linear response approach similar to Refs.[33,34].In essence, we perturb the system as ϕ a → ϕ a + δϕ a with a source j a (t, x) = j 0 a (x)δ(t − t ) at time t and time evolve the perturbation δϕ a according to the linearized equations of motion ∂ t δϕ a = δπ a and We choose to perturb with random Gaussian source fields fulfilling j 0 a (p)(j 0 b (q)) * j = V δ p,q δ ab , where • j denotes average over the perturbations.The retarded part of the spectral function results from linear response theory as Our choice of j 0 a (p) ensures that the canonical commutation relations lim t→t ρ(t, t , p) = 0 and lim t→t ∂ t ρ(t, t , p) = 1 are implemented correctly.
After Fourier transforming with respect to relative time ∆t = t − t , the correlations become with central time τ ≡ (t + t )/2.We used that F (ρ) is even (odd) in ∆t, such that their Fourier transforms are real-valued functions of ω.Strictly speaking, these symmetries only hold for fixed τ .Because of the slow τ -dependence of our observables and for the considered time windows of ∆t, we find the transforms at fixed t to accurately approximate the ones at fixed τ .In practice, we hold t fixed and compute the transforms as where we typically employ ∆t max ≈ 300 − 500.To smoothen the resulting curves, we use standard signal processing techniques by employing the Hann window function and zero-padding.The latter effectively amounts to evaluating (S8) and (S9) at more intermediate frequencies than provided by the usual discrete Fourier transform.We checked that these techniques do not alter the peak structure nor the peak forms considerably, but mainly reduce background ringing.Correlation functions at fixed momentum p are shown in the main text with error bars.To obtain the statistical errors, we first Fourier transform our data for each run and subsequently average over the results.
Because of the effective mass M > m, generated by fluctuations, the exact value of m is not relevant for the infrared physics discussed in this work.In particular, the correlation functions F and ρ plotted as functions of ω − M are independent of m at low momenta p < M .For the N ≥ 2 data shown we used m = 0, but we have also checked that m > 0 does not change the results.For N = 1, the effective mass M decreases as a power law in time for m = 0, whereas it stays approximately constant for sufficiently large nonzero m [7,16,44].To avoid this issue, we use m = 0.5 in all figures with N = 1.However, we have checked that simulations with m = 0 provide the same peak structure and the same properties of the non-Lorentzian peak as for m > 0.

DATA ANALYSIS
To extract properties of the peaks, we perform fits to the spectrum with different fit forms.We distinguish between Lorentzian and non-Lorentzian peaks, where we use a hyperbolic secant form for the latter.The fit forms for these peaks are We employ these functions to fit the positive frequency peaks (ω 0 > 0) of F and ρ.Due to symmetries, the same peaks appear as well at negative frequencies ω 0 → −ω 0 .However, these negative-frequency contributions are irrelevant for our fits since we have γ/ω 0 1 for all peaks and thus, it is sufficient to perform fits in positivefrequency space.
The corresponding fit functions in ∆t are L(∆t; ω 0 , γ) = cos (ω 0 ∆t) e −γ|∆t| , (S13) S(∆t; ω 0 , γ) = cos (ω 0 ∆t) sech(γ∆t) . (S14) In general, we find that fits in ω are easier to perform and more likely to converge than fits in ∆t.However, the latter are useful when the peak width is very narrow, as is the case for the large-N peak in N = 8 (Fig. 2d).To extract the properties of such narrow peaks, we first performed fits in ω space with L(ω; ω 0 , γ) and then used the results as initial values for fits in ∆t with L(∆t; ω 0 , γ).
The properties of the dominating peak for N = 1 were obtained in frequency space.For F , we use (S12) to fit F for p < 0.15.At higher momenta, other peaks start to become visible, as one finds in Fig. 1d.The main additional peaks correspond to a relativistic generalization of Bogoliubov excitations (see next section) with a Lorentzian shape.While additional peaks will be studied in detail in Ref.
[54], the fit function needs to include them at larger momenta, because peaks start to overlap.Hence, to fit F for p ≥ 0.15 we use with the fit parameters ω i , γ i and F i .Here F i is the weight of each peak, and the sum of all weights equals F (t, t, p).However, since S dominates at low momenta one has F nL ≈ F (t, t, p).The insets in Figs.1d and 4a (bottom inset) show the non-Lorentzian part of F (τ, ω, p) (denoted in the plots as F nL ).We obtain this after subtracting the fit results of the Bogoliubov peaks [F 1 and F 2 peaks in (S15)] from the original signal.
For the spectral function ρ, Bogoliubov excitations are always important for N = 1, as can be seen in Fig. 4a.We still find the same non-Lorentzian peak at low momenta, which appears, however, extremely small and hard to detect without signal processing techniques.To extract its properties, we performed a similar fit as above with where the fit parameters are ω i , γ i , A i and µ that turns out to be µ ≈ M .The corresponding results are shown in the left panel of Fig. 2 as blue points.Bogoliubov excitations are also present for other values of N .The corresponding peaks generally appear in two branches with dispersions M ± ω Bog (p) symmetric around M .Therefore, we use the mean value of the two dispersions to calculate M .For N = 8, these excitations are barely visible and we use instead the relativistic dispersion M 2 + p 2 to obtain M .We checked that both methods, where applicable, provide a consistent value for M .
We compute in this section the dispersion of the large-N excitations discussed in the main text.For this we consider the O(N ) scalar theory defined in Eq. (1) for N ≥ 2. The corresponding classical equation of motion reads This equation supports homogeneous solutions of the form where ρ 0 > 0, ϕ 0 is an N -component vector in ϕ-space, and the unit vector e(t) fulfils The effective mass is given by Solutions to Eq. (S19) can be easily constructed as where ω 0 ≡ ±M , |1 = (1, 0, 0, . ..)T , |2 = (0, 1, 0, . ..)T , and R ∈ O(N ) is an arbitrary orthogonal matrix, RR T = 1.These solutions correspond to rotations on an arbitrary but fixed two-dimensional hyperplane in ϕ-space that passes through the origin.
In the following, we study the nature of excitations around solutions of the type (S18).For this we add fluctuations as ϕ(t, x) = ρ 0 + δρ(t, x) e Λ ab θ ab (t,x) e(t). (S22) Here, it is assumed that δρ ρ 0 and |θ ab | 1.The matrices Λ ab correspond to the generators of O(N ) rotations.Expanding this expression to linear order in the fluctuations we obtain ϕ(t, x) = ϕ 0 (t) + δ ϕ(t, x) with δ ϕ(t, These correspond essentially to Pauli iσ y rotations in the plane defined by the axes a, b.Using Eq. (S25), the action of the generators on e and ˙ e can be straightforwardly computed.In this way, we obtain from Eq. (S24) a sum of terms proportional to the vectors e, ˙ e, and |b (b ≥ 3), which are all orthonormal to each other.The resulting linearized equation of motion can then be projected onto the different axes defined by these vectors, leading to a total of N equations, which we analyze in the following.

In-plane excitations
We start by projecting Eq. (S24) onto e and ˙ e.This leads to a system of equations for δρ and θ 12 given by The variables δρ and θ 12 correspond to radial and phase excitations in the plane spanned by ϕ 0 (t) and ∂ t ϕ 0 (t).

(S28)
Here, we defined the effective interaction coefficient These excitations correspond to a Bogoliubov-type phase excitation ω Bog (p) ≡ ω− which becomes linear in p at low momentum, and a massive radial excitation ω rad (p) ≡ ω+ with ω rad (p = 0) = 2 M (M + gρ 0 ).Note that in numerical simulations we look at correlation functions in the original ϕ basis, instead of in the phase and radial variables.Because of this, the dispersions numerically presented in the main text for δ ϕ fluctuations are shifted with respect to the above solutions

FIG. 1 .
FIG.1.Top: a) Typical thermalization scenario for large initial occupancies, the system approaches an attractor during its evolution.b) The distribution function f (t, p) at low momenta close at the attractor for N = 1, 2, 3, 4, 8 at times t = 600, 600, 750, 1000, 2000, respectively.c) Illustration of universality classes for different N .Bottom: Statistical function F (τ, ω, p)/f (τ, p) for N = 1 (left, d) and N = 8 (right, e) at τ = 1000.Red lines correspond to F at fixed p, also shown in the insets.Black dashed lines are fits with Eqs.(5) and (3), respectively.FnL is shown after filtering out other small peaks.A (blue dashed) Lorentzian curve is included for comparison.

FIG. 3 .
FIG. 3. Statistical function F (τ, ω, p)/f (τ, p) for N = 2 (left, a) and N = 3 (right, b).The insets show F (red) for the fixed momentum marked by the vertical red lines.Left: green/yellow dashed lines show the peak positions extracted from fits, as well as the fit results in the inset.Right: a relativistic dispersion (gray) is added for comparison.

FIG. 4 .
FIG. 4. Spectral function |ρ| for N = 1 (left, a) and N = 8 (right, b).The upper insets show F/f (τ, p) and ρ for fixed p marked by vertical blue lines.The lower insets show the validity of Eq. (7) by comparing (ω − M )F and ρ, both normalized to weight one.The ratio of normalizations gives TnL(τ, p).Black/yellow dashed lines are fits to F , ρ.