Strongly anomalous non-thermal fixed point in a quenched two-dimensional Bose gas

Universal scaling behavior in the relaxation dynamics of an isolated two-dimensional Bose gas is studied by means of semi-classical stochastic simulations of the Gross–Pitaevskii model. The system is quenched far out of equilibrium by imprinting vortex defects into an otherwise phase-coherent condensate. A strongly anomalous non-thermal fixed point is identified, associated with a slowed decay of the defects in the case that the dissipative coupling to the thermal background noise is suppressed. At this fixed point, a large anomalous exponent η ≃ − 3 and, related to this, a large dynamical exponent z ≃ 5 are identified. The corresponding power-law decay is found to be consistent with three-vortex-collision induced loss. The article discusses these aspects of non-thermal fixed points in the context of phase-ordering kinetics and coarsening dynamics, thus relating phenomenological and analytical approaches to classifying far-from-equilibrium scaling dynamics with each other. In particular, a close connection between the anomalous scaling exponent η, introduced in a quantum-field theoretic approach, and conservation-law induced scaling in classical phase-ordering kinetics is revealed. Moreover, the relation to superfluid turbulence as well as to driven stationary systems is discussed.


Introduction
The concept of universality has been very successful in characterizing and classifying equilibrium states of matter. It rests upon the observation that symmetries can cause the same physical laws to describe certain properties of very different systems. For example, in the vicinity of a continuous phase transition, certain macroscopic properties of a system become independent of the dynamical details, such as the microscopic interactions between its constituents. In the formulation of renormalization-group theory [1,2], universality manifests itself in scaling laws characterized by a set of universal power-law exponents. These exponents are, to a large extent, determined by the symmetries present in the system and define which class the critical behavior belongs to. As a result, phenomena as diverse as opalescent water under high pressure, critical protein diffusion in cell membranes [3], or early-universe inflationary dynamics [4,5] fall into the same universality class.
Taking a more general perspective, the question arises to what extent universality appears also in the dynamics away from thermal equilibrium. A standard classification scheme of dynamical critical phenomena applies to the linear response of systems driven, in a stochastic way, out of equilibrium [6]. Going further beyond such near-equilibrium settings, the theory of phase-ordering kinetics focuses on universal scaling laws for the coarsening dynamics following a quench across a phase transition [7][8][9][10][11][12]. Closely related dynamical critical phenomena include (wave-)turbulence [13,14], as well as superfluid or quantum turbulence [15,16]. Universal scaling far from equilibrium has recently been analyzed for different types of quantum quenches [17][18][19][20][21][22][23][24][25][26][27][28], see also [29][30][31][32][33] for studies of phase-ordering kinetics in ultracold Bose gases. While coarsening phenomena have partly been associated with the standard dynamical universality classes [11], a rigorous renormalization-group 2 , 1 2 2 where g sets the interaction strength and m r = g is a chemical potential, with mean density ρ. All terms are local in space and time, and we have suppressed the corresponding arguments of the bosonic field operators Φ. We focus on an isolated, homogeneous 2D system of size L×L with periodic boundary conditions. We will also consider the case with coupling to a thermal bath. Details of the description concerning the driving and dissipation relevant to that case are given in appendix C. 1. In the physical situations we are studying in this work, the relevant mode occupation numbers are large. Hence, the time evolution of the quantum many-body system, described by the Hamiltonian(1), can be efficiently simulated with semi-classical numerical techniques. We make use of the Truncated-Wigner approximation [72,73] in which the quantum operator F( ) x t , is replaced by a classical stochastic field variable f ( ) x t , in equation (1) and the time evolution of each realization is determined by the classical Gross-Pitaevskii equation (GPE) of motion, A state is then represented by a statistical ensemble of fields f ( ) x t , and quantum expectation values of an operator O are computed, in a Monte-Carlo-type fashion, via ensemble averages with respect to a positive definite phase-space (Wigner) distribution function [74], Here  is the Wigner function at the initial time t i and O c is the classical (Weyl) symbol of O. The delta distribution together with the functional measure encodes the time evolution of  to any time ¹ t t i . As we are studying homogeneous systems with periodic boundary conditions, one of our natural observables will be the (angle-averaged) single-particle momentum spectrum where the F( ) k t , are the bosonic field operators in momentum space, evaluated at time t. Within the semiclassical approach, defining the initial state is tantamount to defining the statistical properties of the initial distribution of fields, t 0 . The initial distribution  = | t 0 can be assumed to describe Gaussian fluctuations of the fields, separate for every momentum mode f k , * f k . Thus, it is sufficient to define the mean The non-vanishing variance (5) accounts for the vacuum fluctuation of half a particle on average for a noninteracting mode at zero temperature, and is subtracted again before interpreting n as a mode occupation. Given the above initialization of the state the Truncated-Wigner approximation is expected to quantitatively well describe the early and intermediate-time evolution relevant in this work.

Universal time evolution
3.1. Non-thermal fixed points in a superfluid The formation of turbulent ensembles of vortices in an isolated, dilute, superfluid Bose gas [45,46,48] can be associated with the system approaching a non-thermal fixed point [34][35][36] where the time evolution is critically slowed down before the vortices have mutually annihilated and the system eventually thermalizes [38,47]. Turbulent behavior and critical slowing down are signaled by characteristic scaling properties of correlation functions such as the angle-averaged single-particle momentum spectrum defined in equation (4) 5 .
An initial overpopulation [48,75] of momenta at scales on the order of the inverse healing length r = x k g m 2 , induced, e.g., by an instability, can subsequently drive the system to a non-thermal fixed point [76]. This happens because the particles in the overpopulated region are transported in momentum space to modes with lower energy. Close to a non-thermal fixed point, this inverse transport is generically self-similar in space and time, obeying the scaling relation , . 6 0 0 0 At the same time, energy conservation forces a few particles to scatter into higher momentum modes, eventually forming an incoherent thermalized fraction of the gas. The scaling exponents α and β in equation (6) have been calculated for the nonlinear Schrödinger model [38]. The predictions are derived by means of a scaling analysis of Dyson-type field dynamic equations obtained with 2-particle-irreducible effective-action techniques [76].
Assuming a quasiparticle description to apply, these equations lead to effective quantum Boltzmann equations, incorporating a non-perturbative momentum-dependent many-body scattering matrix. A crucial role is taken by the scaling of this scattering matrix. It is obtained within a non-perturbative scheme, summing diagrams to all orders in the bare coupling g [34][35][36][37], equivalent to a next-to-leading order large- approximation for the case of an  ( ) O -symmetric model. The effective quantum Boltzmann equations take into account the scattering of quasiparticle modes the properties of which are encoded by a spectral function During the scaling evolution, this function is assumed to be stationary. It can thus be written as the expectation value of the commutator matrix of the bosonic field operators , defining the frequency dependent response of eigenmode k. At the non-thermal fixed point, the spectral function is assumed to scale as z denotes the dynamical critical exponent. For d-dimensional systems close to their upper critical dimension,  d d up , the anomalous exponent η is expected to be small. Note, however, that η at a non-thermal fixed point does not need to be equivalent to the static anomalous dimension at an equilibrium phase transition. Below, we present an example, in d=2, where η is distinctly different from zero, related to an anomalously small temporal scaling exponent β in equation (6).
For the GP model (1), the exponents α and β read [38] under the condition of particle number conservation, within the region of self-similar transport 6 . Setting the anomalous dimension η to zero, one obtains the set of exponents a = d 2, and b = 1 2 which have been numerically confirmed for the (three-dimensional) Bose gas in [38]. If b > 0 the self-similar build-up (6) of momentum occupations in the infrared reflects a coarsening process. This coarsening can be interpreted in terms of the dynamics and annihilation of vortex defects created through the initial quench [46,47]. The vortices are moving through the superfluid gas, interacting with each other and with background sound waves. This allows for a dilution of the defects through mutual annihilation of vortices with opposite circulation [46,47]. The coarsening evolution implies an algebraic growth law, for the characteristic mean spacing ℓ d between the vortices. Eventually, when the last defects have disappeared, long-range order is established and phase coherence of the Bose-condensed gas is maximized.
In the intermediate-wavelength region of the inverse transport,  ℓ k 1 d , the momentum distribution is found to exhibit characteristic scaling~z -( ) n k k in momentum space [34][35][36][37], e.g.,~n k k 4 in d=2 dimensions and for vanishing η. For  h 2, this power-law behavior is cut off in the infrared at a certain scale l k to ensure non-divergence of the overall particle density ρ, seeequation (10). For the case of vortices in a dilute Bose gas, Î { } d 2, 3 the power law (12), with h = 0, reflects the geometry of the superfluid circulation around the defects, implying a characteristic power-law fall-off of the velocity, f | | v r arg 1 , at distance r perpendicular to the vortex (line). The scaling (12) is also well known (for h = 0) as the Porod tail [11,77] of the momentum distribution of an ensemble of randomly distributed vortices, at > l k k , where the scale l k 1 is related to the mean distance between defects [11,46]. Here, we introduce the d z would apply, e.g., in the case of a self-similar transport conserving energy,  ò definition (12) as a natural extension of Porod's law to include an anomalous dimension η (see section 4.1 for a more detailed discussion).

Preparation of the initial state
In this work, we demonstrate that different non-thermal fixed points can be reached, depending on the particular choice of the initial non-equilibrium state. Universality implies that the dynamical evolution in the vicinity of the fixed point is insensitive to details of the initial state. On the other hand, it is reasonable that the initial state needs to fulfil certain conditions for the system to enter a universal regime of time evolution in the first place. We generate initial states by phase-imprinting vortex defects into a fully phase-coherent Bose gas [50,58]. This offers different parameters, such as the density of vortices, their winding number, and the distribution statistics, to vary the initial conditions.
A fully phase-coherent gas is prepared, given by homogeneous field configurations. In each sample prepared within the Truncated-Wigner scheme, quantum fluctuations are included in the empty modes. We choose, in particular, the mode occupation numbers according to equation (5), for all modes ¹ k 0, while the zero mode is populated with the total number of particles N. Then, the phase pattern q ( ) x t , i of the desired vortex distribution is multiplied into the sampled homogeneous field configurations, , e x pi , i i i . We study two types of initial vortex configurations, regular lattices of non-elementary defects with winding numbers > | | w 1, and uniform random distributions of elementary defects with winding number =  w 1. For both types, we choose distributions with equal numbers of vortices and anti-vortices. Vortex cores in the density f | ( )| x t , i 2 are subsequently formed by means of a short period of evolution according to the imaginary-time GPE, i.e.with t t i in equation (2). In figure 1(a), an example for such a lattice configuration with slightly displaced ( =  w 6)-vortices is shown, depicting the superfluid density field f | ( )| x t , i 2 after imaginary time evolution. The random displacements, in addition to the Truncated-Wigner noise, speed up their decay into elementary vortices. All numerical data is shown in units expressed in terms of the healing length x r = m g 1 2 h (with mass = m 1 2 set in the numerical calculations), see appendix A for more details.

Relaxational production of clustered vortex ensembles
Before analyzing the dynamics of occupation numbers ensuing the phase imprinting, we discuss the spatial dynamics on phenomenological grounds, at the level of single realizations of the order parameter field f ( ) x t , . We focus on an initial arrangement of the vortex defects on a rectangular lattice with alternating signs of the winding number =  w 6, see figure 1(a). During the early time evolution, these defects quickly decay into elementary vortices with = | | w 1 [78]. To speed up this otherwise rather slow decay process, we add initially small random shifts to the regular lattice positions of the defects (see figure 1(a)).
In a first stage, the like-sign elementary vortices form tightly-bound clusters, see figure 1(b), which play an important role in the ensuing dynamics. Each cluster screens the vortex-anti-vortex attraction [79], such that mutual annihilation of vortices and anti-vortices is suppressed. Figure 2 shows examples of the hydrodynamic velocity field , at the same times t as in figure 1. The clustering of like-sign vortices is indicated by (color encoded) extended and strong fluxes.
In later stages, vortices with opposite winding numbers begin to mingle with each other, overcome the screening and start to mutually annihilate. As a consequence, the vortex configuration undergoes a dilution process, see figure 1, panels (c) and (d), which leads to an ordering of the phase field [45][46][47][48]. We discuss the universal character of this ordering process in the next section.
As seen in figures 1(b)-(d), the superfluid density of a single realization develops a considerable amount of small-scale fluctuations. These fluctuations arise in the decay of the non-elementary vortices as well as in the mutual annihilations and are amplified by the nonlinear interaction term in equation (2). They play an important role for the dilution process, as they are known to mediate vortex interactions [29,70].

Anomalous non-thermal fixed point
Near a non-thermal fixed point, the time evolution of the single-particle spectrum ( ) k n t , , equation (4), is expected to correspond to a scaling transformation (6) of a universal scaling function. The occupation spectrum corresponding to the evolution starting from the vortex-lattice state described in the previous section is shown in figure 3. After an initial stage during which the non-elementary vortices decay, a power-law dependence of n on = | | k k builds up in an infrared momentum region. This steep power law levels off at a scale k λ in the infrared. This scale shifts in time towards lower momenta increasing the occupation in the infrared, while the occupation decreases at larger wave numbers. The evolution represents self-similar particle transport towards the infrared, building up a quasi condensate at low wave numbers, similarly as described for the 3D case in [80][81][82].
The self-similar nature of the time evolution of ( ) n k t , , depicted in figure 3, is demonstrated by the scaling collapse according to equation (6). Figure 4 shows that the time evolution of the occupation spectrum can indeed be rescaled to a single curve. The fitting procedure which is described in detail in appendix D yields the scaling exponents The index 'a' is chosen to distinguish the exponents from those obtained later with different initial conditions. In figure 4, the rescaled spectra are shown at five exemplary times between For the scaling analysis we took into account the occupation spectra at 300 logarithmically equidistant times within the same interval.
The scaling evolution shown in figure 3 can be interpreted to reflect turbulent (non-local inverse particle) transport. We note that this inverse transport does not conserve particle number locally, i.e., momentum shell After the initial build-up of a characteristic distribution, a self-similar shift towards smaller momenta is seen during the universal stage of the time evolution. This inverse transport process conserves particle number and, as compared to a turbulent cascade, is non-local in momentum space. The earliest time shown ( x =t 3 10 h 2 2 , gray points) is still within the non-universal stage of evolution right after the decay of the lattice (see figure 1(b)). The dynamics during the later stage of the evolution represents a rescaling in space and time and is analyzed in more detail in figure 4.  per momentum shell, as commonly required for (wave) turbulent cascades within an inertial interval of wave numbers. Nevertheless, in the isolated system we consider, the transport conserves the total particle number.
Our simulations give scaling exponents α and β consistent with the relation a b = d , seeequation (9).
Hence, particle number conservation also holds within the regime of momenta and times governed by the scaling evolution shown in figure 4, whereas a different relation between α and β would apply if the self-similar transport conserved energy or other quantities [38]. Based on the predictions (9) for the scaling exponents α and β we infer a strong anomalous scaling exponent to characterize the non-thermal fixed point reached from our vortex-lattice initial conditions.

(Near-)Gaussian fixed point
In contrast to the above anomalous scaling, characterized by the exponents (13) and (14), also a distinctly different self-similar evolution is possible when starting from different initial states. In the following, we present results for initial states containing random spatial distributions of 2400 elementary defects, with an equal number of vortices and anti-vortices. This corresponds to four times as many vortices as obtained from the decaying non-elementary vortices chosen in the previous section. In figure 5, the time evolution of the occupation spectrum is shown in analogy to figure 3. Figure 6 shows the scaling collapse within the same time frame as figure 4 in the anomalous case. As before, at late times, 10 h 4 2 , the spectrum evolves in a universal self-similar way. The scaling exponents are 1.10 0.08, 0.56 0.08, 15 g g which is again consistent with particle conservation, a b = d , seeequation (9). The crucial difference to what we found in the previous section is that, here, the numerically determined exponents are consistent with the analytical prediction, equation (9), for a vanishing anomalous exponent, indicating a (near) Gaussian nonthermal fixed point: The index 'g' refers to near-Gaussian, in contrast to the strongly anomalous scaling ('a') found in the previous section. We emphasize that our findings are compatible with a small but non-zero anomalous scaling [83]. When we refer to the fixed point as the 'Gaussian' one in the following, it is meant in relation to the much larger value (14) at the anomalous fixed point, not in an absolute sense. This nomenclature includes our anticipation that a small non-zero h g renders the fixed point to be of the Wilson-Fisher type.

Non-thermal fixed points and phase-ordering kinetics
As pointed out in the previous section, if the scaling exponent β is positive, the self-similar build-up (6) of momentum occupations in the infrared reflects a coarsening process. In the following, we discuss in more detail the scaling near the observed non-thermal fixed points in the context of the classical theory of phase-ordering kinetics [7-10, 84, 85] (see [11] for a review). For this, we compare the spatio-temporal scaling relation (6) at the non-thermal fixed point with the scaling forms obtained for the coarsening of ensembles of defects.
Although the evolution of the Bose field, in our settings, is not described by a diffusion equation, as usually is the case in phase-ordering kinetics, we do not expect the effective coarsening dynamics to be fundamentally different in character. As far as scaling laws are determined through the geometric and topological properties of the order-parameter field, the theory of phase-ordering kinetics is expected to apply in principle. In fact, the dynamics at the non-thermal fixed point is found to be in the classical-wave limit of the underlying quantum dynamics. We point that, nevertheless, the scaling exponents resulting for the coarsening dynamics of the unitarily evolving quantum system can be outside the commonly considered universality classes. Not at last, non-thermal fixed points are expected to be the underlying principle for a larger class of universal dynamics phenomena, beyond the closer realm of phase-ordering kinetics.

Porod tails
The theory of phase-ordering kinetics builds on the definition of an order-parameter field y ( ) k t , which, during coarsening, bears nonlinear excitations such as solitary waves, domain walls, and defects. A universal scaling form for the structure factor , , of the defect-bearing order-parameter field can be obtained by means of power counting.
One generically considers an order parameter with  real-valued components, described by an -symmetric model Hamiltonian. For such a field carrying a defect ensemble with a mean defect distance ℓ d , Porod's law [77] (see [8,11] for a generalization to arbitrary  ) predicts a scaling form for the angle-averaged structure factor, in Considering a momentum range between the inverse defect distance and an ultra-violet cut-off set by the size of the defect core, dimensional but otherwise structure-less defect in the order-parameter field. For example, vortices in a 2D Bose condensate are zero-dimensional defects in the order-parameter field f, i.e., we have The scaling (12) of the corresponding bosonic occupation spectrum, for h = 0, is equivalent to the scaling form (17). Hence, the prediction z h = +d 2 for the scaling at a non-thermal fixed point [34][35][36], seeequation (12), is consistent with the Porod law (17) with  = 2, h = 0. Another example are magnetic domains in an order-parameter field obeying a O(1) (Z 2 ) symmetric Hamiltonian such as for the Ising model. The corresponding Porod exponent z = + d 1 also applies to domain walls such as solitons in a ddimensional Bose gas and characterizes the spatial scaling of sound-wave turbulence, see, e.g. [46,86,87]. The rescaled momentum distributions at the non-thermal fixed points studied in the previous section can be fitted, in the infrared momentum region shown as black lines in figures 4 and 6. Note that A is a non-universal parameter which can be fixed by requiring a certain normalisation for the universal scaling function. Thus, the collapsed curve follows a power law The scaling exponent (20) at the Gaussian fixed point (h  0 g ) is considerably smaller than z a , and is consistent with equation (17) for  = 2. The Porod tail with z = 4 g results for dilute ensembles of randomly distributed vortices in d=2 dimensions, seealso [46,47].
For both non-thermal fixed points, the exponent ζ can be understood in terms of an extended Porod law. To this end, we argue that, in general, the Porod law receives an additional correction from the anomalous dimension η, giving The anomalous Porod exponent (21) is consistent with domains in the order-parameter field of an O(1) (Z 2 ) symmetric model, i.e.,  = 1, taking into account the independently determined anomalous exponent h a . Inserting h a from equation (14), we obtain consistent with our numerical result (21). From our simulations, we infer that these domains are distinguished by the two possible orientations of the circulation (vorticity) a cluster of like-sign vortices can take. The conjecture that the relevant order-parameter field is given by the vorticity density w~´v is further supported by the spectral decomposition of the momentum distribution shown in appendix C.3. In figure C3, the red triangles represent the contribution of the divergence-free part of the velocity field to the occupation-number spectrum n(k) shown as black circles. This part corresponds to the vorticity-bearing rotational flow caused mainly by the vortex defects, and its proximity to n(k) demonstrates that it dominates the total spectrum in the region of the Porod-law fall-off, near both, the anomalous and the Gaussian fixed points [46]. This demonstrates that the vorticity plays an important role in the observed dynamics.
We note that the angle-averaged vorticity spectrum scales as ò relative to the singleparticle momentum distribution, seeequation (37) of [46]. Hence, supposing that vorticity represents the order parameter the coarsening of which can be described in terms of a Z 2 -symmetric Landau-Ginzburg type model, one would expect, from our numerical result for z a , equation (21), the respective structure factor to exhibit a Porod tail ∼k −2 . Assuming furthermore the validity of equation (22) for this tail, within a so far not further specified description of phase-ordering kinetics, in d=2 dimensions and for  = 1, corresponding to the Z 2 symmetry, one obtains h  1. A further examination of this is beyond the scope of the present work and will be done elsewhere.
We add the remark that a Z 2 -symmetry breaking clustering transition has been described within an equilibrium formulation in [88,89], with vorticity density as order parameter, see also [90] and our discussion in section 4.6 below.
In this context, it is also interesting that a scaling of n with z = 6 could be associated with the scaling of the angle-averaged radial kinetic energy distribution~~-( ) ( ) E k k n k k 3 3 corresponding to a direct enstrophy cascade [47,[64][65][66]90] known to characterize classical turbulence in a continuously driven 2D incompressible fluid [67]. This power law is steeper than that of a classical Kolmogorov-5/3 law~-( ) E k k 5 3 , corresponding to z = 4.66.
Note, furthermore, that, as seen in figure 3, the power law (21) is found already at very early times, when the vortices still form strong clusters. We recall that an elementary vortex, on scales smaller than the size of its core, gives rise to a power-law momentum spectrum~-( ) n k k 6 , seefigure 7 of [46]. This scaling is here found to survive after the break-up of a non-elementary vortex into elementary vortices, over scales on the order of the cluster size.
We finally note that the η term in equation (22) is also corroborated by the single-particle spectra of soundwave turbulence computed numerically in [46,86]. In this case,  = 1 should also apply because the turbulent transport of compressible (sound-wave) excitations is expected to involve solitary waves which, as mentioned above, represent defects of an O(1)-symmetry breaking order-parameter field.
As was furthermore demonstrated in [41], numerically determined spectra of compressible excitations for = d 1, 2, 3 are consistent with scaling exponents ζ inferred from renormalization-group analyses of the Kardar-Parisi-Zhang (KPZ) equation. Both, the numerically obtained exponents ζ for the compressible excitations, and the exponents determined in the KPZ framework, can be explained by the same non-zero anomalous dimension η. This provides a further support of the anomalous correction in equation (22).

Coarsening as a form of phase-ordering kinetics
Arguments based on classical dissipative dynamics 7 of the order-parameter field [8], or on a renormalisation group approach [11,84], allow to predict the temporal scaling of the defect length as The coarsening dynamics of the order-parameter field is typically described by a diffusion-type equation [11], is the free energy of the order-parameter field ψ defining the model under consideration. The term involving the diffusion constant Γ applies in the case of non-conserved fields, while the first term on the lefthand side accounts for additional conservation laws the field obeys, parametrized by a transport coefficient a [11]. The parameter μ defines the nature of the conservation law which effectively modifies the scaling of the kinetic (typically Laplacian) operator in the equation of motion.
For example, coarsening dynamics near equilibrium, according to the purely dissipative model A within the Hohenberg-Halperin classification [6], has m = 0, while m = 2 applies to the diffusive model B for conserved order parameters.
The coarsening process described by equation (24) is governed by the scaling exponent, see [11] and in particular figure 24 therein for other cases with   m + 2. Inserting the scaling (23) into the spatial scaling form (17), neglecting potential logarithmic corrections, and including a simple infrared cutoff in terms of a constant C results in a scaling form for the structure factor, Assuming that the scaling form (28) and the scaling relation (6) describe the same physical processes, one obtains the scaling relations b b = d and a b = d d between the scaling exponents 8 .
This suggests that the defect coarsening, up to logarithmic corrections, is in one-to-one correspondence with the self-similar transport process near the non-thermal fixed point. For b > 0 d , the transport is directed towards infrared momenta and is subject to particle conservation in the scaling regime (see equations (9) and (10)).

Coarsening at the non-thermal fixed points
With the above results from the theory of phase-ordering kinetics at hand, we can analyze the dynamics near the previously described non-thermal fixed points by evaluating the time evolution of the average defect distance ℓ ( ) t d . Figure 7 shows ℓ ( ) t d for the different initial vortex configurations chosen in our numerics. Blue triangles mark the evolution from the random initial distribution of = ( ) N t 2400 d i elementary vortices and anti-vortices, for which the occupation spectrum is shown in figure 6 and a scaling evolution with exponents (15) is found at late times, indicating the approach of the (near-)Gaussian fixed point. 7 The calculations assume  ( ) O models within the Hohenberg-Halperin classification [6].
The green squares and red circles show the evolution from a vortex-lattice initial state, containing defects with winding number =  w 6. The green data starts from a lattice of 16×16 vortices, the red data from an 8×8-lattice. The evolution of the occupation spectrum corresponding to the red data is shown in figure 4 and gives an anomalously slow scaling evolution, indicating the approach to the anomalous non-thermal fixed point (13). Figure 7 demonstrates that, for both types of initial conditions chosen here, the evolution of the length scales is consistent with scaling behavior,~b ℓ ( ) t t d d , within certain intervals of time. For example, for both, the blue and the red data, scaling is seen for x >  t t 10 h 0 4 2 , simultaneously with the scaling evolution of the occupation spectrum.
The evolution marked by the green squares in figure 7 involves an initial decay of the non-elementary vortices into 1536 elementary defects. At later times, the growth of ℓ ( ) t d (green squares) is also found to approach the Gaussian fixed point. This shows that the system can switch dynamically between the two different types of phaseordering kinetics. As each type is associated with a different non-thermal fixed point, this demonstrates that the system first approaches the anomalous fixed point and only later is attracted by the Gaussian one.
We note that whether and how long the anomalous fixed point is approached appears to be mainly determined by the initial density of defects rather than by their spatial arrangement in the system. Starting with winding numbers larger than one seems to facilitate the approach to the anomalous fixed point. A detailed study of the different conditions determining the attractive basin of the anomalous fixed point is beyond the scope of the present work. As we will show, however, in the next subsection, the effective coupling between the defects and the fluctuating bulk controls the transition from the strongly anomalous to the (near-)Gaussian scaling behavior.
Near both non-thermal fixed points the coarsening exponents are found to be consistent with each other. In the theory of phase ordering kinetics, m > 0 signals that the time-evolving order parameter fulfills additional conservation laws [11]. Formally, in the language of renormalization-group theory, μ takes the role of the scaling exponent η of a wave-function renormalization factor. This corroborates our conjecture, equation (30). Again, we have neglected any possible logarithmic corrections.
We remark that, however, the prediction (9) for β results from a next-to-leading-order large- expansion [38]. Hence, for large  where it is expected to be valid rigorously, the relation (25) applies, and thus, for   2, 3 1 This argument also provides a possible explanation for the independence of the scaling exponent β of  [38].
The energy of the order parameter field per unit defect-core volume is, for   2, mainly given by the contribution from the field around the defect core the scaling of which is fixed by geometry [11]. For example, the velocity field around point defects in d=2 dimensions scales as r 1 in the distance r from the core. Hence, if coarsening is dominated by (radially symmetric) defects, this contribution fixes b = 1 2, for   2, up to the anomalous correction η.
We furthermore note that our results are also compatible with the scaling observed for the scaling dynamics of relativistic  ( ) O models, for  = 2 and 4, giving b  0.5 [38]. As demonstrated in [91] for the case  = 2, charge domains separated by vortex sheets characterize the evolution near the non-thermal fixed point, consistent with the conditions for a coarsening exponent b h = -( ) 1 2 , with h  0. We finally remark that scaling which reflects the occurrence of more intricate topological defects beyond vortices [92], possible in  ( ) O models for  > 2, is expected to require a higher spatial dimensionality,  d 3. In summary, the universal scaling corresponding to the coarsening at the Gaussian non-thermal fixed point is consistent with the class defined by the Hohenberg-Halperin model A. We emphasize that, however, the dynamics at the far-from-equilibrium non-thermal fixed point is not necessarily equivalent to near-equilibrium coarsening at the fixed point describing the corresponding equilibrium phase transition [11]. The coarsening related to the anomalous non-thermal fixed point does not fall, to our knowledge, in any of the known standard classes of dynamical critical phenomena.

Coarsening under coupling to a thermal bath
Our results discussed above suggest that the coarsening dynamics near the Gaussian non-thermal fixed point in the isolated 2D Bose gas is consistent with coarsening. This coarsening occurs according to the Hohenberg-Halperin model A, i.e., within a classical dissipative setting, following quenches across the corresponding thermal phase transition. To investigate this relation further, we present, in the following, analogous results for coarsening in a dilute 2D Bose gas coupled to a thermal bath.
As before, we choose the same type of non-equilibrium initial conditions by setting a weakly irregular rectangular lattice of non-elementary vortices. At time t i , we couple the system to a thermal bath and compute the ensuing time evolution according to the driven-dissipative stochastic Gross-Pitaevskii equation [93]. The temperature is chosen such that, at late times, the system approaches a thermal state without free vortex defects below the Berezinskii-Kosterlitz-Thouless (BKT) transition. Details of the numerical method and results showing the self-similar coarsening evolution of the occupation spectrum are given in appendix C. Figure 8 shows the scaling found in the time evolution of the mean defect distance ℓ d , in analogy to figure 7 for the isolated case. As before, we start from an irregular 8×8 square lattice of non-elementary (anti-)vortices with winding number = | | w 6. We vary the parameters δ and γ, describing the coupling to the bath, seeequation (C.5), in a way that the temperature is held fixed at x = -T 100 h 2 . Figure 8 shows the ensuing time evolution of ℓ ( ) t d for a range of different δ as indicated. For stronger coupling, Gaussian scaling evolution with ℓ ( ) t t d 1 2 is seen during the late-time evolution. As before, in addition to a pure scaling behavior, our numerical results are in principle consistent with weak logarithmic corrections [11,12].
Reducing, however, the coupling δ to the bath, the initial vortex clustering is found to survive for a sufficiently long time such that we are again observing anomalous non-thermal-fixed-point scaling at early times of the evolution, before the system switches to Gaussian scaling. The cross-over to thermal scaling occurs the earlier, the higher the coupling is, at a time t c which approximately scales as d - t c 1.2 0.2 . Hence, as before, our results indicate that the chosen initial vortex ensemble induces the system to be driven closely to the anomalous fixed point at early times of the evolution. We furthermore infer that the late-time Gaussian scaling observed in the isolated system can be associated with a thermal bath built up self-consistently in the system. This implies that the macroscopic parameter controlling the type of fixed point behavior in the isolated system is the effective dissipation to the fluctuating bulk, which builds up depending on the energy in the UV induced through the initial vortex distribution. This leads us to the conjecture that the Gaussian nonthermal fixed point in the isolated system is consistent with the behavior of the coarsening dynamics predicted within the theory of phase-ordering kinetics.

Three-body-collision induced vortex loss
The anomalously slow coarsening behavior can be understood, in leading asymptotic approximation, in terms of the scattering kinetics of vortex defects created in the system. The essential characteristics of this kinetics is seen rather clearly in a motional visualization of the dynamics 9 . As discussed above, the initial states containing non-elementary vortices are seen to exhibit a particularly strong clustering of like-sign defects during the early stage. This observation is corroborated by the results shown in figure 9 for the degree ( ) P t c of clustering at time t. This measure is defined as ) as for the isolated system. Moreover, the crossover time scale from anomalous to Gaussian scaling is controlled by the strength of the coupling δ to the thermal bath. This suggests that also through the different initial conditions chosen in figure 7, a thermal bath is created which eventually forces the system to deviate from the anomalous fixed point. At late times, the data is limited by averaging statistics while at least one vortex-anti-vortex pair is left. in terms of Ripley's  function, where V is the system's volume, is the number of like-sign defects at time t. The θ function counts all pairs of like-sign defects i j , with Euclidean distance d ij below a maximum of ℓ. We chose the integration interval to be bound by x = a 6 h and x = b 150 h . Thus, positive values of P c indicate clustering while for < P 0 c strong anti-clustering is present, i.e., pairing of vortices and anti-vortices. Comparing with figure 7, clustering of like-sign defects appears to be correlated with the anomalously slow coarsening.
Examining the motional visualization (see footnote 8) in more detail one observes that, during the anomalous coarsening, the defects are rather evenly distributed across the volume. While pairs and clusters of vortices with same-sign circulation rotate around each other, pairs of loosely bound vortices and anti-vortices (vortex 'dipoles') travel coherently along the boundaries of the clusters, dragged along by the strong flux. Thereby, loss of vortices appears to be caused by three-body collisions of defects with pairs of both, even-sign and opposite-sign circulation. These collisions allow vortices with opposite circulation to approach each other more closely than a few times the healing length x h . As a result, these vortices can form closely bound pairs moving at a high speed, which are less tied to the large-scale coherent flow. Soon after their formation, these pairs can mutually annihilate by interaction with the background noise.
The motion of point vortices can be described by the Onsager model [79]. This model accounts for an attractive interaction between vortices and anti-vortices which scales with the logarithm of the defect distance. However, the Hamiltonian does not possess a kinetic term which would allow this potential energy to be transferred to kinetic energy and thus enable a direct mutual approach of the defects.
Only at short distances, below a certain threshold of a few healing lengths x h , additional dissipative forces act on the vortices due to the background noise. These forces are not accounted for by the basic Onsager model and allow pairs of vortices to do the final move, viz.directly approach each other and mutually annihilate. During this process, they are still moving at high speed, as a Helmholtz pair, along the direction perpendicular to the dipole-vector.
At larger distances, the only way for a vortex to get close to an opposite-sign vortex and for a dipole to further reduce its binding length, is to scatter with a another vortex or dipole whereby energy can be transferred. As the Onsager potential energy scales logarithmically with distance, the defects taking part in the collision should have similar distances. For recent theoretical studies of such processes, see [47,68,69].
A simple kinetic argument supports our conjecture that, near the anomalous fixed point, three-body scattering prevails while two-body collisions dominate near the Gaussian fixed point: suppose a kinetic equation The two-body collision rate of vortices with mean free traveling time τ can be estimated to scale as t G~-2 1 . That is, within the time where , 34 mfp mfp d on average, a given vortex comes close to another one. Here, l mfp is the mean free path of the vortices, with cross section a few times the healing length and thus scaling as s x h , and v is the mean velocity of the vortex. Consider now the probability that, at the moment of the collision, a further vortex is closer than a small multiple of the healing length to one of the scattering partners. Using our above results, we estimate this probability to scale as x . This is of the order of the fraction of the mean volume per vortex, V N d , which is occupied by a closely bound pair with a linear extent of a few times x h .
As a result, only a fraction x N V h 2 d of all vortex-vortex encounters can be counted as three-body collisions. The resulting rate of vortex loss through three-body scattering is therefore estimated to scale as the product of G 2 and this fraction. Hence the rate is scaling as Inserting τ from equation (34) one obtains the rate to scale as is the standard scaling of a three-body scattering rate between scatterers with concentration scaling as N d .
The crucial difference here is that the velocity v varies with N d as well. The rather uniform distribution of vortices and anti-vortices within the system means that their mean velocity scales with the inverse of the mean pair distance,ṽ l 1 D , which is on the order of the mean defect distance,~l where the constant has units oft 1 , having the scaling solutioñ This is consistent with the increase~l ( ) [ 1 5 of the defect spacing seen in figures 7 and 8. The above kinetics is different from that near the Gaussian fixed point. This can be inferred from the motional visualization (see footnote 8) which shows that vortices are predominantly already paired with each other, with different pairing distance. Clustering of like-sign vortices and concomitant large-scale flows do not appear. As a consequence, dipoles do not have to be formed first, and mutual annihilation follows quickly after two-body scattering of a fast-moving vortex pair with another defect. This is described by a kinetic equation where the rate scales with the mean free traveling time of pairs, , with cross section s x D h , and velocity x- . This is seen in our data, see also previous work in [47].
We conclude that the anomalously slow coarsening process is dominated by the direct interactions between three defects, not obstructed by dissipative interactions with the bath of thermal or near-thermal small-scale fluctuations [29]. In contrast, the main decay mechanism of vortices in the thermally coupled gas is provided by vortex-anti-vortex annihilation after a comparatively fast diffusive mutual approach of the defect pair. This is due to the fact that the drift under the Magnus force, in the dissipative system, exceeds the parallel propagation of pairs of opposite-sign defects according to the Helmholtz vortex law.
The dynamics of the vortex distributions at the anomalous fixed point is likely to have strong relations to 2D classical turbulence. Kinetic-energy spectra and possible relations to classical turbulence were studied in [56, 65, 88-90, 94, 95] exhibiting the relevance of the formation of like-sign vortex clusters for the energy spectra to show Kraichnan-type scaling behavior. In contrast to this, starting from different initial conditions, the defects tend to arrange in randomly distributed bound vortex-anti-vortex pairs showing, at early times, two-body scattering loss, before being slowed down to a three-body-dominated decay later [47].

Relation to scaling behavior in driven stationary systems
Our findings motivate us to ask whether the vortex configurations can be stabilized in a way such that the dynamical scaling behavior at the anomalous non-thermal fixed point can be associated with an appropriate equilibrium fixed point. In the remainder of this chapter we briefly discuss possible interpretations in the context of near-equilibrium dynamical critical phenomena, in particular of vortex glasses in type-II superconductors and a vortex-clustering phase transition in a dilute 2D Bose gas.
The renormalisation-group approach to defect coarsening [11,84] allows for a more intuitive interpretation of the anomalously slow phase-ordering process. Generally, in making a scaling hypothesis for a correlation function, one needs to introduce a dynamical critical exponent z d , defining the relation between the temporal and spatial scaling of the renormalization flow of the order-parameter field 10 . For the structure factor, the scaling hypothesis reads when a regular Abrikosov lattice changes into a phase of disordered, pinned vortices, below the transition to a vortex liquid of more freely propagating defects. The vortices affect the current-voltage characteristics of the superconductor in a way which depends on the density, position, and interactions of the defects. Numerical studies [98] corroborate the measured value of z, and theoretical models for a vortex-glass phase in superconductors have been proposed [99,100]. For a complex Ginzburg-Landau model with a static magnetic field and a Gaussian stochastic mass, van Hove scaling results, , corresponding to a conserved order parameter and an overdamped zero-frequency hydrodynamic mode. A first-order epsilon expansion around the respective upper critical dimension, , similarly as for Ising spin glasses [101]. In d=2 dimensions, this extrapolates to an order-of-magnitude estimate for the dynamical exponent of  z 5.3. In a charge-neutral dilute Bose gas, an Abrikosov lattice can be created by means of rotating-phase laser fields [102,103] and could, possibly, open new perspectives to study dynamical critical behavior in such systems.
For charge-neutral atoms in a closed trap without external fields, a potentially related clustering transition has been reported [88,89]. Onsager's point-vortex model [79] implies a maximum entropy state in which an equal number of elementary vortices and anti-vortices are randomly distributed in the system. It was shown that, continuously driven, such a system can exhibit a dynamical transition from a state described by a positive temperature to one with negative temperature, depending on whether the entropy rises or decreases with increasing energy density. While, for positive temperatures, vortices and anti-vortices mingle and form bound pairs, like-sign vortices are expected to spontaneously cluster with each other at negative temperatures.
Our results indicate that, during the coarsening process, the scaling exponent β switches from the anomalous to the Gaussian value, when the clustering of like-sign defects disappears, see figures 7 and 9. Hence, it is imaginable that in a driven situation, with a vortex ensemble containing a constant number of defects on average and representing a state close to the clustering transition of [88,89], the respective critical dynamics is related to the anomalous fixed point we report here.

Conclusions
In this work, we have presented numerical evidence for a strongly anomalous non-thermal fixed point in a 2D Bose gas. Starting from different kinds of ensembles of elementary and non-elementary vortices, i.e., from a farfrom-equilibrium initial state, the evolving isolated system shows a slow decay of the vortices through mutual annihilation. The decay can be described as coarsening dynamics, showing an anomalously slow self-similar evolution of the single-particle momentum spectrum, equations (6) and (18) (15) and (20), as well as figures 5 and 6. We demonstrate that the self-similar coarsening evolution of the isolated system near the Gaussian fixed point is consistent with the dissipative coarsening when the system is coupled to a thermal bath. This is attributed to the fact that also the isolated system builds up thermal-like fluctuations at short wave lengths. It is a matter of the strength of this coupling to the bath when the Gaussian fixed-point scaling starts to supersede the anomalous one. As a consequence, the observed coarsening represents a type of dynamical critical behavior which appears to be closely related to the equilibrium phase transition at non-vanishing temperatures. In contrast to the entirely thermal state, however, it shows a steep Porod law, equation (20), associated with randomly scattered elementary vortices and anti-vortices on a phase-coherent background [46,48] which are bound to decay.
On the other hand, provided the system is prepared in a suitable initial state, allowing for clustering of likesign vortices, it exhibits the strongly anomalous scaling evolution at early times. The associated anomalous scaling exponent h - 3, seeequation (14), is found to take the same role as the exponent μ characterizing the conservation law governing the coarsening dynamics of the order parameter in phase-ordering kinetics.
Our results indicate that, far from equilibrium, the model possesses a non-thermal fixed point seemingly unrelated to the thermal fixed point. The respective anomalous scaling is seen to last as long as the coupling to thermal noise is sufficiently weak such that three-body collisions of defects represent the dominant loss process. The clustering is reflected in a steeper Porod law at low wave numbers, equation (21), which, as for the Gaussian fixed point, disappears in the thermal state. In both cases, the critical behavior is possible only through strong driving out of equilibrium as accomplished by the initial quench which leads to the coherent hydrodynamic propagation of vortices on the background of a strongly phase-coherent gas.
Here, m is the boson mass, g quantifies the interaction strength in d=2 dimensions, and, in our units,  = 1. Our computations are performed in a computational box of size L 2 on a grid with side length = L N a s s , lattice spacing a s , and periodic boundary conditions. We implement a split-step integration scheme on NVIDIA GTX Titan Black GPUs, in a cluster with 4 host nodes and 4 graphic cards per node. For this, the core functions of the propagation scheme, in particular the fast Fourier transform, are implemented using NVIDIA's CUDA extension for C [104].
For the numerical data presented in this work, we solve equation (A.1) on a grid with N s = 1024 interpolation points and lattice spacing a s = 1 in each direction. Therefore, the only length scale which remains to be fixed is the We use this setup to compute the dynamics of the occupation-number spectrum in analogy to the evolution of the isolated system, seesection 3. Figure C2 shows the approach to the thermal state from a vortex-lattice initial state with initially 8×8 vortices of winding number =  w 6. For this parameter setting, the time evolution of the mean vortex distance, ℓ ( ) t d , after coupling the system to the thermal bath, is shown in figure 8 (see red data points).
The coarsening evolution of the defect scale shows two distinct scaling regimes, for which we analyze the time evolution of the occupation spectrum separately. We find that the occupation spectrum can be rescaled to a universal curve at early times, in the time interval In the upper panel of figure C2, we show the rescaled spectra at five exemplary times within the above interval, together with the non-rescaled spectra in the inset. These exponents are consistent with those near the Gaussian non-thermal fixed point approached in the isolated system, seeequation (15). The infrared momentum powerlaw exponent z =  4.0 0.1 is consistent with the Porod law (20) describing the random vortex distribution near the Gaussian fixed point.
Remarkably, we also find the anomalously slow self-similar evolution of the gas coupled to a thermal bath at early times, with an infrared scaling exponent z =  5.7 0.1. Thus, at early times, the set of scaling exponents is consistent with the scaling at the anomalous fixed point in the isolated system.
Note that the thermally coupled system displays the two types of fixed-point evolution for all values of δ analyzed in figure 8. As discussed in the main text, the coupling strength to the thermal bath changes the characteristic time t c at which the scaling properties change from anomalous to Gaussian (see footnote 8).

C.3. Compressible versus imcompressible excitations
In this appendix we recapitulate the decomposition of the superfluid flow pattern in the system near the nonthermal fixed point, into transverse (incompressible) and longitudinal (compressible) contributions. In this way we can show that the IR scaling is dominated by the incompressible part while in the UV the particles mainly belong to the compressible as well as a quantum pressure components.
To define the decomposition we use the polar representation f r , e x pi , of the field in terms of the density r ( ) t x, and a phase angle q ( ) t x, . This allows to express the particle current * * in terms of the velocity field q =  v . We decompose the kinetic-energy spectrum In this appendix, we discuss the numerical fitting procedure employed in determining the scaling exponents α and β which characterize the self-similar time evolution of the occupation spectrum. Furthermore, we comment on the fits of the universal scaling functions, leading to the spatial scaling exponents ζ. The procedure explained here is used for the data sets presented in figures 4 and 6, and C2 and yields the sets of exponents given in equations (13), (15), (C.6), and (C.7).
The basic strategy for obtaining the temporal scaling exponents is as follows: first, we select a suitable time window within the universal regime, by inspecting the time evolution of a characteristic scale such as ℓ d . Within this window, we compute the occupation spectrum ( ) n k t , j at several times t j , where the times t j are chosen with a logarithmically equidistant spacing. Then, the distance between a chosen reference spectrum ( ) n k t , ref and all rescaled spectra within the window is numerically minimized with respect to α and β. As the occupation spectra are obtained by numerically solving the GPE on a discrete grid, the spectra ( ) n k t , j at each time step are available for fixed, discrete radial momenta k only. Therefore, to be able to evaluate the rescaled spectrum j l j l j , the data is interpolated by means of cubic B-splines, as provided by the Python-SciPy library [111].
In the above procedure, we minimize the quantity summing over all discrete momenta k up to a UV cut-off k Λ and over all rescaled spectra within the scaling window. The cut-off k Λ is set by hand, to exclude the UV tail of any of the spectra from the rescaling procedure. Anticipating positive exponents α and β, this can be achieved by determining k Λ based on the spectrum ( ) n k t , max at the latest time t max within the scaling window. Finally, to minimize the distance (D.1) for a set of occupation-number spectra, we employ the Levenberg-Marquardt algorithm, as provided by the Python-SciPy library [111]. This yields the values for a min and b min for a chosen reference spectrum and a set of times t j . The errors for these values are estimated using the numerically constructed Jacobian of χ (as provided by the SciPy implementation of the algorithm) at the minimum and the residual variance c a b ( ) , min min . To further assess the robustness of the fit, we additionally vary the reference time t ref and the sample times t j which are used in the fitting procedure for a specific scaling window. The values for α and β stated in the main text, together with their errors, result from averaging a min and b min from the Levenberg-Marquardt fits over these additional variations.
Once the optimal values for α and β are determined as described above, the occupation spectra interpolated within the scaling window can be collapsed by rescaling and be averaged. This yields a numerical estimate for the scaling function k ( ) f s . To this estimate, we fit the ansatz (19) for the scaling function, again using the Levenberg-Marquardt algorithm on the distance of logarithms. The normalization constant A is fixed by requiring k =