Non-Thermal Fixed Point in a Holographic Superfluid

We study the far-from-equilibrium dynamics of a (2+1)-dimensional superfluid at finite temperature and chemical potential using its holographic description in terms of a gravitational system in 3+1 dimensions. Starting from various initial conditions corresponding to ensembles of vortex defects we numerically evolve the system to long times. At intermediate times the system exhibits Kolmogorov scaling the emergence of which depends on the choice of initial conditions. We further observe a universal late-time regime in which the occupation spectrum and different length scales of the superfluid exhibit scaling behaviour. We study these scaling laws in view of superfluid turbulence and interpret the universal late-time regime as a non-thermal fixed point of the dynamical evolution. In the holographic superfluid the non-thermal fixed point can be understood as a stationary point of the classical equations of motion of the dual gravitational description.


Introduction
Studies of far-from-equilibrium time evolution of quantum many-body systems have intensified considerably during recent years, driven mainly by new technological possibilities. For example, strongly non-linear dynamics has been observed in ultracold atomic gases [1][2][3][4][5][6], or semi-conductor exciton-polariton superfluids [7][8][9][10]. Moreover, high-energy heavy-ion collision experiments have brought up many questions concerning the thermalisation of the quark-gluon plasma, cf. [11] and references cited therein. Interactions between the constituents of these systems can lead to strong correlations, which render the description of the long-time dynamics intricate and give rise to non-trivial many-body states far from equilibrium. In the case of strong interactions, a quantitative description of dynamical evolution is typically plagued by the absence of suitable approximation techniques, or by technical difficulties such as sign problems when evaluating the dynamics by means of numerical methods. Similarly complicated situations can arise even when a weakly interacting system becomes strongly correlated. In such cases, non-linear excitations can dominate the system's dynamics such as solitary waves or topological defects. In quantum many-body systems the massive appearance of such excitations and their interactions can give rise to quantum turbulence phenomena, i. e., to states the statistical properties of which bear resemblance to correlations observed in classical turbulent systems.
Holographic methods (also known as AdS/CFT or gauge/gravity correspondence) have in recent years opened new vistas on strongly correlated quantum systems. This approach relies on the observation that some strongly coupled quantum theories in D dimensions have a dual (that is equivalent) description in terms of a classical, weakly coupled gravitational theory in D + 1 dimensions. Loosely speaking, the quantum theory can be thought of as living on the D-dimensional boundary of the higher-dimensional space. One therefore often speaks of the 'boundary theory' when addressing the D-dimensional quantum system described in the holographic framework. A finite temperature of the quantum system has a holographic description on the gravity side in terms of a black hole with the corresponding Hawking temperature. The original discovery of the holographic duality [12][13][14] was made for the case of a 4-dimensional gauge theory with maximal supersymmetry. It has subsequently been realised that the method is far more general and can be applied to a variety of physical systems. Applications of the holographic duality by now comprise topics as diverse as ultracold quantum gases and the quark-gluon plasma, for a review see for example [15]. A novel and exciting development is the extension of the holographic duality to various condensed-matter quantum systems, for reviews see for instance [16,17]. The holographic duality allows one to study the dynamical real-time evolution of a strongly correlated many-body quantum system in a genuinely nonperturbative framework. Remarkably, the intricate dynamics of the system -including its far-from-equilibrium behaviour -is entirely captured by a classical gravitational system. Evidently, the holographic description therefore offers the potential to address phenomena in the quantum system which are notoriously difficult to access by other methods. It even carries the promise to discover new phenomena that are unique to strong-coupling situations.
A particularly interesting discovery is that there are gravitational systems which have a dual interpretation in terms of superconductors or superfluids in 2 + 1 dimensions [18][19][20]. Here, the dual gravitational description is in terms of an Abelian Higgs model on and coupled to a (3 + 1)-dimensional spacetime of negative cosmological constant, an Antide Sitter spacetime. The breaking of the Abelian U (1) symmetry at low temperature is associated with the condensation of an order-parameter field and can be interpreted as the emergence of superconductivity or superfluidity.
In the present paper, we will consider a holographic superfluid of this kind at finite temperature and with a chemical potential for the U (1) charge. In this system, vortex excitations exist in the superfluid phase without the presence of an external magnetic field. 1 The system we consider is a (2 + 1)-dimensional relativistic superfluid. Studies of various aspects of that particular holographic superfluid include [21][22][23]. A study of its time-evolution as it relaxes from a far-from-equilibrium initial state, corresponding to an ensemble of vortex defects, was performed in [24]. There, the authors numerically solve the gravitational Einstein-Maxwell-scalar system dual to this superfluid. They have identified a certain regime in the evolution in which the superfluid exhibits Kolmogorov scaling. In the present paper, we perform a similar numerical analysis of the same system. As in [24] we treat the holographic superfluid in the so-called probe limit in which the AdS spacetime in the dual gravitational description is kept fixed. Our numerical methods are sufficiently fast to allow us to investigate two new aspects of the far-from-equilibrium evolution of the superfluid. On the one hand, we follow the system's evolution for a very long time. On the other hand, we study various initial conditions, in particular, we choose random distributions as well as lattices of vortex defects of different densities. This makes it possible to clearly identify the time scales at which the system enters a universal regime.
In particular the investigation of the late-time behaviour of the system leads us to a very interesting observation. Following the propagation and annihilation of the quantum vortices in time, we are able to observe a stage of universal critical dynamics arising in the late-time evolution of the superfluid when the quantum turbulent ensemble relaxes towards equilibrium. More specifically, we observe how the system approaches a nonthermal fixed point, i. e., a far-from-equilibrium field configuration exhibiting universal scaling behaviour [25,26]. We demonstrate how the approach to this fixed point can be related to the dynamics of vortex defects in the order parameter. Non-thermal fixed points were identified, in quantum field theory, as stationary solutions of non-perturbative equations of motion for Green functions [25,27,28]. In the context of non-relativistic Bose gases as well as of relativistic scalar and gauge theories it was shown that superfluid turbulence, related to characteristic distributions of vortex defects [29,30] or more general non-linear excitations [31][32][33], can be interpreted in terms of non-thermal fixed points. General universality classes of such non-thermal fixed points are expected to emerge from a renormalisation-group analysis which includes scaling in space and time [34,35]. The present study is the first analysis of the holographic superfluid in view of the approach to a non-thermal fixed point, comparing universal and non-universal stages of the superfluid's evolution. It opens a new and exciting perspective on time evolution as described in a holographic setting, in particular on the mutual implications of such universal dynamics on both the gravity and the boundary-theory sides.
Our paper is organised as follows. In Sect. 2 we summarise the definition and the properties of the Einstein-Maxwell-scalar gravity model dual to the superfluid, as well as the implementation of the resulting equations within our numerical approach. Sect. 3 contains the details about the different initial conditions considered and summarises basic properties of vortices in superfluids. Further in that section, we present our numerical results on the evolution of the vortex characteristics and of the statistical properties of the ensembles, in particular on occupation number spectra of the boundary theory. In Sect. 4 we finally discuss the holographic perspective on the non-thermal fixed points in the superfluid's evolution. A summary and outlook in Sect. 5 complete the main presentation. Some technical considerations are given in two appendices. Appendix A contains details concerning the numerical implementation, while Appendix B deals with the method for extracting scaling exponents from various spectra.

Holographic superfluid
In this section, we set out the holographic framework for describing the dynamics of a superfluid in (2 + 1) dimensions by means of a gravitational dual in a (3 + 1)-dimensional spacetime.

Gravity model dual to the superfluid
The holographic framework for superfluidity was laid down in [18][19][20]. A scalar field is dynamically coupled to an Abelian gauge theory and gravity on a (3 + 1)-dimensional spacetime. We use the action Here, κ is the Newton constant in four dimensions, and the cosmological constant Λ = −3/L 2 AdS is negative. L AdS sets the curvature scale of the spacetime which arises as a solution of the corresponding Einstein equations. R is the Ricci scalar constructed from the metric g M N , and g is the determinant of that metric. 2 The Lagrangian density L matter accounts for the gauge fields, with the associated field-strength tensor F M N = ∇ M A N − ∇ N A M , and for the scalar field Φ. Here, ∇ M denotes the covariant derivative associated with the Levi-Civita connection. The local U (1) gauge symmetry of the Lagrangian is implemented by upgrading ∇ M to the gauge-covariant derivative According to the holographic dictionary the gauge potential A M in the bulk induces a global U (1) symmetry of the dual field theory. The operator dual to A M is the conserved U (1) current j µ which arises from that symmetry, see App. A.1 for details. Finally, the complex scalar field Φ has mass m and charge q which, by suitable rescaling of the fields, has been pulled out of L matter . If the gravity is taken to be dynamic, q thus quantifies the coupling between the gravity and gauge-matter parts of the model. Solving the gravity model (2.1) allows us to obtain information about the dynamical evolution of a superfluid described by the boundary theory. A superfluid is commonly described by a complex scalar field ψ which, in the symmetry-broken phase, assumes a nonvanishing expectation value ⟨ψ⟩ ̸ = 0 reflecting the presence of a Bose-Einstein condensate. In our holographic model, the field operator ψ is dual to the scalar field Φ, as we explain in detail in App. A.1. The solutions of the equations of motion of the model (2.1) are subject to boundary conditions in the holographic direction. These conditions determine the temperature and the chemical potential for the U (1) charge j 0 in the (2+1)-dimensional boundary theory. In particular, the temperature and chemical potential can be chosen such that the boundary theory is in the symmetry-broken phase with a condensate, ⟨ψ⟩ ̸ = 0 [18]. Holography allows us to compute the time evolution of the quantum expectation value ⟨ψ⟩ starting from various far-from-equilibrium states by solving the classical dynamics of (2.1).
In fact, we also have access to the phase angle of the complex field ⟨ψ⟩ the spatial variation of which encodes information about the superfluid flow. We will use this to construct far-from-equilibrium initial states.
In this article, we consider the so-called probe limit of the action (2.1) in which the back-reaction of the fields Φ and A M on the metric is neglected. This is a good approximation for large scalar charge q, as is clear from the rescaled form of the action (2.1), with 1/q 2 entering as a small pre-factor of the gauge-matter Lagrangian, see for example [19,36]. We can thus treat the gravity and matter parts separately from each other.
Ignoring for the moment the gauge-matter part of the model, the Einstein equations are solved by an AdS 4 spacetime with a planar Schwarzschild black hole. The respective metric reads with the horizon function Here, (t, x) = (t, x, y) are the coordinates of the spacetime on which the boundary field theory is defined, and z ≥ 0 is the coordinate of the holographic direction. As already pointed out, it is often useful to think of the dual field theory dynamics to take place at the boundary z = 0. However, it is the entire bulk information that encodes the boundary physics. The black-hole horizon is situated at z = z h . We use Eddington-Finkelstein coordinates, where a light ray falling into the black hole is affinely parameterised by z, with all other coordinates kept constant. Such coordinates are particularly well-suited for the numerical solution of the system's real-time dynamics [24] and are often employed in holography [37]. Working, then, with a fixed background metric, one is left with the equations of motion for the matter part, i. e., the generally covariant Maxwell and Klein-Gordon equations which are coupled to each other through the bulk electromagnetic current. As the background metric, together with the gauge coupling, allows for spontaneous symmetry breaking in the scalar sector, the problem has thus been reduced to solving a classical Abelian Higgs model on the curved background (2.2).
The holographic model captures important aspects of Tisza's two-fluid model [38], for a review see for instance [39]. More specifically, it has been shown that in the hydrodynamic expansion, to an order which includes only non-dissipative terms, the boundary theory reduces to a relativistic version of the two-fluid model [40]. We point out that our treatment of the holographic model captures more than the effective hydrodynamic limit of the boundary theory. In fact, the holographic description is valid at all scales. Nevertheless, it can still be useful to think of the dynamics of the superfluid in terms of two distinct components. In the probe limit, the presence of the static black hole at z = z h translates, by the AdS/CFT dictionary, to a static heat bath of temperature T = 3/(4πz h ) in the boundary theory. Loosely speaking, this can be viewed as the normal component of the fluid. Similarly, the fields Φ and A M holographically represent the superfluid component. The superfluid component is coupled to the normal component and can dissipate energy and momentum to it. Thus, the model naturally incorporates dissipation to a thermal bath. Further details of the interpretation of the probe limit can be found in [24]. We remark that our formalism is manifestly relativistic. The superfluidity described by the boundary field theory appears in the non-relativistic low-energy limit [24].

Implementation
In this article, we consider the time evolution of the superfluid starting from a far-fromequilibrium situation. For this, we need to numerically solve the full equations of motion for the gauge field A M and the scalar field Φ in the bulk, By the AdS/CFT dictionary, we can extract the expectation value ⟨ψ(x, t)⟩, i. e. the superfluid order parameter, from the asymptotic behaviour of the dual scalar field Φ close to the boundary at z = 0, see App. A.1.
We choose the temperature and chemical potential in the boundary field theory such that the system is in the superfluid phase. Our units are fixed by setting z h = 1 in the metric (2.2), such that the temperature is T = 3/(4π). Then, choosing a chemical potential of µ = 8πT puts the system into the superfluid phase.
As we aim at studying universal aspects of the superfluid, we consider different types of initial conditions, containing topological defects, the details of which are discussed in Sect. 3.1. We take the (x, y) directions to be periodic and use a pseudo-spectral basis for the fields. In order to be able to properly do statistics and suppress finite-volume effects we need to choose our numerical grid sufficiently large. Specifically, we study grid sizes of 352 × 352 as well as 504 × 504 points in the (x, y) plane. The data from these different grid sizes are qualitatively consistent. Since observables computed on the larger grids are considerably less noisy, all data presented in the following were produced on 504 × 504 grids. For this choice, the dimensionless product LT of the extent of the (x, y)-domain and the temperature is LT = 34.4. We use a basis of 32 Chebyshev polynomials and 32 grid points in the holographic direction. To be able to properly assess the genuine late-time behaviour of the system we let it evolve to time t f = 4000 in the aforementioned units, or t f T = 955 in units of temperature. We use an explicit time-stepping scheme for the propagation. We point out that the timestep τ used in our numerics is much smaller than one unit of time, τ ≪ 1. A more detailed discussion of our choice of numerical methods and parameters is given in App. A.
For technical reasons, each of our initial conditions depends on some random data. To get more robust results that do not depend on the definite values of these random data in a specific realisation, we average statistical observables over ten runs for each type of initial condition (random distributions and lattices of vortices, see Sect. 3.1).

Holographic non-equilibrium dynamics
In this section, we first discuss the different types of initial conditions that we employ to induce non-equilibrium behaviour of the superfluid. Specifically, our initial conditions are characterised by ensembles of topological vortex defects. Then, we analyse the evolution of the system in terms of the distribution statistics of the defects. Furthermore, we consider the occupation spectra of the superfluid order parameter.

Initial conditions
For the quantum systems we have in mind, quenches, especially across a or in the vicinity of a phase transition, are being studied intensively, both experimentally and theoretically. A quench in the usual sense involves the rapid change of either a thermodynamic parameter, such as temperature, or a Hamiltonian parameter, for example interaction strength. Temperature quenches across a superfluid phase transition were recently studied in the holographic approach [41,42]. Novel techniques developed for ultracold quantum gases and exciton-polariton superfluids allow to rapidly change parameters of the Hamiltonian. Using these, ensembles of vortex defects can be created in the superfluid [1][2][3][4][5][6][7][8][9][10]. The generic consequence of both types of quenches for superfluids in two spatial dimensions is the nucleation of quantised vortices. This behaviour was also observed in simulations of both relativistic [32,33] and non-relativistic [29-31, 43, 44] Bose systems. Here, we directly prepare ensembles of vortices as quench-like initial conditions for the superfluid's time evolution. For simplicity, we refer to our initial conditions as quenches in the following.
The structure of the core of a vortex, as well as collective properties of the ensemble, are specific characteristics of the superfluid system. The local phase structure of the superfluid order parameter around a vortex defect, however, is determined by topology. The presence of a local vortex with quantisation w ∈ Z\{0} requires that the phase angle ϕ = arg(⟨ψ⟩) of the superfluid order parameter has winding number w around the vortex, i. e., dϕ = 2πw, where the integration contour encircles the particular defect. This property can be used to prepare the initial-time order-parameter field to bear a set of vortices (w > 0) and anti-vortices (w < 0). We do this by 'multiplying' a localised vortex into the phase of the thermalised superfluid, ⟨ψ⟩ → ⟨ψ⟩ · e iwφv with an appropriate polar angle φ v . During the initial propagation of the equations of motion, the appropriate density profile |⟨ψ⟩| 2 around the defect builds up in a short time, providing us with an ensemble of vortices on top of a previously equilibrated system. Technical details can be found in App. A.1. Previous studies of vortex solutions in holographic superfluids include [21,22].
In Fig. 1 we show example realisations of different types of initial conditions we prepare, with a random distribution of vortices and anti-vortices (left column) and a regular lattice distribution (right column). The graphs in the first row show the phase angle field ϕ(x). After imprinting the phase winding of the vortices the system reacts by building up vortex cores and, on a longer time scale, by redistributing the defects, as a consequence of interactions in the superfluid. Thus, the short-time outcome of the respective quench is defined by choosing number, quantisation, and spatial distribution of the vortices. The second row in Fig. 1 shows that indeed vortex cores are built up shortly after the winding phases have been imprinted. This procedure corresponds to the well-established technique of imprinting vortices by 'stirring' a Bose-Einstein condensed dilute atomic gas with the help of a laser [45,46]. Experimental methods have been refined for creating vortices in cold atomic gases [5,6,47,48] and exciton-polariton superfluids [8,9].
In this work, we study two classes of initial vortex distributions: class A consists of random distributions of elementary vortices of winding numbers ±1, while class B comprises regular 12 × 12 lattices of non-elementary vortices with absolute winding numbers |w| > 1, alternating in sign from site to site. For examples from each class see Fig. 1. We observe that a non-elementary vortex of absolute winding number w quickly decays into w elementary vortices of the same sign. Therefore, there are strong correlations built into the initial vortex positions for initial conditions in class B, insofar as after the decay of the non-elementary vortices the like-sign singly quantised vortices are clustered. On the other hand, within the limitations of the finite grid, the initial vortex positions in class A are completely uncorrelated. For each class of initial conditions, we vary the number of vortices by randomly distributing 144, 432, or 720 vortices of either sign in class A, and choosing winding numbers ±2, ±6, or ±10 in class B. This choice implies that after the decay of the non-elementary vortices the total number of elementary vortices is the same in the corresponding cases of both classes. In this way, we vary the initial vortex correlations and the mean separation of vortices, considering both as quench parameters. six different initial conditions are summarised in Tab. 1. Note that in all cases the net vorticity is zero. The vortex phases are imprinted at a time t = t i < 0. After starting the simulation, stable vortex cores develop quickly, typically after ∆t = 5 and ∆t = 10 for quenches of class A and B, respectively. We adjust t i such that at t = 0 the vortex cores are fully formed. For class B (vortex lattices), we furthermore perturb the phase of the bulk scalar field at time t = 0 with random noise to induce variations in the decay pattern. This is illustrated in the upper right panel in Fig. 1. For class A (random distribution) it is not necessary to add phase noise due to the randomness of the vortex positions.

General considerations on vortex dynamics in a superfluid
It will be useful to discuss the real-time dynamics of our holographic superfluid in the framework of an effective description that is well known from other superfluid systems. In this effective picture of quantum turbulence the vortices appear as collective excitations of the order-parameter field. In addition to the vortices the system contains sound waves. These also mediate the effective interaction of the vortex defects. The sound waves can usually be treated in a good approximation as linear perturbations of the order-parameter field.
In this subsection we summarise a few basic properties of the effective description known to characterise vortices and their dynamics in two-dimensional non-relativistic superfluids. For reviews in the context of superfluid helium and cold atomic gases see, e. g., [49][50][51]. The time evolution of an undamped dilute non-relativistic superfluid carrying vortex defects is understood to be well described by the non-linear Schrödinger or Gross-Pitaevskii equation (GPE) [52,53], the classical field equation for the order parameter ⟨ψ⟩ with Schrödinger-type free and quartic interaction parts. The potential field v ∼ ∇ϕ derives from the order parameter's phase angle ϕ = arg(⟨ψ⟩) and describes the local velocity of the superfluid. The velocity field is thus curl free. Vorticity is carried rather by the vortex defects at which the order-parameter field vanishes, permitting a finite circulation of the phase around it.
In studies of two-dimensional turbulence, important observables are the kinetic energy spectrum ∝ ⟨[v(k)] 2 ⟩ and the distribution of vorticity ω = ∂ x v y − ∂ y v x which derives from the local velocity field v = (v x , v y ). As was first discussed in the context of classical fluid dynamics [54][55][56], it is convenient to think of vortices as Coulomb-interacting 'charges' in an effective 'electrostatic' picture where vortices of opposite-sign (like-sign) winding number attract (repel) each other. The dynamics of the GPE vortices is to a good approximation of Hamiltonian character, where the position along one spatial dimension forms the canonical momentum of the position along the other. This forces, in effect, oppositely charged vortices to move in parallel at a fixed distance (Helmholtz law) while equal-sign vortices circulate around each other. The statistics of the vortex distribution contains important information about the temporal and spatial characteristics of the system.
Bogoliubov sound waves form the weak linear excitations of the order parameter field, through which vortices can interact with each other. The interactions with a fluctuating, e. g. thermal, background of excitations of the superfluid causes the vortices to show deviations from the Hamiltonian behaviour. For instance, a sufficiently strong dissipative force can suppress the Helmholtz pair propagation and make oppositely charged vortices move towards each other. This is well known for defect solutions of Ginzburg-Landau equations which represent the generalisation of the GPE with complex parameters [57] and are of relevance, beyond superfluids, for many applications including liquid crystals and biological systems.
We emphasise that not all of the typical properties mentioned above will necessarily be seen in the holographic superfluid that we study here. But we find it helpful to compare our findings to those properties in the following.

Dynamical evolution of the vortex ensembles in the holographic superfluid
We now turn back to our holographic superfluid. In the following we discuss the dynamics of the holographic superfluid for t > 0, after imprinting the initial conditions discussed in Sect. 3.1. The vortices start to move around, subject to an effective interaction. When vortices of opposite sign approach each other sufficiently closely they mutually annihilate. There appears to be a strong suppression of Helmholtz pair propagation: The vortex' and the anti-vortex' trajectories only bend slightly in a common direction before annihilating. This is an indication of strong dissipative effects in the holographic superfluid. Annihilations reduce the total vortex number and change length scales of the vortex distribution. This has important implications for the turbulence properties of the system, which we discuss in Sect. 3.4. At the time t f = 4000 where we end our simulation, only very few vortices are left, typically about 4 to 8 for all types of quenches considered. In our simulations, the non-zero temperature of the black hole is related to dissipation such that sound fluctuations are quickly damped out. The spontaneous generation of vortex-anti-vortex pairs could not be observed. We expect that, after all vortices have annihilated, the excess energy which we initially injected into the superfluid component is completely dissipated into the heat bath. As a consequence, the condensate ⟨ψ(x, t)⟩ would relax to a homogeneous, fully ordered state. Within the finite evolution times, our simulations give indications for this behaviour.
Snapshots of the full time evolution of the system for the two classes of initial conditions are shown in Fig. 2, where the upper (lower) panels correspond to a random, i. e., class-A (class-B, lattice) initial distribution. We plot the superfluid density n(x) = |⟨ψ(x)⟩| 2 at   Note that the apparent initial oscillation of ρ for initial condition B III is not physical, but rather due to uncertainty in the vortex finding algorithm at the high initial densities associated with B III. various times. These snapshots are representative for the different stages encountered in the evolution. 3 In both cases, the main features of the dynamic evolution can be understood as due to the motion and annihilation of vortices. Due to the relatively high initial density of vortices the first stage of the evolution is characterised by a small mean vortex-anti-vortex distance and a high annihilation rate. Each annihilation event releases energy in form of sound waves which is then dissipated into the thermal background [24]. For initial configurations of class B, there is an additional stage in which non-elementary vortices decay. The emerging elementary vortices drift apart for a certain amount of time, thereby expanding like-sign clusters. Naturally, these early stages are highly parameter-dependent and therefore non-universal. The following stage is an evolving dilute vortex gas. Although it is still the annihilation of vortex-anti-vortex pairs which brings the system closer to equilibrium, essential aspects of the time evolution of the system at this stage can be understood from the statistics of the vortex distribution as we will discuss below. Eventually, when only one pair remains, the evolution is governed by the interactions with the sound modes.
Let us now turn to the details of the vortex distribution. At each unit timestep, we determine the position of all vortices and anti-vortices, see App. A.2 for details. In Fig. 3 we show the total vortex density as a function of time (averages are taken over ten runs for each type of quench). For both classes of initial conditions, annihilation proceeds in a non-universal manner, up to a simulation time of t ≃ 400, depending on the particular initial configuration. For initial conditions in class B (lattice), the vortex density is approximately constant at early times. This stage was referred to as the 'drift stage' above. It arises because the like-sign vortex clusters need to expand and dissolve before vortices of opposite sign can encounter each other. Naturally, this takes longer if the initial number of vortices is lower. 4 At time t ≃ 400, starting from any of the initial conditions, a scaling regime is entered, where the vortex density ρ decays algebraically in time, ρ ∼ t −1 . This scaling persists until the end of our simulations at t f = 4000. We observe that the power-law is the same for all tested initial conditions. Hence, starting at t ≃ 400, the vortex density has a universal form. We point out that this decay of vortex density does not exhibit any characteristic time scale. Remarkably, the algebraic decay of the vortex density is characterised by a universal pre-factor. This, together with the value of the exponent, can be explained by a diffusive motion of vortices within a thermal background of sound waves.
To investigate this in more detail we define the following quantities. By l = we denote the mean distance of nearest-neighbour defects of equal sign, while l ̸ = is the mean nearestneighbour distance between vortices and anti-vortices, 5 Here, α = (+, −) denotes the sign of the winding number of the vortex indexed by i α = 1, . . . , N , and x iα is its position. The function nn α (i β ) yields the index of the vortex of sign α nearest to the vortex i β .
In Fig. 4 we show l = (left panel) and l ̸ = (right panel) as functions of time (averages are taken over ten runs for each type of quench). As the inverse of the square root of the vortex density gives the mean vortex distance, the time evolution of the above length scales is directly related to the density evolution shown in Fig. 3. The panels in Fig. 4 show that the nearest-neighbour distances evolve as l = (t) = √ 4D = t and l ̸ = (t) = 4D ̸ = t during the universal late-time period, i. e., diffusively. From our data shown in Fig. 4, we estimate D = ≈ 0. 16   boundary theory for a quantitative comparison with conventional models of superfluidity. Finally, note that the universal regime associated with this scaling behaviour of both l = and l ̸ = is entered at the same time t ≃ 400 for all initial conditions.

Holographic turbulence
So far we have identified at late times t 400 a universal stage of the system's evolution with respect to the effective dynamics of vortices. During this stage, the characteristic length scales of the vortex distribution evolve algebraically, as we have shown in Sect. 3.3. Therefore, their rate of change,l/l, decreases and approaches zero asymptotically. Thus, at late times the system can be considered quasi-stationary. As we will discuss in the following, during this stage the system bears signatures of turbulence. We will, in particular, analyse momentum-space distributions as statistical observables and find them to exhibit algebraic behaviour as it is characteristic for the development of turbulent transport.
While vortices are the building blocks of superfluid turbulence, one needs to study correlation functions of the superfluid order-parameter field in order to obtain a full understanding of the microscopic dynamics of the superfluid. Here, we concentrate on the equal-time two-point correlation function, ⟨ψ * (x, t)ψ(y, t)⟩. Since our system is spatially homogeneous in a statistical sense this quantity can be analysed in relative momentum space, i. e. with respect to the relative coordinate r = x − y. One defines what is known as the (radial) occupation number spectrum as 6 n(k, t) = dΩ k 2π ⟨ψ * (k, t)ψ(k, t)⟩ .
In addition, we average this quantity over a number of runs. This way of approximating the two-point correlation function was also employed in [42]. Although the infrared phenomena we are interested in here are expected to be classical in the statistical sense, it would be very interesting, if numerically demanding, to use the holographic dictionary to extract the full quantum two-point function.
Due to the isotropy of the underlying model we can perform the angular integration without losing information, such that we are left with the radial momentum k = |k|. The standard radial kinetic energy spectrum is defined as and is related to the radial occupation number spectrum via k 3 n(k) = E(k). In general n(k) has a well-defined field-theoretic interpretation.
In his seminal work on turbulence ( [58][59][60], see also [61]), Kolmogorov assumed the existence of an inertial range in momentum space, bounded by two characteristic scales, k in and k diss . At the scale k in , energy is injected into the system. It 'cascades' from momentum shell to momentum shell before it is eventually dissipated into heat at the higher scale k diss . Based on the assumption of such a local transport in momentum space, Kolmogorov found that the kinetic energy spectrum of stationary turbulent flow in an incompressible fluid scales as E(k) ∼ k −5/3 within the inertial range. This result holds also in two spatial dimensions where it is associated with an inverse energy cascade [62]. The essential feature of this so-called Kolmogorov-Richardson cascade is that the system is self-similar, i. e. scale-free, within the inertial range.
Here, we take a more general view on turbulence and analyse the occupation number spectra of the superfluid in terms of scaling laws, n(k) ∼ k −ζ . We extract the scaling exponents ζ by fitting power laws to the spectra after averaging them over ten runs for each initial condition. We estimate the uncertainty of the fitted exponents to be 0.1, for details see App. B. In the following, we relate the scaling seen in the spectra to the statistical properties of the vortex distribution that we discussed in Sect. 3.3. In particular, we want to study whether the universality that emerges in the late-time vortex dynamics is reflected in the occupation number spectra.
In Fig. 5, we show the occupation number spectra n(k) for the initial conditions A II and B II (averages are taken over ten runs for each type of quench). To illustrate their time evolution the spectra are shown at t = 0, 200, 600, and 4000. In both classes, A (random distributions) and B (vortex lattices), the spectra for the parameter sets I and III are very similar to the spectra for the sets II of the respective class. For intermediate and late times we observe inertial ranges in k, the corresponding power laws fitted to the spectra at these times are also shown in the figure. Within the uncertainty in our determination of scaling exponents of occupation spectra (cf. App. B), the scaling exponents at intermediate to late times, t 400, agree even quantitatively for all quench types within each class. Thus, we can discuss the generic features of the spectra at the examples shown in Fig. 5.
First, we note that the initial occupation spectra (corresponding to time t = 0 in Fig. 5) are indeed similar to the far-from-equilibrium initial momentum distributions considered in various classical statistical simulations of dilute Bose gases [29-31, 43, 44]. Similar initial momentum distributions have also been studied in the relativistic case [28,63].
During the early to intermediate stage of the system's evolution, for example at t = 200, the scaling exponents differ between the various initial conditions, indicating that this stage is still non-universal. This could be expected from our analysis of the vortex Occupation number n(k) t 400) to the results reported in [24]. There, the evolution of the holographic superfluid was studied starting from an initial condition very similar to our type B II vortex-lattice configuration. Even though the simulation domain used in [24] was smaller than the one employed here, the initial vortex densities differ by only 10 %. The kinetic energy spectrum was defined in [24], based on hydrodynamic arguments [64], as ϵ kin (k, t) = 1 2 |w(k, t)| 2 k dΩ k . Here, w(k, t) is the Fourier transform of the generalised velocity field w(x, t) = ⟨ψ(x)⟩∇ϕ(x) with ϕ = arg(⟨ψ⟩) the phase of the superfluid order parameter. This definition of the kinetic energy spectrum differs from that of our E(k) in Eq. (3.4). In the cited work, the system was propagated up to time t = 600, and Kolmogorov scaling of the kinetic energy spectrum, ϵ kin (k, t) ∼ k −5/3 , was reported for intermediate evolution times in the range 160 < t < 500. We have analysed our data for quench type B II also in terms of the quantity ϵ kin (k, t). We find, within the uncertainty inherent in the fitting procedure, a scaling law ϵ kin (k, t) ∼ k −1.3 during the intermediate stage of the evolution, for example at t = 200, differing from Kolmogorov scaling.
We now continue with the discussion of our results in Fig. 5. During the evolution, the momentum range where n(k, t) obeys a scaling law gradually grows on its lower end, for all initial conditions. Also, the scaling exponents' absolute values decrease. This can be attributed to the increasing diluteness of the vortex gas. The flow field of a single vortex exhibits a scaling n(k) ∼ k −4 for momenta not resolving the vortex core. The same scaling is expected for randomly distributed vortices and anti-vortices at momenta larger than the mean inverse of the vortex-anti-vortex distance [30]. Hence, the inverse average vortex separation sets a lower cutoff to the scaling regime. As the vortex gas becomes more dilute, the average vortex separation increases, see Fig. 4, so that the lower cutoff of the scaling regime decreases. Furthermore, Fig. 5 shows that the scaling exponents gradually approach the value for single-vortex scaling.
Let us next turn to the late-time evolution. At time t = 600, the scaling exponents of the occupation number spectra are still slightly different for quench classes A and B. They are fitted by ζ ≈ 4.3 and ζ ≈ 4.1 for quench types A II and B II, respectively. At time t = 4000 the scaling exponent is fitted by ζ ≈ 4.1 for both quench types A II and B II. We recall that the uncertainty in the estimation of these scaling exponents is about 0.1, see App. B. As discussed above, within each of the classes A and B of quenches, at late times both the qualitative and quantitative scaling behaviour is unchanged for the alternative choices I and III of the total vortex number. In particular, the scaling exponents in the late stages of the evolution are consistent with 4.1 ζ 4.3 for all initial conditions. To demonstrate this, we plot in Fig. 6 the occupation number spectra for all our quenches (averages are taken over ten runs for each type of quench) at t = 1000 during the late-time stage. At this time, for all initial conditions, the estimates for the scaling exponents are in fact in the narrower range 4.1 ζ 4.2. This indicates the emergence of a universal scaling law for n(k, t). Thus, in the long run, the system loses memory of its initial conditions, confirming our findings in Sect. 3.3. Furthermore, the universal scaling law in the occupation number spectra emerges at late times t 600. This is to be compared with our findings in Sect. 3.3. There, we identified t 400 to roughly mark the onset of the late-time universal scaling behaviour with respect to time.
We note that the late-time scaling exponent still appears to deviate slightly from the single-vortex scaling ζ = 4, indicating that effects of vortex interactions are still relevant in this universal regime. Typically, for all initial conditions, we still have 4 to 8 remaining vortex defects at the final time t f = 4000.

Non-thermal fixed point: the holographic perspective
In the following, we concentrate on the universal late-time behaviour of the holographic superfluid. As we have found in the previous section, the system exhibits two types of universal scaling laws, one in space and one in time, at late times t 600. The first type characterises spatial correlations of excitations and is reflected by the scaling n(k) ∼ k −ζ of the single-particle occupation number spectrum with the exponent in the range 4.1 ζ 4.3. This scaling is observed within an infrared momentum regime which can be viewed as an inertial range. The power-law behaviour terminates in the infrared at a momentum scale that is related to the characteristic length scales of the vortex distribution. The characteristic lengths of the vortex distribution, including the mean vortex separation, follow a scaling law in time, l ∼ t 1/2 , see Sect. 3.3. This leads to the observation that the system's time evolution is slowing down algebraically during the universal stage,l/l ∼ t −1 . Such a behaviour is interpreted as 'critical slowing down' within the context of dynamic critical phenomena well-known from the theory of dynamics near an equilibrium phase transition [65]. Similar scaling features have been observed in numerical calculations of dilute Bose gases far from equilibrium on the basis of a statistical evaluation of Gross-Pitaevskii models [30,43,66]. Within that framework, the scaling features have been interpreted in terms of the more general concept of non-thermal fixed points [25][26][27]. This concept associates stationary points in the time evolution of non-thermally scaling correlation functions with fixed points in a renormalisation-group sense [26,34,35]. In the vicinity of those fixed points, correlation functions assume spatial and temporal scaling behaviour. Non-thermal fixed points thus constitute a generalisation of the concept of critical phenomena near thermal equilibrium. A priori, this situation allows for new universality classes as compared to those in the classification of Halperin and Hohenberg [65].
In the following, we argue that the holographic superfluid indeed exhibits the presence of a non-thermal fixed point. We use observables which have been introduced in [43] for the purpose of numerically identifying non-thermal fixed points in two-dimensional Bose gases. The state of the system can be characterised by two length scales, the mean nearestneighbour vortex-anti-vortex distance l ̸ = , see Eq. (3.2), and the coherence length l C of the superfluid. The latter is defined as with the autocorrelation function g (1) , 7 Here, d 2 x/A denotes the average over the simulation domain A, and dΩ r /(2π) the angular average for the difference vector r. l ̸ = is a measure for the state of the system from the point of view of the effective vortex picture. The coherence length l C measures the degree of spatial correlations of microscopic excitations. The time evolution of both lengths is depicted in Fig. 7 (averages are taken over ten runs for each quench type). Both l ̸ = and l C show universal scaling behaviour in the late stage of the evolution. The scaling exponent of l ̸ = is found to be 0.5, as discussed in Sect. 3.3. The scaling exponent of l C can be narrowed down to a value between 0.5 and 0.6. Thus, it is conceivable that the latter scaling exponent agrees with the scaling exponent of the characteristic length scales of the vortex distribution. This would imply a constancy of ratios of different characteristic length scales typical for the approach to a non-thermal fixed point or a critical point, see for example [67]. It is useful to study the time evolution of the system in a reduced configuration space consisting of the two scales l ̸ = and l C . Fig. 8 shows the trajectories of the system for different initial conditions (averages are taken over ten runs for each quench type) in this configuration space. Here we choose to plot the inverse lengths, 1/l ̸ = and 1/l C , on the axes because in the thermodynamic limit one generically expects characteristic length scales to diverge at critical points. While the inverse coherence length decreases monotonically for all initial conditions, the behaviour of 1/l ̸ = is different for the two classes of initial conditions. Starting from random distributions of vortices (class A), 1/l ̸ = decreases monotonically due to vortex-anti-vortex annihilation. For vortex lattices (class B) 1/l ̸ = increases during the initial drift phase as the vortices of opposite sign move closer to each other. Subsequently, 1/l ̸ = decreases monotonically for all of our initial conditions during the universal stage. Eventually, for all of our initial conditions the system is attracted by the point of maximal coherence and maximal vortex-anti-vortex separation, as marked in the figure. Close to the fixed point, all trajectories meet in a universal curve, indicating that the memory of the initial condition is completely lost and, with that, signalling the universal stage of time evolution. Due to the algebraic slowing-down the system spends a major part of its time evolution near the fixed point. For time t → ∞, all curves would bend over towards (1/l ̸ = , 1/l C ) → (∞, 0) because the last vortex-anti-vortex pair annihilates and the coherence length diverges. 8 This would mark the approach of thermal equilibrium. Our findings indicate that the abstract concept of a non-thermal fixed point in the quantum dynamics can be interpreted in a simple way on the gravity side. Specifically, a fixed point corresponds to a time-independent solution of the classical bulk equations of motion (2.4) and (2.5). There is one unique static solution which is completely stable and corresponds to thermal equilibrium, i. e. the thermal fixed point. But, as for many systems of non-linear partial differential equations, there will be a series of stationary points with different stability characteristics. The intriguing properties of non-thermal fixed points are reproduced by partially stable stationary points, i. e. points which have at least one attractive direction in phase space. This guarantees that the time evolution first approaches such a fixed point -regardless of details of the initial condition -before turning towards thermal equilibrium. Here, we have identified such a partially stable stationary point for the dynamics of the holographic model in the reduced configuration space of Fig. 8.

Summary and outlook
In this work, we have studied the far-from-equilibrium dynamics of a holographic superfluid at finite temperature and chemical potential. In the holographic framework, the classical solution of an Einstein-Maxwell-scalar system in (3+1) dimensions is dual to the quantum dynamics of a (2 + 1)-dimensional superfluid. We have performed a numerical study of the time evolution of the holographic superfluid in its superfluid phase starting from a variety of far-from-equilibrium initial conditions corresponding to quenches of the system. In particular, we have imprinted various kinds of large ensembles of topological defects, i. e. quantised vortices and anti-vortices, in the superfluid, and have evolved the resulting non-equilibrium states for a long time. The evolution can be interpreted in terms of an effective description in which vortices and antivortices interact via sound waves. When vortex-anti-vortex pairs annihilate energy is quickly dissipated into the heat bath.
At early times, the evolution of the system strongly depends on the initial conditions, that is on the distribution and density of vortices and anti-vortices. We find that the system exhibits a new non-equilibrium universality regime in the late-time stages of the evolution. Starting from any of our initial conditions, that regime is entered at times t 600 in the units we have chosen. This universal regime is characterised by a dilute 'gas' of vortex defects. The system remains in the universal regime until the final times t f = 4000 we observe in our simulations. In order to obtain more insight into the dynamics of the universal regime, we have analysed the time evolution of characteristic length scales of the ensemble of topological defects, and have studied spatial correlations of microscopic excitations via the occupation number spectrum. We have observed scaling laws in these observables in the universal regime of the system's evolution and have determined the corresponding scaling exponents. The power-law behaviour of the occupation number spectrum can be related to turbulence. During the universal stage, the occupation number spectrum scales as n(k, t) ∼ k −ζ in an infrared-to-intermediate-momenta inertial range with scaling exponent 4.1 ζ 4.3. This value is close to the scaling n(k, t) ∼ k −4 observed in classical statistical simulations of two-dimensional Bose gases [29,30] where it was shown to be related to a dilute random distribution of vortices and anti-vortices and, in turn, to so-called strong wave turbulence [27].
We have made the interesting observation that the evolution of the system exhibits critical slowing-down during the universal regime. This is a natural feature of non-equilibrium dynamics in the vicinity of non-thermal fixed points. We have presented evidence support-ing this interpretation of the observed late-time universal dynamics of our superfluid. We have hence found the first evidence for the presence of a non-thermal fixed point in the dynamics of a holographic system. This is particularly interesting as the occurrence of a non-thermal fixed point indicates that the dynamics of the corresponding non-equilibrium state is universal and independent of the microscopic details of the system. The same universal dynamics is then expected to occur in a variety of other quantum systems, thus connecting very different fields of physics, see for example [68][69][70]. The observation of a non-thermal fixed point not only gives a new view on the dynamics of the holographic superfluid. The gauge/gravity duality also adds a new dimension to the understanding of non-thermal fixed points, as it translates quantum dynamics in the boundary theory to classical dynamics in the bulk. Thus, the duality might provide a new avenue for an analytic treatment of non-thermal fixed points, as these are dual to stationary solutions of nonlinear partial differential but classical equations of motion.
At an intermediate stage (at about 200 t 400) during the evolution of our system we observe a power-law energy spectrum which is consistent with Kolmogorov 5/3-scaling. Such behaviour has been reported before in [24] for the same system that we study here. However, we find Kolmogorov scaling only for part of the various initial conditions we have chosen. Specifically, Kolmogorov scaling emerges when the evolution starts from random distributions of vortex defects, while this is not the case when we choose vortex lattices as initial configurations, the type of initial conditions used in [24]. At the intermediate times concerned here, the system is still affected by the initial distribution of vortex defects, while soon afterwards it enters the universal regime discussed above which exhibits a different scaling behaviour. Kolmogorov scaling thus seems to occur only in a transient way. However, it appears that details of the initial conditions can influence the system for a sufficiently long time to completely prohibit the emergence of a transient Kolmogorov scaling. This may also depend on the way in which statistical noise is implemented in the numerical simulation for particular sets of quenches. Further study is needed to understand better the conditions for the occurrence of Kolmogorov scaling.
There are many open questions to be addressed in order to obtain a full understanding of the holographic superfluid. It will certainly be very interesting to study how the properties and various observables of the system depend on temperature and chemical potential, both within and across the phase boundaries of the superfluid phase. Another interesting point for further investigation concerns the coupling strength of the holographic superfluid. In general, the holographic duality maps a weakly coupled classical gravity system to a strongly coupled quantum system. It would be useful to investigate what exactly that means in relation to other descriptions of superfluids which have the coupling as an explicit parameter. A detailed comparison of the behaviour of, for instance, typical length or time scales in the holographic superfluid and in semi-classical (Gross-Pitaevskii) descriptions of superfluidity can give insight into this problem.
for useful discussions. C. E. and A. S. are grateful to the Mainz Institute for Theoretical Physics (MITP) for its hospitality and its partial support during the final stages of this work. This work was supported by the Alliance Program of the Helmholtz Association (HA216/EMMI), the Deutsche Forschungsgemeinschaft (GA677/7,8), and the Heidelberg Center for Quantum Dynamics. A. S. acknowledges support in the framework of the cooperation contract between the GSI Helmholtzzentrum für Schwerionenforschung and Heidelberg University.

A Equations of motion and their numerical implementation
In this appendix, we present the explicit form of the equations of motion of our (3 + 1)dimensional holographic model and technical details concerning their numerical solution. The model consists of an Abelian Higgs model coupled to gravity on an AdS 4 spacetime with a Schwarzschild black hole.

A.1 Equations of motion and holographic dictionary
As discussed in Sect. 2.1, we solve the equations of motion of the matter part of the action (2.1) which we quote here again for convenience: We consider a fixed AdS 4 -Schwarzschild background as defined in Eq. (2.2). Our general setup follows [24]. The dual field theory can be tuned to be in the symmetry-broken phase by adjusting the horizon temperature and the chemical potential, if the mass m is suitably chosen. Our choice m 2 = −2/L 2 AdS is within the range of permissible values. For convenience we set L AdS ≡ 1.
In deriving the equations of motion we use the probe limit of large charge q, in which the Einstein equations decouple from the equations for the Abelian Higgs model. q then drops out of these equations. Denoting by ∇ M the metric covariant derivative and by D M = ∇ M − iA M the combined metric and gauge covariant derivative, the equations of motion for A M and Φ take the form with the current We fix the gauge freedom by choosing the axial gauge, A z = 0. Our aim is to describe the dynamical evolution of inhomogeneous solutions of the above equations, describing, e. g., vortex excitations of the dual superfluid. For this we begin by deriving static solutions spatially homogeneous in x and y. Taking the fields A M and Φ independent of the coordinates x, y, and t, we obtain the explicit equations of motion, second order in derivatives with respect to the holographic coordinate z. For the gauge field, they read where i = x, y and the prime denotes a derivative with respect to z. The last equation originates from the dynamic equation for A z and remains as a gauge constraint to ensure the chosen axial gauge. For the scalar field we find The limiting value of the electrostatic component A t of the gauge field at the conformal boundary z = 0 sets the chemical potential in the dual gauge theory, Our units are fixed by setting z h = 1 in the metric (2.2), which also fixes the black-hole temperature. Then, in these units the temperature in the boundary theory is T = 3/(4π), and we further choose µ = 6, such that µ/T = 8π. With these choices the system is in the superfluid phase. At the black-hole horizon z = z h , we need A t (z h ) = 0 for regularity of A M . As we do not want to switch on sources for the spatial parts of the U (1) current dual 9 to A M because these would break isotropy, we impose vanishing of A x , A y at the boundary, and A x (z h ) = 0 = A y (z h ) at the horizon. Close to the boundary, the scalar field Φ behaves as Φ(t, x, z) = η(t, x)z + O(z 2 ) , (A. 10) with the source field η. Since we want the scalar operator ψ to form a condensate due to spontaneous symmetry breaking we choose η = 0. Then, the expectation value ⟨ψ(t, x)⟩ of the operator dual to Φ can be identified with the coefficient of the quadratic term in the expansion, To summarise, the boundary conditions for the gauge fields A µ are where µ is the chemical potential. There is one explicit boundary condition for the scalar field Φ, namely η = 0, which can be expressed as The second boundary condition for Φ is a behavioural one: Φ be regular at the horizon. Physically, this represents the infalling boundary condition [71]. We use these boundary conditions both for the construction of the equilibrium solution and for the full equations of motion. The solution of Eqs. (A.5)-(A.8) yields the holographic representation of the thermal equilibrium configuration of the superfluid. Using (A.10) we extract the order-parameter field ⟨ψ⟩ from the solution. With our choice of the parameters µ and T , and of units, the equilibrium value for the density of the superfluid order parameter is n = |⟨ψ⟩| 2 ≈ 166.8.
Next, let us consider the full equations of motion. To construct the initial conditions, we perturb the homogeneous solution by putting vortices on top of it at an initial time t = t i < 0. For a vortex of winding number w, this is done by imprinting its winding structure onto the bulk scalar field Φ, locally Φ(t i , x, z) → Φ(t i , x, z) · e iwφv for every zslice, where φ v is a polar angle in the xy-plane, centred on the respective vortex. At the vortex positions x v the scalar field has to vanish to remain well-defined after the phase winding is imprinted, so we set Φ(t i , x v , z) = 0 along all z. At all other grid points we leave the absolute value of Φ at the equilibrium value. Starting the simulation, the system very quickly builds up stationary density profiles around the vortex cores, approximately within 5 or 10 units of time for quench classes A and B (see Tab. 1), respectively. Therefore, we choose t i = −5 and −10 for quench classes A and B, respectively. To induce variations in the decay pattern of the vortex lattices (class B), we additionally perturb the phase with noise e i Re(ζ(x)) at time t = 0 (and only there) when the vortex cores are fully formed. ζ(x) is obtained as the inverse discrete Fourier transform of ζ(k). This is in turn constructed by populating the Fourier modes in a disk of radius N/100 about the origin in momentum space with O(1) complex Gaußian noise, where N is the number of grid points along each of the directions x, y. All Fourier modes outside this disk are set to zero.
We thus excite the system to a non-equilibrium state, as discussed in Sect. 3.1. Solving the full set of equations of motion, we follow the subsequent evolution of the superfluid order parameter ⟨ψ⟩. At every time step, we extract it from the bulk field Φ using (A.10). It is convenient for numerics to rescale the scalar field Φ by 1/z and work withΦ ≡ Φ/z. Using the 'lightcone derivative' .14) and with ∇ = (∂ x , ∂ y ) we eventually obtain the following system of equations to solve:

A.2 Numerical methods
For the z-parts of the equations of motion, we use a collocation method with a basis of Chebyshev polynomials on a Gauß-Lobatto grid with 32 points in the holographic direction (see e. g. [72]). After setting the horizon radius z h = 1, we switch to a new coordinatẽ z ∈ [−1, 1] defined by z = 1 2 (z + 1) . (A. 19) With respect toz, we work entirely in real space, implementing ∂z-differentiation via matrix multiplication.
We treat the directions x and y as periodic. This enables us to use discrete Fourier transforms to efficiently compute ∂ x -and ∂ y -derivatives. For the (x, y) grid we choose a grid constant of a = 1/3.5 in the aforementioned units. All data shown in this paper were produced on a 504 × 504 × 32 grid (x, y, z-directions).
To compute the equilibrium configuration of the system we have to solve the boundaryvalue problem defined by (A.5)-(A.8). We treat the nonlinearity of the equations by using the Newton-Kantorovich iteration procedure. Technically, we solve the resulting system of linear equations via an LU decomposition with full pivoting, see e. g. [73].
For the full set of equations of motion we solve the boundary value problems in (A.16), (A.17), and (A.18) for ∇ + A x , ∇ + A y , and ∇ + Φ, undo the shifts (A.14) to get the time derivatives, and use a fourth-order/fifth-order Runge-Kutta-Fehlberg algorithm with adaptive timestep size to propagate the fields one timestep forward. We allow the timestep size τ to vary in the range 0.001 ≤ τ ≤ 0.1 in our aforementioned units. Finally, we use (A.15) to update A t .
During the simulations, in order to investigate the spatial characteristics of the vortex ensembles, we determine the positions of all vortices and anti-vortices in the superfluid order parameter ⟨ψ(x, t)⟩ at every unit timestep. To this end, we iterate over the whole (x, y) grid, measuring the integrated phase of the superfluid order parameter around each elementary plaquette.

A.3 Performance
We implement the algorithm in C++, using the fftw3 library [74] for Fourier transformations and the Eigen library [75] for high-performance linear algebra. We parallelise the numerical code for multicore architectures with OpenMP [76].
Propagating a 352 × 352 grid with 32 points in the holographic direction up to time 600 utilising 4 threads on a regular desktop computer with an Intel i7 processor takes approximately 8 hours.
Running on the Intel Xeon server architecture using 16 threads, a run up to time 4000 on a 504 × 504 grid with 32 points in the holographic direction takes around 120 hours.

B Fitting of scaling exponents of occupation spectra
Scaling laws of occupation number spectra take the simple form n(k) ∼ k −ζ with scaling exponent ζ. In the case of the occupation number spectra discussed in Sect. 3.4, we fit the power law exponents. We employ the Levenburg-Marquardt least-squares fitting algorithm to determine these scaling exponents after choosing a momentum range for the fit by eye.
There are two sources of uncertainty in this: the determination of the momentum range and the uncertainty inherent in estimating the best fit parameters. The fitting algorithm reports an uncertainty of typically 0.03 for the scaling exponents. Fixing the momentum range introduces a larger source of uncertainty. To address this, we vary the endpoints of the momentum range and thus determine a range of estimators for the scaling exponent. As a result of that, we estimate the uncertainty in our determination of scaling exponents of the occupation number spectra to be 0.1. Therefore, we state results for fitted scaling exponents with one decimal digit.
The data on the time dependence of our observables is noisier than the data on occupation number spectra such that using a fitting algorithm is not appropriate in this case. The uncertainty in our determination of scaling exponents with respect to time is difficult to assess.