System size scaling of topological defect creation in a second-order dynamical quantum phase transition

We investigate the system size scaling of the net defect number created by a rapid quench in a second-order quantum phase transition from an O(N) symmetric state to a phase of broken symmetry. Using a controlled mean-field expansion for large N, we find that the net defect number variance in convex volumina scales like the surface area of the sample for short-range correlations. This behaviour follows generally from spatial and internal symmetries. Conversely, if spatial isotropy is broken, e.g., by a lattice, and in addition long-range periodic correlations develop in the broken-symmetry phase, we get the rather counterintuitive result that the scaling strongly depends on the dimension being even or odd: For even dimensions, the net defect number variance scales like the surface area squared, with a prefactor oscillating with the system size, while for odd dimensions, it essentially vanishes.


Introduction
Rapidly sweeping through a second-order (quantum or thermal) symmetry-breaking phase transition in general may trigger the creation of topological defects. This is due to the fact that the system cannot adjust abruptly to the new order, so that (timedependent) correlations are spontaneously generated in the interacting many-body system. Correctly counting the topological defects generated by the sweep amounts to an analysis of these correlations.
The fundamental global requirement for topological defect generation is that the symmetry group of the new phase permits the defects in the sense that the homotopy group of the coset space, formed by the quotient space of unbroken (old) and broken (new) symmetries is nontrivial. The physical consequences of the latter mathematical fact were first used by Kibble [1] in the high-energy context of symmetry breakings from a (grand unified) symmetry group of generally high dimension to one of lower dimension in the earliest stages of the universe. This led to the classification of the symmetry group(s) for which cosmic strings and domain walls can, in principle, be produced by the universe's expansion after the big bang. An implementation of this idea directly accessible to laboratory-scale experimental realizations in condensedmatter systems was devised nine years later by Zurek [2]. Hence the historically developed and by now established name Kibble-Zurek mechanism for this particular symmetry-breaking based scenario of topological defect formation. While Kibble realized that in the early universe relativistic causality alone mandates the appearance of defects according to the nontrivial homotopy groups relevant to the symmetrybreaking scheme at hand, Zurek used scalings of the relaxation time and healing length with the quench time in the vicinity of the critical point to predict topological defect densities. The scaling behaviors derived by Zurek, based on the critical exponents of the phase transition under consideration, therefore led to concretely testable experimental signatures of the general idea of Kibble. The basic physical phenomenon is depicted in Fig. 1 for the simple case of a U (1) order parameter in the plane: It relies on the finite propagation speed communicating the new order and the corresponding correlations between spatially well-separated regions. These finite communication speeds lead to the whole sample not being able to establish the new order homogeneously throughout, thus causing an inhomogeneous order parameter distribution. Order parameter domains form, which, with some appropriate statistical probability, potentially enclose topological defects.
The Kibble-Zurek mechanism may occur in many diverse physical settings, for instance, in nonequilibrium phase transitions during the early universe [3,4], and has been measured in various condensed matter systems like liquid crystals [5] or superconducting films [6]. There has been a tremendously increasing interest in the closely related nonequilibrium quantum many-body phenomena during quantum phase transitions in recent years, cf., e.g., [7,8,9,11,10,12,13,14,15,16,17,18,19,20,21,22], due also to the fact that for controlled dynamical entanglement generation [23] the dynamical evolution of quantum many-body correlations is strongly relevant.
A clean experimental evidence for Kibble-Zurek mediated defect formation in dilute cold atomic gases, which are an experimentally precisely in situ controllable variant of conventional condensed matter systems, has been obtained in spinor as well as scalar Bose gases [24,25]. These dilute gases are particularly suitable to test predictions of the dynamical Kibble-Zurek theory because of their comparatively long re-equilibration times, which effectively permit a real-time study of the evolution of perturbations in the spatial many-body correlations and thus ultimately also a study of the creation process of topological defects. Furthermore, it was shown that the Kibble-Zurek type scaling of the defect number with the sweep rate through the phase transition depends on the inhomogenity of the sample [26,27] (is becoming more strongly dependent on the sweep rate), which is relevant for the intrinsically inhomogenous harmonically trapped gases. While most investigations to date have been carried out for bosonic systems in low dimension (where dimension here refers to real as well as order-parameter dimension), due to the rapidly growing interest in cold atomic gases, more recently fermionic problems in these systems have also been tackled, e.g. in [28,29].
When calculating defect densities from correlation functions, frequently the excitation densities (quasiparticle densities) and the topological defect densities actually created by the quench are assumed to be in a one-to-one correspondence cf., e.g., [18,20,21]. This identification is however misleading in general, because the necessary correlations created by the quench are much stronger for creating topological defects as compared to the correlations for quasiparticle excitations in high dimension [30]. For an illustration of this fact cf. Fig.2. A topological defect at a given vertex enclosed by domains after the quench implies that the direction vectors of all neighbouring domains either point away or towards the vertex (for a monopole configuration of the spin vectors, see also below). In a high-dimensional space, for example, this condition for a defect core to be enclosed by the domains at the vertex such that a winding of the order parameter orientation is induced becomes increasingly difficult to fulfill, given that the direction vectors orient themselves statistically inside the sample (i.e. due to quantum or thermal fluctuations establishing the order Origin of topological defect formation due to the Kibble-Zurek mechanism: Rapidly sweeping through a second-order phase transition statistically generates domains due to the finite propagation speed communicating of the new order in space and establishing the corresponding order parameter correlations. If by chance, at a given vertex, the domains meet such that a (generalized) "phase winding" is produced, a topological defect with a core of the old symmetry (yellow) is formed. In the simple example shown in the figure, a U (1) order parameter in the two-dimensional plane with field direction indicated by the red arrows, this amounts to the order parameter vector orientation winding around the core by 2π, for the planar spin hedgehog depicted. In high dimensions, a rapidly increasing degree of coordination of the order parameter orientation in neighboring domains is necessary to enclose a defect. In the three-dimensional example shown above for the classic monopole configuration, all arrows in the surrounding domains (red) have to point away from the central defect core (yellow), cf. the two-dimensional situation in Fig. 1 which already requires significantly less coordination.
In the following, we will investigate topological defect creation in quantum quenches concentrating on the behaviour of the defect creation rate for high internal and spatial dimension of the system. To this end, we derive essentially exact statements on topological defect production in the large N limit of an O(N ) symmetric model in D = N spatial dimensions. Similar to large-N approaches in condensed matter and field theory we argue that our predictions for the defect production will remain qualitatively valid also for smaller dimensions. More specifically, assuming that there is no critical value of N where the system behaviour changes in a discontinous manner, we expect our results to be qualitatively applicable also to finite N , e.g., N = 3, which are accessible to experimental tests, for example in spinor Bose-Einstein condensates with large spin, cf. the review [31].
In the next Section, we begin by deriving a large N expansion for a general O(N ) model resulting in a bilinear effective action, which is used to obtain the equations of motion describing the instability of direction modes during the phase transition. We then apply these equations to calculate the behavior of the spatial direction correlation functions shortly after the transition. The latter result is used to deduce ab initio the winding number variance for our O(N ) model. Finally, we delineate the behavior of the winding number variance with the size of the sample in dependence on the spatial behavior of the direction correlations (which is, e.g., short-range or long-range periodic).

Large-N Expansion
In order to illustrate the main idea -i.e., how to obtain an effective quadratic action for a N -component field Φ = (Φ 1 , Φ 2 , . . . , Φ N ) in the large N limit, let us consider the (Euclidean) generating functional of the non-linear O(N ) model where the factor 1/N in the second term ensures a well-defined large-N limit. For m 2 > 0, we get an O(N )-symmetric ground state, Φ = 0, whereas m 2 < 0 induces symmetry breaking, Φ = 0. We can formally split the field into modulus χ ≥ 0 and direction n with n 2 = 1 The functional integration measure then becomes where S N −1 = 2π N/2 /Γ(N/2) is the surface area of the unit sphere in N dimensions and thus we get for the generating functional with the remaining integration over χ (note that n ·ṅ = 0 by normalization) In view of the phase-space factor χ N −1 in the measure (4), we get, for N ≫ 1, a strong enhancement of contributions for large χ. On the other hand, the terms in the action in (6) such as χ 2ṅ2 yield a strong suppression at large χ. Therefore, the χ-integral will be dominated by the saddle point of the exponent in (6) given by In the large N limit,ṅ 2 is the sum of many independently fluctuating terms and thus approximately constant (the same is true for n · ∇ 2 n). Therefore, the saddle-point approximation yields a constant field χ, i.e., a classical (mean) field Due to the steep potential for the field χ (in the large N limit), its quantum fluctuations can be neglected and we can approximate it by a constant c-number. In the above formula, we have assumed a suitable ultraviolet regularization resulting in finite expectation values Φ 2 = O(N ) and Φ 2 = O(N ). Consequently, to lowest order in 1/N , we obtain the same results as with a bilinear effective action where the renormalized constants m 2 eff and c 2 eff depend on the shape of the nonlinear potential V . Of course, this line of argument is very similar to the nonlinear σ-model which does also approach a linear theory for N ≫ 1. As one may easily imagine, the same procedure can be applied in more general cases. For example, one may include arbitrary dispersion relations by adding terms like φ · ∇ 4 φ and φ · ∇ 6 φ etc. Furthermore, the prefactors in front of the different terms (such as c 2 ) may depend on φ 2 , similar to the potential V (φ 2 ).
The general idea is always the same: a sum of many (N ≫ 1) independently fluctuating quantities on an equal footing such as φ 2 is approximately a c-number whose fluctuations are suppressed by 1/ √ N (cf. the de Finetti theorem and the law of large numbers). Therefore, we shall use the same linear effective action also in the time-dependent case (with real time t instead of Euclidean time τ ) -and with a possibly time-dependent φ 2 (t), which is due to the growth of the thermal and quantum fluctuations (see below).

Phase transition and correlation function
Having motivated an effectively quadratic action in the large N limit, we can discuss the behaviour of the fluctuations during the phase transition, by using a linear evolution equation of the form where we assumed time-independent coefficients F (−∇ 2 ) and G(−∇ 2 ) initially and after the quench. Note, however, that the validity of such linear equations is not restricted to large N but a linearization is often also possible in lower dimensions, N = O(1), using different expansion parameters instead of 1/N . For instance, for a spinor Bose Einstein condensate, the gas parameter, measuring the diluteness of the system, can be employed, cf. [11]. Before the quench, the field is in a stable configuration, i.e., all excitations have real frequencies. The field operator can be written in terms of initial modeŝ (11) with N orthogonal polarizations, ε k,λ · ε k,λ ′ = δ λλ ′ . The coefficients A k depend on the functions F (k 2 ) and G(k 2 ) in Eq. (10) and are chosen such that the mode operatorsâ † k,λ andâ k,λ creating or annihilating a quasi-particle of wavenumber k and polarization λ obey the usual bosonic commutation relation [â k,λ ,â † k ′ ,λ ′ ] = δ kk ′ δ ′ λλ for a bosonic fieldφ. (In the following, we will take the basis vectors e a labelling the field components as the polarizations ε k,λ .) After quenching the system through a (second-order) phase transition, the O(N )symmetric state is no longer the ground state of the system. As a result, some modes become unstable and thus the (instantaneous) dispersion relation from (9) dips below zero for some wavenumbers, ω 2 k,> < 0, cf. Fig. 3. The time-dependence ofφ can be expressed in terms of instantaneous frequencieŝ with initial mode operatorsâ k,λ and complex coefficients B k and C k depending on the shape of the transition. Note that many of the modes are still stable after the quench and can be viewed as "virtual" short-lived defect-antidefect pairs annihilating immediately after their creation. These stable modes form a time-independent contribution to the correlation function. The field linearization outlined in the previous section (in the large N limit) is still valid for a finite time after the quench -as long as the relative contribution of unstable modes is small compared to the mean field Bearing in mind that we start in the O(N ) symmetric phase, which has φ = 0, we cannot simply insert n = φ /| φ | into the general expression for the winding number given by Eq. (20) below. Instead, in order to address the winding number statistics, we must employ some appropriately defined direction operatorn, see [30]. Since the new phase grows from unstable fluctuations, the basic domain structure with broken O(N ) symmetry and thus also the topological defects are seeded by the field distribution of the unstable modes shortly after the quench. Stable excitations, on the other hand, can be interpreted as virtual defect-antidefect pairs coming only briefly into existence. Especially for large k, these short-lived virtual defect-antidefect pairs would dominate the defect number in a given volume for any fixed time due to the phase space factor k N −1 stemming from the d N k-integral. This corresponds to the strong singularity of the two-point function Since we are not interested in these short-lived pairs, but rather in real long-lived defects originating from unstable fluctuations, we introduce a time-averaged direction operatorn with (temporal) regularization function g(t). This time average supresses all short-lived pairs at large k. Since the spatial behaviour of the n-correlator (16) will be dominated by one wavenumber k * , see Eq. (18) below, and the time-dependence cancels due to the normalization of n, our result is independent of the choice of g(t). The only conditions on g(t) are that its support lies well within the region of applicability of the effective action (9) and that g(t) is smooth enough. More precisely, the spectral width ofg(ω) determines the uncertainty of k = k * ± ∆k which should be small ∆k ≪ k * .
Via the same arguments as in the previous Section, the normalization factor Z = O(N ) can be approximated by a c-number for large N Obviously, the rapidly oscillating parts ofφ, i.e., stable modes, are smoothed out in Eq. (13) and only the exponentially growing excitations remain. Note that the leading order of the normalization factor,Ẑ ≈ Z, would yield exact normalization of the expectation value, n 2 = 1, but only approximate normalization of the operator itself,n 2 ≈ 1. Proper normalization of the direction operator,n 2 = 1, can only be achieved by including operator-valued higher orders ofẐ as well. In view of the mode expansion (12), we get for the time-averaged direction operator where the integration domain in k space runs over all unstable modes. The coefficients C k include the above time-averaging and the exponential growth of the unstable modes, as well as the original coefficients B k and C k , and possibly the initial distribution, e.g., of thermal origin, of the occupation numbers of the different modes. The calculation of the correlation function is now straightforward. Taking the basis vectors e a labelling the components φ a of the field as the polarizations ε k,λ and approximating the normalization factor by its expectation value, we obtain In the case of spatial homogenity and isotropy (as opposed to internal isotropy of the N fields φ a ) after the quench, where no direction is distinguished, we can adopt spherical coordinates and carry out the angular integration where the index ν of the Bessel function J ν is given by ν = N/2 − 1. The remaining k integral is often dominated by one wavenumber k * : In large dimensions N ≫ 1, the phase-space factor k N −1 strongly favors contributions near the largest unstable wavenumber, cf. Fig. 3, whereas for small N = O(1), the dominant wavenumber is usually determined by the largest coefficientC k , i.e., the fastest-growing mode. In either case, as a result of the existence of a dominant wavenumber k * and the ensuing application of the saddle-point integration, the correlation function of the time-averaged direction vector now reads Note again that this is only a part of the direction correlations, generated by the unstable modes, and thus responsible for the formation of topological defects. All stable modes have been averaged out in definition (14) of the direction operatorn.
For large N → ∞, we can use the well-known asymptotic behavior of Bessel functions [32] and obtain Gaussian correlations which fall off on a typical distance √ N /k * . Note that, since the corrections to the normalization factorẐ of the direction operator are suppressed for large N , cf. (14), these Gaussian correlations represent an exact result. The typical oscillatory behavior of Bessel functions is exponentially suppressed for large N → ∞. For finite N , the spatial decay of the correlations is too weak such that the oscillations are retained. However, for finite N , one should also be aware that corrections to the normalization Z would become important, though the leading order of the directional correlations would still be given by expression (18).

Winding numbers in dimension N : Classification of topological defects
As a quantum version of the classical winding number counting the net number of topological defects within a domain M in N dimensions [33], we introduce a winding number operatorN for the direction operatorn of the O(N ) field, cf. [30] where S N −1 = 2π N/2 /Γ(N/2) is the surface area of the N dimensional unit sphere. The operatorN (and similarly the classical winding number) remains unchanged if we subject the fields n to a global internal O(N ) rotation which can be smoothly deformed to the identity transformation. Improper rotations, which are also O(N ) symmetry transformations (but cannot be smoothly deformed to the identity transformation), generally do not leaveN invariant.
To exemplify this, let us consider the simplest nontrivial case, N = 2. Parametrizing the field direction through a single angle,n 1 = cosθ andn 2 = sinθ, we obtainN where dl is a tangent vector along the integration contour C. A rotation of the field directionn ′ = D(φ)n by an angle φ will only add a phase φ to the angle operatorθ and the winding number (21) is unaffected, see also Figure 4. Besides rotations, the group O(N ) also contains improper rotations, i.e., combinations of a rotation with a reflection, which do not conserve winding number. This leads us, apart from the wellknown monopole or spin-vortex configurations depicted in Fig. 4, to field distributions with a nontrivial winding configuration, where the field points away in one direction and towards the defect core in the other, see Fig. 5. This can be generalized to higher dimensions. For example, in three dimensions, any rotation ofn can be parametrized by three Euler angles, i.e., as a sequence of three two-dimensional rotations. Obviously, inversion of an odd number of field direction, e.g., n → −n, cannot be expressed in that way but corresponds to improper rotations again so that the winding number picks up an additional factor −1. Hence, the typical  N = ±1 configurations have the field pointing either away from (+1) or towards (−1) the defect core (in all directions). Other nontrivial field distributions, e.g., those with n pointing away in some and towards the core in all other directions, can be obtained through rotations. The strikingly different appearance of positive and negative defects in two and three dimensions is typical for even and odd number of dimensions N : Field inversion, n → −n, will leave the winding number N → (−1) N N invariant in even dimensions while for N odd, N picks up a factor (−1) N = −1. As we will see in an example below, this has implications for the net defect number created during a quench to the broken-symmetry phase.

Winding number variance
Since, for a fully O(N ) invariant initial state, the probabilities for creating positive and negative defects are equal, the expectation value N of the winding number operator (20) is zero, but the variance N2 is generally non-vanishing and can be used to quantify the net number of defects created during the quench. Due to the totally antisymmetric tensor ε abc... on the right hand side of Eq. (20), each field componentn a must appear exactly once inN. Hence, if different field components are uncorrelated, cf. Eq. (18) the expectation value N 2 must break down into the product of N two-point functions.
The first identity in (22) holds when the system is homogenous, while the second is valid when the system is isotropic as well. This particular form of the two-point function n a (r)n b (r) can be inferred from symmetries: Although the phase with broken O(N ) symmetry is neither homogeneous nor isotropic, the expectation values still are, because we start from an homogeneous and isotropic phase and any realization of the symmetry breaking is equally likely. With (22), there results the following expression for the winding number variance where we exploited the symmetry under exchange of index pairs (a, α) ↔ (b, β), and also changed the order of taking spatial derivatives and expectation values. Due to the total antisymmetry of the integrand, this expression can be regarded as a collection of surface integrals over (N − 1)-forms (in the unprimed indices while keeping the primed indices fixed and vice versa). It becomes therefore possible to apply the generalized Stokes theorem, M dω = ∂M ω, and rewrite (23) as a volume integral over N -forms This expression is generally valid provided different field direction are uncorrelated, n anb ∝ δ ab , cf. (22). The next step is less obvious, but again we can take advantage of the inherent symmetries which enabled us to write the two-point function (22) through a function depending only on the separation L = |r − r ′ |. Using the identity the derivatives under the integral can be evaluated. In view of the total antisymmetry of the Levi-Civita symbols, all terms containing more than two first derivatives of L 2 , e.g., (∂ α L 2 )(∂ ′ α ′ L 2 )(∂ β L 2 )(∂ ′ β ′ L 2 ) cancel after summation. It follows So far, we did not specify any particular coordinates, but in order to simplify this expression, it is advantageous to adopt Cartesian coordinates for which ∂ α ∂ ′ α ′ L 2 = −2δ αα ′ . After performing the summations over α, α ′ etc. and with ∂f /∂L 2 = (1/2L)∂f /∂L, we finally obtain which holds for arbitrary integration manifolds M and any dimension N provided the two-point correlator has the form (22), i.e., the correlator n a (r)n b (r ′ ) depends only on the separation between r and r ′ . If, however, homogenity or isotropy are broken but different directions are still uncorrelated, the most general expression for the winding number variance is given by Eq. (24).

Spherical volumina
The 2N integrals appearing in (27) are, for arbitrary volumina, generally difficult to evaluate directly, but can be simplified drastically by exploiting symmetries of the integration manifold M. In particular, we will demonstrate that for a sphere of radius R -an isotropic manifold defined through M = {r : r 2 < R 2 } -Eq. (27) reduces to a single integration over one remaining variable. To this end, let us instead of the winding number variance N 2 consider its change with radius R of the integration volume. Taking the derivative of (27) with respect to R and using spherical coordinates, one of the radial integrations will break down to yielding a factor R N −1 from the integral measure [note that the integrand of (27) is symmetric in r and r ′ ]. One of the angular integrations d N −1 Ω ′ can be performed due to isotropy of the manifold M and gives the surface area S N −1 of the N dimensional unit sphere. For the radial derivative of N2 , we get The remaining volume integral d N r = J N drd N −1 Ω with the Jacobian J N can be tackled by a coordinate transformation r → L = r − r ′ . Adopting spherical coordinates for L, we see that the (radial) measure L N −1 just cancels the factor 1/L N −1 such that the remaining L dependence is just a partial derivative and the modulus integral dL can be performed easily. Regarding the angular part, we exploit isotropy of M once more and choose the angles φ 1 ...φ N −1 such that the separation L = 2R sin φ 1 depends only on φ 1 . After substituting y = sin φ 1 and performing the integrations over the remaining angles, we then further get Using the identity 1 − y 2 N −3 = −(1/N − 1)y N ∂( 1 − y 2 N −1 /y N −1 )/∂y, Eq. (29) can be rewritten and, after integration by parts, the change of the winding number variance reads In the factors right of the derivative, the integration variable y appears only in the combination Ry. Hence, using the chain rule, ∂/∂y can be replaced by (R/y)∂/∂R and it follows Both sides of this equation are partial derivatives with respect to the manifold radius R and can be formally integrated, where the constant of integration follows by requiring N 2 (R = 0) = 0. Finally with y = sin(θ), we can simplify the result to the more compact form i.e., the winding number variance can be determined through a single angular integration. For short-ranged correlations, where ∂f /∂L falls off fast enough at large a surface scaling [11,30] follows for large R: For sufficiently large R, the main contribution to (32) arises from small θ. By formally extending the integration to infinity and approximating cos θ ≈ 1 and sin θ ≈ θ, one can see that the integral is proportional to 1/R and we obtain scaling of the winding number variance with the surface area of the enclosed volume. Note that expression (32) and thus also the surface scaling for short-ranged correlations (33) hold in any dimension N as long as the direction correlator takes the form (22). In particular, the expression for N2 derived in [11] will be reproduced for N = 2.
We stress that our calculation based on very general principles predicts a surface scaling law for the net defect number variance. This surface law should be contrasted with the volume scaling law obtained from another frequently used and phenomenologically motivated approach [2]. If one would randomly distribute hedgehogs and antihedgehogs in a given N -dimensional space with the same fixed densities, the total number of defects (hedgehogs plus antihedgehogs) inside a given volume would scale with the volume itself. Their difference (i.e., hedgehogs minus antihedgehogs), however, would scale with the square root of that volume (corresponding to a one-dimensional random walk of the winding number) and thus the net defect number variance scales with the volume.
As far as we are aware, the rather generic result obtained here, together with our previous results in [11,30], represent the first proof by direct calculation that (the quantum version of) the Kibble-Zurek mechanism leads to a surface law.

General convex volumina
The surface scaling for short-ranged interactions can be interpreted as a variety of a "confined" phase, where bound defect-antidefect pairs, typically separated by a correlation length, occur. Pairs entirely within or outside the manifold M will not contribute to the net winding number and thus also not to N2 . Only those pairs, where one of the partners is inside and the other outside will contribute. Since it is not clear whether the defect or the antidefect is contained within M, calculation of the winding number would be similar to a random walk with number of steps proportional to surface area. From this simple picture, it can be anticipated that the surface scaling is not restricted to spherical manifolds but might apply in any situation where the surface of the manifold is sufficiently "smooth".
In the following, we will thus consider arbitrary convex manifolds and show that the surface area scaling of the winding number variance is retained. We parametrize the convex manifold in spherical coordinates through a scale factor R and a shape function a(Ω) = O(1), similar to a sphere, where we would have a(Ω) ≡ 1. Obviously, the surface area ∂M ∝ R N −1 so that we need to show N 2 ∝ R N −1 as well. For convex manifolds, the volume element in expression (27) for the winding number variance reads in spherical coordinates where the Jacobian J N (Ω) is taken at r ≡ 1. As for spherical volumina, cf. (28), the scale factor R appears only in the upper bound of the radial integration such that we can use the same trick as before and calculate ∂ N 2 /∂R instead of the variance itself However, in contrast to the spherical case, we cannot simply perform the angular d N −1 Ω integration because of broken isotropy (of the manifold): The separation L = |r − r ′ | generally depends on both azimuthal positions Ω as well as Ω ′ and not only on their relative angles. The coordinate transformation r → L = r − r ′ and subsequent integration in spherical coordinates yields where L max is the distance from a point Ω ′ on the surface ∂M to the opposite boundary of M in directionΩ and is determined through [At this point, it should be mentioned that one needs to be careful regarding the azimuthal positions: Ω ′ and Ω are those of r ′ and r as seen from the predefined centre of the manifold, cf. (35), whereasΩ is the azimuthal position of r as seen from r ′ .] A scaling of the winding number variance with surface area is obtained if the integral in (38) is proportional to 1/R for sufficiently large R -which can indeed be expected for short-ranged correlations (33) because r max and r ′ max both scale with R and thus L max = RΛ max (Ω, Ω ′ ) ∝ R.
We evaluate the scaling behaviour of (38) by choosing suitable coordinates for the integration: In view of the convexity of M, the anglesΩ can be chosen such that for each Ω ′ there exists at least one angle φ 1 ofΩ for which the separation L max = RΛ max changes from zero to order R while keeping all other angles φ 2 to φ N −1 fixed at arbitrary values. The φ 1 integral appearing in (38) will be performed first, it reads Obviously, the integrand vanishes for all φ 1 with Λ max (φ 1 ) = 0. It becomes thus possible to change the domain of integration to . From convexity of the integration manifold M also follows that the maximal separation Λ max cannot have (local) minima in (φ − , φ + ), i.e, Λ max must be monotonically growing in the vicinity of φ − and decreasing near φ + . Furthermore, by demanding short-ranged correlations, (33), ∂f /∂L becomes arbitrarily small for large separations. (See, by contrast, the case of long-range correlations we discuss in the following section). This means that, for sufficiently large R, the main contributions to integral (39) stems from regions where Λ max is small but non-zero, i.e., near φ − and φ + , and we can approximate (39) with small ǫ ≪ 1. For all angles φ ∈ (φ − + ǫ, φ + − ǫ), the distance Λ max > 0 but the integrand is negligibly small due to short-range correlations (33). In the integrals on the right hand side of (40), we can assume that Λ max is strictly monotonic and analytic such that Λ max (φ 1 ) is bijective in (φ − , φ − + ǫ) and similarly in (φ + − ǫ, φ + ). It becomes therefore possible to change the integration variable from φ 1 to Λ max and we finally obtain for the φ 1 integral (39) where we have extended the integration to infinity because, due to short-range correlations, the integrand becomes negligibly small for sufficiently large R. In this limit, the integral yields a constant (independent of system size R but of course still depending on all other angles φ 2 ...φ N −1 as well as Ω ′ ) and the scale factor R appears only as prefactor 1/R. We are thus led to the scaling of (38) with R N −2 for large R and thus, after integration, obtain for the winding number variance i.e., the desired surface scaling result. It should be noted that the R N −2 dependence of ∂ N2 /∂R derived in this section is only valid for sufficiently large R. When integrating ∂ N 2 /∂R to obtain the scaling behavior of the winding number variance, regions with small R, where ∂ N2 /∂R might behave differently, are also included, and potentially contribute subleading orders to the winding number variance, which scales with surface area for large R.

Deviations from the Surface Law
We have obtained the surface scaling (32) for spherical volumes in high dimensions, and generalized this result to arbitrary convex volumes (42). The salient ingredients of the derivation have been the short-range nature of correlations, Eq. (19), and the symmetries of the correlation function (22). By relaxing either or both of these conditions, and depending on the precise functional form of the two-point correlator, it is then possible, as we will now show, to obtain various interesting scaling behaviors with system size. Let us consider broken spatial isotropy of the correlation function (22), e.g., through an externally imposed lattice. (Note that n anb ∝ δ ab is still isotropic in the field directions.) Explicit examples for periodic and short-ranged correlations will be discussed in subsections 5.1 and 5.2. To remain general, we assume that in (22) only the first equality holds, and that the (real space part of) the correlator separates into the product where the correlations g α in different directions can in general have a different functional behavior. We insert this expression into (24), and consider the volume enclosed by a N -dimensional box with size L 1 × L 2 × L 3 .... The second-order partial derivative of the correlator (43) gives where a prime denotes the derivative, e.g., g ′ α = ∂g α (x)/∂x. Inserting the above expression into (24), we see that any terms involving products etc. cancel during summation due to antisymmetry of the ε pseudotensors and we have This correlation function can be evaluated further by accounting for the broken isotropy of f and regarding the integral over each each spatial direction independently. It follows So far, this result is general. We now are going to illustrate it with a few concrete examples for the correlation function.

Periodic Order in the Direction Operator Correlations
As an example for a nonlocalized correlation function, we consider a purely periodic functional behavior, which can, for example, be derived within the following simple model. Using the ansatz for the fluctuating order-parameter field and explicitly break the isotropy of the system by imposing periodicity along all N spatial axes. The integrals in (46) now grow with system size L α for an even number of dimensions (the integrand is strictly positive for even N ), whereas the first term in the sum, g N α (0) − g N α (L α ) = 1 − cos N (kL α ), is a periodic function in the system size L α . Thus the correlation function (48) implies a scaling of the winding number variance (for even N ) with the square of the surface area of integration volume V = L 1 × L 2 × L 3 ...
The factors 1 − cos N (kL α ), however, will still lead to vanishing winding number variance if the size L α = nπ/k matches the periodicity of the correlations in all spatial directions. For an odd number of dimensions, on the other hand, the integrands of (46) oscillate periodically around zero and no scaling of N2 with system size is obtained, This entirely different scaling -surface squared with oscillating prefactor for even dimensions versus essentially no net defects in odd dimensions -follows directly from the behavior of the winding number under field inversion described in Sec. 4.1 (see also Figs. 4 and 5): Due to the strict periodicity of the correlation function (48), the field orientation must be inverted after half a unit cell n(r) = −n(r + e a π/k) in every direction e a , see Fig. 6. Hence, for each defect with topological charge N = +1, there sits another defect with charge N = (−1) N , a distance π/k apart in every lattice direction. This means for even dimensions, one has alternating layers of defects and antidefects, and thus the surface squared scaling with oscillating prefactor results, whereas in odd dimensions, each layer contains positive and negative defects, so that the total winding number essentially vanishes for a box-like integration manifold. (However, if the box-shaped integration manifold were tilted with respect to the lattice, a non-vanishing net defect number could occur in odd dimensions as well.)

Localized (Gaussian) Order
As a second example, we take short-ranged correlations, e.g., a Gaussian correlator similar to (19) derived for spatial isotropy in large N . Now, each of the double integrals in (46) scales with L β , such that the winding number variance is proportional to the surface area of the enclosed volume The above relation implies that the surface scaling found for spatially isotropic shortrange correlations, cf. Eq. (42), is shown to be preserved also for the case of broken spatial isotropy. This should be contrasted with the result of the previous subsection, where it was shown that for long-range periodic order, no scaling proportional to the surface area is obtained.

Conclusion
Starting from a rather general O(N ) symmetric action, we derived a quadratic effective action for the field fluctuations, and studied the instability of modes during a phase transition from the initial O(N ) symmetric ("paramagnetic") phase to a broken symmetry phase in which "ferromagnetic" order develops. Since we are interested only in the exponentially growing modes, which dominate the development of the domains and thus also of topological defects, we introduced a systematic averaging procedure for the fluctuating field operator directions. We have demonstrated that the direction field correlation function can be explicitly evaluated. Our results are very general: The initial state (such as temperature, dispersion relation) enters the coefficients B k and C k in (12) only. Thus, while the initial state affects the field correlator φ a (t, r)φ b (t, r ′ ) , it leaves the direction n-correlator (16) invariant. The external conditions after the transition determine the final dispersion relation which yields the value of k * . Apart from a simple rescaling of this value, our main results for the homogeneous and isotropic case (surface scaling law etc.) thus do not depend on the initial and the final state. This correlator of the field direction then directly determines the statistics of the net defect number created by the quench. A general expression for the winding number variance was derived and was shown to scale with the surface area of the enclosed convex volume in any (large) dimension N . Such a surface area scaling follows from general arguments based on the assumption of short-ranged correlations together with the internal O(N ) and spatial symmetries (homogenity and isotropy). For large enough samples, it should be possible to verify this prediction experimentally. The detection method depends on the specific experimental setup: In Bose-Einstein condensates, the topological defects are phase or spin vortices (in two spatial dimensions) which can be imaged either by absorption (resolving the core as a hole in the condensate [25]), or using polarization-dependent phase-contrast imaging [34]. In superfluid 3 He, vortices can be detected with NMR techniques [35]. Finally, in conventional condensed matter systems with magnetic order like ferromagnets, topological defects can be detected for example with Lorentz transmission electron microscopy [36].
If the correlations are not short-ranged, the scaling law can be very different: E.g., for purely periodic (i.e., long-ranged) correlations on a lattice, the scaling law strongly depends on the parity, i.e., whether the number of dimensions N is even or odd. In even dimensions, the winding number variance basically scales with surface area squared with oscillating prefactor. By contrast, for odd N , the net defect number inside a box of size V = L 1 × L 2 × L 3 ... effectively vanishes. This entirely different behavior can be explained by the strict periodicity of the field directions together with properties of the winding number operator, i.e., how positive and negative defects are aligned on a lattice: When inverting the field direction, n → −n, the winding number N acquires a prefactor (−1) N . Hence, configurations where the field points towards to and those where it points away from the defect core must lie in the same homotopy class (N = +1) for even dimensions and in different classes (N = ±1) in odd dimensions. The strict periodicity of the correlation function, on the other hand, implies a field inversion after half a period (along any of the lattice directions). Thus, one has alternating layers of positive and negative defects for even dimensions, while for odd dimensions positive and negative defects appear in each layer.
Further possible directions of research include the extension of our results for the winding number variance to finite, experimentally accessible values of N . We stress in this context that the replacement of the quantum normalization factorẐ introduced through the time-averaging of the direction vector, by its classical expectation value in the correlator (16), remains the only finite N approximation made in our analysis. This statement holds under the condition that we have a system for which the effective action for the field remains still quadratic for finite N (like in a spinor Bose-Einstein condensate) leading to a dispersion relation with a dominant wavenumber k * > 0. The finite N extension can then be achieved by evaluating Eq. (16) to the next order(s) in 1/N , i.e., by calculating the corrections to the operator-valued normalization factor (14) entering the average quantum direction vector in (15).