text:
# Field theoretic calculation of renormalized-viscosity, renormalized-resistivity, and energy fluxes of magnetohydrodynamic turbulence

## Abstract

A self-consistent renormalization (RG) scheme has been applied to nonhelical magnetohydrodynamic turbulence with zero cross helicity. Kolmogorov’s $`5/3`$ powerlaw has been shown to be a consistent solution for $`dd_c2.2`$. For Kolmogorov’s solution, both renormalized viscosity and resistivity are positive for the whole range of parameters. Various cascade rate and Kolmogorov’s constant for MHD turbulence have been calculated by solving the flux equation to the first order in perturbation series. We find that the magnetic energy cascades forward. The Kolmogorov’s constant for $`d=3`$ does not vary significantly with $`r_A`$ and is found to be close to the constant for fluid turbulence.

PACS numbers: 47.27.Gs, 52.35.Ra, 11.10.Gh

The statistical theory of magnetohydrodynamic (MHD) turbulence is one of the important problems of current research. The quantities of interests in this area are energy spectrum, cascade rates, intermittency exponents etc. In this letter we analytically compute the renormalized-viscosity, renormalized-resistivity, and cascade rates using the field-theoretic techniques.

The incompressible MHD equation in Fourier space is given by

$`\left(i\omega +\nu k^2\right)u_i(\widehat{k})`$ $`=`$ $`{\displaystyle \frac{i}{2}}P_{ijm}^+(𝐤){\displaystyle 𝑑\widehat{p}[u_j(\widehat{p})u_m(\widehat{k}\widehat{p})b_j(\widehat{p})b_m(\widehat{k}\widehat{p})]}`$ (1)
$`\left(i\omega +\lambda k^2\right)b_i(\widehat{k})`$ $`=`$ $`iP_{ijm}^{}(𝐤){\displaystyle 𝑑\widehat{p}[u_j(\widehat{p})b_m(\widehat{k}\widehat{p})]}`$ (2)
$`k_iu_i(𝐤)`$ $`=`$ $`0`$ (3)
$`k_ib_i(𝐤)`$ $`=`$ $`0`$ (4)
where $`𝐮`$ and $`𝐛`$ are the velocity and magnetic field fluctuations respectively, $`\nu `$ and $`\lambda `$ are the viscosity and the resistivity respectively, and $`d`$ is the space dimension. Also,

$`P_{ijm}^+(𝐤)`$ $`=`$ $`k_jP_{im}(𝐤)+k_mP_{ij}(𝐤);`$ (5)
$`P_{im}(𝐤)`$ $`=`$ $`\delta _{im}{\displaystyle \frac{k_ik_m}{k^2}};`$ (6)
$`P_{ijm}^{}(𝐤)`$ $`=`$ $`k_j\delta _{im}k_m\delta _{ij};`$ (7)
$`\widehat{k}=(𝐤,\omega );`$ $`d\widehat{p}=d𝐩d\omega /(2\pi )^{d+1}.`$ (8)
The energy spectra for MHD, $`E^u(k)`$ and $`E^b(k)`$, are still under debate. Kraichnan  and Irosnikov  first gave phenomenology of steady-state, homogeneous, and isotropic MHD turbulence, and proposed that the spectra is proportional to $`k^{3/2}`$. Later Marsch , Matthaeus and Zhou , and Zhou and Matthaeus  proposed an alternate phenomenology in which the energy spectra are proportional to $`k^{5/3}`$, similar to Kolmogorov’s spectrum for fluid turbulence. Current numerical  and theoretical  work support Kolmogorov-like phenomenology for MHD turbulence. In the present paper we show that Kolmogorov’s spectrum ($`k^{5/3}`$) is a consistent solution of renormalization group (RG) equation of MHD turbulence.

Forster et al., DeDominicis and Martin, Fournier and Frisch, Yakhot and Orszag  applied RG technique to fluid turbulence. They considered external forcing and calculated renormalized parameters: viscosity, noise coefficient, and vertex. McComb and his coworkers  instead applied a self-consistent RG procedure; here the energy spectrum was assumed to be Kolmogorov’s powerlaw, and the renormalized viscosity was computed iteratively. For MHD turbulence, Fournier et al., Camargo and Tasso, and Liang and Diamond  employed RG technique in the similar lines as Forster et al. . In this letter we will apply McComb’s self-consistent technique to MHD turbulence. Earlier Verma  had done a self-consistent calculation and showed that the mean magnetic field gets renormalized, and the Kolmogorov’s powerlaw is a consistent solution of MHD RG equation. In this letter we will carry out the renormalization of viscosity and resistivity.

For simplicity of the calculation we assume that the mean magnetic field is absent. This allows us to assume the turbulence to be isotropic to a reasonable approximation. In presence of mean magnetic field, turbulence become anisotropic; this issue has been studied by Sridhar and Goldreich  and Goldreich and Sridhar . In addition to the above assumption, we also take cross helicity $`(2𝐮𝐛)`$, magnetic helicity $`(𝐚𝐛)`$, and kinetic helicity $`(𝐮\omega )`$ to be zero, where $`𝐚`$ is magnetic vector potential, and $`\omega `$ is the vorticity.

In our RG procedure the wavenumber range $`(k_N,k_0)`$ is divided logarithmically into $`N`$ shells. We carry out the elimination of the first shell $`k^>=(k_1,k_0)`$ and obtain the modified MHD equation for $`k^<=(k_N,k_1)`$. This process is continued for all the shells. The shell elimination is performed by ensemble averaging over $`k^>`$ modes . We assume that $`u_i^>(\widehat{k})`$, and $`b_i^>(\widehat{k})`$ have gaussian distributions with zero mean, while $`u_i^<(\widehat{k})`$ and $`b_i^<(\widehat{k})`$ are unaffected by the averaging process. In addition we take

$`u_i^>(\widehat{p})u_j^>(\widehat{q})`$ $`=`$ $`P_{ij}(𝐩)C^{uu}(\widehat{p})\delta (\widehat{p}+\widehat{q})`$ (9)
$`b_i^>(\widehat{p})b_j^>(\widehat{q})`$ $`=`$ $`P_{ij}(𝐩)C^{bb}(\widehat{p})\delta (\widehat{p}+\widehat{q})`$ (10)
Let us denote $`\nu _{(n)}`$ and $`\lambda _{(n)}`$ as the viscosity and resistivity after the elimination of the $`n`$ shell. To first order of perturbation, we obtain

$`\left(i\omega +\nu _{(n)}k^2+\delta \nu _{(n)}k^2\right)u_i^<(\widehat{k})`$ $`=`$ $`{\displaystyle \frac{i}{2}}P_{ijm}^+(𝐤){\displaystyle }d\widehat{p}[u_j^<(\widehat{p})u_m^<(\widehat{k}\widehat{p})`$ (12)
$`b_j^<(\widehat{p})b_m^<(\widehat{k}\widehat{p})]`$
$`\left(i\omega +\lambda _{(n)}k^2+\delta \lambda _{(n)}k^2\right)b_i^<(\widehat{k})`$ $`=`$ $`iP_{ijm}^{}(𝐤){\displaystyle 𝑑\widehat{p}[u_j^<(\widehat{p})b_m^<(\widehat{k}\widehat{p})]}`$ (13)
where

$`\delta \nu _{(n)}(k)`$ $`=`$ $`{\displaystyle \frac{1}{(d1)k^2}}{\displaystyle _{\widehat{p}+\widehat{q}=\widehat{k}}^\mathrm{\Delta }}{\displaystyle \frac{d𝐩}{(2\pi )^d}}[S(k,p,q){\displaystyle \frac{C^{uu}(q)}{\nu _{(n)}(p)p^2+\nu _{(n)}(q)q^2}}`$ (15)
$`S_6(k,p,q){\displaystyle \frac{C^{bb}(q)}{\lambda _{(n)}(p)p^2+\lambda _{(n)}(q)q^2}}]`$
$`\delta \lambda _{(n)}(k)`$ $`=`$ $`{\displaystyle \frac{1}{(d1)k^2}}{\displaystyle _{\widehat{p}+\widehat{q}=\widehat{k}}^\mathrm{\Delta }}{\displaystyle \frac{d𝐩}{(2\pi )^d}}[S_8(k,p,q){\displaystyle \frac{C^{bb}(q)}{\nu _{(n)}(p)p^2+\lambda _{(n)}(q)q^2}}`$ (17)
$`+S_9(k,p,q){\displaystyle \frac{C^{uu}(q)}{\lambda _{(n)}(p)p^2+\nu _{(n)}(q)q^2}}]`$
with $`S_i(k,p,q)`$s as functions of $`k,p`$, and $`q`$. Hence, after the elimination of the $`(n+1)`$th shell, the effective viscosity and resistivity will be $`(\nu ,\lambda )_{(n+1)}(k)=(\nu ,\lambda )_{(n)}(k)+\delta (\nu ,\lambda )_{(n)}(k)`$.

We solve the above equations iteratively. To simplify, we substitute $`C(k)`$ in the Eqs. (15, 17) by one dimensional energy spectrum $`E(k)`$

$`C^{uu,bb}(k)={\displaystyle \frac{2(2\pi )^d}{S_d(d1)}}k^{(d1)}E^{u,b}(k)`$ (18)
where $`S_d`$ is the surface area of $`d`$ dimensional spheres. We assume that $`E^u(k)`$ and $`E^b(k)`$ follow

$$E^u(k)=K^u\mathrm{\Pi }^{2/3}k^{5/3};E^b(k)=E^u(k)/r_A$$
(19)
Regarding $`\nu _{(n)}`$ and $`\lambda _{(n)}`$, we attempt the following form of solution

$`(\nu ,\lambda )_{(n)}(k_nk^{})`$ $`=`$ $`(K^u)^{1/2}\mathrm{\Pi }^{1/3}k_n^{4/3}(\nu ,\lambda )_{(n)}^{}(k^{})`$ (20)
with $`k=k_{n+1}k^{}(k^{}<1)`$ with the expectation that $`\nu _{(n)}^{}(k^{})`$ and $`\lambda _{(n)}^{}(k^{})`$ are universal functions for large $`n`$. We numerically solve for $`\nu _{(n)}^{}(k^{})`$ and $`\lambda _{(n)}^{}(k^{})`$. Our calculations reveal that the solutions of $`\nu _{(n)}^{}(k^{})`$ and $`\lambda _{(n)}^{}(k^{})`$ converge for all $`d>d_c2.2`$. From this observation we can conclude that Kolmogorov’s energy spectrum ($`E(k)k^{5/3}`$) is a consistent solution of the RG equations. Meanwhile, Kraichnan’s $`3/2`$ energy spectrum and $`\nu k^2=\lambda k^2kB_0`$, where $`B_0`$ (a constant) is the magnetic field of the large eddies, do not satisfy the renormalization group equations \[Eqs. (15, 17)\]. Hence $`E(k)k^{3/2}`$ is not a consistent solution of RG equations. Our result regarding the nonexistence of stable RG fixed point for $`d=2`$ is consistent with the RG calculation of Liang and Diamond . Refer to Fig. (1) for illustration of $`\nu _{(n)}^{}(k^{})`$ and $`\lambda _{(n)}^{}(k^{})`$ for $`d=3`$ and $`r_A=1`$.

The values of renormalized parameters for $`d=3`$ and various $`r_A`$ are shown in Table I. For large $`r_A`$, the asymptotic $`\nu ^{}`$ is close to the corresponding value for fluid turbulence, but the asymptotic $`\lambda ^{}`$ is also comparable to $`\nu `$. This implies that in fluid dominated regime, there is a significant magnetic energy flux in addition to the usual Kolmogorov’s flux in fluid modes. As $`r_A`$ is decreased, $`\nu ^{}`$ increases but $`\lambda ^{}`$ decreases. This trend is seen till $`r_A0.25`$, where the RG fixed point with nonzero $`\nu ^{}`$ and $`\lambda ^{}`$ becomes unstable, and the trivial RG fixed point with $`\nu ^{}=\lambda ^{}=0`$ becomes stable. This result suggests an absence of turbulence for $`r_A`$ below 0.25. This is consistent with the fact that MHD equations become linear in the $`r_A0`$ (fully magnetic) limit.

We can proceed further and compute various cascade rates and Kolmogorov’s constant for MHD using the renormalized parameters computed above. To compute these quantities we resort to the energy equations, which are

$`\left({\displaystyle \frac{}{t}}+2\nu k^2\right)C^{uu}(𝐤,t)`$ $`=`$ $`{\displaystyle \frac{1}{(d1)(2\pi )^d\delta (𝐤+𝐤^{})}}{\displaystyle _{𝐤^{}+𝐩+𝐪=\mathrm{𝟎}}}{\displaystyle \frac{d𝐩}{(2\pi )^d}}[S^{uu}(𝐤^{}|𝐩|𝐪)+S^{uu}(𝐤^{}|𝐪|𝐩)`$ (22)
$`+S^{ub}(𝐤^{}|𝐩|𝐪)+S^{ub}(𝐤^{}|𝐪|𝐩)]`$
$`\left({\displaystyle \frac{}{t}}+2\lambda k^2\right)C^{bb}(𝐤,t)`$ $`=`$ $`{\displaystyle \frac{1}{(d1)(2\pi )^d\delta (𝐤+𝐤^{})}}{\displaystyle _{𝐤^{}+𝐩+𝐪=\mathrm{𝟎}}}{\displaystyle \frac{d𝐩}{(2\pi )^d}}[S^{bu}(𝐤^{}|𝐩|𝐪)+S^{bu}(𝐤^{}|𝐪|𝐩)`$ (24)
$`+S^{bb}(𝐤^{}|𝐩|𝐪)+S^{bb}(𝐤^{}|𝐪|𝐩)]`$
where

$$S^{uu}(𝐤^{}|𝐩|𝐪)=\mathrm{}\left([𝐤^{}.𝐮(𝐪)][𝐮(𝐤^{}).𝐮(𝐩)]\right),$$
(25)
$$S^{bb}(𝐤^{}|𝐩|𝐪)=\mathrm{}\left([𝐤^{}.𝐮(𝐪)][𝐛(𝐤^{}).𝐛(𝐩)]\right),$$
(26)
$$S^{ub}(𝐤^{}|𝐩|𝐪)=\mathrm{}\left([𝐤^{}.𝐛(𝐪)][𝐮(𝐤^{}).𝐛(𝐩)]\right),$$
(27)
$$S^{bu}(𝐤^{}|𝐩|𝐪)=S^{ub}(𝐩|𝐤^{}|𝐪)$$
(28)
Here $`\mathrm{}`$ stands for the imaginary part of the argument, and the above integrals have constraints that $`𝐤^{}+𝐩+𝐪=\mathrm{𝟎}`$ ($`𝐤=𝐤^{}`$). The energy equations in the above form have been written by Dar et al. , who interpret the terms $`S(𝐤|𝐩|𝐪)`$ as energy transfer rate from mode $`𝐩`$ (second argument of $`S`$) to $`𝐤`$ (first argument of $`S`$) with mode $`𝐪`$ (third argument of $`S`$) acting as a mediator. This interpretation of energy transfer due to Dar et al.  is consistent with the earlier formalism.

We can derive an expression for the energy transfer rate or energy flux from a wavenumber sphere using $`S(𝐤^{}|𝐩|𝐪)`$. The formula for the energy flux from inside of the $`X`$-sphere ($`X<`$) to outside of the $`Y`$-sphere ($`Y>`$) is

$`\mathrm{\Pi }_{Y>}^{X<}(k_0)`$ $`=`$ $`{\displaystyle _{k>k_0}}{\displaystyle \frac{d𝐤}{(2\pi )^d}}{\displaystyle _{p<k_0}}{\displaystyle \frac{d𝐩}{(2\pi )^d}}S^{YX}(𝐤^{}|𝐩|𝐪)`$ (29)
where $`X`$ and $`Y`$ stand for $`u`$ or $`b`$. In our study we assume that the kinetic energy is forced at small wavenumbers, and the turbulence is steady. We calculate the above fluxes analytically to the leading order in perturbation series using the same procedure as that of Leslie . The flux is calculated using Eq. (29) by taking ensemble average of $`S^{YX}`$. The expression for $`S^{bb}(k|p|q)`$ is

$`S^{bb}(k|p|q)`$ $`=`$ $`{\displaystyle _{\mathrm{}}^t}dt^{}[T_4(k,p,q)G^{bb}(k,tt^{})C^{bb}(p,tt^{})C^{uu}(q,tt^{})`$ (32)
$`+T_8(k,p,q)G^{bb}(p,tt^{})C^{bb}(k,tt^{})C^{uu}(q,tt^{})`$
$`+T_{10}(k,p,q)G^{uu}(q,tt^{})C^{bb}(k,tt^{})C^{bb}(p,tt^{})]`$
where $`T_i(k,p,q)`$ are functions of wavevectors $`k,p`$, and $`q`$. The expressions for other transfer rates $`S^{uu}(k|p|q)`$, $`S^{ub}(k|p|q)`$, and $`S^{bu}(k|p|q)`$ look similar. In the above formulas we substitute Kolmogorov’s spectrum \[Eqs. (19)\] for the energy spectrum, and the following expression for the effective viscosity and resistivity

$`(\nu ,\lambda )(k)`$ $`=`$ $`(K^u)^{1/2}\mathrm{\Pi }^{1/3}k^{4/3}(\nu ^{},\lambda ^{})\text{ for }kk_n`$ (33)
Following the same procedure as Leslie , we obtain the following nondimensional form of the equations

$$\frac{\mathrm{\Pi }_{Y>}^{X<}}{\mathrm{\Pi }}=\frac{4S_{d1}}{(d1)^2S_d}(K^u)^{3/2}_0^1𝑑v\mathrm{ln}(1/v)_{1v}^{1+v}𝑑w(vw)^{d2}(\mathrm{sin}\alpha )^{d3}F_{Y>}^{X<}$$
(34)
where the integrals $`F_{Y>}^{X<}`$ are function of $`v`$, $`w`$, $`\nu ^{}`$, and $`\lambda ^{}`$. After a bit of manipulation we can obtain $`\mathrm{\Pi }_{Y>}^{X<}/\mathrm{\Pi }`$ and the constant $`K^u`$. In addition we can also obtain the Kolmogrov’s constant $`K`$ for total energy

$$E(k)=K\mathrm{\Pi }^{2/3}k^{5/3}$$
(35)
using $`K=K^u(1+r_A^1)`$. The values of $`\mathrm{\Pi }_{Y>}^{X<}/\mathrm{\Pi }`$ and $`K`$ for $`d=3`$ and various $`r_A`$ are listed in Table I.

The entries in the Table I show that the cascade rates $`\mathrm{\Pi }_{b>}^{u<}`$, $`\mathrm{\Pi }_{b>}^{b<}`$, $`\mathrm{\Pi }_{b<}^{u<}`$, $`\mathrm{\Pi }_{u>}^{b<}`$are approximately of the same order for $`r_A`$ between 0.5 and 1, but the flux $`\mathrm{\Pi }_{u>}^{u<}`$is rather small. The sign of $`\mathrm{\Pi }_{b>}^{b<}`$is positive, indicating that the ME cascades forward, that is from large length-scales to small length-scales. The magnetic energy thus appearing at small length-scales will be lost due to resistive dissipation, and the large-scale magnetic field is maintained by the $`\mathrm{\Pi }_{b<}^{u<}`$flux. The Kolmogorov’s constant $`K`$ is approximately constant and is close to 1.6, same as that for fluid turbulence ($`r_A=\mathrm{}`$), all $`r_A`$ greater than 0.5.

To summarize, we employed a self-consistent RG scheme to MHD turbulence and found that Kolmogorov’s 5/3 powerlaw is a consistent solution of RG equations for $`dd_c2.2`$. For Kolmogorov’s solution, the renormalized viscosity and resistivity have been calculated, and they are found to be positive. For $`d=3`$, variation of $`\nu ^{}`$ and $`\lambda ^{}`$ with $`r_A`$ shows some interesting features. For large $`r_A`$, $`\nu ^{}`$ is same as that for fluid turbulence, but $`\lambda ^{}`$ is also nonzero, in fact larger than $`\nu ^{}`$. As $`r_A`$ is decreased, $`\nu ^{}`$ increases but $`\lambda ^{}`$ decreases until $`r_A0.25`$ at which value turbulence disappears.

Using the flux equations we have obtained various fluxes and Kolmogorov’s constant $`K`$. For $`d=3`$, $`K`$ does not vary significantly with the variation of $`r_A`$, and it is close to $`K`$ for fluid turbulence. We find that the cascade rate from magnetic-sphere to outside magnetic-sphere ($`\mathrm{\Pi }_{b>}^{b<}`$) is positive, a result consistent with the numerical results of Dar et al. .

In this paper we have restricted ourselves to nonhelical turbulence. Helical MHD turbulence is very important specially in the light of enhancement of magnetic energy (dynamo). However, the physics of helical turbulence is more complex with the appearance of inverse cascade of magnetic helicity etc. The field-theoretic analysis for this case will be taken up later. Recent studies show that the mean magnetic field has a strong effect on energy spectrum, and it induces anisotropy. A full-fledge field theory calculation in the presence of mean magnetic field is also necessary for a clearer picture of MD turbulence

The author thanks J. K. Bhattacharjee for very valuable discussions and ideas. He also thank G. Dar and V. Eswaran for many insights from numerical results.


format: text

text:
# Incoherent Interplane Conductivity of 𝜅-(BEDT-TTF)2Cu[N(CN)2]Br

\[

## Abstract

The interplane optical spectrum of the organic superconductor $`\kappa `$-(BEDT-TTF)2Cu\[N(CN)2\]Br was investigated in the frequency range from 40 to 40,000 cm-1. The optical conductivity was obtained by Kramers-Kronig analysis of the reflectance. The absence of a Drude peak at low frequency is consistent with incoherent conductivity but in apparent contradiction to the metallic temperature dependence of the DC resistivity. We set an upper limit to the interplane transfer integral of $`t_b^2/t_{ac}10^7`$ eV. A model of defect-assisted interplane transport can account for this discrepancy. We also assign the phonon lines in the conductivity to the asymmetric modes of the ET molecule.

\]

Of the $`\kappa `$-type (BEDT-TTF)-based organic superconductors $`\kappa `$-(BEDT-TTF)2Cu\[N(CN)2\]Br has the highest superconducting transition temperature at ambient pressure, $`T_c=12.5`$ K. BEDT-TTF, or bis(ethylenedithio)tetrathiafulvalene and hereafter further abbreviated to ET, is a platelike molecule that, in the $`\kappa `$-type ET-based superconductors, forms conducting layers of orthogonally arranged dimers separated by insulating layers of anions. The resulting quasi-two-dimensional electronic structure makes this class of organic superconductors of great interest due to their potential similarity to the high temperature cuprate superconductors. One important issue is the nature of interlayer transport: is it coherent giving rise to a three dimensional Fermi liquid at sufficiently low temperature, or does it remain incoherent to the lowest temperatures? This issue is well illustrated by two two-dimensional superconductors, Sr2RuO4 which shows coherent interplane transport at low temperature and Bi2Sr2CaCu2O8+δ (Bi-2212) where the interplane transport is incoherent at all temperatures.

The first, Sr2RuO4, is isostructural with the cuprate superconductor La2-xSrxCuO4, and has a superconducting $`T_c`$ of 1.5 K. Although highly anisotropic, it shows coherent interplane transport. At low temperature Sr2RuO4 shows metallic resistivity both in the in-plane and interplane directions, although the resistivity in the interplane direction is two orders of magnitude higher than that of the in-plane direction and has a metallic temperature dependence only below 100 K. The in-plane optical conductivity shows a Drude-like peak, and a weak Drude component can also be seen in the interlayer direction.

The cuprate superconductor Bi-2212 shows an in-plane DC resistivity that, in optimally doped Bi-2212, varies linearly with temperature up to several hundred Kelvin, while in the interplane direction the resistivity is four orders of magnitude higher and increases with decreasing temperature. The real part of the optical conductivity of Bi-2212 shows a Drude-like peak centred at zero frequency in the in-plane direction, but no such peak is seen in the interplane direction. Instead, the interplane conductivity consists entirely of phonon lines. The weak residual interplane transport is due to an incoherent tunneling mechanism.

The transport properties of $`\kappa `$-(ET)2Cu\[N(CN)2\]Br are also highly anisotropic. Within the conducting $`ac`$-planes, the resistivity shows an unusual broad peak near 100 K, but has a clearly metallic behaviour at low temperature with resistivity decreasing with decreasing temperature. The interplane $`b`$-axis resistivity has qualitatively the same metallic temperature dependence as the in-plane resistivity, but is three orders if magnitude higher. Its value of 1 $`\mathrm{\Omega }`$cm at 15 K, which depends somewhat on how quickly the sample is cooled, is six orders of magnitude higher than what is seen in good metals. Thus, in the interplane direction, there is an apparent contradiction between the temperature dependence of the resistivity that suggests coherent transport, and the magnitude of the resistivity that suggests incoherent transport.

One criterion for evaluating the coherence of the interplane conductivity is the Ioffe-Regel-Mott minimum metallic conductivity. Modified for the open Fermi surfaces of these highly anisotropic systems, it gives a lower limit on conductivity which corresponds to a coherence length comparable to the size of a unit cell in the interlayer direction:

$$\sigma _{min}\sqrt{\frac{\sigma _{}}{\sigma _{}}}\frac{e^2}{2\pi ^2\mathrm{}}\frac{l_{}}{l_1l_2}$$
(1)
where $`\sigma _{}`$ and $`l_{}`$ are the conductivity and lattice constant in the interplane direction and $`\sigma _{}`$, $`l_1`$ and $`l_2`$ are the conductivity and lattice constants in the planes. Table I lists these parameters for Sr2RuO4, Bi-2212 and $`\kappa `$-(ET)2Cu\[N(CN)2\]Br. The conductivity anisotropy ratio should be evaluated at $`\sigma _{}=\sigma _{min}`$. For Sr2RuO4 this gives coherent conductivity below 52 K which is indeed the temperature where the resistivity starts to deviate from the low-temperature $`T^2`$ dependence. For Bi-2212 and $`\kappa `$-(ET)2Cu\[N(CN)2\]Br, however, $`\sigma _{min}>\sigma _{}`$ even near $`T_c`$ where the anisotropy is largest. This suggests that the low temperature interplane conductivity of Sr2RuO4 is coherent but becomes incoherent above 52 K due to thermal fluctuations, while both Bi-2212 and $`\kappa `$-(ET)2Cu\[N(CN)2\]Br remain incoherent at all temperatures.

Optical conductivity provides another method of investigating the coherence of interplane transport where in a coherent system the optical conductivity $`\sigma _1`$ shows a Drude peak centred at zero frequency with a maximum corresponding to the DC conductivity $`\sigma _{DC}`$ and a width equal to the scattering rate $`\mathrm{\Gamma }`$ of the free carriers:

$$\sigma _1(\omega )=\frac{\sigma _{DC}\mathrm{\Gamma }^2}{\omega ^2+\mathrm{\Gamma }^2}$$
(2)
The question of the coherence of the interplane transport in the ET-based organic superconductors can be addressed with measurements of low temperature far-infrared optical conductivity in the interplane direction. To date however, infrared studies have focused on the in-plane properties,  and the few interplane measurements which have been made  were at room temperature and above 600 cm-1. This is due to the difficulty of growing crystals with large faces perpendicular to the conducting planes. Recently, however, high quality crystals of sufficient size for far-infrared interplane measurements have become available. The interplane reflectance measurements of $`\kappa `$-(ET)2Cu\[N(CN)2\]Br presented in this paper are the first such measurements on any ET based superconductor.

The single crystals of $`\kappa `$-(ET)2Cu\[N(CN)2\]Br were synthesized by the electrocrystallization technique described elsewhere. Typical crystal sizes were $`1.5\times 1.5\times 1.5`$ mm with faces as large as 1 mm2 parallel to the interlayer $`b`$-axis. Polarized reflectance measurements between 40 and 8000 cm-1 were performed on these as-grown faces with a Michelson interferometer using three different detectors. A grating spectrometer with three additional detectors was used to make measurements at 300 K for the rest of the range up to 40,000 cm-1 (5 eV).

The reflectance of $`\kappa `$-(ET)2Cu\[N(CN)2\]Br with the light polarized in the interplane direction is shown in Fig. 1 for four temperatures above the superconducting transition temperature. The reflectance is approximately 0.15 over the entire range with several sharp phonon peaks at low frequencies and some broader interband-like features at higher frequencies.

It should also be pointed out that there was some sample-to-sample variation of the interplane reflectance, particularly the temperature dependence of the background reflectance above 200 cm-1. Fig. 2 is a comparison of the sample shown in Fig. 1 (upper panel) with another sample showing a much stronger temperature dependence (lower panel). The phonon lines below 200 cm-1 also appear to be stronger in this second sample although the optical conductivity should be calculated for a true comparison. Unfortunately we were unable to measure reflectance above 800 cm-1 for this second sample, and we have no explanation for the variation. Complete data sets need to be collected on more samples to properly investigate this phenomenon. It is interesting to note that the feature has the appearance of a gap, and evidence for a pseudogap in $`\kappa `$-(ET)2Cu\[N(CN)2\]Br at $`T^{}=50`$ K in ESR and 13C NMR measurements has been reported. This would give $`2\mathrm{\Delta }/k_BT^{}=5.7`$ which is not far from the value 4.3 reported for the high temperature cuprate superconductors.

The calculation of the optical conductivity requires extrapolation of the reflectance to all frequencies for the Kramers-Kronig analysis. The spectra were extrapolated to high frequencies using power law extrapolations: $`\omega ^1`$ above 5 eV and $`\omega ^4`$ above 62 eV. Below our lowest measured point at 25 cm-1 we used a Drude extrapolation, although an insulating extrapolation gives very similar results. We estimate our experimental uncertainty of the reflectance to be $`\pm 0.005`$. Combined with uncertainties due to the extrapolations, this gives an uncertainty in our optical conductivity of $`\pm 8\%`$ between 200 and 2000 cm-1. Outside this range the uncertainty rises, reaching $`\pm 40\%`$ at 25 and 5000 cm-1. The resolution of the spectra is 2 cm-1 up to 200 cm-1, 4 cm-1 up to 680 cm-1, and 15 cm-1 up to 8000 cm-1.

The real part of the interplane optical conductivity is shown in Fig. 3. The value of the DC conductivity from the work of Su et al. at 15 K is shown at 1 $`(\mathrm{\Omega }cm)^1`$ on the vertical axis. Clearly there is no sign of the usual Drude peak that accompanies coherent conductivity. Instead the conductivity is dominated by sharp phonon lines on a background due to interband-like features in the mid-infrared. This confirms that the interlayer conductivity in $`\kappa `$-(ET)2Cu\[N(CN)2\]Br is incoherent.

A fit to Lorentz oscillators which will be discussed later shows no sign of any free carrier component to the optical conductivity above 30 cm-1. This is similar to what has been reported for Bi-2212. It is difficult to estimate an upper limit for the plasma frequency of a hypothetical Drude peak, but assuming the width is equal to the in-plane width reported by Eldridge and Kornelsen of $`\mathrm{\Gamma }=20`$ cm-1 we estimate $`\omega _p=50`$ cm-1 for the plasma frequency of the Drude component. Combining the in-plane scattering rate $`\mathrm{\Gamma }=20`$ cm-1 with the DC resistivity gives a somewhat lower value of $`\omega _p=25`$ cm-1. Using the 50 cm-1 value and a simplified tight-binding model we can estimate the interplane transfer integral using

$$\omega _p^2=\frac{e^2}{ϵ_0\mathrm{}^2}\underset{BZ}{}f(E_𝐤)\frac{^2E}{k_\mu ^2}$$
(3)
where $`f(E_𝐤)`$ is the Fermi-Dirac occupation number and the derivative is to be taken in the direction of the field. We use the simplified tight-binding band

$$E_𝐤=2t_a\mathrm{cos}(k_ad_a)2t_b\mathrm{cos}(k_bd_b)2t_c\mathrm{cos}(k_cd_c)$$
(4)
where $`d`$ are the ET molecular repeat distances and $`t`$ are the average transfer integrals along the various directions. We assume an open Fermi surface in the interplane direction with $`t_bt_{ac}`$ and ignore the in-plane anisotropy to get

$$\omega _{pb}\frac{e^2}{ϵ_0\mathrm{}^2}\frac{b^2/V_m}{\sqrt{\pi }\mathrm{sin}(\sqrt{\pi })}\frac{t_b^2}{t_{ac}}$$
(5)
where $`b`$ is the interplane lattice constant and $`V_m`$ the volume per ET molecule. This gives $`t_b^2/t_{ac}10^7`$ eV which gives $`t_b10^4`$ eV if $`t_{ac}10^1`$ eV as has been estimated for (ET)2I3. Of course, this analysis assumes a coherent component to the conductivity which seems unlikely given the earlier discussion of minimum metallic conductivity.

We now return to the problem of the temperature dependence of the interplane DC conductivity which shows a clear metallic character and in fact follows rather accurately the in-plane conductivity although it is a factor of 1000 smaller. Since the interplane conductivity is below the minimum metallic limit it is highly likely that its metallic temperature dependence is due to some special process that makes it mirror the in-plane conductivity. There are several models that can do this and we focus here on a model originally proposed for Bi-2212 by Martin et al., where the conducting planes are connected by a random network of shorts approximated as a regular array of links of resistance $`R_{lb}`$ distance $`\xi `$ apart along the conducting plane. Given the resistance $`R_{lac}`$ along the plane between the shorts and the interplane distance $`d`$, the apparent interplane resistivity is given by

$$\rho _b=(R_{lb}+R_{lac})\times \frac{(\xi /2)^2}{d}$$
(6)
and assuming $`R_{lac}=\rho _{ac}/d`$ we get for the anisotropy

$$\frac{\rho _b}{\rho _{ac}}=\left(1+\frac{dR_{lb}}{\rho _{ac}}\right)\times \left(\frac{\xi }{2d}\right)^2$$
(7)
Buravov’s data above 50 K are nearly temperature-independent with $`\rho _b/\rho _{ac}1000`$, while below 50 K the anisotropy rises dramatically. Fig. 4 compares the interplane resistivity and anisotropy of Buravov et al. with a least squares fit using Eq. 7 assuming a temperature-independent $`R_{lb}`$. The fit, which reproduces both the temperature independence of the anisotropy at high temperatures and the dramatic increase in anisotropy at low temperatures, gives $`R_{lb}=138\pm 9`$ k$`\mathrm{\Omega }`$ and $`\xi =802\pm 7`$ Å if we take $`d=15`$ Å. The model explains the rise of the anisotropy below 50 K as due to the rapid drop of $`R_{lac}`$ relative to the temperature-independent $`R_{lb}`$, whereas models where the proportionality between in-plane and interplane resistivity is built in predict a constant anisotropy, independent of temperature. The observed sample-to-sample variation shown in Fig. 2 is also consistent with a process controlled by defects.

The model of a meandering current will fail at high frequency. To estimate the characteristic maximum frequency we model the current path as a transmission line of series resistors $`R_{lac}`$ shunted to ground by capacitors $`C_s`$ where we set $`C_s2\pi \xi `$. The impedance of such a transmission line is

$$Z=\frac{R_{lac}}{2}\left(1+\sqrt{1+\frac{4i}{\omega R_{lac}C_s}}\right)$$
(8)
Thus above a frequency given by $`\omega _0=(R_{lac}C_s)^1`$ the capacitors short the AC current and the in-plane fields do not have time to build up. With the parameters determined above we find that $`\omega _010`$ cm-1. A test of the model would be a reduction of the interplane conductivity from its DC value to a much lower value at this frequency. Dressel et al.  find in the millimeter wave range (1-3 cm-1) conductivities that agree with DC values consistent with our picture. Unfortunately the strong phonon background discussed in the next paragraphs prevents us from giving an accurate value of the interplane electronic conductivity in the far infrared range other than an upper limit of about 1 ($`\mathrm{\Omega }`$cm)-1.

It was possible to fit the conductivity with a series of Lorentzian oscillators according to

$$\sigma _1(\omega )=\frac{1}{4\pi }\underset{i}{}\frac{\omega _{pi}^2\omega ^2\mathrm{\Gamma }_i}{(\omega _{0i}^2\omega ^2)^2+\omega ^2\mathrm{\Gamma }_i^2}$$
(9)
51 oscillators were used in the fit along with a weak Drude term to account for phonons below the measurement range and one strong oscillator at high frequency to provide the background tail of the mid-infrared features. It was necessary to add some asymmetry to the strongest oscillators at 492 and 1450 cm -1, but all other oscillators were symmetric. Table II lists the frequency $`\omega _0`$, plasma frequency $`\omega _p`$ and width $`\mathrm{\Gamma }`$ of the oscillators at all four temperatures. Most of the low frequency lines seem to be lattice modes as their frequencies increase with decreasing temperature. Some of these lines were previously reported in a powder absorption experiment and are marked with an asterisk. The rest of the lines are related to internal vibrations of the ET molecule, and these mode assignments were made by comparison with Eldridge et al. who assign the normal modes of the ET molecule and relate these ET modes to observed lines in in-plane infrared and Raman spectra of $`\kappa `$-(ET)2Cu\[N(CN)2\]Br.

In general, infrared spectra are sensitive to asymmetric (ungerade) modes while Raman spectra are sensitive to symmetric (gerade) modes. Since the ET molecule consists of two mirrored halves joined by a single C=C bond, the infrared and Raman spectra of ET contain similar sets of lines.  Each vibration of atoms in one half of the molecule can be in phase or out of phase with an identical vibration in the other half producing symmetric/asymmetric pairs of modes. Since the two halves are nearly independent, the members of each mode pair have nearly the same energy. This argument does not apply to modes involving the central C=C bond.

Like the ET spectra, the in-plane infrared and Raman spectra of $`\kappa `$-(ET)2Cu\[N(CN)2\]Br also have similar sets of lines, however in this case both sets of lines are ascribed to the symmetric modes of ET. The reason the symmetric ET modes are infrared-active in $`\kappa `$-(ET)2Cu\[N(CN)2\]Br is the arrangement of ET molecules into dimers so that the members of a dimer can vibrate out of phase transferring charge back and forth and producing the dipole moment required for infrared activity. The Raman line positions in $`\kappa `$-(ET)2Cu\[N(CN)2\]Br are shifted from the ET line positions due to the charge on the ET molecules, and the infrared line positions are further shifted by their coupling to the charge transfer between members of a dimer. To assist them in their mode assignments, Eldridge et al. also cause further shifts using several isotopic substitutions.

We have assigned the lines in our interplane infrared spectra to the asymmetric modes of the ET molecule. The mechanism that makes the symmetric modes infrared-active in the in-plane measurements does not apply to our interplane measurements since the electric vector of the light is perpendicular to the direction of charge transfer between members of a dimer. Fig. 5 shows a comparison of our $`\sigma _1`$ to the ET powder absorption spectrum of Eldridge et al. For the most part, our mode assignments are based on this comparison since we do not have many measurements on isotopic analogs to assist us. We do however have a partial reflectance spectrum of a $`\kappa `$-(13C(2)-ET)2Cu\[N(CN)2\]Br crystal in which the two central carbon atoms have been substituted with 13C. Fig. 6 shows the 23 cm-1 shift of a line at 791.2 cm-1 confirming its identity as the $`\nu _{32}`$(B1u) mode. Since this mode involves the central C=C bond, its symmetric counterpart is at a very different frequency and is seen in the in-plane infrared and Raman spectra near 450 cm-1 while no line is seen near 791 cm-1.

In general, all of the lines become narrower as temperature decreases as expected, however the plasma frequencies, or line strengths, also increase significantly in some cases. In particular, the line at 81 cm-1, which is quite strong at 15 K, has nearly disappeared above 200 K. An increase in line strength with decreasing temperature in organic conductors has generally been interpreted as a signature of charge-density-wave fluctuations.  However, temperature-dependent phonon intensities have also been associated with spin-density-wave transitions.

In summary we have measured the interplane reflectance of the quasi-two-dimensional organic superconductor $`\kappa `$-(ET)2Cu\[N(CN)2\]Br and calculated the optical conductivity using Kramers-Kronig relations. We find strong evidence of incoherent transport as has been reported for the high temperature superconductors, and estimate an upper limit for the free carrier plasma frequency of 50 cm-1 from which we derive an upper limit for the interplane transfer integral of $`10^4`$ eV. We have used a defect model to explain the crossover from a temperature-independent resistivity anisotropy at high temperatures to a rapidly increasing anisotropy at low temperatures. We have also fit the phonon lines in the conductivity to a series of Lorentzian oscillators and assigned these to the asymmetric modes of the ET molecule.

The work at McMaster University was supported by the Natural Sciences and Engineering Research Council of Canada. Work at Argonne National Laboratory is sponsored by the U.S. Department of Energy, Office of Basic Energy Sciences, Division of Materials Sciences, under Contract W-31-109-ENG-38.


text:
# Phase transition dynamics in the hot Abelian Higgs model

## I Introduction

Whilst there are many useful techniques for studying the equilibrium properties of finite-temperature field theories, understanding the non-equilibrium dynamics is a much harder task. Nevertheless, it would be essential for many fields of physics, for instance cosmology, heavy ion physics and condensed matter physics. In all these fields, new empirical data will be available in near future, which would allow the theories to be tested, but the complexity and the non-equilibrium nature of the phenomena make it difficult to derive theoretical predictions that could be compared with the data.

One fairly generic consequence of phase transitions is formation of topological defects . If the phase transition is associated with a spontaneous breakdown of a global symmetry, this process is well understood. The correlation length of the order parameter cannot keep up with its equilibrium value, which diverges at the transition point. The direction of the symmetry breaking must therefore be uncorrelated at long distances, and at places where these correlated domains meet, topological defects are formed. This is called the Kibble-Zurek (KZ) mechanism (see e.g.  for a review).

If the symmetry that gets broken is a local gauge invariance, the above argument cannot be used directly, because the direction of the order parameter is not a gauge invariant quantity. We studied this recently in the context of the Abelian Higgs model , and pointed out that the thermal fluctuations of the magnetic field lead to another mechanism that forms topological defects. The argument was based on fairly generic assumptions, but leads to some concrete predictions that were confirmed in numerical simulations.

The aim of this paper is to study in more detail the dynamics of the Abelian Higgs model during the phase transition from the Coulomb phase to the Higgs phase. In particular, we concentrate on those degrees of freedom that are relevant for defect formation. This allows us to test the scenario of Ref.  on a more quantitative level.

The theory considered in Ref.  was classical, and although the same arguments apply to the quantum theory as well, the details of the dynamics are different. The full quantum field theory cannot be simulated in practice, but one can argue that the dynamics of the relevant long-wavelength degrees of freedom are classical . By integrating out the short-wavelength fluctuations perturbatively, one obtains a classical effective theory with non-local interactions , to which we refer as the hard-thermal-loop (HTL) improved theory. In order to understand how the quantum effects change the dynamics, we simulate this HTL improved theory using the method developed in Ref. .

The structure of the paper is the following: in Sect. II we present both the classical and HTL improved Abelian Higgs models. In Sect. III we discuss defect formation in the model, comparing the mechanism presented in  with the Kibble-Zurek scenario. In Sects. IV and V we describe our numerical simulations and present the results. Conclusions are given in Sect. VI and technical details of the HTL improved equations of motion and the lattice formulation in the two Appendices.

## II Abelian Higgs model

The Abelian Higgs model is defined by the Lagrangian

$$=\frac{1}{4}F_{\mu \nu }F^{\mu \nu }+|D_\mu \varphi |^2m^2|\varphi |^2\lambda |\varphi |^4,$$
(1)
where $`D_\mu =_\mu +ieA_\mu `$ and $`F_{\mu \nu }=_\mu A_\nu _\nu A_\mu `$.

A particularly interesting feature of this theory is the existence of Nielsen-Olesen vortex solutions . These string-like topological defects are characterized by a zero of the Higgs field at the centre of the vortex around which the Higgs phase angle has a non-zero winding number

$$n_C=_C𝑑\stackrel{}{r}\gamma (\stackrel{}{r})0.$$
(2)
Here $`C`$ is a closed path around the vortex and $`\gamma `$ is the Higgs phase angle, i.e., $`\varphi =|\varphi |\mathrm{exp}(i\gamma )`$.

At finite temperature, perturbation theory is plagued by infrared divergences , which can be partly cured by a resummation of the perturbative expansion, but even the resummed expansion breaks down near the transition. Static equilibrium quantities, such as the phase diagram of the theory, can still be calculated reliably with non-perturbative Monte Carlo simulations.

The model has a phase transition between the high-temperature Coulomb phase and the low-temperature Higgs phase at $`T^2=T_c^212(m^2)/(3e^2+4\lambda )`$. In the perturbative regime ($`\lambda e^2`$), the transition is of first order , and at larger $`\lambda `$ it becomes continuous . There are no local order parameters, but a number of non-local ones: the photon mass and the vortex tension are non-zero in the broken phase and vanish at the transition . Therefore the transition is not a smooth crossover like the electroweak phase transition .

Monte Carlo simulations cannot be used for real-time quantities in the quantum theory, because the necessary path integral is not Euclidean but consists of a complicated path in complex time . However, we can utilize the fact that modes with different momenta behave in very different ways . The soft, long-wavelength modes ($`kT`$) have large occupation numbers, and they can be approximated very well by a classical theory. This makes numerical simulations feasible, because the time-evolution of a classical field theory can be found simply by solving the equations of motion numerically.

### A Classical theory at finite temperature

Classical field theory at finite temperature is ultraviolet divergent, and thus the results depend on the lattice spacing $`\delta x`$. Divergences like these are generic to all low-energy effective theories, and are exactly cancelled by corrections the high-momentum modes induce to the effective Lagrangian. If one is only interested in static equilibrium quantities, these corrections can be calculated in the limit of high temperature and small lattice spacing $`\delta x`$ ,

$$m_T^2=m^2+(3e^2+4\lambda )\left(\frac{T^2}{12}\frac{3.176T}{4\pi \delta x}\right).$$
(3)
The approach we use in Sect. IV is to take this correction into account and solve the classical equations of motion

$`_\mu F^{\mu \nu }`$ $`=`$ $`2e\mathrm{Im}\varphi ^{}D^\nu \varphi ,`$ (4)
$`D_\mu D^\mu \varphi `$ $`=`$ $`m_T^2\varphi 2\lambda (\varphi ^{}\varphi )\varphi .`$ (5)
However, as was pointed out in Ref. , these equations do not reproduce the real-time dynamics of the quantum theory correctly. Thus, the results are not reliable on a quantitative level, but they can still give a reasonably good qualitative picture of the dynamics, and provide a non-trivial test for the scenario presented in Section III.

### B HTL improved theory

If the couplings are small, the system is close to thermal equilibrium and $`T\delta x1`$, it is possible to construct a classical theory which approximates the dynamics of the original quantum theory to leading-order accuracy in the coupling constants .

Near the phase transition, $`m^2e^2T^2`$, and we can use the high temperature approximation $`Tm`$ in our loop integrals, provided $`e`$ is small. We calculate the one-loop corrections from the hard modes to the self-energies of $`\varphi `$ and $`A_i`$, and resum them into the effective Lagrangian

$`_{\text{HTL}}`$ $`=`$ $`{\displaystyle \frac{1}{4}}F_{\mu \nu }F^{\mu \nu }{\displaystyle \frac{1}{4}}m_D^2{\displaystyle \frac{d\mathrm{\Omega }}{4\pi }F^{\mu \alpha }\frac{v_\alpha v^\beta }{(v)^2}F_{\mu \beta }}`$ (7)
$`+|D_\mu \varphi |^2m_T^2|\varphi |^2\lambda |\varphi |^4,`$
where $`m_T^2`$ is given by Eq. (3), and the integration is taken over the unit sphere of velocities $`v=(1,\stackrel{}{v})`$, $`\stackrel{}{v}^2=1`$. The Debye mass has the value $`m_D^2=\frac{1}{3}e^2T^2+\delta m_D^2`$, where $`\delta m_D^2`$ is a counterterm that cancels the UV divergences and is discussed in more detail in Appendix B 2.

All the degrees of freedom remaining in this effective theory are classical, and therefore it can be treated as a classical theory. The time evolution of the fields will then be determined by the equations of motion \[cf. Eq. (5)\]:

$`_\mu F^{\mu \nu }`$ $`=`$ $`m_D^2{\displaystyle \frac{d\mathrm{\Omega }}{4\pi }\frac{v^\nu v^i}{v}E^i}2e\mathrm{Im}\varphi ^{}D^\nu \varphi ,`$ (8)
$`D_\mu D^\mu \varphi `$ $`=`$ $`m_T^2\varphi 2\lambda (\varphi ^{}\varphi )\varphi .`$ (9)
As such, the equations of motion (9) are not well suited for our purposes. Firstly, it is not obvious how to find the corresponding Hamiltonian, which is necessary for preparing the initial configurations. And secondly, the equations of motion are non-local both in space and time.

In Ref. , a convenient local formulation was presented along the lines of Ref. . The latter work introduces a new local field $`W(t,\stackrel{}{x},\stackrel{}{v})`$, representing the departure from the equilibrium distribution function for hard particles with velocity $`\stackrel{}{v}`$. Since the hard particles move at the speed of light, $`\stackrel{}{v}`$ consists of two free coordinates, making $`W`$ a six-dimensional field. In our formulation, part of the velocity dependence decouples from the dynamics of the soft modes, and we can describe exactly the same dynamics with two five-dimensional fields, $`\stackrel{}{f}(t,\stackrel{}{x},z)`$ and $`\theta (t,\stackrel{}{x},z)`$, where $`z[0,1]`$ is the cosine of the angle between the hard mode velocity and the gradient of the distribution function. The details of this formulation are given in Appendix A.

In order to simulate the model numerically, we define the theory on a periodic spatial lattice, and the $`z`$ dependence of $`\stackrel{}{f}`$ and $`\theta `$ is discretized by expressing them as sums over a finite number of Legendre modes, whose order ranges from $`0`$ to $`N_{\mathrm{max}}`$ . The resulting equations of motion are detailed in Appendix A. We merely note here that, as shown in Ref. , for a given value of $`N_{\mathrm{max}}`$ and a given momentum $`k`$, the approximation breaks down at times

$$t>t_c(k)=4N_{\mathrm{max}}/k.$$
(10)
Therefore one can strictly speaking only measure correlators up to $`\mathrm{\Delta }tN_{\mathrm{max}}\delta x`$, but since the modes that are in equilibrium will remain in equilibrium even beyond that, one can essentially trust the results until $`\mathrm{\Delta }tt_c(k)`$, where $`k`$ is the relevant momentum scale.

In Sect. V, we study the dynamics of the model using the HTL improved equations of motion, and describe how these corrections change the results from the classical case.

## III Defect formation

In cosmology, temperature of the universe decreases as a result of its expansion, and this leads to phase transitions. If we write the equations of motion in conformal coordinates, the effect of the expansion can be absorbed completely into a varying zero temperature Higgs mass term, $`m^2=m^2(t)`$.

In order to see this, we perform the conformal rescaling $`dx^\mu adx^\mu `$, $`v^\mu av^\mu `$, $`A_\mu a^1A_\mu `$ and $`\varphi a^1\varphi `$, where $`a=a(t)`$ is the scale factor of the Universe. If we can neglect the expansion rate in comparison to the microscopic timescales of the theory, which are given by $`m_T^1`$ and $`m_D^1`$, the effect of this rescaling on the action is to replace the mass terms $`m_T^2`$ and $`m_D^2`$ by $`(m_Ta)^2`$ and $`(m_Da)^2`$ respectively. The hard modes are assumed to be ultra-relativistic, so they will stay close to a thermal distribution with temperature $`T=\overline{T}a^1`$. Hence $`(m_Da)^2`$ stays constant, as do the rescaled thermal corrections to the Higgs mass, leaving the only time dependence in the parameter $`m^2(t)=m^2a^2(t)`$. In the following, we will assume that the phase transition is triggered in this way, by a mass term decreasing below a critical value, but we believe that the qualitative features would not be very different if some of the other parameters were changing at the same time.

When the system enters the Higgs phase, Nielsen-Olesen vortices are formed. In the limit $`e0`$, this can be understood in terms of the Kibble mechanism . When $`m^2`$ approaches its critical value, the Higgs correlation length $`\xi `$ grows, and if the system could remain in equilibrium it would eventually diverge at the transition point. However, $`\xi `$ cannot grow arbitrarily fast, because at the very least it is constrained by the finite speed of light. Therefore it remains finite if the transition takes place in finite time. At the transition point, the system consists of correlated domains of size $`\widehat{\xi }`$ determined by the maximum correlation length reached. In each of these domains, the phase angle of the Higgs field is chosen independently of all others, and this gives rise to frustrations, vortices, where these domains meet. Up to a numerical factor, the number of vortices piercing a unit area is

$$n=N/A\widehat{\xi }^2.$$
(11)
In practice, the maximal rate of change of the correlation length may be well below the speed of light, and this argument can indeed be made more precise by considering the dynamics of the system in more detail .

In Ref. , we argued that this picture is not complete if $`e>0`$. After all, the phase angle of the Higgs field is not gauge invariant, and therefore arguments based on it cannot apply. Nevertheless, if the amplitude of the magnetic field is small, we can fix the gauge in which $`A_i0`$, and then we can use the above picture in this gauge. In these cases the above Kibble-Zurek scenario should work. However, if the initial state is at a non-zero temperature, magnetic field is never exactly zero, and it must be taken into account.

In the symmetric phase, the thermodynamics of the gauge field is described by ordinary electrodynamics. The energy of any magnetic field configuration is approximately given by

$$H_{\mathrm{EM}}[\stackrel{}{B}(\stackrel{}{x})]=\frac{1}{2}d^3x\stackrel{}{B}(\stackrel{}{x})^2,$$
(12)
and the probability with which the thermal fluctuations can generate that configuration is proportional to $`\mathrm{exp}(H_{EM}/T)`$. Let us consider a circular loop $`C`$ that bounds a surface $`S`$. The magnetic flux through the surface is

$$\mathrm{\Phi }_S=𝑑\stackrel{}{S}\stackrel{}{B}.$$
(13)
Although this is zero on the average, it is non-zero in almost all configurations, and therefore the typical value, given by $`\sqrt{\mathrm{\Phi }_S^2}`$ is non-zero. We can estimate its value by calculating the energy of the field configuration that minimizes the energy for a given value of $`\mathrm{\Phi }_S`$. This minimal configuration is simply a magnetic dipole, and its energy is

$$E_{\mathrm{min}}(\mathrm{\Phi }_S)\mathrm{\Phi }_S^2/R,$$
(14)
where $`R`$ is the radius of the loop $`C`$. Thermal fluctuations can create this configuration if $`T>E_{\mathrm{min}}(\mathrm{\Phi }_S)`$, and solving this for $`\mathrm{\Phi }_S`$ shows that the typical flux through the loop is

$$\mathrm{\Phi }_S\sqrt{TR}.$$
(15)
When the system enters the Higgs phase, magnetic flux is confined into flux tubes, which costs energy. Therefore the dynamics tries to decrease the magnetic flux, but it can only do so at short distances, because the range of the interaction between the flux tubes decreases rapidly.

We can be more specific in Fourier space. We can write the two-point correlator of the magnetic flux density $`\stackrel{}{B}`$ as

$$B_i(\stackrel{}{k})B_j(\stackrel{}{k}^{})=\left(\delta _{ij}\frac{k_ik_j}{k^2}\right)(2\pi )^3\delta (\stackrel{}{k}+\stackrel{}{k}^{})G(k).$$
(16)
In the symmetric phase, different Fourier modes $`\stackrel{}{B}(\stackrel{}{k})`$ behave as independent oscillators in thermal bath, and thus each of them has the same amplitude $`G_0(k)=T`$. When the system enters the broken phase, magnetic field becomes massive, and the equilibrium distribution changes into

$$G(k)=T\frac{k^2}{k^2+m_\gamma ^2}.$$
(17)
In order to stay in equilibrium, the amplitude of the long-wavelength modes must drop rapidly, but the time scale $`\tau (k)`$ of the dynamics of these modes is very slow. It depends on the individual system, but in general $`lim_{k0}\tau (k)=\mathrm{}.`$ Thus the modes with $`k`$ less than some critical value $`\widehat{k}`$ cannot remain in equilibrium. If we know $`\tau (k)`$, we can calculate $`\widehat{k}`$ from the condition

$$\left|\frac{d\mathrm{ln}G(\widehat{k})}{dt}\right|\frac{1}{\tau (\widehat{k})}.$$
(18)
The consequence of the above process is that there will be long-wavelength magnetic fields present even in the Higgs phase. In particular, at distances less than $`\widehat{\xi }=2\pi /\widehat{k}`$ it looks like there was a uniform external magnetic field. We can estimate that its amplitude $`B_{\mathrm{avg}}`$ is

$$B_{\mathrm{avg}}^2_0^{\widehat{k}}\frac{d^3k}{(2\pi )^3}B_i(\stackrel{}{k})_0^{\widehat{k}}\frac{d^3k^{}}{(2\pi )^3}B_i(\stackrel{}{k}^{})_0^{\widehat{k}}\frac{d^3k}{(2\pi )^3}G_0(k)T\widehat{k}^3.$$
(19)
This magnetic flux must be confined into vortices, and because each vortex carries one flux quantum $`\mathrm{\Phi }_0=2\pi /e`$, the number density of vortices per unit area is

$$n\frac{e}{2\pi }B_{\mathrm{avg}}\frac{e}{2\pi }T^{1/2}\widehat{k}^{3/2}.$$
(20)
Even if we do not know $`\tau (k)`$, we can still make some concrete predictions based on this scenario. For instance, the spatial distribution of vortices turns out to be completely different than in the KZ scenario. At distances shorter than $`\widehat{\xi }=2\pi /\widehat{k}`$, the magnetic field points in the same direction, which means that the vortices tend to be aligned, while in the KZ scenario they prefer to be anti-aligned. This prediction was confirmed in the simulations in Ref. .

## IV Classical simulations

### A Equilibrium

In order to test the scenario presented in Sect. III, we carried out a number of numerical simulations. Let us first discuss the simulations of the classical model (5).

Because we are interested in transitions that start from close to thermal equilibrium, we will first have to thermalize the system. This means preparing an ensemble of field configurations with the canonical equilibrium distribution $`\mathrm{exp}(\beta H)`$, where $`\beta =1/T`$ and the Hamiltonian $`H`$ is

$$H=d^3x\left[\frac{1}{2}\stackrel{}{E}^2+\frac{1}{2}(\stackrel{}{}\times \stackrel{}{A})^2+\pi ^{}\pi +(D_i\varphi )^{}(D_i\varphi )+m_T^2\varphi ^{}\varphi +\lambda (\varphi ^{}\varphi )^2\right].$$
(21)
Here $`\pi =_0\varphi `$ is the canonical momentum of $`\varphi `$, and in the temporal gauge ($`A_0=0`$), the electric field is simply $`\stackrel{}{E}=_0\stackrel{}{A}`$. In addition, the fields must satisfy the Gauss law as an extra constraint

$$\stackrel{}{}\stackrel{}{E}=2e\mathrm{Im}\varphi ^{}\pi .$$
(22)
Because of the constraint (22), a straightforward Metropolis algorithm would not work very well. Instead, we used a hybrid algorithm, in which we thermalized the component of $`\pi `$ orthogonal to Eq. (22) with a heat bath algorithm, and performed a number of Metropolis updates to the gauge field $`\stackrel{}{A}`$. Because $`\stackrel{}{A}`$ does not appear in Eq. (22), it leaves the constraint unchanged. Then we evolved the system with the equations of motion

$`_0\stackrel{}{A}`$ $`=`$ $`\stackrel{}{E},`$ (23)
$`_0\varphi `$ $`=`$ $`\pi ,`$ (24)
$`_0\stackrel{}{E}`$ $`=`$ $`\stackrel{}{}\times \stackrel{}{}\times \stackrel{}{A}+2e\mathrm{Im}\varphi ^{}D_i\varphi ,`$ (25)
$`_0\pi `$ $`=`$ $`D_iD_i\varphi m_T^2\varphi 2\lambda (\varphi ^{}\varphi )\varphi ,`$ (26)
which, again, leave the Gauss law unchanged. We repeated this procedure a number of times so that the system thermalized.

In our simulations, we used the couplings $`e=0.3`$ and $`\lambda =0.18`$. The lattice size was $`120\times 120\times 20`$ (the reason for choosing one short dimension is discussed in Sec. IV B), the lattice spacing was $`\delta x=6T^1`$ and the time step was $`\delta t=0.3T^1`$. The details of the lattice implementation are given in Appendix B 1.

Since $`\lambda >e^2`$, the phase transition is continuous. In order to determine the location of the transition point, and to test the accuracy of the tree-level result (17), we measured the correlator (16) at various values of $`m^2`$, starting deep in broken phase. For thermalization to each value of $`m^2`$, we used 24 hybrid Monte Carlo cycles each consisting of 400 Metropolis sweeps and time evolution for $`\mathrm{\Delta }t=600T^1`$. We carried out the measurement in nine independently thermalized configurations, measuring the average correlator in an unperturbed run of length $`\mathrm{\Delta }t=12000T^1`$ (and $`\mathrm{\Delta }t=72000T^1`$ for $`m^2=0.083T^2`$). The results are shown in Fig. LABEL:fig:spatcorra. The solid lines show that the agreement with Eq. (17) is excellent, and we find that the phase transition takes place at $`m^20.083T^2`$.

In the symmetric phase, the tree-level result (17) corresponds to a constant $`G(k)=T`$, but the data measured at $`m^2=0.083T^2`$ (see Fig. LABEL:fig:spatcorrb), clearly turns down at small momenta. It is customary to parameterize the corrections to the tree-level result by introducing the static photon self energy $`\mathrm{\Pi }_T`$, defined by

$$G(k)=T\frac{k^2}{k^2+\mathrm{\Pi }_T(k)}.$$
(27)
Our measurements seem to contradict the results of Kraemmer et al. and Blaizot et al. , who showed that after a resummation, the lowest-order term of $`\mathrm{\Pi }_T`$ is proportional to $`k^2`$. Such a quadratic term would only change the overall normalization of the correlator and any higher-order terms would only modify the high-$`k`$ end of the spectrum. However, they considered specifically the case in which the zero-temperature Higgs mass vanishes, i.e., $`m^2=0`$. In that case, the thermally generated effective mass for the Higgs field $`M^2`$, which is approximately equal to $`m_T^2`$, is always $`e^2T^2`$ \[see Eq. (3)\] and acts as an infrared regulator in the loop integral, making $`\mathrm{\Pi }_T(k)`$ an analytic function of $`k^2`$. Because one can show that the constant term is forbidden, the lowest-order term must be $`O(k^2)`$.

In our case, this argument breaks down because we are studying the system at the transition point where the effective Higgs mass $`M`$ becomes very small. In the momentum range $`M|\stackrel{}{k}|eT`$, the perturbative behaviour of $`\mathrm{\Pi }_T(k)`$ at one loop order is

$$\mathrm{\Pi }_T(k)=\frac{e^2T}{16}k=0.005625Tk.$$
(28)
This is shown as a solid line in Fig. LABEL:fig:spatcorrb, and agrees fairly well with the data. Assuming that perturbation theory is applicable, the discrepancy at very low $`k`$ can be explained by a non-zero effective Higgs mass, which arises because we are not exactly at the transition point. With $`M>0`$, the one-loop result is

$$\mathrm{\Pi }_T(k)=\frac{e^2MT}{4\pi }\left[\frac{4M^2+k^2}{2kM}\mathrm{arctan}\frac{k}{2M}1\right].$$
(29)
The best fit gives $`M=(0.0040\pm 0.0005)T`$, and is shown in Fig. LABEL:fig:spatcorrb as a short-dashed line. It is curious how well it agrees with the measurements, because one would expect perturbation theory to break down at low momenta near the phase transition.

Assuming that the Higgs mass vanishes at the transition point, all that remains is Eq. (28). The existence of a linear term in the self energy $`\mathrm{\Pi }_T`$ is surprising, because it implies magnetic screening. This screening is not as strong as it would be if there was a constant term, in which case the correlator would fall exponentially in the coordinate space. With a linear term, the low-momentum behaviour of the transverse gauge field correlator is $`k^1`$, and consequently, the long-distance behaviour in the coordinate space is $`r^2`$ instead of the usual $`r^1`$.

We also measured real-time correlators in the same simulations. In each of the nine independent configurations, we measured the correlator $`G(t,k)`$ and took the average of the results. Two examples, measured for $`k=0.026T`$ at $`m^2=0.083T^2`$ and $`m^2=0.086T^2`$, are shown in Fig. 2a. At all values of $`m^2`$ and $`k`$ we measured, the data were well described by the function

$$G_{\mathrm{fit}}(t)=a_0\mathrm{exp}(\gamma _Lt)+a_1\mathrm{exp}(\gamma _pt)\mathrm{cos}(\omega _pt+\delta ),$$
(30)
where $`a_0`$, $`a_1`$, $`\gamma _L`$, $`\gamma _p`$, $`\omega _p`$ and $`\delta `$ are free parameters. Physically, $`\gamma _L`$ is the Landau damping rate, $`\gamma _p`$ is the plasmon damping rate, and $`\omega _p`$ is the plasmon frequency. We fitted this function to the data at each point and estimated the errors using the bootstrap method. The results are shown in Figs. 2b–d.

We can compare these results with perturbative calculations, but we must keep in mind that since these simulations were carried out in a classical lattice theory, the result is not the same as in quantum theory in continuum. In the symmetric phase, the plasmon frequency should behave as

$$\omega _p(\stackrel{}{k})=\sqrt{\stackrel{}{k}^2+m_p^2},$$
(31)
where $`m_p`$ is the plasmon mass, which has been calculated in classical lattice perturbation theory in Refs. ,

$$m_p^20.086e^2\frac{T}{\delta x}0.0013T^2.$$
(32)
We have plotted this curve in Fig. 2b, and it agrees very nicely with the measured frequencies at $`m^2=0.083T^2`$. In the Higgs phase, the photon becomes massive due to the Higgs mechanism and this increases $`m_p`$.

The plasmon damping rate has not been calculated perturbatively for the classical lattice theory. However, a calculation has been carried out in the HTL-improved theory by Evans and Pearson , who showed that it was peaked just below the phase transition. Our classical results also show a rising trend as the transition is approached, although we cannot directly compare the numerical values.

The Landau damping rate in the symmetric phase has been calculated perturbatively in Ref. ,

$$\gamma _L\frac{7\delta x}{e^2T}k^3470\frac{k^3}{T^2}.$$
(33)
However, in our results the dependence on $`\stackrel{}{k}`$ is milder and we find the best fit with $`\gamma _L15T^{1.1}k^{2.1}`$ above the transition point. Below the transition point, the exponential decay rate $`\gamma _L`$ becomes very large so that at late times, the correlator simply oscillates around zero (see Fig. 2a). In Fig. 2d, this behaviour can be seen in the values of the Landau damping rate.

In the defect formation scenario, the freeze-out occurs when Eq. (18) ceases to be satisfied. This happens most easily at the transition point, and there the real-time correlators are still the same as in the symmetric phase. Therefore, it is the symmetric phase correlators that determine the relevant time scale $`\tau (k)`$, and Fig. 2 shows clearly that for the long-wavelength modes, the lowest time scale is that of Landau damping, $`\tau (k)\gamma _L(k)^1`$.

Substituting $`\tau (k)=\gamma _L(k)^1`$ to Eq. (18), we find

$$\frac{1}{\widehat{k}^2}\frac{\delta m^2}{\tau _Q}\gamma _L15T^{1.1}\widehat{k}^{2.1},$$
(34)
which implies

$$\widehat{k}\left(\frac{0.0061}{T\tau _Q}\right)^{1/4.1}T0.29T(T\tau _Q)^{0.24}\tau _Q^{0.24}.$$
(35)
### B Phase transition

We simulated the phase transition by thermalizing a number of configurations with $`m^2=m_0^2=0.044T^2`$, and solving numerically the equations of motion (23), varying the mass term with time according to

$$m^2(t)=m_0^2\delta m^2\left(\frac{4}{3\pi }\mathrm{arctan}(t/\tau _Q1)+\frac{1}{3}\right),$$
(36)
where $`\delta m^2=0.089T^2`$. This form of $`m^2(t)`$ has the advantage that long time after the transition, the system is in thermal equilibrium, making it possible to compare the final states.

The vortices produced in the transition are closed loops, and in general they will soon shrink to a point and disappear, which makes it very difficult to even define what we mean by the final vortex number. However, on a periodic lattice some of the vortex loops can be non-contractible, i.e. wind around the lattice in some direction. These loops can only disappear if they annihilate with another vortex of opposite direction, but this process is very slow, since the interactions between the vortices are exponentially suppressed at long distances. Therefore we chose one of the lattice dimensions much shorter than the other two (but still much longer than the microscopic length scales such as the Debye screening length $`m_D^1`$.) A long time after the transition, when the system is deep in the broken phase, these vortices still remain and are well-defined macroscopic objects, and they can be counted without any ambiguities using the gauge-invariant lattice definition for the winding number .

Because the scenario of defect formation discussed here is based on the assumption that the distribution of the magnetic field freezes in the transition, we can carry out a very simple test for the scenario by quenching the system through the transition and measuring the spatial correlator. We used $`\tau _Q=300T^1`$ and stopped the quench at $`m^2=0.089T^2`$ to carry out the measurement. The spatial correlator is shown in Fig. LABEL:fig:spatcorra as gray squares, and indeed, resembles much more the symmetric phase correlator than the equilibrium correlator at the same value of $`m^2`$. We also measured the real-time correlator, and the corresponding time scales are shown in Figs. 2b–c. The plasmon frequency and decay rate do not differ significantly from their equilibrium values, but Landau damping gets extremely slow. This means that once a mode has fallen out of equilibrium at the transition point, it takes a very long time before it thermalizes.

When one of the dimensions is shorter than $`\widehat{\xi }=2\pi /\widehat{k}`$, the prediction (20) changes into

$$n\frac{e}{2\pi }T^{1/2}L_z^{1/2}\widehat{k}^1.$$
(37)
We tested this prediction by simulating the time evolution with different values of $`\tau _Q`$ and measuring the vortex number a long time after the transition at $`t=\tau _Q+2400T^1`$. The results were published in Ref. , and are shown in Fig. 3. Each point is an average over around 15 runs starting from different initial configurations.

Combining the result $`\widehat{k}\tau _Q^{0.24}`$ from Eq. (35) with Eq. (37), we expect $`Nc\tau _Q^{0.24}`$. If we determine the constant $`c`$ from the best fit to the data at $`\tau _Q>200T^1`$, we find $`c=39.5\pm 1.0`$ and $`\chi ^2=7.6/13`$ dof. This is shown in Fig. 3 as the dashed line. If we also leave the exponent as a free parameter, we find $`N=(43.9\pm 7.9)(T\tau _Q)^{0.255\pm 0.026}`$, with $`\chi ^2=7.2/12`$ dof.

If we calculate precisely the prediction (37) using Eq. (35), we find $`c650`$ for the above fit parameter. However, it is not surprising that it differs from the measured value by a factor of $`15`$, not only because have we neglected factors or $`2\pi `$ and other numerical factors, but also because in our estimates we have assumed that the mechanism is ideally efficient and that all the magnetic flux is converted into vortices.

From the cosmological point of view, we are more interested in the area density of vortices in a fully three-dimensional case than in a thin box. Because the topology of the system does not prevent vortex loops from shrinking, the resulting network is not stable, unless it is stabilized by the expansion of the universe. We can therefore only use Eq. (20) as an estimate for the area density of vortices immediately after the transition. Comparing Eqs. (20) and (37), we find

$$n_{3\mathrm{D}}=\left(\frac{2\pi }{e}\right)^{1/2}T^{1/4}L_z^{3/4}n_{2\mathrm{D}}^{3/2},$$
(38)
where $`n_{2\mathrm{D}}`$ is given by Eq. (37) and $`n_{3\mathrm{D}}`$ is the three-dimensional area density given by Eq. (20). In our case, $`L_z=120T^1`$ and $`A=5.2\times 10^5T^2`$, and we find

$$n_{3\mathrm{D}}166T^1n_{2\mathrm{D}}^{3/2}1.3\times 10^4(T\tau _Q)^{0.38}.$$
(39)
## V HTL simulations

### A Equilibrium

As discussed in Section II, the classical theory discussed above does not describe the dynamics of the quantum theory correctly. Therefore we also carried out the same simulations with the HTL improved theory (7). In the same way as in Section IV A, we first studied the equilibrium properties of the theory. We used the same couplings $`e=0.3`$ and $`\lambda =0.18`$ as in the classical case, and the same lattice. With these parameters, the Debye mass has the value $`m_D^2=0.03T^2`$. The number of Legendre modes was $`N_{\mathrm{max}}=4`$. The only effect of the HTL corrections to the thermodynamics is to give an extra contribution to electric screening, so we expect that the phase diagrams of the theories are practically the same. Therefore we only carried out the equilibrium measurements at the transition point, $`m^2=0.083T^2`$.

Using the Hamiltonian (A24), we first prepared the thermal initial conditions with a Monte Carlo algorithm. This can be done in two steps:

1. Generate the soft mode configuration. The part of the Hamiltonian that involves the hard modes is Gaussian and therefore they can be integrated out exactly. This results in a theory with only the soft fields $`\varphi `$ and $`A_i`$ and their canonical momenta. The Gauss law leads to an extra Debye screening term

$$\delta H_{\mathrm{Debye}}=\frac{1}{2m_D^2}\left(\stackrel{}{}\stackrel{}{E}2e\mathrm{Im}\varphi ^{}\pi \right).$$
(40)
in the Hamiltonian.<sup>*</sup><sup>*</sup>*In principle, the canonical momenta could also be integrated out at this stage, resulting in a theory with an extra neutral scalar field $`A_0`$. This would lead to an effective theory that is equivalent with dimensional reduction at one-loop order . We used a Metropolis algorithm for $`\varphi `$ and a heat bath algorithm for all other fields, and carried out around 10000 thermalization sweeps.
2. Generate the hard modes in the background of the soft modes generated in step 1. Since the Hamiltonian is Gaussian it could in principle be diagonalized and the field values could be taken directly from a normal distribution. However, we did this only for $`\theta `$ and its momentum $`\mathrm{\Pi }`$. Note that the lowest Legendre mode of $`\mathrm{\Pi }`$ is fixed by the Gauss law,

$$\mathrm{\Pi }^{(0)}=\frac{1}{m_D}\left(\stackrel{}{}\stackrel{}{E}2e\mathrm{Im}\varphi ^{}\pi \right).$$
(41)
The field $`\stackrel{}{f}`$ and its momentum $`\stackrel{}{F}`$ were generated using a heat bath algorithm with around 5000 sweeps.

Again, we evolved the configurations taken from the thermal ensemble using the equations of motion (A17) for the time $`\mathrm{\Delta }t=12000T^1`$, measuring the equal-time and real-time correlators. The equal-time correlator is shown as open circles in Fig. LABEL:fig:spatcorr, and it agrees fairly well with the corresponding classical correlator. This time, the fit to the perturbative one-loop result (29) favours $`M=0`$, which suggests that because of the presence of the hard modes, the transition point is at slightly larger $`m^2`$ than in the classical case, and that $`m^2=0.083T^2`$ is very close to it.

Because perturbation theory cannot be trusted at low momenta, we do not adopt the perturbative result (28), but instead we simply assume that the self energy is linear at small momenta

$$\mathrm{\Pi }_T(k)=k_{\mathrm{np}}k,$$
(42)
and determine the coefficient $`k_{\mathrm{np}}`$ from the best fit to the data, which gives

$$k_{\mathrm{np}}=(0.0071\pm 0.0005)T,$$
(43)
and use this in our later estimates.

The hard modes have a significant effect on the real-time correlators, shown in Figs. 2b–d. Because they mimic the effect of the hard modes in the continuum quantum theory, we should now be able to compare the results with the standard perturbative calculations. For the plasmon mass, the perturbative result is

$$m_p=\frac{1}{\sqrt{3}}m_D=0.1T.$$
(44)
The dashed line in Fig. 2b shows the corresponding curve $`\omega _p=\sqrt{k^2+m_p^2}`$. The measured values are slightly below this curve, but the agreement is still very good.

The continuum plasmon damping rate at zero momentum has been computed perturbatively in the HTL approximation in Refs. ,

$$\gamma _p(T)=\frac{e^2T}{24\pi }A(T/T_c,\lambda /e^2).$$
(45)
The function $`A(T/T_c,\lambda /e^2)`$ takes the value 1 at the critical point, peaks at value $`1`$ just below it, and vanishes above it when $`m_p>2M`$. In our case, the damping rate at the transition point would be

$$\gamma _p(T)0.0012T,$$
(46)
which is shown in Fig. 2c as a dashed line and is lower than value the measured at $`m^2=0.083T^2`$ by roughly a factor of three. It has been observed earlier , that also in the SU(2) theory, the perturbative calculation underestimates the plasmon damping rate significantly. We can also note that the damping rate is lower than in the classical theory.

The Landau damping rate can be obtained directly from the HTL improved Lagrangian and is

$$\gamma _L=\frac{4k^2}{\pi m_D^2}(k+k_{\mathrm{np}})42.4\frac{k^2}{T^2}(k+k_{\mathrm{np}}),$$
(47)
where we have taken into account the linear contribution (42) to the self energy. This curve is shown as the dashed line in Fig. 2d, and agrees well with the measured values. Again, we can conclude that the relevant time scale is that of Landau damping. If we can ignore the linear correction, we find \[cf. Eq. (35)\]

$$\widehat{k}\left(\frac{\pi }{4}m_D^2\frac{\delta m^2}{\tau _Q}\right)^{0.2}0.29T(T\tau _Q)^{0.2}\tau _Q^{0.2},$$
(48)
and in slowest transitions, we have

$$\widehat{k}\left(\frac{\pi }{4}\frac{m_D^2}{k_{\mathrm{np}}}\frac{\delta m^2}{\tau _Q}\right)^{0.25}0.74T(T\tau _Q)^{0.25}\tau _Q^{0.25}.$$
(49)
However, this corresponds to $`\tau _Q10^8T^1`$, which is much larger than the values of $`\tau _Q`$ we are able to use in out simulations.

### B Phase transition

As in Sect. IV B, we studied the non-equilibrium dynamics of the phase transition by starting from thermal configurations at $`m^2=0.044T^2`$ and evolving the system with a time-dependent mass term (36). We used two different values for $`N_{\mathrm{max}}`$, $`4`$ and $`16`$.

In Fig. 3, we show how the final vortex number, measured at $`t=2\tau _Q+2400T^1`$ depends on the quench rate $`\tau _Q`$. Each data point is an average over 20–30 runs. In fast transitions, $`\tau _Q<1000T^1`$, the results for $`N_{\mathrm{max}}=4`$ and $`16`$ agree, but in slower transitions there is a statistically significant difference. This can be understood in terms of Eq. (10): For small $`\tau _Q`$, $`N_{\mathrm{max}}=4`$ is sufficient because by the time the approximation breaks down, the system is already so deep in the Higgs phase that the vortex number cannot change any more. When $`\tau _Q`$ gets larger, eventually a point is reached at which the breakdown occurs so early that it would have an effect on the final state. This may have happened even for $`N_{\mathrm{max}}=16`$ in the slowest transitions with $`\tau _Q=6000T^1`$.

Combining Eq. (48) with Eq. (37), we find $`N\tau _Q^{0.2}`$. A fit to the $`N_{\mathrm{max}}=16`$ data at $`\tau _Q>100T^1`$ with $`N=c(T\tau _Q)^{0.2}`$ gives $`c=45.5\pm 0.9`$ with $`\chi ^2=6.8/6`$ dof, which shows that the results are compatible with the prediction. The fit is shown in Fig. 3 as a solid line. The prediction of Eq. (37) is $`c660`$, which is greater than the measured value by a factor of $`15`$. This is the same factor found in the classical case, which further supports our scenario. If we keep the exponent as a free parameter, we find $`N=(45.9\pm 4.6)(T\tau _Q)^{0.201\pm 0.015}`$, with $`\chi ^2=6.8/5`$ dof.

Again, we relate the results to three dimensions using Eq. (38), and find

$$n_{3\mathrm{D}}1.4\times 10^4(\tau _QT)^{0.30}.$$
(50)
The conjectured non-perturbative behaviour at low momenta (42) would change this scaling law in very slow transitions. Eq. (20) would become $`n\widehat{k}^2`$, and Eq. (49) would therefore imply $`n\tau _Q^{0.5}`$.

## VI Conclusions

In this paper we have presented a thorough study of the dynamics of the high temperature phase transition in the Abelian Higgs model, using both classical and HTL improved approximations. The aims were twofold. The first was to measure the equilibrium properties of the theory, and in particular find the relevant timescale for the decay of the long wavelength modes of the gauge field, which were identified in  as crucial to the understanding of vortex formation.

In Ref. , we performed numerical simulations of the Abelian Higgs model phase transition in real time by using the classical theory and changing the mass parameter of the Higgs field over a characteristic quench time $`\tau _Q`$. In this paper, our second aim was to do the same quenches with the HTL improved theory, and to compare the resulting scaling law for the number of vortices $`N`$ formed as a function of the quench time.

Our measurements of the equilibrium correlators show that perturbation theory gives a reasonable account of their behaviour, except perhaps for the plasmon damping rate. A new and unexpected result is that near the transition, the equal-time correlator exhibits power-law magnetic screening, with a coefficient that is similar but not equal to the perturbative one. In the HTL improved simulations, the numerical values of the plasmon mass and the Landau damping rate agree well with the perturbative values in the Coulomb phase. For the plasmon damping rate, however, the measured value was significantly higher than the perturbative estimate. As might be expected, the agreement is not as good in the classical theory: although the plasmon mass agrees well, the dependence of the Landau damping rate $`\gamma _L`$ on wavenumber $`k`$ is $`\gamma _Lk^{2.1}`$ rather than the expected $`k^3`$.

We found significant differences between the scaling laws for the classical simulations, and for HTL improved quenches with $`N_{\mathrm{max}}=4`$ and $`N_{\mathrm{max}}=16`$. The difference between the classical scaling law and the HTL improved one with $`N_{\mathrm{max}}=16`$ can be ascribed to the discrepancy in the Landau damping rate, and lends force to the contention made in  that it is the balancing of the cooling rate with the Landau damping rate which decides the length scale above which the fields fall out of equilibrium. The difference between $`N_{\mathrm{max}}=4`$ and $`N_{\mathrm{max}}=16`$ can be understood as stemming from a lack of phase space in the hard modes, which causes the HTL approximation to break down at times greater than $`t_c(k)=4N_{\mathrm{max}}/k`$.

Our simulations have been carried out to leading order in the couplings $`e`$ and $`\lambda `$, and hence do not include the effect of high momentum transfer scattering between the hard modes. However, the hard modes can still scatter by exchanging soft modes, we do not expect this to change the dynamics qualitatively on the relatively short time scales which we have been able to study. In very slow transitions, it may be that hard mode scattering and also non-perturbative effects in the photon self-energy start to become important and and the predicted scaling law $`n\tau _Q^{0.3}`$ ceases to be valid.

###### Acknowledgements.

We would like to thank Tim Evans, Edmond Iancu and Mikko Laine for useful discussions. Part of this work was conducted on the SGI Origin platform using COSMOS Consortium facilities, funded by HEFCE, PPARC and SGI. We also acknowledge computing support from the Sussex High Performance Computing Initiative. AR was supported by PPARC and also partly by the University of Helsinki.

## A HTL improved equations of motion

In Sect. II the local formulation of the HTL improved Abelian Higgs model was described. In this Appendix we detail the resulting equations of motion for the fields $`\varphi `$ and $`\stackrel{}{A}`$ and Legendre modes $`\theta `$ and $`\stackrel{}{f}`$, which encode the effect of high momentum ($`k>T`$) particles.

These fields satisfy the equations of motion

$`_0^2\stackrel{}{A}`$ $`=`$ $`\stackrel{}{}\times \stackrel{}{}\times \stackrel{}{A}2e\mathrm{Im}\varphi ^{}\stackrel{}{D}\varphi `$ (A2)
$`+m_D{\displaystyle _0^1}𝑑zz^2\left(\stackrel{}{}\theta m_D\stackrel{}{A}+\sqrt{{\displaystyle \frac{1z^2}{2z^2}}}\stackrel{}{}\times \stackrel{}{f}\right),`$
$`_0^2\stackrel{}{f}(z)`$ $`=`$ $`z^2\stackrel{}{}\times \stackrel{}{}\times \stackrel{}{f}+m_Dz\sqrt{{\displaystyle \frac{1z^2}{2}}}\stackrel{}{}\times \stackrel{}{A},`$ (A3)
$`_0^2\theta (z)`$ $`=`$ $`z^2\stackrel{}{}\left(\stackrel{}{}\theta m_D\stackrel{}{A}\right).`$ (A4)
Since these equations for $`\stackrel{}{f}`$ and $`\theta `$ are linear, it is easy to solve them and show that the dynamics of $`\varphi `$ and $`A_i`$ is identical to the original non-local theory (9).

Not only is this reformulation of the theory local, but the equations of motion are in a canonical form and we can therefore write down the corresponding Hamiltonian

$`H`$ $`=`$ $`{\displaystyle }d^3x{\displaystyle _0^1}dz[{\displaystyle \frac{1}{2}}\stackrel{}{E}^2+{\displaystyle \frac{1}{2}}(\stackrel{}{}\times \stackrel{}{A})^2+\pi ^{}\pi +(D_i\varphi )^{}(D_i\varphi )+m_T^2\varphi ^{}\varphi +\lambda (\varphi ^{}\varphi )^2`$ (A6)
$`+{\displaystyle \frac{1}{2}}\stackrel{}{F}^2+{\displaystyle \frac{1}{2}}\mathrm{\Pi }^2+{\displaystyle \frac{z^2}{2}}(\stackrel{}{}\times \stackrel{}{f})^2+{\displaystyle \frac{z^2}{2}}(\stackrel{}{}\theta m_D\stackrel{}{A})^2m_Dz\sqrt{{\displaystyle \frac{1z^2}{2}}}\stackrel{}{f}\stackrel{}{}\times \stackrel{}{A}],`$
where $`\stackrel{}{F}=_0\stackrel{}{f}`$ and $`\mathrm{\Pi }=_0\theta `$ are the canonical momenta of $`\stackrel{}{f}`$ and $`\theta `$, respectively. We also need two extra conditions, namely the transverseness of $`\stackrel{}{f}`$ and Gauss’s law

$`\stackrel{}{}\stackrel{}{f}`$ $`=`$ $`\stackrel{}{}\stackrel{}{F}=0,`$ (A7)
$`\stackrel{}{}\stackrel{}{E}`$ $`=`$ $`m_D{\displaystyle _0^1}𝑑z\mathrm{\Pi }(z)+2e\mathrm{Im}\varphi ^{}\pi .`$ (A8)
The $`z`$ dependence of $`\stackrel{}{f}`$ and $`\theta `$ is discretized  by introducing the Legendre modes

$`\stackrel{}{f}^{(n)}`$ $`=`$ $`{\displaystyle _0^1}𝑑zz\sqrt{{\displaystyle \frac{2}{1z^2}}}P_{2n}(z)\stackrel{}{f}(z),\theta ^{(n)}={\displaystyle _0^1}𝑑zP_{2n}(z)\theta (z),`$ (A9)
$`\stackrel{}{F}^{(n)}`$ $`=`$ $`{\displaystyle _0^1}{\displaystyle \frac{dz}{z}}\sqrt{{\displaystyle \frac{2}{1z^2}}}P_{2n}(z)\stackrel{}{F}(z),\mathrm{\Pi }^{(n)}={\displaystyle _0^1}𝑑zP_{2n}(z)\mathrm{\Pi }(z).`$ (A10)
Note that we have used slightly different definitions from Ref. , in order to write the Hamiltonian in a practical form.

In terms of these modes, the equations of motion become

$`_0\stackrel{}{A}`$ $`=`$ $`\stackrel{}{E},`$ (A11)
$`_0\stackrel{}{f}^{(n)}`$ $`=`$ $`C_n^+\stackrel{}{F}^{(n+1)}+C_n^0\stackrel{}{F}^{(n)}+C_n^{}\stackrel{}{F}^{(n1)},`$ (A12)
$`_0\theta ^{(n)}`$ $`=`$ $`\mathrm{\Pi }^{(n)},`$ (A13)
$`_0\stackrel{}{E}`$ $`=`$ $`\stackrel{}{}\times \stackrel{}{}\times \stackrel{}{A}+2e\mathrm{Im}\varphi ^{}\stackrel{}{D}\varphi +{\displaystyle \frac{1}{3}}m_D^2\stackrel{}{A}`$ (A15)
$`{\displaystyle \frac{1}{3}}m_D\left(\stackrel{}{}\theta ^{(0)}+2\stackrel{}{}\theta ^{(1)}+\stackrel{}{}\times \stackrel{}{f}^{(0)}\stackrel{}{}\times \stackrel{}{f}^{(1)}\right),`$
$`_0\stackrel{}{F}^{(n)}`$ $`=`$ $`\stackrel{}{}\times \stackrel{}{}\times \stackrel{}{f}^{(n)}+m_D\stackrel{}{}\times \stackrel{}{A}\delta _{n,0},`$ (A16)
$`_0\mathrm{\Pi }^{(n)}`$ $`=`$ $`C_n^+\stackrel{}{}^2\theta ^{(n+1)}+C_n^0\stackrel{}{}^2\theta ^{(n)}+C_n^{}\stackrel{}{}^2\theta ^{(n1)}m_D\stackrel{}{}\stackrel{}{A}\left({\displaystyle \frac{1}{3}}\delta _{n,0}+{\displaystyle \frac{2}{15}}\delta _{n,1}\right),`$ (A17)
where

$$C_n^+=\frac{(2n+1)(2n+2)}{(4n+1)(4n+3)},C_n^0=\frac{1}{4n+1}\left(\frac{(2n+1)^2}{4n+3}+\frac{4n^2}{4n1}\right),C_n^{}=\frac{2n(2n1)}{(4n+1)(4n1)}.$$
(A18)
The Hamiltonian becomes

$`H`$ $`=`$ $`{\displaystyle }d^3x\{{\displaystyle \frac{1}{2}}\stackrel{}{E}^2+{\displaystyle \frac{1}{2}}(\stackrel{}{}\times \stackrel{}{A})^2+\pi ^{}\pi +(D_i\varphi )^{}(D_i\varphi )+m_T^2\varphi ^{}\varphi +\lambda (\varphi ^{}\varphi )^2`$ (A24)
$`+{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}[{\displaystyle \frac{1}{4}}{\displaystyle \frac{1}{4n+3}}((2n+1)\stackrel{}{F}^{(n)}+(2n+2)\stackrel{}{F}^{(n+1)})^2`$
$`{\displaystyle \frac{4n+1}{4}}\left(C_n^+\stackrel{}{F}^{(n+1)}+C_n^0\stackrel{}{F}^{(n)}+C_n^{}\stackrel{}{F}^{(n1)}\right)^2`$
$`+{\displaystyle \frac{4n+1}{4}}\left(\stackrel{}{}\times \stackrel{}{f}^{(n)}\right)^2{\displaystyle \frac{1}{4}}{\displaystyle \frac{1}{4n+3}}\left((2n+1)\stackrel{}{}\times \stackrel{}{f}^{(n)}+(2n+2)\stackrel{}{}\times \stackrel{}{f}^{(n+1)}\right)^2`$
$`+{\displaystyle \frac{4n+1}{2}}\left(\mathrm{\Pi }^{(n)}\right)^2+{\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{4n+3}}\left((2n+1)\stackrel{}{}\theta ^{(n)}+(2n+2)\stackrel{}{}\theta ^{(n+1)}\right)^2`$
$`+{\displaystyle \frac{m_D}{3}}\stackrel{}{}\stackrel{}{A}(\theta ^{(0)}+2\theta ^{(1)}){\displaystyle \frac{m_D}{3}}\stackrel{}{}\times \stackrel{}{A}(\stackrel{}{f}^{(0)}\stackrel{}{f}^{(1)})+{\displaystyle \frac{1}{6}}m_D^2\stackrel{}{A}^2]\},`$
## B Lattice discretization

### 1 Classical theory

In order to carry out numerical simulations described in Section IV, we discretize the Hamiltonian (21) and the equations of motion (23) in the standard leap-frog fashion. The Higgs field $`\varphi `$ was defined at the lattice sites, and its canonical momentum at temporal links between time slices. The gauge field $`\stackrel{}{A}`$ was represented by real numbers defined at links between lattice sites, and the electric field $`\stackrel{}{E}`$ by temporal plaquettes that connect these links. Therefore $`\pi _{(t,\stackrel{}{x})}`$ is actually defined at the point $`(t+\delta t/2,\stackrel{}{x})`$, $`A_{i,(t,\stackrel{}{x})}`$ at $`(t,\stackrel{}{x}+\widehat{ı}/2)`$ and $`E_{i,(t,\stackrel{}{x})}`$ at $`(t+\delta t/2,\stackrel{}{x}+\widehat{ı}/2)`$. Here $`\widehat{ı}`$ is a vector of length $`\delta x`$ in the $`i`$ direction.

The lattice version of the Hamiltonian (21) is

$`H`$ $`=`$ $`{\displaystyle \underset{\stackrel{}{x}}{}}\delta x^3[{\displaystyle \frac{1}{2}}{\displaystyle \underset{i}{}}E_i+{\displaystyle \frac{1}{2}}{\displaystyle \underset{i}{}}\left(ϵ_{ijk}\mathrm{\Delta }_j^+A_k\right)^2+\pi ^{}\pi `$ (B2)
$`{\displaystyle \frac{2}{\delta x^2}}{\displaystyle \underset{i}{}}\mathrm{Re}\varphi _{(\stackrel{}{x})}^{}U_{i,(\stackrel{}{x})}\varphi _{(\stackrel{}{x}+\widehat{ı})}+(m_T^2+{\displaystyle \frac{6}{\delta x^2}})|\varphi |^2+\lambda |\varphi |^4].`$
where

$`U_i`$ $`=`$ $`\mathrm{exp}\left(ie\delta xA_i\right),`$ (B3)
$`\mathrm{\Delta }_i^\pm \varphi _{(\stackrel{}{x})}`$ $`=`$ $`\pm \delta x^1\left(\varphi _{(\stackrel{}{x}\pm \widehat{ı})}\varphi _{(\stackrel{}{x})}\right).`$ (B4)
We also define the lattice version of the covariant derivative

$`D_i^+\varphi _{(\stackrel{}{x})}`$ $`=`$ $`\delta x^1\left(U_{i,(\stackrel{}{x})}\varphi _{(\stackrel{}{x}+\widehat{ı})}\varphi _{(\stackrel{}{x})}\right),`$ (B5)
$`D_i^{}\varphi _{(\stackrel{}{x})}`$ $`=`$ $`\delta x^1\left(\varphi _{(\stackrel{}{x})}U_{i,(\stackrel{}{x}\widehat{ı})}^{}\varphi _{(\stackrel{}{x}\widehat{ı})}\right).`$ (B6)
The value of the bare lattice mass $`m_T^2`$ is given by Eq. (3) and was chosen in such a way that the Hamiltonian (B2) describes the thermodynamics of the finite-temperature theory with renormalized mass $`m^2`$ correctly .

The discretized equations of motion are

$`\mathrm{\Delta }_tE_{i,(t,\stackrel{}{x})}`$ $`=`$ $`ϵ_{ijk}ϵ_{klm}\mathrm{\Delta }_j^{}\mathrm{\Delta }_l^+A_{i,(t,\stackrel{}{x})}2e\mathrm{Im}\varphi _{(t,\stackrel{}{x})}^{}D_i^+\varphi _{(t,\stackrel{}{x})},`$ (B7)
$`\mathrm{\Delta }_t\pi _{(t,\stackrel{}{x})}`$ $`=`$ $`D_i^{}D_i^+\varphi _{(t,\stackrel{}{x})}m_T^2\varphi _{(t,\stackrel{}{x})}2\lambda |\varphi _{(t,\stackrel{}{x})}|^2\varphi _{(t,\stackrel{}{x})},`$ (B8)
$`\mathrm{\Delta }_tA_{i,(t+\delta t,\stackrel{}{x})}`$ $`=`$ $`E_{i,(t,\stackrel{}{x})},`$ (B9)
$`\mathrm{\Delta }\varphi _{(t+\delta t,\stackrel{}{x})}`$ $`=`$ $`\pi _{(t,\stackrel{}{x})},`$ (B10)
where $`\mathrm{\Delta }_t\varphi _{(t)}=\delta t^1[\varphi _{(t)}\varphi _{(t\delta t)}]`$ etc.

The lattice version of the Gauss law (22) is

$$\underset{i}{}\mathrm{\Delta }_i^{}E_{i,(t,\stackrel{}{x})}=2e\mathrm{Im}\varphi _{(t,\stackrel{}{x})}^{}\pi _{(t,\stackrel{}{x})}.$$
(B11)
This is an extra constraint the initial field configuration must satisfy.

### 2 HTL theory

In the HTL simulations described in Section V, the soft modes were discretized in the same way as in the classical case. The extra field $`\theta `$ is defined at lattice sites and $`f_i`$ at the plaquettes. We denote by $`f_{i,(t,\stackrel{}{x})}`$ the field value at $`(t,\stackrel{}{x}+\widehat{x}/2+\widehat{y}/2+\widehat{z}/2\widehat{ı}/2)`$.

The lattice version of the HTL-improved Hamiltonian (A24) is

$$H_{\mathrm{HTL}}=H+\delta x^3\underset{\stackrel{}{x}}{}\left[_F+_f+_\mathrm{\Pi }+_\theta \right],$$
(B12)
where

$`_F`$ $`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{4}}[{\displaystyle \frac{1}{4n+3}}((2n+1)\stackrel{}{F}^{(n)}+(2n+2)\stackrel{}{F}^{(n+1)})^2`$ (B14)
$`(4n+1)(C_n^+\stackrel{}{F}^{(n+1)}+C_n^0\stackrel{}{F}^{(n)}+C_n^{}\stackrel{}{F}^{(n1)})^2],`$
$`_f`$ $`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}[{\displaystyle \frac{4n+1}{4}}\left(ϵ_{ijk}\mathrm{\Delta }_j^{}f_k^{(n)}\right)^2{\displaystyle \frac{1}{4}}{\displaystyle \frac{1}{4n+3}}((2n+1)ϵ_{ijk}\mathrm{\Delta }_j^{}f_k^{(n)}`$ (B16)
$`+(2n+2)ϵ_{ijk}\mathrm{\Delta }_j^{}f_k^{(n+1)})^2]{\displaystyle \frac{m_D}{3}}ϵ_{ijk}\mathrm{\Delta }_j^+A_k(f_k^{(0)}f_k^{(1)})+{\displaystyle \frac{1}{6}}m_D^2A_i^2,`$
$`_\mathrm{\Pi }`$ $`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{4n+1}{2}}\left(\mathrm{\Pi }^{(n)}\right)^2,`$ (B17)
$`_\theta `$ $`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{4n+3}}\left((2n+1)\mathrm{\Delta }_i^+\theta ^{(n)}+(2n+2)\mathrm{\Delta }_i^+\theta ^{(n+1)}\right)^2+{\displaystyle \frac{m_D}{3}}\mathrm{\Delta }_i^{}A_i\left(\theta ^{(0)}+2\theta ^{(1)}\right).`$ (B18)
The Gauss law (A8) can be written in the form

$$\mathrm{\Pi }_{(t,\stackrel{}{x})}^{(0)}=\frac{1}{m_D}\left(\underset{i}{}\mathrm{\Delta }_i^{}E_{i,(t,\stackrel{}{x})}2e\mathrm{Im}\varphi _{(t,\stackrel{}{x})}^{}\pi _{(t,\stackrel{}{x})}\right),$$
(B19)
and therefore we can eliminate $`\mathrm{\Pi }^{(0)}`$ from the Hamiltonian.

The bare Higgs mass has the value given in Eq. (3), and the Debye mass is

$$m_D^2=\frac{1}{3}e^2T^22e^2\frac{3.176T}{4\pi \delta x}.$$
(B20)
These counterterms were calculated by matching static correlators and in the absence of Lorenz invariance, they do not remove all the ultraviolet divergences from real-time quantities. However, since our lattice spacing is relatively large, this leads only to small errors.

The discretized equations of motion are

$`\mathrm{\Delta }_tE_{i,(t,\stackrel{}{x})}`$ $`=`$ $`ϵ_{ijk}ϵ_{klm}\mathrm{\Delta }_j^{}\mathrm{\Delta }_l^+A_{i,(t,\stackrel{}{x})}2e\mathrm{Im}\varphi _{(t,\stackrel{}{x})}^{}D^+\varphi _{i,(t,\stackrel{}{x})}`$ (B22)
$`{\displaystyle \frac{m_D}{3}}\left(\mathrm{\Delta }_i^+\theta _{(t,\stackrel{}{x})}^{(0)}+2\mathrm{\Delta }_i^+\theta _{(t,\stackrel{}{x})}^{(1)}+ϵ_{ijk}\mathrm{\Delta }_j^{}f_{k,(t,\stackrel{}{x})}^{(0)}ϵ_{ijk}\mathrm{\Delta }_j^{}f_{k,(t,\stackrel{}{x})}^{(1)}\right),`$
$`\mathrm{\Delta }_t\pi _{(t,\stackrel{}{x})}`$ $`=`$ $`D_i^{}D_i^+\varphi _{(t,\stackrel{}{x})}m_T^2\varphi _{(t,\stackrel{}{x})}2\lambda |\varphi _{(t,\stackrel{}{x})}|^2\varphi _{(t,\stackrel{}{x})},`$ (B23)
$`\mathrm{\Delta }_tF_{i,(t,\stackrel{}{x})}^{(n)}`$ $`=`$ $`ϵ_{ijk}ϵ_{klm}\mathrm{\Delta }_j^+\mathrm{\Delta }_l^{}f_{m,(t,\stackrel{}{x})}^{(n)}+\delta _{n,0}m_Dϵ_{ijk}\mathrm{\Delta }_j^+A_{k,(t,\stackrel{}{x})},`$ (B24)
$`\mathrm{\Delta }_t\mathrm{\Pi }_{(t,\stackrel{}{x})}^{(n)}`$ $`=`$ $`C_n^+\mathrm{\Delta }_i^+\mathrm{\Delta }_i^{}\theta _{(t,\stackrel{}{x})}^{(n+1)}+C_n^0\mathrm{\Delta }_i^+\mathrm{\Delta }_i^{}\theta _{(t,\stackrel{}{x})}^{(n)}+C_n^{}\mathrm{\Delta }_i^+\mathrm{\Delta }_i^{}\theta _{(t,\stackrel{}{x})}^{(n1)}`$ (B26)
$`m_D\mathrm{\Delta }_i^{}A_{i,(t,\stackrel{}{x})}\left({\displaystyle \frac{1}{3}}\delta _{n,0}+{\displaystyle \frac{2}{15}}\delta _{n,1}\right),`$
$`\mathrm{\Delta }_tA_{i,(t+\delta t,\stackrel{}{x})}`$ $`=`$ $`E_{i,(t,\stackrel{}{x})},`$ (B27)
$`\mathrm{\Delta }_t\varphi _{(t+\delta t,\stackrel{}{x})}`$ $`=`$ $`\pi _{(t,\stackrel{}{x})},`$ (B28)
$`\mathrm{\Delta }_tf_{i,(t+\delta t,\stackrel{}{x})}^{(n)}`$ $`=`$ $`C_n^+F_{i,(t,\stackrel{}{x})}^{(n+1)}+C_n^0F_{i,(t,\stackrel{}{x})}^{(n)}+C_n^{}F_{i,(t,\stackrel{}{x})}^{(n1)},`$ (B29)
$`\mathrm{\Delta }_t\theta _{(t+\delta t,\stackrel{}{x})}^{(n)}`$ $`=`$ $`\mathrm{\Pi }_{(t,\stackrel{}{x})}^{(n)}.`$ (B30)