Dynamical simulations of electroweak baryogenesis with fermions

We perform real-time numerical lattice simulations of a one-family version of the Standard Model. We model the quantum fermions using the ensemble method and treat the bosonic scalar and non-abelian gauge fields classically. Our main interest is electroweak baryogenesis, and we test the approach by considering Standard Model baryon number violation through the chiral anomaly, a truly quantum phenomenon. We find that the method correctly reproduces the anomaly, and perform the first dynamical simulations of electroweak baryon number violation including fermions.


Introduction
Since the conditions for generating a matter-antimatter asymmetry were laid out by Sakharov [1], many scenarios for baryogenesis have been proposed, based on high-energy processes in the early Universe (for a review, see for instance [2]). Here we report on progress of a technique that can be used to examine various models of electroweak-scale baryogenesis [3][4][5], whose common feature is that the source of baryon number violation is the electroweak anomaly [6].
All the key symmetry violations (CP, C and baryon number) occur in the Standard Model, and the requirement of departure from thermal equilibrium is determined by early-Universe dynamics and the electroweak symmetry breaking transition. In the Standard Model (SM), CP violation is provided by the CKM matrix of quark mixings [7,8] and possibly by the PMNS matrix of neutrino mixings [9,10], although this is less well measured experimentally. Beyond the SM, CP-violation can be provided for instance through Majorana masses, multiple Higgs field with complex coupling and/or higher-dimensional effective terms.
In electroweak baryogenesis the baryon asymmetry is expected to occur at the electroweak scale of ∼100GeV, at which time the Hubble expansion rate is around 10 −5 eV -very slow by the standards of microscopic electroweak processes. Given that a thermal electroweak phase transition in the Standard Model cannot be first order for current Higgs-mass bounds [11,12] , we are left with the problem of finding dynamics that can cause sufficient departure from thermal equilibrium. Such a possibility is raised by models inspired by preheating at the end of inflation, where the phase transition is non-thermal [13][14][15][16][17][18][19][20], or in extensions of the Standard Model with a more complicated scalar sector, where the phase transition can be strongly first order [21].
Whilst it is nowadays a relatively straightforward task to simulate the classical dynamics of the bosonic degrees of freedom of the Standard Model and minimal extensions, using techniques from lattice gauge theory, the fermions are more problematic. And given that all the CP violation of the Standard Model is in the fermion sector, and the baryon number itself is a fermion quantum number these cannot merely be ignored. One approach is to "integrate out" the fermions, and treat their dynamics as new terms in an effective action. This approach has been taken in a number of studies [22][23][24][25] (see also [26][27][28]).
Here we present a numerical method that evolves the fermions along with the bosons and, when extended to three families of fermions, can include CP violation from the fermion sector directly.
The technique that we use to simulate the dynamical fermions was first laid out in [29], with a sample study of the effect of fermions on the evolution of oscillons in 2+1 dimensions. The idea is to model the quantum averages |...| by ensemble averages ... e , and so by simulating an ensemble of fermion realizations, one is effectively sampling the evolution of the whole set of fermion mode functions. The utility of this is that the number of mode functions varies as n 3 x in three-dimensions (here n x is the number of lattice sites in one direction), which is prohibitively expensive numerically, whereas the number of modes required to give a reasonable sample is expected to vary as n x [29] (see also the 3+1D simulations in a scalar-fermion theory in [30]).
One check that the method is viable for reproducing the quantum anomaly has already been performed in an axial Yang-Mills-Higgs-fermion system in 1+1 dimensions [31], where the baryon-number violating processes where directly simulated. The results of the ensemble method were in complete agreement with the evolution using the full set of mode functions [32,33]. While there is no numerical gain in 1+1 dimensions by using the ensemble method, the results of [31] gives confidence in the technique, and provided motivation for the present 3+1 study, which would not be practical using the full set of mode functions, certainly for lattices of the size used in [24].
In this paper we use the fermion ensemble method to simulate the electroweak quantum anomaly directly. In this first study we shall be using only a single family of Standard Model fermions 1 which is enough to establish the viability of the method by direct reproduction of the anomaly. We shall also be omitting the U(1) hypercharge and SU(3) colour of the Standard Model, as this does not affect the anomaly process.
The structure of the paper is as follows: We start in section 2 with a presentation of the continuum model that we aim to simulate and in section 3 discuss how we treat the fermions and bosons numerically. In an equilibrium environment, the key baryon number violating process is the sphaleron transition. In section 4 we test the real-time fermions in the background of handmade sphalerons. We then describe how to do full boson-fermion dynamics in the context of a fast electroweak transition in 5, where we demonstrate the approach with sample simulations. We conclude in section 6. In order to keep a flow to the paper we have included a number of the technical details in the appendices, including a list of conventions in App. A, the lattice implementation in App. B, the fermion initial conditons in App. C and the details of the handmade sphalerons in App. D

The SU(2)-Higgs model with chiral fermions in 3+1 dimensions
Our field content is one family of the Standard Model (extended to include a right-handed neutrino), but without U(1) hypercharge, SU(3) gluons and without colour degrees of freedom. Therefore we have a complex Higgs doublet φ; SU(2) gauge fields W a µ ; a lefthanded SU(2)-doublet quark field q L = (u L , d L ); two right-handed SU(2)-singlet quark fields u R , d R ; a left-handed SU(2)-doublet lepton field l L = (ν L , e L ); and two right-handed SU(2)-singlet lepton fields e R , ν R . The continuum action is then written as where the set of different components is given by The absence of a U(1) hypercharge means that one is allowed more Yukawa terms than the equivalent Standard Model action, but we simply set to zero those that would not be allowed had we included the Standard Model U(1).
The charges of the fields are determined by the covariant derivatives, which are 8) and the SU(2) field-strength is defined by As well as the conserved current following from the presence of the gauge symmetry, there is a baryon current and a lepton current coming from the global symmetries q → exp(iαγ 5 )q and l → exp(iαγ 5 )l, Classically, these currents are conserved. Quantum mechanically however, one finds that [34] (2.14) and n f is the number of fermion families, which we have equal to one for our simulation. We then find that the baryon number, N f = d 3 xj 0 (b) , is related to the Chern-Simons number, N CS = d 3 xK 0 , as follows The main aim of this paper is to confirm that this truly quantum mechanical relation holds in our numerical simulations, and to apply the method to a fast electroweak transition. We will also calculate the average Higgs field 16) and the Higgs winding number, For convenience: Partial charge conjugation transformation On a lattice it is is more convenient to redefine the Fermi fields so that their kinetic terms are just the standard Dirac form, i.e. the fermion-gauge interactions are vector-like rather than chiral. This may be achieved by the following definitions [32,33,35,36], We then find thatl where some integration by parts has been done on the first three equations, leaving us with The end result is that whereas we before had two left-handed doublets and four righthanded singlets, these are now collected into one full Dirac doublet and two Dirac singlets. The singlets only interact via the Yukawa term. It is straightforward to check that in this formulation one has (2.30)

Modelling the bosons and fermions
The continuum model above is discretized on the lattice as described in App. B. For the bosonic fields the classical Hamiltonian equations are derived, with the understanding that fermion bilinears are to be represented by quantum averages. The fermion evolution equations are also found in the usual way by variation of the action, and amount to a linear Dirac equation in a time-and space-dependent bosonic background. One issue to address are the lattice fermion doublers. For every physical fermion mode, the standard lattice discretization prescription generates a set of fifteen doublers, which count as real degrees of freedom at finite lattice spacing [37]; in the context of the anomaly these are particularly troublesome. Since they contribute with opposite signs, the doublers cancel the anomaly from the physical fermions, and set the total anomaly to zero [38] (see also Fig. 3).
Fortunately, a simple approach exists to breaking the symmetry between these doublers, namely adding a Wilson term to the lattice action (B.13). As was the case in 1+1 dimensions [31][32][33], we find that it is sufficient to add a spatial Wilson term to cancel the space-like doublers, and to not initialize the time-like doublers. For small enough timestep in the numerical evolution, the time-like doublers stay unexcited, at least for the duration of our simulations (see also section 4).
The next step is to replace quantum averages by ensemble averages, as discussed in detail in [29,31]. The process starts by noting that the continuum canonical anticommutation relations for continuum fermions are and we may expand the field operator as in terms of the creation and annihilation operators. The anti-commutation relation are then equivalent to imposing and we may calculate bi-linears such as where µ is the fermion mass; it is these kinds of quantities that we aim to reproduce using ensemble averages.
To model the quantum averages we introduce two ensembles of fermions, M(ale) and F(emale), according to the mode expansion and where the exact same random numbers ξ, ζ are used in a given Male and Female pair. By then requiring that the variables ξ and η satisfy the ensemble average relations is understood to refer to the lattice version of the Dirac delta), we may for example calculate which leads to the simple prescription (3.10) For more details of the method see [29]. The mode expansion described above constitutes the initial conditions for our fermions, which then amounts to generating sets of η k , ξ k and inserting them into eq. (3.6). In the case of a general Yukawa interaction and CKM mixing matrix, one has to take care to diagonalize the fermions first to initialize the mass eigenmodes; this is explained in more detail in App. C. For the simulations presented here, we will restrict ourselves to Yukawa couplings proportional to the identity (see also App. B)

Sphaleron transitions
Now that we have the full prescription for evolving the system of bosons and fermions, we wish to test whether the method is reliable for the anomaly dynamics driven by sphaleron processes. We are mainly interested in whether the numerical fermions react in the correct way as the SU (2) gauge field undergoes a change in winding. To that end we shall start by evolving the Higgs and gauge fields by hand, rather than through their equations of motion, so we need to know how to construct a winding event. This is the method used in 1+1 dimensions [31][32][33], and here we extend the construction to the SU(2) sphaleron following the techniques of [39] (see also [40] for a nice exposition). The detailed implementation can be found in App. D. The end result is a family of configurations starting and ending in a pure-gauge vacuum of the gauge-Higgs theory, which takes the fields through a sphaleron transition. In this way there is a relative winding number between the vacua, both in terms of winding in the Higgs field, and the Chern-Simons number of the gauge field. On this bosonic background we evolve the fermions according to their dynamical evolution equations, and test that they respond in the expected way and acquire a non-zero fermion number in accordance For the moment, we will set the Yukawa couplings to zero. We generate a doublesphaleron trajectory on a n 3 x = 32 3 lattice with m H /m W = 2, a timestep a t /a x = 0.05 a lattice spacing am H = 0.42 and a Wilson coefficient r = 1.0. Fig. 1 shows Chern-Simons number N cs in time (black) and the combination (B + L)/2 increasing the fermion ensemble. The coloured lines show the ensemble averaged anomaly for N q = 640, 1280, 2560 and our largest ensemble N q = 10240. We see that the anomaly is very well reproduced at the ten percent level for these ensemble sizes, and that this can be improved upon in a straightforward way. For a 32 3 lattice, numerical efficiency is only gained over the "non-statistical" approach for ensemble sizes less that 32 3 = 32768.
The next test is to vary the lattice spacing (for fixed number of lattice points), shown in Fig. 2. We see that the anomaly is well reproduced for all lattice spacings, although the am H = 0.105 is a little low, possibly because the physical volume is quite small n x am H = 3.36. Also the largest lattice spacing am H = 0.42 is about 5 percent low, which may in turn be put down to discretization errors. We note that the anomaly is expected to be one of the "hardest" observables to get right, as it does not involve any volume averaging (such as for instance particle number in a homogeneous background would [30]). Also, as can be deduced from the level-crossing picture of the anomaly, all momentum modes contribute to the anomaly, both in the IR and the UV. Hence the statistical averaging has to be accurate for all modes independently, to get the correct anomaly out. It is not possible to leave out some of the UV modes, as if one were only interested in the IR physics. In this light, we consider it remarkable that our simulations reproduce the anomaly so well. In Fig. 3 we show how the anomaly depends on the coefficient of the Wilson term. As mentioned, for r = 0 lattice doublers are expected to pair up and cancel out the anomaly, whereas for r = 1, the anomaly is manifest. We indeed see that the anomaly is identically zero for r = 0, but also that once r > 0.25 the anomaly is largely unaffected by r (within statistical errors). This is as we would have expected. Also note that the time-like doublers are not cancelled by a Wilson term, but that these are not excited enough to cause problems.
We then consider the case where scalars are directly coupled to the fermions, i.e. the fermions are massive. As shown in Fig. 4, at finite Yukawa coupling convergence slows down with increasing coupling. This is completely analogous to the situation in 1+1 dimensions [31], and the solution is also the same: increase the size of the ensemble. N q = 10240 is again sufficient for λ yuk < 0.1, corresponding to a fermion mass of λ yuk v/ √ 2 ≃ 17 GeV, i.e. all Standard Model quarks except the top. But to include the top, λ yuk ≃ 1, we need many more realizations. We did a test at the smaller lattice size of n 3 x = 16 3 with N q = 8 × 10240 and indeed found that the anomaly is still present, shown in Fig. 5. Such simulations suffer from having small physical volume, and the fermion number is a little low. We also explicitly checked that decreasing the timestep and/or increasing the Wilson coefficient r did not improve convergence. [39] We can conclude from these tests that the fermion ensemble method works even for very subtle quantum observables like the anomaly, and that we need N q = O(10 4 ) to get a measurement at the ten percent accuracy level, at zero and small Yukawa couplings, up to a fermion mass of a few GeV. A larger ensemble is required to consistently include heavy quarks, i.e. the top. Lattice issues like lattice spacing dependence and fermion doublers are well under control. We will now proceed to also make the bosonic fields dynamical and include the fermion back-reaction.

Including fermion backreaction: tachyonic electroweak transition
In the tests with handmade sphalerons presented above, the fermions were just "along for the ride", to see how well the fermion equations of motion behave in a time-dependent bosonic background. But for baryogenesis, we need the bosonic fields to be dynamical, too, and for the fermions to back-react. Indeed, in the SM CP-violation is a fermion back-reaction effect. As a testing ground for this, we consider a cold tachyonic electroweak transition, where the Higgs field experiences a fast quench. The Higgs potential flips instantaneously as The Higgs field is generated by Monte-Carlo sampling of a classical ensemble, reproducing the quantum correlation functions at zero temperature in the symmetric phase [18,19,23]. The fermion fields are generated as an ensemble of N q vacuum realizations at φ = 0, as above. The initial gauge fields are then found by setting A µ (t = 0) = 0, and solving the lattice Gauss law in the background of the generated Higgs and averaged fermion fields. The dynamics start at the moment when the Higgs potential is flipped, and the lowmomentum modes k < µ experience a tachyonic (or spinodal) instability, and start growing exponentially. This leads to a strongly out-of equilibrium electroweak transition, where large particle numbers, effective diffusion of Chern-Simons number and in the presence of CP-violation, baryogenesis can occur [13,20,23,24]. In the present work, we have only implemented one generation of Standard Model fermions, and so we are not able to include CP-violation via the CKM matrix. We will leave that for future work. As in [23] we will monitor the volume-averaged Higgs field, the Chern-Simons number and the Higgs winding number. But rather than infer the baryon and lepton number through the anomaly equation, we now calculate it directly. We note that the anomaly equation states that changes in the Chern-Simons and Baryon/Lepton number obey ∆B = ∆L = ∆N cs → ∆N cs = ∆B + ∆L 2 .
(5.2) Fig. 6 shows a tachyonic run where the fermions do not back-react (dashed) and one where they do (full). Yukawa couplings are zero, and N q = 5120. We see that the Higgs field (black) "falls off the hill" and oscillates around its broken phase minimum. Meanwhile, the Chern-Simons number (red) and Higgs winding number (green) bounce around out of equilibrium, but eventually settle down in a common minimum. To check our numerics, we also monitored Gauss law, which throughout is conserved to computer accuracy. As for the sphalerons above, we see that the fermions (blue) reproduce the anomaly rather well, however, the gauge and Higgs fields are now dynamical.
What is particularly interesting is the effect of having fermions affecting the bosonic dynamics. At early times, the evolution of the bosonic observables is unchanged. But then around m H t = 15, the trajectories of both fermion number and Chern-Simons number diverge, depending on whether the bosons feel the fermions, eventually ending up at different minima, and hence a different number of fermions.
Finally, in Fig. 7 we show a run with small enough λ yuk = 1/128, that N q = 10240 gives reliable convergence. This corresponds to a 1.4 GeV fermion, i.e. slightly heavier than the charm quark. As for the sphaleron runs, we see that the agreement is only approximate, but clearly the anomaly is reproduced, improveable with larger statistics.

Conclusions
While there have been many simulations of the bosonic fields in the Standard Model, this is the first that has included the fermionic sector, rather than integrating them out and/or treating them as passive observers to the gauge-scalar dynamics. This has been achieved by applying the technique pioneered by Borsányi and Hindmarsh [29], where the many mode-functions of the fermion fields, required for a full quantum treatment, are replaced by an average over a statistical ensemble of Male and Female fields.
We have performed a detailed study of the lattice parameters, and found that the fermion number does indeed reproduce the quantum anomaly, with the fermion number tracking the Chern-Simons number through sphaleron transitions. We also simulated a fully dynamical fast electroweak transition including fermions, and found that the backreaction of fermions can indeed significantly alter the evolution of the gauge fields.
From a technical point of view, the numerical effort is significant, in that implementing N q fermion realizations on a n 3 x lattice requires n 3 x × N q × 1.6 kB, so that for our single generation, no-colour full run of Fig. 6 we use approximately 270 GB. It is however important to note that the anomaly is probably one of the most sensitive observables, since it requires every single lattice mode to be very precisely reproduced. This means that for studying the (IR) dynamics of electroweak baryogenesis, we can probably make do with a somewhat smaller N q . Alternatively, one may consider returning to the non-ensemble approach [32,33], where one is however hampered by the need for large lattices to correctly simulating tachyonic transitions.
The next step is now to study dynamics both for baryogenesis and electroweak preheating; including all three generations, mixing and CP-violation for both the quark (CKM) and lepton (PMNS) sector as relevant for the most minimal version of Cold Electroweak Baryogenesis; but also for electroweak baryogenesis in more general setups like Standard Model + second Higgs doublet or Standard Model + scalar singlet. Demonstrating that this is a realistic aim (also in practice) is perhaps the most significant result of the present work.
with in the Weyl representation The charge-conjugate scalar doublet is given by the complex conjugate

B Lattice action
The theory is discretized on a lattice with spatial lattice spacing a x , timelike spacing a t = dta x . We find it convenient to work with the following rescaled and re-defined variables thus making them dimensionless. Then we have the link variables, and corresponding electric fields We need the derivatives, where we will not include the lattice spacings (they are compiled into lattice coefficients, see below) Note that we are using the compact formulation of the gauge action, and have chosen to evolve using temporal gauge W a 0 = 0 (U 0 = I). In order to remove the spatial fermion doubler modes we will employ Wilson fermions in space, for which we need the Wilson term [32,33,41] where r is a constant parameter. We will also introduce the quantities The action is then given by To which we add the Yukawa terms, now written in terms of the matrix Φ, and where we for simplicity impose that (B.18) and in addition set Then we have, using the addition to the action The equations of motion follow by variation. For the Higgs field where δΦ = δφ 0 I + δφ a iσ a , (B.23) and are ensemble average over fermion realizations. For the gauge field we have Tr iσ a Φ y (U y,n Φ y+n ) † − a t aβ t g Ψ y γ n iσ a U y,n Ψ y+n +Ψ y+n γ n U † y,n iσ a Ψ y and we have Gauss law Finally, the fermion fields evolve according to The currents that are particularly relevant to us are the baryon and lepton currents which, when converted to male and female fermions, becomes

C Fermion initial conditions
In order to specify the initial conditions we need some mode functions for the fermion fields. The ones we use are those that are vacuum fluctuations for a Higgs field taking the value φ = (0, v/ √ 2),φ = (v/ √ 2, 0). The lattice equations of motion then read with M andM some general Yukawa mixing matrices. This is a set of linear equations, so we perform a Fourier decomposition and look at positive frequency modes and find that by defining now definem = m(k)/a andṽ = v/ √ 2 and make the simplification that the Yukawa couplings are related by so that defining allows us to diagonalize the equations by introducing the canonically normalized Fermi fields Ω i , , (C.10) , (C.11) , (C.12) . (C.13) The positive frequency modes of these Fermi fields then satisfy − iγ µ s µ U 1 = (m + |g u | 2 )U 1 , (C.14) −iγ µ s µ U 2 = (m − |g u | 2 )U 2 ,

D Handmade Sphalerons
The basic idea [39] is that the vacuum manifold of the Higgs scalar is a three-sphere, S 3 vac , and one wants to construct some other three-sphere in order to have non-trivial maps between them. At spatial infinity the Higgs field is in the vacuum, but spatial infinity is only a two-sphere, S 2 ∞ , and the maps S 2 ∞ → S 3 vac are topologically trivial. However we can construct a three-sphere by considering a loop in configuration space, S 1 loop , and forming the smash product S 2 ∞ ∧ S 1 loop , which i s a three-sphere, and so one may construct nontrivial maps. This loop in configuration space means we have a family of configurations parametrized by an angle.
Of course, we also need an ansatz for the gauge field, and this is achieved by constructing the unitary matrix The gauge field is in the vacuum at infinity, and so we take the asymptotic ansatz to be with the full radial dependence of the family of gauge fields being given by W j (Γ, r, θ, ϕ) = f (r)W j(∞) (Γ, r, θ, ϕ), (D. 6) where f (r) is a smooth monotonic function satisfying f (0) = 0, f (∞) = 1.