Spontaneous vortex formation on a superconducting film

We carry out numerical simulations of a temperature quench in a Hamiltonian system that models a superconducting film in three-dimensional (3D) space. Theoretical arguments suggest that vortices are spontaneously formed by two separate mechanisms: the Kibble–Zurek mechanism, which is based on the dynamics of the order parameter, and the flux trapping mechanism, which is based on the dynamics of the electromagnetic field. We find that fast quenches are in agreement with the predictions of the flux trapping scenario, but the results from slow quenches appear to point to the Kibble–Zurek mechanism. Our results demonstrate the crucial role the electromagnetic field plays in this phenomenon, making it very different from vortex formation in fully 2D systems or in superfluids.


Introduction
When a superconductor is rapidly cooled from the normal to the superconducting phase in the absence of an external magnetic field, vortices form. This phenomenon has been observed in several experiments [1]- [8] and has analogues in other condensed matter systems such as superfluids [9]- [11] and liquid crystals [12]- [15].
There has been considerable discussion in the literature of the mechanism by which vortices form and in particular of what controls their number density. In the case of superfluids and liquid crystals, the order parameter is electrically neutral, and defect formation can be understood in terms of the so-called Kibble-Zurek mechanism [16]- [18]. This was first proposed in connection with the possible formation of defects at phase transitions in the very early universe. When a system described by a complex order parameter ψ is cooled rapidly through its critical temperature T c , there is spontaneous breaking of the U(1) phase symmetry: ψ becomes nonzero and has to choose a random phase. In a large system undergoing a rapid phase transition, this choice is made independently in widely separated regions. Inevitably there will be loops around which the phase varies by 2π (or a multiple). Such a loop must be threaded by at least one vortex. So we expect a random tangle of vortices to form. One can then estimate the number density of vortices some time after the transition by comparing the rate of temperature quench with the rate at which phase information can propagate [18].
Analogous ideas have also been tested in a nonlinear optical system [19] and in a fluid undergoing the conduction-convection transition [20]. The situation is more complicated however for superconductors, because the Cooper pairs described by the order parameter ψ are electrically charged. In this case the electromagnetic field plays a crucial role. It has been shown that there is another distinct defect-formation mechanism involving fluctuations of the magnetic field rather than the scalar order-parameter field [21,22]: before the transition, when the system is still in the normal phase, there are long-wavelength thermal magnetic field fluctuations. In the 3 DEUTSCHE PHYSIKALISCHE GESELLSCHAFT superconducting phase, such magnetic fields are forbidden by the Meissner effect, but because of causality, they are unable to decay if the phase transition takes place in a finite time. Therefore, they survive even after the transition, and the only way they can do it is by forming vortices each of which carries one flux quantum. One way to understand how this happens is that in the presence of a nontrivial electromagnetic vector potential, the minimum of energy corresponds to a nonzero gradient of the complex phase of the order parameter. Magnetic field, as the curl of the vector potential, will therefore produce frustrations in the complex phase in the form of vortices.
This flux-trapping mechanism leads to a very different pattern of defect formation. In particular, in a thin film, if flux-trapping is dominant, one may expect to find clusters containing only vortices, or only antivortices, whereas the Kibble-Zurek mechanism would suggest a much more even distribution. Possible confirmation of this second mechanism has been provided by an experiment on thin-film superconducting rings [6].
Understanding the mechanisms by which vortices form will also have implications well beyond condensed matter physics. In particular, similar mechanisms may have produced topological defects such as cosmic strings, magnetic monopoles and domain walls in the early universe [16,23], either in cosmological phase transitions [17] or in brane collisions [24]. In cosmology, a neutral order parameter corresponds to a global symmetry and a charged order parameter to a local gauge symmetry. In the latter case, which is more common, the vortices are Nielsen-Olesen cosmic strings [25]. The role of the Cooper pairs is played by the Higgs field, which obtains a nonzero vacuum expectation value at the transition.
In this paper, we investigate vortex formation in the case of a two-dimensional (2D) film in 3D space, which is the typical set-up in actual superconductor experiments [1]. We demonstrate that the electromagnetic field plays a crucial role in this phenomenon, as it extends to the empty space outside the film and gives rise to long-range interactions. In the case of a neutral order parameter, a 2D film can be studied using a fully 2D simulation, but the 3D nature of the electromagnetic field means that the same is not true for superconductors [26].
Because our set-up captures some of the peculiarities of fully 3D systems, it is also more relevant for cosmological applications than previous 2D studies [27]. A different borderline case between two and three spatial dimensions has been studied in the past in [21,28]. This paper is organized as follows. In section 2, we describe the set-up of our scenario and the Hamiltonian approach we follow. In section 3, we describe the expected evolution of the two-point correlation function of the magnetic field G(k) and explain its role in the flux-trapping mechanism. In particular, we obtain predictions for the number density and spatial distribution of vortices. In section 4, we describe the details of numerical simulations we carried out to test these predictions. In section 5, we present the results of the simulations and compare them with the predictions of the flux-trapping scenario derived in section 3. We summarize our conclusions in section 6. A comparison of our approach with the time-dependent Ginzburg-Landau (TDGL) equations is included in appendix A. Appendix B describes the method we used to discretize the equations of motion in lattice simulations.

Set-up
The system we shall study consists of a 2D film in a 3D space. We consider a hypothetical experiment similar to the one reported in [1]. The whole system is initially in thermal equilibrium 4 DEUTSCHE PHYSIKALISCHE GESELLSCHAFT at a temperature above T c , so that the film is in the normal phase. Then the system is gradually cooled through the phase transition to the superconducting phase.After the transition, we examine the number and spatial distribution of the produced vortices.

Hamiltonian model
Numerical studies more or less like this have been carried out in the past, but they have generally ignored the electromagnetic field outside the film. It is also a common practice in theoretical studies of superconductor dynamics to use the TDGL approach [29]- [33], in which dissipative equations of motion are coupled to a noise term that models a thermal heat bath. In reality, the dynamics of the electromagnetic field outside the film is, of course, described by non-dissipative Maxwell's equations, and since we are interested in precisely the behaviour of the electromagnetic field, we have to take this into account. To be consistent, we will therefore use a toy model with Hamiltonian dynamics instead. The differences between these two approaches are discussed in appendix A.
Our starting point in designing our model is the Ginzburg-Landau free energy density, where ψ is the order parameter, A is the electromagnetic vector potential (gauge field) and B = ∇ × A is the magnetic field. The parameters e and m are the electric charge and mass of the Cooper pair, and α and β are phenomenological parameters. In the Ginzburg-Landau picture, α is temperature dependent and is negative in the superconducting phase. To simplify our expressions, we will work in natural units c =h = µ 0 = 0 = k B = 1 from now on. We now note that, in addition to realistic superconductors, the Ginzburg-Landau free energy (1) also describes the phase transition of a simple field theory with Hamiltonian (in the temporal gauge) where we have rescaled the field ψ → √ 2mψ and used the covariant derivative notation D = ∇ + ieA. We have also included the canonical momenta E = −∂A/∂t (which corresponds to the electric field) and π =ψ * in a Lorentz invariant way.
The thermal equilibrium state of theory (2) at temperature T is given by the canonical ensemble, and as usual, the free energy is The sum in this expression is over all possible field configurations. More explicitly, the partition function can be written as where the delta function imposes the Gauss law at microscopic level. The nonzero temperature breaks the Lorentz invariance, but it is important that the underlying model is still invariant because we are interested in effects that are due to a finite signal propagation speed.

DEUTSCHE PHYSIKALISCHE GESELLSCHAFT
As a function of temperature T , our model has two different phases, which correspond to the normal and superconducting phases, and a phase transition between them at some critical temperature T c . The order of the phase transition depends on the parameter β, but for our choice it is of second order. Because our model has the same symmetries as the Ginzburg-Landau free energy (1), one expects on general grounds that near the phase transition the free energy of our model is well approximated by an expression of the form (1).
We now have a relativistic microscopic model that has the same Ginzburg-Landau effective description as a real superconductor. Taking this point of view literally, we assume that the parameters in the Hamiltonian are temperature-independent. Using the Hamiltonian (2) we can derive the equations of motion for ψ and A. These equations are Lorentz invariant and conserve energy. The only place where the temperature enters is in the initial conditions, which are drawn from the canonical ensemble described by (4). In particular, there is no damping or noise term in the equations of motion, because we are describing the system at the microscopic level, and there is no need for them because the thermal effects that they are normally introduced to approximate are included explicitly.

Superconductor film
However, instead of the fully 3D system described by (2), we want to describe a superconductor film in 3D space. It is easy to modify the Hamiltonian to take this into account by setting the order parameter ψ to zero everywhere outside the film. Denoting the film thickness by l and ignoring the derivatives of ψ in the orthogonal direction, we find the Hamiltonian To distinguish them from three-vectors denoted by bold symbols X, we denote two-vectors on the film by symbols with an arrow on top X. The covariant derivative on the film is therefore Dψ = ( ∇ + ie A)ψ, where the 2D vector potential A(x, y) is obtained from A by restricting to the plane z = 0, and x-and y-directions. In other words, A(x, y) consists of the x and y components of A(x, y, 0). The equations of motion are obtained from the Hamiltonian (5) using Hamilton's equations, which giveψ where j is the electric current. Because the only charged degree of freedom is ψ, the current is confined on the plane z = 0, i.e. j = δ(z/ l)(j ψ x , j ψ y , 0), where the electric current due to the field ψ is j ψ = 2eImψ * Dψ. The latter one of equations (6) is, of course, equivalent to Maxwell's equations.
The equations of motion (6) conserve the electric charge defined as ρ = 2eImψ * ψ . If the Gauss law ∇ · E = ρ is satisfied by the initial conditions, it will remain satisfied at later times.
The parameters we used in the simulations are given in table 1. For comparison, the experiment in [1] was carried out using a YBa 2 Cu 3 O 7−δ film with thickness 100 nm. In natural units the Cooper pair charge is e ≈ 0.6. The penetration depth and correlation length of YBa 2 Cu 3 O 7−δ are λ 0 ≈ 150 nm and ξ 0 ≈ 2 nm [34]. In units of nanometres, we find that these correspond to α ≈ −0.25 and β ≈ 4000. One can see that our simulations correspond to a 6 DEUTSCHE PHYSIKALISCHE GESELLSCHAFT thinner film and a shorter penetration depth λ ≈ 4. The reason for this choice is that numerical computations are easier if all the characteristic length scales are comparable.

Quench
Our aim is to simulate the dynamics of our theory (5) starting in thermal equilibrium at T > T c and gradually cooling the system down so that it undergoes a phase transition to the superconducting phase. Because the equations of motion (6) are Hamiltonian, the total energy of the system is conserved. In particular, that means that if the system is initially in thermal equilibrium at a given temperature T , it will remain at the same temperature, i.e. the field configurations will have the same probability distribution (4) at later times.
To cool the system down, we therefore have to modify the equations of motion. There are various ways in which this can be done. For instance, in [21], the quadratic term α was slowly varying, becoming more and more negative, and this triggered the phase transition by changing the critical temperature T c .
In this work, we choose a simpler approach and modify the equations of motion by adding small damping terms to the equations of motion, The only purpose of these damping terms is to cool down the system, and they should not be confused with the damping terms in the TDGL approach, which are introduced to model thermal effects. In particular, we do not introduce a corresponding noise term. For simplicity, we choose the damping rates in both (7) and (8) to have equal values. This has the convenient consequence that the Gauss law remains exactly satisfied at all times, which means that we can still think of our model as a microscopic theory. As long as the explicit damping term σ introduced in (7) and (8) is much smaller than the effective damping generated by thermal effects, we do not expect them to have any other significant effects than cooling the system down.
As we will discuss in more detail in the next section, our model will remain close to thermal equilibrium but with a slowly decreasing effective temperature, provided that the damping rate σ is small enough. When the effective temperature reaches the critical value T c , the system undergoes a phase transition into the superconducting phase.

Flux trapping
We will now examine analytically how the flux-trapping scenario, which was outlined in the introduction section, takes place in our model, and make predictions for numerical simulations. In the scenario, the vortices are formed by the magnetic field orthogonal to the film, To predict the number and spatial distribution of the vortices, we will therefore need to know the statistical properties of B at the time of the transition.

Evolution of the magnetic field
Initially the film is in the normal phase, ψ ≈ 0, and therefore it does not have a significant effect on the dynamics of the electromagnetic field. If we can ignore the presence of the film, the dynamics is linear, we simply have free electromagnetic waves with a small damping term. The Fourier modes of the orthogonal magnetic field B with wavenumber k greater than σ will therefore oscillate with a decreasing amplitude Again, to the extent that the film can be ignored, the equilibrium distribution given by (4) is approximately Gaussian in the normal phase. The statistical properties of the magnetic field are therefore to be completely characterized by the two-point function where the brackets · · · indicate a statistical ensemble average. In thermal equilibrium at temperature T , a straightforward calculation gives [26] Equation (9) implies that G(k) ∝ exp(−σt), and together with (11) this means that modes with | k| σ retain the equilibrium distribution with an exponentially decreasing effective temperature When the temperature reaches the critical value T c , the order parameter ψ becomes nonzero and the dynamics becomes nonlinear. If the system were to stay in equilibrium, it would now be in the superconducting phase and therefore repel magnetic fields. Because of the nonzero temperature, the magnetic field would not be exactly zero, but its two-point function would be given by [26] G(k) ≈ Tk 2 where is a length scale, which can be thought of as a magnetic screening length. Equation (11) shows that = 0 in the normal phase, so it starts from zero at T c and grows exponentially as the temperature decreases. At any finite temperature, thermal fluctuations with wavelength less than are suppressed relative to the normal phase, but longer wavelengths are unaffected. This is in stark contrast with fully 2D or 3D systems in which longer wavelengths are suppressed first, and means that on a film there is no sharp transition. Instead, the apparent critical temperature depends on the length scale.
Even though it is impossible to solve the nonlinear equations of motion analytically, we can generally say that the amplitude of a Fourier mode B( k) with a given wavenumber k cannot decay arbitrarily fast. Moreover, the longer wavelengths react slower. For instance, if the dynamics of the long-wavelength fluctuations is diffusive, the fastest possible decay rate is γ max (k) = k 2 /D, where D is the diffusion constant. It is therefore unavoidable that if the transition takes place in a finite time, the longest wavelengths are too slow to react, and freeze out. This means that there is a critical wavenumber k c so that modes with k < k c still have approximately their initial 8 DEUTSCHE PHYSIKALISCHE GESELLSCHAFT amplitude after the transition. Schematically, one can write the final two-point function as where T c is the critical temperature and θ(x) is a step function. The value of k c is determined by the dynamics, but to keep the discussion general, we do not attempt to calculate it. The survival of the long-wavelength fluctuations means that at distances much shorter than the critical wavelength λ c ≈ 1/k c , there is effectively a uniform magnetic field. If λ c is longer than the Pearl length, i.e. the size of a vortex, and the total flux integrated over the area λ c is greater than one flux quantum, this magnetic field must form Abrikosov vortices [21], just like a uniform magnetic field does. This mechanism of vortex formation is called flux trapping.
The number density of vortices per unit area produced by flux trapping is given by the typical (root-mean-square) magnetic field strength divided by the flux quantum 0 = 2π/e [26]: where Using (13), we find Moreover, because the frozen-out magnetic field consists of modes with wavelengths longer than 2π/k c , there are clusters of size 2π/k c of equal-sign vortices. We can estimate the number of vortices in each cluster by assuming that they are disks of radius π/k c : For the estimate in (14) to be valid, N cl has to be greater than one. Otherwise, one will have to consider the coupled dynamics of both the electromagnetic field and the order parameter. It seems likely that vortex formation can then be described by the Kibble-Zurek mechanism [16]- [18], although the electromagnetic field may still play a role. We should also note that the friction term σ causes modes with | k| σ to freeze out even in the absence of any critical dynamics. The solution of the linearized equations of motion with thermal initial conditions gives

DEUTSCHE PHYSIKALISCHE GESELLSCHAFT
where k 2 = k 2 + k 2 z andσ 2 = σ 2 − 4k 2 , and Erfc is the complementary error integral. The second line is valid at late times when | k| σ. This gives rise to an apparent critical wavenumber k c = √ σ/2t. We are mainly interested in the physically more relevant case in which the freeze-out is caused by the critical dynamics. This restricts us to relatively low damping rates σ. Nevertheless, simulations with higher σ are also useful, because with an exact two-point function (18) we can test in high detail the second part of the flux-trapping scenario, namely how the magnetic field turns into vortices. This is discussed in the next subsection.

Flux quantization
Let us now assume that at the freeze-out, the two-point function is given by some function G(k), and calculate how many vortices should be formed. From (14) we know that In the unphysical friction-dominated case, we can combine the integrations in (18) and (19) into one integral, where the approximate equality is valid at late times. In principle, t should be chosen to be the time the vortices form. The number of vortices per cluster is approximately In the physically more relevant case of smaller σ, the two-point function G(k) cannot be calculated analytically. Instead of the sharp step function of (13), we now assume that it can be approximated by a Gaussian function, where T 2 and k 2 are constants, which correspond to the freeze-out temperature and the critical wavenumber, respectively. The freeze-out temperature T 2 is expected to be below T c , but not necessarily equal to it, since (12) shows that long-wavelength modes do not feel the Meissner effect until the screening length is comparable with their wavelength. This suggests that T 2 is determined by the temperature at which (T 2 ) = 1/k 2 . Using (22), (19) yields

DEUTSCHE PHYSIKALISCHE GESELLSCHAFT
In analogy with (17), the number of vortices per cluster is roughly and in principle this number should be much greater than one for (19), and consequently also these estimates, to be valid. In a numerical simulation, we can measure not only n but also the two-point function G(k) after the transition. Comparing the two gives us a strong test of the validity of the scenario. However, (18) and (22) are idealized descriptions of what the two-point function G(k) should look like at the time of the freeze-out, when the magnetic field is still smooth over long distances. When the vortices form, the magnetic flux is quantized and this changes the two-point function.
To compare the measured two-point function with predictions, we have to understand how it changes due to this effect.
Let us first assume that we have N vortices and N antivortices. Each of which carries one flux quantum, so that the magnetic field has a delta function at the centre of the vortex. Given the spatial distribution of the vortices, we ask what the two-point function of the magnetic field would be. The total magnetic field is where x ± i are the positions of the vortices and antivortices. The two-point function is where the two-vortex correlator G ij ( x − y) is and the brackets indicate integration over the positions x ± i . As long as i = j, G ij ( x − y) = G 0 ( x − y) is independent of i and j. On the other hand, for i = j, we have The two last terms give contributions of order N/A 2 , which we will ignore, but the two first terms give delta functions, so where A is the area of the film.
In total, we have in Fourier space,

DEUTSCHE PHYSIKALISCHE GESELLSCHAFT
The continuous-field limit is basically equivalent to taking N → ∞, and 2N/A = n is the number density of vortices, so we find where G cont (k) is the two-point function given by the continuous magnetic field. In a realistic vortex, the magnetic field is not completely localized, which means that the delta functions spread out over a finite distance. If we write the Fourier transform of the magnetic field profile of a vortex as G 1 (k), we have Thus, if the two-point function of the field at the time of the transition was G cont (k), it becomes (32) when the field gets trapped into vortices. Therefore we would expect to find (32) if we measure the two-point function after the transition. In principle, G 1 (k) could be calculated by finding the static vortex solution, but we simply parameterize it by an exponential, which seems to give a fairly good fit to the numerical data. We use the Gaussian ansatz in (22) for G cont (k), so that the whole two-point function is where T 1 should be proportional to the number of vortices and k 1 inversely proportional to the size of a vortex. The constants T 2 and k 2 have the same interpretation as in (22). We will use this form for fitting our numerical results. In particular, k 1 should not depend on how fast the quench was, because it is a property of the static vortex solution.
The linearized result in (18) is not of exactly the same form as (22), and therefore the parameters T 2 and k 2 will not be exactly the same as their theoretical values in the fast-quench limit, but they should still be close to them.

Simulations
We tested our predictions by solving the fully nonlinear equations of motion (B.5) numerically using the initial conditions given by the partition function (4). To do this, we defined the system on a 3D lattice, one slice of which corresponds to the film. The details of the lattice implementation are presented in appendix B.
In the preparation of the thermal initial conditions we employed the fact that the energy functional E tot in (B.7) is quadratic and diagonal in E and π. Therefore, their probability distribution is Gaussian. The only complication is the Gauss law (B.6), without which the field values at different points would be uncorrelated.
We started by generating random, uncorrelated values for the component of π that is parallel to ψ, i.e. π ψ = ψ(Reψ * π)/(ψ * ψ). This component does not appear in the Gauss law (B.6), and therefore it has the simple Gaussian probability distribution Then, we went through all plaquettes and considered a change of the electric field at all links around the plaquette by the same amount, Such a change does not change the divergence of E and therefore does not affect the Gauss law constraint. It would change the energy by an amount that is at most quadratic in , where the constants Q and L depend on the field values at the neighbouring links. Therefore the probability distribution for is Gaussian, It is straightforward to draw the value of from this distribution. Next we evolved the field configuration for time t th using the equations of motion (B.5), with zero conductivity. We repeated this cycle of randomization and evolution steps N th times, monitoring the evolution of the two-point function G(k). When it reached the equilibrium form given by the analytical calculation, we considered the system to have thermalized.
We then solved the equations of motion (B.5) with the initial conditions produced by the thermalization algorithm. We used a number of different values for the conductivity σ to test the dependence on the cooling rate. For each value of σ, we repeated the run several times using different initial conditions taken from the same thermal ensemble, carrying out measurements separately for each initial condition and averaging the results.
The values of the other parameters are shown in table 1. In terms of the zero-temperature coherence length and penetration depth, they correspond to ξ 0 = 2 and λ 0 = 2, which are very convenient values for the numerical simulations. The Ginzburg parameter has the value κ = 1, which means that our system is a type II superconductor.

Results
Our main aim was to test the flux-trapping theory, and therefore we measured the two-point function G(k) of the magnetic field fluctuations on the film. In figure 1, we show the time evolution of G(k) for σ = 1.0. One can see that, as expected, the long-wavelength modes freeze out to a high amplitude, whereas short wavelengths are exponentially suppressed. The plot also shows how the amplitude at around 0.3 k 0.8 increases after t ≈ 10, when the system enters the superconducting phase and the two-point function changes from (22) to (33) because of flux quantization.
To make it possible to compare the results for different values of σ, we chose to carry out our measurements at time t = 20/σ, so that the effective temperature is the same in each case. Figure 2, shows the measured two-point function at that time for selected values of σ. They are compared with the results for linearized friction-dominated freeze-out, calculated numerically from the first line of (18). The plot shows that the linearized result works very well in the fastest quench (σ = 5). For σ 1, one can see a clear discrepancy, which is a sign that the nonlinear effects due to the critical dynamics have become important. The difference is also partly due to the contribution G 1 (k) from flux quantization.
We fitted the two-point functions with the ansatz in (33). In figure 3, we show k 1 as a function of σ. Its value is independent of σ in slow quenches, as it should be if it is to be interpreted as the vortex size. At higher σ, the vortex density becomes higher and non-equilibrium effects more important, and therefore it is not surprising that k 1 starts to grow. This can be interpreted as the size of the vortex getting smaller as more of them are packed in the same area.
To improve the accuracy of the fit parameters in subsequent fits, we fixed k 1 to 0.1. Some examples of these fits are shown in figure 4. The plateau at very short wavelengths, which corresponds to modes still in thermal equilibrium, was excluded from the fits.
In figure 5, we show the fit parameters T 1 , T 2 and k 2 as functions of σ. The errors were estimated using the bootstrap method, and contain only the statistical error. We did not attempt to estimate the systematic error due to the choice of the fitting function.
In the (unphysical) high-σ limit, the two-point function should be given by (18), and to the extent that it can be approximated by a Gaussian, we expect T 2 ≈ T = 10 and k 2 ≈ √ σ/2t,  Figure 3. The exponential cutoff scale k 1 as a function of σ. For σ → 0 it approaches a constant value k 1 = 0.104(1), which characterizes the size of a static isolated vortex solution. At high σ, vortex density is high and the system still away from equilibrium, and the vortices are not well described by the static solution. In this regime, we fit k 1 = 0.0140(6)σ 1.60(3) . In all other fits, we fixed k 1 = 0.1. (filled circles) and T 2 (empty diamonds) of the exponential and Gaussian terms in (33) respectively. The dotted line shows the power-law fit T 1 = 9.76(9)σ 1.218 (6) to the first six data points. The dashed line shows a fit T 2 = 2.6(5)σ 0.43 (6) to the first four data points. (b) The critical wavenumber k 2 . The solid line shows the theoretical prediction k 2 = √ σ/2t, and the dashed and dotted lines are power-law fits k 2 = 0.186(5)σ 0.318 (8) and k 2 = 0.115(7)σ 0.95(3) of the first five and last six data points, respectively. which for our choice of t gives k 2 ≈ σ/ √ 40. This is confirmed by the measurements. In this limit the values of T 1 show significant scattering, but this is understandable, because G(k) is dominated by the Gaussian term and we had also fixed k 1 to a very different value from the best fit.
The deviation from the prediction k 2 ≈ √ σ/2t is a sign of nonlinear critical dynamics. Figure 6 shows the number of vortices plus antivortices measured at time σt = 20. The agreement with the analytical prediction (20) shown by the dashed line is good for fast quenches σ 1, apart from a constant factor of about 2. This indicates that vortices are predominantly formed by flux trapping.
For σ 1, we can fit the data very well with a power law N = 772(18)σ 1.20(2) . This is compatible with the expectation that T 1 ∝ N. We have also plotted the theoretical prediction in (23) calculated using the fitted parameter values. As the plot shows, it gives the correct order of magnitude, although it does not agree perfectly for slow quenches. The prediction seems to suggest a different power-law behaviour from the observed one, and may be indicating that in these slow quenches vortex formation is dominated by the Kibble-Zurek mechanism.
As figure 1 shows, the fluctuations have frozen out, but it is possible that their amplitude is not high enough to produce the characteristic signatures of flux trapping. Indeed, figure 7 shows the predicted number N cl of vortices in a cluster calculated using (24) and the fit parameters T 2 and k 2 . The predictions of the flux-trapping mechanism in section 3.2 assume it to be greater than one, but in slow quenches it is just below one. It is interesting to note that N cl calculated from T 2 and k 2 depends very weakly on the cooling rate σ. This suggests that the strength of the flux-trapping mechanism would not depend on the cooling rate but only on the parameters of the model, or in the case of real superconductors, the material and the thickness of the film. This is a consequence of the film geometry, since previous studies of fully 2 and 3D systems have suggested that flux trapping always dominates in slow quenches. The reason for this can be traced back to (12) and the fact that the effective critical temperature depends on the length scale, which makes T 2 strongly dependent on σ.
It is possible to increase the predicted N cl and therefore the expected strength of the fluxtrapping mechanism by using different parameter values, and we carried out one set of runs with T = 1000. The clustering is indeed evident in the snapshot of the magnetic field at time t = 400 in a run with σ = 0.25 shown in figure 8. In the same figure, we have also plotted the quantity, n C (r) which measures the average winding number around a circle of radius r centred at a vortex [21]. Clustering is indicated by values greater than one, and the maximum value gives a measure of the typical number of vortices in a cluster.

Conclusions
Fast quenches, with σ 1, are unphysical because the relevant degrees of freedom are overdamped even outside the film. However, because the freeze-out of the magnetic field can be calculated analytically, they allow us to see in detail how the frozen-out field turns into vortices. The number of produced vortices agrees well with the prediction of the flux-trapping theory. This provides strong support for the flux-trapping theory even in more realistic cases in which analytical calculation is not possible.
Although overdamped dynamics is often used in studies of superconductors and vortex formation [27], the underdamped case with smaller σ is more realistic. In fact, one could carry The curve in the inset shows the quantity n C (r), the average net number of vortices within a circle of radius r around a vortex, and confirms that the clustering that is evident in the field configuration is statistically significant. out the simulation without introducing an explicit damping term at all [21,28], because the nonlinear Hamiltonian dynamics gives rise to effective damping on the film. However, even in our case, the relevant wavelengths of the electromagnetic field are underdamped in slow quenches.
It is much more difficult to describe slow quenches theoretically, and therefore the predictions are not as robust. In spite of this, our results show a reasonably good agreement between theory and simulation, although the scaling of the vortex number with the cooling rate σ does not seem to agree. This may be an indication that the amplitude of trapped fluctuations is so low that vortex formation can be better described by the Kibble-Zurek mechanism. Our results also suggest that unlike in fully 2 or 3D set-ups, the flux-trapping mechanism does not seem to become more dominant in slow quenches, at least within the range of values we used.
Simulations with a higher initial temperature support this picture, and show clear vortex clustering even in slow quenches. This is an unmistakable sign of flux trapping, and demonstrates that it is the dominant mechanism in this case.
To understand better how vortices form when the amplitude is not high enough for flux trapping, it would be important to derive theoretical predictions of the Kibble-Zurek scenario and compare them with our results. In particular, one must understand how a weak background magnetic field biases the Kibble-Zurek mechanism.
Our results show that under the right conditions, the electromagnetic field plays an important role in defect formation, also in the experimentally relevant geometrical set-up of a 2D film in 3D space. It remains to be seen if these conditions are satisfied in actual superconductor experiments. Likewise, the cosmological consequences of this remain largely unexplored.

Appendix A. Comparison with TDGL
Our model is a physically consistent toy model with all the characteristic properties of a superconducting film, but it is not meant to be a quantitative model of any real superconducting material. It resembles in many ways the TDGL approach [29]- [33], which has become standard in non-equilibrium studies of superconductors, but it is not identical. TDGL is a phenomenological model of the macroscopic dynamics, whereas our model is defined in terms of microscopic degrees of freedom. This guarantees that our model is consistent with fundamental principles such as charge and energy conservation and special relativity and that our results are not dependent on any implicit assumptions about the dynamics. Nevertheless, as we show in this appendix, the two models have the same effective behaviour at long distances.
The full TDGL equation derived in [30] contains both first and second-order time derivatives, but it is common to drop the second-order terms to obtain a purely diffusive system of equations, (A.1) These equations have sometimes been used to study the dynamics of superconductor films. However, the dynamics of the electromagnetic field outside the film should be non-dissipative, and it would therefore be more realistic to use standard second-order Maxwell equations for the vector potential, as we have done in this paper.
In the TDGL approach, thermal fluctuations are introduced by adding noise and damping terms, which leads to a Langevin equation. At late times, the system then approaches a thermal equilibrium state, which is described by a partition function that resembles (4), with temperature determined by the amplitude of the noise. The only difference is that the partition function of TDGL, or indeed that of the TDGL approach to critical behaviour of superconductors [37,36], does not have the delta function δ (∇ · E − ρ) in the integrand. This delta function describes the long-range Coulomb interactions between electric charges, and the TDGL approach is therefore only valid at distances longer than the Debye screening length. In real superconductors, this is a good approximation, because the screening length is very short.
In our model, we have included the electric field, but it turns out to be screened by thermal fluctuations. In fact, the calculation presented in [26] for the magnetic field is just as valid for the screening of the electric field. We also checked this numerically. As shown in figure A.1, the electric charges are screened, and therefore our model has the same static properties at long distances as the TDGL model.
Before the addition of the small damping terms in (7) and (8), the microscopic dynamics of our model is non-dissipative, but the macroscopic dynamics is dissipative, just like in the TDGL model. Irrespectively of the initial configuration, the system will approach the equilibrium state (4) at late times. Any perturbation away from equilibrium will therefore disappear with time. This can be seen as the damping of the unequal-time correlation function of the charge density in figure A.2. Because the system stays in equilibrium, the fluctuation-dissipation relation is satisfied, with thermal noise arising from interactions with short-wavelength modes. There is, therefore, no need to add an explicit noise term to the equations of motion. . This can be interpreted as the total charge in a slab of width 2y around a unit charge at the origin. Without screening, this would be independent of y, but the data agrees with power-law screening of the type discussed in [26] with screening length l D ≈ 10 (dashed line). The equilibrium unequal-time correlator ρ − k (0)ρ k (t) of a Fourier mode ρ k = d 2 x exp(i k · x)ρ( x) of the charge density. This corresponds to the time evolution of a linear perturbation from the equilibrium state. The curve, which corresponds to wavelength λ = 128, is well fitted by an exponentially decaying sine wave, showing that interaction with microscopic degrees of freedom gives rise to effective dissipation for long-wavelength modes.
Thus, we have shown that, in spite of the apparent differences, our models has the same effective dynamics as the TDGL model. The agreement is only qualitative, however, so our simulations cannot make quantitative predictions about real superconductors. One should, therefore, carry out a similar study using the TDGL approach and realistic parameter values. In this paper, in contrast, our aim is only to test the general validity of the scenario, and a simple toy model is more suitable for the purpose.
In the TDGL-Langevin approach one can decrease the temperature directly by decreasing the amplitude the noise term. In our approach, as in the real world, this is not possible because temperature is not an external parameter, but rather a quantity that characterizes the state of the system. Therefore, we modified the equations of motion by adding small damping terms (7), (8), which leave the dynamics unchanged at microscopic time scales, but slowly cool the system down over long periods of time.