Cold Baryogenesis from first principles in the Two-Higgs Doublet model with Fermions

We present a first-principles numerical computation of the baryon asymmetry in electroweak-scale baryogenesis. For the scenario of Cold Baryogenesis, we consider a one fermion-family reduced CP-violating two Higgs-doublet model, including a classical SU(2)-gauge/two-Higgs sector coupled to one quantum left-handed fermion doublet and two right-handed singlets. Separately, the C(CP) breaking of the two-Higgs potential and the C and P breaking of the gauge-fermion interactions do not provide a baryon asymmetry. Only when combined does baryogenesis occur. Through large-scale computer simulations, we compute the asymmetry for one particularly favourable scalar potential. The numerical signal is at the boundary of what is numerically discernible with the available computer resources, but we tentatively find an asymmetry of $|\eta|\leq 3.5\times 10^{-7}$.


Introduction
Electroweak baryogenesis is still the subject of significant scientific attention, being a very elegant and testable set of scenarios to explain the baryon asymmetry of the Universe. The central element is the non-perturbative violation of baryon and lepton number in the electroweak sector of the Standard Model, which when combined with C-and CP-breaking interactions out of equilibrium allows for baryogenesis [1]. Since the Standard Model and its simplest extensions (multiple scalar fields, massive leptons) break C and CP, there only remains to compute the asymmetry. This turns out to be very challenging to do from first principles, since it involves non-perturbative quantum dynamics of in particular fermions out of thermal equilibrium, and because the final asymmetry one is trying to compute is very small [2] η = n B n γ (6.0 ± 0.1) × 10 −10 . (1.1) Traditionally, the approach has been to split the computation up into a number of equilibrium and out-of-equilibrium quantities, which may each be treated using dedicated techniques. The state-of-the-art is quite advanced, and often involves considering effectively bosonic systems obtained through integrating out the fermions and/or dimensional reduction [3]. This has allowed a quite precise determination of the sphaleron rate [4]; the phase diagram of the Standard Model and its extensions [5,6]; the bubble nucleation rate (in case of a first-order transition) [7]; the dynamics of bubbles and their observational signatures [8,9]; the effective CP-violation in the Standard Model [10,11]; the interaction of a bubble wall with fermions (see for instance [12][13][14]). Putting these quantities together gives a handle on the baryon asymmetry produced, and some are even no-go results, such as the non-existence of a first-order phase transition and the very strong suppression of Standard Model CP-violation. Historically, these no-go results removed the need for actually calculating a number for the asymmetry, apart from order-of magnitude estimates based on maximally favourable parameters. The need for a strong phase transition is in the Cold Baryogenesis scenario replaced by the requirement of a cold tachyonic transition at the end of inflation [15][16][17]. A compelling feature of this scenario is the option of numerically computing the final asymmetry from first principles, simply because the mechanism is very simple: the particle creation process of reheating during an electroweak spinodal transition is asymmetric in particles and antiparticles. Such simulations have been performed for more than a decade [18][19][20][21][22][23][24], showing that the scenario is viable, given a cold initial condition. The simulations were however done in purely bosonic versions of the Standard Model or extensions with a singlet [25] or a doublet scalar field [26].
For electroweak baryogenesis, it is crucial to note that C-, CP-but also P-breaking is necessary to generate a baryon asymmetry. This follows from the anomaly equation relating the Chern-Simon number of the SU(2) gauge field and the baryon and lepton numbers B(t) − B(0) = n f [N cs (t) − N cs (0)] = L(t) − L(0), (1.2) and the observation that Chern-Simons number is odd under P but even under C. What this means in practice is that it is not enough to break C, thereby breaking CP if P is conserved. And it is also not enough to break P and thereby break CP if C is conserved. One needs to break C and P and CP, and this is indeed achieved in the Standard Model through the left-handed coupling to gauge fields (C and P broken maximally, CP conserved), and the complex phase in the CKM matrix (CP broken, by breaking C). Alternatively, in the two-Higgs doublet model (ignoring the CKM matrix), the left-handed coupling to gauge fields again provides C and P breaking, and in addition, complex parameters in the Higgs-Higgs potential break C and thereby CP. And only when combining them is an asymmetry produced.
The results of different combinations of the breaking of discrete symmetries was demonstrated in [19] for the bosonic part of the Standard Model with an effective CP and Pviolating operator. One finds non-zero N cs , and the fermion number was inferred through the anomaly equation. The subsequent task is then to compute the coefficient of such an operator from integrating out the fermions in the full theory, including both C, P and CP-breaking. This turns out to be a substantial calculation, yielding a slightly different set of effective operators [10,27,28].
In the bosonic part of the two-Higgs doublet model, it was seen that simply adding the C-breaking potential does not produce a net Chern-Simons number; only when in addition including an effective C-and P-breaking, but CP-conserving operator is an asymmetry produced. In [33], the operator was considered. Again, it remains to compute the coefficient δ C/P of this operator, this time from the CP-even fermion sector 1 . This has not been attempted beyond the simplest estimates [24], but may in principle be done using the methods of [10,28]. The obvious, but as it turns out far from straightforward, alternative is to include the fermions in the real-time dynamics. Fermions are inherently quantum mechanical, but may be formulated on a classical bosonic background [29]. The numerical solution is computationally extremely demanding, but possible using what is known as "ensemble" fermions [30]. Crucially, the anomaly equation carries through and the fermion back reaction on the bosons is reliably described [31,32]. The C-and P-breaking of the gauge-fermion coupling may now be directly introduced, and in combination with the C-breaking of the two-Higgs potential, an asymmetry should be created. The questions we wish to address are then the following: • Is an asymmetry created in a cold spinodal transition?
• Is this asymmetry large enough that we can see it numerically on the lattice?
• Is it comparable to or ideally larger than the observed asymmetry?
• Can we connect these results to previous work [33], for instance to provide an estimate of the coefficient δ C/P ?
• What is the numerical effort involved, and is it realistic to sweep a multidimensional experimentally allowed parameter space?
Many of these questions will depend on the speed of the spinodal transition, the numerical effort available and the parameter range adopted for the Higgs-Higgs potential. But for moderately favourable choices, the answers are: yes, perhaps, yes, yes and unlikely. In Section 2 we will introduce the reduced Standard Model that we will consider, including only one generation of fermions. We will set out the observables, numerical implementation, initial condition and the Higgs potential. In Section 3 we explain how to compute the baryon asymmetry and describe our numerical results. We conclude in Section 4. A number of technical details are relegated to a set of Appendices.

The reduced Standard Model
We will consider a reduced version of the Two-Higgs Doublet Standard Model, where the SU(3) and U(1) interactions are ignored, and we only keep one generation of fermions, including a left-handed quark SU(2) doublet, q L = (u L , d L ) T , a left-handed lepton SU (2) doublet, l L = (e L , ν L ) T and right-handed singlets u R , d R and ν R , e R . Although we use a notation suggesting the first generation of the SM (u, d, e, ν e ), we expect the main effect of CP-violation to come from the heaviest fermions, in practice the third generation (t, b, τ, ν τ ). We emphasise that although for SM CP-violation it is crucial to have 3 generations of fermions (allowing for a complex phase in the CKM matrix to be physical), for the 2HDM CP-violation considered here, this is not necessary. Adding up to three families is straightforward, but require three times as much computational time.

Continuum action
Having made these simplifications, the continuum action reads where the SU(2) gauge covariant derivatives are and we have ordinary partial derivatives for the right-handed singlet fermion fields u R , e R , d R and ν R .

Higgs potential
The two-Higgs scalar potential is given by This is the standard parametrisation [34], where we have chosen to put the couplings λ 6,7 to zero for simplicity. Allowing λ 5 and/or µ 2 12 to be complex provides for C-violation at tree-level. If in addition Im(µ 2 12 ) = 2|v 1 ||v 2 |Im(λ 5 ), both Higgs vevs v 1,2 can be chosen real. Otherwise the vevs will have a relative phase Arg(v † 2 v 1 ). We will choose this argument to be maximal, Not all sets of complex λ 5 and µ 12 lead to C-violation. A useful parametrisation is provided by [35], which guides our choices below. We wish to maximise the effects of CP-violation, by tuning the parameters of the Higgs potential (the C-and P-breaking in the gaugefermion coupling is already maximal). In addition we wish to have well separated Higgs mode masses, and we choose for the neutral Higgs bosons We have also fixed the total Higgs vev, so that |v 1 | = 110 GeV, |v 2 | = 220 GeV, tan β = 2, (2.7) with a relative phase of π/2 and After a number of tests at small lattice volumes, we settled on which amount to the angles mixing angles of the three neutral Higgs modes [35] (α 1 , α 2 , α 3 ) = (0.314, −1.49, 2.58), (2.11) The potential breaks C, due to the non-zero values of Arg(µ 2 12 ) = 0.45π, Arg(λ 5 ) = −0.14π. (2.12) In Figure 1, we show a contour plot of the potential in the |v 1 |, |v 2 | plane. The potential energy in the symmetric phase is (2.13)

Yukawa couplings
The Yukawa coupling terms can in all generality be parametrized as (2.14) For simplicity, we will use a single coupling constant, There is no issue in principle using general Yukawa couplings, although it may become a little cumbersome. We take It was demonstrated in [31] that the ensemble averaging procedure is under control for such large Yukawa couplings, which correspond to fermion masses of about 11/ √ 2 8 and 22/ √ 2 16 GeV, depending on which Higgs field they are coupled to, easily larger than all the quark masses except for the top mass.

Asymmetric observables
The baryon and lepton numbers are the spatial integrals over the zero-component of the baryon and lepton currents We will assign baryon number 1 to the quarks to emulate summing over SU(3) colour. And lepton number 1 to the lepton field. Fermions coupled chirally to an SU(2) gauge field will experience a quantum anomaly, so that these currents obey for each fermion doublet. The Chern-Simons current is given by and we therefore have that under a change of Chern-Simons number N CS over time, we pick up a change in baryon and lepton number of (2.23) The Chern-Simons number changes continuously between one gauge vacuum and the next, and is an integer in the infinite series of (almost) degenerate gauge vacua. In contrast, the Higgs field winding number is always an integer, changes discontinuously halfway between vacua, and coincides with Chern-Simons number in the vacua. Because it is an integer, the winding number turns out to be a better observable to count baryons, since if at the end of a simulation one has a well-defined N W , one knows that because the sphaleron barriers are now up, N cs and therefore B and L will eventually relax to the same integer value. We also note that there are two winding numbers N 1,2 W , one for each Higgs field, and they both coincide with the Chern-Simons number in the vacuum. Because of lattice artefacts, these are in fact not numerically exactly integers, and the jumps not strictly discontinuous. In addition to the expectation of the Higgs fields and the various energy components, our prime observables will be Chern-Simons number, the two winding numbers and the fermion number of the quark and lepton fields.

CP-symmetric initial conditions
We will consider an instantaneous quench of the Higgs potential at zero temperature. This in practice means that we evolve the initial configurations using the potential (2.1.1) for times t > 0 but that we generate these initial conditions in the vacuum corresponding to the quadratic potential The vacuum is defined through the "half" method of [18], which is to have field correlators and momentum correlators obey where π i k and φ i k are the momentum space variables for each real component of the Higgs fields (8 in total), and ω 2 k = M 2 + k 2 , with M 2 the eigenvalues of the mass matrix. The understanding being, that if m 12 is non-zero, one should first diagonalise in field space, generate the configuration in this basis and then rotate back to the φ 1,2 basis. This prescription is the natural analog of the quench in [22], where in a single-Higgs model, the coefficient of the quadratic term "flips" instantaneously.
For each randomly generated configuration φ i , π i , one also generates an entire ensemble of N f = 4000 fermion configurations (u, d, e, ν separately, with "male" and "female" copies [30]). Since the initial state is the symmetric vacuum ( φ 1 = φ 2 = 0), the fermions are taken to be massless initially. Setting A µ = 0 initially, we solve Gauss law for the SU (2) electric field E i in the background of the scalar-fermion fields.
The behaviour of each of our fields under CP-transformations is and In order to reduce statistical noise, we will construct an explicitly CP-symmetric ensemble of initial conditions, where for every randomly generated initial configuration of φ 1,2 and u, d, e, ν and E i , we include its CP-transformed. The CP-symmetric ensemble ensures, that when C-symmetry is turned off in the Higgs potential, the ensemble averages of the CP-asymmetric observables N 1,2 W and N cs are identically zero. We have checked numerically that this is indeed the case (see also [31,32]). This also means that there is no spurious CP-violation from counter terms or discretisation.

Results
We implement the action on a space-time lattice, and derive the classical equations of motion for the bosonic fields and the linear operator equations for the fermions (see appendix through the "ensemble" fermion method [30]. For details on the numerical implementation of this method, we refer to [31], and to the appendices in the present paper.
The fermion expectation values are formally divergent in the continuum limit. This can be resolved by introducing counter terms in the bosonic equations of motion. This is described in Appendix B.
For each random initial realisation of the Higgs field and fermion ensemble, we solve the evolution equations in temporal gauge A 0 = 0. We use a spatial lattice of size (Lv) 3 = (n x av) 3 with N x = 32 and av = 1.2 denoting the Higgs vev in lattice units. MPI-parallelised runs each take about 8 hours on 80 cpu's. The total cpu-time used for 400 CP-conjugate pairs is therefore approximately 500.000 hours.
3.1 Which observable to use: N W , N cs or N f ?
In Figure 2, we show the average Higgs field, the winding numbers, the Chern-Simons number and the fermion number for a single configuration. As was demonstrated in [31] the fermion number follows the Chern-Simons number in accordance with the anomaly equation, but the statistical noise of the observable N f is substantial. This noise originates in the UV of the lattice operator, and is largely white noise that may be reduced by increasing the fermion ensemble. What is important for our purposes here is that the back-reaction on the bosonic degrees of freedom has converged, and this is achieved to a reasonable degree for N f 4000 [31].
The bosonic fields evolve in a plausible way: The two Higgs fields roll off the potential barrier at φ 1,2 = 0, and oscillate and damp asymptotically to their respective vacuum expectation values, which is here normalised to unity (blue and black). The Higgs transition is over by vt 10, but the topological observables do not settle until vt 20, i.e. after twice as long. Chern-Simons number (red) changes from zero to (in this case) 4 in the interval 10 < vt < 20, following the two Higgs winding numbers, but lagging a little behind (yellow and green). Winding number changes most readily when the average Higgs field is small, and we indeed see a failed attempt at a winding number transition around the first Higgs minimum at vt 7.
Both winding numbers and Chern-Simons number suffer some amount of lattice corrections, which tend to reduce the vacuum asymptotic values to somewhat below integer. Figure 3 shows a histogram of final winding and Chern-Simons numbers, over all the configurations simulated. We see a clear set of peaks, corresponding to integer values. In fact the values in the figure have been rescaled by 1/0.73 to compensate for lattice artefacts. The overall distribution of integers looks roughly Gaussian, and although the difference is not large, the winding numbers tend to be somewhat more peaked than Chern-Simons number for each integer. Because of the statistics noise on the fermion number and the lattice artefacts, we choose to use N 1,2 W as the cleanest observables. As we will see below, since we are looking for complete integer flips, this choice is convenient, but N cs could be used as well.

Looking for flips
Next, we demonstrate that "flips" do occur, configurations where the observables and the CP conjugate configuration observables do not add up to zero after the transition. First, Figure 4 shows a pair of configurations that average out to zero in both winding number and Chern-Simons number. The observables for the CP-conjugate pair have been given the opposite sign for comparison. We see that not only is there no flip, the observables are also closely similar, almost to the point that they are indiscernible. The effect of C, P and CP violation stays quite small.  N cs (N) the kind of events we are looking for, and these will stay as a permanent asymmetry in the fermions. In the present case the flip is of +1 in the total winding number. Closer inspection shows that the mismatch begins to accumulate around vt = 20, having been vanishing for earlier times. Also, the transition seems to be controlled by the evolution of N cs , which first splits away from N 1,2 W in both configurations, but for one the winding number moves to rejoin the Chern-Simons number; in the other Chern-Simons number moves to rejoin the winding number.
This sort of event is reminiscent of the early work of [24]. A configuration with mismatched N cs and N W is gauge equivalent to a N cs → 0, N W → N W − N cs localised texture. One may then ask under which circumstances such a texture decays by the winding number going away, rather than the gauge field moving to screen the topological charge, i.e. change Chern-Simons number. In [24] it was argued that there is a typical physical size of the texture, deciding one or the other, and that this critical size is shifted in the presence of CP-violation. A similar behaviour was observed in [36].
In our setup, C and P violation from the fermions acts on the gauge fields and C violation enters in the scalar dynamics from the potential. We may therefore identify our flips with events where a texture is created at random, and its unwinding is just biased enough by C/P and C violation that it goes opposite ways for the CP-conjugate pair.
In Figure 6, we display the final values of N 2 W for all 400 configurations (black) and their CP-conjugate (white), rescaled by the lattice artefact correction 1/0.73. Firstly, we notice how precisely most of them overlap, so that the black symbols are practically hidden. But a few (16) do not match, and they correspond to flips, 6 times +1 and 10 times −1.

Averages
We conclude by showing in Figure 7 the observables in time, averaged over the entire ensemble of bosonic realisations. The green bands correspond to one standard deviation. We see that the average Higgs fields φ 2 1,2 are very well behaved statistically, and describe a strongly damped oscillation, approaching the equilibrium expectation value. Because of finite temperature corrections, this is slightly smaller than the vevs, to which the Higgs fields are normalised. Although the Chern-Simons number N cs seems to have a broad statistical uncertainty, the average evolution looks strikingly similar to the purely bosonic simulations [19,21,22,26,33] (up to an overall sign). First a semi-exponential growth until vt 20, then a dip to a final value of the opposite sign. In [19], the initial rise was explained in terms of a linear treatment of the effective CP-violation, and the general transfer of energy to the gauge field during the spinodal preheating. That this initial behaviour should carry over from using effective bosonic operators to the full fermion dynamics is perhaps not surprising. What is maybe more interesting is that also the qualitative "bouncing" to the opposite side occurs with the dynamical fermions providing the C-and P-violation. A similar behaviour is observed for the winding numbers. Having discovered and asymmetry of 6-10=-4 in 400 pairs, we would expect final winding and Chern-Simons numbers of −4/800 −0.005. This does not follow directly from these averaged observables, without including the lattice artefact correction 0.73, which itself is only approximate. Including this, the integer counting corresponds to N 1,2 W 0.004 at the final time. The naive statistical error bars are as large as the signal, and must therefore be taken as inconclusive. We do suspect that the naive standard deviation is an over-estimate of the statistical error, and therefore count these results as suggestive.
The fermion number is badly behaved, but seems to agree between different bosonic realisations. This shows that the lattice fermion number observable is indeed dominated by UV modes, that are all initialised identically irrespective of the random background bosonic field.

Computing the asymmetry
Our simulations are performed on a lattice of physical size v 3 V = v 3 L 3 = n 3 x ×(av) 3 . Given N + − N − flips in N tot pairs of configurations, this constitutes a baryon asymmetry of The total energy available in the simulation is V 0 = (226 GeV) 4 , which should be distributed over all the active degrees of freedom Then the baryon to photon ratio is

(3.3)
Using g * = 28 (in the reduced SM simulated here, g * 100 in the full Standard Model), v = 246 GeV, V 1/4 0 226 GeV on a lattice n x = 32, av = 1.2, we find where the one standard deviation error applies for N tot 1. The observed asymmetry of 6.0 × 10 −10 corresponds to one uncancelled flip in 10 5 configurations. This is numerically a hopeless task to find, since only fully completed flips count. An overall shift close to a minimum in the Chern-Simons number potential is not meaningful. We were able to simulate 400 such configuration pairs, so that one uncancelled flip amounts to η 8.6 × 10 −8 . However, the statistical error is large if there are many flips cancelled out ( 4 Conclusion:

4.1
Can we see a net asymmetry?
In a ensemble of 400 pairs of realisations, we found N + = 6 flips with +1 difference in N W ; and N − = 10 with -1. We therefore have about 600 times as big as the observed asymmetry; but also consistent with zero, if these simple statistical error estimates are to be relied upon. Taking the average value as physical, one may readily compare to the results of [33], obtained by replacing the fermion degrees of freedom by an effective C/P breaking bosonic operator of the form The expectation is, that to leading order in a gradient expansion, this effective term arises upon integrating out the fermions (see however [10,27,28]). In this way one could in principle compute the coefficient δ C/P . In the present work including fermions, this term and all higher order terms are in principle included at one-loop 3 . Our values of (α 1 , α 2 , α 3 ) correspond exactly to the the lowest left black point in the bottom right panel of Figure 1 in [33]. The simulations in that paper were done at µ = Im(µ 2 12 ) = 100 GeV (half of the value in the present work). Typical asymmetries quoted in [33] are |η| = (1 − 2) × 10 −4 , with a maximal value of Under the assumption that the precise value of µ is not decisive, and to the extent that the fermion degrees of freedom can be represented by such a bosonic operator, we conclude that in this scenario, the effective coefficient is of the order of This is much smaller than one, but not nearly as small as for Standard Model CP-violation at finite temperature [10], which seems to totally rule out electroweak baryogenesis in the Standard Model. This is likely due to the large Yukawa couplings and the fact that CPviolation in the SM kicks in at higher loop order than in the present case of a tree-level CP breaking effect combined with a tree-level C and P breaking effect.

Outlook
In conclusion, we have successfully developed a method to numerically compute the baryon asymmetry from first principles in a viable reduced Standard Model. The numerical effort is significant, and with our available numerical resources, we were only barely able to see a signal. But clearly, the effect of CP-violation is present, and using an explicitly CPsymmetric ensemble of realisations allows us to see this. Using a fermion ensemble of N = 4000 for each bosonic realisation is a cautious choice, where we may be confident that the fermion back-reaction on the bosonic evolution is well converged. This back-reaction is important, since it provides the breaking of P absent in the Higgs potential and necessary for the asymmetry to be generated. We tentatively find an asymmetry of ≈ 600 times the observed one, but with a very large statistical uncertainty. An increase in computing time of a factor 10 is necessary to pin this down convincingly. We are currently seeking to acquire such resources for future work.
Obvious applications and extensions of this method are to consider other experimentally allowed regions of two-Higgs doublet parameter space; including all three families of fermions with the physical values of the Yukawa couplings; including the mechanism responsible for the cold spinodal transition (low-scale inflaton, singlet Higgs field); and the interaction of fermions with a bubble wall in a first order phase transition, relevant for "hot" electroweak baryogenesis. The inclusion of U(1) and SU(3) gauge fields may also be envisaged, but then the transformation from a chiral to a vector theory no longer works, and one will have to in another way circumvent the issue of chiral fermions on the lattice.
Acknowledgments: ZGM would like to thank CSC for financial support. PMS thanks STFC for financial support. AT is supported by the Villum Kann Rasmussen Foundation. The numerical work was performed in part at Dirac facility COSMOS (UK) and the Notur facility Abel (Norway).

A Lattice implementation
Putting chiral fermions on the lattice is notoriously difficult, but because of the pseudoreality of SU(2), it is possible to bypass this problem by defining (in the continuum) a new set of fermion fields, in the following way: This amounts to a charge conjugation on the left-handed component. Now, instead of having four right-handed singlets and two left-handed doublets, we have a full doublet (left-and right-handed) and two full singlets. This turns the reduced Standard Model into a vector theory rather than a chiral one, and the lattice implementation becomes similar to 2-colour QCD.
The action now reads: with the 2-Higgs scalar potential The Yukawa interactions are now somewhat more complicated, with in general in terms of a number of projectors For the redefined fields, CP-transformation amounts to The lattice parameters are related to the continuum ones and the lattice spacings as where a t = adt is the lattice spacing in the time direction, with a the spatial spacing. We will refer to dt as the "time step".
It follows that the energy density is given by (for the G(auge), H(iggs) and F(ermion) components, respectively) summing to the total energy density. Because of the field transformation, we compute the baryon and lepton currents of the original theory as the chiral current in the transformed theory We now specialise to the Yukawa interactions with a single coupling y uk , This translates to the lattice theory as Our choice is simply meant as a simplification, and may readily be generalised using the general form (A.9). We have carefully considered accidental symmetries of the Yukawa term, to ensure that CP is in fact broken. Fermions enter in the bosonic equations of motion as two-point correlators. These are formally divergent, and we introduce counter terms ct 11,12 in the following way We do not consider weaker logarithmic divergences that appear in the gauge equations of motion. S C should be added to (A.5-A.7).
From the complete lattice action, we derive the equations of motion by straightforward variation. For the fermions, linear equations in the bosonic field background: For the gauge field Ψ y γ n iσ a U y,n Ψ y+n +Ψ y+n γ n U † y,n iσ a Ψ y , And for the two Higgs fields In the bosonic equations of motion, all fermion bilinears are replaced by quantum expectation values over the fermion ensemble fields,.

B Counterterms and Renormalization
The counterterms are chosen so as to compensate the backreaction terms of the fermion fields in vacuum, as these terms appear in the bosonic potentials and kinetic terms. Appearing as relevant or marginal operators, these counterterms include terms quadratically or logarithmically divergent in the continuum limit. Since we do not intend to take this limit in our simulations, the counter terms simply provide a subtraction of a finite (yet generically large) number. In practice, the most severe (quadratic) divergences are our main concern A.23.
When the Higgs fields are constant, as in the vacuum, the integration over the fermion modes yields where m D and m S are the masses of the doublets and singlets respectively, and K denotes the matrix Notice in the backreaction of fermion fields, there is a symmetry between Φ 1 and Φ 2 , owing to the fact that K is invariant under the exchange of Φ 1 and Φ 2 , and up and down simultaneously, while the latter exchange does not affect the determinant of the matrix. Therefore, the simplicity of the counterterms A.23 can be attributed to the symmetry in Yukawa interactions A.21.
To obtain the expression of ct 11 and ct 12 , we consider two cases: (a) Φ 2 = Φ 1 = φ y uk I, The series on the right-hand side contains a constant term which is the result of the free massless fermion, as one can see by simply setting φ = 0 on the left-hand side. Thus, instead of implementing the equation above, we consider its derivative where the constant term is removed,  where H refers to the Hamiltonian  Considering the continuum theory, one may immediately notice that ct 11 is quadratically divergent, as expected; while the divergence of ct 12 is logarithmic, or null for massless fermions. Conveniently, the solutions may be applied directly to Wilson fermion, in which case m D = 1 2 ar D (p 2 1 + p 2 2 + p 2 3 ) and m S = 1 2 ar S (p 2 1 + p 2 2 + p 2 3 ), with Wilson parameters r D and r S for doublets and singlets separately.