Magnetic field production at a first-order electroweak phase transition

TL;DR

This study simulates the generation of magnetic seed fields during a first-order electroweak phase transition, revealing that about 10% of latent heat converts into magnetic energy, mainly during the Higgs oscillation phase.

Contribution

It provides a numerical analysis of magnetic field production during a first-order electroweak phase transition using lattice simulations, highlighting the energy conversion and spectrum evolution.

Findings

Approximately 10% of latent heat becomes magnetic energy.

Magnetic field spectrum peaks shift to larger scales during transition.

Magnetic fields are mainly generated during the Higgs oscillation stage.

Abstract

We study the generation of magnetic field seeds during a first-order electroweak phase transition, by numerically evolving the classical equations of motion of the bosonic electroweak theory on the lattice. The onset of the transition is implemented by the random nucleation of bubbles with an arbitrarily oriented Higgs field in the broken phase. We find that about 10% of the latent heat is converted into magnetic energy, with most of the magnetic fields being generated in the last stage of the phase transition when the Higgs oscillates around the true vacuum. The energy spectrum of the magnetic field has a peak that shifts towards larger length scales as the phase transition unfolds. By the end of our runs the peak wavelength is of the order of the bubble percolation scale, or about a third of our lattice size.

Tables

Table 1. Table 1: The magnetic energy density measured in units of m H 4 superscript subscript 𝑚 𝐻 4 m_{H}^{4} and the conversion rate ρ B / ρ total subscript 𝜌 𝐵 subscript 𝜌 total \rho_{B}/\rho_{\text{total}} for p B = 10 − 4 subscript 𝑝 𝐵 superscript 10 4 p_{B}=10^{-4} .

$\gamma$	${\rho }_B/m_H^4$	${\rho }_B/{\rho }_{\text{total}}$
$0.00$	$0.0160$	$0.0331$
$2.34\times {10}^{-5}$	$0.0162$	$0.0335$
$1.00\times {10}^{-4}$	$0.0160$	$0.0332$
$1.00\times {10}^{-3}$	$0.0139$	$0.0301$
$1.00\times {10}^{-2}$	$0.00842$	$0.0233$

Table 2. Table 2: Magnetic energy density at t ∼ T stop similar-to 𝑡 subscript 𝑇 stop t\sim T_{\text{stop}} and at the end of our simulation.

$\gamma$	$p_B$	${\rho }_{B,BC}/m_H^4$	${\rho }_{B,H{O}^{\ast }}/m_H^4$
$2.34\times {10}^{-5}$	${10}^{-4}$	$0.016$	$0.062$
$2.34\times {10}^{-5}$	${10}^{-6}$	$0.016$	$0.060$
$2.34\times {10}^{-5}$	${10}^{-7}$	$0.015$	$0.060$
$1.00\times {10}^{-2}$	${10}^{-6}$	$0.0075$	0.037

Equations

$$\mathcal{L}=-\frac{1}{4}{W}_{\mu \nu }^a{W}^{a\mu \nu }-\frac{1}{4}{B}_{\mu \nu }{B}^{\mu \nu }+|D_{\mu }\Phi {|}^2-\lambda (|\Phi {|}^2-{\eta }^2{)}^2,$$

$$D_{\mu }={\partial }_{\mu }-\frac{i}{2}g{\sigma }^a{W}_{\mu }^a-\frac{i}{2}{g}^{\prime }{B}_{\mu }.$$

$${\partial }_0^2\Phi =$$

$${\partial }_0^2{B}_i=$$

$${\partial }_0^2{W}_i^a=$$

$${\partial }_0{\partial }_j{B}_j$$

$${\partial }_0{\partial }_j{W}_j^a$$

$$\ddot{r}-r{\dot{\theta }}^2+c\dot{r}+r$$

$$r^2\dot{\theta }$$

$$\ddot{\mathbf{r}}+c\mathbf{r}\frac{d}{dt}\ln r+\mathbf{r}=0,$$

$${\partial }_0^2{\Phi }_a$$

$$a=1,2,3,4,$$

$$\text{Im}[{\Phi }^{†}\Xi ]$$

$$\text{Im}[{\Phi }^{†}{\sigma }^i\Xi ]$$

$$\Xi \propto \Phi ,$$

$$0=4\pi {\int }_0^{\infty }{r}^{2}dr\left[({\partial }_{r}f{)}^2-2\lambda {\eta }^2{f}^2+\lambda f^4\right].$$

$${\partial }_{r}f=-\sqrt{\lambda }(2{\eta }^2-f^2{)}^{1/2}f.$$

$${\partial }_{\rho }F=-(1-F^2{)}^{1/2}F$$

$$F(\rho )=\frac{2C{e}^{-\rho }}{1+C^2{e}^{-2\rho }}.$$

$$f(r)=\eta \frac{[1+(\sqrt{2}-1{)}^2]e^{-m_{H}r/\sqrt{2}}}{1+(\sqrt{2}-1{)}^2{e}^{-\sqrt{2}{m}_{H}r}}$$

$$A_{\mu }=\sin {\theta }_w{n}^a{W}_{\mu }^a+\cos {\theta }_w{B}_{\mu },$$

$$n^a\equiv -\frac{{\Phi }^{†}{\sigma }^a\Phi }{{\eta }^2}$$

$$A_{\mu \nu }$$

$$-i\frac{2}{g{\eta }^2}\sin {\theta }_w\left[(D_{\mu }\Phi {)}^{†}(D_{\nu }\Phi )-(D_{\nu }\Phi {)}^{†}(D_{\mu }\Phi )\right].$$

$$⟨B_i^{\ast }(\mathbf{k},t)B_j({\mathbf{k}}^{\prime },t)⟩=(2\pi {)}^3{\delta }^{(3)}(\mathbf{k}-{\mathbf{k}}^{\prime })F_{ij}(\mathbf{k},t),$$

$$B_i\left(\mathbf{k},t\right)$$

$$B_i\left(\mathbf{x},t\right)$$

$$\frac{F_{ij}(\mathbf{k},t)}{(2\pi {)}^3}=({\delta }_{ij}-{\hat{k}}_i{\hat{k}}_j)\frac{E_M(k,t)}{4\pi k^2}+i{ϵ}_{ijl}{k}_l\frac{H_M(k,t)}{8\pi k^2}.$$

$${\rho }_B\left(t\right)=\frac{1}{2}⟨{\mathbf{B}}^2(\mathbf{x},t)⟩={\int }_0^{\infty }{E}_M\left(k,t\right)dk.$$

$$k_{\text{mean}}\left(t\right)=\frac{{\int }_0^{\infty }k{E}_M\left(k,t\right)dk}{{\int }_0^{\infty }{E}_M\left(k,t\right)dk},$$

$$k_i=2\pi K_i^{\prime }/L,$$

$$K^{\prime}_{i}=\left\{\begin{array}[]{ll}K_{i}&K_{i}\leq N/2\\ K_{i}-N&K_{i}>N/2.\\ \end{array}\right.$$

$${\lambda }_k=2\pi /k=N\Delta x/K^{\prime },$$

$${\rho }_B(t)=\sum _{K_c^{\prime }}{E}_M(K_c^{\prime },t)\Delta K^{\prime },$$

Full text

I Introduction

Magnetic fields are pervasive in the Universe. Micro-gauss fields coherent on scales up to ten kpc have been detected in nearby spiral galaxies, such as the Milky Way and in higher redshift galaxies (e.g. Kronberg (1994)). Similar stochastic magnetic fields are found in clusters of galaxies. Such fields are believed to result from the dynamo amplification of weak magnetic field seeds, whose origin remains a mystery (see e.g. Kandus et al. (2011); Durrer and Neronov (2013); Subramanian (2016)). Recent observational evidence that even the intergalactic medium in voids is pervaded by a weak (below femto-gauss strength) magnetic field Neronov and Vovk (2010); Tavecchio et al. (2010); Tashiro et al. (2014); Chen et al. (2015a, b); Biteau (2018), points to a primordial origin of the field seed, since it is difficult to account for a field that fills the volume in the void regions by astrophysical processes in the late universe Furlanetto and Loeb (2001); Bertone et al. (2006) (but see Durrive and Langer (2015); Beck et al. (2013)). Furthermore, the study of the diffuse gamma-ray background provides clues for the existence of a nonzero helical component in the intergalactic magnetic field Tashiro et al. (2014); Chen et al. (2015a); Duplessis and Vachaspati (2017). Interestingly, while most magnetic field generation mechanisms discussed so far produce nonhelical fields, physical processes associated with electroweak baryogenesis imply fields that are helical Cornwall (1997); Vachaspati (2001). This is important, since conservation of magnetic helicity plays an important role in the evolution of primordial fields and leads, via inverse cascade, to a larger coherence scale than for nonhelical fields Christensson et al. (2001).

The production of helical fields is related to the violation of baryon plus lepton (B+L) number during the electroweak phase transition (EWPT). To change B+L requires a change in the Chern-Simons number of the electroweak gauge fields, which proceeds through “sphaleron” configurations. The decay of the sphaleron releases helical magnetic fields, and this is borne out by both numerical simulations Copi et al. (2008) and analytical arguments Chu et al. (2011). Generically, helical magnetic fields are produced at the electroweak phase transition if the Chern-Simons number changes during the phase transition Zhang et al. (2017). Inhomogeneities in the Higgs field can also give rise to the generation of helical magnetic fields as shown in Diaz-Gil et al. (2008a, b) in the context of low-scale electroweak hybrid inflation Garcia-Bellido et al. (2003).

The study of the properties of cosmological magnetic fields can, thus, open a window to the early universe and very high energy particle interactions. To investigate the epoch of magnetogenesis and identify the mechanism at the origin of the primordial seeds we need to make full use of the statistical properties of the magnetic fields matching theoretical predictions with observations Caprini and Gabici (2015). These are encoded in the power spectrum which, for Gaussian random fields, entirely describes the statistical properties of the fields. In this respect, the helicity spectrum can be recovered from observations of the diffuse gamma-ray background Chen et al. (2015a).

On the theoretical front, it has been shown that a measurement of the sign of the helicity can distinguish fields produced during electroweak baryogenesis from those generated at an earlier epoch associated with leptogenesis Long et al. (2014). However, a detailed study of the statistical properties of the magnetic fields is required to narrow the gap between theoretical studies of the microscopic properties and the macroscopic properties that can be experimentally observed (field strength and coherence scale), as was done for the case of low-scale hybrid inflation in Diaz-Gil et al. (2008a, b).

In this paper, we study the dynamics of a first-order EWPT by numerically evolving the classical bosonic electroweak theory, and we examine the properties of the generated magnetic field. In the Standard Model (SM), the EWPT is first order only if the mass of the Higgs boson lies below $m_{h}\lesssim 70$ GeV Bochkarev and Shaposhnikov (1987); Kajantie et al. (1996). Since the experimentally observed mass violates this bound Aad et al. (2012); Chatrchyan et al. (2012), we work under the assumption that there is physics beyond the SM that influences the character of the phase transition. Although the details of the particular SM extension are not important for our purpose, we have in mind models that make electroweak baryogenesis viable (see e.g. Riotto and Trodden (1999); Quiros (1999); Dine and Kusenko (2003); Morrissey and Ramsey-Musolf (2012) for reviews). In that case, the properties of the cosmological magnetic fields can be related to the observed baryon asymmetry of the universe Vachaspati (2001). Furthermore, a cosmological first-order phase transition can also be a source of gravitational waves (GWs) Witten (1984); Hogan (1986). In this respect, large-scale numerical simulations of a scalar field theory on the lattice have recently become available to verify the production of GW radiation during a first-order phase transition Cutting et al. (2018); Ellis et al. (2019); Kharzeev et al. (2019). Our results bring another handle to probe the universe at the time of the EWPT: the observation of cosmological magnetic fields.

When the phase transition occurs, the Higgs field will leave the symmetric phase and gradually settle down around the true vacuum, $|\Phi|=\eta$. In order to mimic this behavior on the lattice, we introduce a phenomenological damping term for the Higgs boson. The modified electroweak evolution equations preserve gauge invariance and satisfy the Gauss constraints. They are reviewed in Sec. II, where we also outline the calculation of the magnetic field spectrum from our lattice simulations.

Magnetic fields are produced as a result of nonvanishing gradients of the Higgs field Vachaspati (2001). In a first-order phase transition, bubbles are randomly nucleated as a result of quantum tunneling, and subsequently expand and collide with each other, producing an out-of-equilibrium environment. Magnetic fields, will be produced when bubbles collide. We argue that the details of the short period of bubble nucleation are not of major concern, and we can safely mimic the quantum tunneling by a simple random procedure, controlled by a parameter $p_{B}$, the nucleation probability. The specifics of our numerical implementation are presented in Sec. III.

To gain intuition about the general features of the magnetic fields generated from bubble collisions, we first consider in Sec. IV a set of constrained simulations with a regular array of bubbles so that we can control the size of the colliding bubbles. In Sec. V we allow for bubbles to nucleate at random locations in the unbroken phase and we present the spectrum of the magnetic field induced during a first-order EWPT. At the end of the bubble collision stage, we find that the field strength is about $10^{23}$ gauss, with a comoving correlation length of about 2 cm. If the fields are nonhelical, the field strength will evolve to $\sim 10^{-15}$ gauss at present Banerjee and Jedamzik (2004). Finally, we discuss and summarize our results in Sec. VI.

II Theoretical Formulation

II.1 Classical equations and Higgs damping

We consider the classical bosonic electroweak theory, which includes the Higgs doublet $\Phi$, the SU(2)-valued gauge fields $W_{\mu}^{a}$ and the U(1) hypercharge field $B_{\mu}$. The Lagrangian is given by

[TABLE]

where

[TABLE]

Since we will evolve the field equations using the standard Wilsonian approach for lattice gauge fields Ambjorn et al. (1991); Moore (1996); Tranberg and Smit (2003); Garcia-Bellido et al. (2004), it is advantageous to use the temporal gauge, $W_{0}^{a}=B_{0}=0$, which allows a simple identification of the canonical momentum. Our implementation of the lattice equations can be found in Rajantie et al. (2001); Zhang et al. (2017), with the addition of a Higgs damping term that we now discuss. The classical equations of motion (EOMs) in the continuum are given by:

[TABLE]

The solutions are subject to the Gauss constraints:

[TABLE]

The term proportional to $\gamma$ in Eq. (3) gives rise to a linear damping of the magnitude of the Higgs field. Its purpose is to attenuate the radial oscillations of the Higgs field within a reasonable simulation time. The term also simulates the effects of energy losses due to Higgs decays into fermion-antifermion pairs. In principle, one could include separate damping terms for each component of the Higgs field (e.g. to account for the decays of gauge bosons). However, this is harder to implement since we want to respect gauge invariance and the conservation of electric charge.

The form of the damping term in Eq. (3) is motivated by considering the analogy of a 2D simple harmonic oscillator with linear damping in the radial direction. In polar coordinates, the EOMs are

[TABLE]

where, $r=|{\bm{r}}|$, ${\bm{r}}=\left(x,y\right)$. The second equation is the conservation of angular momentum and corresponds to the conservation of charge in the field theory. If the two equations are transformed into Cartesian coordinates, we find

[TABLE]

which suggests the modification of Eq. (3).

To check that the damping term is not in conflict with the Gauss constraints, let us write a general additional term for the Higgs EOM,

[TABLE]

where $\Xi=\left(\Xi_{1}+i\Xi_{2},\Xi_{3}+i\Xi 4\right)^{T}$ is a doublet that causes damping. The Gauss constraints should be satisfied throughout the evolution assuming that they are initially satisfied, which requires that the time derivatives of Eqs. (6) and (7) should vanish. Differentiating the two equations, and using the modified EOM Eq. (11) for the Higgs field, together with Eqs. (4) and (5), to replace the second-order time derivatives, we obtain

[TABLE]

The solution to these equations, which ensures the Gauss constraints, is

[TABLE]

and this gives the damping term in Eq. (3).

To simulate a first-order electroweak phase transition we also need to include the quantum nucleation of broken phase bubbles, which we now discuss.

II.2 Bubble profile

As mentioned above, the SM does not admit a first-order electroweak phase transition. However, we expect that for our purposes the Higgs potential as given in Eq. (1) captures the dynamics of the phase transition in an extension of the SM that allows for a viable electroweak baryogenesis. The equations of motion, with the initial condition $\Phi=\dot{\Phi}=0$, are then supplemented by bubble nucleation events that occur randomly in regions where the symmetry is unbroken.

In numerical simulations, we use a simple method for randomly nucleating bubbles of the broken phase to mimic the tunneling effect. The sites of the lattice can be numbered as a sequence $S$ from 1 to $N^{3}$. At each time step of the simulation, we first randomly shuffle this sequence, and then we select sites from this shuffled sequence $S_{\text{rand}}$, with probability $p_{B}$ that controls the rate of nucleation on the lattice. The selected sites are stored $S_{\text{select}}$, which contains $\sim p_{B}N^{3}$ elements. We scan each site $s_{i}$ in $S_{\text{select}}$ and if the Higgs field in all lattice sites within a radius $r_{0}$ from $s_{i}$ is still in the symmetric phase, then a new bubble centered at $s_{i}$ is nucleated at this time step, otherwise, the site $s_{i}$ is skipped. The random shuffle procedure of all sites guarantees that the nucleation procedure is unbiased in any direction.

At every nucleation event we set up a bubble with a spherically symmetric profile $|\Phi|=f(r)$, determined by demanding that the nucleation process conserves energy. This requires that the energy change due to bubble nucleation should vanish. Therefore,

[TABLE]

This equation can be satisfied by choosing

[TABLE]

Rescaling $\rho=\sqrt{2\lambda}\,\eta r$ and $F=f/(\sqrt{2}\,\eta)$ gives

[TABLE]

and the solution is

[TABLE]

The integration constant, $C$, is fixed by requiring that the center of the bubble be in the true vacuum, so $f(0)=\eta$ ($F(0)=1/\sqrt{2}$). This leads to $C=\sqrt{2}\pm 1$, and we choose the smaller value so as to have a gentler bubble profile. The final solution for the bubble profile is

[TABLE]

where $m_{H}=2\sqrt{\lambda}\,\eta$ is the Higgs mass. 111One issue with this bubble profile function is that it has a kink at $r=0$. For example, along the $x$ axis and close to the origin it behaves as $\exp(-|x|)$. This is not a problem numerically as finite differences will not resolve the kink.

In this way we have fixed the bubble profile using the (stronger) requirement that bubble nucleation conserve energy locally. The direction of the Higgs within the bubble is assumed to be uniform. Since the vacuum manifold – zeros of the Higgs potential – defines a three sphere, the direction of the Higgs is chosen by randomly selecting a point with a uniform distribution on the three sphere. Different bubbles will have different Higgs field orientations.

II.3 Definition of electromagnetic field

Once the Higgs field leaves the symmetric phase ($\Phi=0$), we can define the electromagnetic field as follows:

[TABLE]

where,

[TABLE]

indicates the direction of the Higgs field, and $\theta_{w}$ is the weak mixing angle. The corresponding field strength is constructed as ’t Hooft (1974); Vachaspati (1991)

[TABLE]

II.4 Magnetic energy spectrum

We assume that the magnetic field after production is a statistically homogeneous and isotropic, Gaussian-distributed vector field. The field can then be described in terms of the equal time correlation function, which we write, following the conventions in Brandenburg et al. (2017, 2018), as

[TABLE]

where $B_{i}({\bm{k}},t)$ is the Fourier transform of $B_{i}({\bm{x}},t)$ with the convention

[TABLE]

The spectrum, $F_{ij}({\bm{k}},t)$, can be divided into a symmetric (nonhelical) part and an antisymmetric (helical) part,

[TABLE]

The mean magnetic energy density can be written as

[TABLE]

The average wave number

[TABLE]

provides a characteristic of the energy distribution in a given field configuration.

We implement discretized versions of the expressions above on a three-dimensional lattice containing $N$ nodes separated a distance $\Delta x$ along each spatial dimension of length $L=N\Delta x$. Every lattice point is labeled by a triplet of integers, ${\bm{X}}$, each ranging from [math] to $N-1$: $X_{i}\in\left\{0,1,...,N-1\right\}$ , for $i=1,2,3$. The corresponding Fourier space is also described by triplets of integers, $K_{i}$, of the same form. The largest wavelength along a particular direction, corresponding to the smallest momentum, is of order $L$, while the smallest wavelength that can be effectively described is determined by the spacing $\Delta x$. Allowing for modes traveling in the positive and negative directions, a given physical wave number corresponds to a triplet $K^{\prime}_{i}$:

[TABLE]

where,

[TABLE]

For the physical wavelength associated to a momentum with magnitude $k$ we have

[TABLE]

where the magnitude, $K^{\prime}=|{\bf K^{\prime}}|$, ranges from 0 to $\sqrt{3}N/2$ on a cubic lattice. We divide this range into bins of size $\Delta K^{\prime}$, each with center value $K^{\prime}_{c}$.

Then Eq. (27) can be approximated by

[TABLE]

where,

[TABLE]

Also, the discrete Fourier transform is given by

[TABLE]

III Numerical Simulation

As mentioned above, we follow the strategy in Rajantie et al. (2001); Zhang et al. (2017) to evolve the electroweak EOMs on the lattice. Our code is based on the LATfield2222http://github.com/daverio/LATfield2 library Daverio et al. (2015), and the linear algebra operations are performed with the help of the Eigen333http://eigen.tuxfamily.org library Guennebaud et al. (2010). Our simulations use periodic boundary conditions and the dimensionless constants entering the EOMs are fixed to their physical values: $g=0.65$, $\sin^{2}\theta_{w}=0.22$, $g^{\prime}=g\tan\theta_{w}$ and $\lambda=0.129$. The spatial and time spacing are chosen to be $\Delta x=0.25$, $\Delta t=\Delta x/4=0.0625$, respectively. The dimensionful vacuum expectation value of the Higgs, denoted by $\eta$, is $174.13\leavevmode\nobreak\ {\rm GeV}$. In our numerical code we set $\eta=1$, so that $\eta\Delta x=0.25$, and then $m_{H}\Delta x=2\sqrt{\lambda}\eta\Delta x=0.18$ where $m_{H}$ is the mass of the Higgs. This choice of lattice spacing gives us enough resolution to ensure that we capture all the dynamics. For instance, since $m_{H}\Delta x=0.18$, momenta of order $m_{H}$ are well resolved. The bulk of our simulations is performed on a lattice with size $N=256$, although we use a larger lattice for several runs in Sec. IV. We denote by $T$ the (integer) time step number, and the physical time $t$ is $t=T\Delta t$.

The bubble profile function, Eq. (19), does not have any free parameters and its tail has infinite extent, which we truncate on the lattice as follows. We define the symmetric phase to correspond to locations where $|\Phi|\leq 0.01\eta$. With this prescription, the “size”, $r_{0}$, of the bubble turns out to be $\eta r_{0}=9.0$ ($m_{H}r_{0}\approx 6.5$), since the profile in Eq. (19) falls below $0.01\eta$ for $r>r_{0}$. With our lattice parameters, this gives $r_{0}$ to be $36\Delta x$. We use this value to prevent the nucleation of new bubbles within existing ones: a bubble can only be nucleated at a particular site if all lattice points within a distance $r_{0}$ are still in the symmetric phase ($|\Phi|\leq 0.01\eta$). Once a bubble is nucleated, it will expand and collide with other bubbles if there are any in the vicinity. The expansion of a single bubble is shown in Fig. 1, while Fig. 2 shows the evolution and collision of several randomly generated bubbles.

Two additional inputs required for our runs are the Higgs damping $\gamma$, defined in Eq. (3), and the nucleation probability $p_{B}$, which determines the probability of bubble nucleation per lattice site per time step. These two parameters cannot be determined within the model we are considering, and thus we will compare the results by varying the two parameters. We consider several values in the range $0\leq\gamma\leq 0.01$, including the experimentally measured decay width of Higgs boson, $\gamma\sim\Gamma_{\rm Higgs}\sim 4.07\times 10^{-3}\text{ GeV}$ Patrignani et al. (2016), which corresponds to $\gamma\sim 2.34\times 10^{-5}$ in our lattice units. $p_{B}$ is chosen to be in the range $10^{-8}\leq p_{B}\leq 10^{-3}$ in our simulations.

Since we are concerned with the generation of magnetic fields during the electroweak phase transition, we need a criterion to determine when the phase transition is completed. Our strategy is to compute the minimum $|\Phi|^{2}$ among all the lattice sites at each time step. To avoid spurious fluctuations, we work with the ten-step moving average of $|\Phi|_{\text{min}}^{2}$, denoted as $|\Phi|_{\text{MA10}}^{2}$, and we stop the simulation at the first time step $T_{\text{stop}}$ when $|\Phi|_{\text{MA10}}^{2}>0.25\eta^{2}$. In this manner, we ensure that the Higgs field is away from the symmetric phase.

One caveat of our formalism is that our field equations do not include the effects of other charged particles that might be present or generated at the time of the phase transition.

IV Test runs with nonrandom bubble distributions

A single expanding bubble in our analysis does not generate magnetic fields. This can be verified from the field equations since the gauge field currents [right-hand sides of Eqs. (4) and (5)] vanish for a spherically symmetric expanding bubble. Hence, the simplest setup where we can observe the generation of electromagnetic fields involves the collision of two bubbles.

Accordingly, we nucleate two bubbles along the $z$ axis at $T=0$ and let them expand and collide. The initial radius $r_{0}$ of each bubble is fixed, $r_{0}=36\Delta x$, but the initial orientations of the Higgs field inside the bubbles are random. As shown in Fig. 3, the two bubbles initially expand freely before they collide. Once the collision occurs, the process of magnetic field generation can start. Specifically, as we can see from the bottom panel in Fig. 3, the magnetic field is generated at the intersection of the two bubbles. For this two-bubble configuration, a ring-shaped magnetic field will be produced, at least initially. Let us also emphasize that, since we are using periodic boundary conditions, the two bubbles actually collide twice along the $z$ axis during the expansion as can be seen in Fig. 3.

To better understand the general features of the magnetic field resulting from multiple bubble collisions, we start with a constrained simulation in which we control the initial separation of the colliding bubbles. To this effect, we consider a regular array of bubbles with centers separated by a fixed distance $r_{s}$, which are all simultaneously nucleated at $T=0$. No additional bubbles are nucleated at $T>0$. In this case, the percolation size $r_{p}$ of the bubbles can be estimated by $r_{p}\approx r_{s}$. We consider different values for the initial separation, $r_{s}/\Delta x=80,96,112,128$, and two lattice sizes, $N=256,400$. The magnetic field is calculated using Eq. (22), and then Fourier transformed to obtain the power spectrum. We show the time evolution, in units of $m_{H}t$, of the mean wave number, $k_{\text{mean}}\Delta x/2\pi$, in Fig. 4, and the peak of the spectrum, $k_{\text{peak}}\Delta x/2\pi$, is displayed in Fig. 5. To smooth out abrupt jumps of the peak location, a ten-step moving average is taken for the latter. In both cases, we adjust the starting time in the plot to coincide with the instant when the bubbles first collide, allowing for a meaningful comparison between runs with different parameters. Assuming that the bubble is expanding at roughly the speed of light, the starting time of bubble collision can be estimated as $t=\left(r_{s}-2r_{0}\right)/2$.

From Fig. 5, we observe that the peak of the spectrum is located around $k\Delta x/2\pi\approx 0.02$ independently of the initial bubble separation $r_{s}$ and of the lattice size $N$. This value corresponds to a wavelength of $\lambda_{k}\approx 50\Delta x$, or $m_{H}\lambda_{k}\approx 9$. On the other hand, the mean wave number of the spectrum does not show a clear dependence on either the initial bubble separation $r_{s}$ or the lattice size $N$. Nevertheless, curves with larger $r_{s}$ in Fig. 4 show relatively larger fluctuations. This is because for a given lattice size, larger values of $r_{s}$ lead to a fewer number of bubbles. For instance, for $N=256$ and $r_{s}=128$ there are only 4 bubbles in the lattice, which might result in significant statistical fluctuations.

V Random bubble nucleation

Having gained an intuition for the magnetic field resulting from multiple bubble collisions, we are now in a position to model accurately a first-order EWPT by allowing for bubbles to nucleate at random locations in the unbroken phase. We perform 30 simulations on a lattice of size $N=256$, varying the nucleation probability over 6 different values, $p_{B}=10^{-3},10^{-4},10^{-5},10^{-6},10^{-7},10^{-8}$, and considering also 5 different values for the Higgs damping, $\gamma=0,2.34\times 10^{-5},10^{-4},10^{-3},10^{-2}$, where the second one is equal to the experimentally measured decay width. As can be seen from our results (Fig. 6), for $p_{B}\gtrsim 10^{-4}$ the total number of bubbles reaches the maximum number of bubbles that the lattice can accommodate almost instantaneously and the percolation size does not depend on $p_{B}$. Conversely, for $p_{B}\lesssim 1/2N^{4}\approx 10^{-10}$, the nucleation probability is so small that only one bubble may be generated before the phase transition is completed in our lattice. In practice, already for $p_{B}=10^{-8}$ the results show a large variance because our lattice contains too few bubbles. To explore configurations with smaller $p_{B}$ would require a larger lattice. In the rest of the paper, we shall focus on configurations with $10^{-7}\leq p_{B}\leq 10^{-4}$.

To better understand the results of the simulations, we divide the evolution into three stages:

Free expansion (FE stage): the time before the nucleated bubbles collide with each other. During this stage no magnetic field is generated. 2. 2.

Bubble collision (BC stage): this stage starts when the bubbles collide with each other, and the generation of magnetic field starts. At this point, the broken phase does not yet extend to the whole lattice. 3. 3.

Higgs oscillation (HO stage): the Higgs field has completely left the symmetric phase but is still oscillating, and the generation of magnetic fields continues.

Although we cannot draw clear lines between the three stages, the shift from Stage 1 (FE stage) to Stage 2 (BC stage) is signaled by the onset of magnetic field generation. This typically occurs when $m_{H}t\apprle 20$ in the simulations considered here, although, it can generically depend on the value of $p_{B}$. The boundary between Stage 2 (BC stage) and Stage 3 (HO stage) is roughly given by $T_{\text{stop}}$, which determines when the Higgs field at each lattice site has left the symmetric phase. This transition occurs when $m_{H}t\sim 200$ with some slight dependence on $p_{B}$ for the range of values that we consider.

We have carried out a set of “bubble-collision stage” simulations that focus on the magnetic field generation during the phase transition. These simulations cover the first two stages, and the magnetic fields generated up to $T_{\text{stop}}$ are analyzed. In addition, we have also performed “Higgs-oscillation stage” simulations that are run until well into Stage 3, when the Higgs has completely left the symmetric phase but is still oscillating and producing magnetic fields. We discuss them in turn.

V.1 Bubble-collision stage simulation

In this section the results of the 30 parameter combinations defined at the beginning of Sec.V will be compared during Stages 1 and 2 of the phase transition, i.e. for $T\leq T_{\text{stop}}$. As mentioned above, the stopping point $T_{\text{stop}}$ of the simulations is chosen to be when $|\Phi|_{\text{MA10}}^{2}=0.25\eta^{2}$.

In Fig. 7, various contributions to the energy density, as well as the energy density of the generated magnetic field, are shown as a function of time. Two different parameter combinations, $\gamma=2.34\times 10^{-5}$ & $p_{B}=10^{-4}$ and $\gamma=0.01$ & $p_{B}=10^{-6}$, are displayed. The total energy is not conserved due to the presence of the Higgs damping term. Among the different contributions to the energy density, it is the Higgs potential energy the one that shows a significant and sustained decrease. This is in agreement with the expectation that the Higgs damping term is working properly only on the Higgs radial degree of freedom. The simulations start with vanishing gauge fields. However, once bubbles are nucleated and start to collide with each other, energy is transferred to the gauge fields, and magnetic fields are generated as well. Although, in the period considered here, the energy in the magnetic field is only a few percent of the total energy, it does not stop increasing by the end of these runs as can be observed in Fig. 8, where we plot the magnetic energy versus time for a selection of $p_{B}$ and $\gamma$ values. Fig. 8a shows that increasing the damping reduces the magnetic field. This is so because the condition $|\Phi|_{\text{MA10}}^{2}\geq 0.25\eta^{2}$ is met earlier due to the faster attenuation of the Higgs field, and the duration of stages 1 and 2 is decreased. Also, for larger damping, a larger proportion of the total energy gets dissipated, thus reducing the energy available for magnetic field generation. This is clearly seen for $\gamma\apprge 10^{-4}$; while for $\gamma\apprle 10^{-4}$, we found the effects of Higgs damping become negligible. Fig. 8b shows that increasing $p_{B}$ increases both the magnitude and the rate $d\rho_{B}(t)/dt$ of magnetic energy density generated. Furthermore, we notice that for smaller $p_{B}$, the onset of magnetic field generation (i.e. the beginning of bubble collision) occurs later. For $p_{B}\apprge 10^{-4}$, bubble nucleation, as well as magnetic field generation, is saturated on the lattice, and the corresponding curves that fall in this range are similar to each other.

The dependence of the magnetic energy with the Higgs damping $\gamma$ is shown in Fig. 9 for different values of $p_{B}$. When the damping falls below $\gamma\apprle 10^{-4}$, magnetic field production reaches “saturation” and there is no dependence on $\gamma$. Indeed, choosing the damping equal to the SM Higgs width makes almost no difference compared to the case with no damping at all.

Also from Fig. 9, we deduce that the ratio $\rho_{B}/\rho_{\text{total}}$ of the magnetic energy to the total energy at the end of the BC stage, i.e. at $T=T_{\rm stop}$, turns out to be $\approx 2\%$-$4\%$. More specifically, for $\gamma\apprle 10^{-4}$, this ratio saturates to the value $\sim 3.3\%$, while for $\gamma\apprge 10^{-4}$ it decreases with increasing Higgs damping, approaching $\sim 2.3\%$ when $\gamma=0.01$. These results are summarized in Table 1.

A more detailed picture of the generated magnetic field can be obtained from Fig. 10, which shows the field spectrum $E_{M}(k)$, defined in Eq. (31), for various values of $\gamma$ at fixed $p_{B}=10^{-4}$. The dependence on spectrum with $p_{B}$, for fixed $\gamma=2.34\times 10^{-5}$, is displayed in Fig. 11. In both cases, the spectrum is shown as a function of the physical wave number $k\Delta x/2\pi=K^{\prime}/N$, defined in Eq. (29). The width of each bin is $\Delta k\Delta x/2\pi=\Delta K^{\prime}/N$. The spectrum is normalized in such a way that the area under the curve is equal to the magnetic energy density. One caveat is that, the nonspherical geometry of the lattice may lead to the underestimation of the spectrum for $k\Delta x/2\pi\apprge 0.5$. The vanishing tail for large $k$ shows numerical artifacts from $k\sim 1/\Delta x$ are well controlled. It is clear from Figs. 10 and 11 that certain features of the spectrum are largely independent of $p_{B}$ and $\gamma$. Specifically, the peak of the spectra lies at $k\Delta x/2\pi=0.018$. Thus, the dominant wavelength is $\lambda_{k}=2\pi/k=\Delta x/0.018\approx 56\Delta x$, or $m_{H}\lambda_{k}\approx 56m_{H}\Delta x\approx 10$ ($\lambda_{k}\approx 10^{-15}$ cm).

V.2 Higgs-oscillation stage simulations

At HO stage, the Higgs field has left the symmetric phase but it is still oscillating around the minimum of the potential. As a result, magnetic fields can continue to be generated and this is what we set out to study in this section.

To this effect we select a representative subset of the configurations considered before with the following set of parameters:

• (i)

$\gamma=2.34\times 10^{-5}$, $p_{B}=10^{-4}$

• (ii)

$\gamma=2.34\times 10^{-5}$, $p_{B}=10^{-6}$

• (iii)

$\gamma=2.34\times 10^{-5}$, $p_{B}=10^{-7}$

• (iv)

$\gamma=1.00\times 10^{-2}$, $p_{B}=10^{-6}$.

The evolution is followed for 100,000 time steps, $m_{H}t\approx 4500$, and a snapshot of the configuration is saved every 200 time steps. The outcome of this calculation is shown in Fig. 12. The upper-left plot displays $k_{\text{mean}}(t)\Delta x/2\pi$. The upper-right plot shows $k_{\text{peak}}(t)\Delta x/2\pi$, the mode of the spectrum, where, as before, the ten-step average is used. Finally, the lower plot depicts the magnetic energy density, $\rho_{B}(t)/m_{H}^{4}$.

First, we notice that the magnetic energy density keeps increasing within the time range of the simulations.

For example, for the configuration with $p_{B}=10^{-4}$ and $\gamma=2.34\times 10^{-5}$ (red curve), before $m_{H}t\sim 1000$, the magnetic energy density $\rho_{B}$ grows roughly linearly with time; although the rate slows down after that, $\rho_{B}$ keeps on increasing throughout the simulation time range. Keeping in mind that $m_{H}t_{\text{stop}}\sim 150$ for this configuration, it is clear that magnetic energy is indeed generated in the HO stage. In fact, the magnetic energy density at $m_{H}t\sim 150$ is $\rho_{B}/m_{H}^{4}\sim 0.016$, and it grows by a factor of $\sim 4$ to reach $\rho_{B}/m_{H}^{4}\sim 0.062$ at the end of our simulation. Hence, the generation of magnetic fields in the HO stage dominates over that in the BC stage.

The plot of $k_{\rm mean}$ and $k_{\rm peak}$ in Fig. 12 shows that the magnetic energy has power on length scales that are much larger than the particle physics scale $m_{H}^{-1}\approx 6\Delta x$. For example, when $\gamma=2.34\times 10^{-5}$, independent of $p_{B}$, $k_{\text{mean}}\Delta x/2\pi$ converges to $\sim 0.04$, equivalent to a wavelength of $m_{H}\lambda_{k}\approx 4.2$. The power spectrum of the magnetic field peaks at even larger length scales. From the plot of $k_{\rm peak}$ we see that the peak moves to larger length scales with time and at the end of our run, $k_{\text{peak}}\Delta x/2\pi\approx 0.011$ for all parameters. (The plot is jagged because of binning effects.) This corresponds to a wavelength of $\lambda_{k}=2\pi/k=\Delta x/0.011\approx 91\Delta x$ ($m_{H}\lambda_{k}=15.2$).

In Fig. 13 we show the energy spectrum of the magnetic fields at the end of our simulation for $\gamma=2.34\times 10^{-5}$ and $p_{B}=10^{-6}$. A peak is clearly seen in Fig. 13 and its location is largely independent of the parameters we varied in this paper. We conducted several runs on large lattices to test if the peak is due to finite lattice size and always found the peak indicating the same wavelength, independent of the lattice size. Further study is needed to determine what parameters control the location and height of this peak.

VI Discussion and Conclusion

We have simulated the classical dynamics of the bosonic electroweak theory to study the generation of magnetic fields during the EWPT, assuming that physics beyond the SM yields a first-order transition. Bubbles with the Higgs in the broken phase are randomly nucleated in regions where the Higgs is still in the symmetric phase, with the nucleation rate controlled by a parameter $p_{B}$ that we vary over a range. To account for the energy damped into fermions due to Higgs decays that will be present in the full theory we add a damping term to the scalar EOM that acts on the magnitude, $|\Phi|$, in a gauge invariant way.

We found it useful to divide the phase transition into three stages: free expansion stage (FE), bubble collision stage (BC), and Higgs oscillation stage (HO). During the FE stage broken phase bubbles are nucleated and expand in the symmetric phase, but bubble collisions have not started yet and no magnetic fields are generated. Once bubbles start crossing, the energy density in magnetic fields grows to $\sim 3\%$, or $B\sim 10^{23}$ gauss, for $\gamma$ equal to the observed Higgs decay width. Assuming the fields are nonhelical, the field strength will evolve to $\sim 10^{-15}$ gauss at present times Banerjee and Jedamzik (2004). At this point the spectrum of the magnetic field has a peak at $k\Delta x/2\pi=0.018$, or equivalently, the peak wavelength is $m_{H}\lambda_{k}=10$ (or $\lambda_{k}\approx 2$ cm in comoving scale). This is larger than the initial size of the nucleated bubbles ($m_{H}r_{0}=6.5$) but not by much. We do not see a clear dependence of the peak location on either $\gamma$ or $p_{B}$. These findings are consistent with the results obtained in similar simulations (e.g. Diaz-Gil et al. (2008b); Stevens et al. (2012)).

Similar behavior was found in Diaz-Gil et al. (2008b, a), where the generation of magnetic fields in low-scale electroweak hybrid inflation was studied on the lattice. In this context, electroweak symmetry breaking also occurs via the nucleation and growth of Higgs bubbles and the system eventually enters a regime where magnetic fields with energy density $\rho_{B}/\rho_{\text{total}}\sim 0.01$ were found. Furthermore, as in our scenario (see Fig. 8), $\rho_{B}$ was also observed to grow linearly with time. Let us emphasize however, that we consider initial conditions appropriate for a first-order phase transition with random bubble nucleation. Unlike the scenario in Diaz-Gil et al. (2008b, a), our magnetic field is initially zero and is entirely dynamically generated.

A detailed characterization of the magnetic helicity is left for future work. Nevertheless, let us point out that our equations of motion do not include an explicit $CP$-violating term so we expect the average helicity to vanish. Nevertheless, as observed in Diaz-Gil et al. (2008b, a), the dispersion is expected to be nonzero leading to nonvanishing helical magnetic susceptibility.

After the BC stage, the evolution enters the HO stage in which the Higgs oscillates around its true vacuum for a very long time. Magnetic energy is seen to continuously increase during the HO stage, even at the end of our simulation runs. Due to limitations in computation power, we are not able to see the asymptotic value of magnetic energy density from our simulations. However, it is clear that most of the final energy in the magnetic field is produced during the HO stage, exceeding that produced in the BC stage by a factor of $\sim 4$ (for typical values of $\gamma$ and $p_{B}$) and still growing at the end of our runs. When decomposing the magnetic energy in Fourier space, we find the spectrum to peak at $k_{\text{peak}}\Delta x/2\pi=0.011$, or $\lambda_{k}=91\Delta x$, comparable to the bubble size at percolation which, is approximately 0.33 of the lattice size (see Fig. 6). Our simulations suggest that the peak location is not sensitive to the damping $\gamma$ nor the bubble nucleation probability $p_{B}$. This is consistent with Fig. 6 in which we see that the percolation size is not sensitive to these parameters, in the range that we have considered.

To summarize, using numerical simulations we find that a first-order EWPT generates a significant amount of magnetic fields. While magnetic field generation has not stopped by the end of our simulations, we find that $\sim 10\%$ of the electroweak false vacuum energy is converted into magnetic fields. The energy spectrum of the magnetic field has a peak that shifts towards larger length scales. By the end of the BC stage, the peak wavelength is of the order of the bubble percolation scale, and it shifts to a longer wavelength in the HO stage.

Acknowledgements.

F.F. thanks the Institut de Física d’Altes Energies and the Universitat Autònoma de Barcelona for their hospitality. T.V. thanks the Institute for Advanced Study for their hospitality. F.F. and Y.Z. were supported in part by the U.S. Department of Energy, Office of High Energy Physics, under Award Number No. DE-SC0017987 at Washington University. T.V. is supported by the U.S. Department of Energy, Office of High Energy Physics, under Award No. DE-SC0018330 at Arizona State University.