Electroweak Baryogenesis with Vector-like Leptons and Scalar Singlets

We investigate the viability of electroweak baryogenesis in a model with a first order electroweak phase transition induced by the addition of two gauge singlet scalars. A vector-like lepton doublet is introduced in order to provide CP violating interactions with the singlets and Standard Model leptons, and the asymmetry generation dynamics are examined using the vacuum expectation value insertion approximation. We find that such a model is readily capable of generating sufficient baryon asymmetry while satisfying electron electric dipole moment and collider phenomenology constraints.

1 Introduction Throughout the observable universe there is significantly more matter than antimatter, and establishing the dynamical origin of this asymmetry remains an open problem in physics. The baryon asymmetry of the universe, as determined by measurements of the cosmic microwave background made by the Planck experiment [1], indicates a baryon-to-entropy ratio Y B = n B s = (8.66 ± 0.04) · 10 −11 . (1.1) Any mechanism that might explain the origin of this asymmetry must satisfy the three Sakharov conditions [2]: There must be baryon number violating processes; violation of both charge conjugation (C) invariance and its combination with parity invariance (CP); and either out-of-equilibrium dynamics or violation of CPT conservation. A range of mechanisms have been proposed that can satisfy these criteria and thus account for the observed asymmetry. Electroweak baryogenesis (EWBG) is one such mechanism, in which the asymmetry is generated during a first-order electroweak phase transition (EWPT). The nucleation and expansion of the bubbles of broken electroweak symmetry during the early universe satisfy the out-of-equilibrium criterion, while CP-violating (CPV) interactions with the bubble wall catalyse asymmetry generation through electroweak sphaleron processes. In order for EWBG to successfully explain the observed asymmetry we require a strongly firstorder (SFO) electroweak phase transition, or SFOEWPT. However, the Standard Model (SM) instead features a crossover transition [3] which cannot provide the needed out of equilibrium dynamics. Additionally, even with an SFOEWPT, the amount of CPV present in the SM is not enough to yield the required asymmetry [4] Thus, the viability of electroweak baryogenesis requires the presence of Beyond the Standard Model (BSM) physics to generate both a strongly first order EWPT and to provide additional sources of CP-violation. Since electroweak symmetry breaking occurs at a temperature T ∼ 100 GeV, the corresponding BSM mass scale must not be too much higher than the electroweak scale. Such extensions to the SM should thus be testable at current or future colliders [5][6][7], or via high sensitivity fundamental symmetry tests such as electric dipole moment (EDM) searches, which place stringent constraints on additional sources of CPV [8].
While it is straightforward to obtain a strongly first-order electroweak phase transition by extending the SM scalar sector (the simple addition of a gauge singlet scalar often suffices [5,6,[9][10][11][12]), introducing new CP-violating phases in Higgs interactions while avoiding EDM constraints can be challenging. These constraints can be avoided in models where the EDMs arise at two-loop order, or where the CPV interactions are flavor non-diagonal, vector-like, or "partially secluded" from the SM 1 . In the latter instance, the introduction of additional scalars charged under SU (2) can result in electroweak symmetry breaking (EWSB) which occurs via multiple, successive phase transitions [13,14]. In these models EWBG may occur during the first EWSB transition, with a subsequent transition resulting to the usual SM Higgs-phase. This class of models have the benefit that the new CPV interactions can be hidden in the new scalar sector and hence avoid the EDM constraints [15].
One can also consider an EWBG scenario which combines elements of the above possibilities. In what follows, we consider a model involving two scalar singlets, S i , and a vector-like lepton doublet, ψ, with the same quantum numbers as the SM left handed lepton doublets. The vector-like nature of the new leptons allows for CP-violating Yukawainteractions with the gauge singlet scalars and SM leptons, while also avoiding any contri- 1 By partial seclusion we refer to scenarios where the CPV asymmetries are generated in species that carry no SM gauge quantum numbers (so there no EDM effects). These asymmetries get transferred to the SM via particle number changing reactions bution to anomalies. The scalar singlets catalyse a single-step strongly first order EWPT, while their interactions with the vector-like and SM leptons provide the necessary source of CP-violation. The EDM constraints for this scenario are somewhat relaxed as: • The electron EDM d e arises at two-loops due to restricting the CPV interactions to the third generation.
• The scalar singlets couple to the first generation leptons only via mixing with the SM Higgs boson (a form of "partial seclusion"), so the two-loop EDM involves a small mixing angle.
• Finally, the new CPV terms in the weak currents which contribute to EDMs in similar models [16] do not contribute due to the lack of light right handed neutrinos.
A model similar to the one outlined in this paper has recently been considered in Ref. [17]. In that work, a complex scalar singlet interacts with vector-like quark singlets, and CPV arises spontaneously from the vacuum expectation value (VEV) of the complex scalar singlet, which changes during the electroweak phase transition. This, and other similar models, often involve vector-like fermions which have large Yukawa couplings to the SM Higgs [18][19][20][21][22]. This is desirable in the context of electroweak baryogenesis since large Yukawa couplings can provide large CP violating interactions necessary for EWBG, and may also play a role in acquiring a strongly first order EWPT [18,19]. However, such models lead to significant corrections to Higgs boson properties, leading to important constraints from on the diphoton branching ratio and production processes.
In contrast, the scenario considered here easily evades these constraints as the SM Higgs coupling to the vector-like leptons is O(y τ ) (or even smaller) and plays no significant role in the asymmetry generation process. Both classes of models have to contend with the non-observation of direct vector-like lepton production at colliders, which implies a lower bound on their masses. This lower bound cannot be much larger than the temperature at which the electroweak phase transition takes place (T ∼ 100 GeV) in order for EWBG to remain viable. An in-depth examination of the collider signatures of a minimal vector-like lepton model has been undertaken in Ref. [23]. Using 8 TeV ATLAS data, Ref. [23] places a lower bound of ∼ 270 GeV on the masses of vector-like lepton doublets that decay via mixing with the SM leptons. However, the situation in our scenario is more complicated, as the existence of the new scalars modifies the decay chain of the vector-like leptons. We will also consider lepton universality, lepton flavor violation, and electron EDM constraints. We find that that such a model can readily generate the observed asymmetry while avoiding all current bounds, though future searches for vector-like leptons at colliders will begin placing severe constraints.
The layout of the remainder of this paper is as follows. In Section 2 we define our model along with a discussion of the interactions, mass mixing, and new CPV phases. Section 3 introduces some benchmark parameters and discusses various observables and constraints. In Section 4 we discuss EWBG methodology, how we derive and solve the transport equations, and our treatment of the EWPT. We conclude and discuss future directions in Section 5. Technical details relating to the electron EDM and transport equations are given in Appendices A and B, respectively.

Model
We extend the SM by introducing two real gauge singlet scalars S 1 and S 2 , and a vector-like lepton doublet ψ. These fields transform under the SM gauge group SU(3)×SU(2)×U(1) as (2.1) We use the notation where the s i and h are the dynamical fields, G 0 and G + are the Goldstone bosons, and v X is the VEV of particle X. As the ψ have the same quantum numbers as the SM leptons, they can in principle mix and interact with all of the SM lepton generations. Interactions with all of the lepton generations would face strict constraints from lepton universality, flavor violation, and electron EDMs. Furthermore, multi-generation interactions are not necessary for successful electroweak baryogenesis. Hence, for simplicity we only consider interactions between the ψ and the third-generation of SM leptons. The Lagrangian then includes the following Yukawa interactions terms, where L 3 = ν τ , τ − L T is the SM third generation left-handed lepton doublet. These Yukawa interactions result in mixing between the ψ and SM leptons, and provide the CPV necessary for EWBG. There are also two possible Dirac mass terms The M ψL mass mixing term can be removed via a redefinition of the ψ L and L 3 fields, If all of the Yukawa couplings are suitably redefined this transformation leaves the rest of the theory unchanged.

Scalar Potential
We now examine the scalar potential and outline how its parameters are related to the zero temperature VEVs and scalar masses. The scalar potential is where the sums go over the two scalar singlets, and the couplings are symmetric in all indices. As the S i are gauge singlets they can be shifted by a constant (S i → S i + c), and if the couplings are suitably redefined, this will leave the form of the Lagrangian unchanged. This shift freedom is frequently used to eliminate the linear terms a i S i . However, in order to simplify the resulting mass matrices that need to be diagonalised, we instead use this freedom to relocate the global minimum of the effective potential at zero temperature to a point where v S i = 0. This choice implies which we use to fix a i . We then require that the global minimum resembles includes a SM-like Higgs with mass and VEV, and µ 2 H such that we obtain the correct Higgs mass after diagonalisation of the mass matrix. Similar to the shift, we can also always perform a rotation of the singlets. We choose to rotate such that the scalar singlet mass mixing terms disappear, which fixes a 12 such that The remainder of the scalar potential couplings are then free parameters, subject to the requirement that the SM Higgs-phase is indeed the global minimum of the effective potential at zero temperature. We are interested in the scenario where this potential enables baryogenesis by having a strongly first-order phase transition during which both of the scalar singlets have changing VEVs. In this paper we focus on the transport and asymmetry generation dynamics and do not perform a thorough scan of the scalar potential parameter space. We utilise a single EWPT benchmark point which fixes these scalar potential parameters and satisfies the SM-like global minimum requirement (see Section 4.2). A more detailed study of the phase transition is left to future study.

New CPV phases
We now briefly discuss the new CPV phases that appear in our model and how we parameterise them. Some of the phases in the Yukawa couplings can be removed by rephasing fields via the transformation (2.10) With these rephasings we find that the Lagrangian is unchanged if we take such that we have the following invariant phases, (2.12) We will find that δ 1 appears within the EWBG calculation as the primary source of the CP-violation necessary to generate the asymmetry. The other phases will contribute to electron EDM constraints but not significantly impact baryogenesis. From here on we choose to re-phase our couplings such that the following conditions are satisfied which is equivalent to setting (2.14) The choice of this specific phase will simplify some of the expressions in the following sections.

Mass and Mixing Matrices
We now consider the diagonalisation of the mass matrices. The neutral scalar and lepton mass terms are given by, These mass matrices can be diagonalised by a redefinition of the fields such that The M N matrix is already diagonal and real, such that no redefinitions are necessary. However, the form of the weak currents after diagonalisation can be simplified by choosing to re-phase the neutral vector-likes N L/R = N L/R e iδ 4 . In the limit where the scalar mass mixing terms v H b i are small relative to the scalar mass differences, the elements of P are well-approximated by, For the fermions, in the case where M ψ in M E is much larger than any other term, we can approximate the mixing matrices via and δ 4 is given in eq. (2.14). The τ mass defined in eq. (2.17b) will be shifted slightly from its SM relation with Higgs Y τ Yukawa and is given by, (2.21) As we want the physically measured m τ , we fix Y τ as a function of M ψ and Y ψ to give the correct value. If |Y ψ | ∼ Y τ and M ψ > 200 GeV the correction to Y τ will be orders of magnitude smaller than the current uncertainty in Y τ inferred from Higgs decay measurements [24]. The mass of the charged component of ψ will also be shifted slightly, This will lead to a small mass splitting between the neutral and charged vector-like leptons which will be further enhanced by radiative corrections; however, we will still generally refer to the masses of the ψ components as simply M ψ . Using these mixing matrices one is able to write out the new forms of the weak current and Yukawa interactions after diagonalisation. These are given in Appendix A, where we have employed the notation of Ref. [25], whose electron EDM expressions will be used in the next section. As the mixing angles are small, the diagonalisation matrices were all selected such that the mass basis states are labelled by their primary components, e.g., τ will be the physical lepton that consists primarily of τ L and τ R . Throughout the phenomenology section we will refer to the mass basis and simply drop the primes.

Phenomenology
In order to quantitatively investigate the phenomenology of this model we consider the four benchmark scenarios defined in Tables 1 and 2. The resulting mixing parameters are provided in Table 3. The scalar potential parameters of benchmark points A and C were obtained from a scan for a suitable electroweak phase transition, as discussed in Section 4. Benchmark B is a variant of A with slightly different Yukawa couplings which modify the decays of the new particles leading to significant effects on the collider phenomenology. Benchmark D is a variant of C with no mixing between the physical Higgs and the new scalars. We will study the constraints from Higgs physics, direct collider searches sensitive to the vector-like leptons, and precision lepton phenomenology.

Vector-Like Lepton Collider Phenomenology
The charged and neutral components of the vector-like lepton doublet can be produced at colliders via Drell-Yan processes. In minimal vector-like lepton models they decay into a SM lepton (a τ or ν τ ) plus either a SM Higgs or weak gauge boson via the mixing induced by the SM Higgs Yukawa Y ψ . This results in a lower bound on minimal vector-like lepton doublet masses M ψ 270 GeV [23] (using 8 TeV data). However, in our scenario if the gauge singlet scalars are lighter than the vector-like leptons, the leptons may instead decay via the singlet scalars. This is the case for all of our benchmark points. This is motivated by the requirement that both of the scalars have changing VEVs during the electroweak phase transition, which generally requires them to have non-negligible negative mass terms. When combined with perturbativity bounds on the b ii couplings, this leads to an upper bound on the zero temperature scalar masses m 2 We thus expect that if we have the desired phase transition the scalars will generally be lighter than the vector-like leptons. 0.511 0.962 b 22 1.36 0.036 0.05 Table 1: Benchmark point parameters. Parameters not listed here are the same for all the benchmark points, and are given in Table 2.  Table 3: Mixing parameters for the benchmark points as defined in Section 2.3. Benchmark points A and B have the same mixing angles. θ L,R are the same across all the benchmark points as they only depend on M ψ and Y ψ .
The dominant decay modes of the vector-like leptons are depend on the relative sizes of the Yukuwa couplings with the scalar singlets λ Lψi , and the Yukawa coupling with the SM Higgs Y ψ ,and the mixing angles θ L and θ R . If Y ψ dominates, then the main decay modes of the new leptons will be as for the minimal vector-like lepton doublet model. On the other hand, if Y ψ and the associated mixings are small relative to λ Lψi then ψ will mainly decay into a SM lepton and a singlet s i . In the benchmark points we consider this is the case, as can be seen in Table 4, which shows the main branching ratios of the vector-like leptons and scalar singlets for benchmarks A and D.
The decays of the scalar singlets depend on whether the singlets mix with the SM Higgs boson. If they do, then they inherit their decays from the Higgs decay channels, with the Higgs mass set to that of the appropriate singlet. This is the case for benchmark points A, B and C . Since the scalar singlets for these points have masses between 127 and 200 GeV their decays are mainly into bb, W W ( * ) and ZZ ( * ) . For benchmark point D the singlet mixing with the SM Higgs is turned off and the singlets decay into τ + τ − pairs through the λ Lψi coupling and E-τ mixing. We have checked that singlet lifetimes in all cases do not lead to displaced vertices at the LHC or constraints from BBN. Therefore, in E + E − pair production the final state will always contain a τ + τ − pair and the decay products of two singlets. As the singlets frequently decay into W W ( * ) /ZZ ( * ) , which may in-turn decay leptonically, we expect multilepton searches featuring signal regions with multiple τ leptons to be the most sensitive to our scenario.
To examine in detail the constraints arising from collider analyses we have implemented our model in FeynRules 2.3.29 [26,27]. Using the UFO format [28] we import our model into MadGraph5_aMC@NLO 2.6.1 [29,30] in order to generate events at parton level, before showering them using Pythia 8.230 [31,32]. These are then fed into CheckMATE 2.0.26 [33][34][35][36][37], which includes an interface with Delphes 3.4.1 [38], in order to determine the most relevant ATLAS or CMS analyses and constraints on our model. QCD corrections to the pair production of SU(2) triplet leptons in type-III seesaw models at 13 TeV leads to K-factors ranging from 1.17-1.25 for lepton masses ranging from 100-1000 GeV [39]. Based on this result we scale the ψ production cross sections passed to CheckMATE using a K-factor of 1.2.
For the benchmark points listed in Table 1 we varied M ψ from 200 to 980 GeV in steps of 20 GeV and generated two million events of pp → EE, EN , N N for each parameter point. CheckMATE evaluates the expected number of signal events S and associated error ∆S per signal region for each of the implemented analyses and then evaluates the ratio of the 90% confidence lower limit on the number of expected signal events S − 1.64∆S and the 95% confidence upper limit on experimentally observed number of events S 95 . This ratio, provides a conservative quantification of the constraints on given model, with any point having r > 1 considered excluded [33]. The resulting r-values are shown in Fig. 1. As expected, the most important analyses are generally ATLAS or CMS searches for charginos and neutralinos with multilepton final states, specifically the searches in Refs. [40][41][42] which involve 36fb −1 of data taken at 13 TeV. The exclusion limits on benchmark D, where the singlets have no mixing with the SM Higgs boson, are noticeably stronger than the rest. In this case the singlets always decay to τ + τ − , so that each signal event involves 4 − 6 τ leptons. This is strongly constrained by signal regions of the above CMS and ATLAS analyses involving hadronic τ leptons, missing energy, and either 0 or 1 pairs of opposite sign-same flavour leptons.
On the other hand, when mixing with the SM Higgs is present, the ψ mass can be brought as low as 350-500 GeV. The lowest mass bound of the considered benchmark Benchmark point A Particle Total Width Decay Products Branching Fraction Total Width decay products Branching fraction Table 4: List of particle widths, decay modes, and branching fractions as calculated by MadGraph5 (via MadWidth) for benchmark points A and D. Benchmark points A and D were selected in order to illustrate the effect of s i -h mixing, which is not present in benchmark point D.
points is given by benchmark point A. The relative sizes of the Yukawa couplings in benchmark A cause the ψ to decay primarily to just the lightest scalar singlet s 2 . This decays mostly to bb, leading to more jets than leptons and thus weaker constraints on the vector-like lepton masses.
We have also calculated the production cross-sections for scalar singlets, which may in principle be constrained by Higgs boson searches and measurements. We find that Higgs searches have no impact for our benchmark points, due to the very small mixing between the scalars and the SM Higgs, as given in Table 3. The singlet production cross-sections in standard Higgs production channels are ∼ 10 −3 fb for benchmarks A and B, and hence unobservable even at the HL-LHC. The cross section for benchmark C is larger at ∼ 0.1fb, but this leads to no constraints from currently available searches. In principle the scalar mixing could be larger, which could lead to significant constraints that require a more indepth analysis. See, e.g., Refs. [43,44] and references therein. However, a larger mixing does not reduce the lower bound on the ψ masses and is thus uninteresting in the context of the viability for this model to generate a baryon asymmetry. We have used small mixing terms for simplicity.
Singlet pair production is also possible. In the SM, Higgs pair production is due to a triangle and box diagram, which interfere destructively. Singlet pair production via an intermediate off-shell SM Higgs (via the triangle diagram) is not suppressed by the singlet-Higgs mixing angle. However, the box diagram is suppressed by two powers of the mixing angle and is thus negligible in our model. We therefore expect this leads to an O(1) in the dihiggs production cross-section at the LHC. Given that this process is not expected to be observed until the end of the high-luminosity LHC run, this is unconstrained by current LHC searches [45]. We refer the reader to Refs. [43,44,[46][47][48] for discussions on singlet pair production at current and future colliders.
Finally, the vector-like leptons also lead to a correction to the Higgs diphoton decay rate. However, as we consider the scenario where the ψ Higgs Yukawa is of the order of the τ Yukawa, the correction is negligible compared to the SM contributions and current experimental accuracy.
We conclude that scalar mixing is critical for a low mass bound and that M ψ 500 GeV is likely a strong enough requirement to satisfy current collider bounds, though the specific lower bound will depend on λ Lψ1 λ Lψ2 and may in some cases be larger. Figure 2: Electron EDM contribution arising from new interactions. CPV comes from either phases in the weak currents or phases in the new scalar-fermion interactions which arise through scalar mixing. Mirrored versions of the second diagram (i.e., Z ↔ φ i , E − ↔ E + ) also contribute.

Electron Electric Dipole Moment
The electron EDM d e provides a strong constraint on new CP violation. For large phases this will lead to constraints of the sizes of the Yukawa couplings, scalar mixing, and ψ-L mixing. Current experimental measurements place the upper bound on the electron EDM [49] at, We consider contributions to the electron EDM via the two two-loop Barr-Zee style diagrams shown in Fig. 2. The relevant EDM formulae used to obtain numerical results are provided in Appendix A.
The diagram in Fig. 2a could produce EDM contributions due to new CPV interactions in the charged weak currents, as is the case in vector-like quark models [16]. This diagram relies on a relative phase between the left-and right-handed charged weak current inter- where M L,R are the left-or right-handed weak current mixing matrices that are given in Appendix A.1. However, in our scenario, these matrices are real, such that there is no EDM contribution from this diagram. This can can be attributed to the lack of a right handed ν τ for the N R to mix with. We do, however, obtain nonzero EDM contributions driven by scalar mixing, as in diagram in Fig. 2b. These contributions are similar to those that appear in supersymmetric models with chargino and neutralino loops and thus we make use of the formulae provided in Ref. [25] to evaluate these contributions.
Substituting the scalar and ψ masses into the EDM loop functions yields expressions for the EDM contributions as a function of phases, Yukawa interactions, and mixing angles, though the latter are not independent of the masses. Doing so for benchmark A or B and retaining only terms first order in mixing angles we obtain  where P ij are the components of the scalar mixing matrix defined in Section 2.3, with numerical values shown in Table 3. The predicted values are O(10 −2−3 ) less than the current bound, but could be probed at future EDM experiments. Full numerical electron EDM results for benchmarks A-C are given in Table 5. The EDM is much larger in benchmark C as it has significantly larger b i couplings which lead to scalar mixing. As benchmark D does not feature scalar mixing there is no two-loop electron EDM contribution. The contribution of the CPV phases δ 1 and δ 4 to the EDM are suppressed relative to the δ 2 and δ 3 contributions by an additional factor of the L-ψ mixing angles θ L/R . This is because the diagram shown in Fig. 2b that involves the Yukawa couplings λ Lψi and Y ψ associated with these phases requires non-zero L ↔ ψ flavor changing neutral currents. These are proportional to the mixing angles θ L/R . This is not the case for the λ ψψi couplings associated with the δ 1 and δ 4 phases.
The suppression of the λ Lψi contribution is significant as it is these couplings that are the primary source of CP-violation for baryogenesis. Hence, while greater precision in electron EDM measurements will place constraints on the scalar mixing matrix P and λ ψψi couplings, this will not directly constrain the parameters most relevant to EWBG. Furthermore, as the suppressing mixing angles are directly dependent on the Yukawa coupling Y ψ as shown in eq. (2.20), and as Y ψ does not play a significant role in EWBG, we can further seclude the CP-violation relevant to baryogenesis by taking Y ψ =0. We would then require at least 3-loop diagrams to contribute to the electron EDM.

Lepton Phenomenology
We will briefly discuss corrections to lepton universality and Z/W couplings. Because the new particles couple only with third generation leptons, we do not obtain a contribution to lepton flavour violating decays. However, the new particles modify the couplings of the weak gauge bosons to τ and ν τ . The weak currents and Yukawa interactions are provided in Appendix A.1. Due to the mixing between τ and E the Zτ R τ R and W ± τ L ν τ interactions are modified, where c w and s w are the sine and cosine of the weak mixing angle, and P L,R are the standard left-and right-projection operators. There is no correction to the Zτ L τ L vertex. The corrections to the weak-current couplings lead to small deviations from lepton universality.
As an example, the τ decays will be modified due to the non-unitarity of the standard 3 × 3 Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix U PMNS . The τ → νν decay rate Γ τ is proportional to where the unitarity of the PMNS matrix is recovered in the θ L → 0 limit. Note that these universality-violating effects are second order in the mixing angles θ L and θ R . Using the mixing angles in the benchmark points listed in Table 3, the corresponding universalitybreaking corrections will be of O 10 −7 . There will also be a contribution to the 1-loop Zτ τ coupling from vertex corrections of the types shown in Fig. 3. While superficially divergent, the divergences cancel after including diagrams with Goldstone bosons and accounting for the τ wave-function renormalisation. Evaluating the relevant loop diagrams using Package-X [50] and inserting the benchmark couplings and mixing angles, we find that this vertex correction is ∼ 0.01 This correction is many orders of magnitude smaller than those found in vector-like quark models [51], as the latter feature factors of m t /M ψ rather than m τ /M ψ , and generally consider much larger Yukawa couplings than those considered here. Corrections to the weak gauge couplings of this order of magnitude are easily small enough to avoid existing constraints on lepton universality and Z-coupling measurements [52,53].

Electroweak Phase Transition and Baryogenesis
This section is separated into three parts: an overview of EWBG, a discussion about the EWPT and bubble nucleation, and a discussion of the transport dynamics around the moving bubble wall. Our methodology for EWBG closely follows that of Refs. [15,54].

Overview of EWBG
Electroweak sphalerons are anomalous processes that violate B + L while preserving B − L. Prior to electroweak symmetry breaking, the rate for these processes is relatively rapid, being proportional to the temperature T . The presence of a non-zero chemical potential, µ L , for any fermions carrying SU(2) L charge can act to bias the electroweak sphalerons to produce a net B+L charge density. If the EWSB transition is sufficiently out-of-equilibrium, the net B + L charge will be preserved in regions of broken electroweak symmetry before minimisation of Gibbs free energy drives B + L → 0.
A first order electroweak phase transition can provide the necessary conditions for this B + L production and preservation if the effects of CP-violation generate a sufficiently large left-handed fermion chemical potential and if the degree of sphaleron "quenching" is sufficiently strong in regions of broken electroweak symmetry. The transition proceeds via the nucleation and expansion of the bubbles of broken electroweak symmetry. In the present scenario, the VEVs of the SM Higgs doublet and the scalar singlets evolve, varying across the bubble walls. The corresponding CP-violating Yukawa interactions generate a non-zero µ L , as described in detail below. This left-handed charge diffuses ahead of the expanding bubbles, thereby biasing the electroweak sphalerons into net B + L generation. The expanding bubbles capture the resulting asymmetry, preserving it if the sphaleron transitions are sufficiently quenched inside the bubbles.
The degree of preservation is governed by the broken phase EW sphaleron rate Γ WS that is exponentially suppressed as where E WS is the weak sphaleron energy and where the prefactor A is a function of T . The value of E WS can be related to a T -dependent energy scalev(T ) associated with electroweak symmetry breaking. For a discussion of the relationship betweenv(T ) and the scalar field VEV and the associated issue of gauge invariance, see Ref. [55]. It is conventional to write where g is the SU(2) L gauge coupling and B is a computable function of g and the other couplings in the gauge-Higgs sector. Note that both A and B will vary depending on the representation of the EW symmetry breaking scalar [56][57][58]. Preservation of the B + L asymmetry requires a sufficiently large E WS /T at the the bubble nucleation temperature T N , which typically lies just below the critical temperature of the transition T C . A firstorder phase transition that satisfies this requirement is dubbed "strong", a characterisation that is usually translated into the requirement that For recent discussions of sphaleron rate computations in perturbation theory and the associated theoretical uncertainties see, e.g., Refs. [55,[57][58][59]. The dynamics of the expanding bubbles, together with the CPV transport dynamics, constitute the crucial elements for generating the left-handed number density that catalyses B + L generation. In what follows, we concentrate on the transport dynamics but note here the importance of the bubble expansion rate, characterised by the wall velocity w. On the one hand, one must have w > 0 in order to generate any asymmetry (see below). On the other hand, the expansion must be sufficiently slow to allow the left-handed number density to diffuse ahead of the advancing wall and "seed" the EW sphalerons before they are quenched in the bubble interior. For a detailed discussion of the bubble dynamics and wall velocity, we refer the reader to Refs. [60][61][62].
For purposes of analysing the CPV transport dynamics, it is conventional to treat the bubble wall as a flat plane moving with a constant w. Under this assumption, the primary inputs needed are the specific CPV interactions, T , w, and the bubble wall profile as a function of the direction normal to the wall. As a result, the procedure of evaluating the asymmetry generated by a given model during the EWPT can generally be split into two steps: verifying the presence of a strongly first-order electroweak phase transition and obtaining the moving bubble wall information, and then solving a system of outof-equilibrium transport equations to compute the final baryon asymmetry. For a more thorough review of EWBG see [63,64].

Electroweak Phase Transition
In this section we outline our treatment of the electroweak and the selection of the benchmark points listed in Section 3. Our goal is not to perform a comprehensive study of the model parameter space that yields a strongly first-order EWPT but rather to identify illustrative parameter choices for use in the treatment of the transport dynamics.
To that end, we will employ the high-temperature effective potential which is acquired by adding thermal mass terms to the tree level potential, where V 0 is the tree level potential given in eq. (2.6), and the thermal masses δm 2 X are computable in terms of the benchmark model parameters and are listed in Table 9 within Appendix B.1. This treatment has the advantage of being directly gauge-independent, numerically fast to evaluate, and still yields results comparable to a more thorough gaugeindependent treatment, as seen in Ref. [59]. Note that, in this context, the scalar field VEVs are manifestly gauge invariant, and we will identify these VEVs with the corresponding scalē v introduced above.
In order to find benchmark points that yield the desired phase transition, we performed a random scan over the scalar potential parameter space and utilised a modified version of the CosmoTransitions package [65] to determine T C , the nucleation bubble wall profiles, and the bubble nucleation probability density per unit time [66], where S is the three-dimensional Euclidean bubble action. The expected number of bubbles in a Hubble volume is then given by [67], The corresponding nucleation temperature is defined by the criterion where g is the effective number of relativistic degrees of freedom and m pl is the Planck mass. This T N is then used throughout the remainder of the EWBG calculation.
The results of using CosmoTransitions to trace the minima of the effective potential for benchmark point A are shown in Fig. 4. As the temperature decreases benchmark point A exhibits two second-order phase transitions where the scalar singlets gain their VEVs. Eventually there is a strongly first-order electroweak where the Higgs gains a VEV and the scalar singlets lose theirs. Such a phase transition will result in a potentially observable gravitational wave signal [68][69][70][71][72]. It should be noted that if the preceding phase transitions where the scalar singlets gained their VEVs were also first order, these transitions would also produce gravitational waves. When combined with the gravitational waves from the Higgs phase transition this would lead to a multi-peaked gravitational wave power spectrum that is characteristic of multi-step phase transitions [67,73].
Determining the wall velocity and moving bubble wall profiles is a complicated task requiring solving a set of integro-differential equations [60,61]. However both nucleation and moving wall profiles are generally well approximated by hyperbolic tangent profiles with a characteristic wall width L, where z is the distance from the planar bubble wall and z > 0 corresponds to the inside of the bubble where the SM Higgs has a non-zero VEV. In the limit where the wall width is small relative to the ratio of the diffusion scale D and reaction rates Γ in the thermal bath L D/Γ, the generated asymmetry is not strongly-dependent on the wall width. The reaction rates and CPV source term can then be well approximated by step and delta functions, respectively. However, as we discuss in more detail in Section 4.3, the magnitude of the CPV source term is sensitive to the relative magnitudes of the various VEVs within the bubble wall, and thus an accurate expression for the moving wall profiles is still necessary. For simplicity, we fit tanh-functions to the nucleation profiles and scale the profiles such that the wall width of the Higgs VEV is representative of values typically found in studies of the moving wall, L H = 10/T [60]. These studies find wall velocities w ∼ 0.1-0.3. For our benchmark we take w = 0.1. Recent progress has been made towards generalising moving wall dynamics for minimal extensions of the SM [61,62]. However, the results of Ref. [62] do not apply straightforwardly to models with multiple scalars, so we proceed with the approach outlined above and leave a detailed study to further work.
The radial nucleation bubble wall profiles obtained by CosmoTransitions for benchmark A are shown in Fig. 5. The rescaled parameters for the profile function defined in eq. (4.10) are given in Table 6 along with the nucleation temperature. The differences in the phase transitions between benchmarks A and B are negligible, and similarly for C and D. Hence, only results for A and C are given in the table. In all cases, it is essential that the profiles of the two singlet scalars differ from each other in the region where one or both are varying.  Table 6: EWPT Temperature and rescaled VEV-profile parameters for eq. (4.10). The parameters for benchmarks B and D are approximately equivalent to those of A and C, respectively.

Quantum Transport Equations
The baryon asymmetry is generated by electroweak sphalerons acting on the net density of fermions charged under SU(2) that are diffusing ahead of the bubble wall. Computing this number density requires solving a set of transport equations of the form where J λ i and µ i are the number density current and chemical potential of particle species i. The Γ ij... are equilibration rates that arise due to various interactions in the plasma, while S CPV i is a CPV source term arising from interactions with the bubble wall. The number density and chemical potential for species i are related as where n f (ω, µ) is the fermion distribution function and g i is the number of degrees of freedom associated with the particle species. Treating the bubble as a flat plane moving with speed w, employing the diffusion approximation J = −D ∇n, and using eq. (4.12), we can re-write eq. (4.11) as We use the diffusion constants D i derived by Ref. [74], which are provided in Appendix B.1 along with the thermal masses and widths. The fermion diffusion constants were derived for massless particles and arise from their gauge interactions. Even though the ψ are massive, for now we simply assume D ψ ≈ D L i and will briefly investigate the dependence of the final asymmetry on D ψ . Calculating the rates Γ ijk and CPV source terms S CPV requires the use of out-ofequilibrium finite-temperature field theory. We will use rates and source terms calculated using the VEV-insertion approximation (VIA) [54,75,76], which itself relies upon the Schwinger-Keldysh closed time path framework. The VIA formalism utilises self-energy diagrams, such has the one shown in Fig. 6, to compute the reaction rates and source terms for the species in the thermal bath through the use of VEV-insertions. A more thorough treatment necessitates a VEV-resummation, resulting in a spatially varying mass matrix and flavour oscillations. The VIA formalism may overestimate the generated asymmetry, and the results ought to be considered illustrative of what may be expected from such a model. For a more thorough review of the theoretical issues see Ref. [64]. We do not give detailed derivations of the rates and source terms in this paper. Instead we provide an overview of rates involved along with the full system of transport equations, and then give analytic formulae and numerical results for the rates and source terms in Appendix B.3.
When performing the VIA derivations, the type of diagram shown in Fig. 6a gives rise to two equilibration rates and a CPV source term. These equilibration rates, which we denote as Γ ± AB S i ,H (µ A ± µ B ), are proportional to the product of the VEVs of S i or H, and hence have a spatial dependence across the bubble wall. The CPV source term arising from the VEV insertion, denoted S CP V AB , is non-zero only within the bubble wall when the VEVs are changing. The absorptive parts of Yukawa-loops, shown in Fig. 6b, give rise to equilibration rates, Γ AB S i ,H (µ A − µ B ± µ S i ,H ). These rates also acquire some spatial dependence as the masses of the particles in the diagram will vary across the bubble wall. Additionally, strong and EW sphalerons will also act on the number densities. The effect of the sphalerons is to introduce terms proportional to where Γ WS is the EW sphaleron rate, Γ ss is the strong sphaleron rate, u i and d i are the three generations of right-handed SU(2) singlet quarks, Q i are the three generations of lefthanded quark doublets, and L i are the three generations of left-handed SM lepton doublets.
We approximate the sphaleron rates as [77][78][79], with κ WS ≈ 20, κ SS ≈ 14, and Γ WS exponentially suppressed within the bubble such that it is effectively zero as far as the transport dynamics are concerned. In summary we have have a combination of VEV-insertion, Yukawa-loop, and sphaleron equilibration rates driving the number densities back to zero while a CPV term, which is non-zero only within bubble wall, provides the source for the number densities that diffuse into and ahead of the bubble. EW sphalerons acting outside of the bubble may generate a non-zero net baryon density. We can simplify our system of transport equations by reducing the number of species to be considered via some equilibrium considerations. Due to rapid weak interactions we take the components of the SU(2) doublets to be in equilibrium, such that for the third generation leptons we have with similar relations holding for other doublets. Figure 6: Self-energy diagrams that contribute terms to the transport equations, including; (a) VEV-Insertions with spatially varying VEVs that will provide a CPV source term S CP V AB and an equilibration rate Γ ± AB S i ,H acting on the fermions, and (b) Yukawa-loops which just provide an equilibration rate Γ AB S i ,H acting on all of the particles in the diagram.
At this order our treatment will only provide a CPV source term for the SM third generation leptons and ψ via their scalar singlet Yukawas. Non-zero densities for the other particles are introduced only via EW sphalerons, or via the τ or ψ Yukawa couplings to the Higgs. As the τ Yukawa is relatively small, and as the benchmarks under consideration have an even smaller ψ Yukawa, we can simplify our system of equations by neglecting these Yukawa couplings as their resulting reaction rates are 2-3 orders of magnitudes smaller than the relaxation rates arising from the λ Lψ couplings. Then as the strong sphaleron rate equilibrates left and right handed quarks faster than the weak sphalerons produce them, and as the sphalerons act equally on all generations, we can use It should be noted that while the top quark Yukawa interactions are fast, they drive µ Q 3 − µ d 3 − µ H → 0. Hence they will not change these relations if there is no source for µ H . The first part of eq. (4.14) then reduces to When Γ WS is much smaller than all other relevant reaction rates, it is reasonable to decouple the weak sphaleron rate from the system of transport equations and compute the resulting baryon asymmetry in two steps: first, compute the left-handed chiral charge; second, use the latter to compute the B + L asymmetry from the sphaleron rate equation [15,54,80].
In the present case, however, we will find that the weak sphaleron rate (∼ 5 · 10 −4 GeV) may be comparable to the L 3 and ψ relaxation rates (∼ 10 −3 GeV in benchmark A). Thus we will instead include the weak sphaleron effects into our system of transport equations. Using the previous approximations, we obtain the following set of coupled transport equations: The relevant equilibration rates and CPV source terms are given by, The Λ 0,± L 3 ψ are numerically evaluated integrals that appear when evaluating the rates arising from the VEV-insertion diagrams. Similarly the I F function is a numerically evaluated integral arising from the Yukawa-loop processes that are dependant on the thermal masses of the particles in the loop. The thermal masses are given in Appendix B.1, and the formulae for Λ ±,0 and I F are provided in Appendix B.3 along with numerical values for the reaction rates in benchmark A. We also consider a more complex system of transport equations, outlined in Appendix B.2, which is suitable when the ψ Yukawa is large and we can no longer use some of the earlier simplifying assumptions.
Note that although we have four CPV phases δ 1 -δ 4 , we only have one CPV term which is proportional to Im λ Lψ1 λ * Lψ2 = |λ Lψ1 λ Lψ2 | sin δ 1 . (4.21) The δ 4 phase does not contribute since introducing its associated Yukawa coupling Y ψ , needs a Higgs-VEV-insertion. The only way to construct self-energy diagrams of the type seen in Fig. 6a is to have two such insertions, leading to CPV sources proportional to terms like Im Y ψ Y * ψ = 0. Regarding the absence of the δ 2 and δ 3 phases, note that if the vector-like leptons had a Dirac mass that was much smaller than the thermal masses, then we could treat ψ L and ψ R as independent degrees of freedom [81]. Then the δ 2 and δ 3 phases would source an asymmetry in n ψ L and n ψ R . However for our large Dirac mass these sources have have no significant effect on our transport equations.
This system of four linearly independent ODEs are solved numerically using a technique similar to the one outlined in appendix C of Ref. [15]. The approach outlined there treats the VEV and mass dependent rates as a step function across the bubble wall. We go one step further and also treat the source term as a delta function, which introduces negligible error for the bubble wall lengths that we consider. A more thorough treatment only becomes important once the scale of the wall width L approaches that of the typical diffusion-reaction distance, when we no longer have L D/Γ. Once the transport equations are solved,  Table 7: Baryon asymmetry generated by each benchmark point using the full and approximate system of transport equations.

EWBG Results
The resulting baryon asymmetries generated by each of the benchmark points are listed in Table 7, for both the simplified and full system of transport equations as given in Section 4.3 and Appendix B.2, respectively. The number densities obtained by solving the simplified system of transport equations for benchmark A are shown in Fig. 7, where the bubble wall is located at z = 0, with z > 0 corresponding to the interior of the bubble. The number density profiles obtained for benchmark A are illustrative of those obtained for the other benchmarks. The effects of varying some of benchmark A's parameters (λ Lψi , M ψ and Y ψ ), the bubble velocity w and the diffusion rate D ψ are illustrated in Fig. 8. As the simplifying assumptions no longer hold when Y ψ becomes larger, we have used the full system of transport equations for the contour plots. As the source term is proportional to |λ Lψ1 λ Lψ2 |, we would expect the contours of constant baryon asymmetry to be hyperbolic in the λ Lψ1 -λ Lψ2 plane, which is indeed seen in Fig. 8a. The small deviations from this relationship arise due to the equilibration rates arising from the Yukawa-loop self-energies. In general the resulting asymmetry can vary by several orders of magnitude depending on the sizes of these Yukawa couplings. However, as discussed in Section 3.1, the relative size of these couplings can have a significant effect on the decay of the ψ, and thus has a non-trivial impact on the collider phenomenology. Figure 8b shows the dependence of the final asymmetry on the bubble wall velocity. A thorough treatment of the phase transition dynamics is vital for a reliable EWBG calculation. However, there is significant freedom in the size of the Yukawa interactions to boost the amount of asymmetry generated. Hence, assuming the bubbles are not relativistic, EWBG ought to remain viable even for fast-wall situations. Figure 8c shows that variations in D ψ or Y ψ results in a relatively small change in the final asymmetry. The error introduced by utilising the massless fermion derivation of the diffusion constant (D ψ = D L ) is thus likely negligible compared to other factors affecting the asymmetry generation. The asymmetry is much more strongly dependent on D L as the left-handed leptons diffusing ahead of the wall are what directly bias the sphalerons, and D ψ only acts to effectively modify the rate with which n L 3 equilibrates back to zero.

Conclusion
We have investigated the capability of vector-like leptons and scalar singlets to generate the observed baryon asymmetry, with an emphasis on the phenomenology of such a model. The singlets are necessary to induce a strongly first order electroweak phase transition. The vector-like leptons introduce the CP-violating interactions with the singlets and SM leptons needed to generate asymmetry during the electroweak phase transition. Unlike some other electroweak baryogenesis models featuring vector-like fermions [18,19,21], this model shares a trait with multi-step EWPT models [15] in that the necessary CPV does not come from interactions with the SM Higgs but from the new scalar content.
From collider constraints we found that the presence of the new scalars will generally lead to a larger mass bound than minimal vector-like lepton models unless mixing is introduced in the scalar sector. The scalar mixing then provides the primary contribution to the electron EDM. However, the generated electron EDM is still a couple of orders of magnitude below current experimental bounds. Additionally, the CPV phase which enables EWBG contributes only at second order in the mixing angles. The capability for baryon asymmetry generation may begin to be significantly restricted if the vector-like lepton mass bounds from collider phenomenology continue to increase.
The model examined in this paper, which is complementary to the one in Ref. [17], indicates that scalar singlet plus vector-like fermion models can readily generate the observed baryon asymmetry during the EWPT. There are a number of areas in which further work could be pursued, including more thorough investigations of the moving wall dynamics, going beyond the high-temperature approximation for a more rigorous treatment of the EWPT, and resolution of uncertainties in the accuracy of the VEV-insertion approximation formalism. A better understanding of these issues would lead to a more accurate determination of the allowed parameter space of our model.

A.1 Mass Basis Weak Currents and Yukawa Interactions
As mentioned in the main text, the notation we use for the weak currents and Yukawa interactions are selected to match the notation of Ref. [25], such that it is straightforward to apply their general electron EDM formulae if one makes the following replacements, We define a matrix M such that after moving to the mass basis the charged weak currents become, where P L,R are the standard left-or right-projection operators. Similarly for the neutral currents we define a matrix G such that, For the Yukawa interactions we define matrices D i such that, Using the mixing angles and phases defined in Section 2, these matrices are given by, In the expression for D i provided above we have dropped any terms second order in mixing angles P i,j =i and θ L/R , though the numerical calculations included these terms.

A.2 EDM Formulae
Using the notation introduced in Section A.1 and the results from Refs. [16,25,83], the primary contributions to the electron EDM are given by Where we have used, Substituting the matrices defined in eq. (A.5) one finds that in our case d W − W − e = 0. To first order in the mixing angles the only nonzero contributions to the EDM come from terms proportional to (D L,R φ i ) 11 , which correspond to Barr-Zee style diagrams involving an E + -E − loop.

B Thermal Properties and Transport Equation Functions
B.1 Thermal Properties The thermal mass formulae and diffusion constants used are given in Tables 9 and 8 respectively. For massless fermions the thermal widths at zero momentum are given by [84] Particle Q, d, u L τ H D i /T 6 100 380 110 Table 8: Diffusion constants for the transport equations as derived in Ref. [74].
Particle   where C F denotes the quadratic Casimir invariant. In the limit of masses heavy compared to the temperature one instead has, which differs from the massless case by a factor of a half. The resulting widths for the relevant particles are given in Table 10.

B.2 Full System of Transport Equations
If the τ and ψ Yukawa couplings are not negligible, the system of transport equations becomes significantly more complicated. We must then consider n H , n Q 3 , n u 3 in more detail. Neglecting the Yukawa interactions of the lighter fermion, instead of eqs. (4.17a) we will get µ L 1 = µ L 2 , µ Q 1 = µ Q 2 , µ u 1 = µ u 2 = µ d i . (B. 3) The full system of quantum transport equations to be solved is then,

B.3 Transport Equation Formulae
The relaxation rates and source terms used in eqs. (4.20) and (B.5) were derived using the VEV insertion approximation as outlined in Ref. [54]. For the SM fermions, which are massless prior to EW symmetry breaking, we make use of the interacting thermal propagator that has particle (P) and hole (H) poles at energies k 0 = E X P,H (k) with residues Z X P,H (k), where we use X and Y to denote some SM fermion species. As the ψ Dirac mass is significantly larger than the thermal mass, we can safely utilise the free, non-interacting, thermal propagator. Following Ref. [54], we introduce the thermal widths listed in Table 10 by taking ω k − i → E(k) = ω k − iΓ. The rates and source terms are then given by, (B.6e) The formulae E X P,H (k) and Z X P,H (k) are derived in Ref. [85]. Due to the large vector-like mass the primary contributions to the integrals in Λ 0,± Xψ come from regions where the momenta are much larger than the thermal masses of the SM fermion. In this large momentum limit we have Z X P ≈ 1, Z X H ≈ 0 so that these integrals are in agreement with those used in Refs. [17,20,80].
The Yukawa-loop self energy function I F is taken from Ref. [76], In the event where I F = 0 due to the mass-thresholds we utilise the approximation for the four-body rates introduced in Refs. [15,76]. For our benchmarks this is only relevant for Y QT H . Due to the insignificance of quark and Higgs densities this approximation is expected to introduce negligible error. In units of GeV, the rates evaluated for benchmark point A are given in Table 11. In this benchmark the CPV source across the bubble wall is, dzS CP V Lψ (z) = −1.98 · 10 −2 GeV 3 (B.8) Rate Γ(z < 0) (GeV) Γ(z > 0) (GeV) Γ +