Universal far-from-equilibrium Dynamics of a Holographic Superconductor

Symmetry breaking phase transitions are an example of non-equilibrium processes that require real time treatment, a major challenge in strongly coupled systems without long-lived quasiparticles. Holographic duality provides such an approach by mapping strongly coupled field theories in D dimensions into weakly coupled quantum gravity in D+1 anti-de Sitter spacetime. Here, we use holographic duality to study formation of topological defects -- winding numbers -- in the course of a superconducting transition in a strongly coupled theory in a 1D ring. When the system undergoes the transition on a given quench time, the condensate builds up with a delay that can be deduced using the Kibble-Zurek mechanism from the quench time and the universality class of the theory, as determined from the quasinormal mode spectrum of the dual model. Typical winding numbers deposited in the ring exhibit a universal fractional power law dependence on the quench time, also predicted by the Kibble-Zurek Mechanism.

Non-equilibrium quantum phenomena are of wide importance across several disciplines of physics. Despite their fundamental relevance, few widely applicable principles are known for field theories far from equilibrium [1][2][3]. Gauge-gravity duality is a powerful tool in this respect, as it provides a well-defined first-principles framework to study strongly-coupled field theories in the non-equilibrium setting. The theories we consider are strongly correlated in the sense that there exists no weakly-coupled quasi-particle picture upon which one could base a perturbative treatment. This makes understanding their dynamics, even near equilibrium, an extremely challenging problem. Consequently progress has been confined mostly to cases in which the dynamics is integrable [4][5][6][7].
Holography [8] is a powerful tool enabling us to explicitly analyze such theories without having to rely on integrability: One maps the quantum-field theory setup of interest to a 'dual' gravity problem in asymptotically Anti de Sitter (AdS) space-time, which can be solved in great detail. This duality can reveal significant new insights, for example it provides new examples of interaction driven localization transitions [9], as well as modelling finite-temperature transport near the superfluid/insulator critical point [10]. This article extends the framework in yet another direction, namely into the regime of non equilibrium dynamics across phase transitions. We analyze the breaking of a continuous U (1) symmetry, by studying the evolution of black holes in AdS, giving a dual description of superconductors.
Previous work on dynamics of holographic superfluids analyzed the condensation process by perturbing an uncondensed initial state below criticality [11], as well as the nonequilibrium phase diagram resulting from sudden quenches of the superfluid phase [12]. Both of these studies focus on spatially homogeneous configurations, as do [13,14], which present scaling results for the timescale of the breakdown of adiabaticity during finite-rate quenches.
Subsequent work, relaxing the constraint of spatial homogeneity, explored vortex turbulence of holographic superfluids [15] (see also [16] for recent work on anisotropic quenches). In this study we explore a canonical non-equilibrium paradigm associated with broken symmetries, leading to universal scaling results in the formation of defects when a critical point is crossed at a finite rate [17,18]. In this context, inhomogeneous configurations of the order parameter are known to be crucial during the formation of topological defects. We now describe some general features of the physics of defect formation during second-order phase transitions.

Note:
Independently, Chesler, Liu and Garcia-Garcia [19] have explored the Kibble-Zurek scaling via the dual black hole quasi-normal mode spectrum and numerical analysis in two spatial dimensions.

A. Universal Scaling of Winding Number
Phase transitions were traditionally studied as equilibrium phenomena. Thus, in the broken symmetry phase, the whole system was assumed to make the same selection of the broken symmetry vacuum. However, as noted in the cosmological context, where a sequence of symmetry breaking phase transitions is thought to have resulted in the familiar fundamental forces, rapid cooling of the post Big Bang Universe combined with relativistic causality makes it impossible to coordinate symmetry breaking outside of the Hubble horizon. As a consequence, distinct domains of the Universe will choose broken symmetry vacua on their own. The resulting mosaic of broken symmetry vacua will -in the course of the subsequent phase ordering -attempt to smooth out. As Kibble pointed out [17], these disparate choices can crystalize into topological defects that may have significant consequences for the subsequent evolution of the Universe.
As noted by one of us [18], systems undergoing second order phase transitions cannot ever follow a sequence of instantaneous equilibria. This is because of the critical slowing down in the vicinity of the critical point: Relaxation time τ of the order parameter diverges as a function of the dimensionless distance from e.g., the critical temperature T c : so that the "reflexes" of the system are characterized by the universal power-law, where z and ν are the dynamic and correlation length critical exponents, respectively. An arbitrary cooling ramp can be linearized around the critical point As a result of critical slowing down, the system loses the ability to adjust to the change even when it happens slowly, on any finite quench timescale τ Q . The instantt when the system can no longer keep up with the change of happens when its relaxation time becomes comparable with /˙ , the rate of change of . This leads [18] to the equation: Using the critical slowing down scaling relation, Eq. (2), one obtains: This time scale allows us to split the crossing of the phase transition into a sequence of three stages in which the dynamics is first adiabatic (t < −t ), then effectively frozen during the interval (−t,t), and finally adiabatic again deep in the broken-symmetry side of the transition (t >t ). The value ofˆ corresponding to the time scale separating the frozen and adiabatic stages,ˆ enters into the sonic horizon estimate, the analog of the causal horizon in the early universe.
As a result, the characteristic sizeξ of the domains that can coordinate the choice of broken symmetries exhibits a universal power-law dependence [18] on the quench time: The density of topological defects can then be estimated by recognizing that aξ-sized fragment of defect can be expected within a volumeξ-sized domain. This reasoning is expected to yield correct scaling of the density of defects, but only an order of magnitude estimate of the prefactor. The scenario just described is often referred to as the 'Kibble-Zurek mechanism' and we will refer to it throughout this paper as 'KZM'.
Here we consider the setting where the phase transition happens in a ring of circumference C. In view of our above discussion one can expect that the broken symmetry will be chosen independently in sections of sizeξ, so there will be C/ξ fragments of the ring that select broken symmetry independently [18]. Consequently, phase mismatch resulting from a random walk of phase with C/ξ steps needed to circumnavigate C is expected to lead to: ∆Θ ≈ Ĉ ξ . The net "phase distance" travelled ∆Θ will then have to settle to a multiple of 2π, It follows that the dispersion of the values of W will scale with the quench rate as [20]: This universal power-law is expected to hold whenever C/ξ 1, away from the onset of adiabatic dynamics. The prefactor has to be taken with the proverbial "grain of salt". As seen in previous numerical experiments, such KZM calculations yield correct scalings, but tend to overestimate the density of defects [20][21][22][23] as well as the typical winding numbers [24,25].

II. A SUPERCONDUCTING RING IN HOLOGRAPHY
We study a one-dimensional superconducting ring of circumference C. The quantity of interest is the winding number where φ is the angle along the ring. It follows from gauge invariance that in equilibrium we have (dΘ − A) = 0 so that a non-vanishing winding number is accompanied by a nonvanishing line integral of the vector potential around the loop. It should be noted that W as defined above is gauge invariant under single-valued gauge transformations.
We consider the so-called 'bottom-up' holographic superconductor model [26,27] in threedimensional AdS [28], denoted from here on AdS 3 and work in the probe limit. This amounts to neglecting the effect of the charged components of the system on the neutral ones and vice versa [15,26]. From the gravity point of view this corresponds to neglecting the bulk gravitational backreaction of the charged scalar and the Maxwell field, but keeping the backreaction of the scalar and Maxwell field on each other. This means that we consistently study the dynamics of a Maxwell field A coupled to a scalar field ψ on a fixed gravitational background. We choose the scalar field to have vanishing mass, which means that it is dual to a classically marginal operator in the boundary field theory. We do not expect our results to differ qualitatively for other choices of the mass, as long as it remains in the range corresponding to a marginal or relevant operator (we refer here to the RG scaling dimension in the UV). The details of our bulk action and equations of motion can be found in the Methods section. Maxwell theory without the charged scalar in AdS 3 and its dual field theory have been studied previously in [29][30][31], taking advantage of the fact that one may conveniently dualize the bulk vector field to a bulk scalar.
Studying the field theory at finite temperature corresponds to studying the bulk system in the background of a black hole, in our case the three-dimensional BTZ black hole [32], which is characterized by the parameter z h , its horizon radius. For more details see Fig. 1.
With this data we have the Hawking temperature Via the AdS/CF T correspondence this is directly translated into the temperature of the dual field theory [33]. This theory is therefore studied in a thermal ensemble at temperature T H , which is an external parameter in our dynamics. In order to fully specify the ensemble (in the equilibrium case) as well as the dynamics we must augment the equations of motion with suitable boundary conditions at the UV boundary z = 0. We impose in addition to the finite temperature T a fixed charge density ρ corresponding to Neumann boundary conditions on the bulk gauge field. We give Dirichlet boundary conditions for the scalar field. From this it follows that the field theory is studied in the absence of a source for the order parameter density. This condition ensures that any symmetry breaking occurring will be spontaneous. In order to compare our full non-equilibrium results to the prediction of KZM, we first need to determine the universality class of the field theory just defined, as well as the microscopic parameters τ 0 and ξ 0 . As is usual in holography, this boils down to computing the so-called quasinormal modes of the bulk black hole [34].
III. NEAR-EQUILIBRIUM UNIVERSALITY: ν, z AND τ 0 , ξ 0 In order to determine the critical exponents z, ν, as well as the microscopic parameters τ 0 and ξ 0 , it is sufficient to study the holographic superconductor near equilibrium. We will study the non-compact case, since the finite size of the ring in our simulations does not affect the results for local correlations, so long as ξ is much smaller than the circumference. Since ξ formally diverges near T c , this assumption will break down for extremely slow 'quenches', but we have not found that our simulations entered this regime for the range of τ Q under study.
Adapting the analysis of [35] to the present situation, the correlation function of the order parameter field O near, but slightly above, T c takes the momentum-space form where the last expression holds for small ω, k and we introduced the parameters c and 1/ξ 2 .
We have defined the dynamical susceptibility of the order parameter χ(ω, k) in the first equality. One obtains the correlation function from the susceptibility χ(ω, k) by a Fourier transform We can find the relaxation time by looking at the Fourier transform of the zero-momentum response, χ(ω, k = 0), which at late times takes the form and τ 0 = 0.89. We have the approximate relation ξ < 0 = √ 2ξ > 0 between ξ 0 below and above the transition. This is the "opposite" of the usual Landau-Ginzburg mean-field relation Thus, while the critical exponents have the "mean field" values, equilibrium correlations lengths clearly show that the theory we are dealing with is outside of the the Landau-Ginzburg paradigm.
The equal time correlation function, following from the Fourier transform of the static correlation function χ(ω = 0, k), falls off like From the fact that the relaxation time is proportional to the square of ξ it follows that the dynamical critical exponent z = 2. This leaves us to determine the correlation length critical exponent ν. For this we must study how the correlation length ξ diverges near the critical In holography the poles of two-point functions of boundary operators coincide with the quasinormal modes of the bulk fields dual to the operators appearing in the correlation function [36]. As we saw above, in order to extract ξ and τ it is sufficient to study the static susceptibility χ(ω = 0, k) and the dynamical susceptibility χ(ω, k = 0) separately. As shown in Fig. 3, we find that the spectrum of quasinormal modes of χ(ω, k = 0) contains poles only in the lower half complex plane (as demanded by stability), while χ(ω = 0, k) has two series of poles along the positive and negative imaginary axis. This structure is also apparent from the correlation function (12). In each case the relevant poles are the ones closest to the real axis, governing the exponential decay at long time scales, i.e. τ , and the falloff of correlations at large distances, i.e. ξ respectively. We find from the quasinormal mode analysis that 1/ξ 2 ∼ | |, which implies ν = 1/2, consistent with the results of [35]. By studying the motion of the leading poles as T → T c , we can also deduce that τ 0 = 2.02±0.01 and ξ 0 = 0.28 ± 0.02 if the critical point is approached from above and τ 0 = 0.89 ± 0.01 and ξ 0 = 0.39 ± 0.01 if it is approached from below. For more details, we refer the reader to

IV. FAR FROM EQUILIBRIUM DYNAMICS
We simulate cooling of the superconducting ring by numerically evolving the bulk equations (Eqs. (20) in Methods) in the black-hole background with changing temperature (11).
We implement a piece-wise linear protocol starting at an initial temperature of T i with corresponding i , and cooling the system at finite rate according to (t) = −t/τ Q through the critical point, until the temperature T f with corresponding f is reached. We implement the temperature ramp by changing the dimensionless ratio T C = C/2πz h , while holding the density ρC fixed. For a precise definition of these quantities the Reader may consult the Methods section. In order to allow the system to break the symmetry dynamically, one needs to add fluctuations to the classical Einstein equations, which encode the field theory in the large-N limit. To achieve this, we introduce noise into the evolution in a manner consistent with the fluctuation-dissipation theorem for the horizon temperature T H . Thus we sample the boundary values, ψ(z, t, φ) z=z UV , of the scalar field (its real and imaginary parts) from a Wiener process, which ensures that their average values vanish but give rise to a non-vanishing correlation. We treat the amplitude of the noise as a phenomenological parameter, but it would be enlightening to determine its precise form in the future, by deriving the relevant fluctuation-dissipation relation from bulk quantum gravity (for work in this direction, see [37][38][39]). Before the linearly decreasing temperature ramp we allow the system to thermalize for a fixed amount of time. In order to solve the nonlinear partial differential equations determining the evolution of the bulk fields, we use a Chebyshev grid in the holographic ('radial') direction z and a Fourier decomposition in the periodic direction φ. We use a fixed step to evolve in time. Each quench is started at some >ˆ for a given τ Q and is stopped at a time t L , when the condensate, averaged over the ring, has reached a fixed fraction of the equilibrium value at (t L ). At this point the winding number W is recorded together with the lag time t L . Both of these quantities are predicted in the KZM scenario to follow the universal scaling relations (4) and (9) once averaged over noise realizations. The winding number W is extracted from the discretized φ direction (denote the grid points by As in the case of the continuum definition above, the sum as a whole is gauge invariant under single-valued gauge transformations. We find that it stabilizes soon after the ordered phase is formed (see Fig. 4). Winding number becomes a good observable by the time we stop the evolution to extract both t L and W , when it is frozen in at integer values. From this point on it is no longer susceptible to noise or late-time equilibration dynamics of the ordered phase, as can be seen in Fig. 4. An example condensation process for τ Q = 12 is shown in Fig. 2, where boundary and bulk physics leading to a W = 3 configuration is illustrated.
We find very good agreement between the full simulations and KZM predictions based on the equilibrium critical exponents deduced from our independent quasinormal mode analysis.
From ν = 1/2 and z = 2 it follows that |W | ∝ τ −1/8 Q and t L ∝ τ 1/2 Q , while the simulation results in The quoted uncertainties give a single standard deviation from the fitted value. The discrepancy in the value of t L is likely a result of the fact that the lag time does not vanish at very small quench times, but rather saturates to a finite value, so that the simple scaling form is no longer a very good fit for the rapid quenches at the fast end of our window of τ Q .
In summary, we find good agreement with universal KZM values for the scaling of t L as well as the dispersion of winding number σ (W ). Matching the prefactors in (7)  Evidently the numbers in our simulation are in rather good agreement, compared to past simulations where mismatches by factors of O(10) were not uncommon. Provided that the numerics is good enough, the degree of agreement depends on the microscopic dynamics.
Furthermore our value for τ 0 is in line with the mismatch encountered in previous work [24].

V. DISCUSSION AND OUTLOOK
The work reported here constitutes the first demonstration of universal KZ scaling of defects spontaneously created in the far-from equilibrium dynamics of a strongly-correlated field theory without quasiparticles. Holography allowed us to map the time-dependence of this system to the evolution of a set of nonlinear partial differential equations together with stochastically sampled boundary conditions. Our setup is comparable to a holographic version of the Gross-Pitaevskii equation, applicable for systems without quasiparticles.
We introduced stochasticity by sampling boundary conditions from a Wiener process characterized by a phenomenological parameter α. This is natural, because it constitutes only a minimal extension of the usual holographic dictionary, which tells us how to map bulk to boundary quantities, but has the limitation that there is a free parameter. While we found that the universal results concerning the scaling law as a function of the cooling rate are very robust to changes in α, and the parameters τ sim 0 and ξ sim 0 of our study show only weak dependence, it would nevertheless be desirable to derive it from first principles. This means deriving the fluctuation-dissipation relation associated with the thermal Hawking radiation of the bulk black hole, which will force us to take into account quantum effects in the bulk gravity model along the lines of [37][38][39].
Our scaling results are consistent with exponents obtained in the mean-field approximation of the boundary theory, even though the theory is not described by a simple Landau-Ginzburg (LG) model, as evidenced e.g. by a violation of the standard LG relation between correlation length above and below the transition (see Fig. 3 and its caption). The meanfield like scaling results are a consequence of working in the classical gravity limit, which, however, does manage to capture some features beyond mean-field. A tempting next goal would be to investigate quench dynamics in theories that do not result in mean-field scalings. A much more ambitious (but perhaps not completely unrealistic [40]) goal would be to study strongly coupled theories directly relevant to condensed matter physics.
An open question left unanswered by the present investigation is the origin of the saturation we observed. There are several plausible explanations, and we hope to return to a detailed investigation of the reasons for the saturation in the model we have presented here in forthcoming work.
Another interesting future project would be to study the effect of the noise in more detail.
In addition to the first principles derivation of the properties of the Wiener process we have already noted, two distinct phenomenological tacks can be considered. We have already begun to investigate the effect of the amplitude of the Wiener process (our parameter α) applied at the boundary. Preliminary conclusion is that -at least within the range we have explored -the quantities we have followed are insensitive to α, with the dependence likely to be sublogarithmic. More detailed characterization of this dependence would be desirable.
Furthermore, one can consider noise applied throughout the AdS interior (rather than only on the boundary). This raises the possibility of seeding topological defects in the interior, and may pose the question of the relation between them and the behavior of the field on the boundary.

Acknowledgements:
We would like to thank Tarek

A. Bulk action and Details of the Holographic Mapping
The full action for the AdS 3 holographic superconductor reads where D µ ψ = ∇ µ ψ − iψ and we choose the case m = 0. We denote bulk indices running over z, t, φ by x µ and boundary indices, running over t, x, with x i . Taking q → ∞ has the result of decoupling the gravity part from the matter part [26], tantamount to the 'probe limit' discussed in this article. One fixes a solution of the gravity equations, to be discussed shortly, and treats the dynamics of the matter fields separately. Hence the effective bulk action from which the equations of motion follow is simply Maxwell-scalar theory in a nontrivial background, that is, the equations of motion are for the current J µ = i(ψ∇ µ ψ * − ψ * ∇ µ ψ). Covariant derivatives are taken with respect to a fixed background metric, which we take to be the metric of the BTZ black hole, written in ingoing Eddington-Finkelstein coordinates The length sets the AdS curvature scale; its role is to fix the number of degrees of freedom 'N ' of the dual theory, which must be large for classical gravity to apply. At leading order in large N this quantity scales out of the equations we solve. The metric function takes the For the scalar field one finds the asymptotic expansion The interpretation of ψ (0) (t, φ) is that it sources the symmetry-breaking operator O(t, φ), while ψ (2) (t, φ) gives its expectation value Therefore the requirement that the symmetry be broken spontaneously means that we must set the source ψ (0) (t, φ) to zero for all time. This translates into a homogenous Dirichlet boundary condition on the field ψ(t, z, φ) in the UV, The vector field in AdS 3 is more subtle [29]. Its asymptotic behavior is given by where the vector j µ (t, φ) is an external current in the boundary theory, not to be confused with the bulk current J µ (t, z, φ). We introduced the scale Λ to make the argument of the logarithm dimensionless. We shall see that our chosen boundary condition on A µ is independent of this scale. Note that in the normal phase the solution is j µ log(z/z h ) and Λ = z h is enforced by regularity at the horizon. It is convenient to work in axial gauge so that A z = 0 and thus the current j µ (t, φ) has components only in the field theory directions, j i (t, φ).
Note that this is not the same as choosing axial gauge in the Schwarzschild like coordinate system (D1), in which the bulk metric is diagonal. In the latter choice of coordinates, the equations of motion imply the equation Thus the current is conserved in the absence of a source for the operator O(t, φ) (We have denoted the current evaluated in the Schwarzschild like coordinates as j i S ). From the point of view of the boundary field theory this is simply the Ward identity for the one-point function of the current. Operationally, the conservation condition follows from the z−component of the Maxwell equations, which gives rise to a constraint. Turning to the gauge field, the boundary condition is where Π µ A is momentum conjugate to A µ with respect to z slicing This is a Neumann boundary condition. Thus we are free to fix j i (t, φ) subject to the conservation condition, and leave a i (t, φ) free to fluctuate. As stated above, the scale Λ drops out from the boundary condition (27). This choice corresponds to a dual vector operator of dimension ∆ = 1, the right dimension for a dynamical gauge field in the boundary theory.
Indeed the residual gauge transformation preserving axial gauge (λ res (t, z, φ) = λ(t, φ)) acts on this as a standard field-theory gauge transformation so that a i (t, φ) is a bona-fide fluctuating gauge field. We fix a constant background charge density, so that the current has the only non-vanishing component, j t = ρ.

B. Details on Numerics
The simulations in this article were performed on a pseudo-spectral spatial grid comprised of 21 Chebyshev points in the radial direction and 111 plane waves in the angular direction of the boundary. Writing the system in ingoing Eddington-Finkelstein coordinates allows us to employ a characteristic evolution scheme that is linear in time derivatives. In order to propagate in time we used forward Euler as well as explicit fixed-step RK4 integration with mutually agreeing results.
For each run we start the system at the initial time slice in the normal phase, defined by setting to zero the field ψ and giving the gauge field a non-trivial time component j t = ρ.
We then evolve forward in time, updating the noise according to the rule (16), implemented as a discrete Wiener process. Since the scaling quantity of interest pertains to the standard deviation of winding number σ(W ), we average over O(10 2 ) noise realizations for each value of τ Q .

Appendix A: Holographic Renormalization
In axial gauge, and setting ψ (0) = 0, the equations of motion (20) have solutions with asymptotic behavior We regularize the on-shell action by introducing a finite cutoff z = , so that we now have wheren is the outward pointing unit normal to the boundary and γ = det(γ ij ) the determinant of the induced metric. The omitted terms depend only on the boundary values j i , a i , ψ (2) and vanish in the limit → 0. The divergent terms can be cancelled by adding a counterterm As shown in [30] this counterterm becomes a manifestly local boundary quantity when written in terms of a dual scalar field S, obtained from the two-form gauge field F via F = dS.

Appendix B: Evolution Scheme
For the simulations reported on in this article we employed a characteristic evolution scheme for the Maxwell and scalar fields. After gauge fixing A z = 0, denoting A t (t, z, φ) = T (t, z, φ), A x (t, z, φ) = X(t, a, φ) and writing ψ(t, z, φ) = a(t, z, φ) + ib(t, z, φ) we have four evolution equations and one constraint equation. In the numerical evolution it proves convenient to work with rescaled fields, rather than the 'bare' ones appearing in the action (19). We furthermore subtract the leading log terms. Thus definingT = z (T + ρ log z) and X = zX we find the equations where the fields Φ i are all fields other thanT . The sources S i and S T depend non-linearly on the fields as well as their spatial derivatives. In addition, the radial (z) component of the Maxwell equation gives the constraint equation After the rescaling, the boundary conditions onT andX are noŵ Note that Eqs. (B1) are linear equations for the time derivatives suggesting the following evolution scheme. Assume that at time t = t n the values of the fields a, b, X are known.
We can now use the constraint equation to solve forT (t n , z, φ) and then solve the linear equations for Φ i t to obtain the time derivatives of a, b,X. We can then use our favorite time evolution scheme to obtain the values of a, b,X on the next time slice t n+1 . We have used explicit RK4 integration as well as simple forward Euler with good results. For 111 Fourier modes in the spatial boundary direction, the step size was taken to be 0.003/N z , where N z is the number of grid points in the radial direction. The scaling analysis of Figs. 5 was obtained on a grid of 21 points in the radial direction and 111 Fourier modes in the annular direction choosing a characteristic size of the ring C = 50 . It would be desirable to repeat the analysis for a larger ring with higher spatial resolution, which is likely to require more significant computing power, especially if higher statistics on noise realizations is desired.
The linear equations for ∂ t Φ i n on a given time slice can then be solved efficiently and in parallel for each spatial grid point φ j . For the above-mentioned spatial grid we choose a fixed time step of ∆t = 0.003/N z . In order to achieve stable long-time evolution we filter out high-frequency modes in the φ direction using the Orszag 2/3 rule every five time steps. This can be achieved very efficiently using FFT and iFFT. Due to the presence of logarithmic terms in the asymptotics of the fields one does not expect spectral accuracy, and this is borne out in preliminary convergence tests. We found that our evolution scheme is stable even in the presence of stochastic boundary conditions simulating thermal noise and gives accurate results (for example by comparing the dynamical solutions at late time to the equilibrium results obtained in the standard way from solving ODEs). We now turn to a detailed description of our implementation of stochastic boundary conditions. Under classical evolution, corresponding to the classical large N limit of the dual field theory, fluctuations are suppressed and the symmetry of the order parameter cannot be dynamically broken. Said differently, even though below T c the phase with O = 0 is unstable, there is nothing in the classical evolution equations to push the order parameter off its precarious perch on top of the potential. Taking our lead from the literature on BEC dynamics using the stochastic Gross-Pitaevksii equation [20] we add fluctuations by sampling our boundary conditions from a thermal noise distribution, that is from a Wiener process.
Thus the boundary conditions for a(t, z, φ) and b(t, z, φ) should be modified. Instead of imposing the strict Dirichlet boundary condition (24), we only impose this condition on where · denotes noise average and ψ (0) = a (0) + ib (0) . Usually one determines α from a fluctuation-dissipation relation as α = 2ηT c , where η is a damping parameter. In this work we treat α as a phenomenological parameter. We found that varying α, even by several orders of magnitude has little influence on the scaling results, but we do find a weak dependence of the absolute magnitude of t L on α. Clearly it is desirable in future to determine α from a first-principles holographic calculation.
In practice the noisy boundary condition is realized by sampling each spatial boundary point from an independent normal distribution of zero mean and unit variance N (0, 1). That is, we set for each boundary point φ i . Note that one cannot choose both the boundary value of j i and ψ (0) independently, since they are constrained by the current Ward identity (26). This equation is consistent only for ψ (0) = 0, so in what sense can one choose the noisy boundary condition above? This can be seen by considering the noise as a small perturbation j i + δj i and ψ (0) + δψ (0) satisfying By solving the z component of the Maxwell equation near the boundary we automatically impose this constraint at each time step, that is the noise fluctuations (C2) also imply a fluctuation in the current δj i , such that (C3) is satisfied.
Note in particular that (C2) implies that the phase is a random variable. In our numerics we find that it is sufficient to update the noise boundary condition at larger intervals, say every 100 time steps. We have also run the simulations updating the noise at every time step, as well as every 10 time steps, again with no apparent impact on the results -provided the noise amplitude is adjusted in accordance with (C2). Clearly there is a lower limit on the sampling frequency for which this statement is true -consider, e.g., the extreme case of only updating the noise once or twice during the entire simulation. However our results show no detectable dependence on sampling frequency, which implies that we stayed far away from this lower limit throughout. by a factor of ten. The inset shows a blowup of the 'knee' region where the condensation process happens. We see that this process is largely insensitive to the noise amplitude.
Since the equations are linear in ω, but quadratic in k, these admit an expansion [35] in small frequency and momentum α = α (0,0) (z) + ωα (1,0) (z) + k 2 α (0,2) (z) + · · · , β = β (0,0) (z) + ωβ (1,0) (z) + k 2 β (0,2) (z) + · · · , where for z → z UV . In this expression a superscript s denotes source behavior, while v denotes expectation value. Thus Green functions of the dual operator have the small ω, k behavior as claimed in (12) above. More precisely, the Green functions of the operators dual to α and β mix, but the eigenvalues of the matrix of correlation functions will be of the functional form (D6). This establishes analytically that the dynamical critical exponent is z = 2. In order to determine ξ, and thus ξ 0 and ν, we solve the system of equations (D3) numerically, using a finite difference discretization. Results of this analysis are described in more detail in Sec. III above, and summarized in Fig. 3.

Broken Phase
We now describe the more involved calculation in the broken phase, where for ω = 0 and k = 0 more mode mixing occurs. If we were to go beyond the probe approximation, the sound channel [41], holographically encoded as a scalar fluctuation of the metric, will also contribute. This means that in general the condensate-condensate two points function below T c no longer takes the simple form (12).
Below the critical temperature, i.e., when the background contains a nontrivial scalar field, we have to take into account the modes e −iωt+ikx {a t (z), a z (z), a x (z), α(z), β(z)}, where a µ (z) are perturbations of the gauge field and α(z) and β(z) are perturbations of the real and imaginary parts of the complex scalar field respectively. Not all of these modes are physical, since we can generate an infinitesimal perturbation by acting on the background with a gauge transformation e −iωt+ikx λ(z). This generates the perturbation modes {δa t , δa x , δa z } = e −iωt+ikx {−iωλ(z), −ikλ(z), λ (z)} in the gauge sector, as well as a perturbation of the imaginary part of the scalar δβ = e −iωt+ikx a 0 (z)λ. We cannot generate a real part of the scalar perturbation in this way, because the background value a 0 (z) is purely real. Thus we have three gauge-invariant perturbations