Ensemble fermions for electroweak dynamics and the fermion preheating temperature

We refine the implementation of ensemble fermions for the electroweak sector of the Standard Model, introduced previously. We consider the behavior of different observables as the size of the ensemble is increased and show that the dynamics converges for ensemble sizes small enough that simulations of the entire electroweak sector become numerically tractable. We apply the method to the computation of the effective preheating temperature during a fast electroweak transition, relevant for Cold Electroweak Baryogenesis. We find that this temperature is never below 20 GeV, and this in combination with the early results convincingly rules out Standard Model CP-violation as the origin of the baryon asymmetry of the Universe.


Introduction
The non-equilibrium dynamics of the electroweak sector of the Standard Model and its extensions is crucial for the understanding of baryogenesis and leptogenesis in the early Universe. A large body of work exists on equilibrium quantities including the sphaleron rate (see [3] for recent results) and electroweak phase diagram [4], based on a dimensionally reduced version of the theory, and out-of-equilibrium dynamics has been studied using the classical approximation for the bosonic degrees of freedom (see for instance [5][6][7][8][9][10][11][12]).
However, a complete understanding of, in particular, fermion production and the baryon asymmetry requires us to include fermionic degrees of freedom, and these are inherently quantum mechanical. It has been known for some time how to combined quantum fermions with classical bosonic fields out of equilibrium [13][14][15][16][17], using either the complete set of quantum modes, or a statistical ensemble approach. Recently this method was implemented on a lattice for the reduced Standard Model; including only SU(2) gauge fields, a Higgs field and a single family of mass-degenerate quarks and leptons [1,18].
It was demonstrated that the chiral anomaly is indeed generated in this formalism, and that a sizable number (tens of thousands) of realizations are required for the ensemble fermion number to converge, in particular for large values of the fermions-scalar Yukawa coupling. In the Standard Model, the top quark has a coupling of approximately 1, requiring another factor of ten in ensemble size. Although it was highly encouraging that the anomaly is reproduced, the massive numerical effort threatened to make the method unmanageable for realistic systems.
In this work, we point out that the slow convergence of the fermion number is irrelevant to the dynamics of the fields, but is a feature of the observable. We refine the analysis and show that all other observables of interest converge at much smaller ensemble size, and in particular that the back reaction of fermions onto the bosonic fields converges faster. In practice, this means that simulations are reliable with of order 2000-3000 fermion realisations, which is certainly numerically tractable.
We also introduce a way to compute the particle spectrum and extract an effective temperature from the fermions, using as a testing ground a fast electroweak quench transition. As a direct application of this, we calculate this temperature as a function of the quench rate. This number enters in recent computations of effective CP-violation in Cold Electroweak Baryogenesis [2,9,[19][20][21][22].
The structure of the paper is as follows: In section 2 we set up the reduced Standard Model and in section 3 the ensemble fermion method. In section 4 we investigate the behavior and convergence properties of global variables such as the Chern-Simons number, Higgs winding number, fermion number and the energy. Section 5 discusses the measurement of the particle number and the effective temperature, and in section 6 we perform the full simulations to find the quench time dependence, further averaging over an ensemble of bosonic fields. We conclude in section 7.

The reduced Standard Model
We will consider a simplified Standard Model including the SU(2) gauge field W a µ coupled to the Higgs doublet φ and one family of fermions; a left-handed quark doublet q L = (u L , d L ), and two right handed singlets u R , d R , and similarly a left-handed lepton doublet l L = (ν L , e L ) and two right handed singlets e R , ν R , including a right-handed neutrino. As a consequence, we ignore hypercharge and gluonic gauge fields as well as the second and third family of fermions. The continuum action is then written as with the components 3) The covariant derivatives are 8) and the SU(2) field-strength is defined by The Higgs mass is taken to be 125 GeV. It is convenient to redefine the Fermi fields so the kinetic terms have the standard Dirac form, with vector-like gauge-fermion interactions with = iσ 2 , and C = diag(− , ). It follows that leaving us with Whereas we before had two left-handed doublets and four right-handed singlets, these are now collected into one full Dirac doublet and two Dirac singlets. The latter only interact via the Yukawa term. For simplicity, we will restrict ourselves to the case which corresponds to all fermions (quarks, charged lepton, neutrino) having the same mass. The global symmetry q → exp(iα)q and l → exp(iα)l implies classical conservation of the currents 19) and in terms of the redefined fields, we find that (2.20) At the quantum level, the baryon and lepton currents are no longer conserved due to the chiral anomaly, and and n f is the number of fermion families, here taken to be one. The baryon number, , is related to the Chern-Simons number, N CS = d 3 xK 0 , as It was demonstrated in [1,18] that this relation is reproduced on a lattice using ensemble fermions as described in the following.

Ensemble fermions
The continuum model above is discretized on a lattice of V /a 3 = N 3 sites, as described in [1]. The dynamics of bosonic fields is assumed to be classical, and their time evolution follows from a simple variation of the action. Similarly, since the fermions are bi-linear in the action, variation with respect to the fields leads to linear fermion equations of motion involving the bosonic fields, and these equations are solved simultaneously with the bosonic field equations. We deal with the fermion doublers by including a Wilson term in the spatial directions, and by not initialising the lattice time-like doublers. With a small enough timestep, it then takes a long time before these get excited, longer than the duration of our simulations (see also Appendix A).
In the bosonic equations of motion, fermions enter through bilinear quantum correlators. These can be computed by expanding the fermions on momentum modes in terms of time-independent creation/annihilation operators [13,14], these operators in turn encode the initial condition. Alternatively one may replace quantum averages by ensemble averages, as discussed in [15,18]. Since fermion observables can be extracted from the time-ordered Green function, and the fermion back reaction on the bosons is through the equal-time correlation functions, it will be sufficient to give the fermion back reaction, currents and the energy, and demonstrate that these local two-point functions can be well represented by the method.
The equal-time correlation functions for the fermion can be written as, and we notice that We may expand the field operator as in terms of the annihilation and creation operators b and d † . The fermion anti-commutation relations correspond to so the equal-time correlation function can be written out explicitly in the vacuum We introduce two ensembles of fermions, M(ale) and F(emale), and where the exact same random numbers ξ, ζ are used in a given Male and Female pair. We then require that the variables ξ and η satisfy the ensemble average relations 1 , we may calculate Generating sets of η k , ξ k and inserting them into eq. (3.7) provides the initial condition for the fermion ensemble fields, and these are solved in position (x) space and averaged over to generate the bilinears at each time-step, which are in turn fed into the bosonic equations of motion. The number of realizations in the fermion ensemble is denoted N q , and we must ensure convergence of the physical observables as N q is increased. In [1], it was found that for small values of λ yuk , N q 10000 gave convergence of the fermion number observable N f .
The fermion evolution is unitary, and will conserve inner products for each gender, Male or Female where the field is discretized and normalized by the lattice spacing (see [1]). There exist local forms for arbitrary x and p. We will use them to monitor the temporal doubler (see also Appendix A).

Global observables
We will compute a number of different observables to check the convergence of the simulation and determine N qc , the smallest N q for which the observables can be said to have converged. N qc should be robust to different initial conditions that bring different dynamics.
We define the average Higgs field squared, It is scaled to be unity when the Higgs field is in the broken phase vacuum. The total energy is conserved up to the order O(a t /a), with energy density components defined from the lagrangian The contribution e C includes the counterterms needed to formally keep the theory finite (see Appendix A). The physical fermion energy is obtained by summing over e F , e Y , e W and e C , and normalized by subtracting its initial value. Then around a well-defined vacuum, the fermion energy can be interpreted as the total energy of all excited particles. We test the dynamics, convergence and observables in a fast electroweak transition, where the Higgs potential undergoes an instantaneous quench +µ 2 → −µ 2 . In Fig. 1, we show the convergence of a set of global observables (Higgs field, N CS , N W and energy). We make a comparison between N q = 9600, 2400 and 1440, with lattice size N = 32, Yukawa coupling λ yuk = 0.03 and the lattice spacing am H = 0.42. We see that convergence is achieved for these observables with N q 1000 − 2000 and we take as a conservative choice N qc = 2400. We found that this was safe also for larger lattices of N = 48. Convergence is less good for larger Yukawa couplings, but improve for smaller couplings. We note that λ yuk = 0.03 corresponds to a fermion mass of 5.1 GeV, so that all Standard Model fermions except the top correspond to smaller values of the coupling.
In Fig. 2 we show one special observable, the baryon number (or fermion number), which has not converged yet, even at N q = 4800. This is consistent with the findings of but we now see that it is a result of the observable being badly behaved, not the dynamics. In particular, from Fig. 1 we can conclude that the fermion back reaction on the bosonic field has also converged at N qc , and as a measure of the baryon production, we can simply use bosonic operators, N CS and N W . This is with the understanding that had we increased N q for a factor of 10 or more, N f would converge to this same number [1]. A further example is to look at the individual energy components, shown in Fig. 3 for two different values of the Wilson coefficient r w . When r w = 0 (left plot), the energy transfer to the fermions is much faster than when r w = 0.5 (right plot). Apparently the Wilson term is making the doubler modes more heavy, so that they can no longer be excited. And so even when not looking specifically for the effect of the baryon anomaly, it may be worthwhile including the Wilson term in the dynamics. In our simulations, the Wilson parameter r w is fixed to 0.5.
The lattice spacing dependence is shown in Fig. 4, showing the different energy components on two different lattices with the same physical volume, but different lattice spacing (N = 32, am H = 0.63 and N = 48, am H = 0.42). We see that by choosing the counterterms carefully, the lattice spacing dependence is nicely under control. We note that the random bosonic initial condition is different in the two cases (since there are a different  number of modes present), and so the agreement is quite non-trivial.

Spectrum
The fermion particle number can be extracted from the two-point correlation functions. If the fermion field is close to thermal equilibrium, fitting of the particle number with the Fermi-Dirac distribution will give the (effective) temperature and chemical potential of the system. We define the correlation function in momentum space through a Wigner transform and averaging over the space volume, which is equivalent to computing |Tψ(p, t)ψ(p, t)| . For future use, we define: A gauge-invariant correlation function for the gauge doublet can be defined as where U (y, x; t) is the gauge link connecting x and y at the time t. For the gauge singlet field, this is not needed D ξ,χ (x, y; t) = 0|Tξ(x, t)ξ(y, t)|0 , 0|Tχ(x, t)χ(y, t)|0 .
which is obviously true for the setup (3.3) and (3.5). By assuming the on-shell condition ω(p, t) = p 2 + m 2 (p, t), we have enough freedom to measure the effective energy ω(p, t) and the average particle number N av (p, t) ≡ (N (p, t) +N (−p, t))/2 simultaneously. In the unitary gauge, the mass eigenstates Ω 1,2 can be written in terms of the weak eigenstates as where Ψ u , Ψ d are the up and down parts of the gauge doublet.
Having computed the particle number N av (ω), we can proceed to extract the effective temperature and chemical potential, by fitting to the form If the system is properly equilibrated, the low-momentum range should indeed show this form (Fig. 5).
To estimate the uncertainty from the gauge non-invariance, we have compared results from different choices of correlators in Fig. 6 (left). The singlet (field ξ) is simply the correlator (5.4) and doublet (field Ψ u ) is (5.3), but omitting the link variable. We compare these to the temperature of an mass eigenstate (Ω 1 ) in the unitary gauge (5.8). We see that the doublet temperature deviates from the other two, with the unitary gauge one the smoothest. The right-hand plot shows the effective chemical potential, which only at later times becomes positive.
In Fig. 7 we show a comparison between difference ensemble sizes N q demonstrating that our choice of N qc = 2400 also holds for these derived observables. We can also compare results for different values of the Yukawa coupling, Fig. 8. The values λ yuk = 0.1, 0.03 and 0 correspond to fermion masses ∼ 17, 5.1, and 0 GeV respectively, and we see that the temperature is largely independent of the choice of mass. This suggests that most of the energy transfer comes down through the gauge field coupling, rather than directly through the coupling to the Higgs. This conclusion may not hold for the top quark mass, which is high above the effective temperature. The dependence on the Wilson term coefficient was found to be comparable (not shown).

Fermion temperature in Cold Electroweak Baryogenesis
Cold Electroweak Baryogenesis assumes that the Universe came out of inflation around the electroweak scale, triggering a fast electroweak quench. The out-of-equilibrium dynamics of this tachyonic transition is then responsible for Baryogenesis, without a need for a first order phase transition with bubble nucleation [9,[20][21][22]. It is known that the observed asymmetry can be generated in the presence of a generic additional source of CP-violation [10,11]. But extensive work has also gone into the possibility that the CP-violation already present in the Standard Model through the complex phase in the CKM matrix may be sufficient [2,19,[23][24][25].
The status is that a number of effective bosonic operators arise upon integrating out the fermions at finite temperature, and the coefficient functions are strongly suppressed with temperature. Bosonic simulations suggest, that the coefficient should be of order 10 −6 in some normalization [11], corresponding to an effective temperature in the fermions of around 1 GeV. Although an equilibrium computation does not necessarily represent the out-of-equilibrium state during the tachyonic transition, it certainly gives the best estimate currently available. With the techniques outlined above, we are now in a situation to compute this effective temperature. It is also known that generating the asymmetry requires the quench transition to be fairly fast [12]. This is in order for the dynamics to be violent enough that non-perturbative effects such as baryon number violation can take place. Finally, it is known that the asymmetry is created during the first or second period of the Higgs field oscillation.
We now introduce a quench time as in [12] to parametrize the flip of the mass parameters in the Higgs potential We then simulate the transition for different values of τ Q , in each case computing the effective temperature averaged during the first and second period of the Higgs field oscillation. We use N q = 2400, N = 48, l yuk = 0.03, am H = 0.63 and we in addition average over 8 realizations of the bosonic fields. An example of such an averaged effective temperatures and chemical potential is shown in Fig. 9, at a quench time of m H τ Q = 30. 1st Higgs minimum 2nd Higgs minimum the late time limit Figure 10: The first and second Higgs minimum correspond to the first and second time the Higgs field rolling back to the minumum of the average Higgs field squared. Each point is the statistical result of 8 runs, with the error bar stands for σ confidential interval. The statistical error for the late time limit is smaller than 0.2 GeV.
We finally show the effective temperatures at the first and second minimum of the Higgs field, as well as the late time limit, as a function of quench time (Fig. 10). The effective temperature oscillates in time, and the results are therefore sensitive to the exact time assigned to the Higgs minima. But the picture is clear: the temperature decreases in time, and it also decreases with quench time. For each of the measurement times (first, second minimum, late time), we can fit the temperature by a straight line. extrapolating these to large quench times, we find that the required 1 GeV may be reached for m H τ Q 150. This is an order of magnitude slower than the requirements for a violent out-of-equilibrium phase during the transition [12]. We conclude that Standard Model CP-violation cannot be strong enough to generate the observed baryon asymmetry, also not in the context of Cold Electroweak Baryogenesis.

Conclusion
In summary, we have carefully developed the real-time lattice implementation of the electroweak sector of the (reduced) Standard Model introduced in [1]. It includes the SU(2) gauge and Higgs fields, treated classically, and one generation of fermions, treated quantum mechanically, using the ensemble fermion method. All the fermions are taken to have the same mass, around 5 GeV.
We found that the ensemble method converges for or order 2000 realisations of the fermions, for the dynamics and all observables except the fermion/baryon number. This is a reduction of the numerical effort by an order of magnitude, compared to what was anticipated in [1], making the method viable for modelling the entire electroweak sector. Including the full fermion spectrum, in particular the massive top quark will make convegence slower. Our largest lattices (48 3 , N q = 2400) fill 135GB memory, which is certainly tractable by modern supercomputers. Triple this for the full SM fermion content.
As an application of the method, we considered the tachyonic preheating mechanism after hybrid inflation, in the context of Cold Electroweak Baryogenesis. Standard Model CP-violation gives rise to effective higher order bosonic interactions, the coefficients of which are strongly temperature dependent. In order for SM CP-violation to be sufficient for generating the observed baryon asymmetry, we need a transient stage during preheating, where the fermions have a temperature in the region of 1 GeV. We find that the transient effective temperature in the IR is 10 to 20 times higher than this, lower for slower quenches. Extrapolation to very slow quenches potentially leads to a lower preheating temperature, but we come in conflict with the need for a far-from-equilibrium state also needed for successful baryogenesis.
We conclude that SM CP-violation is insufficient for baryogenesis, also in the Cold scenario. One caveat is that we have not included the top quark; however, since the top mass is about 5 times the largest temperature encountered here, only non-perturbative processes during the initial rolling down of the Higgs field can source it. We believe that this effect is unlikely to bridge the gap in temperature. Similarly, the reheating temperature scales with the number of relativistic degrees of freedom as (g * ) −1/4 , and so including all the SM degrees of freedom (a factor of about 4 larger compared to what we have here, depending on whether W and Z are taken to be relativistic) also does not reconcile the measured temperature with 1 GeV.
The approach to fermion dynamics considered here has now matured to the point where one may apply it to a number of phenomena, including magnetic field and gravitational wave creation during the electroweak phase transition, preheating dynamics in extensions of the SM as well as high temperature baryogenesis mechanisms where bubble walls sweep through a hot plasma. A first step is the implementation of the full three generations of the SM, with physical masses and mixings, and extensions of the Higgs sector, projects presently under way Acknowledgments: The numerical simulations were implemented and performed on the COSMOS supercomputer, part of the DiRAC HPC Facility which is funded by STFC and BIS. ZGM wishes to thank Shuang-Yong Zhou and Paul Tognarelli for useful discussions.

A Fermion doublers and the Wilson term
In our simulation, we include only the spatial Wilson term to reduce the effect of fermion doublers, The temporal doubler is suppressed by choosing the initial condition carefully, so that the doubler starts out un-excited. For a long enough simulation, the doubler mode will return, but we find that this happens on a timescale much longer than the simulations presented here.
Considering a single mode of the U part (with U one of the eigenspinors), which follows the difference equation where the t is discretized into integer values, and on the first step U k (1) = U k e ik.x . The eigenspinor U k is the solution of where the dispersion relation reads sin 2 k 0 = i sin 2 k i + m 2 k . So here the recurrence relation is The general solution is The second term on the right-hand side is the doubler and will switch the sign from step to step. To remove it, one needs to select the initial condition that U k (2) = e −ik 0 U k (1), and similarly, V k (2) = e ik 0 V k (1) for the V part. With initial conditions chosen carefully for each mode, the temporal doubler will not be excited in the beginning. During the simulation, we keep track of the physical and doubler parts by measuring three closest time steps, If there are only physical modes, and on the contrary, with only doubling modes, In practical simulations, A.8 is accurate up to the order of 10 −5 .

B Two-point functions on the lattice
For the Dirac field, the free two-point function has the form: The The above perturbation respects the gauge symmetry, in the sense that the Ward identity is fulfilled in the form i µ (1 − e ikµ )Γ µ (p + k, p) = −iG −1 0 (p + k) + iG −1 0 (p). (B.5)

C Counterterms
We introduce counterterms for the back reaction of the quantum fermions onto the classical bosonic fields. The full action S C is chosen to be, The first term describes the screening effect when Z 3 < 1, and will reduce the coupling constant. In the equation of the motion for the Higgs field, the cancellation is where G(x, y) and G (x, y) are Green functions for different mass eigenstates in the weak theory. For the gauge field 4(Z 3 − 1) g 2 (E a n (x) − E a n (x − 0) − 4(Z 3 − 1) g 2 m D ab m Tr iσ b U x,m U x+m,n U † x+n,m U † x,n = Ψ x γ n iσ a U x,n Ψ x+n +Ψ x+n γ n U † x,n iσ a Ψ x − r w Ψ x iσ a U x,n Ψ x+n −Ψ x+n U † x,n iσ a Ψ x = iTr[iσ a U n (x)G(x + n, x)γ n ] + iTr[U † n (x)iσ a G(x, x + n)γ n ] − ir w Tr[iσ a U n (x)G(x + n, x)] +ir w Tr[U † n (x)iσ a G(x, x + n)]. This leaves three dimensional discretised lattice sums, which are computed numerically. We choose the counterterm for the potential where ct 1 depends on the lattice spacing quadratically, and ct 2 depends on the logarithm of the lattice spacing. φ can be set to be constant in B.2 and C.2 to get the ct 1 and ct 2 quickly. To obtain Z 3 and Z φ , one may select the field to have a particular momentum to simplify the calculation. With our choice of lattice parameters, 0.98 < Z 3 < 0.99, and Z φ ∼ 1.