Quantum Formation of Topological Defects

We consider quantum phase transitions with global symmetry breakings that result in the formation of topological defects. We evaluate the number densities of kinks, vortices, and monopoles that are produced in $d=1,2,3$ spatial dimensions respectively and find that they scale as $t^{-d/2}$ and evolve towards attractor solutions that are independent of the quench timescale. For $d=1$ our results apply in the region of parameters $\lambda \tau/m \ll 1$ where $\lambda$ is the quartic self-interaction of the order parameter, $\tau$ is the quench timescale, and $m$ the mass parameter.


I. INTRODUCTION
The formation of topological defects during a quantum phase transition is a novel process in which the quantum vacuum spontaneously breaks up into classical objects. In a thermal phase transition, the formation of defects is also a transition from a collection of particles above the critical temperature to a collection of a complex of particles (defects) and new excitations at low temperatures. It is no surprise that there has been so much theoretical and experimental  interest in understanding details of defect formation.
The number density of defects formed during a phase transition is sensitive to the rates at which external parameters are changed to pass through the phase transition. The leading theoretical framework for estimating the number density of defects is the "Kibble-Zurek" analysis [1][2][3][4][5][6]. Numerical simulations have further strengthened the model [24][25][26][27][28][29][30]. However predictions of the Kibble-Zurek model have not yet gained universal confirmation, with most experiments in systems involving 4 He, liquid crystals, superconductors, superfluids in agreement [9][10][11][12][13][14][15][16][17][18][19][20][21][22][23] and others in disagreement [11,18,31] with the predictions. In particular, the appearance of vortices in 4 He was claimed in [13] but was retracted in [11] since it was found that the vortices in the former case was an externally induced artifact. Overall, the analysis of the phenomenon of defect formation in various systems is an ongoing field of research and has broad implications.
In the present work we follow our analysis of [32] and solve for the number density of defects (kinks, vortices and monopoles) formed during a quantum phase transition. The analysis is rigorous and without recourse to approximation but the quantum field theory models we consider are "free", the only interaction being with external parameters that drive the phase transition. These models provide us with zeroth order solvable problems in different dimensions that we fully analyze. Even with these minimal interactions, the analysis is highly nontrivial and in part has to be done numerically. We discuss how other interactions may be included in the analysis using perturbation theory and under what conditions we expect the zeroth order approximation to be accurate.
We are generally interested in Poincaré invariant fieldtheoretic models in d + 1 spacetime dimensions, featuring an internal (global) O(d) symmetry which is spontaneously broken during a quantum phase transition. In particular we will be considering d real scalar fields Φ 1 , . . . , Φ d assembled in an O(d)-multiplet Φ ≡ (Φ 1 , . . . , Φ d ) T whose dynamics are given by the Lagrangian density 1 Here the potential V β is O(d)-invariant and depends on a (possibly time-dependent) external parameter β. We assume that V β is such that the vacuum manifold is O(d)symmetric for β < 0 and O(d − 1)-symmetric for β > 0. In other words, as the parameter β increases from negative to positive values, the system transitions from a higher symmetry phase to a lower symmetry one, and the average vacuum field configuration starts exhibiting topological defects. These defects then annihilate with one another and eventually disappear. It is precisely this dynamics of formation and annihilation of topological defects that we are concerned with in this paper. In fact, our main purpose will be to determine the number density of topological defects as a function of time and its dependence on the external parameter β, using a combination of analytical and numerical methods. For concreteness we will take β to be the (timedependent) mass squared of the field, so that where 1 In this paper we use a mostly plus signature for the Minkowski metric and natural units, = c = 1. and λ, m, τ are positive parameters. In particular, the quench parameter τ is a time scale quantifying the rate of change of the potential during the phase transition. It is clear that for t −τ , the vacuum manifold reduces to the null field configuration Φ = 0 and is therefore O(d)-symmetric, while for t τ it includes all field configurations on the O(d − 1)-symmetric hypersphere given by In Fig. 1 we sketch the potential (2) at a few different times.
It is well known that these models have topological defects -kinks (d = 1) in one spatial dimension, vortices (d = 2) in two spatial dimensions, and monopoles (d = 3) in three spatial dimensions [33]. In each of these cases the vacuum manifold described by (4) has non-trivial topology: for d = 1 it is 2 points, for d = 2 it is a circle, and for d = 3 it is a two-sphere. The defect locations are described by zeros of Φ even in the symmetry broken phase. The zeros are trapped due to the non-trivial topology of the vacuum manifold. We realize that the topology persists even if we set λ = 0 and the problem of defect formation simplifies. Then the λ = 0 problem can be thought of as the zeroth order problem. We discuss the λ = 0 problem for d = 1 in greater detail in Sec. V where we find that λ dependent corrections are small if λτ /m 1. The overall strategy will be to regulate both the IR and UV behaviors of the field theory by working in a finite box of size L d (with periodic boundary conditions) and discretize space on a N d point lattice (with lattice spacing a = L/N ). Then we can determine an exact expression for the field probability density functional as a function of a finite number of quantities that can be computed numerically. We then find the average expectation value of a judiciously constructed quantum operator that counts the number density of zeros of the field multiplet Φ in the limit of the finite resolution imposed by the lattice. We finally take both the continuum limit N → ∞, a → 0, and the infinite volume limit, L → ∞ (in this exact order), to recover the full field theory result. Up to spurious zeros due to vacuum fluctuations that can consistently be discarded, this accurately gives the number density of topological defects. The case of a sudden phase transition (τ = 0) can be treated analytically but the general case will be treated numerically.
The paper is structured as follows. In Sections II and III we fully describe the average dynamics of kink (d = 1) and vortex (d = 2) condensation respectively. In Section IV we extend these results to 3 and higher dimensions. In Sec. V we discuss how the previous results constitute only the zeroth order approximation in a perturbative expansion in λ and estimate the next-order corrections. We end in Section VI by emphasizing the importance of our results and contrasting them with previous work done on the subject.

II. ONE DIMENSION: KINKS
One of the challenges in finding the number density of kinks is to first define a kink in the quantum field theory given by (1) with d = 1, where we denote the singlecomponent scalar field Φ by φ. This can be done using the Mandelstam "kink operator" [34], which is a twocomponent fermionic operator,χ. A key property ofχ is that it satisfies the equal time commutation relations, where η is a real number. If |s is an eigenstate ofφ(t, y) such thatφ|s = 0 (for all y), then we find that the state |s ≡χ(t, x)|s satisfieŝ Hence the operatorχ has created a step in the value of φ at x by an amount η. If φ = 0 and φ = η are two possible vacuum expectation values of φ,χ would have created a kink that interpolates between two vacua. The number density of χ quanta would then correspond to the number density of kinks. Unfortunately the relation between χ and φ is quite complicated -χ involves exponentials of φ andφ and other quantum field theory subtleties -and we do not have a clear way to utilize the Mandelstam operator. Instead we find it useful to work entirely with the φ field, simply defining the kink to be a jump in the value of φ as further discussed in Sec. II B. Our definition of the kink operator is also helpful in the case of vortices and monopoles for d = 2, 3 as these objects correspond to intersections of domain walls i.e. kinks extended to higher dimensions.

A. Setup and quantization
We start by treating the d = 1 case in detail. The relevant Lagrangian density for the real scalar field φ is thus Clearly, for t < 0 the model has a unique vacuum φ = 0 while for t > 0 it has two degenerate vacua at φ = ±m/ √ λ corresponding to the two minima of the doublewell potential. It is well-known that in the t τ limit (where m 2 (t) ≈ −m 2 ), there exist static classical kink and anti-kink solutions given by These solutions are non-perturbative and topologically non-trivial: they interpolate between the two vacua over a spatial scale ∼ 1/m. Of course, Poincaré invariance allows the construction of displaced or even "dynamical" kinks from the above solutions but, whatever the frame, they will always be characterized by their topological charge In fact, a kink always has positive topological charge since the field undergoes a negative to positive sign change, while an anti-kink has the exact opposite property.
Multi-kink and anti-kink solutions can be constructed as well, but these will not be static anymore since the kinks and anti-kinks will attract each other and they will eventually annihilate. If separations are large and the different kinks and anti-kinks are initially at rest, such configurations will however be approximately static. Even though the topological charge of such field configurations does not inform us about the number of kinks or antikinks involved (since the topological charge is a binary valued quantity), one can however in principle recognize the presence of individual kinks and anti-kinks in a given field configuration by focusing on the points where the field changes sign: a negative-to-positive sign change will be identified as a kink while a positive-to-negative one will be identified as an anti-kink. Of course this is only part of the picture because not every sign change should be counted as a kink or anti-kink especially if it occurs on time and distance scales shorter than the characteristic width of 1/m. We will discuss this subtlety in Sec. II B.
We are interested in the production of kinks during a quantum phase transition and in particular in how their average number density scales with time. As we have discussed in Sec. I, we will first be analyzing the λ = 0 case and the Lagrangian density we will work with will thus be We now need to quantize this model. We start by assuming that the volume (or length since d = 1) of space is finite of size L and that the field obeys periodic boundary conditions. (We can alternatively think of space as a circle of length L.) We then discretize space on a lattice consisting of N points separated by a distance a = L/N . At each lattice point x j ≡ ja, we define the discretized field φ j ≡ φ(x j ) and the full Lagrangian of the discretized theory reads where we have assembled the discretized fields in a column vector φ ≡ (φ 1 , . . . , φ N ) T and the matrix Ω 2 is defined by [ Introducing the canonically conjugate momentum fields, π j ≡ aφ j , and assembling them in a column vector π ≡ (π 1 , . . . , π N ) T , we can promote both the φ j s and π j s to operators satisfying canonical commutation relations [φ j ,π l ] = iδ jl . The quantum Hamiltonian of the discretized theory [35] then readŝ where hats denote operator valued quantities. It is apparent from (13) that the discretized theory describes the quantum dynamics of a set of N quadratically coupled, simple harmonic oscillators. We are interested in how the (unique) quantum vacuum before the phase transition (at a time t 0 −τ when the potential is upright and m 2 (t 0 ) ≈ m 2 ) is destabilized by the quench and evolves into a more complicated state featuring dynamical kinks and anti-kinks. To understand the dynamics of this process we need to solve the functional Schrödinger equation associated with (13), where the wave functional Ψ[φ 1 , . . . , φ N ; t] is such that |Ψ| 2 gives the probability density of a given field configuration at time t, and the Laplacian operator is defined by One can easily check that the wave functional for the vacuum state at t = t 0 is where and fractional powers of the positive-definite matrix Ω 2 (t 0 ) are unambiguously defined in the standard way. For instance Ω where O is the orthogonal matrix diagonalizing Ω 2 , and λ j are the (positive) eigenvalues of Ω 2 . Given this initial condition, the solution for the wave functional at time t will be given by and and we can write Indeed, using (20), (21) and (22), it is easy to check that this expression yields a symmetric matrix sinceŻZ −1 − (ŻZ −1 ) T is a conserved quantity which vanishes at time t 0 . We can now write the probability density functional as To simplify this expression we first use the fact that where K ≡ ZZ † is a real positive definite symmetric matrix; indeed, using (20), (21) and (22), it is easy to check that ZZ † − Z * Z T is a conserved quantity which vanishes at time t 0 . Next, we make use of another conserved quantity which can also be verified via (20), (21) and (22), to show that Finally, plugging (25) and (27) into (24) yields a simplified (and manifestly normalized) expression for the probability density functional This expression (along with (20), (21) and (22)) contains all the information that we will need in order to determine the average number density of kinks in the lattice. Note that K is a time-dependent matrix, whose timedependence is given by that of the matrix Z.
Before going any further, we mention a separate interpretation of the matrix Z. Working in Heisenberg picture with respect to time t 0 , we can define creation and annihilation operators at time t 0 bŷ Notice that we have used column vector notation here but that the dagger refers to the adjoint operation on the Hilbert space only: it does not turn column vectors into row vectors. Then we can expand the Heisenberg picture discretized field operators at time t as followŝ Eqs. (20), (21) and (22) ensure that the Heisenberg equations as well as the proper initial conditions at t 0 are verified. Now it is easy to see that the matrix K is simply the covariance matrix of the discretized field since, using (31), Here the Heisenberg picture vacuum |0 is timeindependent and defined by the wave functional (16). In principle we now have all the ingredients needed to discuss the quantum production of kinks during the phase transition. Indeed equation (28) along with the N 2 complex linear ordinary differential equations (20) fully determine the quantum dynamics of the field configuration. However, it turns out that not all components of the matrix Z are relevant and we can reduce the number of differential equations that need to be solved. It can be shown that the matrix Z is circulant [36] i.e., its matrix elements Z jl only depend on j − l (mod N ). We can therefore diagonalize it via the discrete Fourier transform: This allows us to recast (20), (21) and (22) in terms of the complex mode functions c n (t) thus obtaining and Rewriting the dynamical equations in terms of mode coefficients provides an enormous computational gain: we now only have to solve N equations instead of N 2 . Additionally, as we will shortly see, mode coefficients are particularly well suited to discussing problems related to the N → ∞ limit and divergences related to vacuum fluctuations of the quantum field. We can achieve further simplification by writing the mode functions in trigonometric form where ρ n and θ n are respectively the modulus and argument of the complex number c n , and making use of the conserved quantity (26), which in this representation takes the form of a conserved angular momentum, Then (34) reduces to a set of N real (but non-linear) ordinary differential equations, with initial conditions Even though working in terms of modes is computationally advantageous, kinks are configurations (field zeros) in physical space. Thus we have to straddle the two descriptions as in the following Sec. II B.

B. Average kink number density
We are now in a position to tackle the problem of kink production during the phase transition. As mentioned earlier, since kinks and anti-kinks occur at zeros of the field configuration we first introduce a quantum operator n Z that gives the number density of zeros in a given field configuration: More precisely, such an operator is sensitive to the number of sign changes that occur between adjacent points on the lattice. We should stress that this is only accurate up to the finite resolution given by the lattice spacing a. It may in fact undercount the number of zeros of the actual continuous field configuration (if there are multiple sign changes within a lattice spacing). We expect however that, as N becomes large enough, this operator will become more and more accurate. This assumption is reasonable as long as we can find a way to disregard high frequency noise-like fluctuations due to the quantum vacuum thus only counting "true" kinks and anti-kinks. We now calculate the vacuum expectation value of this operator or, in Heisenberg picture Given that we know the probability density functional explicitly for the Schrödinger picture time-dependent state we can write wherë sgn Äφ Introducing the permutation (shift) matrix and performing the change of variables φ → P 1−j φ, we can rewrite (45) as sgn Äφ As mentioned earlier, Z is a circulant matrix and, consequently, it has to be polynomial in P . Therefore, the matrix K = ZZ † is also circulant and K −1 is seen to commute with P . This implies that, and the average number density of zeros simply reduces to Let it be mentioned here that the circulant property of the covariance matrix K is the mathematical counterpart of the fact that the system has translational invariance (which is maintained at a discretized level by our choice of periodic boundary conditions). In other words, it is a consequence of the fact that two-point correlation functions φ(x)φ(y) only depend on the relative position |x − y|.
We now need to evaluate (48) more explicitly. We start by writing where the sum runs over the four quadrants in the (φ 1 , φ 2 ) plane (denoted by Roman numerals). We then decompose K −1 into suitably sized blocks, where A and C are real symmetric matrices of respective We also assume that C is invertible, which will be true generically. This allows us to rewrite the bilinear in the exponent in (50) as Using we can perform the Gaussian integral over φ 3 , . . . , φ N and obtain sgn Äφ is the so-called Schur complement of C. The left-over two-dimensional quadrant integrals can also be carried out. For the first quadrant, for example, sgn(φ 1 φ 2 ) = +1 and we can write where in going from the first to the second line we used the change of variables φ 1 → sφ 2 .
The integrals over the remaining three quadrants are readily obtained from I I as follows. To begin with, the change of variables φ 1 → −φ 1 and φ 2 → −φ 2 makes it clear that I III = I I , and I II = I IV . Furthermore, notice that the change of variables φ 1 → −φ 1 (leaving φ 2 unchanged) on I II has the same effect (up to an overall sign) as changing A 12 into −A 12 in (57). We thus obtain and all the four integrals I Q appearing in (54) are accounted for. We can achieve further simplification by taking advantage of the properties of the matrix A . In particular since we have and (54) collapses tö But we can go even further. Indeed inverting (59), shows that A −1 coincides with the upper-left 2 × 2 block of the matrix K. More explicitly we can write, where, using (33) and the reality of K, Thus and (61) becomes Finally we obtain the average number density of zeros Recall however that we are interested in the average number density of kinks which may differ from the number density of zeros as given in (68) because the latter includes zeros due to vacuum fluctuations of the quantum field. The difference between the two quantities is most clear long before the phase transition, where the field is in its unique vacuum and its expectation value vanishes everywhere on the lattice. However the field fluctuates about zero and there is a non-zero average number density of zeros. This is to be contrasted with the average number density of kinks which should be exactly zero before the phase transition. Moreover the average number density of zeros is expected to be highly sensitive to the number of lattice points N since the finer the resolution, the more zeros can be identified. This is again different for the average number density of kinks which are supposed to be extended objects whose separation is set by the correlation length of the field fluctuations. We therefore need a systematic procedure to eliminate the spurious zeros from the result in (68). One way is to restrict the sums in (64), (65) to those modes c n (t) that are not oscillating [37], in other words to indices n verifying It is indeed the presence of such unstable modes that allows for the production of the non-perturbative kink and anti-kink solutions. Then the formula for the average number density of kinks, n K , is obtained by restricting the modes that enter (68), giving us where nowᾱ These equations only apply for t ≥ 0 when the modes start to become unstable. For t < 0, there are only fluctuating modes and we set n K = 0. We will discuss the difference between n Z and n K in Sec. II D.
After the phase transition and as long as the lattice spacing a is small enough, a < 2/ |m 2 (t)| for all times t > 0, we can introduce n c (t), the time-dependent critical value of n that separates unstable modes from modes that oscillate, where denotes the integer part function. Then n c (t) < N/2 and (71), (72) can be rewritten in a more explicit wayᾱ Here we have identified c −n with c N −n for concision, and exploited the symmetry c N −n = c n (valid for 1 ≤ n ≤ N − 1) which can be checked directly via (34), (35), (36).
Since the ratioβ/ᾱ belongs to the interval [0, 1] one can also rewrite (70) as Before diving into analytical and numerical estimates of n K we need to discuss the continuum and infinite volume limits of our discretized theory. We start with the continuum limit. Keeping L fixed, and noticing that, for all N , n c (t) ≤ mL/4, we can safely take the N → ∞ limit in expressions involving n/N . In particular and Then the expression for the ratioβ/ᾱ reads Now, it is clear that this expression is of the form 1 − 2x 2 with x ∈ [0, 1] and therefore we may use the identity to simplify (76) and obtain This is the expression of the continuum limit (N → ∞) average number density of kinks. The main property of this expression is that is does not depend on N anymore. Indeed, although the system's dynamics is governed by an infinite number of mode functions, only a finite number appears in the formula; it is only those modes with n ≤ mL/2π that trigger the instabilities required for the production of kinks. This means that the result is stable in the UV limit and does not depend on the resolution of our discretization. Physically, the contribution of the vacuum fluctuations of the quantum field has been discarded.
Let us now end this section by discussing the infinite volume (or length since we are working in one spatial dimension) limit L → ∞. This is readily done by noticing that the finite size of the spatial dimension is responsible for the discreteness of the wave vectors corresponding to different modes. As L increases, however, these wave vectors become more and more numerous and densely packed until they form a continuum spanning the entire interval [−m, m]. At this point, it is convenient to switch notations and index any relevant quantities by k n instead of just n. Then and The average kink number density can therefore be rewritten In the L → ∞ limit the sums over k n become integrals over k, so that Here we have tacitly introduced the infinite volume mode functions c k verifying and used their k → −k symmetry properties. We now have all the tools required to perform simple analytical estimates of the average kink number density.

C. Analytical estimate
In the limit of a sudden phase transition (τ = 0), we can solve (34) exactly since m 2 (t) = −m 2 Θ(t) (where Θ is the Heaviside function). In fact one may then choose the initial time to be t 0 = 0 − and solve the differential equationsc with initial conditions Since the time dependence of the frequency has disappeared, the above differential equations can be solved analytically. This yields the unstable mode functions c n (t) involved in the formula for the average number density of kinks (70) i.e. those verifying |n| ≤ N sin −1 (ma/2)/π. More precisely we have where κ n = » m 2 − 4 a 2 sin 2 πn N . Taking first the continuum limit N → ∞ while keeping L fixed we obtain, for |n| ≤ mL/2π, where we labelled the mode functions by k n = 2πn/L as in the previous section. In the L → ∞ limit, the discrete variable k n becomes continuous and we can write an analytical formula for the average kink number density as in (86): With this expression in hand we can immediately estimate a few important quantities. First of all, we can predict the late time behavior of the average kink number density. Indeed for large t the integrals simplify considerably and it is easy to see that they are dominated by values of k m. We can then estimate (93) to be Using Eq. (93), one can also estimate the maximum number density of kinks that are produced after the phase transition. In fact, taking a time derivative of (93), it is easy to convince oneself that this maximum occurs at t = 0 + , in other words, immediately after the phase transition. Moreover its value can be computed exactly to be Both the power law for the asymptotic behavior of the average kink number density and the maximum number of kinks value will be numerically confirmed in the following subsection. Our analytic results agree with previous work on sudden phase transitions in thermal quenches studied in [38,39] using different techniques.

D. Numerical results
We now discuss our numerical results for the time evolution of n K for different values of the quench parameter τ . In principle this involves solving the complex differential equations (34) with initial conditions (35) and (36), for the unstable mode functions c n (t) -those with |n| ≤ n c (t). We can then directly evaluate the average number density of kinks using (70). However, since this formula only involves |c n (t)| = ρ n (t), considerable computational gain can be achieved by instead solving the real differential equations (39) with initial conditions (40), (41).
It turns out that this system of ordinary differential equations presents a major computational difficulty caused by the fact that ρ n (t) grows exponentially for |n| ≤ n c (t). Therefore the numerical evolution is limited to short time periods after the phase transition beyond which the numbers involved become extremely large and results cannot be trusted. One way to get around this problem is to factor out the exponential growth, i.e. the zero mode ρ 0 (t) = ρ N (t), from the other modes and evolve it separately. So we write, for n = 1, . . . , N − 1. With this redefinition it can be shown that the differential equation (39) now becomes, and its corresponding initial conditions are given by, Recall here that and ω 2 . Furthermore, one can also efficiently solve for ρ 0 (t) by introducing the auxiliary function q(t) = ln ρ 0 (t), verifying with initial conditions, By going to the q(t) variable we avoid the exponential growth of ρ 0 (t). Thus, both the differential equation for r n (t) (97) and its corresponding initial conditions can be rewritten in terms of this auxiliary function: with initial conditions, In summary, the numerically efficient way to study the dynamics of kink formation in our model, is to solve (101) and (103) with respective initial conditions (102) and (104). The computational problem we had is indeed resolved since we managed to eliminate the exponential growth of ρ n (t) by suitable function redefinitions. The numerics can now be trusted for much longer periods of time. In our numerical work we work in units where m = 1 and pick t 0 = −200. To get accurate results we choose large L and N . Most of our results are for L = 6400 and N = 12800. The evolution of the average number density of kinks n K for different quench time scales τ is shown in Fig. 2. The different curves exhibit the same qualitative behavior: immediately after the phase transition (t = 0) the average number density of kinks increases from 0 to a maximum value (n K ) max within a time t max , and this is followed by a gradual decrease that asymp-totically converges to a power law. Physically this corresponds to the production of a random distribution of kinks and anti-kinks during the phase transition, followed by their mutual annihilation over time. Noticeably, the asymptotic behavior of the average kink number density is independent of the quench time scale: at late times the plots for different values of τ converge to the same function that falls off as t −1/2 . (We have also cross-checked this result by calculating the correlation length ξ(t) of field fluctuations and plotting 1/ξ ∼ n K as a function of time.) This scaling law agrees with the analytical estimate of Eq. (94) and shows that the τ = 0 solution is a universal attractor. To analyze the rate at which the kink densities for different values of τ converge, we plot ∆n K (t, τ 1 , τ 2 ) ≡ n K (t, τ 1 ) − n K (t, τ 2 ) versus t in Fig. 3. We observe that at late times these differences fall off as t −3/2 . We can therefore conclude that where C K ≈ 0.22 is a constant of proportionality which is independent of the quench time scale τ . This agrees well with the analytical estimate found in Eq. (94): 1/(π √ 2) ≈ 0.225. We can explicitly check, as shown in Fig. 4, that our results are independent of both L and N as long as they are sufficiently large and a = L/N is sufficiently small. In Fig. 5 we have also plotted n Z and n K for different values of N . Although the curves depend on N (or are UV sensitive) near the phase transition, the late time behaviors are universal. This is to be expected since unstable modes grow exponentially and dominate the sums in (74) and (75). Thus our technique of restricting the mode sums to differentiate between field zeros and kinks is reasonable and gets rid of the artifacts arising due to finite N . The plots of (n K ) max versus τ , and t max versus τ , are shown in Fig. 6 and Fig. 7 respectively. From these we note that the faster the phase transition (smaller quench time τ ), the more kinks and anti-kinks are produced and the faster their maximum number density is attained. In Fig. 6 we see that the maximum density of kinks (n K ) max flattens, i.e. it becomes a constant as quench time scales approach zero. The value of (n K ) max for which this happens is seen to be approximately 0.175. This agrees remarkably well with the analytical estimate in Eq. (95).

III. TWO DIMENSIONS: VORTICES
The analysis done in Section II can be generalized to the d = 2 case. We will be considering a two-dimensional complex scalar field Φ whose dynamics are described by the Lagrangian density This theory is known to possess solitonic solutions called vortices, characterized by a topological charge known as the winding number. Assuming a vortex field configuration Φ(x, y) = r(x, y)e iθ(x,y) = φ(x, y)+iψ(x, y) centered at a point (x 0 , y 0 ), the winding number is given by where C is any closed loop around (x 0 , y 0 ). Generically a non-zero winding number along a closed loop implies the existence of a vortex configuration and the vanishing of the field somewhere within the bounded region. Therefore, as in the case of kinks, vortices are to be found among zeros of Φ.
To study the production of vortices during the quantum phase transition we will thus do a similar analysis to the one we did for kinks. We start by setting λ to zero and express the Lagrangian density in terms of the two real scalar fields φ and ψ, respectively defined as the real and imaginary part of the complex field Φ: This is a model for two non-interacting real scalar fields in two spatial dimensions. In order to apply the methods outlined in Section II, we need to discretize this model. We first compactify both spatial dimensions by assuming periodic boundary conditions, φ(x+L, y) = φ(x, y +L) = φ(x, y) (and similarly for ψ). Space is thus seen to be a 2-torus of area L 2 . We then discretize it on a regular square lattice consisting of N 2 points separated by a distance a = L/N along both the x and y directions. Now for each lattice point (x j , y l ) ≡ (ja, la) we can define the discretized fields φ jl ≡ φ(x j , y l ) and ψ jl ≡ ψ(x j , y l ). Writing the discretized Lagrangian and quantizing it can be done analogously to the one-dimensional case, with the understanding that any vectors and matrices are now N 2 and N 2 × N 2 dimensional respectively. For example, the vector of discretized field values of φ is given by More generally, any N 2 ×N 2 matrix A will be represented by a two-dimensional array of matrix elements A ij,kl arranged in the following way:  (20), (21) and (22) as long as the matrix elements of Ω 2 are given by otherwise .
(110) It is then easy to write the probability density functional as in Eq. (28), where the matrix K is still related to Z via K = ZZ † . We can be even more explicit by realizing that the matrix Z(t) is once again real and circulant, i.e., the matrix elements of Z, Z pq,rs depend only on p − r (mod N ) and q − s (mod N ). We can therefore again diagonalize Z using the discrete Fourier transform: Using equations (20), (21) and (22), the complex mode functions c n,n (t) verifÿ and Note that c n,n = c n ,n which immediately implies that Z pq,rs = Z qp,sr and again we assume the initial time t 0 to be such that t 0 −τ . This follows from the rotational symmetry of the system.

A. Average vortex number density
To find the vortex number density, we first need a quantum operator that counts the number of zeros n Z of the complex field Φ (as in Sec. II B), or in other words, coincident zeroes of both the fields φ and ψ. Since space is discretized, such an operator necessarily yields a coarsegrained estimate of the actual number of zeros of a given field configuration. As the number of lattice points N 2 increases so does the operator's resolution: while certain "zeros" cease to be counted, new ones are revealed. In the limit where N → ∞ we expect divergences, just as in the kink case, and we will return to this point later on.
We think of the vortex as the intersection of a domain wall of φ -for our purposes, a domain wall is a curve on which φ = 0 -with a domain wall of ψ. Then, as shown in Fig. 8, there could be a situation where a φ domain wall enters a plaquette through one edge and leaves through the opposite edge, while a ψ domain wall passes through the plaquette in the orthogonal direction. Then the two domain walls must intersect, leading to coincident zeros that correspond to a vortex within that plaquette. Other possibilities include the case where the φ wall enters the plaquette from the lower edge but leaves from the right edge in Fig. 8 while the ψ domain wall goes through as shown or bends to exit from the top edge. It is ambiguous whether a coincident zero exists in these other cases but the ambiguity is minimized as the lattice resolution is increased (N → ∞). Hence we can count zeros of Φ in the large N limit by counting the plaquettes in which φ and ψ domain walls enter across orthogonal edges.
Then, motivated by the discussion in Sec. II B, we can define the number density of zeros of Φ bŷ We can now write down the vacuum expectation value of the operatorn Z : n Z ≡ 0|n Z (t)|0 . Using the fact that the fields φ and ψ are independent and that, consequently, the probability density functional factorizes as in Eq. (111), we first notice thaẗ sgn Äφ ijφi,j+1 Then, using the fact that the matrix K −1 is circulant and, moreover, symmetric under interchange of its first (or last) two indices -properties that are inherited from Z, we can establish thaẗ sgn Äφ Physically, this set of equalities is a manifestation of the translational and rotational invariance of the system. It is also clear that, φ being a dummy variable in the integral of Eq. (117), These properties thus allow us to write the average number of zeros of the field in a very simple form: From this point on, the computation of the average number of zeros follows along the same lines as in Sec. II, and we obtain where α and β are now defined as β ≡ K 11,12 = N n,n =1 |c n,n | 2 cos(2πn /N ) . (123) Here, once again, we have used the reality of the matrix Z.
Eq. (121) gives us the number density of field zeros but we are interested in counting the number density of vortices. We have already discussed how quantum fluctuations can induce a non-zero number density of zeros of the field even in the absence of spontaneous symmetry breaking. We thus need to eliminate such spurious zeros by restricting the sums in (122) and (123) to the mode functions c n,n (t) that are non-oscillating. In this case we include the modes corresponding to n and n verifying, (124) The average number density of vortices formed after the phase transition is finally given by, where,ᾱ β ≡ ω (n,n ) 2 ≤0 |c n,n | 2 cos(2πn /N ) .
Similarly to the case of kinks (see discussion in Sec. II B), this result only makes sense after the phase transition; it is ill-defined before. As might be intuitively expected, the average number density of vortices is obtained, up to a combinatorics factor due to the φ ↔ ψ symmetry, by squaring the average number density of kinks. In the next subsections we will see that this intuition is supported by both analytical and numerical estimates of the asymptotic dynamics of the problem. Analogously to Sec. II B, Eq. (125) can be further simplified by first taking the continuum limit N → ∞ (at fixed volume L) to obtain where the sums run over pairs of integers (n, n ) ∈ Z 2 verifying ω (n,n ) 2 and it is understood that c −n,n ≡ c N −n,n , c n,−n ≡ c n,N −n for 0 ≤ n, n ≤ N − 1. Relabelling the mode functions by the discrete two-dimensional wave vector and taking the large L limit, Eq. (128) can be recast as Here we have once again introduced the infinite volume mode functions c k -labelled by a continuum of twodimensional wave vectors k ≡ (k x , k y ) -verifying and defined k ≡ | k| and k c (t) = |m 2 (t)| as in Sec. II B.
Noticing that c k only depends 2 on k and going to polar coordinates, we can turn the double integrals in (131) into single integrals to finally obtain the continuum, infinite volume limit of the average vortex number density: (133)

B. Analytical estimate
Just as we did in the case of kinks in Sec. II C, we can also compute the average number density of vortices at late times in the limit of a sudden phase transition (τ = 0). This can be achieved once again by exactly solving the differential equations for the mode coefficients c n,n (t). As we saw in Sec. III these differential equations are as follows, with initial conditions The solution to these equations can be obtained analytically. In fact they look very similar to the ones we obtained in the kinks case but now involve two indices instead of just one. This gives the unstable mode functions c n,n (t) involved in the formula for the average number density of vortices: cosh (κ n,n t) where Now, taking first the continuum limit N → ∞ while keeping L fixed we obtain, for n 2 + n 2 ≤ (mL/2π) 2 , where we have relabelled the mode functions by k n,n = (k (n) x , k (n ) y ) and recall that k = 2πn /L. In the limit L → ∞, the discrete variables k n,n become continuous and, as in Eq. (133), we can write an analytical formula for the average number density of vortices: Using this equation, we can once again estimate the late time behavior of the average number of vortices. In the limit, k, k m, we have (141) As mentioned below Eq. (127), the vortex number density is obtained by squaring the kink number density and multiplying by the combinatorial factor 2! due to the exchange symmetry φ ↔ ψ.
Furthermore, like in the case of kinks, the maximum number density of vortices can be estimated using Eq. (140). This maximum is reached immediately after the phase transition, at time t = 0 + and is found to be (142)

C. Numerical results
We use numerical techniques to solve (113) and then calculate the average vortex number density using (125). For reasons discussed earlier, the parameters L and N that we choose for our numerical simulations need to be sufficiently large to accurately describe the continuum infinite volume limit. We choose, L = 2000 and N = 4000. As in the case of kinks, the results are insensitive to the UV and IR cutoffs. In practice, because of the order N 2 computational complexity of the problem and the exponential growth of the magnitudes of mode functions, we directly solve for ρ n,n = |c n,n | and factor out the zero mode to improve the numerical accuracy (see Sec. II D for details). In Fig. 9 we show the average vortex number density for different quench parameters τ as a function of time. As in the kink case, the plots of n V vs. t for different τ converge to the same function and decay as t −1 as we expect from the analytical estimate in (141). The result also agrees with the intuition that a vortex corresponds to the intersection of two independent domain walls. Fig. 9 also shows that immediately after the phase transition, n V increases from zero to some maximum value (n V ) max in a time t max . As time goes on n V starts to decay. At very early times, after the phase transition, randomly distributed vortices of positive and negative winding number are produced, but then the system starts relaxing, the vortices-antivortices start annihilating, and the dynamics reaches its scaling regime. We can also plot the differences of vortex number densities for different values of τ as we did in the case of kinks: ∆n V (t, τ 1 , τ 2 ) ≡ n V (t, τ 1 ) − n V (t, τ 2 ). This is shown in Fig. 10 which shows that ∆n V (t, τ 1 , τ 2 ) decays as t −2 at late times. We thus deduce the asymptotic form, where, C V is some constant of proportionality which is independent of the quench time scale τ . Numerically, we find C V ≈ 0.092. This is again in reasonable agreement with the value we calculated analytically for a sudden phase transition (τ = 0) in Eq. (141), more precisely 1/π 2 ≈ 0.101. The plots of (n V ) max versus τ , and t max versus τ are shown in Fig. 11 and Fig. 12 respectively. The intuitive understanding that a faster phase transition (smaller quench time scale τ ) leads to greater and more rapid vortex production is confirmed by these plots. Moreover, from Fig. 11 we see that the maximum number density of vortices (n V ) max flattens as the quench time scale τ approaches zero. This happens for a value (n V ) max ≈ 0.0483 which is once again in good agreement with our analytical result in Eq. (142).
As a final remark, comparing Fig. 6 to Fig. 11 shows us right away that for the same quench time-scales τ , the maximum vortex number density (n V ) max is much lower than the maximum kink number density (n K ) max . For example, in the limiting case of τ → 0, (n V ) max ≈ 0.050 while (n K ) max ≈ 0.175. This is again to be expected since the formation of a vortex requires the simultaneous vanishing of two fields, which is less probable than the vanishing of a single field necessary for the formation of a kink in one dimension.

IV. HIGHER DIMENSIONS: MONOPOLES
Having worked out the details of the d = 1 and d = 2 cases, it is easy to see that the methods described in the previous sections directly generalize to higher dimensions. Without going into the details of a rigorous proof, the average number density of zero-dimensional topological defects n D formed in the d dimensional field theory discussed in Sec. I is given by for late times. The factor of d! arises because of permutation symmetry. To get a monopole in d dimensions we need coincident zeros of d fields in a cell of the lattice. As in Sec. III, the point Φ = 0 corresponds to the intersection of d orthogonal domain walls. The d! permutations of the wall positions preserves the Φ = 0 point which leads to the d! prefactor in (144). In Fig. 13 we show numerical results for d = 3 for the monopole number density as a function of time, obtaining the first term on the right-hand side of (144). In Fig. 14 we provide evidence for the second term on the righthand side of (144).

V. EFFECT OF SELF-INTERACTIONS
A key question is to understand the range of parameters for which our results are a good approximation even when λ = 0. We will address this in the context of the model in one spatial dimension given in (7). We check for self-consistency of our solution and examine the conditions under which it breaks down.
Our solution for the wavefunction is a Gaussian at all times and so φ 4 = 3 φ 2 2 . With λ = 0, the evolution of the wavefunctional, Ψ[φ, t], will depend on λ. As long as Ψ can be approximated by a Gaussian centered at φ = 0 we can use the Hartree approximation (e.g. [40]) to write λφ 4 as 3λ φ 2 φ 2 . Taking into account mass renormalization at lowest order in λ we obtain an effective mass squared m eff 2 (t), m eff 2 (t) = m 2 (t) + where the "in" subscript refers to evaluation at the initial time (t → −∞). The mass counterterm 3λ φ 2 in /2 is chosen such that the effective mass equals m at the initial time. Therefore, in the Hartree approximation, the effects of interactions are negligible if the λ dependent corrections to m 2 are small and m eff 2 (t) ≈ m 2 (t), or, 3λ The condition in (146) will fail in two circumstances. First, around the time of the phase transition, t ∼ 0, when m 2 ∼ −m 2 t/τ (see (3)) is very small; second, at late times, when φ 2 grows large. We can make these statements more precise by noticing that (146) is strongly violated whenever the function becomes negative. It turns out that, generically, f λ,τ has three zeros that we denote as t 1 , t 2 , t 3 and it is negative on the intervals [t 1 , t 2 ] and [t 3 , ∞) (see Fig. 15 for a qualitative sketch of f λ,τ ). The late time violation is not important for us as long as by that time all the kinks have already been formed. Moreover their mutual interactions are exponentially suppressed on distances longer than 1/m in d = 1, and they can be completely neglected given that the average separation of the kinks is larger than ∼ (n K ) −1 max ∼ 6/m. On the other hand, the early time violation in the interval [t 1 , t 2 ] can be important as it might interfere with kink production and change the maximum kink number density. We can thus deduce three necessary conditions for the kink number density in the λ = 0 model to be a good approximation to that in the λ = 0 case: (i) The duration of early time violation of (146) needs to be finite i.e. t 2 < ∞.
(ii) All the kinks need to have been produced by the time the late time violation of (146) sets in i.e. t max < t 3 (iii) The duration of the early time violation of (146) needs to be much smaller than the fastest timescales of variation of the wave-functional i.e. ∆t ≡ t 2 − t 1 1/m.
We have swept the (λ, τ ) parameter space to determine the regions where the above conditions are verified. This has been done numerically by approximating f λ,τ via f λ,τ (t) ≈ 2|m 2 | − 3λ N n=1 |c n (t)| 2 − |c n (t 0 )| 2 , (148) and determining the corresponding values of t 1 , t 2 , t 3 for a wide range of values of λ and τ . The results are shown in Fig. 16 where we used the same numerical parameters as in Sec. II D. The regions shaded in red, orange and pink are excluded by the necessary conditions (i), (ii) and (iii) respectively. Alternatively we expect the λ = 0 model to be accurate inside the green region. Remarkably, the λτ /m = 1 curve lies deep inside this region which indicates that λτ /m 1 is a sufficient condition for the approximation to be valid.

VI. CONCLUSIONS
In this work we carried out a thorough analysis of the dynamics of topological defect formation in a quantum field theory where the only interactions are with external parameters that induce a quantum phase transition. We thus worked in the limit where self-interactions can be neglected. Results for the number density of kinks in one spatial dimension are summarized in Fig. 2, for vortices in two spatial dimensions in Fig. 9, and for monopoles in three spatial dimensions in Fig. 13. These results indicate that the number density of topological defects in d spatial dimensions scales as t −d/2 and does not depend on the quench time scale, in the late time limit. Moreover, we showed that the sudden phase transition analytical result is a universal attractor. These novel results stand in contrast to the Kibble-Zurek prediction for a thermal phase transition.
We have also discussed the limit within which our results can be expected to be a good approximation for a more realistic theory where self-interactions are not explicitly set to zero. In the case of kinks (d = 1) we found the condition λτ /m 1 where λ is the self-interaction coupling strength, to be a sufficient condition for our results to hold. This condition can be generalized on dimensional grounds to be λm d−2 τ 1 in d spatial dimensions.