Universal Statistics of Vortices in a Newborn Holographic Superconductor: Beyond the Kibble-Zurek Mechanism

Traversing a continuous phase transition at a finite rate leads to the breakdown of adiabatic dynamics and the formation of topological defects, as predicted by the celebrated Kibble-Zurek mechanism (KZM). We investigate universal signatures beyond the KZM, by characterizing the distribution of vortices generated in a thermal quench leading to the formation of a holographic superconductor. The full counting statistics of vortices is described by a binomial distribution, in which the mean value is dictated by the KZM and higher-order cumulants share the universal power-law scaling with the quench time. Extreme events associated with large fluctuations no longer exhibit a power-law behavior with the quench time and are characterized by a universal form of the Weibull distribution for different quench rates.


Introduction
The dynamics across a continuous phase transition is a paradigmatic scenario of spontaneous symmetry breaking in which adiabaticity inextricably breaks down. In any finite time scale, a quench from the high-symmetry phase to the lower-symmetry phase is governed by critical slowing down and the effective freezing of the system. Facing a degenerate manifold, causally separated regions of the system make disparate choices of the broken symmetry that result in the formation of topological defects [1,2]. The characterization of the latter depends on the topology of the vacuum manifold. In the formation of superconductors and superfluids, entailing U (1) symmetry breaking, vortices with quantized flux appear [3,4]. Spontaneously formation of vortices was observed in experiments with neutron-irradiated superfluid Helium [5,6].
The mean number of topological defects generated in the course of a phase transition is predicted by the KZM to follow a universal power law with the rate at which the phase transition is crossed. The verification of this prediction has been the subject of a long-time quest [7]. The validity of the KZM is not only supported by theoretical models and numerical simulations but has been established in a variety of experimental platforms ranging from colloids [8] to quantum simulators [9,10].
In strongly coupled systems, the validity of KZM cannot be taken for granted. A natural framework to account for strong coupling is provided by holography [11,12]. In this context, the spontaneous current formation in a superconducting ring [13][14][15] is well described by the KZM prediction [2,16]. However, deviations from KZM have been predicted and the power-law scaling of the mean number of topological defects is expected to be modified by logarithmic terms of the quench rate [17]. In the laboratory, pioneering experiments on the spontaneous vortex formation in the light of the KZM were restricted to the weakly interacting regime, accessible with Bose-Einstein condensates [18][19][20] and ferroelectrics [21]. Remarkably, recent progress has allowed probing the strongly-interacting regime using a unitary Fermi gas [22]. The measured density of defects was found to be compatible with the KZM scaling laws as in the weakly interacting case. Theoretical [17] and experimental [22] results on vortex formation at strong coupling are thus in conflict.
The predictive power of the KZM is restricted to the average number of topological defects. Spatial correlations between topological defects have been discussed in the framework of the Halperin-Liu-Mazenko theory [23,24]. Fluctuations of the number of topological defects have recently been explored in spin chains [10,25,26] and one-dimensional φ 4 theory [27]. These studies have unveiled signatures of universality in the full counting statistics of topological defects that lie beyond the scope of the KZM.
In this work, we explore the statistics of vortices in a newborn holographic superconductor in (2+1) dimensions and show that the full counting statistics of vortices is universal. The mean density is shown to follow the KZM power-law prediction. Fluctuations beyond the mean are probed by low-order cumulants of the vortex number distribution, which are found to exhibit a universal power-law scaling with the quench time. The vortex number distribution is well described by a binomial distribution restricted to even outcomes by the topology of the system, making it possible to probe rare events far away from equilibrium. Large deviations away from the mean vortex number no longer exhibit a power-law behavior and the corresponding extreme value statistics is characterized by a Weibull distribution.

Formation of a newborn holographic superconductor
We simulate the superconducting transition from a normal metal to a holographic type-II superconductor in two spatial dimensions by implementing a thermal quench in a finite time τ Q . This results in the spontaneous formation of vortices that are pinned. The system is described making use of the gauge-gravity duality and numerical simulations involving a (3+1) dimensional gravity. In this setup, a thermal quench can be effectively simulated by changing the charge density in the boundary of the black hole, see section 3 [28]. The phase transition is continuous and of second-order. Critical slowing down in the proximity of the critical point leads to vortex formation. To characterize the resulting vortex number distribution, we consider a homogeneous system, free from external potentials that can alter the KZM scaling [29].
According to the KZM, traversing the phase transition at finite rate leads to the formation of domains of characteristic length scaleξ = ξ 0 (τ Q /τ 0 ) ν 1+zν , where ν and z are the correlation-length and dynamical critical exponents of the continuous phase transition, and ξ 0 and τ 0 are microscopic constants. Within such domains, the superconductor phase is chosen coherently. According to the geodesic rule [3], when multiple domains merge at a point, there is a chance that the quantized circulation of the superconductor phase φ around that point is non-zero and a multiple 2π. This configuration can lead to the formation of a vortex. Typical values of its vorticity V = 1 2π dγ∇φ ∈ Z are V = ±1. The number of vortices induced by the thermal quench is thus proportional to n ∝ A/ξ 2 where A is the area of the superconductor. As a result, the mean vortex number scales with the quench rate following the universal power law n ∝ (τ 0 /τ Q ) 2ν 1+zν , which is the key prediction of the KZM. This universal scaling quantifies the intuition that fast quenches result in a high number of vortices, the number of which decreases as the rate of the transition is reduced. As the transition is thermal, the number of vortices is not deterministic but it constitutes a stochastic variable described by a probability distribution.
Characterizing the full counting statistics of spontaneously formed vortices across the phase transition is our central goal. As we shall see, fluctuations away from the mean value are universal. Specifically, we show that cumulants of the distribution share with the mean value a universal power-law behavior in the limit of slow quenches, required for scaling theory to hold. In addition, knowledge of the exact vortex number distribution allows us to characterize extreme events associated with large deviations from the mean value.

Holographic setup
We work in the black brane background in Eddington-Finkelstein coordinates The horizon is z h while z = 0 is the boundary in which the field theory lives. In the probe limit we adopt the Abelian-Higgs Lagrangian density as where D µ = ∇ µ − iA µ is the covariant derivative, A µ is the U (1) gauge field, F µν = ∂ µ A ν − ∂ ν A µ is the gauge field strength and Ψ is the scalar field. We take the ansatz as Ψ = Ψ(t, z, x, y), A t = A t (t, z, x, y), A x = A x (t, z, x, y), A y = A y (t, z, x, y) and A z = 0. The equations of motions (EoMs) of these fields can be obtained readily from the Lagrangian density. The explicit forms of the EoMs can be found in the Appendix A. Near the boundary z → 0 the expansion of the fields are (we have set L = 1 and m 2 = −2), From gauge-gravity duality, the coefficients a t , a i (i = x, y) and Ψ 0 can be interpreted as chemical potential, gauge field velocities, and scalar operator source in the boundary, respectively. Their corresponding conjugate variables are achieved by varying the renormalized on-shell action S ren with respect to these coefficients. To get a finite S ren , some counter terms should be added. According to holographic renormalization [30], the counter terms of the scalar field is S ct = d 3 x √ −γΨ * Ψ, where γ is the reduced metric on the z → 0 boundary. We have imposed Neumann boundary conditions for gauge fields as z → 0 in order to get dynamical gauge fields in the boundary [31,32]. Therefore, a surface term where n µ is the normal vector perpendicular to the boundary, should also be added to have a well-defined variation. Therefore, the expectation value of the order parameter on the boundary, O = Ψ 1 , can be obtained by varying the finite renormalized action S ren with respect to the source term Ψ 0 . The conservation equation of the charge density and current on the boundary is To study spontaneously symmetry breaking of the order parameter, we set Ψ 0 = 0 in the z → 0 boundary. Besides, the Neumann boundary conditions for gauge fields can be imposed from the above conservation equations. In the spatial (x, y)-directions we impose the periodic boundary conditions for all the fields. At the horizon, we set the time component of gauge fields as A t (z h ) = 0, and the regular finite boundary conditions for other fields.
To drive the system out of equilibrium, we quench the charge density ρ on the boundary while fixing the temperature of the black hole, which was commonly implemented in holographic superconductor settings [28]. The mass dimension of black hole temperature T is one, while the mass dimension of the charge density ρ is two. Thus, T / √ ρ is a dimensionless quantity. Therefore, decreasing the temperature is equivalent to increasing the charge density. In order to have a linear quench of temperature across the critical point, one can quench the charge density ρ as

Numerical Scheme
Before quench, we thermalize the system by adding random seeds into the system in the normal state. The random seeds of the fields are added in the bulk by satisfying the statistical distributions s(t, x i ) = 0 and s(t, In principle, h cannot be too large since the seeds serve as fluctuations to thermalize the system. The system is quenched linearly from T i = 1.4T c to T f = 0.8T c with T c the critical temperature. We evolve the system by using the fourthorder Runge-Kutta method with a time step ∆t = 0.02. In the radial direction z, we use the Chebyshev pseudo-spectral method with 20 grids. Since along (x, y)-directions periodic boundary conditions are imposed, we adopt the Fourier decomposition in the (x, y)directions. The size along (x, y) is 50 × 50, and the number of grids are 201 × 201. Filtering methods are implemented following the rule that the uppermost one-third Fourier modes are removed [33]. We count the number of vortices as the average order parameter just arrived at its equilibrium value. Due to the large dimensions of the system (3+1-dim), the time cost in large quench time is considerable. For instance of τ Q = 4000, each trajectory of the simulation will cost more than two hours. Therefore, collecting statistics with 1000 trajectories required about three months of running time. Due to the large consumption of time, we have limited the number of trajectories for each quench time to 1000.

Vortex counting statistics
We consider the formation of a newborn superconductor in two spatial dimensions with periodic boundary conditions in the (x, y) spatial directions. The topology of the system is that of a torus T 2 with zero Euler characteristics χ(T 2 ) = 0. As a result of the Poincaré-Hopf theorem [34], the total vorticity of the superconductor equals χ and vanishes identically. By contrast, in the case of open boundary conditions, the vortex number would be unrestricted by the topology of the system and could take both even and odd values. The regime of parameters studied is such that no vortex with vorticity other than V = ±1 is observed, and the number of vortices and anti-vortices is thus balanced, see Figure 1.
We focus on the total number of vortices, regardless of their vorticity. By numerically solving the stochastic dynamics in the holographic setting, an ensemble of realizations is used to collect statistics and build a histogram for different values of the vortex number. For increasingly fast quenches (τ Q = 4000, 2000, 1000 and 20) the distribution shifts to higher mean values, while simultaneously broadening. By contrast, at the onset of adiabatic dynamics the distribution narrows down and becomes increasingly asymmetric, as the fully adiabatic regime with n = 0 acquires a significant probability P (0); see Fig. 2. The vortex number statistics is found to be precisely described by a binomial distribution restricted to even outcomes, which we refer to as an even-binomial (EB) distribution in the following. The later results from assuming that the probability for vortex formation at the merging point between adjacent domains occurs with probability p, while no vortex is formed at such location with probability (1 − p). With only two possible outcomes this process can be described as a single Bernoulli trial. The total number of vortices is the result of a number of trials N ∼ A/ξ 2 . We next assume that vortex formation at different locations unfolds as a result of uncorrelated events that can be described by N independently and identically distributed (iid) random variables. The probability to observe a given vortex number n is thus  [35] guarantees that the statistics becomes even-Poissonian (EP) where λ = N p is a parameter and the mean reads n = λ tanh(λ). Additional supporting evidence of the excellent agreement between the P (n) values extracted from the histogram and those predicted by the distribution (4.2) is provided in the Appendix B for different values of n.
To characterize universal signatures in the full counting statistics, we describe the scaling of the low-order cumulants of the distribution as a function of the quench time. Given the Fourier transform of the distributionP (w) = E[e iwn ], where w is variable conjugated to n, cumulants κ p are defined through the expansion logP (w) = ∞ p=1 κ p (iw) p /p!. The first cumulant equals the average vortex number n and is predicted by the KZM. The second cumulant equals the variance κ 2 = Var(n) while the third one is related to the skewness Skew(n) of the distribution through the identity κ 3 = Skew(n)κ 3/2 2 . In the Poissonian limit for large average vortex numbers, all cumulants approach the first, thus inheriting the universal power-law scaling as a function of the quench time dictated by KZM. This prediction is explicitly verified in Fig. 3 where the first three cumulants are plotted as a function of the quench time using about 1000 trajectories. This sampling size is limited by the computational cost (see Methods). The first two cumulants show saturation at a plateau for fast quenches, followed by a universal power-law behavior for longer values of the quench time. The scaling of the first cumulant is dictated by the KZM to follow a power-law n ∝ τ −1/2 Q for mean-field critical exponents ν = 1/2 and z = 2 in two spatial dimensions. A fit to the data shows that n ∝ τ −0.518±0.0243 Q . The even Poissonian distribution for a large mean number of vortices predicts a power-law scaling of the vortex number variance κ 2 = Var(n) ∝ τ −1/2 Q , i.e., equal to the KZM scaling for the mean. The large fluctuations observed in the third cumulant are expected for the number of trajectories considered and their suppression would require one to increase the number of trajectories by one to two orders of magnitude.

Large fluctuations
In what follows we turn our attention to extreme statistics associated with rare events, that can be estimated efficiently with the available sample size. Fluctuations far away from the mean vortex number can be expected to be sensitive to defect-defect interactions. Their characterization can be achieved by adding the contribution from the tails of the distribution. In the even-Poissonian limit, the cumulative probability for P EP (n ≤ r) and P EP (n ≥ r) are in which 1 F 2 is a hypergeometric function. For the values of r = 6, 10, 20 and 30, the cumulative probability is shown in Fig. 4 as a function of the mean number with an excellent agreement between the numerical data and the prediction for the even Poissonian distribution. As n is predicted by the KZM, it can be fitted to power the scaling of n ≈ 333.3 × τ −0.518 Q in the range of quench times τ Q ∈ [1000, 4000] (see Fig. 3), together with Eqs. (5.1) and (5.2), to quantify the quench time dependence of rare events, shown in Fig. 5.
For further characterization, we resort to large deviation theory and characterize the distribution of the maxima in long sequences of realizations. According to the Fisher-Tippett-Gnedenko theorem [36], the extreme values of the iid variables satisfy the generalized extreme value (GEV) distribution, with location parameter µ, scale parameter σ and shape parameter ξ. The GEV distribution includes as limiting cases of the Weibull (ξ < 0), Gumbel (ξ → 0) and Fréchet (ξ > 0) laws.  To determine the relevant GEV distribution, we use the block maxima method. Data from the stochastic trajectories is partitioned in blocks. In each block, the maximum vortex number n M is found and from the ensemble of blocks the probability Prob(n M < x) is determined. We analyze the statistics of large vortex number deviations for slow quenches in the universal scaling regime in Fig. 6, that show the probability density function (PDF) and the cumulative distribution function (CDF) for their GEV distributions in different groups. Both PDF and CDF are shown as a function of the variable y = x−µ σ . Specifically, for τ Q = 1000, data is partitioned into 100 groups (top row in Fig. 6) and 200 groups (bottom row in Fig. 6). The corresponding parameters are µ = 13.056, σ = 1.666, ξ = −0.149 and µ = 11.892, σ = 1.943, ξ = −0.184, respectively. Thus, GEV is described by the Weibull distribution, which has reflecting the existence of an upper bound to the vortex number. As shown in the Appendix E, the validity of the Weibull distribution further extends away from the universal scaling regime and characterizes deviations also in the saturation regime associated with fast quench times.

Discussion
The characterization of the full distribution of topological defects generated across a phase transition is expected to have wide applications, ranging from condensed matter to quantum simulation and computation, and cosmology. Experimental efforts to date have focused on the universal power-law dependence of the mean number of defects with the quench time, which is successfully predicted by KZM. Our results show that fluctuations away from the mean exhibit universality but are no longer captured by the KZM power-law scaling. The full counting statistics of vortices in a strongly coupled superconductor follows a universal binomial distribution. The distribution of maxima in a sequence of realizations is captured for the Weibull law in large deviation theory. The dependence of the tails of the distribution with the quench time dictates the suppression of topological effects, near or far from the onset of adiabatic dynamics. This dependence is thus crucial to analyze rare events associated with profusion or absence of topological defects. The complete suppression of topological defects is sought after in the preparation of novel phases of matter in quantum simulators and finite-time quantum annealing and quantum optimization. It is also of relevance in a cosmological setting, as the observation of cosmic strings predicted by the KZM remains elusive.

A Equations of motion for a holographic superconductor
In the probe limit, the equations of motions for Ψ and A µ read, in which (.) represents imaginary part. In explicit form, these equations read in which (.) represents imaginary part. In explicit form, these equations read where Φ = Ψ/z. The above five equations are not independent, and their L.H.S. satisfy the following constraint equation Therefore, there are four independent equations for four fields, Φ, A t , A x and A y . This also implies that our choice of the gauge A z = 0 is viable for the setup of the system.

B Properties of binomial and Poissonian distributions restricted to even outcomes
The even-binomial (EB) distribution is obtained by restricting to even outcomes the binomial distribution and is where N is the number of domains with broken symmetry, p is the success probability to form a vortex, n represents a given number of vortex and belongs to non-negative even integers, and A = 1+(1−2p) N 2 is the normalization constant. The first three cumulants of EB distribution are These cumulants satisfy a recursion relation such as κ q+1 = p(1 − p) dκq dp . In the limit of N → ∞ and keeping the parameter N p = λ finite, we get the even-Poisson (EP) distribution C Full counting statistics of vortices: numerical simulations for a holographic superconductor and even-Poissonian distribution In Fig. 8, we compare the numerical results of the probability distribution P (n = 0, 2, 4, . . . , 30) of a given vortex number n with respect to theoretical value of the EP distribution, and the (unrestricted) Poisson distribution, which is a limit of (unrestricted) Binomial distribution with N → ∞ and keeping the parameter N p = λ finite. The agreement between the numerical simulations and the EP distribution is excellent provided sufficient statistics.  The probability P (n = 0) to observe no vortices at all is a rare event away from the adiabatic limit and can be estimated from the EP distribution. Using Eq. (B.5), it reads P EP (n = 0) = sech(λ). (C.1) Numerical results for P (n = 0) with respect to n are shown in Fig. 9 with P (n = 0) ≈ 3.998 × sech( n ). The value of the prefactor is found to decrease with an increasing number of trajectories of simulations. This implies that the error between the numerical results and the theoretical EP distribution is due to the limited sampling statistics in the simulations. By increasing it, the numerical results are expected to approach the theoretical EP distribution. In the main text, we show the cumulative probability of even-Poissonian distribution for P EP (n ≤ r) and P EP (n ≥ r) as where 1 F 2 is a hypergeometric function and · is the floor function. In Figure 10, the numerical results of the cumulative probabilities are shown to match very well the theoretical predictions for a broad range of rates at which the transition is crossed, ranging from slow quenches with τ Q = 3000 to the fast quench limit with τ Q = 20.

D.1 Chernoff bound
The Chernoff bound [37] can be used to derive exponentially decreasing bounds on the tail distributions of vortex numbers. In its looser form, the Chernoff bound can be written as In Fig. 11 we plot bound of lower tail and upper tail for δ = 1, 2, 3. The numerical results satisfy the Chernoff bound very well. The latter can thus be used to capture the dependence on the quench time of the large fluctuations away from the mean, associated with the tails of the distribution.

E Extreme value distribution of vortex numbers
According to the Fisher-Tippett-Gnedenko theorem [36], the extreme maximal values of the independently and identically distributed (iid) variables satisfy the generalized extreme value (GEV) distribution, in which, µ is the location parameter, σ is the scale parameter and ξ is the shape parameter. Note that ξ(x − µ)/σ + 1 > 0 and zero otherwise. The above GEV distribution function G(x; µ, σ, ξ) is the cumulative density function (CDF), whose corresponding probability density function (PDF) can be written as If ξ < 0, the GEV distribution is called Weibull distribution which is upper bounded. if ξ = 0, the GEV distribution is called Gumbel distribution which has a light tail. Finally, if ξ > 0, the GEV distribution is called Fréchet distribution which has a heavy tail and a lower bound.
In practice, to analyze the extreme value distributions for iid variables, it is customary to separate the data into several groups (or blocks), and then proceed to identify the maximum in each group. The final list of maxima will tend to satisfy the above GEV distribution. This method is called 'Block Maxima' method, and we adopt it to study the maximum value distributions for the vortex numbers in numerical simulations. There are some arbitrary choices in the partition of the data. We partition the data into more than 100 groups, which is sufficient for the observed vortex-number maxima distribution to be identified with the GEV.
The extreme maximal value distributions of the vortex number for a slow quench (such as τ Q = 1000) are shown in the main text. Here, we show the PDF and CDF of the maximum values for a fast quench (τ Q = 20) in Fig. 12. Specifically, for τ Q = 20 we have 11655 numerical data of the vortex number. They are partitioned into 777 groups (top row in Fig. 12) and 111 groups (bottom row in Fig. 12). Both PDF and CDF are shown as a function of the variable y = x−µ σ . In the top row of Fig. 12, the parameters are µ = 32.841, σ = 2.382, ξ = −0.113; while in the bottom row the parameters are µ = 36.558, σ = 2.114, ξ = −0.090. The analysis of the data based on these partitions shows that the GEV distribution of the maxima for τ Q = 20 belongs to Weibull distribution, which means there is an upper bound. Figure 12. Probability density function (PDF) and cumulative distribution function (CDF) of extreme maximum values of the vortex number spontaneously generated in the limit of fast quenches with τ Q = 20. Independently of the partitioning size, he data is well described by a Weibull distribution with an upper bound.