Real-time Fermions for Baryogenesis Simulations

We study how to numerically simulate quantum fermions out of thermal equilibrium, in the context of electroweak baryogenesis. We find that by combining the lattice implementation of Aarts and Smit [1] with the"low cost"fermions of Borsanyi and Hindmarsh [2], we are able to describe the dynamics of a classical bosonic system coupled to quantum fermions, that correctly reproduces anomalous baryon number violation. To demonstrate the method, we apply it to the 1+1 dimensional axial U(1) model, and perform simulations of a fast symmetry breaking transition. Compared to solving all the quantum mode equations as in [1], we find that this statistical approach may lead to a significant gain in computational time, when applied to 3+1 dimensional physics.


Introduction
In electroweak baryogenesis [3] the source of baryon number non-conservation is the quantum anomaly of fermions chirally coupled to the Standard Model SU(2) gauge field. When the gauge field evolves in such a way that its Chern-Simons number changes, the fermion, and hence B(aryon) and L(epton), number changes as where n f = 3 is the number of fermion generations in the Standard Model. The question of successful baryogenesis thus reduces to whether a permanent change of Chern-Simons number can take place in the early Universe, presumably under the influence of CP-violation and the back-reaction of the fermions. Various models of baryogenesis have been proposed, of which the most popular (and most developed) is "hot" electroweak baryogenesis [4], where walls of bubbles nucleated in a first order phase transition interact in a CP-violating manner with the fermions in the hot plasma. In this way a net left-right fermion asymmetry is generated inside and outside the bubbles, and equilibrium gauge dynamics (sphaleron transitions) convert this asymmetry into a baryon asymmetry.
The rate of sphaleron transitions can reliably be calculated in thermal equilibrium using sophisticated Monte-Carlo methods [5,6]. In such a setup, fermions can be included in terms of effective couplings for the bosonic theory, for instance through dimensional reduction [7].
An alternative scenario is "Cold" electroweak baryogenesis [8,9,10,11,12,13], where the electroweak phase transition does not involve bubble nucleation, but instead a fast quench of the Higgs potential. Here, baryon number violating processes are not equilibrium Sphaleron transitions, but complicated out-of-equilibrium field dynamics.
Numerical real-time simulations of electroweak baryogenesis have until now neglected dynamical fermions. Instead, purely bosonic systems are evolved and baryon number has simply been assumed to follow the gauge field Chern-Simons number in accordance with the anomaly equation, ignoring fermionic backreaction.
One case where this is certainly not allowed is for minimal electroweak baryogenesis, since CP-violation in the Standard Model originates from the fermion mass matrix. A possible approach employed in [14,8,12] is to integrate out the fermions in the path integral or in perturbation theory, thus recovering CP-violation effects in terms of a series of higher-dimensional bosonic terms.
The current understanding that Standard Model CP-violation is strongly suppressed at high temperatures, and therefore insufficient for successful baryogenesis follows from such a computation (see for instance [14,15]). In contrast, at low temperatures relevant for "Cold" baryogenesis, recent calculations have shown that the suppression is absent [16,17,18,19], and direct numerical simulations have in turn indicated that Standard Model CP-violation may in fact be large enough to accommodate the observed asymmetry [13,20].
A possible caveat to this procedure is that it is based on a gradient expansion in the gauge and Higgs fields, which may not be valid during electroweak symmetry breaking. And so although the work in [13,20] is very encouraging indeed, it would be even better not having to integrate out the fermions, but include them directly in real-time simulations of the transition. In this way, the CKM matrix and CP-violation could be included from first principles.
In [1] Aarts and Smit showed how to implement quantum fermions in real-time, coupled to classical bosonic gauge and scalar fields. The method involves a proper lattice discretization in Minkowski space, and the realisation that since fermions are bilinear in the action, the field operators can be expanded into mode functions, in terms of time-independent creation-annihilation operators. These mode functions can then be solved in the classical bosonic background, with the back-reaction on the bosonic fields defined as the quantum averages over the creation-annihilation operators for some given initial state.
In practice, the problem is that for every momentum mode k (equal to the number of lattice sites n D x , where D is the number of spatial dimensions), one needs to solve a separate real-time field equation (the mode function equation) for which the numerical effort is also proportional to the number of lattice sites. Hence the total numerical problem scales as n 2D x , and quickly becomes unmanageable for large three-dimensional lattices. Large lattices are often required in baryogenesis simulations to accommodate extended objects such as sphalerons and for having enough infrared modes for a fast quench to be correctly reproduced. Some time ago [2], Borsanyi and Hindmarsh showed how to replace the n D x mode equations by an ensemble of fermion field realisations, approximating the quantum fermion expectation values through a statistical averaging procedure. In the context of a scalar-fermion theory, they showed that one can significantly reduce the numerical effort, at least in three dimensions. This is because the number of random realisations in the ensemble N q can be much smaller than n D x . In this work, we will implement the "low cost fermion" or "fermion ensemble" method of Borsanyi and Hindmarsh to the 1+1 dimensional axial-U(1) model with fermions of Aarts and Smit. This will act as a toy model for the electroweak part of the Standard Model, and will provide a testing ground for the method. In particular, we will investigate whether this method correctly reproduces the anomaly equation, charge conservation and the correct dynamics, and determine how large the fermion ensemble needs to be to get reliable results. We also want to understand when it is correct to neglect fermion backreaction for the boson dynamics.
The paper is structured as follows: In section 2, we will introduce the model, discretize it on the lattice (section B), and derive the equations of motion. In section 3 we introduce an adapted version of the "Male" and "Female" fermion fields [2] required to generate the fermion correlators with c-number fields. In section 4 we describe the numerical setup and the results, and we conclude in section 5.

The Axial-U(1)-Higgs-fermion model in 1+1 dimensions
We will consider the 1+1 dimensional Abelian-Higgs model, coupled axially to fermions. The action reads in the continuum: in terms of the components and with the definitions The action is invariant under gauge transformations of the form if we take q = 1/2. In lattice simulations it is more convenient to work with vector gauge symmetry, rather than axial, and so noting that the left and right chiral components have opposite charge, it is therefore natural to charge-conjugate one of them [1], where C is the charge-conjugation matrix given in Appendix A. Upon doing this, the action in the new variables (but omitting the primes) reads It is no longer axially coupled, but the Yukawa interaction has become a Majorana term, and the gauge symmetry has become vector-like with the continuum equations of motion being We have introduced the gauge currents There is one further symmetry of this system, the one that this work is principally interested in, and it is the global U(1) symmetry, ψ → exp(−iωγ 5 )ψ. This symmetry has an associated current which is precisely the fermion current in the original 1 theory, and classically conserved if one naively applies the equations of motion. Quantum mechanically, however, it is the subject of an anomaly 15) and this allows us to relate the total fermion number, Q(t) = dx j 0 5 , to the Chern-Simons number, There is one further number that is worth mentioning, the winding number of the Higgs field. When the Higgs field is away from zero, it takes values on a circle parametrized by its phase θ, φ(x) = |φ(x)|e iθ(x) . Using this phase we may define a Higgs winding number, describing the number of times the field winds around this circle on a given spatial section, In a vacuum state we know that the covariant derivative of the Higgs field vanishes, and that its modulus is constant, in which case we have that ∂ x θ = A x , leading to the sum of the Higgs winding and Chern-Simons numbers vanishing in the vacuum. For the numerical work, we discretize the Abelian-Higgs-fermion model on a 1+1 dimensional lattice of size L = a 1 n x at the level of the action, and derive lattice equations of motion as described in Appendix B.

Bosons and fermions
We are interested in the time-evolution of this system, and we will adopt the approach of [1], where the dynamics of bosonic and fermionic degrees of freedom are treated differently. The gauge and scalar fields are evolved using the classical equations of motion described in the previous section. Classical dynamics is an excellent approximation to the quantum dynamics for processes dominated by infrared physics and for fields with large occupation numbers. The fermions are treated completely quantum-mechanically, in the sense of solving the quantum equation of motion (2.11) in the classical bosonic background, in terms of field operators. Since the fermions are bi-linear in the action, the equation of motion is linear, and the field can in all generality be expanded in terms of a set of mode functions and time-independent creation-annihilation operators (see below).
This leaves the question of the back-reaction of the fermions on the classical bosonic fields. Following [1] again, we interpret the fermionic terms in the gauge and scalar equations of motion as expectation values of the corresponding operators, evaluated in some state encoded in the expectation values of the creation-annihilation operators. These states are time-independent, and amount to specifying an initial condition. The time-evolution is in the mode functions only.
We then take one step further by representing these creation-annihilation operators by a set of random numbers, thereby generating an ensemble of fermion field-realisation [2]. These can each be evolved in the same bosonic background, and the field expectation values are then replaced by simple averages over the ensemble. The point is to note that the number of field realisations (N q ) in the ensemble can be much smaller than the number of mode functions (n D x ), and the statistical approach can therefore be much cheaper in terms of computational effort.

Boson initialisation
We will consider two setups for the bosonic fields. The first (in section 4.1) is to by hand set the gauge-Higgs evolution to be a sequence of sphaleron transitions, thus forcing the Chern-Simons number to change (as in [1]). The fermion fields evolve dynamically in the background of these handmade sphalerons. We will use this setup to test the ability of the ensemble to capture the anomaly, and to find out how large the ensemble needs to be.
When considering the non-perturbative field dynamics (in sections 4.2 and 4.3), we instead initialise the bosonic fields by setting A µ (x, t = 0) = 0, ∂ 0 φ(x, t = 0) = 0 and introducing random noise for the scalar field, φ(x, t = 0) in terms of random numbers φ 1,2 k , with the correlator The gauge field momenta ∂ 0 A 1 (x, t = 0) are found by solving the Gauss constraint (2.12) with fermion sources. As described in [21,22] this initialisation represents 2 an initial quantum vacuum before Higgs symmetry breaking, V ini = λv 2 φ * φ. In the subsequent evolution, momentum modes k 2 < λv 2 will grow exponentially, and from some time on they can be described using classical dynamics. The fermions do not grow, and are still treated quantum mechanically.
The amount of growth of the scalar modes is determined by the (in 1+1 dimensions dimensionless) parameter v. This can be seen in various ways. The growth lasts until backreaction from self-interactions kick in, i.e. when φ 2 ≃ v 2 . For a given mode, we have where initially, n k = 0. Classical dynamics is a good approximation once the mode has grown so much that n k + 1/2 ≫ 1/2. Hence large v allows for classicality. Another way of phrasing this is to note that once φ ≃ v, the scalar-gauge interaction and the effect of the scalar on the fermions goes as ev and Gv, respectively, whereas back-reaction of fermions on bosons is e and G. Hence for large v, fermion effects are relatively smaller (the fields have relatively smaller amplitude).
In the following, we will employ v = 64 and v = 8. Since only modes with k 2 < λv 2 are unstable, only they will be classical, and these are therefore the only bosonic modes we initialise.

Fermion mode expansion
Now we need to know how to set up the initial conditions for the fermion field, and for this we will be using the usual mode expansion. There is a slight complication, however, due to the fact that the fermion equation of motion is not linear in ψ, but involves both ψ and ψ ⋆ (2.11). This leads to the real and imaginary components having different equations of motion, particularly when the lattice Wilson term is included, and so it is convenient to write the Dirac spinor as a combination of two Majorana spinors, which in our conventions (Appendix A) just means breaking ψ into real and imaginary parts.
and it is these components that we write as a mode expansion in terms of the constant spinors U and V (given in Appendix A) and a set of creationannihilation operators b † k , b k . We then note that the fields Ψ 1,2 are canonically normalized, and that their conjugate momenta are iΨ T 1,2 , so the canonical anti-commutations relations are which may be achieved by imposing In the equations of motion for the bosonic fields we require the quantum expectation value of fermion bilinears, and so we follow [2] in constructing the two-point functions leading to where we take b k | = 0. We note that although the fields are real, the two-point function is imaginary, D ⋆ αβ (x, y) = −D αβ (x, y). The observation of [2] is that we can construct a bi-linear of classical spinor fields, for which the ensemble average two-point function matches (3.10). This allows us to simulate the quantum backreaction of fermion fields using ensemble averages of classical spinor fields; this is what we shall now do.

Fermion ensemble, Male and Female
If we were to simply evolve an ensemble of fermions, where we draw the initial conditions of each realization from a sample with the appropriate distribution and then take the ensemble average Ψ(x)Ψ(y) to mimic the quantum two-point function, we cannot reproduce (3.10). However, if one introduces two "genders" of fermions, male and female, and writes their mode expansion as then we find that taking so we now have an explicit way of replacing quantum averages, |X| , with ensemble averages, X . This leads to us evolving rather than the equations of motion appearing in Appendix B. The fermion gaugecurrent is also modified in this prescription, with the requirement of its conservation leading to Furthermore, we need a representative of the anomalous current, Because of cancellation between lattice doublers, this quantity is conserved for vanishing Wilson term (r 1 = 0, see Appendix B), but for r 1 = 1 the current is anomalous, as we will see below.  We first want to check that our approach of replacing mode functions by a random ensemble still leads to correct dynamics and that the anomaly equation is satisfied, as in [1]. We also want to determine how large the ensemble needs to be to get statistically reliable results for the anomaly and the dynamics.

Hand-made Sphaleron transitions
An elegant way of doing this is to by hand set the gauge-Higgs field evolution to be a series of sphaleron transitions, thereby continuously changing the Chern-Simons number in a controlled way. The explicit expression for the bosonic fields can be found in [1]. The important point is that sphaleron transitions take place at half-integer values of t/t 0 , and the fields are in vacuum at integer values. We choose the timescale t 0 so that the transitions are slow enough that the fermions, which are evolved using the equations of motion in the sphaleron-vacuum background, do not lag too much behind, m A t 0 = 4. From the point of view of the fermion, the evolution is almost adiabatic, and no additional spurious particle creation takes place. Only the particles associated with the anomaly contribute. At the sphaleron configuration, the Higgs field length vanishes, and Higgs winding changes discontinuously from one integer to the next. In the vacuum, N w = −N cs , and we will always plot −N w . Fig. 1 shows the evolution of Chern-Simons number N cs , Higgs winding number −N W and the fermion number N f . The parameters used were m A L = evL = 25.6, G = 0, N q = 90, n x = 128, v = 2, timestep dt = 0.05. The agreement between Chern-Simons number and fermion number is remarkably precise, even for such a small ensemble.
As was pointed out in [1], fermion number is periodic on a finite lattice with period 2n x , and so for large N cs the agreement will fail. Fig. 2 shows fermion number for very large Chern-Simons number at different values of the lattice size n x (volume is fixed evL = 6.4). These show the lattice behaviour, and can indeed be rescaled by n x (along both axes) to end up on top of each other (inset). This is exactly as in [1], and means that sufficiently large lattices can accommodate any Chern-Simons number. Notice that the ensemble is still N q = 90.  Although the small ensemble very convincingly reproduces the anomaly when looked at by eye, it is only prudent to investigate the statistical precision. This is shown in Figs. 3 and 4 for ensembles of 10, 30, 90, 270, 810 and 2430 random realisations, respectively. The left-hand plot is on a n x = 32 lattice, and we see that the agreement is always fairly good, even for N q = 10. Looking closer (inset), we do see that the curves converge, and in fact converge to a value slightly off N cs . This is the finite volume effect as described before. In the right-hand plot with n x = 128, this discrepancy is gone and increasing the ensemble, fermion number converges to the Chern-Simons number value. We conclude that convergence in N q is achieved at the few-percent level for N q = O(1000).
As reported in [1], including the Yukawa coupling G makes the lattice artefacts   stronger. Fig. 5 shows the fermion number in the hand-made sphaleron background for n x = 32, evL = 6.4 with varying G. The anomaly holds until N cs ≃ 5, after which the finite size effects kick in, stronger with increasing G. However, increasing the lattice size again ameliorates the situation, as shown in Fig. 6, where G/e = 0.1 is kept constant, and the lattice discretization is made finer (constant volume evL = 6.4 and n x increasing 3 ). like to be able to run baryogenesis simulations for timescales m A t = evt = O(100), on a large enough lattice to fit in the appropriate physics m A L = evL ≫ 1, while having the anomaly correctly reproduced for realistic values of the Chern-Simons number, say N cs = O(10). And we also need to include a non-zero Yukawa coupling, at least for physics around the electroweak transition. The question is whether we can find a combination of evL, n x , G and N q that can accommodate this. Fig. 7 shows a run on a much larger lattice, evL = 51.2, v = 2, with Yukawa cou-pling G/e = 0.1, in a range of N cs = 0 − 10. We first note that the larger volume makes the anomaly agree less well than for the runs in Fig. 6, which did reasonably well until N cs = 5. This is the case both for n x = 256 and for n x = 512, with half the lattice spacing (not shown). However, this is just adding up of statistical fluctuations, and can be compensated for by increasing the ensemble size. At N q = 810 the agreement is again convincingly reproduced. The finite volume effect is not apparent at these lattice sizes.  In a fast-quench symmetry breaking transition, Higgs field modes with k 2 < λv 2 will be unstable ("tachyonic") and grow exponentially. This drives the gauge field to also grow until non-linear backreaction begins to dominate, stop the growth and eventually leads to thermalisation. In the presence of CP-violation, such a transition may lead to a baryon asymmetry (see for instance [13]).

Non-equilibrium dynamics
For our choice of initial conditions, the initial gauge field is driven by the fermion ensemble fluctuations and the initial scalar field. Our goal is that for a given scalar field configuration, the evolution should be independent of N q , so that the statistics reliably reproduce the fermion state. We need a large enough ensemble to have statistical fluctuations under control. We will set λ e 2 = 1 4 , n x = 256, 1) and vary N q . Fig. 8 shows the Higgs field (black line), Higgs winding number (green), Chern-Simons number (red) and fermion number (blue) in such a transition. G = 0, v = 64 and N q = 2430. The Higgs field "falls off the hill" as expected, and performs oscillation around its finite temperature minimum. Meanwhile, the Chern-Simons number grows and oscillates near the integer-value Higgs winding number. The anomaly is so well obeyed that fermion number is essentially indistinguishable from Chern-Simons number. The Higgs winding bounces around in the beginning, an effect of the Higgs field length being small, and the winding number therefore ill-defined. But once symmetry breaking gets going, winding number is stable, integer and consistent with the Chern-Simons number.
We illustrate the convergence of the dynamics with increasing N q in Fig. 9, where tachyonic transitions are performed for v = 64 for different sizes of the ensemble. As expected, we see convergence in N q , but also that the required ensemble is O(1000), to get agreement at this value of v and for these times. In fact, the Chern-Simons number is very sensitive to fluctuations in the fermion source. This is at least partly because in 1+1 dimensions, a U(1) gauge field only has one dynamical degree of freedom (i.e. up to gauge transformations), which is precisely the Chern-Simons number. This means that in the sea of fermion degrees of freedom, the single degree-of-freedom gauge field can easily be bounced around. These issues are specific for 1+1 dimensions, and we will proceed with v = 64 and N q = 2430, for which convergence is under control at least for m A t < 100, and qualitatively correct for m A t < 150. This will suffice for the present work, but can be improved depending on the level of precision required.
Since the Yukawa coupling is absent, the fermions are massless throughout. Also, the fermion and boson total charges are individually conserved at the level of O(10 −13 ), and Gauss law is conserved at a (relative) level of O(10 −8 ).
For non-zero Yukawa coupling the fermions acquire masses as the Higgs transition proceeds (in addition to the gauge and Higgs fields). As we saw in section 4.1, the Yukawa coupling introduces stronger lattice artefact, but these could be cured by using larger ensembles. In Fig. 10 and 11 we show the evolution and convergence in N q of a simulation with G/e = 0.1, v = 64. We see that although there is a clear effect of non-zero G, convergence still holds by increasing the ensemble to a few thousand, and fermion number (dashed lines) follows Chern-Simons number fairly well.

Application: Cold Electroweak Baryogenesis in 1+1 dimensions
In the minimal model of electroweak baryogenesis, CP-violation is provided through the CKM fermion mass matrix. For hot baryogenesis, this effect is much too small to account for the asymmetry and a separate source of CP-violation is required. The   situation is less clear for "cold" baryogenesis (see for instance [13]). For illustration, we will postpone this issue, and simply introduce C(P) violation 4 in our 1+1 dimensional model through a bosonic term in the action (exactly as in [21]), which amounts to an addition to the bosonic equations of motion (2.12), (2.10) of The conservation of Gauss law and the anomaly and the convergence in N q is unaltered by this addition, and the C(P) violation is not obvious from a given random realisation of the bosonic fields. We now need to also average over an ensemble of bosonic realisations, each with a separate ensemble of fermion fields. This represents a quantum initial state of the Higgs fields coupled to fermions initially in the vacuum. Fig. 12 shows the scalar-ensemble averaged observables, Higgs field, Chern-Simons number, fermion number and winding number, at κ = 0.04, v = 64 and N q = 2430. We average over a set of 64 random realisations plus the corresponding C(P)-conjugate configurations. This makes the ensemble explicitly C(P) symmetric, and the asymmetry will be identically zero for κ = 0. This procedure is similar to the one employed in [24,25]. We see that an asymmetry is indeed generated in Chern-Simons, winding and fermion numbers as the transition proceeds. The anomaly is very well obeyed until times m A t ≃ 100, where fermion number goes a little low. We checked that this is indeed due to the configurations with relatively large N cs = 8 − 10, for which the lattice artefact makes a small deviation.
From Fig. 13, a blow-up at early times, we see that because Higgs winding can only take place in the presence of a local Higgs field zero, the asymmetry is created when the average Higgs field is low in its oscillation. When it is high, winding number is essentially constant. Also, since the C(P)-violating force is proportional to|φ 2 |, the gauge field picks up speed in-between Higgs extrema. The net effect is a "pumping" behaviour, of the baryon asymmetry as Higgs symmetry breaking proceeds. As the Higgs approaches a uniform vev, and oscillations damp out, the asymmetry creation gradually stops. This is very similar to the case in the 3+1 dimensional SU(2)-Higgs model [12,24,25], where the gauge field is much more complicated dynamically than the model considered here.
One important question is to what extent fermions can be ignored dynamically, compared to the bosonic fields. If not, many simulations of baryogenesis may need  For v = 8 however, the fermions begin to influence the evolution of the gauge field, even when the bosonic fields are subject to a tachyonic instability and therefore grow large. For the Standard Model in 3+1 dimensions, there is no v-ambiguity, and it will be crucial, to what extent back-reaction is important. In particular, CP-violation itself is a backreaction effect, which will dynamically generate effective terms similar to the bosonic C(P)-violation used here. Finally, to illustrate the type of calculations that are possible, we show how the asymmetry depends on κ (Fig. 15). For small enough κ, the dependence is nicely linear. Also remember that because our scalar ensemble is explicitly C-symmetric, N f (κ = 0) = 0.

Conclusion
By combining the methods of [1] and [2], we have demonstrated how to do first-principle numerical simulations of bosonic scalar-gauge systems with quantum fermions, in a numerically efficient manner. Although there is no gain in numerical effort in the specific 1+1 dimensional toy model considered here (compared to [1]), in the physically relevant case of 3+1 dimensions, we expect a significant decrease in the required computing time.
As an example, including fermions in the simulations of [12], an ensemble of N q = 2430 should be compared to the 90 3 (= 729, 000) mode functions otherwise required, a gainfactor of 300. Either way, fermions are numerically challenging, but for the setup of [12], there would then be no need for the resource-consuming CP-violating term. Except for the scalar ensemble averages in section 4.3, all the simulations presented here were done on a normal desktop computer in less than 24 hours in total.
We found that Gauss' law, fermion back-reaction as well as the baryon anomaly are well reproduced in terms of a statistical ensemble of O(1000) fermion field realisations. As described in [2] implementation of the fermion correlators, including the anti-commutativity of the fermionic operators requires a doubling of the fields into "male" and "female" (not to be confused with the standard lattice fermion doublers), adapted to the system at hand. In addition, the usual lattice doubler problem has to be addressed; in the present case we found that Wilson fermions in space and a small timestep was sufficient to keep the doublers sufficiently decoupled that they stayed unexcited for the timescales required here. Failure to do this leads to an exact cancelling out of the baryon anomaly.
The method requires careful consideration of the interplay between lattice size, ensemble size and the size of couplings. In particular, the Yukawa coupling introduces additional lattice artefacts, which have to be compensated for, and we have demonstrated how to do this.
The upshot is that fermions are included completely in the dynamics, since they are bi-linear in the action, and so at least in cases where gauge fields are dominated by large particle numbers and long wavelength, this approach provides a very reliable description of the full field dynamics. The obvious application of the method is (electroweak) baryogenesis, where baryon number violating processes are classical in nature, whereas the CP-violation 5 and the actual baryon number are carried by the fermion degrees of freedom. This applies both to "hot" and "cold" baryogenesis.
Simulations of the early stages of the heavy-ion collisions are also within the scope of the work presented here, since it involves very large (boosted) gluons fields coupled to (sea and valence) quarks. The valence quarks source the gauge field, which then evolves and may in turn source the emission of fermions. With the method presented here, fermions may be included in the dynamics completely.
The obvious next step is to implement ensemble fermions in 3+1 dimensions, coupled to SU(2)-Higgs bosonic fields as in the Standard Model, where the gauge-fermion interaction is chiral rather than axial. Including the Standard Model CP-violation via the CKM matrix will require all three fermion generations, and represents a significant numerical challenge; with the method described here, this numerical effort can be reduced by one or even two orders of magnitude. This is work in progress.

A. Conventions
We use the metric signature (−, +), and for the Dirac algebra, we employ the Weyl-Majorana representation with explicitly

B. Lattice equations
On the lattice, we define the gauge link variables (a µ , µ = 0, 1 are the lattice spacings), We are using the non-compact formulation of the gauge action, with A µ the basic gauge field variable. We define the derivatives We will deal with the spatial fermion doublers by including a Wilson term The lattice action then becomes, This immediately gives the lattice equations of motion and the currents, and so we have Gauss' law from The chiral current is Chern-Simons number and fermion number is trivially adapted to the lattice, and following [23], we write φ(x) = |φ(x)|e iθ(x) , and then define the integer lattice winding number as

C. Fermion doublers
Doublers contribute to the anomaly with the opposite sign to the non-doubler modes, and the anomaly will then average out if we do not remove the doublers from the dynamics. Fig. 16 shows the anomaly in a tachyonic transition for "naive" fermions r 1 = 0 and with a Wilson term r 1 > 0. The anomaly disappears for the naive fermions, when doublers are allowed to get excited. By adding a Wilson term in space, but not time, the fermion equation will still lead to temporal doublers, so where we thought we were evolving a single Fermi-field ψ, we are actually evolving two Fermi-fields, which we call ψ + and ψ − . To see this we define Now, if the lattice time derivative is evaluated on an even t slice, then it actually only uses fields evaluated on the preceeding odd t slice, and the following odd t slice, in which case one finds the equation of motion (B.11) becomes Similarly, on odd t slices we have 0 = γ µD µ ψ + (x) − a 1 r 1 2 D 1 D ′ 1 ψ + (x) + Gφ(x)ψ + * (x) (C.3) showing that both ψ + and ψ − satisfy the fermion equation of motion, and that we are actually evolving two fermi degrees of freedom. In Fig. 17 we see a set of runs where we vary the timestep. We only initialise the non-doubler modes, and we see that for all timesteps that give a stable integration of the equations of motion, the time-like doublers stay un-excited. In all simulations in the main paper, we use dt = 0.05, the smallest time-step presented here, and we see no doubler effects.

D. Spinors
We now construct the basis spinors required in the mode expansion of the fermion operators, which we do by setting A 1 = 0 and taking |φ| 2 = v 2 /2 so that for Ψ 1 , m f = Gv/ Similarly, the negative frequency solutions are found by x , k 0 = ω 1 > 0, (D.5) but because the field is Majorana, and therefore real, we immediately have V = U * . In order to calculate the two-point function we need the following identity, wherek µ = (ω 1 , s k ). To get the mode functions for Ψ 2 , simply make the replacement m f → −m f .