Topological defects as relics of spontaneous symmetry breaking from black hole physics

Formation and evolution of topological defects in course of non-equilibrium symmetry breaking phase transitions is of wide interest in many areas of physics, from cosmology through condensed matter to low temperature physics. Its study in strongly coupled systems, in absence of quasiparticles, is especially challenging. We investigate breaking of U(1) symmetry and the resulting spontaneous formation of vortices in a (2 + 1)-dimensional holographic superconductor employing gauge/gravity duality, a ‘first-principles’ approach to study strongly coupled systems. Magnetic fluxons with quantized fluxes are seen emerging in the post-transition superconducting phase. As expected in type II superconductors, they are trapped in the cores of the order parameter vortices. The dependence of the density of these topological defects on the quench time, the dispersion of the typical winding numbers, and the vortex-vortex correlations are consistent with predictions of the Kibble-Zurek mechanism.


Introduction
Critical dynamics in strongly coupled non-equilibrium phase transitions is one of the most interesting and important problems in modern physics [1]. The conventional quasiparticle approaches do not apply in this case. We study the critical dynamics of the superconducting phase transition, focusing on the formation and evolution of topological defects, by utilizing the AdS/CFT correspondence [2]. Generation of topological defects is expected in such transitions and can be used to test the Kibble-Zurek mechanism (KZM) [3][4][5][6].
The basic idea of KZM is that, as a system approaches the critical point starting in the symmetric phase, its dynamics undergoes critical slowing down reflected in the divergence of the relaxation time. As a consequence, different domains of the system cannot communicate with each other, and select how to break the symmetry independently. The dimension of these domains is limited by the size of the sonic horizon -by how far the relevant sound can propagate in the near-critical time interval. These independent choices of the broken symmetry lead to irreconcilable differences -formation of topological defects can be expected, and has been observed. Numerical simulations [7][8][9][10][11][12][13] as well as experiments in liquid crystals [14][15][16], 3 He superfluids [17,18], Josephson junctions [19][20][21][22], thin-film superconductors [23][24][25] and quantum optics [26] have shown results consistent with KZM (for reviews, see [27,28]).
Gauge/gravity duality has been employed to study strongly coupled systems to bypass the difficulties caused by the absence of quasiparticles [29][30][31][32], for a review see [33]. Previous studies on KZM in holographic superfluids and 1D superconducting loop can be found in [34,35]. In these two holographic studies, scaling laws between number density JHEP03(2021)136 of defects and the quench time were found to match KZM in slow quench regime. Other papers on holographic KZM are [36,37].
We use the AdS/CFT correspondence to examine the U(1) gauge symmetry breaking in the (2+1)-dimensional boundary system, and study the non-equilibrium critical dynamics in a strongly coupled field theory. The spontaneous generation of magnetic fluxons (with quantized fluxes) is observed for the first time using holographic numerical simulation, figure 1. Spatial distributions and other characteristics of these topological defects are investigated. This includes their density as a function of the quench time, the dispersion of the typical winding numbers in the resulting superconductor, the flux distribution, and the correlation between the charges and locations of topological defects. We conclude that the resulting defects exhibit the short-range and nearest neighbor vortex-antivortex correlations consistent with the predictions from KZM.
Previous work on holographic KZM did not realize the spontaneously generated magnetic fluxons. For example, in [34] the authors studied the holographic KZM for vortices in a superfluid without any magnetic fluxons. In [35] the authors studied the holographic KZM in a spatial one dimensional ring which obviously did not have the magnetic fluxons. The authors in [36,37] analytically investigated the KZM scaling laws nearby the critical point of the phase transition, without any numerical simulations of the vortices, not to mention the magnetic fluxons. The key difference of our paper to previous work is that we realize the spontaneously generated magnetic fluxons in the superconductor system. The advantage of the existence of the magnetic fluxons is that the vortices will exhibit the "pinning effect" [38], which means the vortices will move very slowly in space. Thus, the vortices will be very hard to meet and annihilate. Therefore, the slowing down of the movement will be beneficial to counting the number of vortices, measuring the typical winding numbers and evaluating the vortex-antivortex correlations, which are the tasks in our paper.

Results
As observed by Kibble [3], in the early Universe symmetry breaking phase transitions induce local choices of broken symmetry that cannot be coordinated between domains larger than the distance travelled by light since Big Bang till the transition. The resulting a patchwork of domains can then lead to topological defects.
In the laboratory setting speed of light is no longer relevant. However, critical slowing down that occurs in course of second order phase transitions limits the size of the domains. Equilibrium critical exponents can then predict scaling of the defect density as well as other excitations as a function of the quench time [5].

Magnetic fluxons from symmetry breaking
Requirements of minimal free energy and periodicity of the phase θ of the complex scalar field Ψ imply the quantization of the magnetic fluxes generated from U(1) symmetry breaking [38], i.e., Φ = 2πZ, where Z is an integer. Typical configurations of magnetic fluxons and order parameter vortices generated from U(1) symmetry breaking after quench (with JHEP03(2021)136  figure 1(b) has Φ ≈ 2π, which demonstrates the quantization of magnetic flux. All ten positive magnetic fluxons shown in this plot have the average Φ ≈ (1.9864±0.0067)π; while the other ten negative magnetic fluxons have average Φ ≈ (−1.9834 ± 0.0093)π. Therefore, the net magnetic flux of the whole system vanishes within numerical errors, consistent with the fact that there is no external magnetic field imposed for the system. Thus, all magnetic fluxons are spontaneously generated from U(1) symmetry breaking due to KZM. 2 At the bottom of figure 1(c), we show the density plot of the order parameter vortices and the streamlines for the angle of the scalar field phase θ. The arrows indicate the directions of the phase angles. One can read out the positive or negative vorticity of the vortices from the streamlines. Consequently, we see that the locations of positive vortices correspond to positive magnetic fluxons, and vice versa.

Kibble-Zurek mechanism
Near the critical point of a second order phase transition, both the relaxation time τ and the correlation length ξ diverge as, where = 1 − T /T c is the reduced dimensionless temperature (or, more generally, dimensionless distance from the critical point), while ν and z are spatial and dynamical critical exponents. One can usually assume that traverses the critical point approximately linearly in time t Above τ Q is the quench time.
KZM recognizes that at the instantt before the critical point, when the rate of change imposed by the quench is comparable to the system's relaxation time τ , the order parameter will cease to follow or even approximate its equilibrium value. Thus, one obtains Time −t marks the beginning of the non-adiabatic evolution and +t its end. What happens inbetween is occasionally described as a "freeze-out" but a more accurate picture, especially in the cases when the order parameter is underdamped, is based on the idea of "sonic horizon" -distinct parts of the system can still evolve and influence one another, but they can coordinate their choices of broken symmetry only at distances given by how far the relevant sound (associated with the perturbations of the order parameter) can propagate during the time interval (−t, +t) [5]. 1 We actually quench the charge density on the boundary by fixing the temperature, which is equivalently to quench the temperature. Please see the details in the Methods (appendix A). 2 In fact we made 100 times of simulations for each τQ, and the differences of the magnetic fluxes for different runs are very tiny, which are in the range of the order 10 −2 π ∼ 10 −3 π.

JHEP03(2021)136
Now, from eq. (2.1) and eq. (2.3), one can obtain the corresponding correlation length, The new, post-transition order parameter will randomly select how to break symmetry in domains that are this far apart, as they have no time to communicate with one another. This sonic horizon argument [5] leads to domains having the size ∼ξ -same scaling (although different pre-factors [39,40]) as these given by the "freeze-out".
For vortices in two-dimensional space the estimated number density of point defects is; The same scaling laws (although, again, with somewhat different prefactors) follow from the arguments based on the "sonic horizon" paradigm [5,28,39,40]. Equations (2.3), (2.4) and (2.5) can be used to test the validity of KZM in laboratory experiments and in numerical simulations.
Another important feature of KZM is the spatial distribution of the charges of topological defects. KZM predicts that the random choices of the locally broken symmetrye.g., phase of the superfluid wavefunction -will ultimately lead to anticorrelated charges of topological defects. This is in contrast to the possibility that topological charges are distributed at random.
The scalings of the typical winding number W subtended by a loop C with circumference C = 2πr can be used to distinguish between these two alternatives [41]. The winding number W inside C is W ≈ n + − n − where n + and n − are the numbers positively and negatively charged vortices inside C. If vortices were distributed at random with randomly assigned topological charges, the dispersion of typical winding numbers W 2 would be proportional to the square root of total number of defects, n ≡ n + + n − inside the loop. Therefore, it would scale as a square root of the area A C inside C, i.e., However, according to KZM [41], the broken symmetry of the local order parameter -e.g., phases of the superfluid wavefunction -rather than topological charge of defects is distributed at random. Therefore, W is determined by the winding of the phase along the loop C. The random choices of the phase inξ-sized domains along C imply that the accumulated typical winding number should be proportional to C/ξ as long as C >ξ, where C/ξ is the number of domains with the size given by the correlation lengthξ. Therefore, the dispersion of W predicted by KZM should scale as, Furthermore, in the range W 2 1 the absolute value of W has the same scaling limit as W 2 . More precisely, the relation is [41], in which the prefactor π/2 comes from the Gaussian approximation to the distribution of W. The above scalings need to be adjusted if the magnitude of W is smaller than 1. This happens as the radius r of the contour is smaller than the correlation length r <ξ. In this case W is proportional to the probability of finding one vortex inside C. Thus, |W| ≈ p + + p − , where p +/− is the probability to find a positive/negative vortex. In this case W only has three possible values +1, 0 and −1, thus we can deduce |W| = W 2 . Therefore, in the limit of |W| 1 we arrive at (2.8) in which A C is the area surrounded by the contour C.

Dynamics of symmetry breaking and nascent topological defects
Growth of the average absolute value of order parameter |O(t)| from t = 0 (T = T c ) to the final equilibrium state is seen in figure 2(a). The instantaneous dynamical values of |O(t)| remain negligible, i.e. close to what it was in the symmetric vacuum, and, hence, lags behind the instantaneous equilibrium values. For instance for quench with τ Q = 190, its instantaneous value remains negligible until the lag time t L /τ Q ∼ 0.263, and then begins to grow rapidly reaching the approximate equilibrium value at t/τ Q ∼ 0.368. This behavior, with t L larger than but proportional tot, is predicted by KZM, and was reported in the past [12,35].
which belongs to type II superconductors [38].

Number density of topological defects and "freeze-out" time
We count the vortex number density n as the average order parameter saturates to its equilibrium value. Scalings between n and τ Q are exhibited in figure 3(a), with the final equilibrium temperature T f = 0.64T c and the size of the boundary (x, y) = (100, 100). In the slow quench regime (large τ Q ), the scaling relation is fitted as n ≈ (836.0430±1.1844)× τ −0.5126±0.0230 Q , where the uncertainties give standard deviations. The quasi-normal modes (QNMs) analysis in the supplementary material indicates a mean-field theory with ν ≈ 1/2 and z ≈ 2. Thus, the exponent in the scaling between n ∼ τ Q in eq. (2.5) is roughly −1/2, which is in good agreement with the above numerical results. By contrast, in the fast quench regime (small τ Q which is beyond the scope of KZM) the vortex number is approximately constant. This is consistent with the previous results [27,34,35].
The timet ∼ √ τ Q in the symmetry breaking phase is the instant that the system leaves the phase of the quench when the dynamics of the systems cannot keep up with the changes imposed by the quench and enters the adiabatic region [28]. Following [12,35] we define the lag time t L as the time when the order parameter begins to grow rapidly. Lag time t L reflects the "freeze-out" timet [12,35]. Scalings between t L and τ Q are shown in figure 3

Typical winding number W
The winding number W is phase accumulated along the contour C. In the broken symmetry phase it is dominated by the net number of topological defects inside C. According to KZM, random choices of the phase of the order parameter determine the distribution of the topological defects. As a consequence, defect charges are not distributed randomly, but, rather, anticorrelated. This leads to predictions about the distribution of winding numbers as a function of the circumference of C [41].   Figure 4(b) exhibits the relations between W and the average total vortex number n inside the contour C. Because of n ∼ A C /ξ 2 ∝ r 2 /ξ 2 , the doublings of powers for scalings W ∼ n and W ∼ r are explicitly shown by comparing figure 4(b) and figure 4(a). In the limit r <ξ, it is usually only possible to find one vortex (either positive or negative) inside C, thus |W| ≈ W 2 ≈ n inside a small C. This is demonstrated as wellln( |W| ) ≈ ln( W 2 ) ≈ ln( n ) in figure 4(b) when n is small. Therefore, the universal scalings of typical winding number W in our holographic study are consistent with predictions from KZM [16,41,42].

Vortex-vortex correlation function with polarity
Correlation lengthξ can be estimated from the above W. We see that, |W| = W 2 as r <ξ, however, this equality is violated as r ξ . Figure 4(a) shows that |W| departs from W 2 at around 1 < ln(r) < 2, thus we can estimate 2.7 <ξ < 7.3. Following [24,25], one can evaluateξ from the vortex-vortex correlation function G(r) with vortex polarities as well. G(r) is defined as G(r) ≡ n(r)n(0) , with n(r) = +1/ − 1 at the location of a positive/negative vortex, and 0 elsewhere. In practice, G(r) can be evaluated by summing over all charges of vortices (with ± polarities) sitting at the circumference a contour C, whose center is at a positive vortex. Meanwhile, one can also define the net vortex number n c (r) by summing over vortex charges inside the above contour. 3 Figure 5(a) shows that n c = 1 at r = 0, which is obvious from its definition. Away from r = 0, n c decreases to zero, which demonstrates the short-range vortex-antivortex correlations. From the definition of G(r), one can set G(0) = 0. Negative minimum of G(r) in figure 5(b) also reflects the short-range, nearest neighbour vortex-antivortex correlations between vortices, and the position of the minimum isξ. Fitting G(r) ≈ ar 2 × e −br 2 to the theory [43,44], we get a ≈ −0.0412 and b ≈ 0.0703. Therefore, the correlation length can be estimated asξ ≈ 1/ √ b ≈ 3.7, which lies in the range 2.7 <ξ < 7.3 of our previous estimate. There is another length scale -mean vortex separation r av = A/ n , in which A is the area of the system [24]. Thus, r av ≈ 12.9099 as τ Q = 30 from figure 3(a) (with A = 100 × 100, n ≈ 60). Consequently,ξ ≈ 0.2921r av which is comparable to the experimental resultsξ ≈ 0.35r av in [24], where the authors studied the distributions of the magnetic flux quanta from KZM in a 2D superconducting film.

Conclusions
Taking advantage of the AdS/CFT correspondence we have simulated quench-induced symmetry breaking in the transition from the normal to superconducting phase of the strongly coupled holographic field theory. We have observed formation of topological defectsfluxons with quantized fluxes trapped within the vortices, with properties consistent with type II superconductor. Their densities accord with the predictions of KZM, as does their distribution. In particular, they are anticorrelated. This is related to the distribution of the winding numbers of the phase of the condensate. They provide evidence that it is the condensate phase (i.e., post-transition choice of the broken symmetry) that is random, which in turn results in the anticorrelation of the topological charges. We have also observed that the lag time -the instant at which the order parameter begins to grow rapidly -scales as the freezeout timet, central to KZM.

JHEP03(2021)136
where D = ∇ − iA. (Throughout this paper, we work in the units with e = c = = k B = 1.) The ansatz we take is Ψ = Ψ(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. Then the equations of motion (EoM) read, The asymptotic expansions of fields near the boundary z → 0 are (we have set In the numerics we have scaled L = 1. From AdS/CFT correspondence, a t , a i (i = x, y) and Ψ 0 are interpreted as the chemical potential, gauge field velocity and source of scalar operators on the boundary, respectively. Their corresponding conjugate variables can be achieved by varying the renormalized on-shell action S ren. with respect to these source terms. From holographic renormalization [46], we can add the counter terms of the scalar fields S c.t. = d 3 x √ −γΨ * Ψ into the divergent on-shell action, where γ is the determinant of the reduced metric on the z → 0 boundary. In order to get the dynamical gauge fields JHEP03(2021)136 in the boundary, we impose the Neumann boundary conditions for the gauge fields as z → 0 [47,48]. Thus, the surface term S surf. = d 3 x √ −γn µ F µν A ν near the boundary should be added as well in order to have a well-defined variation, where n µ is the normal vector perpendicular to the boundary. Hence, we obtain the finite renormalized on-shell action S ren. . Therefore, the expectation value of the order parameter, O = Ψ 1 , can be obtained by varying S ren. with respect to Ψ 0 . Expanding the z-component of Maxwell equations near boundary, we get ∂ t b t + ∂ i J i = 0, which is exactly a conservation equation of the charge density and current on the boundary, since from the variation of S ren. one can easily deduce that b t = −ρ with ρ the charge density and On the boundary, we set Ψ 0 = 0 in order to have spontaneously broken symmetry of the order parameter. The Neumann boundary conditions for the gauge fields can be imposed from the above conservation equations. Therefore, dynamical gauge fields on the boundary can be computed and lead to the spontaneous formation of magnetic vortices. Moreover, we impose the periodic boundary conditions for all the fields on the spatial boundary along (x, y)-directions. Near the horizon we set A t (z h ) = 0 and the regular finite boundary conditions for other fields.
From the dimension analysis, we know that the temperature of the black hole T has mass dimension one, while the mass dimension of the charge density ρ has mass dimension two. Therefore, T / √ ρ is dimensionless. From holographic superconductor [45], decreasing the temperature is equivalent to increasing the charge density. Therefore, in order to linearize the temperature near the critical point like eq. (2.2), we can actually quench the charge density ρ as ρ(t) = ρ c (1 − t/τ Q ) −2 with ρ c ≈ 4.06, which is the critical charge density from the equilibrium state of the holographic superconductor. For all quench rates in our paper, we quench the temperature from the initial temperature T i = 1.30T c to the final temperature T f = 0.64T c (Except in figure 1, the final temperature is T f = 0.9T c ).
After the quench, we maintain the temperature at T f until the system arrives at the final equilibrium state.
Numerical schemes. We thermalize the system thoroughly before quench in order to make an equilibrium initial state. In order to thermalize the system initially, different from putting the seeds on the boundary in [34,35], we put the random seeds of the fields in the bulk by satisfying the statistical distributions s(t, x i ) = 0 and s(t, x i )s(t , x j ) = hδ(t − t )δ(x i − x j ), with the amplitude h = 10 −3 . 4 System evolves by using the fourth order Runge-Kutta method with time step ∆t = 0.025. In the radial direction z, we use the Chebyshev pseudo-spectral method with 21 grids. Since in the (x, y)-directions, all the fields are periodic, we use the Fourier decomposition along (x, y)-directions with the spatial spacing ∆x = ∆y = 0.25. The codes are indeed very robust to the grids in time and spatial directions. We choose these grids in our paper considering both the consumptions of time and the good quality of the results. Filtering of the high momentum modes are implemented following the "2/3's rule" that the uppermost one third Fourier modes are removed [49].