{'fraction_non_alphanumeric': 0.03462569832402235, 'fraction_numerical': 0.009474860335195531, 'mean_word_length': 4.225115301535408, 'pattern_counts': {'":': 0, '<': 16, '<?xml version=': 0, '>': 20, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 15, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Topological Defects from First Order Gauge Theory Phase Transitions

We investigate the mechanism by which topological defects form in first order phase transitions with a charged order parameter. We show how thick superconductor vortices and heavy cosmic strings form by trapping of magnetic flux. In an external magnetic field, intermediate objects such as strips and membranes of magnetic flux and chains of single winding defects are produced. At non-zero temperature, a variety of spontaneous defects of different winding numbers arise. In cosmology, our results mean that the magnetic flux thermal fluctuations get trapped in a primordial multi-tension string network. The mechanism may also apply to the production of cosmic-like strings in brane collisions. In a thin type-I superconductor film, flux strips are found to be meta-stable while thick vortices are stable up to some critical value of the winding number which increases with the thickness of the film. In addition, a non-dissipative Josephson-like current is obtained across the strips of quantized magnetic flux.


Introduction
The Kibble mechanism [1] has proved successful in explaining the formation of topological defects in phase transitions of field theories which experience spontaneous breakdown of a global symmetry. It applies to both first and second order phase transitions. The vacuum (minimum energy configuration) is degenerate and the system can fall into any vacuum state in different spatial locations. In first order phase transitions, where bubbles nucleate spontaneously, it has been tested theoretically in numerical simulations and analytical calculations [2,3,4]. Experimentally it has been verified in nematic liquid crystals through the spatial distribution of disclinations [5,6,7].
In second order phase transitions Kibble's argument was completed by Zurek [8] on what is known as Kibble-Zurek mechanism (KZ). It applies to the generation of topological defects as the system in question cools down in a finite time. When the transition takes place, non-equilibrium phenomena make the coherence length finite. This allows the system to choose different vacuum states in different locations. The mechanism has been tested experimentally in the prediction of the vorticity pattern of superfluid 3 He [9,10].
In cosmology, the analogues of disclination lines and superfluid vortices are global cosmic strings. They may have been generated in the early stages of the universe [11,12,13]. If they exist, Kibble mechanism is expected to explain their origin.
The order parameters which describe the phases of liquid crystals, 3 He and field theories with global cosmic strings are neutral. In the contrary, the order parameter which describes the symmetry of a gauge theory carries a charge and is coupled to a gauge field. KZ mechanism is based on the dynamics of a neutral order parameter, and therefore it is insufficient for describing generally either the formation of local cosmic strings in cosmological phase transitions of gauge theories or the origin of vortex lines in superconductors [14].
It is easy to see this in a simple condensed matter example. When a type-II superconductor is slowly cooled from the normal to the superconducting phase in an external magnetic field, an Abrikosov lattice consisting of vortices of the same sign is formed. This cannot be explained by Kibble-Zurek mechanism because it predicts the same number of vortices and anti-vortices. In the same way, Kibble mechanism is not sufficient to explain the formation of intermediate structures such as strips of magnetic flux in type-I superconductors under an external magnetic field [15,16].
For cosmology and particle physics prospects, gauge theories are more relevant than global theories. Grand Unified Theories (GUTs) are gauge theories which predict generically the spontaneous formation of local cosmic strings, monopoles and other topological defects. In particular, cosmic strings are tubes of quantized magnetic flux. The value of the tension of the strings increases with their flux. The phenomenological effects of these strings are mainly gravitational and depend quantitatively on their tension. Different cosmological observations generally attribute different upper bounds to the tension of cosmic strings [17]. These are average values which could actually correspond to a distribution of windings as a result of the evolution of a primordial multi-tension string network [18]. Kibble mechanism, on the contrary, gives rise to a unique tension because all strings have unit winding.
In a more theoretical context, recent findings suggest that cosmic strings may be a generic prediction of super-string theory [19] and could form in brane collision scenarios [20]. In some brane-inflation models like, for instance, the K 2 LM 2 T model [21], Dstrings are produced as a result of the collision of D-branes. The D-strings are intended to be cosmic strings and the D-branes contain a net magnetic flux. The resultant Dstring network confines that flux in an analogous manner that a vortex network does in a superconductor when an external magnetic field is applied. Therefore, a different mechanism to KZ must be responsible for the formation of those strings.
As it was shown in [14,22], generally in gauge theories KZ mechanism is only a good approximation in those cases where the dynamics of the gauge field is negligible. In the opposite limit, to which the examples above correspond, KZ mechanism is not enough and a different explanation is needed for the formation of topological defects. In second order phase transitions an explanation in terms of the non-equilibrium dynamics of the gauge field has already been found [14,23,24]. Numerical simulations show the formation of clusters of equal sign vortices.
On the other hand, the possibility for the production of high winding topological defects was already suggested in [25]. It was first argued in [26] and then discussed in more detail in [22] that they could form as a result of trapping of thermal fluctuations in first order phase transitions. We will investigate in this paper the mechanism by which thick vortices and heavy cosmic strings indeed form when the underlying symmetry is the U(1) gauge symmetry and the phase transition is strong first order. It is in this case the joint dynamics of the gauge field together with the complex phase of the order parameter -not necessarily out of equilibrium -that gives rise to trapping of magnetic flux.
The reader is warned that throughout this paper we commonly use the term vortex to refer to both cosmic strings and vortices regardless of the number of spatial dimensions we deal with. Only where the dimensionality is of some relevance we will stress the distinction. Also, we will refer to the regions in a minimum energy state either as 'vacuum regions' or as 'superconducting regions'. Likewise, we will refer to the regions in the symmetric phase either as 'symmetric phase regions' or as 'normal phase regions'. This paper is organized as follows. Section 2 describes the model and the setup of the numerical simulations. In section 3 we present the mechanism by which high-winding topological defects form in gauge theories. Section 4 is devoted to the mechanism of formation of strips of magnetic flux and chains of single-winding defects. In section 5 we explain the different properties of the vortex lines in a superconductor film with respect to the vortices in the fully two-dimensional theory. In section 6 we argue that spontaneous vortices and cosmic strings form in the thermal ensemble. Section 7 describes some mathematical aspects of the U(1) gauge theory which lead to the derivation of Josephson-like effects as well as to the reinterpretation of the so-called geodesic rule.

The model
Let us consider the Ginzburg-Landau (G-L) free energy for superconductors near the phase transition point: where m is the effective mass of a Cooper pair and µ 0 is the magnetic permeability of free space. Let us also consider the Abelian Higgs model, which can be seen as an idealized relativistic extension of the above G-L theory. In natural units, ‡ its Lagrangian is Greek indices run from 0 to the number of space dimensions. We will use the indices 0 and t to denote time components and both Latin indices and x, y, z for the spatial components. In the following tensorial formulae only the spatial components apply to (1). The charged complex scalar field φ ≡ |φ| exp (iθ) is generically referred to as 'the order parameter'. It is interpreted as the quantum wave function of Cooper pairs in the G-L theory and as the Higgs field in a particle physics context. The vector field A µ is a U(1) gauge field intended either as the electromagnetic field of the usual Maxwell theory in the context of superconductivity or as a U(1) field appearing in the spectrum of a grand unified theory (not necessarily the photon field). F µν is the Maxwell field strength and D µ is the covariant derivative operator: In the G-L theory the value of e is twice the electron charge in agreement with the interpretation of a Cooper pair as a two-electron bound state. The Lagrangian (2) and the free energy (1) are invariant under the U(1) gauge transformations θ → θ + Λ, where Λ is a real scalar function of space and time.
In the context of Grand Unified Theories, the model described by the Lagrangian (2) is meant to be just an idealized toy model. When implementing radiative corrections, additional couplings of A µ and φ to other fields must be considered. Also first order time derivatives in the effective equations of motion for both fields, φ and A µ , emerge as considering non-zero momenta corrections of order one.
Concerning dynamical phenomena in superconductivity, the time-dependent Ginzburg-Landau equations (TDGL) were derived in [27,28,29] as effective macroscopic equations from the microscopic BCS theory. A first order time derivative of φ emerges as a result of the coupling of the electrons in the superconducting phase to the electrons in the normal phase. The equation for A i also contains a first order time derivative as a result of the finite conductivity in the normal phase. This is an ohmic current. Thus, the dynamics of both fields, φ and A i , is approximately diffusive in a superconducting film. The TDGL equation for φ does also include a 'small' second order time derivative. The equation for A i must contain a second order time derivative as well. This is necessary if the gauge field has to satisfy Maxwell equations out of a superconducting film [24]. Therefore, we can explain with analogous equations of motion the dynamics of both a ‡ Natural units, c = = µ 0 = ǫ 0 = k B = 1, are used hereafter.
superconductor -described by the time-dependent G-L theory -and a particle physics phase transition -described by an effective field theory. In the temporal gauge, A 0 = 0, those equations read where σ is the damping rate, ǫ ijk is the three dimensional Levi-Civita tensor and summation over spatial indices applies if repeated in different location on pairs of tensor symbols. The electric current Following the original phenomenological Ginzburg-Landau theory in which the free energy can be expanded analytically in powers of |φ| near the phase transition and the fact that many GUTs predict a radiatively corrected effective potential of a similar kind, we will work with a scalar potential of the form where m H , κ and λ are in general positive definite temperature-dependent parameters. A first order phase transition takes place at some critical temperature T c . Below T c the coupling parameters must satisfy the inequalities 9κ 2 − 32m 2 H λ > 0, κ 2 − 4m 2 H λ > 0 for the existence of a meta-stable minimum at φ = 0. The true minimum is located at |φ| min = (3κ + 9κ 2 − 32λm 2 H )/8λ. In what follows we will refer to |φ| min as the vacuum expectation value (v.e.v.) and to the set of absolute minimum energy local states as vacuum manifold. Those states are characterized by the conditions |φ| = |φ| min , The transition takes place by means of spontaneous nucleation of bubbles of true vacuum. These bubbles appear as a result of the growth of very special and infrequent fluctuations of φ [30,31]. Spontaneous symmetry breaking (SSB) occurs as the order parameter φ chooses randomly a point α on the U(1) vacuum manifold such that φ vac = |φ| min e iα ( figure 1(b)). That way, spatial inhomogeneity in the manner the symmetry breaks down arises naturally in a first order phase transition. This is a requirement for defect formation.
Nucleation of bubbles and expulsion of magnetic field out of them drive the system out of equilibrium. These processes produce latent heat. Restoration of equilibrium involves readjustment of the magnetic flux and steady bubble expansion. It is achieved by means of the dissipation provided by the damping terms which, at the same time, remove the latent heat.
The damping rates -conductivity in the case of the gauge field -determine the time scale of the dynamics of the fields -the steady bubble expansion rate in the case of the scalar field. In principle no a priori assumption exists as regarding the damping rates. However they are important phenomenologically since their values determine whether the fields evolve in adiabatic, diffusive or dissipative regimes. For the sake of simplicity, we will assume that conductivity is uniform in the symmetric (normal) phase and we will assign the same damping rate σ to both fields. In doing so we guarantee conservation of charge by fulfilling Gauss' law ∇ · E = −q [24]. § 2.2. The gauge choice and the gauge-invariant field D µ θ Beyond tree-level, spontaneous breakdown of a local symmetry has been proved to be impossible without gauge fixing [32]. In the mean-field theory approximation a similar situation holds. Spontaneous symmetry breaking manifests itself in the gauge fixing as we will illustrate in section 4.
The fact that topological defects do form implies that a magnetic field gets confined around a discrete number of zeros of φ. It results in a multiply-connected vacuum region. Along any curve passing through one of the zeros of φ the complex phase θ presents a discontinuity of value π. These discontinuities cannot be gauged out because they are physical and so gauge-invariant.
Deep in a simply connected domain of broken phase like a bubble, it is always possible to choose the unitary gauge or, more generally, θ(x) = α uniform. In fact, because the magnetic flux is expelled by Meissner effect, it is natural to choose A µ = 0 there. Conversely, in the symmetric phase region between bubbles where φ = 0, θ is ill-defined and so is ∂ µ θ. As the bubbles of broken phase expand and collide they may form a multiply-connected broken phase region. The gauge θ(x) = α can only be maintained if this region is not to contain topological defects. If topological defects do form, the unitary gauge is not possible because, for any closed loop C around a defect, The above situations illustrate the fact that neither ∂ µ θ nor A µ alone are good quantities to be tracked in order to explain the formation of topological defects in gauge theories for they are not gauge-invariant quantities. In fact, the relevant degrees of freedom to be considered in the region where topological defects are to be formed are in the gauge-invariant combination of the Goldstone modes and the gauge field

Setup of the numerical simulations
In the numerical simulations presented in this paper we deal mainly with two kinds of scenarios. In the first one we show how thick vortices can be obtained when an external magnetic flux gets trapped in a multi-bubble collision. At the same time, chains of § Violation of conservation of charge is not actually a problem in superconductors phenomenology. In fact the TDGL equations violate the gauge symmetry explicitly if one considers the system they describe as isolated. See section 7.1 for further comments on this definition.
single-winding vortices form at the junction of each pair of bubbles. In the second scenario we show how chains of vortices and flux strips form in the gap between two parallel superconducting blocks in the presence of a net magnetic flux.

Multi-bubble collision picture
Nucleation of bubbles is implemented by setting-up disconnected spherical bubbles of broken phase. Their initial profile |φ|(r) is approximately that of the bounce solution in the thin bubble wall approximation [30,31]. Spontaneous symmetry breaking consists of choosing randomly a uniform complex phase value in each bubble: where r is the space-coordinate vector, m 2 h is given by (58) in section 5, N b is the number of bubbles, α i is the initially uniform complex phase of bubble i, r i stands for the position vector of the center of bubble i and R N is the initial radius of the bubbles. R N is large enough so that bubbles expand. When letting the system evolve with the equations of motion, the bubbles grow until colliding. The distance between bubbles is large enough so that the bubble-wall profile takes its actual form much before they collide. The bubble expansion rate is controlled by the damping rate σ appearing in (6). Soon after bubbles begin to grow the expansion rate v approaches v ≈ exp (−σ).
When a uniform magnetic flux is present, nucleation of bubbles requires the expulsion of the magnetic flux out of them before they start expanding. We first set up a uniform magnetic field B parallel to the y-axis in the gauge A To expel flux we let the gauge field evolve according to (7) while the scalar field configuration is kept fixed at its original form (11). That way, it is by means of dissipation that the system relaxes until getting to the minimum energy configuration compatible with both the nucleated bubbles with profile (11) and the presence of a given magnetic flux in equilibrium. The optimal damping rate is given by σ relax ≈ 1/l, where l is the typical length of the lattice. The time it takes for the magnetic field to relax is approximately t relax ≈ l. After relaxation, both φ and A i are evolved according to equations (6) and (7). As it was mentioned above, in order to satisfy conservation of charge, the same damping rate σ is used in both equations.
The approach presented here corresponds to the mean-field zero-temperature effective theory. However, when a uniform external magnetic field is applied, this field resembles a long wave-length magnetic field thermal fluctuation.

Collision of two-superconducting blocks with magnetic flux in the junction
This set-up is meant to mimic the collision between two big bubbles when a magnetic flux gets trapped in the junction. It is a good approximation as long as the radius of the bubbles is much larger than the width of the bubble-wall. The cubic lattice has dimensions l x × l y × l z , where 0 ≤ x ≤ l x , 0 ≤ y ≤ l y and 0 ≤ z ≤ l z . The collision axis is the x-axis, orthogonal to the walls of the blocks. The gap between the blocks has width d and its middle point is localized at x = l x /2. A uniform magnetic field parallel to the y-axis is confined in the gap. The total flux is quantized with value 2π e N w . Periodic boundary conditions are imposed in the directions z and y and the axial gauge A x = 0 is chosen: where B = 2πNw e d lz and (α 2 − α 1 ) is the spontaneous gauge-invariant phase difference (see below) between both blocks.
In the fully two dimensional simulations, l y = 0. In the superconductor film model, the film thickness µ satisfies µ ≪ l x , l y , l z . In both models, we refer to the strip between blocks of width d along the x-axis, length l z along the z-axis and thickness µ -which is zero in the fully two dimensional model-as 'junction'. We will show that it is in many aspects similar to a Josephson junction. In the fully three dimensional model, as the blocks extend initially all along the y-axis, translational invariance holds at any time in the y direction and the problem is actually two dimensional. We refer to the x direction as 'transverse direction across the junction' and to the z direction as 'direction along the junction'. In the collision of two three-dimensional bubbles along the x axis, rotational invariance in the y, z plane rather than y-invariance holds if no magnetic field is applied. We are however interested in the situation where a magnetic field is applied in the y direction. Since we take a planar approximation in which the bubbles are big enough so that they can be considered as infinite blocks in the z, y plane, the problem becomes effectively two dimensional. Therefore, we will maintain the nomenclature we use in the fully two dimensional model as referring to the directions along (i.e. the z direction) and transverse to (i.e. the x direction) the junction.

The superconductor film model (SC)
The superconductor film simulations require extra conditions to account for the finite film thickness effects. In the first place, spherical bubbles are replaced with cylinders of length equal to the thickness of the film µ. More importnatly, vanishing of the variation of the free energy (1), δF , with respect to arbitrary variations of the scalar field φ * in the superconductor film not only implies (6). An additional surface term with the form Ω δφ * (r)D i φ(r) ds i (r) (12) must vanish as well. Ω denotes the boundary of the film and ds i is the differential vector normal to the surface at point r. For arbitrary δφ * , in order for (12) to vanish, it is necessary that the component of D i φ normal to the film goes to zero at any point of the boundary. On the other hand, because the electromagnetic field does not couple to any charge particle out of the film, purely Maxwell equations hold in the free space and no dissipative term (i.e. σ = 0) is to be included there. That is, if the film is placed in the xz-plane with center at y = 0 and thickness µ, the additional superconductor film conditions read: We will be mainly interested in the stability of flux strips and thick vortices. The setup for the formation of flux strips is the same as in the previous subsection for the confinement of magnetic flux in between two superconducting blocks.
Formation of thick vortices of winding N w > 1 is achieved by confining initially a flux of value 2π e N w in a square of side d equal to the typical diameter of the vortex, that is d ≈ √ 2Nw e|φ| min . If the center of the square is placed at (x = l x /2, y = 0, z = l z /2), the initial conditions are: where B = 2πNw e d 2 . As the system evolves, a stable, meta-stable or unstable thick vortex form. The lattice size is big enough so that finite size effects are negligible. Periodic boundary conditions are imposed on the direction y orthogonal to the film. Some simulations in this paper differ slightly from the ones presented here. We will describe the corresponding setups in each particular case. A description of the lattice discretization method is found in Appendix B. The values of the parameters used in each simulation are compiled in table C1.

Global Theory
First suggested by Kibble in [1], the formation of topological defects in a scalar field theory with a global continuous symmetry can be explained by means of the spontaneous breakdown of the symmetry and the so-called geodesic rule.
In a first order phase transition of the Abelian global theory, spontaneous symmetry breaking takes place as bubbles nucleate at arbitrary points of the U(1) vacuum manifold. Since in the global theory any vacuum configuration requires ∂ µ θ = 0, the value of the complex phase of the scalar field in each bubble remains uniform as the bubble expands. When two bubbles coalesce the above condition cannot be satisfied at the contact area unless both bubbles nucleated with the same complex phase value. The geodesic rule states that the complex phase θ of the scalar field must interpolate across the junctions between pairs of bubbles in such a way that it traces the shortest path in the U(1) manifold. When several bubbles nucleated at different points {α i } of the vacuum manifold meet (figure 1(a)), the complex phase interpolates between each (α i , α j )-pair as shown in figure 1(b). In the case of a three-bubble collision the condition for the formation of a single winding topological defect is This condition is fulfilled with probability 1/4. As the complex phase interpolates according to the geodesic rule, the gradient energy ∂ i φ∂ i φ * dx across the junctions between bubbles, a < x < b, is minimized: where ∂ i |φ|∂ i |φ| is the bubble-wall tension and |φ| 2 ∂ i θ∂ i θ is the current self-interaction energy density. ¶ More precisely, what is minimized by following the geodesic rule is the term associated to ∂ i θ, that is, the current self-interaction energy density. ¶ The one-dimensional integrations in (15) are actually either linear or superficial densities of energy in two or three spatial dimensions respectively. For brevity, we will refer to them as energy.

Gauge Theory
In the gauge theory ∂ µ θ is not a gauge-invariant quantity anymore. The analogue of equation (15) is where the current self-interaction energy density reads ε(x) ≡ |φ| 2 D i θD i θ. When two bubbles collide, a different rule to the geodesic follows: The phase gradient D x θ and the gauge-invariant phase difference γ ab ≡ b a D x θ dx across the junction a < x < b between any pair of bubbles decrease in time until vanishing. In doing so, A i evolves (and so does ∂ i θ) satisfying the conservation of flux requirement. That way, the magnetic flux is expelled, a gradient in the complex phase, ∂ x θ, is generated and the current self-interaction energy gets to its minimum value. At the same time the order parameter takes the v.e.v. across the junction.
In a multiply-connected broken phase region the above rule applies to any given pair of vacuum sectors along each loop C enclosing a simply-connected symmetric phase area. Thus, C D i θ dx i = 0 holds when C is contained in the surrounding vacuum region.
In absence of magnetic field the geodesic rule is a good approximation. In fact it is equivalent to the gauge rule above if the gauge A 0 = 0, A x (t = 0) = 0 is taken, where x stands for the direction along the path traversing from one to another bubble. Hence the formation of a single vortex between the three bubbles of figure 2 can be explained by Kibble mechanism. When the three bubbles of figure   in absence of an external magnetic field can be found in section 7.3.
The difference between the gauge rule and the geodesic rule in the gauge theory manifests itself in the presence of an external magnetic field. Let us consider for instance the phase transition of an infinite superconductor film in an external magnetic field. Because the total flux is conserved, the resultant vortex pattern must consist of a lattice of vortices of the same sign. On the contrary, the randomness of the complex phase pattern obtained by SSB would give rise to the same number of positive and negative vortices according to Kibble's argument. Hereafter we will refer to the above rule in gauge theories as geodesic-Meissner rule.
A typical scenario for the formation of a high winding topological defect is a multibubble collision in presence of a magnetic flux. As the bubbles of figure 4 grow, D i θ goes to zero inside them and the magnetic flux is expelled. That is the Meissner effect. Ultimately, when the bubbles meet and the geodesic-Meissner rule is applied to each two-bubble junction, the flux in the symmetric phase region enclosed by the bubbles gets trapped and a thick vortex forms at the center of the collision. The main contribution to the total flux in the vortex comes from the trapping of part of the magnetic flux which is at the moment bubbles nucleate, t N , inside the triangle with contour curve C whose vertices are the centers of the bubbles ( figure 3(a)). An upper bound is in fact is the contour curve of the area enclosed by the bubbles when they meet at some time t > t N (figure 3(a)). The precise amount depends on kinematic effects, i.e. on the bubble expansion rate and the dynamics of the magnetic field. The upper bound Φ max F.T. is saturated when either the bubbles expand at a rate close to the speed of light or the dynamics of the magnetic field is diffusive. In both cases the magnetic flux expelled from the bubbles gets stuck to the bubble walls (within a layer of thickness ∼ 1/σ in the latter case) all the way until they meet. This effect appears signaled by the rings in bright green color surrounding the bubbles of figure 5. In the opposite limit, when the bubbles expand adiabatically and the magnetic field escapes from bubbles at the speed of light, Φ min F.T. is a better approximation. In general, the precise value of Φ F.T. depends on the bubble expansion rate and the damping rate of the magnetic field. In addition, order-one numerical factors are determined by the actual geometry of the collision. The quantitative study will be reported in a separate article [34].
The snapshots of the two-dimensional simulation in figure 4 correspond to a soft bubble collision in an external magnetic field. After the three bubbles collide gently, i.e. v ≪ 1, a thick vortex gets formed at the very center.
because σ is chosen so that the dynamics of |φ| is close to adiabatic while the dynamics of the magnetic field is diffusive. It is easy to calculate its value by simply counting the times the complex phase increases 2π units around the vortex.
The contribution of Kibble mechanism through SSB is negligible in the presence of an intense magnetic field, that is, for Φ F.T. ≫ 2π/e. Let us consider a fixed value for the initial magnetic field in the setup of figure 4(a). If one chooses randomly the initial {θ i } values of the complex phases of the three bubbles the winding number of the vortex that gets formed at the center of the collision takes just one of the three possible values e 2π Φ F.T. ± 1. Only if the contributions due to Kibble mechanism and Flux Trapping were additive, according to the probability for condition (14) to be satisfied, the probability P of getting each value would be: However we will show at the end of section 4 that Kibble mechanism and Flux Trapping are not actually additive. If the collision is violent, i.e. v ∼ 1, as it is the case in figure 5, the flux trapped corresponds roughly to Φ max F.T. . The non-equilibrium dynamics is relevant, reason why the flux trapped is not only confined in the form of a thick vortex at the very center of the collision. On the contrary, chains of single-winding vortices form at the junctions.

Magnetic flux strips and single vortex chains
By causality, when two bubbles collide in an external magnetic field, a finite amount of magnetic flux gets trapped at the junction. Let us assume in first approximation that small structure effects -i.e. effects at scales comparable to the bubble wall thickness ∼ m −1 h -can be neglected in the dynamics of the bubbles. That way the bubbles expand isotropically at constant rate v at all times, even after their walls overlap. Let us assume the less favorable situation for magnetic flux to get trapped. That is, let us assume that the magnetic flux consists of massless magnetic particles which are expelled from the surface of the bubbles at the speed of light as if they were light-rays. That way we are neglecting several effects. First, we are assuming no interaction between the magnetic field and the scalar field at small scales of the order of ∼ m −1 h . Second, we are ignoring the long range magnetic field interaction, which we will show later on is present in the superconducting model. Third, we are neglecting any diffraction effect. All these effects contribute positively to the trapping of magnetic flux otherwise. As a consequence, the estimate of the flux trapped that we proceed to calculate disregarding these effects is a lower bound.
Let us take two identical bubbles separated a distance D, nucleated at time t = 0 Figure 6. Trapping of a 'magnetic ray' by causality. A 'magnetic particle' is expelled at the speed of light c from the surface of the left bubble at (x 1 (t 0 ), z 1 (t 0 )). It gets trapped by causality at some point (x = D/2, z) along the junction between the two bubbles. The centers of the bubbles are separated a distance D and both bubbles expand at constant rate v. and expanding at constant rate v as shown in figure 6. A magnetic particle abandons the surface of the bubble on the left at height z 1 (t 0 ) from the collision axis when the bubble radius is r 0 = vt 0 . The ray experiences several reflections as it strikes several times on each bubble. These reflections may, or may not, follow Snell's law. In the case they do, the process is referred to in [34] as magnetic optic approach. Nevertheless, in the worst of the cases the exit angle is arbitrarily large at each reflection. However, the distance the magnetic particle travels in a time interval ∆t = t − t 0 > 0 cannot exceed c(t − t 0 ). If t > D/2v, bubbles meet and overlap during that time interval such that the vertex of the overlap is located at height z vertex 2 (t). Because the speed of the vertex along the vertical axis z is initially superluminal, it is possible that for some time interval T f , where (x 1 (t 0 ), z 1 (t 0 )) stands for the coordinates of the point from which the ray is first expelled and (D/2, z vertex 2 (t f )) stands for the coordinates of the vertex. If inequality (18) is satisfied for some t f ∈ T then the light-ray must have got trapped at some point of the junction (x = 0) where bubbles overlap.
The vertex height, z vertex 2 , evolves in time as The maximum height z 2 at which a photon expelled at (x 1 (t 0 ), z 1 (t 0 )) can get along the junction at time t > t 0 is From (20) one reads that the rate at which z 2 grows is initially infinite. However, from the analysis of its second derivative, it tends asymptotically to c faster than the first derivative of (19) tends to v. As c ≥ v, inequality (18) may get saturated at time t M ax f for some pair ( Therefore, the flux in the differential area z M ax 1 (x 1 )dx gets trapped. We have solved the equations above numerically discretizing the x-axis with space step δx = D/2N. We search for z M ax Starting with a uniform initial magnetic field B 0 , the flux trapped by causality is given by Figure 7 shows the causality bound estimate (curve (b)) in comparison to the actual data (curve (a)) of the simulation of a two bubble collision in an initially uniform magnetic field B 0 . Quantization of the magnetic flux has been taken into account in the causality bound computation. A damping rate σ is used to slow down the bubble expansion rate v through (6). However, no damping term is included in the equation for A i . That way the dynamics of A i is neither dissipative nor diffusive. Having shown that magnetic flux gets necessarily trapped at the junction by kinematical reasons (modulo quantization effects), we aim to explain how it gets confined in the form of vortex chains like those in figure 5(d). This involves the study of small structure effects, that is, the dynamics of |φ| and D i θ at the junctions.
Let us simplify the problem by considering that the collision area of two big bubbles can be approximated by two parallel superconducting blocks separated by a gap of width d in the transverse direction x in which a net magnetic flux is placed. For the sake of simplicity we will work with the fully two-dimensional theory. Periodic boundary conditions (pbc) are imposed in the direction z along the junction. The total flux in the gap is quantized. Its value is 2π e N w , where N w is an integer. The reason for imposing pbc is to confine the flux as it would happen when two big bubbles coalesce (by dynamical reasons in the latter case).
Given two disconnected broken phase blocks we introduce a uniform magnetic flux in between. With no loss of generality we will choose the gauge A x = 0 at t = 0 as shown in figure 8. In this gauge the field A i is continuous and so are ∂ i θ and D i θ in the  areas of broken phase. Furthermore, the gauge-invariant phase difference γ ab (z) across the junction when the blocks meet at time t = 0 is equivalent to the complex phase difference ∆θ ab (z). That is Therefore, according to the gauge choice in figure 8, where (α B − α A ) is the spontaneous phase difference and −(N w /R) z is the induced phase difference. + After fixing the gauge there exists still the freedom in choosing ∆α ab ≡ (α B − α A ). This freedom reflects the invariance of the system along the junction (z-invariance) when the blocks are separated. The equality γ ab = ∆α ab holds in absence of magnetic field. Thus, the arbitrary choice of ∆α ab corresponds to the spontaneous breakdown of the U(1) symmetry as announced in section 2.2. The spontaneous difference ∆α ab is physical and manifests itself as blocks overlap and the z-invariance gets broken. At the moment the blocks get in contact, the physical local gauge-invariant quantities are not translation invariant anymore along the junction due to the presence of the magnetic field (see Appendix A for a proof). The scalar field profile falls off exponentially within a typical length equal to the inverse of its mass 1/m h . Consequently, it is necessary that the size of the gap between the blocks d is less than d s ≈ 2/m h so that the blocks overlap and interact with each other.
If d ≫ d s , the blocks behave as fronts of big bubbles whose expansion is driven by the potential V (|φ|) and the bubble-wall tension ∂ i |φ| ∂ i |φ|. The z-invariance of |φ| and B along the junction persists and no current traverses the gap. Consequently no vortex can form. On the other hand the trapped flux gives rise to a repulsive force which, if d > d s , can balance the expansion of the blocks towards each other. If one neglects the bubble-wall tension in first approximation it is easy to calculate the critical value for the magnetic field necessary to prevent the blocks from approaching each other beyond some distance greater than d s . At zero-temperature the equilibrium configuration corresponds to the minimum energy configuration compatible with a net magnetic flux in the junction. In the fully two-dimensional theory the magnetic field B does not present long-range interaction in the symmetric phase and the equilibrium configuration corresponds to a magnetic flux density uniformly distributed in the junction. Therefore, for a fixed value of the magnetic flux, the magnetic energy depends only on the size of the junction. The equilibrium distance d eq between the two superconducting blocks of figure 8 is the distance at which the expansion of the blocks is compensated by the + See section 7.1 for further explanation about this decomposition. magnetic pressure, where 2πR is the length of the blocks along the junction, Φ E.M. is the total flux and V min stands for |V (|φ| min )|. Requiring d eq = d s we calculate the minimum value for the magnetic flux so that the blocks never overlap and the flux strip remains straight, that is, If the value of Φ E.M. in the gap is less than Φ strip,min the small scale structure effects are relevant and the dynamics of both |φ| and D i θ cooperate towards the formation of topological defects. A detailed analysis of the generation of Josephson-like currents can be found in [33]. In what follows we just make use of the essential results there.
At the moment the blocks get in contact, γ ab (z) = This implies that D x θ(z) is a function of γ ab (z) and the equations of motion incorporate non-local effects and non-linear terms in D x θ. In particular, a Josephson-like coupling shows up (see also sections 7.2 and 7.3). The quantity to be studied is γ ab (z) which, right at the moment the blocks get in contact, follows a perturbed sine-Gordon equation in 1 + 1 dimensions [35], where |φ| 0 is the z-independent factor of the norm of the scalar field at the time the blocks meet, i.e. |φ| 2 |x|<χ/2 ≃ |φ| 2 0 × f (z). In [33] |φ| 2 0 is referred to as 'reduced Cooper pairs density'. For χ ≪ d, it is Therefore, initially, the vortex formation dynamics corresponds to that of onedimensional sine-Gordon solitons. At this stage the current self-interaction energy density depends on z as At |φ| 0 fixed, equation (24) implies that ε ab gets maximum at those points z n satisfying γ ab (z n ) = (2n + 1)π, which are where the linear flux density is expressed in terms of the winding number, dΦ E.M. d z = Nw eR . When the blocks get in contact, φ| x<0 and φ| x>0 interfere at the middle of the gap x = 0 giving rise to an interference pattern along the junction as Such an interference is constructive at the points z n ′ such that γ ab (z n ′ ) = 2πn ′ (equal colours according to the color gradient legend in figure 8). On the contrary it is totally destructive if (30) is satisfied. As a result N w minima of |φ| turn up at the points {(x = 0, z = z n )} when the blocks first meet. (see figure A1(b)).
In [33] it is shown that |φ| 0 is not z-invariant anymore either as |φ| follows the equation of motion (88). We will just use here a minimum energy principle to argue for the formation of vortices. In view of formula (31), the potential (29) (27) must be intended as a good starting point to evaluate the initial evolution of γ ab . Afterwards, additional non-linear terms enter and the problem becomes two-dimensional in space.
It might still happen that the magnetic repulsion stops the vortex formation process. Unit-winding vortices are rotationally symmetric with a radius r v that scales with the London penetration length λ L which is equivalent to the inverse of the photon mass [36]: Therefore the inequality πR N w r v must be satisfied in order to accommodate N w single vortices in a chain of length 2πR. This condition imposes an additional constrain on the maximum value of the flux trapped for a chain of single-winding vortices to form: In the intermediate situation the magnetic flux lies in the range Φ strip,min . In that case the N w minima of |φ| do form but single vortices do not. Instead, periodic flux 'ripples' appear along the junction between blocks. Figure 9 shows the three typical cases. The stability/meta-stability of the configurations 9(b) and 9(c) depends on the dimensionality of the theory. That is, flux 'ripples' and chains of singlewinding vortices are not stable against condensation into thicker vortices or heavier cosmic strings. However they can be meta-stable for high enough flux density. That meta-stability is stronger in the superconducting model where a long-range magnetic interaction is present (see next section). Figure 10(a) shows a membrane of magnetic flux in between two superconducting blocks. The flux membrane contains nine winding units of flux. They manifest in the form of 'ripples' in figure 10(b). The equilibrium of such configuration is just maintained by the periodicity in the direction z along the junction. However, when we introduce a small perturbation δ|φ| ≪ |φ| min at z = 0, the periodicity is broken and the structure becomes unstable. The vortices gather in groups of winding three (figure 10(c)) and three thick vortices of winding three get formed ( figure 10(d)).
The spontaneous phase difference α B − α A is indeed physical. It gives rise to a shift in the location of the vortices along the junction according to (24). In a two-  bubble collision, whether a vortex forms at either side with respect to the collision axis depends on this shift. In figure 11 two bubbles 'sandwich' a strip of magnetic flux as they coalesce. The relative shift difference of value 3π/4 between the cases (a) − (b) and (c) − (d) gives rise to the production of one more vortex along the junction in figure  11(b) with respect to figure 11(d). This is the main effect that SSB has in the total There is also another way in which SSB and flux trapping mechanism interact. Let us consider the current traversing the gap between the two bubbles in the presence of a magnetic flux. The total current follows an interference pattern (see formula (99), section 7.2) and oscillates as a function of the total flux. Its amplitude decreases notably as the flux increases and it is zero if the total flux in the junction takes integral values of 2π/e. As a result of that current a fluxon ring would develop around the collision axis giving rise to an extra contribution to the magnetic flux in the thick vortex (similar to what happens in figure 2 in absence of magnetic field). However, because that current is strongly suppressed for high values of the magnetic flux in the junction, we expect that fluxon to be very weak. Consequently its contribution to the flux in the thick vortex is negligible. For instance, in the collision of the three bubbles of figure 4, the probability for the final winding in the thick vortex to be e 2π Φ F.T. must be greater than the value 3/4 in (17) due to this interference effect. We conclude that he contributions from flux trapping and KZ mechanism are not additive.

The superconductor film model
The model which describes a superconducting film differs from the fully two-dimensional G-L model in the fact that the order parameter φ is confined in a thin film while the electromagnetic potential A i lives in the three dimensional bulk space. The magnetic field is not a pseudoscalar but a three-component pseudovector. It propagates in the bulk and induces a long-range interaction on the magnetic flux density in the film which is absent in the fully two-dimensional case [23].
Let us take a thin superconductor film of thickness µ, placed at y = 0 and parallel to the xz-plane. If a net magnetic flux is present, it gets confined in areas of symmetric phase. For a stationary configuration, the electromagnetic energy reduces to the magnetic energy confined in a volume V whose cross section at y = 0 is the symmetric phase area on the film, stat. is quadratic in B i , its most general form in terms of the film magnetic flux density B y ( r) reads where two-dimensional vectors with components x and z in the film are labeled with arrows on the top, V mag. ( r, r ′ ) is the induced magnetostatic potential interaction and the subindex sym. stands for 'symmetric phase area'. The stationary configuration corresponds to that for which (33) is minimum. That is, B y ( r) such that subject to the presence of a fixed magnetic flux Φ E.M. = sym. d 2 r B y ( r). We first calculate V mag. ( r, r ′ ) from the stationary Ampere's law, Because the current lives in the superconducting film on the plane y = 0, only the components x, z are non-zero in the vector equation (35). On the other hand, because the magnetic field is orthogonal to the film at y = 0 and the flux gets trapped in the symmetric phase regions, Ampere's law reads where Π y=0 sym. stands for 'symmetric phase area in the plane y = 0'. As a first approximation we will consider the film infinitely thin, that is, with zero thickness, µ = 0. That way, the y-component of the magnetic field on the film can be written in Fourier space as It is more convenient to write (36) in function of the gauge field A i (r), In the Coulomb gauge, We solve for A i (r) using the Green's function of the operator ∇ 2 , that is, G ∇ 2 ij such that In Fourier space,G The gauge field can be computed as and likewise the magnetic field, making use of (41) and expanding B y ( r) in Fourier modes, where ǫ ijk stands for the Levi-Civita tensor, the subscript y is fixed, m, n take labels {x, z} and i, j, k run over {x, y, z}. Next, we define the 'pseudo-Green's function' and its corresponding Fourier transform Therefore, we can write the magnetic energy as where † denotes adjoint. From (33), one reads immediately Again, we make use of the Fourier transforms to write, in the infinitely thin film limit µ → 0, After the integration of q, p, s (which yields a product of delta functions) we get to where k 4 is the 4 th power of the norm of the three-component wave vector k; r, b, m, d take labels {x, z} and i, j, a run over {x, y, z}. By applying the identity ǫ ijk ǫ klm = δ l i δ m j − δ m i δ l j and integrating the mode k y , it is easy to get to Comparing (51) with (47,48), we identify It is remarkable that V µ=0 mag. ( r 1 ′′ , r 2 ′′ ) equals the inverse of the two-point function B y ( r 1 ′′ )B y ( r 2 ′′ ) in the high temperature field theory, normalized to k B T = 1 [23]. The reason being that the long wavelength spectrum of fluctuations of the gauge field at high temperature -that is, the Rayleigh-Jeans spectrum of balckbody radiation-equals the Green's function of the stationary Ampere's equation (39) in the symmetric phase. Should we take into account the finiteness of the thickness of the film, µ, on the magnetic flux densityB y ( k), we would integrate those modes for which k y ≤ 1/µ. We may use a Gaussian function to suppress shorter wavelength modes. Replacing the delta function δ(r y ) in the integrand of (37) and thereafter with the Gaussian distribution 1 √ 2πµ exp (−r 2 y /2µ 2 ), we define the film magnetic flux density as This way, analog of (52) with µ > 0 reads where Erfc stands for 'complementary error function'. If the typical length of the symmetric phase area is much larger than the film thickness, the long wavelength limit | k| ≪ 1/µ can be taken in good approximation such that Erfc(µ| k|) ≈ 1 and the approximate expression (52) is recovered. We end up with the result that, in a superconductor film, vortices present a long-range interaction associated to their magnetic flux. Also, the Coulombian-like interaction V mag. ( r, r ′ ) implies that the magnetic flux density on the symmetric phase areas of the film behaves as if it were an electric charge density on the surface of an ordinary metal conductor.
When a magnetic flux is trapped in a symmetric phase region, the stationary configuration corresponds to that for which (33) is minimum subject to the existence of a total magnetic flux Φ E.M. . That is, B y ( r) such that δ δB y ( r) sym.

Stability and meta-stability of thick vortices
For the purpose of studying the dynamics of vortices, there are in principle two length scales to be considered, the coherence length ξ and the London penetration depth λ L . ξ is equal to the inverse of the mass of the order parameter at |φ| min , that is, ξ = 1/m h , where For the potential in (8), m 2 h = m 2 H − 3κ|φ| min + 6λ|φ| 2 min . Notice that even though m h is of the order of m H they differ by an order one numerical factor. λ L is the inverse of the photon mass, λ L = 1/m γ = ( √ 2e|φ| min ) −1 . The Ginzburg parameter, κ G ≡ λ L /ξ = m h /m γ , depends on both the value of |φ| min and the curvature of the effective potential m h at |φ| min . Finite temperature radiative corrections [37] and the presence of an external magnetic field [38] may give rise to variations in κ G . In a type-I superconductor, where λ L < ξ, there exists an attractive vortex-vortex force at distances shorter than ξ associated to the order parameter. At distances shorter than λ L there is also a magnetic force interaction which is repulsive for equal sign vortices. If no other force exists, two vortices of the same sign separated by a distance less than 1/m h can gather to give rise to a stable thick vortex. That is the case in the fully two and three-dimensional theories [39].
In a real type-I superconductor film, the induced long-range repulsive force (53) between either equal sign vortices or equal sign flux densities, modifies the picture above. In a thick film, µ ≫ ξ, the long-range Coulombian interaction is dominant at very long distances. Therefore, a meta-stable cluster of vortices of the same sign form. If the film thickness is of the order of ξ and still greater than λ L , it happens that high winding vortices are stable up to some critical value of the winding number. For higher winding they can be meta-stable . Finally, if the film is thin, µ λ L , the Coulombian repulsive interaction acts at distances r λ L and so only single-winding vortices are stable and an Abrikosov lattice form in an external magnetic field . The last case corresponds to the situation in which a thin type-I bulk superconductor film behaves as type-II in a perpendicular external magnetic field [40,41,42].
Let us concentrate on the case where both stable and meta-stable thick vortices can form. The winding number of the stable vortices cannot be arbitrarily high. Beyond some critical winding value, they break up into smaller pieces. The argument for this to happen is as follows. The total energy of a stationary winding N w vortex contains three terms, where C(µ/r v ) is a geometrical factor which depends on the exact form of V mag. ( r, r ′ ) through µ/r v -formula (53) is only exact in the thin film limit µ/r v ≈ 0 -, V = πV min and T stands for 2π times the tension of the vortex wall. Stability requires Considering the film thickness fixed and the vortex radius large in comparison to both µ and λ L so that formula (53) is a good approximation and the energy associated to the vortex wall tension is negligible (i.e. T r v ≪ CN 2 w r −1 v , µV r 2 v ), the first and fourth terms in (60) can be ignored and Consequently, and multiple winding vortices are not stable. Let us illustrate this phenomenon by analyzing the stationary configuration of the magnetic flux density in a thick vortex. For sufficiently large values of N w , the vortex radius scales as r v ∝ N 2/3 w and the magnetic energy is dominated by the long-range interaction of the flux density in (53). In that case the minimum-energy equilibrium configuration of the magnetic flux Φ E.M. = 2π e N w in a vortex approaches the distribution of an electric charge placed on an ordinary conducting surface of circular shape and radius r v [43]. Solving for B y ( r) in (57), where r is the distance to the center of the disc. Formula (63) gets maximum along the periphery of the circle rather than at the center as would be the case in a fullytwo dimensional vortex. As N w is increased, r v goes to infinity while B y | r=0 tends to zero. Such a tendency maximizes the interaction energy between the magnetic field and the scalar particles through (16). The magnetic flux density B y goes to zero at r = 0 where φ is zero while B y increases as r goes to r v where |φ| r=rv ≈ |φ| min , making the configuration energetically unfavourable. However, given a pair of values C(µ/r v ), µ in (60), the energy associated to the vortex wall tension may compensate the energy difference in the inequality (62) for small values of N w and r v . If that is the case, inequality (62) can be inverted. Figures 12 and  13 show this phenomenon in numerical simulations. In figure 12 the magnetic flux is 6π/e and λ L is approximately one in lattice spacing units. When the superconducting film is thick (µ = 21 in lattice spacing units) the magnetic flux profile presents a maximum at the very center and falls off exponentially towards the periphery of the vortex. As the thickness decreases, the profile of |φ| does not change sensibly. However the magnetic flux profile does approach the form of formula (63) and the maximum magnetic flux density migrates towards the periphery of the vortex. Finally, for a thickness µ = 1 ≈ λ L the winding-three vortex splits up in three single-winding pieces. In figure 12 only the profile of one of them is plotted. While the winding-three vortex is stable for µ = 21, it is meta-stable for µ = 9, 5, 3 and unstable for µ = 1. Figure  13 shows how the energy per winding number of thick vortices E(N w )/N w scales with the winding number N w in comparison to the energy of a single-winding vortex E(1). For µ = 1 only single-windings are stable. For 15 > µ > 1 only single-windings are still stable but numerical simulations show that winding two and three vortices are meta-stable. It is difficult to test how much meta-stable higher windings are because the number of possible combinations and geometries in which a thick vortex can decay increases quickly with its winding. For µ = 15 a winding-two vortex is stable against decay into single-winding pieces. For µ = 21, winding-four vortices are stable and, at least up to N w = 6, thick vortices do not split up in single-winding pieces. A winding-six is however meta-stable and can decay into a pair of stable winding-three vortices.
We expect that, by the mechanism explained in previous sections, a variety    of stable and meta-stable vortices form during the phase transition of a type-I superconductor film in an external magnetic field. Because of the inhomogeneous nature of a first order phase transition, the primordial lattice can be very different to the one corresponding to the minimum energy configuration. Its evolution in the first stages must be dominated by the long-range magnetic interaction. It is also very likely that, by the enormous amount of meta-stable intermediate states, the evolution of the vortex lattice gets stuck to some meta-stable configuration. It is only for very thin films that the only stable vortices are single-winding and the final configuration is a stable Abrikosov lattice. The Abrikosov lattice configuration is unique up to the fixing of its center of mass. The center of mass consists of two degrees of freedom which correspond to the two-dimensional analogue of the spontaneous phase difference ∆α ab that determines the position of the single-winding vortices in a finite vortex chain. Thus, they are also vestiges of the spontaneous breakdown of the gauge symmetry.

Meta-stability of flux strips
The strips of magnetic flux described in section 4 are generally unstable in the fully two and three-dimensional models. For thick superconducting films like the one in figure  10 where µ ≈ 3λ L they may be unstable as well. They can however be meta-stable for sufficiently thin superconducting films and strong enough magnetic field. That is so even if their width is much less than 2/m h and short distance effects become relevant. Numerical simulations show that, for a fixed value of the linear magnetic flux density higher than Φ chain,max E.M.
/2πR (2πR is the length of the strip), the resultant strip is metastable for a film thickness less than some critical value. Conversely, given a fixed value for the film thickness, there exists a critical magnetic flux linear density for which the strip is meta-stable. As an example let us take the same setup of figure 10 for a superconductor film of finite thickness µ. Using the parameter values in Appendix C, if µ is set to 3 in lattice spacing units, the resultant strip is unstable for N w = 9 and three thick vortices of winding three form as in figure 10. However, if we increase the total magnetic flux up to 12 2π e , the flux strip configuration turns out to be meta-stable. Likewise, fixing the value of the total magnetic flux at 12 2π e , instability takes place when increasing the thickness up to µ = 7. As a result, two thick vortices of winding six form. Figure 14 shows both configurations. Similar structures of intermediate states have been observed experimentally [15,16]. However the ones described in [15] and [16] do not arise as a consequence of a phase transition in the presence of an external magnetic field but as a result of the restoration of the symmetry in a type-I superconducting film as a uniform external auxiliary magnetic field H is increased.
In order to test the meta-stability/instability of a flux strip, we proceed as follows. We first let the system relax to some relative minimum energy state compatible with periodic boundary conditions along the z-axis. That gives rise to a non-straight strip like the one in figure 10(b). Both the electric field and the scalar field momentum are zero at that moment. Afterwards we implement small perturbations δ|φ| ≪ |φ| min randomly along the strip, once at a time. The size of the fluctuations in the transverse direction x vary from one lattice spacing up to the width of the strip. In the longitudinal direction z the size takes values from one to πR. If the perturbation does not grow the system goes back to the original strip structure. That a strip is meta-stable and not stable can be easily seen by comparing its energy with the one in which the same amount of flux is confined in separated stable vortices. The latter is energetically favourable.

Formation of spontaneous topological defects
In previous sections we have shown how local vortices and cosmic strings form in the collision of bubbles by trapping of an external magnetic flux. A similar mechanism applies to the formation of topological defects when, instead of an external magnetic field, magnetic flux is provided in the form of thermal fluctuations. One can think of the magnetic field long wavelength fluctuations as a uniform external magnetic field. The actual magnetic field fluctuations and nucleation of the actual bubbles need of a Monte Carlo simulation like the one performed in [18]. The analytical computation is a difficult task. In the finite-temperature field theory a first order phase transition effective potential of the form (8) comes from radiative corrections to the zero-temperature potential. In the particular case of the U(1)-invariant renormalizable theory, The first order transition effective potential has been calculated in the type-I regime (λ/e 2 < 1) in several contexts. Some of them are concerned with quantum corrections [44] while others deal with thermal perturbations [45,46,47,13,48]. Because e 2 > λ, radiative corrections due to the gauge field are dominant. The nucleation rate and the bubble structure are affected by the interaction of the fluctuation of both the scalar field and the gauge field [49,50]. In particular, this implies that the bounce solution and the nucleation rate as calculated in [30,31] are only approximations. We will therefore simplify the problem by considering that the nucleation rate Γ, the bubble expansion rate v and the typical bubble nucleation radius R N are free parameters which must be ultimately given by the actual field theory. Thus, the finite-temperature vacuum manifold consists of the set of minima of a potential like −λ(φφ * − η 2 ) 2 radiatively corrected. The equilibrium state corresponds to the Gibbs ensemble. Thermal fluctuations of the gauge field are strongly suppressed by Meissner effect in the broken phase where the photon is heavy. They are however present in the large areas of symmetric phase between bubbles. Therefore, with regards to the production of spontaneous magnetic field, we will concentrate on the long-wavelength fluctuations of A i .

Estimation of the typical winding number of thermal topological defects
We will review here the arguments and predictions presented in [22] and [18].
Let us assume that the dynamics of the whole system is adiabatic. Therefore, bubbles expand slowly enough such that the magnetic field can always be considered in thermal equilibrium in the symmetric phase areas between bubbles. Thus, the spectrum of the wavelength fluctuations corresponds to the long wavelength spectrum obtained in the thermal ensemble where the saddle point approximation is accurate enough given that the fluctuations are nearly gaussian. The spectrum is so the two-point correlation function of the gauge field in the symmetric phase. We have already calculated the functional form of that spectrum in the context of the zero-temperature field theory in section 5. In the thermal ensemble it is the Boltzmann factor exp (−F/T ) what determines the weight of a configuration. Therefore, the spectrum of fluctuations is the Green's function of the equation of A i for an extreme of the function F/T , that is: ∇ 2 n G ij = T δ n (r − r ′ )δ ij . * The label n denotes the number of space dimensions and so the functional form of the magnetic field fluctuations depends on the spatial dimensionality: Formulas (65), (66) and (67) correspond to the fully two-dimensional, superconducting and fully three-dimensional models respectively. The superconducting one (66) was derived in [23] following similar steps to those we followed in the computation of (53) in the context of the zero-temperature theory. Again, top-arrows label two-dimensional vectors living in a superconductor film which is embedded in a three dimensional space and B y stands for the component of the magnetic field orthogonal to the film. Generally, we will write the formulas above as where the label m denotes the model to which the function corresponds and the indices i, j and the integer n take their values accordingly. As bubbles expand the magnetic flux associated to those fluctuations is expelled by Meissner effect. When several bubbles meet, a simply connected region of symmetric phase arises at the center of the collision. The magnetic flux in that region gets trapped by the mechanism explained in section 3.2 and a thick vortex or heavy cosmic string form. The typical amount of flux trapped can be estimated under adiabatic conditions. As the free energy F m is quadratic in the magnetic field components B i (r), the magnetic energy F m mag /T = 1 2T B i (r)B i (r)d n r can be written in thermal equilibrium as where (G m ) −1 ij (r, r ′ ) is the inverse of G m ij (k) defined in (68) and so the two-point interaction potential. In the fully two-dimensional model the thermal fluctuations are not correlated, (G 2D ) −1 (r, r ′ ) = δ 2 (r − r ′ ). In that case, the equilibrium value of (69) can be exactly computed on a surface of size S in function of the total magnetic flux Φ E.M. . That is
In the superconducting model (SC) and the fully-three dimensional model, because the thermal fluctuations are correlated, the integral in (69) does depend on the shape of the surface S. Nevertheless the magnetic energy is a quadratic function of the total flux in both cases: The factors of proportionality which are missing are geometrical coefficients of order one.
As the free energy is quadratic in the magnetic flux in thermal equilibrium, the thermal fluctuations of Φ E.M. adjust to the gaussian distribution ∼ exp (− 2T S ) in twodimensions and to ∼ exp (− 2T √ S ) in three dimensions. Correspondingly, the typical value of the thermal flux in a surface of area S is the root-mean-square of Φ E.M. : Φ rms,3D E.M.
The typical size S of a region of symmetric phase enclosed by bubbles can be estimated using dimensional arguments. The characteristic length of that region must be of the order of the bubble radius. If the transition is strong first order, the radius is determined by the nucleation rate Γ and the expansion rate v. S ∼ (v/Γ) 2/3 holds in the fully two-dimensional and the superconducting models while S ∼ (v/Γ) 1/2 does in the fully three-dimensional model. Making use of the estimations above we can compute the typical winding number of vortices and cosmic strings in each case: A variety of windings is expected in the actual distribution. The properties of the winding numbers distribution function are to be reported in [34].

Flux trapping vs. Kibble mechanism
As it was mentioned in section 3.2 (see also section 7.3), Kibble mechanism can be used in the form explained there to explain the formation of local topological defects in absence of magnetic field. On the other hand, because the typical winding predicted by Kibble mechanism is of the order of the square root of the number of bubbles involved in a typical collision, it must be of order one. In particular, if the width of the bubble wall is negligible, it is exactly one. That is so because the probability for more than three bubbles to collide at a single point is zero. Therefore the trapping of thermal flux will be dominant over Kibble mechanism in the formation of vortices as N w,typ. > 1. In order to estimate when that is so we have to find out a relation between the parameters involved in the equations (75), (76) and (77). First, we will assume that the electromagnetic coupling e is of order one. That is a reasonable assumption in both GUTs and electromagnetism. † Then, we have to compare the value of the temperature T to the value of some power of v Γ which is equivalent to some power of the radius of the bubbles. As mentioned at the beginning of this section, Γ is difficult to estimate and so is v. We will take a more conservative approach. If the transition is not weak Γ is exponentially suppressed either by exp [−F (φ T )/T ] if the transition is temperature driven or by exp [−S E (φ 0 )] if it is a quantum transition. S E (φ 0 ) stands for the Euclidean action of the bounce solution at zero-temperature while F (φ T ) does for the free energy of the bounce solution at finite-temperature T . Also, the size of the bubbles when they meet must be larger than or equal to the initial nucleation radius at the transition temperature R N (T ). Following Linde's arguments [31], If the transition is temperature driven, R N ≥ T −1 . Therefore, the above conditions together with our estimates in equations (75), (76) and (77) imply that N w > 1 is guaranteed if the transition is strong first order and temperature-driven. For the theory in (64) that is the case in the couplings range e 2 > λ > e 4 . Otherwise not only the nucleation process but also the bubble expansion rate must be considered.

The gauge-invariant phase difference
In section 2.2 it was argued that it is D µ θ the field to be tracked in the gauge theory in order to describe the formation of topological defects. The definition of D µ θ in (10) can be inferred from the application of the differential operator D µ onto the product exp (iθ) · |φ|. Following the Leibnitz's rule, On the other hand, the LHS of (80) can also be written as By expanding the RHS of (81) and using the fact that the covariant derivative of a neutral function like |φ| is equivalent to its ordinary derivative in (80), identity (10) holds.
The local observables can be written in terms of the gauge-invariant quantities |φ| and D µ θ as: charge density q e = 2e|φ| 2 D 0 θ, where in this case arrows denote spatial components vectors. Likewise, as ∂ µ ∂ ν θ = ∂ ν ∂ µ θ, it is possible to write down the equations of motion in terms of |φ| and D µ θ alone: where in the last equation we can identify In the global theory D µ θ reduces to ∂ µ θ and the phase difference b a ∂ x θ dx = α B − α A ≡ ∆α ab is meaningful. In the gauge theory, when no magnetic field is present and the gauge A x = 0 can be chosen, γ ab = b a D x θ dx = ∆α ab as well. In this case, the only reason for the existence of non-zero D µ θ is the spatial overlap between two different vacuum domains A and B as it occurs in the junction between two colliding bubbles. In absence of magnetic field, D µ θ can be interpreted as an induced spontaneous connection as a result of the coupling of the light ('fast') degrees of freedom living in the junction to the heavy ('slow') degrees of freedom living deep in the broken phase where m γ = √ 2e|φ| min . Even in the actual gauge theory this spontaneous connection cannot be gauged away. One can nevertheless perform the gauge transformation in such a way that D x θ remains invariant and so does γ ab . If a magnetic field B i does exist there is an extra divergenceless piece of connection induced by B i . It is always possible to write D i θ as a sum of an exact derivative ∂ i θ plus a divergenceless component eA i in the Coulomb gauge: As dDθ = edA = eF , ∂ i θ is the locally flat component of the connection while A i carries the curvature in the gauge (91). In summary, the local connection D i θ is composed of two parts: •A spontaneous, flat component ∂ i θ which, in the Coulomb gauge, carries the longitudinal modes of D i θ. Its origin is in the breakdown of the U(1) symmetry. In the global theory, its interaction with the neutral scalar particles gives rise to the geodesic rule.
•An induced, curved component eA i which, in the Coulomb gauge, carries the transverse modes of D i θ. Its origin is in a net magnetic field. In the gauge theory, its interaction with the charged scalar particles gives rise to the Aharonov-Bohm effect.
Both the local quantities ∂ i θ and eA i and their integrals ∂ i θ dx i and e A i dx i come to be respectively equal (up to a sign) deep in the broken phase where D i θ = 0. As one approaches the center of a vortex of winding N w at r = 0, D i θ is not zero any more for r r v ≈ √ N w λ L , where the magnetic flux is confined. The holonomy of any loop C of radius r enclosing the center of the vortex, C ∂ i θ dx i = 2πN w , does however remain unaltered as r goes to zero by continuity. It is the discrepancy between the holonomy and the magnetic flux in the circle of radius r v what gives rise to the closed steady screening currents in the core of the vortex.

DC Josephson-like effect across bubbles junctions and flux strips
To our knowledge Weinberg was the first to show the periodicity of the Josephson currents using field theory arguments [51]. This will be our starting point to study the electric current traversing the junction between two bubbles as they collide.
We will consider stationary situations and so neglect time oscillations, i.e. D 0 θ ≈ 0, ∂ t |φ| ≈ 0. Gauge invariance implies that, assuming no gradients along the junction between two superconducting pieces separated by a gap of width |b − a|, the value of the current self-interaction energy through the junction is an even function of where {c m } are constants which incorporate the degrees of freedom of |φ| in a < x < b.
If stationary conditions hold, the current density is [51] j x ab = −e ∂ε ab ∂γ ab .
If the current flow is parametrised with the coordinate x, a uniform steady current turns up in a < x < b [33] if linear superposition of the scalar fields is valid, that is, if 2/m h ≪ |b − a|. Uniformity of j x along the integration path a < x < b implies that the main contribution to formulas (92,93) comes from the terms m = 0 and m = 1. Uniformity of the reduced charge-carriers density |φ| 2 0 implies that c 0 ≈ c 1 ≈ −2|φ| 2 0 /χ. Therefore, under these conditions, (92) and (93) read respectively ‡ where |φ| 0 is given in (28). The potential (94) is a Josephson coupling and (95) is a direct Josephson-like current though its origin is not in quantum-tunneling [52]. Its amplitude, is instead a function of the reduced charge-carriers density and the length scale χ.
In the actual collision of two superconducting bubbles or blocks, steady and ‡ A more precise computation of this equation can be found in [33]. uniformity conditions approximately hold at the middle of the junction. The charge carriers density |φ| 2 is roughly uniform at the very middle of the gap, x = 0, where |φ| presents a minimum and so ∂ x |φ| x=0 = 0 ( figure 17(b)). The non-zero value of D x θ concentrates also in a very narrow strip along the junction with center at x = 0 and width χ. § When ∂ z |φ| is negligible the current flows across the junction in the xdirection. Figure 15(a) shows the shape of the transverse current density j x ab in an actual simulation. It is in good agreement with expression (95). The current is measured at the middle of the gap (x = 0) between the two superconducting blocks of figure 8. There is a magnetic flux of value 6π/e uniformly distributed in the gap so that the gradient of γ ab is uniform in the direction z along the junction. Figure 15(b) shows the corresponding current self-interaction energy as a function of γ ab normalized to 2 χ |φ| 2 0 = 1. An interesting scenario in which an actual steady superconducting current can be found is a meta-stable flux strip. In general, neither the density of charge carriers nor the magnetic flux are uniform in a meta-stable strip like the one in figure 14 with µ = 3. The gradients of |φ| and D x θ across the strip make equations (94,95) inapplicable. Instead, the general expression (92) must be used. In addition, the gradients along the strip must be incorporated in the coefficients of the sum. Therefore, ε ab (z) takes the general form Correspondingly, the transverse current density is The functions {C m (z)} are not completely arbitrary. All the physical quantities are periodic along a meta-stable strip. If the length of the strip is L z and it contains approximately N w winding units of magnetic flux, their period is ∆z ≈ L z /N w . Consequently, both the functions {C m (z)} and γ ab (z) have period ∆z along the junction (modulo 2π in the case of γ ab ). In terms of the total magnetic flux Φ, ∆z = 2πLz eΦ . A first consequence of the current density distribution (98) along a flux strip is an interference pattern of the transverse total current different to that in an actual Josephson junction. Taking into account the symmetry and periodicity properties of the functions in j x ab (z), it can be formally integrated along the junction of a meta-stable strip in a superconductor film of thickness µ and total flux Φ. That is where the coefficients C ′ m are constants and ∆α ab stands for the spontaneous phase difference. The thickness µ is taken large enough such that j x ab is approximately uniform in the direction normal to the film. Figure 16 shows the functional form of the first three harmonics.
The meta-stability of a flux strip requires the gap width |b − a| to be relatively broad. Broad means in the first place that |b − a| is greater than 1/m h so that the edges of the strip are approximately straight and the linear magnetic flux density is nearly uniform. Broad also implies that χ ∂ x |φ| dx ≪ χ −1 χ |φ| dx, χ ≪ |b − a|. If these conditions are satisfied the density of charge carriers is nearly uniform in the range of interest. Therefore the first harmonics m = 0, m = 1 are the dominant terms in the sums (97-99).
Because the total current only vanishes for integral values of the total flux in units of 2π/e, there must be a net current passing through the gap during the process of formation of a meta-stable flux strip. That total current stops when the flux in the junction gets quantized. It is worth noting that the remaining current density j x ab (z) along the strip does not radiate because of topological reasons. It is the flux quantization rather than the existence of a mass gap what prevents this current from dissipating.
The agreement between the lattice simulation and the Josephson-like current in (95) deserves some comments. At the moment the blocks get in contact is of order π and, as we mentioned in section 4, the dynamics of γ ab is governed by a sine-Gordon equation. In the numerical simulation, this apparent problem is overcome naturally with the use of the non-compact U(1) lattice formulation [53]. As described in Appendix B, it is the use of the link field U i = exp (ieδxA i ) that incorporates the nonlinear effects. Hence, the Josephson-like coupling − 2

The geodesic rule in gauge theories
In [54] it was argued that the geodesic rule in a gauge theory holds by equilibrium arguments, i.e. as a result of the minimization of the energy. In [55] it was shown that if no magnetic field is present the geodesic rule holds in the gauge theory by dynamical arguments, i.e. as a result of the equations of motion. A more detailed study of the bubble collision dynamics was carried out in [26] in terms of gauge-invariant quantities. We will show with the aid of numerical simulations that both energetic [54] and dynamical [26] arguments hold and agree with the geodesic-magnetic rule presented in section 3.2. The global geodesic rule as proved in [55] is only satisfied if some particular prescription for the gauge choice can be given (see below).
Let us consider two superconducting bubbles A and B facing each other. Once more, if the radius of the bubbles is much larger than the bubble-wall thickness the contact area can be well approximated by two superconducting blocks with a gap in between. Block A is nucleated at the point α A while block B is at the point α B of the vacuum manifold. The gap in between is initially in the symmetric phase and no magnetic flux at all is in it. Let us take as an example ∆α ab = α B − α A = 3π/4 and let the two blocks expand and coalesce. Figure 17  γ ab and ∆θ ab across the junction a < x < b. In agreement with the geodesic-magnetic rule, γ ab decreases in time until vanishing. In the contrary, ∆θ ab remains constant and equal to 3π/4. The complex phase θ does interpolate between 0 and 3π/4 as stated by the geodesic rule because the gauge A µ = 0 has been chosen initially. However, at the moment the blocks meet we could have made the following gauge transformation everywhere: The complex phase would have been uniform in both the blocks and the junction. In that case the evolution of γ ab would have still been the one shown in figure 17(a) while ∆θ ab would have remained equal to zero. This shows that the physical modes are those of D x θ, neither the Goldstone modes ∂ x θ nor the longitudinal modes of the gauge field A x . When the temporal gauge A 0 is chosen, the evolution of ∂ x θ and A x just depends on the initial fixing of the remaining gauge degrees of freedom. Thus, the geodesic rule as applied in the global theory is a gauge artifact in the gauge theory. It holds in the gauge theory when the gauge choice A µ = 0 can be made [55], which is only possible in the absence of a magnetic field. One can however consider that in the absence of a magnetic field topological defects form by Kibble mechanism. After all, in the gauge theory it is also the spontaneous breakdown of the symmetry that gives rise to the D i θ modes.
In energetic terms what is happening is that, as bubbles coalesce, the current tends to vanish as γ ab goes to zero. The current could also vanish if γ ab evolved towards π. However π corresponds to a maximum and so unstable point of the current selfinteraction energy ε ab ( figure 15(b)). The current self-interaction potential ε ab does decrease as γ ab evolves toward its minimum at γ ab = 0.
In dynamical terms, the Maxwell equation (89) becomes a generalized Klein-Gordon equation for D µ θ in the gauge ∂ ν D ν θ = const., Integrating the resultant equation for D x θ across the junction (using the rough approximation |φ| ∼ uniform), the equation for γ ab is The study of the time-evolution of γ ab was carried out in detail in [26] taking |φ| ≈ |φ| min , ∂ x |φ| ≈ 0. However, because the typical phase difference between any two bubbles in a multi-bubble collision is of the order of π for a vortex to form, (103) is not a good approximation. It is instead the sine-Gordon equation (27) what should be used in the first stages of the collision. It is only when γ ab relaxes to γ ab ≪ π that (103) is valid and γ ab goes roughly as ∼ (α B − α A )e − 2e 2 |φ| 2 σ t . The numerical simulation of figure 17(a) so confirms. ‡ Figure 18 shows the evolution of the magnetic fluxon ring generated around the ‡ In the simulation of figure 17(a) the dynamics of γ ab is overdamped so that oscillations in time are suppressed.
collision axis when two bubbles meet with α B − α A = 3π/4. Because the collision there is soft, the fluxon can spread out in the directions orthogonal to the collision axis. That is the reason for its constant asymptotic value.

Conclusions and outlook
In the Abelian Higgs model we have shown that it is the local dynamics of the gaugeinvariant field D i θ = ∂ i θ + eA i and the global dynamics of the gauge-invariant phase difference γ = D i θ dx i that determine the formation of vortices in superconductors and hypothetical cosmic strings in grand unified theories.
In the mean-field theory approximation we have found that, when a first order phase transition takes place in an external magnetic field, a variety of objects of different winding form. These are meta-stable strips and membranes of magnetic flux, chains of single-winding vortices, thick vortices and heavy strings. Thick vortices turn out to be stable in thin superconducting films up to some critical winding value. The cause of the stability and/or meta-stability of vortices and strips is a Coulomb-like long-range magnetic interaction induced by the gauge field which propagates out of the film. It might determine the order of the phase transition of a thin superconducting film. A Monte Carlo simulation similar to the one in [18] with an external magnetic field is needed to investigate this matter.
The spontaneous breakdown of the U(1) gauge symmetry manifests in the location of the center of mass of single-winding vortex chains and finite size Abrikosov lattices.
Under steady conditions we have derived a direct Josephson-like current traversing the meta-stable flux strips. For that calculation we have invoked non-local effects within a length scale χ which we identify with the correlation length of D i θ.
Thermal long-wavelength fluctuations of the gauge fields get confined in a variety of winding topological defects in a first order phase transition. In the particular case of the fully three-dimensional Abelian Higgs model, a multi-tension cosmic string network arise [18]. The typical winding of the strings is ∼ e 2π T 1/2 (v/Γ) 1/8 if the transition is strong first order, e 2 > λ > e 4 . For a thin superconducting film the prediction is ∼ e 2π T 1/2 (v/Γ) 1/6 . Even though the latter is not easy to test experimentally, it is our belief that trapping of an external magnetic flux in the form of strips, vortex chains and thick vortices can be achieved experimentally by inducing nucleation of bubbles in a type-I superconducting film sample.
In the context of brane-world cosmology and extra dimensions phenomenology, it would be interesting to study the hypothetical non-Abelian currents of grand unified theories traversing a 3 + 1 dimensional brane (with or without fluxes) in between two bulk disconnected -or weakly overlapped -vacuum domains.
If the collision is violent, the highly non-equilibrium dynamics of the scalar field can lead to the formation of unit-winding vortex-anti vortex pairs [56] giving rise to an apparent violation of the geodesic rule.

Acknowledgments
I would like to thank very specially Arttu Rajantie for his useful advice, discussion and corrections throughout all this paper. I thank Tom Girard for his comments on section 5 and discussion about experimental issues in superconductivity. I am also grateful to the referee for pointing out an important error in Section 4.
Appendix A. Breakdown of the translational invariance of the setup in figure 8 Let us consider the fully two-dimensional theory and take two parallel superconducting blocks separated by a finite gap as in figure 8. The blocks are semi-infinite in the transverse direction x and infinite along z. A non-zero magnetic field exists in the gap, x ∈ [−d/2, d/2], and it is invariant along the z-axis. Before the blocks get in contact, neglecting the tension of the walls, the initial equilibrium configuration for the magnetic flux corresponds to that of a uniform magnetic field B. If either the walls tension were considered or we were working in the framework of the superconducting model, translation invariance along the junction would also hold but the magnetic field strength would show a gradient across the gap. We aim to show that, regardless of the choice of gauge: a) local translation invariance along the junction does hold as long as a finite width gap in which |φ| = 0 ∀z exists between the two blocks, b) local translation invariance breaks down as the superconducting blocks get in contact. Condition a) is necessary for both SSB and for the stability of the straight strip in figure 9(a). Condition b) is necessary for the formation of the structures presented in figures 9(b) and 9(c).

Proof of part a).
The initial conditions of the local physical quantities are: 1) Translation-invariance along the z direction: The above equations also imply: ∂ z j i | t=0 = 0, at x such that ∂ z |φ| x,t=0 = 0, ∂ z D i θ| t=0 = 0, at x such that ∂ z |φ| x,t=0 = 0 and |φ| > 0. (c) t = 0. Initial profile of the square of the current density j 2 ≡ j 2 x + j 2 z . Invariance along z holds. The current density concentrates at the edges of the blocks where the magnetic field initially migrates due to small structure effects. (d) t = 30. z-invariance breaks down as blocks meet. j 2 gets concentrated around the minima of |φ|.

Appendix B. Lattice disretization
To write down the discretized equations of motion, we define the forward and backward derivatives where δx is the lattice spacing,î is a unit vector in the i direction and vectors are denoted by bold-faced letters. Similarly, we define the time derivative where δt is the time step. Using the link variable . Figure  ,where l x × l z × l y is the lattice size in lattice spacing units; (pbc) stands for 'periodic boundary conditions'; δx is the lattice spacing; δt is the time step; σ is the damping rate; B ini is the initial magnetic flux density; µ is the film thickness in the superconducting model; d stands for bubble centers separation; d s is the initial square side; R N is the initial bubble radius; |b − a| is the gap size between blocks or bubbles. The value for σ in fig.no. 7 is 0.0 in the equation of the gauge field and is variable in the equation for φ as it determines the bubble expansion rate.