text:
# Quantum tunneling of superconducting string currents.

## I Introduction

Our universe may contain cosmic strings which are relics of a symmetry breaking phase transition at early times. (For a review, see .) In many models these strings can be superconducting, and the resulting current can have important effects on the string dynamics. In particular, if the current is large enough, the pressure due to the charge carriers can oppose the tension of the string and lead to a stable “vorton” state .

The string current depends on the winding number density of the phase of a condensate field; in a loop the total winding number will be an integer, and thus classically cannot change as long as this phase is well-defined. When the current becomes large enough, there is the possibility of a ”quench”, where the condensate is driven to 0, and thus the phase can unwind and current can be lost. However, even below this threshold current, there is the possibility that current can be lost through quantum tunneling effects. It is those effects which we consider here.

The simplest field theory model that exhibits bosonic superconductivity is of the form,

$$=\overline{D_\mu \sigma }D^\mu \sigma +_\mu \overline{\varphi }^\mu \varphi \frac{1}{4}f_{\mu \nu }f^{\mu \nu }V(\sigma ,\varphi ),$$
(1)
where

$`f_{\mu \nu }`$ $`=`$ $`_\mu Z_\nu _\nu Z_\mu ,`$ (2)
$`D_\mu \sigma `$ $`=`$ $`(_\mu ieZ_\mu )\sigma .`$ (3)
This model has a local $`U(1)_{local}`$ symmetry associated with the gauge field $`Z^\mu `$ and a global symmetry $`U(1)_{global}`$ of the phase transformations of the field $`\varphi `$. In order to have a superconducting string model, we should also impose that the interaction potential between the two complex scalar fields, $`\varphi `$ and $`\sigma `$ is such that the broken $`U(1)_{local}`$ symmetry gives rise to a cosmic string configuration and the $`U(1)_{global}`$ becomes broken only inside the string. A bosonic condensate will then be formed on the string, with massless Goldstone bosons playing the role of “charge carriers”. Note however that in this model the condensate field $`\varphi `$ is not coupled to any gauge field, so the current is neutral.

Low energy excitations of the condensate along a straight string in the $`z`$ direction will have the form

$$\varphi (x,y,z,t)=\varphi (x,y)\mathrm{exp}[i\theta (t,z)].$$
(4)
Using this ansatz, we can write an effective action for the field $`\theta `$ by inserting Eq. (4) in the four dimensional action to get,

$$S=\mathrm{\Sigma }𝑑z𝑑t[(_t\theta )^2(_z\theta )^2]$$
(5)
where

$$\mathrm{\Sigma }=𝑑x𝑑y|\varphi |^2.$$
(6)
We can now define the neutral current associated with the field living on the string worldsheet by§§§Here and below, Latin indices take values 0 and 1, Greek indices take values from 0 to 3, and the signature of the metric is $`g_{00}>0`$ both in four dimensions and on the string worldsheet.

$$J_a=_a\theta ,$$
(7)
so the current is a measure of the winding number per unit length of the phase of the bosonic condensate on the string.

If we imagine the straight string as part of a large loop, we can write the total winding number

$$N=\frac{1}{2\pi }_\sigma \theta d\sigma ,$$
(8)
which is classically conserved. Nevertheless, this is not the whole story. It was already noticed by Witten in his original paper on superconducting strings that quantum effects would be able to decrease the current. He argued that since the field would have to unwind in order to decrease the value of N, the process would have to involve some kind of localized region on the string worldsheet where the absolute value of the field $`\varphi `$ would vanish.

Several authors have estimated the decay rate of the current following similar arguments. In this paper we shall study the quantum tunneling of the current by analytically computing the worldsheet instanton for this process and the corresponding decay rate in a loop with a uniform current.

## II Worldsheet Instanton

As a first approximation we will study the process on a straight string and will neglect any backreaction effect on the profile of the string. Therefore, we will only be concerned with the evolution of the field $`\varphi `$ in a fixed background. Also, we will not consider situations where the unwinding of the field takes place differently at different locations in the plane transverse to the direction of the string. This enables us to integrate over transverse directions, and use only the effective Euclidean action for the field $`\theta `$,

$$S_E=\mathrm{\Sigma }𝑑z𝑑\tau [(_\tau \theta )^2+(_z\theta )^2]$$
(9)
The Euclidean field equation for $`\theta `$ is then

$$_a_a\theta =0.$$
(10)
The limits of validity of this approximation will be indicated below.

Let us first assume that the current is spacelike,

$$J_aJ^a<0.$$
(11)
Then we can choose a frame of reference where the charge density vanishes, $`J_0=0`$ (pure current). We are interested in studying the quantum tunneling process from an initial state of homogeneous pure current of winding $`N`$,

$$\theta =\frac{2\pi N}{L}z,$$
(12)
where $`L`$ is the length of the loop, to a final state of winding $`N1`$. The corresponding initial current is therefore,

$`J_0`$ $`=`$ $`0,`$ (13)
$`J_z`$ $`=`$ $`{\displaystyle \frac{2\pi N}{L}}.`$ (14)
Semiclassically, the tunneling is described by a Euclidean instanton solution $`\theta (\tau ,z)`$ which approaches the initial configuration (12) at $`\tau \pm \mathrm{}`$. The field configuration immediately after tunneling is given by $`\theta (0,z)`$. Since the winding along the $`z`$-axis $`\tau =0`$ should be one unit less than that at $`\tau \pm \mathrm{}`$, it is clear that the field $`\theta (\tau ,z)`$ should have a vortex below and an anti-vortex above the $`z`$-axis (or vice versa). The corresponding solution of Eq. (10) is

$$\theta (\tau ,z)=\mathrm{arctan}\left(\frac{z}{\tau a}\right)\mathrm{arctan}\left(\frac{z}{\tau +a}\right)+\frac{2\pi N}{L}z,$$
(15)
where $`2a`$ is the vortex separation. This solution is well known from hydrodynamics , where $`_a\theta `$ corresponds to the velocity of an incompressible fluid. The solution describes a vortex-antivortex pair on the background of a fluid of constant velocity. The vortex separation $`a`$ is chosen so that the two vortices are in equilibrium due to the balance between their mutual attraction and the oppositely directed Magnus forces exerted by the background. This gives

$$a=\frac{L}{4\pi N}.$$
(16)
The equilibrium is unstable, as it should be, since the instanton solution should have a negative mode  .

For our purposes, instead of the variable $`\theta `$ it will be convenient to introduce the dual field $`\chi `$ defined by

$$_a\chi =ϵ_{ab}_b\theta .$$
(17)
Like $`\theta `$, $`\chi `$ obeys Laplace’s equation, $`_a_a\chi =0`$ away from the vortex cores, but a vortex in $`\theta `$ becomes a source for $`\chi `$. This can be used to write the equation for $`\chi `$ in the whole space,

$$_a_a\chi =2\pi \rho ,$$
(18)
where the “charge density” $`\rho `$ is given by

$$\rho =\stackrel{~}{\delta }(𝐱𝐱_+)\stackrel{~}{\delta }(𝐱𝐱_{}).$$
(19)
Here, $`𝐱`$ is a vector in the $`(\tau ,z)`$-plane, $`𝐱_\pm =(\pm a,0)`$ are the locations of the vortex centers, and $`\stackrel{~}{\delta }(𝐱)`$ is a delta-like function normalized so that

$$\stackrel{~}{\delta }(𝐱)d^2x=1.$$
(20)
The width of the function $`\stackrel{~}{\delta }(𝐱)`$ is comparable to that of the vortex core and its precise shape will not be important for our considerations. We thus have an analogy with 2-dimensional electrostatics. The problem to be solved is that of two opposite charges of magnitude $`1`$ in a constant external field of magnitude

$$𝐄=(2\pi N/L,0).$$
(21)
The corresponding solution of (18) outside the sources is

$$\chi =\mathrm{ln}\frac{|𝐱𝐱_+|}{|𝐱𝐱_{}|}+𝐄𝐱.$$
(22)
The mutual attraction of the charges must be balanced by the force of the external field acting on them in opposite directions, which occurs when $`E=1/2a`$, in agreement with (16). See Fig 1.

The Euclidean action corresponding to the field equation (18) can be written as

$$S_E=\mathrm{\Sigma }d^2x[(_a\chi )^2+4\pi \rho \chi ].$$
(23)
The overall factor has been chosen so that the actions (9) and (23) coincide when $`\chi `$ is a solution of (18) and $`\theta `$ is related to $`\chi `$ by Eq. (17).

The current decay rate is given by

$$\mathrm{\Gamma }=Ae^B,$$
(24)
where

$$B=S_ES_E^{(0)},$$
(25)
$`S_E`$ is the instanton action and $`S_E^{(0)}`$ is the action for the constant current configuration (12). The semiclassical approximation applies when $`B1`$. In this regime, the tunneling rate is determined mainly by $`B`$, and the value of the prefactor $`A`$ is not very important.

### A Evolution subsequent to the decay.

It is easy to check that our solution (15) conserves energy and momentum as well as charge on the worldsheet. We can also analytically continue it into Lorentzian time to obtain the evolution of the current after the tunneling process has occurred. By calculating the charge density on the Lorentzian worldsheet after the current decay, we can see (Fig. 2) that this solution represents the propagation of two lightlike trains of opposite charge travelling in opposite directions along the string. They leave behind a region of the same winding number density as the initial conditions, but if we now calculate the total winding along the string we see that

$$\frac{1}{2\pi }_\sigma \theta d\sigma =(N1).$$
(26)
This is because the two travelling waves have each negative winding number density, in other words they have the opposite winding to the original current.

## III Decay rate

In order to calculate the decay rate for this process we have to compute the Euclidean action for the instanton, or in terms of the electrostatic analogy, the difference in energy between the pure background configuration and the equilibrium solution described above. However, both energies are infinite, due to the infinite two-dimensional volume occupied by the background electric field. Moreover, one finds that a formal large-distance cutoff gives different answers, depending on the details of the cutoff procedure. We have therefore introduced a “physical” cutoff by replacing the homogeneous electric field with a field produced by a finite distribution of charges, $`\rho _{ext}(𝐱)`$. The solution (22) is then replaced by

$$\chi =\mathrm{ln}\frac{|𝐱𝐱_+|}{|𝐱𝐱_{}|}+\chi _{ext}(𝐱)\stackrel{~}{\chi }(𝐱)+\chi _{ext}(𝐱),$$
(27)
where $`\chi _{ext}(𝐱)`$ is the solution of (18) with the source $`\rho _{ext}(𝐱)`$ and it is assumed that $`\chi _{ext}(𝐱)`$ vanishes at $`|𝐱|\mathrm{}`$. It is assumed also that the external charges are located at distances $`|𝐱|a`$, and in the end we take the limit when these charges are removed to infinity, while keeping the electric field near the origin fixed and equal to (21),

$$d^2x\rho _{ext}(𝐱)\frac{𝐱}{r^2}=𝐄.$$
(28)
We also have to deal with the short-distance logarithmic divergence in the self-energy of the two charges. This divergence is an artifact of our approximation treating the vortices as pointlike objects. It is easily fixed by replacing $`\mathrm{ln}0`$ by $`\mathrm{ln}\delta `$, where $`\delta `$ is the characteristic radius of the vortex core. Although $`\delta `$ is defined only by the order of magnitude, the result is not sensitive to the precise definition, due to the logarithmic character of the divergence.

After integration by parts in Eq. (23), the energy difference between the instanton configuration and the background can be expressed as

$$B=2\pi \mathrm{\Sigma }\left(\rho \chi d^2x\rho _{ext}\chi _{ext}d^2x\right)=2\pi \mathrm{\Sigma }[\chi (𝐱_+)\chi (𝐱_{})]2\pi \mathrm{\Sigma }\rho _{ext}(𝐱)\stackrel{~}{\chi }(𝐱).$$
(29)
Since the support of $`\rho _{ext}(𝐱)`$ is at $`|𝐱|a`$, we can use the large-distance expansion of $`\stackrel{~}{\chi }`$,

$$\stackrel{~}{\chi }\frac{2𝐱𝐚}{r^2},$$
(30)
in the last term of (29). With the aid of Eq. (28), this gives

$$B=4\pi \mathrm{\Sigma }[\mathrm{ln}(2a/\delta )2Ea],$$
(31)
where $`\delta `$ is the vortex core radius, as explained above. Note that the result (31) is insensitive to the details of the external charge distribution $`\rho _{ext}(𝐱)`$.

The action (31) is maximized for $`a_0=1/2E`$, in agreement with (16). It is easily verified that $`d^2S/da^2(a_0)<0`$, indicating that $`a=a_0`$ is a local maximum of the action. In terms of the electrostatic analogy, this indicates that $`a=a_0`$ is an unstable equilibrium, as it should be for an instanton. Substituting $`a=a_0`$ in (31) we obtain finally

$$B=4\pi \mathrm{\Sigma }[\mathrm{ln}(L/2\pi N\delta )1],$$
(32)
and for the tunneling rate

$$\mathrm{\Gamma }A\left(\frac{2\pi N\delta }{L}\right)^K,$$
(33)
where $`K=4\pi \mathrm{\Sigma }`$.

The prefactor $`A`$ for tunneling in $`(1+1)`$ dimensions is given by $`A=BD`$, where $`D`$ is a determinantal factor defined in Ref.  and the factor $`B`$ comes from including the zero modes contribution to the determinant.

The exact calculation of this prefactor is a complicated task and we do not attempt it here. Our main interest here is in the dependence of the prefactor on the worlsheet current. In the Appendix we show that we should expect this prefactor to scale as,

$$AB\delta ^2\left(\frac{L}{2\pi N\delta }\right)^2.$$
(34)
Analytic estimates and numerical calculations show  that, if the fields $`\sigma `$ and $`\varphi `$ are weakly coupled to each other and to themselves, then $`\mathrm{\Sigma }1`$ in a large portion of the parameter space of the model (1). In this case, $`K10`$, and it follows from Eqs. (33) and (34) that the current decay is strongly suppressed, as long as the condensate winding does not get very large, namely, $`NL/2\pi \delta `$.

We finally discuss the validity of the approximations we made in deriving the tunneling rate (33). First, we assumed that the vortex cores give a negligible contribution to the instanton action, compared to the contribution coming from the region

$$\delta |𝐱𝐱_\pm |a.$$
(35)
This is justified if the logarithm in Eq. (31) is large,

$$\mathrm{ln}(L/2\pi N\delta )>1,$$
(36)
that is, when the current is very small. We expect the vortex radius $`\delta `$ to be comparable to the thickness of the condensate; then Eq. (36) says that the wavelength of the condensate should be large compared to its thickness.

We also assumed that the instanton solution can be factorized as in (4), with $`\theta `$ independent of the coordinates $`(x,y)`$ in the plane transverse to the string. This approximation is likely to break down in the vicinity of the vortex cores, where $`\theta `$ varies on a length scale comparable to the string thickness, but we expect it to be valid in the range (35) which, assuming (36), gives the dominant contribution to the action. Eq. (36) is thus the only condition we need to impose to justify the approximations made above.

## IV Charged current instanton.

We would like now to extend the previous analysis to the case where the condensate is coupled to the electromagnetic field $`A_\mu `$. The simplest way to achieve this is to gauge the $`U(1)_{global}`$ symmetry we have been discussing. In this case the action for the model takes the form,

$$S_{el}=d^4x\overline{D_\mu \sigma }D^\mu \sigma \frac{1}{4}f_{\mu \nu }f^{\mu \nu }+\overline{𝒟_a\varphi }𝒟^a\varphi V(\sigma ,\varphi )\frac{1}{4}F_{\mu \nu }F^{\mu \nu }$$
(37)
where

$`F_{\mu \nu }`$ $`=`$ $`_\mu A_\nu _\nu A_\mu ,`$ (38)
$`𝒟_\mu \varphi `$ $`=`$ $`(_\mu ieA_\mu )\varphi .`$ (39)
As in the previous case, we are interested in the quantum mechanical decay of a spacelike current of the form

$$\varphi =\varphi (x,y)\mathrm{exp}[i\theta (\tau ,z)].$$
(40)
In terms of the angular variable the effective action becomes

$$S_{el}=\mathrm{\Sigma }_{\text{eff}}𝑑z𝑑\tau (_a\theta +eA_a)^2\frac{1}{4}d^4xF_{\mu \nu }F^{\mu \nu }.$$
(41)
where $`\mathrm{\Sigma }_{\text{eff}}`$ can be written in terms of the parameters in the original $`U(1)_{local}\times U(1)_{local}`$ Lagrangian. This expression is correct if the scalar field condensate and the vector field remain approximately constant in the transverse direction to the string. This is ensured if the penetration depth of the vector field in the core of the string is large compared to the condensate thickness.

The important difference in this case is that the electromagnetic field is not confined to the worldsheet, so our instanton is somewhat more complicated. Adding the electromagnetic field makes the Euclidean equations of motion for the worldsheet look like those of a two dimensional superconducting sheet embedded in a four dimensional Euclidean space. Our initial condition at $`\tau =\mathrm{}`$ is a sheet carrying a constant planar current with its accompanying magnetic field. The instanton includes a vortex-antivortex pair in equilibrium under the action of the Magnus force, the electromagnetic force due to the long-range magnetic field of the vortices, and (to a smaller extent) the force due to the scalar field distribution on the sheet.

When this work was in progress we learned that this type of instanton has been discussed by Duan in the context of a real superconducting wire. The instanton solution cannot be found analytically in this case, but remarkably the instanton action, (excluding the contribution from the vortex cores) can still be found exactly,

$$B=\frac{4\pi \mathrm{\Sigma }}{\alpha }\left[\mathrm{ln}\left[1+\alpha \mathrm{ln}\left(\frac{2a}{\delta }\right)\right]a\left(\frac{4\pi N\alpha }{L}\right)\right]$$
(42)
where $`\alpha =\frac{\mathrm{\Sigma }e^2}{\pi }`$, and $`N`$, $`L`$, and $`\delta `$ have the same meaning as before. Note that this reduces to our Eq. (31) in the neutral limit $`e0`$. The action is maximized at

$$a_0\left(\frac{L}{4\pi N}\right)\frac{1}{1+\alpha \mathrm{ln}\left(\frac{L}{2\pi N\delta }\right)}.$$
(43)
and finally it becomes

$$B\frac{4\pi \mathrm{\Sigma }}{\alpha }\left[\mathrm{ln}\left[1+\alpha \mathrm{ln}\left(\frac{L}{2\pi N\delta }\right)\right]\frac{\alpha }{1+\alpha \mathrm{ln}\left(\frac{L}{2\pi N\delta }\right)}\right]$$
(44)
As we discussed at the end of Section III, we are interested in the range of parameters where $`\mathrm{\Sigma }>100`$ and $`\mathrm{ln}(L/2\pi N\delta )>1`$. With $`e^210^2`$ this gives $`\alpha >1`$ and Eq. (44) reduces to

$$B\frac{4\pi ^2}{e^2}\mathrm{ln}\left[\alpha \mathrm{ln}\left(\frac{L}{2\pi N\delta }\right)\right]$$
(45)
Following the discussion for the neutral case we see that the tunneling rate for the charged case is

$$\mathrm{\Gamma }A\left(\alpha \mathrm{ln}\left(\frac{L}{2\pi N\delta }\right)\right)^{(2\pi /e)^2}$$
(46)
which clearly indicates that this process is very much supressed in the electromagnetic case as well.

## V Discussion

Our result (33) was derived in the frame where the global charge density vanishes, $`J_0=0`$. A generalization for an arbitrary frame is easily obtained if we notice that the rate $`\mathrm{\Gamma }`$ is Lorentz invariant, so we only need to express it in terms of the invariant quantity $`J_aJ^a`$. With the aid of Eqs. (33), (34) and assuming (36), we obtain

$$\mathrm{\Gamma }K\delta ^2X^{(K2)/2}\mathrm{ln}(1/X),$$
(47)
where

$$X=|J_aJ^a\delta ^2|.$$
(48)
It is interesting to note that the rate (47) vanishes in the chiral limit $`J_aJ^a=0`$. This is not surprising: since $`\mathrm{\Gamma }`$ is an invariant, it can only depend on the invariant combination $`J_aJ^a`$ and should therefore vanish for a chiral current. We expect this conclusion to be valid for strings with charged as well as neutral currents. The stability of chiral vortons with respect to this type of decay reinforces the argument that cosmic string theories able to produce chiral vortons will have serious constraints from their overproduction in the course of the history of the string network.

For a timelike current, $`J_aJ^a>0`$, instantons of the type we discussed here do not exist. In this case, one can go to the Lorentz frame where the current vanishes and there is only a non-zero global charge density. Charge is conserved on the worldsheet and cannot therefore be reduced by worldsheet instantons or any other processes operating on the string worldsheet. It is possible that when the charge density gets high enough, charges can escape from the string by tunneling into the external space. We shall not consider such tuneling processes in the present paper.

## VI Acknowledgments

We would like to thank Jaume Garriga for very helpful conversations. This work was supported in part by funding provided by the National Science Foundation. J.J.B-P was supported in part by the Relativity Group PPARC Rolling Grant.

## VII Appendix

In this Appendix we will estimate the dependence of the decay rate prefactor on the uniform current on the worldsheet. The usual procedure for this type of calculation involves the computation of the imaginary part of the energy coming from the contribution around the one-instanton solution of the Euclidean equations of motion. The decay rate per unit volume per unit time is given by,

$$\mathrm{\Gamma }=\underset{i=1,2}{}\left(\frac{N_i}{2\pi }\right)^{1/2}\frac{\left|Det^{1/2}S_{}^{\prime \prime }{}_{E}{}^{}\right|}{Det^{1/2}S_{}^{\prime \prime }{}_{E}{}^{(0)}}e^{(S_ES_E^{(0)})}$$
(49)
where $`S_E^{\prime \prime }`$ and $`S_{}^{\prime \prime }{}_{E}{}^{0}`$ denote the second order perturbation expansion of the instanton and the background Euclidean action respectively. In our notation we have

$$S_ES_E^{(0)}=B.$$
(50)
Furthermore, the first term, proportional to $`N_i`$, is included to take into account the normalization of the zero modes, and in general it can be shown to be proportional to the instanton action $`B`$. Consequently we do not have to consider those in the following calculation. This is signaled by the expression $`Det^{}`$ which denotes the determinant of the fluctuations once the zero modes have been excluded. In our system the zero modes correspond to the collective motion of the vortex pair in the $`(\tau ,z)`$ plane. It is clear that translating them together does not alter the action, so this kind of fluctuations must correspond to zero modes.

Among the modes included in $`Det^{}`$ is the negative eigenvalue mode of the fluctuations around the instanton solution. As we discussed earlier, this negative mode corresponds to the fluctuation that changes the distance between the vortices in the direction perpendicular to the current flow.

There is also another type of fluctuation which corresponds to a relative motion of the vortices in the direction parallel to the current. This is a positive mode.

If the current is very small, in other words, if the distance between the two vortices is very large compared to the size of the vortex core, it could be expected that the only fluctuations which would not cancel with the ones in the background are the fluctuations associated with motions of the vortices. In this case the effective one-loop contribution to the prefactor can be written as the following integral,

$$\mathrm{\Gamma }Be^B\frac{dz_rd\tau _r}{\mu ^4}e^{\beta (\frac{z_r}{a_o})^2}e^{\beta (\frac{\tau _ra_o}{a_o})^2},$$
(51)
where the factor of $`B`$ is the usual normalization factor for the two zero modes in 1+1 dimensions. This expression singles out the contribution from the modes that are related to the position of the vortices by defining $`z_r`$ and $`\tau _r`$ as the relative distance between the vortices on the worlsheet coordinates. Furthermore, $`a_0`$ denotes the distance along the Euclidean time direction that maximizes the instanton action, $`\beta =4\pi \mathrm{\Sigma }`$ and $`\mu `$ is a normalization factor of order $`|\varphi |\delta ^1`$ summarizing the result of doing the small fluctuation integral around the one vortex configuration.

Performing the gaussian integrals we arrive at

$$\mathrm{\Gamma }Be^B\delta ^2\left(\frac{L}{2\pi N\delta }\right)^2$$
(52)


format: text

text:
# Geometrical explanation and scaling of dynamical heterogeneities in glass forming systems

## Abstract

We show how dynamical heterogeneities in glass forming systems emerge as a consequence of the existence of dynamical constraints, and we offer an interpretation of the glass transition as an entropy crisis in trajectory space (space–time) rather than in configuration space. To illustrate our general ideas, we analyze the one dimensional ($`d=1`$) Fredrickson–Andersen and East models. Dynamics of such dynamically constrained systems are shown to be isomorphic to the statics of $`d+1`$ dimensional dense mixtures of polydisperse non–interpenetrating domains. The domains coincide with arrested regions in trajectory space.

PACS numbers: 64.70.Pf, 75.10.Hk, 05.70.Ln

A glass forming system, like a supercooled liquid, exhibits a precipitous onset of slowness. As temperature is decreased in these systems, typically in a range of a few decades, relaxation times and viscosities increase by several orders of magnitude, eventually surpassing experimentally accessible times. For practical purposes, these systems effectively freeze at the glass transition temperature $`T_g`$. For reviews see . Interestingly, this dynamical arrest carries no evident static structural signature of growing length scales. Rather, experiments and simulations show that supercooled liquids are dynamically heterogeneous . Molecules in one region of the liquid translate or rotate several orders of magnitude faster or slower than those in a neighboring region. The spatial extent of these dynamical heterogeneities is mesoscopic, and the time scale of the slowest domains increases with decreasing temperature at least as fast as the relaxation time of the system. Such structural behavior seems beyond description with homogeneous methods like mode coupling  and mean field theories , and it is widely neglected in analytical treatments (see, however, ). Nevertheless, we show here that for a broadly applicable mechanism of dynamical arrest, these heterogeneities are intrinsic to the nature of glass forming systems.

Our central result is that dynamical heterogeneities are a manifestation of the existence of nontrivial structure in the trajectories of glassy systems. This structure associated with dynamics is independent of any specific static properties. Instead, the nontrivial dynamical structure is a consequence of local dynamical rules that significantly restrict the size of accessible trajectory space. For example, consider a highly compressed (or supercooled) glass former. Atoms in most regions of space are jammed, making mobility possible in only a relatively low fraction of spatial regions. These rare regions are those that are already unjammed, or those that may be close in space to an unjammed region. In the evolution of such a system, one therefore expects a clustering of mobile regions and thus a mesoscopic demixing of mobile and static regions. Macroscopic demixing is not expected because dynamics should conserve a canonical distribution. This picture of the origin of dynamic heterogeneity is in accord with the idea that glassiness is not necessarily a consequence of either disorder or frustration in the static interactions but of the existence of effective constraints on the dynamics of the system .

The simplest microscopic models that illustrate this view are the Fredrickson–Andersen (FA) model  and the East model . They consist of a chain of Ising spins $`\sigma _i=\pm 1`$ $`(i=1\mathrm{}N)`$, with trivial Hamiltonian $`H=_i\sigma _i`$, and single spin flip dynamics subject to local kinetic constraints. In the FA model, a spin can flip if either of its nearest neighbors is in the up state . In the East model, a spin can flip only if its nearest neighbor to the right is up . The equilibrium behavior of both models is that of an uncorrelated spin system, with $`c=1/\left(1+e^{1/T}\right)`$ being the average concentration of up spins at temperature $`T`$. The competition between decreasing the energy and the need for facilitating spins leads to a glassy slow down at low temperatures. The relaxation times go as $`\tau e^{3/T}`$ in the FA model, and as $`\tau e^{1/(T^2\mathrm{ln}2)}`$ in the East model. These models therefore correspond to strong and a fragile glass formers, respectively. See  for a short review.

Let us first establish that both the FA and East models display dynamical heterogeneities. Consider the coarse grained spins $`s_i(t;\mathrm{\Delta }t)(\mathrm{\Delta }t)^1_0^{\mathrm{\Delta }t}𝑑t^{}\sigma _i(t+t^{})`$. Fast and slow spins in the time window of width $`\mathrm{\Delta }t`$ at time $`t`$ will correspond to low and high values of $`s_i^2(t;\mathrm{\Delta }t)`$, respectively. The spatial distribution of $`s_i^2(t;\mathrm{\Delta }t)`$ will determine the extent to which the system is dynamically heterogeneous. The heterogeneity made evident with this field depends on the coarse graining time. For $`\mathrm{\Delta }t`$ very short, only the trivial uncorrelated static structure is probed. The same is true for $`\mathrm{\Delta }t`$ much larger than the relaxation time where ergodicity is restored. For intermediate values of $`\mathrm{\Delta }t`$, however, we may expect to see spatial structure in the $`s_i^2(t;\mathrm{\Delta }t)`$. In Fig. 1 we show that this is indeed the case by plotting the structure factor, $`S(k)`$, for the field $`s_i^2(t;\mathrm{\Delta }t)`$ calculated from simulations of the FA model (left panel) and East model (right panel). We have used coarse graining times of about a fifth of the relaxation time, and simulations were performed using continuous time Monte Carlo  for systems of $`10^5`$ spins averaged over $`10^3`$ samples. For both models, the existence of a correlation length $`\xi _{\mathrm{\Delta }t}`$ is clear. Moreover, $`\xi _{\mathrm{\Delta }t}`$ grows, albeit slowly, with decreasing temperature and thus increasing relaxation time. For the FA model the structure factor decays approximately as $`S(k)k^2`$ for large $`k`$. In the East model, at intermediate $`k`$, the structure factor goes as $`S(k)k^{\mathrm{ln}3/\mathrm{ln}2}`$. For larger momenta, the structure factor is oscillatory, a feature which becomes more pronounced with decreasing temperature.

In order to understand these behaviors, we study the generating functional for the trajectories of these systems. Consider a trajectory $`\{\stackrel{}{\sigma }_t\}(\stackrel{}{\sigma }_0,\stackrel{}{\sigma }_1,\stackrel{}{\sigma }_2,\mathrm{},\stackrel{}{\sigma }_𝒯)`$, where $`\stackrel{}{\sigma }_t(\sigma _{1t},\sigma _{2t},\mathrm{},\sigma _{Nt})`$ denotes the configuration of the $`N`$ spin system at time $`t`$. Each trajectory has a probability $`P(\{\stackrel{}{\sigma }_t\})`$, so that the total probability of going from $`\stackrel{}{\sigma }_0`$ to $`\stackrel{}{\sigma }_𝒯`$ is the weighted sum of such trajectories, $`Z_{\stackrel{}{\sigma }_0,\stackrel{}{\sigma }_𝒯}_{\{\stackrel{}{\sigma }_t\}}P(\{\stackrel{}{\sigma }_t\})`$. This partition sum, $`Z_{\stackrel{}{\sigma }_0,\stackrel{}{\sigma }_𝒯}`$, can be written as

$$Z_{\stackrel{}{\sigma }_0,\stackrel{}{\sigma }_𝒯}=\underset{\{\stackrel{}{\sigma }_t\}}{}e^{S_0[\{\stackrel{}{\sigma }_t\}]}\left[\underset{it}{}\delta (h_{it})\right]e^{\mathrm{\Delta }S[\{\stackrel{}{\sigma }_t\}]}.$$
(1)
The first factor in the summand of (1) corresponds to the probability of a trajectory in the absence of dynamical constraints. Its action is just that of $`N`$ non–interacting ferromagnetic Ising chains in uniform fields, $`S_0(\{\stackrel{}{\sigma }_t\})=_{i=1}^N_{t=0}^{𝒯1}(J_+\sigma _{it}\sigma _{it+1}`$ $`+`$ $`J_+\sigma _{it}`$ $`+`$ $`J_+\sigma _{it+1}`$ $`+`$ $`J_{+++})`$, where $`J_{\mu \nu \epsilon }\{\mathrm{ln}[1\gamma \delta t(1c)]`$ $`+`$ $`\mu \mathrm{ln}(1\gamma \delta tc)`$ $`+`$ $`\nu \mathrm{ln}[\gamma \delta t(1c)]`$ $`+`$ $`\epsilon \mathrm{ln}(\gamma \delta tc)\}/4`$, $`\gamma `$ is a microscopic rate, and $`\delta t`$ the time step. The second factor in the summand of (1) embodies the kinetic constraint. Only trajectories which satisfy the condition $`h_{it}=0`$ at all space–time points are allowed. If $`h_{it}`$ is chosen to be a positive number when the constraint is not satisfied, then $`_{it}\delta (h_{it})=\delta ()=e^\lambda `$, where $`_{it}h_{it}`$ and $`\lambda \mathrm{}`$. The infinite coupling constant means that the space of allowed trajectories is that of the ground states of $``$. The respective forms of $``$ for the FA and East models are

$`_{\mathrm{FA}}`$ $`=`$ $`{\displaystyle \underset{it}{}}\left(1\sigma _{it}\sigma _{it+1}\right)\left(1\sigma _{i+1t}\right)\left(1\sigma _{i1t}\right)/8,`$ (2)
$`_{\mathrm{East}}`$ $`=`$ $`{\displaystyle \underset{it}{}}\left(1\sigma _{it}\sigma _{it+1}\right)\left(1\sigma _{i+1t}\right)/4,`$ (3)
so that the $`i`$–th term is zero unless the spin $`i`$ flips and the facilitating spins are in the down state. The third factor in the summand of Eq. (1) accounts for preservation of norm and detailed balance. Trajectories which are allowed by the constraints have different probabilities than they would in the unconstrained case. In particular,

$`\mathrm{\Delta }S_{\mathrm{FA}}`$ $`=`$ $`{\displaystyle \underset{it}{}}\left(1\sigma _{i+1t}\right)\left(1\sigma _{i1t}\right)\left(\stackrel{~}{J}_{}\sigma _{it}+\stackrel{~}{J}_+\right)/4,`$ (4)
$`\mathrm{\Delta }S_{\mathrm{East}}`$ $`=`$ $`{\displaystyle \underset{it}{}}\left(1\sigma _{i+1t}\right)\left(\stackrel{~}{J}_{}\sigma _{it}+\stackrel{~}{J}_+\right)/2,`$ (5)
where $`\stackrel{~}{J}_\pm J_{\pm 00}+\mathrm{ln}[1\gamma \delta t(1c)e^\lambda ]/4\pm \mathrm{ln}(1\gamma \delta tce^\lambda )/4`$. These expressions describe a proper dynamics for any $`\lambda `$. When $`\lambda =0`$, the expressions reduce to the trivial model of independent spins with unconstrained dynamics. When $`\lambda \mathrm{}`$, they correspond to the FA and East models.

$``$ and $`\mathrm{\Delta }S`$ introduce competing spatial interactions in the space of trajectories. Those of $``$ are strong, and they are ferromagnetic in the sense that they favor the clustering of like spins. In contrast, those of $`\mathrm{\Delta }S`$ are weak and antiferromagnetic. Moreover, the scaling with distance of interactions in $``$ is different than that for those in $`\mathrm{\Delta }S`$, so that nontrivial structure in the space of trajectories can be expected. In Fig. 2, we show samples of equilibrium trajectories for the unconstrained case (top), the FA (bottom left) and East models (bottom right), at $`T=1.0`$. The difference between the constrained and unconstrained dynamics is striking. The trajectories in both the FA and East models display an extensive number of domains of down spins. These domains are the origin of the dynamical heterogeneities. Spins within these domains do not change, so when coarse grained in time, they correspond to slow regions. To the extent that lowering temperature and thus energy in an atomic system coincides with decreasing facilitating spin concentration $`c`$, the structures observed in Fig. 2 coincide with the correlation observed in simulation of supercooled liquids . In particular, slow and fast dynamical heterogeneities correlate with regions of low and high energy, respectively. Since each column of the pictures in Fig. 2 is an equilibrium configuration of the noninteracting $`H`$, the structure seen in the trajectories is purely dynamical.

Down spins must form closed domains in a trajectory as a consequence of the local and causal nature of the dynamical constraints. This fact is apparent from the illustration in Fig. 3. A spin is able to flip only if an appropriate neighboring spin is in the up state at the same time. It therefore will have flipped up previously and/or it will flip down later. As such, a well defined closed boundary must exist between regions of up and down spins. These boundaries are formed out of segments with shapes like those depicted on the top of Fig. 3, with the possibility of different slopes. In the FA model, the dynamical constraint is spatially symmetric, and four possible kinds of boundary segments are allowed, with the restriction that all the segments of the first (second) kind must be below those of the third (fourth) kind. It follows that spin down domains must form semi convex polygons like that pictured at the bottom left of Fig. 3. Namely, any spatial line cuts the boundary only twice. In the case of the East model, only the third and fourth boundary segments are allowed, which restricts the domains to shapes like that pictured at the bottom right of Fig. 3. Since slow dynamical heterogeneities correspond to spatial projections of the compact space–time down spin domains, they are necessarily compact, while the converse is true for fast regions. This observation is in agreement with what is found with computer simulations of atomistic models .

The geometrical construction described above implies that spin down domains cannot penetrate each other. Therefore, the set of trajectories maps to the configuration space of a two dimensional mixture of polydisperse non–interpenetrating objects of all the possible shapes allowed by the dynamical constraints. This observation motivates a description of the dynamics in terms of $`\rho (l,t)`$, the density of domains of typical spatial size (height) $`l`$ and extension in time (length) $`t`$. The partition function for trajectories is then $`ZD\rho (l,t)\mathrm{exp}\left(\mathrm{\Omega }[\rho (l,t)]\right)`$, where $`\mathrm{\Omega }\left[\rho (l,t)\right]`$ is the free energy functional of this density. Whatever estimate is used for $`\mathrm{\Omega }\left[\rho (l,t)\right]`$, it is crucial that it enforces the constraint

$$𝑑t\rho (l,t)=\rho _{\mathrm{eq}}(l)=ce^{cl}.$$
(6)
This condition ensures that any spatial cut of the equilibrium trajectories is an equilibrium configuration. For a density functional theory, it provides an approximation to the effects of the antiferromagnetic interactions of $`\mathrm{\Delta }S`$.

In general, $`\rho (l,t)=\rho (l|t)p(t)`$, where $`\rho (l|t)`$ is the probability density of $`l`$ conditioned on $`t`$, and $`p(t)`$ is the probability density of domains of length $`t`$. In the case of the FA model, a domain is bounded by the random walks of two up spins between successive encounters, and $`\rho (l|t)`$ can be obtained by standard arguments

$$\rho (l|t)l^2(Dt)^{3/2}\mathrm{exp}(l^2/Dt)$$
(7)
with $`Dc`$. If the domains were isolated, $`p(t)`$ would be given by the probability of first return of a random walker, which goes as $`p(t)t^{3/2}`$ for large $`t`$, leading to $`t\mathrm{}`$, and to the formation of unbounded domains. This result, however, corresponds to only the ferromagnetic part of the interactions in trajectory space. It is frustrated by the condition (6). By combining (7) with (6), we have instead $`p(t)(\tau t)^{1/2}e^{\sqrt{t/\tau }}`$, with $`\tau D^1c^2e^{3/T}`$. This result for $`\tau `$ is precisely the relaxation time in the FA model. The forms of $`\rho (l|t)`$ and $`p(t)`$ are verified in simulations of the FA model . We also see that in terms of the scaling variables $`l^{}=cl`$ and $`t^{}=t/\tau `$, the density $`\rho (l^{},t^{})`$ is independent of temperature. An illustration of these scaling relations is given in Fig. 4. The typical height of domains grows with the cubic root of the relaxation time, thus explaining the slow increase of the correlation length $`\xi _{\mathrm{\Delta }t}`$ with decreasing temperature in the structure factors of Fig. 1. Slow growth and domains being mesoscopic rather than macroscopic are consistent with experimental observations . Moreover, since the trajectories are extensive in interfaces, the form of the structure factor for large $`k`$ corresponds to a spatial projection of Porod’s law $`S(k)k^{(d+1)}`$ .

In the case of the East model, typical distances and times are related at low temperatures by a $`T`$ dependent dynamic exponent $`z(T)1/(T\mathrm{ln}2)`$ . We therefore expect the conditional probabilities to be $`\rho (l|t)\mathrm{exp}[l^{z(T)}/t]`$. The condition (6) then leads to a stretched exponential form for the persistence function, $`_t^{\mathrm{}}𝑑t^{}p(t^{})\mathrm{exp}[(t/\tau )^{\beta (T)}]`$, where $`\beta (T)1/z(T)`$ and $`\tau c^{z(T)}`$, in accordance with previous results . In this case, $`\xi _{\mathrm{\Delta }t}`$ increases even more slowly with decreasing $`T`$ than in the FA model. The geometry of the domains also explains the differences in the FA and East structure factors shown in Fig. 1. In contrast to the FA case, where neighboring domains can be compressed to the point of contact, the boundary between domains in the East model cannot be formed just by a single line of up spins in space-time. Instead, to decrease energy and thus increase probability, the boundary is wet by smaller domains, ideally like in Fig. 5. This structure has fractal dimension $`d_f=\mathrm{ln}3/\mathrm{ln}2`$, and gives rise to the behavior $`S(k)k^{d_f}`$ at intermediate momenta. At larger $`k`$, $`S(k)`$ probes the granular structure of the boundary, and thus exhibits the oscillatory behavior shown in Fig.1 .

The ideas presented here in detail for the FA and East models generalize to higher dimensions provided the dynamical constraints are causal and local, so that excitations favor the creation of neighboring excitations. In general, dynamical constraints will appear only below a crossover temperature $`T_x`$ or above a corresponding packing fraction. This can be incorporated into our description by allowing the coupling $`\lambda `$ to be a function of temperature or packing fraction, increasing with decreasing $`T`$ or increasing packing. The temperature $`T_x`$ will then coincide with the so-called “landscape” temperature.

The domains that thus appear for $`T<T_x`$ and grow with further decrease in $`T`$ can lead to a broken symmetry. Consider, for example, the partition function $`Z_T(q,𝒯)`$ for trajectories between configurational fluctuations with overlap $`q`$ at time difference $`𝒯`$, $`Z_T(q,𝒯)_{\{\stackrel{}{\sigma }_0,\stackrel{}{\sigma }_𝒯\}}\widehat{\rho }(\stackrel{}{\sigma }_0)Z_{\stackrel{}{\sigma }_0,\stackrel{}{\sigma }_𝒯}\delta \left(qN^1_i\delta \sigma _{i0}\delta \sigma _{i𝒯}\right)`$, where $`\delta \sigma _{it}\sigma _{it}\sigma `$, and $`\widehat{\rho }`$ is Boltzmann distribution. Since $`Z_T(q,𝒯)`$ is proportional to the number of such trajectories, it is natural to write it as $`Z_T(q,𝒯)\mathrm{exp}[N\omega _T(q,𝒯)]`$, where $`\omega _T(q,𝒯)`$ is the entropy density in trajectory space. It is approximately the entropy of mixing of the slow domains in trajectory space. At values of $`T`$ for which the temporal extension of these domains is much smaller than $`𝒯`$, $`\omega _T(q,𝒯)`$ is extensive in time, and $`Z_T(q,𝒯)`$ is peaked at $`q=0`$. In this case, correlation functions are exponential.

Non-zero overlap is probable, i.e., $`Z_T(q,𝒯)`$ is peaked at finite $`q`$, only when $`𝒯`$ is comparable to or smaller than the length of typical slow domains. As temperature is decreased, the size of these domains increases, constricting the available trajectory space. At a low enough temperature, we may therefore expect $`\omega _T(q,𝒯)`$ to become sub-extensive, leading to a probable finite $`q`$ throughout a relatively large range of $`𝒯`$. This change in behavior of $`\omega _T(q,𝒯)`$ signals a dramatic change in the number of available trajectories and a corresponding onset of dynamical arrest. In this picture, therefore, the glass transition coincides with an entropy crisis in trajectory space, rather than in configuration space.

This work was supported by the National Science Foundation, the Glasstone Fund, and Merton College, Oxford.