Electroweak Baryogenesis with Vector-like Leptons and Scalar Singlets
Nicole F. Bell Matthew J. Dolan Leon S. Friedrich Michael J. Ramsey-Musolf Raymond R. Volkas


ACFI-T19-02




a]Nicole F. Bell,
a]Matthew J. Dolan,
a,1]Leon S. Friedrich,Corresponding author
b,c,d]Michael J. Ramsey-Musolf,
a]and Raymond R. Volkas,
[a]ARC Centre of Excellence for Particle Physics at the Terascale,
  School of Physics, The University of Melbourne, Victoria 3010, Australia
[b]Amherst Center for Fundamental Interactions, Department of
  Physics, University of Massachusetts Amherst, Amherst, MA 01003, USA and California Institute of Technology, Pasadena, CA 91135 USA
 [c]Tsung-Dao Lee Institute and  School of Physics and Astronomy, Shanghai Jiao Tong University, 800 Dongchuan Road, Shanghai, 200240 China
[d]Kellogg Radiation Laboratory, California Institute of Technology, Pasadena, CA 91125 USA
n.bell@unimelb.edu.au matthew.dolan@unimelb.edu.au
leon.friedrich@unimelb.edu.au
mjrm@physics.umass.edu
raymondv@unimelb.edu.au





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.

Electroweak Baryogenesis with Vector-like Leptons and Scalar Singlets
    [
    December 30, 2023
=====================================================================



§ 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 <cit.>, indicates a baryon-to-entropy ratio

    Y_B = n_B/s = (8.66 ± 0.04) · 10^-11.

Any mechanism that might explain the origin of this asymmetry must satisfy the
three Sakharov conditions <cit.>: 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 first-order (SFO) electroweak phase
transition, or SFOEWPT. However, the Standard Model (SM) instead features a
crossover transition <cit.> 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 <cit.>

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, 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 <cit.>,
or via high sensitivity fundamental symmetry tests such as electric dipole
moment (EDM) searches, which place stringent constraints on additional sources
of CPV <cit.>.

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 <cit.>),
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[By
  partial seclusion we refer to scenarios where the CPV asymmetries are
  generated in species that carry no SM gauge quantum numbers (such that there
  are no EDM effects). These asymmetries get transferred to the SM via particle
  number changing reactions.]. 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 <cit.>. 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 <cit.>.

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 Yukawa-interactions with the gauge singlet
scalars and SM leptons, while also avoiding any contribution 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 <cit.> 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. <cit.>. 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 <cit.>.
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 <cit.>. 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 𝒪(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) 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. <cit.>. Using
8 ATLAS data, Ref. <cit.> places a lower bound of ∼ 270
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 <ref> we define our model along with a discussion of the
interactions, mass mixing, and new CPV phases. Section <ref>
introduces some benchmark parameters and discusses various observables and
constraints. In Section <ref> 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 <ref>. Technical
details relating to the electron EDM and transport equations are given in
Appendices <ref> and <ref>, 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

    S_i ∼ (1,1,0),   ψ∼(1,2,-1/2)  .

We use the notation

    ψ = [   N; E^- ],   
      H = [                      G^+; 1/√(2) (v_H + h + i G^0) ],  
      S_i =  v_S_i + s_i ,

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,

    - ℒ^Y_ψ =
         ∑_i [ λ_L ψ i L̅_3ψ_R S_i   + λ_ψψ i ψ̅_Lψ_R S_i  ]   + Y_ψ ψ̅_L  (i σ_2 H^*)  τ_R   + h.c.,

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

    - ℒ^M_ψ = M_ψψ̅_L ψ_R + M_L ψL̅_3 ψ_R + h.c.


The M_L ψ mass mixing term can be removed via a redefinition of the
ψ_L and L_3 fields,

    [ ψ_L; L_3 ]   →[ ψ_L'; L_3' ]
        =
        𝒰·[ ψ_L; L_3 ],     𝒰 = 
        [   e^ i ϕ_1 cosθ    e^ i ϕ_2sinθ; - e^ i ϕ_1 sinθ   e^ i ϕ_2 cosθ ] 
        , 
    sinθ   = | M_L ψ|/√(|M_ψ|^2 + |M_L ψ|^2  ),     ϕ_1= Arg(M_ψ),     ϕ_2= Arg(M_L ψ),

  such that eq. <ref> becomes,

    - ℒ^M_ψ = √(|M_ψ|^2 + |M_L ψ|^2  ) ψ̅'_L ψ_R' + h.c.  =   M_ψ' ψ̅'_L ψ_R' + h.c.

We can then perform a redefinition of the Yukawa couplings,

    [ Y_ψ Y_τ ]→[ Y_ψ' Y_τ' ]
      =
      [ Y_ψ Y_τ ]·𝒰^†  ,

and similarly for λ_L ψ i and λ_ψψ i. Aside from
removing M_L ψ this transformation leaves the form of the
Lagrangian unchanged and we drop the primes.




 §.§ 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

    V_0  =    -    μ_H^2 H^† H  + λ_H (H^† H)^2 
        
    
        +    ∑_i a_i S_i  - 1/2∑_i,j a_ij^2 S_i S_j   +  1/3∑_i,j,k a_i j k S_i S_j S_k +  1/4∑_i,j,k,l a_i j k l S_i S_j S_k S_l  
        
    
        +    1/2∑_i,j  b_ij S_i S_j H^† H  + ∑_i b_i S_i H^† H,

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

    . ∂ V/∂ v_S_i. = a_i + 1/2 v_H^2 b_i = 0,

which we use to fix a_i. We then require that the global minimum resembles
that of the SM, resulting in a SM-like Higgs with mass and VEV,
  
    (v_H, v_S_1, v_S_2) = (246, 0 , 0 ),      m_h = 125.

This fixes λ_H = μ_H^2/2 v_H^2 and μ_H^2 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

    a_12 = - 1/2 b_12 v_H^2.


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 <ref>). 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

    ψ_R →ψ_R e^i α,        ψ_L →ψ_L e^i β,
        
        L_3→ L_3 e^i γ,          τ_R→τ_R e^i δ.
 With these rephasings we find that the Lagrangian is unchanged if we take

    λ_L ψ i    →λ_L ψ i  e^i (γ-α),      Y_τ   → Y_τ e^i (γ - δ), 
        
    λ_ψψ i   →λ_ψψ i e^i (β-α) ,
            Y_ψ   → Y_ψ e^i (β -δ),
    
        M_ψ   → M_ψ e^i (β-α),
 such that we have the following invariant phases,

    δ_1    = Arg(λ_L ψ 1λ_L ψ 2^*),   δ_2    = Arg(λ_ψψ 1 M^*_ψ), 
    δ_3     = Arg(λ_ψψ 2 M^*_ψ),   δ_4    = Arg(Y_τ^* Y_ψλ_L ψ 2 M^*_ψ).


We will find that δ_1–δ_3 will appear within the EWBG
calculation as the sources of the CP-violation necessary to generate the
asymmetry. δ_4 can contribute to electron EDM constraints but will not
directly contribute to the asymmetry generation.

From here on we choose to re-phase our
couplings such that the following conditions are satisfied

    Arg( M_ψ)  = Arg( Y_τ)   =  Arg( λ_L ψ 2)   =   0

which is equivalent to setting

    δ_1    = Arg(λ_L ψ 1),  
           δ_2    =Arg(λ_ψψ 1),  
    δ_3    =Arg(λ_ψψ 2) , 
           δ_4    =Arg(Y_ψ) .

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,
 
  
  
    1/2ϕ^T ℳ^2_ϕϕ + (  ℰ̅_L ℳ_ℰℰ_R + 𝒩̅_L ℳ_𝒩𝒩_R  + h.c.),

  where,
  
    ϕ   =
          [   h; s_1; s_2 ]
        ,    ℰ   =
                          [ E; τ ]
        ,    𝒩_L    =
                            [ N_L; ν_τ ]
        ,    𝒩_R    =
                            [ N_R;   0 ]
        ,

  
    ℳ^2_ϕ
         =
          [       - μ_H^2 + 3 λ v_H^2                   v_H b_1                   v_H b_2;                   v_H b_1 - a_11^2 + 1/2 b_11 v_H^2                         0;                   v_H b_2                         0 - a_22^2 + 1/2 b_22 v_H^2 ],

    
    ℳ_ℰ   =
                         [          M_ψ v_H Y_ψ/√(2);            0 v_H Y_τ/√(2) ],
                          ℳ_𝒩   =
                                       [ M_ψ   0;   0   0 ].


These mass matrices can be diagonalised by a redefinition of the fields

    ϕ' = 𝒫ϕ
      ,     ℰ'_L/R = 𝒰^†_L/Rℰ_L/R
      ,

such that

  
  
    𝒫ℳ^2_ϕ𝒫^T    =  Diag(m^2_h,m^2_s_1,m^2_s_2),
        
    𝒰^†_Lℳ_ℰ𝒰_R   =  Diag(m_E, m_τ) .


The ℳ_𝒩 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 𝒫 are well-approximated by,


    𝒫_i i   ≈ 1, 
    𝒫_1 i   ≈ - 𝒫_i 1≈(M_ϕ^2)_1 i/(M_ϕ^2)_11 - (M_ϕ^2)_ii , 
    𝒫_2 3   ≈(M_ϕ^2)_1 3 (M_ϕ^2)_1 2/( (M_ϕ^2)_22 - (M_ϕ^2)_33)((M_ϕ^2)_11 - (M_ϕ^2)_22) , 
    𝒫_3 2   ≈𝒫_2 3(M_ϕ^2)_11 - (M_ϕ^2)_22/ (M_ϕ^2)_11 - (M_ϕ^2)_33 .




For the fermions, in the case where M_ψ in ℳ_ℰ is
much larger than any other term, we can approximate the mixing matrices via

    𝒰_L,R
      =
      [  e^i δ_4cosθ_L,R -e^i δ_4sinθ_L,R;         sinθ_L,R         cosθ_L,R ],

where

    θ_R ≈v_H | Y_ψ|/√(2) M_ψ,     θ_L ≈ v_H^2 Y_τ| Y_ψ|/2 M_ψ^2,

and δ_4 is given in eq. (<ref>). The τ mass defined
in eq. (<ref>) will be shifted slightly from its SM
relation with Higgs Y_τ Yukawa and is given by,

    m_τ≈Y_τ v_H/√(2)(  1 - v_H^2 | Y_ψ|^2  /4 M_ψ^2).

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 the correction to Y_τ will be orders of
magnitude smaller than the current uncertainty in Y_τ inferred from Higgs
decay measurements <cit.>. The mass of the charged component of
ψ will also be shifted slightly,

    m_E ≈ M_ψ(  1 + 
         v_H^2 | Y_ψ|^2 /4 M^2_ψ).

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 <ref>, where we have employed the notation of
Ref. <cit.>, 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 <ref>
and <ref>. The resulting mixing parameters are provided in
Table <ref>. 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 <ref>. 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 <cit.> (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_s_i^2 ≈ - a_ii +
1/2 b_ii v_H^2 ≲ - a_ii + 2 π v_H^2. We thus expect that
if we have the desired phase transition the scalars will generally be lighter
than the vector-like leptons.



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 <ref>, 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 their decays are mainly
into bb̅, W W^(*) and Z Z^(*). 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 WW^(*)/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  2.3.29 <cit.>. Using
the UFO format <cit.> we import our model into
 2.6.1 <cit.> in order to
generate events at parton level, before showering them using
 8.230 <cit.>. These are then fed into
 2.0.26 <cit.>, which
includes an interface with  3.4.1 <cit.>, 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 leads to K-factors ranging from 1.17–1.25 for lepton masses ranging from 100–1000 <cit.>. Based on this result we scale the ψ production cross sections passed to  using a K-factor of 1.2. 

For the benchmark points listed in Table <ref> we varied
M_ψ from 200 to 980 in steps of 20 and generated two million
events of pp→ℰℰ,ℰ𝒩,𝒩𝒩 for each parameter
point.  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,

    r = S - 1.64 Δ S/S_95,

provides a conservative quantification of the constraints on a given model, with
any point having r>1 considered excluded <cit.>. The resulting
r-values are shown in Fig. <ref>. 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. <cit.> 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. The lowest mass bound of the considered
benchmark 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 <ref>. 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.1 fb, 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 in-depth
analysis. See, e.g., Refs. <cit.> 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 that this leads to an 𝒪(1) correction to
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 <cit.>. We refer the
reader to
Refs. <cit.>
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 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.





 §.§ 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 <cit.> at,

    | d_e | < d_e^expt = 1.1 · 10^-29e·cm (90% CL).

We consider contributions to the
electron EDM via the two two-loop Barr-Zee style diagrams shown in Fig. <ref>. The relevant EDM formulae used to obtain numerical results
are provided in Appendix <ref>.




The diagram in Fig. <ref> could produce EDM contributions due to new
CPV interactions in the charged weak currents, as is the case in vector-like
quark models <cit.>. This diagram relies on a relative phase between the
left- and right-handed charged weak current interactions with d_e^W^- W^- ∝Im[M^L (M^R)^†], where M^L,R
are the left- or right-handed weak current mixing matrices that are given in
Appendix <ref>. However, in our scenario, these matrices are real,
such that there is no EDM contribution from this diagram. This 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. <ref>. 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. <cit.> 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


    ∑_i d_e^γϕ_i   ≈
                                 -( 56.2  |λ_ψψ 1 | 𝒫_1 2 sinδ_2
                               + 21.2  |λ_ψψ 2| 𝒫_1 3 sinδ_3 ) · 10^-29e·cm ,
                               
    ∑_i d_e^Z ϕ_i   ≈(7.90   |λ_ψψ 1 | 𝒫_1 2 sinδ_2
                               + 7.60   |λ_ψψ 2| 𝒫_1 3 sinδ_3) · 10^-29e·cm,


where 𝒫_ij are the components of the scalar mixing matrix defined
in Section <ref>, with numerical values given in
Table <ref>.

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. <ref> 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 δ_2 and
δ_3 phases.


  Full numerical electron EDM results for benchmarks A–C are given in
  Table <ref>. The predicted values for benchmark A and B are
  two orders of magnitude less than the current bound, whereas benchmark C is
  on the edge of being excluded at 90% confidence level. The EDM is much
  larger in benchmark C as it has significantly larger b_i couplings which
  lead to scalar mixing. No results are given for benchmark D as it does not
  feature scalar mixing, such that there is no two-loop electron EDM
  contribution and a higher order calculation is required. An increase in the
  sensitivity of EDM experiments by two orders of magnitude would place severe
  constraints on the scalar mixing matrix and the λ_ψψ i couplings. 

















 §.§ 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 <ref>. Due to the mixing between τ and E the Z
τ_R τ_R and W^±τ_L ν_τ interactions are modified,

    j^μ_Z⊃s_w^2/c_wτ̅γ^μℙ_R τ    → 1/c_w( s_w^2 + 1/2sin^2 θ_R ) τ̅γ^μℙ_R τ ,                                                                                         
    
      j^μ_W^-⊃1/√(2)τ̅γ^μℙ_L ν_τ    → 1/√(2)cosθ_L  τ̅γ^μℙ_L ν_τ  ,

where c_w and s_w are the sine and cosine of the weak mixing angle, and
ℙ_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

    Γ_τ∝∑_i=1^3 | (U_PMNS)_3 i|^2  ≈ 1 - θ_L^2  ,

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 <ref>, the corresponding universality-breaking corrections will be of
𝒪( 10^-7).

There will also be a contribution to the 1-loop Zττ coupling from vertex
corrections of the types shown in Fig. <ref>. 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 <cit.> and inserting
the benchmark couplings and mixing angles, we find that this vertex correction
is ∼ 0.01 | Y_ψ|^2/16 π^2∼𝒪(10^-9). This correction is many orders of magnitude
smaller than those found in vector-like quark models <cit.>, 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 <cit.>.









§ 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. <cit.>. 



 §.§ 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

    Γ_WS(T) = A(T) e^-  E_WS/T

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 scale v̅(T) associated with electroweak symmetry
breaking. For a discussion of the relationship between v̅(T) and the
scalar field VEV and the associated issue of gauge invariance, see
Ref. <cit.>. It is conventional to write

    E_WS =   4π B/g v̅(T)   ,

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 <cit.>.
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
first-order phase transition that satisfies this requirement is dubbed
“strong”, a characterisation that is usually translated into the requirement
that

    v̅(T_C)/T_C≳𝒪(1)    .

For recent discussions of sphaleron rate computations in perturbation theory and
the associated theoretical uncertainties see, e.g.,
Refs. <cit.>.


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. <cit.>.


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
out-of-equilibrium transport equations to compute the final baryon asymmetry.
For a more thorough review of EWBG
see <cit.>.



 §.§ Electroweak Phase Transition



In this section we outline our treatment of the electroweak and the selection of
the benchmark points listed in Section <ref>. 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,

    V^high-T = V_0 + δ m^2_H  H^† H T^2 + 1/2δ m^2_S_1 S_1^2 T^2  + 1/2δ m^2_S_2 S_2^2 T^2,

where V_0 is the tree level potential given in eq. (<ref>), and
the thermal masses δ m^2_X are computable in terms of the benchmark
model parameters and are listed in Table <ref> within
Appendix <ref>. This treatment has the advantage of being directly
gauge-independent, numerically fast to evaluate, and still yields results
comparable to a more thorough gauge-independent treatment, as seen in
Ref. <cit.>. Note that, in this context, the scalar field VEVs are
manifestly gauge invariant, and we will identify these VEVs with the
corresponding scale 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 
package <cit.> to determine T_C, the nucleation bubble wall
profiles, and the bubble nucleation probability density per unit
time <cit.>,

    Γ_nuc = T^4 ( S/2 π T)^3/2 e^- S/T  ,

where S is the three-dimensional Euclidean bubble action.
The expected number of bubbles in a Hubble volume is then given by <cit.>,

    ⟨ N(T)  ⟩ = ∫_T^T_c T/TΓ_nuc/H(T)^4.

The corresponding nucleation temperature is defined by the criterion

    ⟨ N(T_N)  ⟩≳ 1
      
        Γ_nuc(T_N) ≳ H(T_N)^4 = ( 4 π^3 g_⋆ T^4/45 m^2_pl)^2 
        S/T_N - 3/2lnS/T_N≲  173.8 - 2 ln g_⋆ - 4 lnT_N/GeV,

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  to trace the minima of the
effective potential for benchmark point A are shown in
Fig. <ref>. 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 <cit.>. 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 <cit.>.



Determining the wall velocity and moving bubble wall profiles is a complicated task requiring
solving a set of integro-differential equations <cit.>.
However both nucleation and moving wall profiles are generally well approximated
by hyperbolic tangent  profiles with a characteristic wall width L,

    v_i(z) = Δ v_i/2( 1 ±tanh( z -  δ z_i/L_i) ),      i ∈{H,S_1,S_2},

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.


The finite width of the wall can be neglected in the limit where the width
is small relative to the distances that particles diffuse ahead of the wall. In
the thin wall limit the reaction rates and CPV source terms can be well approximated by
step and delta functions, respectively. The typical reaction-diffusion distances
for particle species are given by √(D/Γ), where D is a
diffusion constant arising due to scattering interactions with the plasma, and
Γ is a dominant relaxation rate that drives a species back towards
equilibrium. In our model the particles responsible for biasing sphalerons to
generate the asymmetry are leptons, which diffuse much further ahead of the
wall than heavy quarks do. This is in part due to larger
diffusion constants as a result of leptons not being charged under SU(3), and in part
due to smaller relaxation rates as a result of smaller Yukawa couplings. Hence,
we find L ≪√(D/Γ), and will neglect the finite width of the
wall.

However, as will be shown in Section <ref>, the magnitude of the CPV
source term is sensitive to the relative magnitudes of the VEVs within the
bubble wall, and thus an accurate expression for the moving wall profiles is
still necessary to evaluate the source term. For this purpose we utilise the
nucleation bubble wall profiles as obtained by . The
profiles for benchmarks A and C are shown in
Fig. <ref> along with fitted tanh-functions.
Parameters for the fitted profile functions, as defined in
eq. (<ref>), are given in Table <ref> along with
the nucleation temperatures. The differences in the phase transitions between
benchmarks A and B are negligible, and similarly for C and D. Studies of moving wall dynamics typically find wall velocities w
∼ 0.1–0.3 <cit.>,
so for our benchmarks we will take w = 0.1.







Note that in general the nucleation VEV profiles differ from the moving wall
profiles. In particular, the wall widths of moving walls tend to be smaller than
the nucleation widths by a factor of 2–4 <cit.>.
However as mentioned earlier, we ignore the finite width of the wall and will
ignore this scaling. It should be noted that other studies of EWBG have found
that the asymmetry scales as 1/L. This is a result of the formalism
used to derive the CPV source terms. As will be discussed in the next section,
our CPV source term is proportional to the first derivatives of the VEVs. Studies using
the semi-classical force approximations instead find source terms proportional to the
second derivative <cit.>. This would
introduce an additional factor of 1/L and explains the discrepancy.






The presence of non-zero singlet VEVs in the unbroken phase will introduce a
mass term that leads to mixing between the ψ and SM leptons.
The singlet VEVs modify the Dirac mass terms,

    M_ψ→ M_ψ + λ_ψψ 1 v_S_1 +
      λ_ψψ 2 v_S_2, 
        
      M_L ψ→ M_L ψ + λ_L ψ 1 v_S_1 +
      λ_L ψ 2 v_S_2 ,

where M_ψ L is the mass term that we removed via a field redefinition in
Section <ref>. This new mixing term can once again be removed in the
singlet phase via the same transformation, eq. (<ref>), with
M_ψ and M_ψ L appropriately modified. We perform such a rotation
and work in the mass basis within the singlet phase. As a result, the ψ and L_3
Yukawa couplings that will appear in the following sections will be primed to
denote the fact that they correspond to the transformed Yukawas, as given by
eq. (<ref>), not the zero temperature Yukawas. A benefit of
the singlet VEVs is that if the λ_ψψ i couplings have the right
phases, they can drive the ψ to be lighter in the early universe than they
are at zero temperature, which results in an enhanced asymmetry. The benchmark
λ_ψψ i and their phases were selected such that this is the
case.






 §.§ 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

    ∂_λ J^λ_i (x) = - ∑_ΓT^2/6Γ_i j …(x) (μ_i - μ_j ±… ) + S_i^CPV(x)   ,

where J^λ_i and μ_i are the number density current and chemical
potential of particle species i. The Γ_i j
  … are  equilibration rates that arise due to various interactions
in the plasma, while S_i^CPV is a CPV source term arising from
interactions with the bubble wall. The number density and chemical potential for species i are related as

    n_i    = N_i - N̅_i   = g_i ∫^3 k/(2 π)^3[ n_f(ω_i(k), μ_i) - n_f(ω_i(k), -μ_i) ]
            , 
       ≈k_i(m_i/T) T^2/6μ_i  +  𝒪(μ_i^3/T^3) 
            , 
    
      k_i( m_i/T)    = g_i 6/π^2∫_m_i/T^∞ x x e^x/(e^x + 1)^2√(x^2 - m_i^2/T^2),

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. (<ref>), we can re-write eq. (<ref>) as

    w ∂/∂ z n_i(z) - D_i  ∂^2/∂^2 zn_i(z) = - ∑_ΓΓ_i j …(z) (n_i/k_i - n_j/k_j±…) + S_i^CPV(z)   .

We use the diffusion constants D_i derived by Ref. <cit.>,
which are provided in Appendix <ref> 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-of-equilibrium finite-temperature field theory. We will
use rates and source terms calculated using the VEV-insertion approximation
(VIA) <cit.>, 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. <ref>, 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 is an estimate of the contribution to the source due to flavor mixing that arises at first order in a gradient expansion across the bubble wall. 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. <cit.>. For a discussion of terms arising at second order in gradients, we refer the reader to Ref. <cit.>. We also refer the reader to Refs. <cit.> and <cit.> for discussion of these issues in the context of flavor mixing involving scalar fields. 

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 <ref>.

When performing the VIA derivations, the type of diagram shown in
Fig. <ref> gives rise to two equilibration rates and a CPV source term.
These equilibration rates, which we denote as Γ^±_A B 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_A B^CPV, is non-zero only
within the bubble wall when the VEVs are changing. The absorptive parts of
Yukawa-loops, shown in Fig. <ref>, give rise to equilibration rates,
Γ_A B 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

    T^2/6Γ_WS∑_i ( 3 μ_Q_i + μ_L_i)      and    T^2/6Γ_SS∑_i ( 2 μ_Q_i - μ_u_i - μ_d_i),

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 left-handed quark
doublets, and L_i are the three generations of left-handed SM lepton doublets.
We approximate the sphaleron rates as <cit.>,

    Γ_WS  = 6 κ_WS  T α_w^5 ,     Γ_SS =  6 κ_SS  T α_s^4,

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

    μ_τ_L = μ_ν_τ = μ_L_3,    n_L_3 = n_τ_L + n_ν_τ ,    k_L_3 = k_τ_L + k_ν_τ,

with similar relations holding for other doublets. 


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


    μ_ℓ = μ_L_1= μ_L_2,        k_ℓ = k_L_1 + k_L_2,        n_ℓ = n_L_1 + n_L_2,
        
    μ_q = μ_u_i = μ_d_i = μ_Q_i,        k_q = ∑_i k_Q_i +  k_u_i + k_d_i,        n_q = ∑_i n_Q_i + n_u_i + n_d_i.


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. (<ref>) then reduces to

    Γ_WS( 9 n_q/k_q  + 2 n_ℓ/k_ℓ + n_L_3/k_L_3).

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 <cit.>. In the present case, however, we
will find that the weak sphaleron rate (∼ 5 · 10^-4) may be
comparable to the L_3 and ψ relaxation rates (∼ 10^-3–10^-2 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:



    w ∂/∂ zn_ψ(z) - D_ψ∂^2/∂^2 zn_ψ 
      =    
      - ( Γ_L_3ψ S_i ^- + Γ_L_3ψ  S_1  + Γ_L_3ψ  S_2)
        ( n_ψ/k_ψ  - n_L_3/k_L_3)
      
       - Γ_L_3 ψ S_i ^+( n_ψ/k_ψ + n_L_3/k_L_3)
           + S_L_3 ψ^CPV,
    
           
           
      w ∂/∂ zn_L_3 - D_L_3∂^2/∂^2 zn_L_3  =   ( Γ_L_3ψ S_i ^- + Γ_L_3ψ  S_1  + Γ_L_3ψ  S_2)
        ( n_ψ/k_ψ  - n_L_3/k_L_3)
      
       + Γ_L_3 ψ S_i ^+( n_ψ/k_ψ + n_L_3/k_L_3)
       - Γ_WS( 9 n_q/k_q + 2 n_ℓ/k_ℓ + n_L_3/k_L_3)
           - S^CPV_L_3 ψ ,
             
           
           
      w ∂/∂ zn_ℓ - D_ℓ∂^2/∂^2 zn_ℓ  =     
       - 2 Γ_WS( 9 n_q/k_q + 2 n_ℓ/k_ℓ + n_L_3/k_L_3)
         ,
    
         w ∂/∂ z n_q - D_q∂^2/∂^2 zn_q  =     
       - 9 Γ_WS( 9 n_q/k_q + 2 n_ℓ/k_ℓ + n_L_3/k_L_3)
    .


The relevant equilibration rates and CPV source terms are given by,



    S_L_3 ψ^CPV(z)    
     = - 2 Im[λ'_L ψ 1λ'^*_L ψ 2]  ( v_S_1(z)  / tv_S_2(z) -  v_S_2(z) / t v_S_1(z) ) 
                 Λ^0_L_3 ψ,
      
    Γ_L_3 ψ S_i ^±   = 12/T^2|λ'_L ψ 1  v_S_1(z) + λ'_L ψ 2 v_S_2 (z) |^2   Λ^±_L_3 ψ,
      
    Γ_L_3 ψ S_i    =  12/T^2|λ'_L ψ i|^2  I_F(m_ψ, m_L, m_S_i),


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 <ref>, and the formulae for
Λ^±,0 and I_F are provided in Appendix <ref> along with
numerical values for the reaction rates in benchmark A. We also consider a
more complex system of transport equations, outlined in
Appendix <ref>, which is suitable when the ψ Yukawa is large
and we can no longer use some of the earlier simplifying assumptions.


  Note that due to the mixing induced by the singlet VEVs the couplings leading
  to a nonzero CPV source term S^CPV∝Im[λ'_L ψ 1λ'^*_L ψ 2 ], are a non-trivial function of the λ_L ψ
    i, λ_ψψ i, and their associated phases
  δ_1–δ_3.

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. <ref> is to have two such
insertions, leading to CPV sources proportional to terms like
Im[Y'_ψY'^*_ψ] =
0. 








This system of four linearly independent ODEs are solved numerically using a
technique similar to the one outlined in appendix C of Ref. <cit.>. 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, the asymmetry is given by dividing the
baryon-density deep inside of the bubble by the entropy density s,

    Y_B =  n_q(+∞)/3 s  , 
          
      s  = 2 π^2 g_⋆/45 T^3   ,      g_⋆=108.75.






 §.§ EWBG Results





The resulting baryon asymmetries generated by each of the benchmark points are
listed in Table <ref>, for both the simplified and full
system of transport equations as given in Section <ref> and
Appendix <ref>, respectively. The number densities obtained by
solving the simplified system of transport equations for benchmark A are shown
in Fig. <ref>, 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, λ_ψψ i, M_ψ and Y_ψ), the bubble velocity
  w and the diffusion rate D_ψ are illustrated in
  Fig. <ref>. As the simplifying assumptions used to obtain the
  simplified transport equations no longer hold
  when Y_ψ becomes large, we have used the full system of transport
  equations for the contour plots.










  Figures <ref> and <ref> show the
  dependence of the asymmetry on the scalar Yukawa couplings. If |λ_L ψ i|≪|λ_ψψ i|, the role of
  the λ_L ψ i is to generate mixing between the ψ and L_3
  and produce large λ'_Łψ i which lead to a non zero CPV source
  term as given in eq. (<ref>). The λ_ψψ i play a
  dual role of both reducing the mass of the ψ in the singlet phase, and
  contributing to λ'_L ψ i and hence the CPV source term. Note
  that in Fig. <ref> there is a ring where the generated
  asymmetry crosses through zero. This occurs due to a cancellation of
  the imaginary parts of the Λ_L_3 ψ^0 function, which appears in
  the expression for the CPV source term eq. (<ref>) and is given in eq. (<ref>).
  
  
  
  
  
  
  

  Figure <ref> shows the dependence of the final asymmetry on the
  bubble wall velocity and ψ mass. The decrease in the asymmetry when w
  becomes large is due to the inability of the particles to efficiently diffuse
  ahead of the wall, while the decrease for low w is a result of the CPV
  source term being proportional to the velocity. Meanwhile decreasing the
  ψ mass increases the resulting asymmetry, with a resonance occurring when
  M_ψ∼ m_L_3.
  
  

  
  
  Figure <ref> 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 the L_3 equilibrates.
  
  
  
  
  
  







§ 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 <cit.>, this model shares
a trait with multi-step EWPT models <cit.> 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. 

  An increase in the lower bound on the ψ masses reduces the capability for
  asymmetry generation. This effect can be offset slightly by using large
  λ_ψψ i Yukawa couplings and selecting their phases to make
  the vector-like fermions lighter in the early universe. However, the Yukawa
  couplings considered in our benchmarks are already uncomfortably large, a
  factor of three larger than the top quark Yukawa. Additional, it is these
  couplings and their phases that provide the dominant EDM contribution. Hence
  increased collider constraints alongside more precise electron EDM
  measurements will be sufficient to rule out this model.



The model examined in this paper, which is complementary to the one in
Ref. <cit.>, 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.



We thank Phillip Basler, Christopher Lee, David Morrisey, Hiren
Patel and Gaham White for helpful discussions. Feynman diagrams were drawn using
the TikZ-Feynman package <cit.>. This work was supported by the
Australian Research Council. Michael J. Ramsey-Musolf was supported in part
under U.S. Department of Energy Contract DE-SC0011095.







§ ELECTRON EDM FORMULAE


  
  

 §.§ 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. <cit.>, such
that it is straightforward to apply their general electron EDM formulae if one
makes the following replacements,

    h^0, H^0    →ϕ_i,   χ^-_a    →ℰ_a,
     
        c_e^ϕ_i   →𝒫_i1,   χ^0_i    →𝒩_i.

We define a matrix M such that after moving to the mass basis the
charged weak currents become,

    j_W^-^μ    ⊃ 1/√(2)τ̅γ^μℙ_L ν_τ +  1/√(2)E̅γ^μ  N  
        
       
        =  
        -  ℰ̅'̅γ^μ( M^L ℙ_L + M^R ℙ_R )  𝒩',

where ℙ_L,R are the standard left- or right-projection operators.
Similarly for the neutral currents we define a matrix G such that,

    c_w j_Z^μ    ⊃ τ̅γ^μ( s_w^2 - 1/2ℙ_L )  τ +  (s_w^2 - 1/2) E̅γ^μ  E
        
       
        =  ℰ̅'̅γ^μ( G^R ℙ_R + G^L ℙ_L )   ℰ' .

For the Yukawa interactions we define matrices D_i such that,  

    h/√(2)( Y_ψE̅ + Y_ττ̅)
         ℙ_R τ +  ∑_i s_i
        (  λ_ψψ iE̅  + λ_L ψ iτ̅)
        ℙ_R E
          + h.c.
         = 
        ϕ_i e/√(2) s_wℰ̅'̅( D_i^R ℙ_R + D_i^L ℙ_L  ) ℰ',

Using the mixing angles and phases defined in Section <ref>, these
matrices are given by,

  
  
    M^L      = -1/√(2)[  cosθ_L  sinθ_L; -sinθ_L  cosθ_L;         ]
    ,     
    M^R      = -1/√(2)[  cosθ_R       0; -sinθ_R       0;         ]
    ,
    
    G^L      = (s_w^2 - 1/2)
    [ 1 0; 0 1;   ]
    ,     
    G^R      =    
    [ s_w^2 - 1/2 cos^2 θ_R        1/4sin (2 θ_R);        1/4sin (2 θ_R)   s_w^2 -1/2sin^2 θ_R;                       ],


    D^R_ϕ_1   ≈s_w/e[    √(2)|λ_ψψ 1| e^i δ_2𝒫_12+√(2)|λ_ψψ 2| e^i δ_3𝒫_13+| Y_ψ|θ_R                                                  | Y_ψ|+θ_LY_τ; √(2)|λ_L ψ 1| e^i (δ_1+δ_4)𝒫_12+√(2) e^i δ_4λ_L ψ 2𝒫_13+θ_RY_τ                                                  Y_τ-| Y_ψ|θ_L;                                                                ]
    ,
    
    D^R_ϕ_2   ≈s_w/e[ √(2)(e^i δ_2|λ_ψψ 1|+|λ_ψψ 2| e^i δ_3𝒫_23+|λ_L ψ 1| e^i (δ_1+δ_4)θ_L)                                   -| Y_ψ|𝒫_12-√(2)|λ_ψψ 1| e^i δ_2θ_R;   √(2)(e^i (δ_1+δ_4)|λ_L ψ 1|+e^i δ_4λ_L ψ 2𝒫_23-|λ_ψψ 1| e^i δ_2θ_L)                               -√(2)|λ_L ψ 1| e^i (δ_1+δ_4)θ_R-𝒫_12Y_τ;                                                                       ]
    ,
    
    D^R_ϕ_3   ≈s_w/e[         √(2)(e^i δ_3|λ_ψψ 2|-|λ_ψψ 1| e^i δ_2𝒫_23+e^i δ_4λ_L ψ 2θ_L)                                  -| Y_ψ|𝒫_13-√(2)|λ_ψψ 2| e^i δ_3θ_R; √(2)(e^i δ_4λ_L ψ 2-|λ_L ψ 1| e^i (δ_1+δ_4)𝒫_23-|λ_ψψ 2| e^i δ_3θ_L)                                      -√(2) e^i δ_4λ_L ψ 2θ_R-𝒫_13Y_τ;                                                                      ], 
    
        D_ϕ_i^L    = (D_ϕ_i^R)^† .


In the expression for D_i provided above we have dropped any terms second order in mixing angles
𝒫_i, j ≠ i and θ_L/R, though the numerical calculations
included these terms. 










 §.§ EDM Formulae


Using the notation introduced in Section <ref> and the results from Refs. <cit.>, the primary contributions to the electron EDM are given by

    d_e = ∑_i=1^3( d_e^γϕ_i + d_e^Z ϕ_i)  + d_e^W^- W^-,


    d_e^γϕ_i
      =  e α^2 𝒫_i1/8 √(2)π^2 s_w^2 m_e/m_W m^2_ϕ_i∑_a=1^2Im(D_ϕ_i, a a^R) M_ℰ_a∫_0^1  x j( 0 , r_ℰ_a h_1/x(1-x))/x (1-x),


    d_e^Z ϕ_i   = e α^2 𝒫_i1(s_w^2 - 1/2)/8 √(2)π^2 c_w^2 s_w^4 m_e/m_W m^2_ϕ_i∑_a,b=1^2Im(G^R_a b D^R_ϕ_i b a - G^L_a b D^L_ϕ_i b a) m_ℰ_b
           × ∫_0^1  x 1/xj( r_Z ϕ_i , x r_ℰ_a h_1 + (1-x) r_ℰ_a ϕ_i/x(1-x)),


    d_e^W^- W^-   =
        -∑_a,i=1^2 e α^2 m_e m_ℰ_a m_𝒩_i/8 π^2 c_w^4 m_W^4Im( M^L *_a i M^R_a i) 
        
           × ∫_0^1  z_1 ∫_0^1- z_1 z_2  Z log( K_a i/Z)/(K_a i - Z )^2  - 1/K_a i - Z.

Where we have used,

    j(x,y)    = x log(x)/x-1  -  y log(y)/y-1/x-y
        ,    
        Z    = (z_1 + z_2)(1-z_1- z_2),
        
    
        K_a i   = r_𝒩_i W + z(r_ℰ_a W - r_𝒩_i W)
        ,      r_x y    =  m^2_x/m^2_y.

Substituting the matrices defined in
eq. (<ref>) one finds that in our case
d_e^W^- W^-=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.









§ THERMAL PROPERTIES AND TRANSPORT EQUATION FUNCTIONS





 §.§ Thermal Properties



The thermal mass formulae and diffusion constants used are given in
Tables <ref> and <ref> respectively.


For massless fermions the thermal widths at zero momentum are given by <cit.>

    Γ≈∑_i g_i^2 T C_F,i/4 π,

where C_F denotes the quadratic Casimir invariant. In the limit of masses
heavy compared to the temperature one instead has,

    Γ≈∑_i g_i^2 T C_F,i/8 π,

which differs from the massless case by a factor of a half. The resulting widths for the relevant
particles are given in Table <ref>.






 §.§ 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. (<ref>) we will get  

    μ_L_1 = μ_L_2,   μ_Q_1= μ_Q_2,   μ_u_1 = μ_u_2 = μ_d_i .

The full system of quantum
transport equations to be solved is then,



    ∂_μ J_ψ^μ =    
       - Γ_ψτ_R H(   n_ψ/k_ψ -  n_τ_R/k_τ_R - n_H/k_H)
      - ( Γ_L_3 ψ S_i ^- + Γ_L_3ψ  S_1  + Γ_L_3ψ  S_2)
        ( n_ψ/k_ψ  - n_L_3/k_L_3)
      
       - Γ_L_3 ψ S_i ^+( n_ψ/k_ψ + n_L_3/k_L_3)
            - Γ_ψτ_R H^- (   n_ψ/k_ψ -  n_τ_R/k_τ_R)  - Γ_ψτ_R H^+ ( n_ψ/k_ψ + n_τ_R/k_τ_R)
           + S_L_3 ψ^CPV,
    ∂_μ J_L_3^μ =   Γ_L_3ψ  S_i ^+( n_ψ/k_ψ + n_L_3/k_L_3)
      + ( Γ_ψL_3 S_i ^ - + Γ_ψL_3 S_1  + Γ_ψL_3 S_2)
        ( n_ψ/k_ψ  - n_L_3/k_L_3)
      
              - Γ_L_3τ_R H(   n_L_3/k_L_3 -  n_τ_R/k_τ_R - n_H/k_H) - Γ_L_3 τ_R H^- (   n_L_3/k_L_3 -  n_τ_R/k_τ_R)
       
      - Γ_L_3 τ_R H^+ ( n_L/k_L + n_τ_R/k_τ_R)
       - Γ_WS∑_j ( 3 n_Q_j/k_Q_j + n_L_j/k_L_j)
           - S^CPV_L_3 ψ ,
    ∂_μ J_L_1,2^μ =    - Γ_WS∑_j ( 3 n_Q_j/k_Q_j + n_L_j/k_L_j), 
    ∂_μ J^μ_τ_R =   Γ_ψτ_R H(   n_ψ/k_ψ -  n_τ_R/k_τ_R - n_H/k_H)+ Γ_ψτ_R H ^- (   n_ψ/k_ψ -  n_τ_R/k_τ_R)
      + Γ_ψτ_R H^+ ( n_ψ/k_ψ + n_τ_R/k_τ_R)
      
       
          +  Γ_L_3τ_R H(   n_L_3/k_L_3        -  n_τ_R/k_τ_R - n_H/k_H) + Γ_L_3 τ_R H ^- (   n_L_3/k_L_3        -  n_τ_R/k_τ_R)
      + Γ_L_3 τ_R H ^+ ( n_L/k_L + n_τ_R/k_τ_R)
                                   ,
    ∂_μ J^μ_Q_1,2       =    
     - 2 Γ_SS∑_j ( 2 n_Q_j/k_Q_j - n_u_j/k_u_j - n_d_j/k_d_j)
                    - 3 Γ_WS∑_j ( 3 n_Q_j/k_Q_j + n_L_j/k_L_j)
                    , 
    ∂_μ J^μ_Q_3       =      - 2 Γ_SS∑_j ( 2 n_Q_j/k_Q_j - n_u_j/k_u_j - n_d_j/k_d_j)
                    - 3 Γ_WS∑_j ( 3 n_Q_j/k_Q_j + n_L_j/k_L_j) 
       
        - Γ_Q_3u_3 H(  n_Q_3/k_Q_3 - n_u_3/k_u_3 + n_H/k_H) 
        - Γ_Q_3 u_3 H ^- (  n_Q_3/k_Q_3 - n_u_3/k_u_3) 
      - Γ_Q_3 u_3 H ^+ ( n_Q_3/k_Q_3 + n_u_3/k_u_3)
                    , 
    ∂_μ J^μ_d_1,2,3       =   ∂_μ J^μ_u_1,2  =
     Γ_SS∑_j ( 2 n_Q_j/k_Q_j - n_u_j/k_u_j - n_d_j/k_d_j)  , 
    ∂_μ J^μ_u_3       =   ∂_μ J^μ_d_1  
      - ∂_μ J^μ_Q_3       +  ∂_μ J^μ_Q_1, 
    ∂_μ J^μ_H       =      
        - Γ_Q_3u_3 H(  n_Q_3/k_Q_3 - n_u_3/k_u_3 + n_H/k_H) 
           + Γ_ψτ_R H(   n_ψ/k_ψ -  n_τ_R/k_τ_R - n_H/k_H) 
      
       
           + Γ_L_3τ_R H(   n_L_3/k_L_3 -  n_τ_R/k_τ_R - n_H/k_H) 
    .


Here n_H = n_H^++n_H^0 are the Higgs densities, and the scalar singlets
do not appear as they have zero chemical potential. Some of the relevant equilibration rates and CPV source terms are given in eqs. (<ref>) while the remainder are given here,



    Γ_ψτ_R H^±(z)    = 12/T^2| Y'_ψ v_H(z) |^2 Λ^±_ψτ_R,
      
    Γ_L_3 τ_R H^±(z)    = 12/T^2| Y'_τ v_H(z) |^2 Λ^±_L_3 τ_R,
      
    Γ_Q_3 u_3 H^±(z)    = 12/T^2| Y_u_3 v_H(z) |^2 Λ^±_Q_3 u_3,
      
    Γ_L_3 τ_R H    = 12/T^2| Y'_τ|^2 I_F(m_L, m_τ_R, m_H),
      
    Γ_ψτ_R H    = 12/T^2| Y'_ψ|^2 I_F(m_ψ, m_τ_R, m_H),
      
    Γ_Q_3 u_3 H    = 12/T^2| Y_u_3|^2 I_F(m_Q_3, m_u_3, m_H),







 §.§ Transport Equation Formulae



The relaxation rates and source terms used in eqs. (<ref>) and
(<ref>) were derived using the VEV insertion approximation as outlined
in Ref. <cit.>. 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 = ℰ^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. <cit.>, we introduce the thermal widths
listed in Table <ref> by taking ω_k - i ϵ→ℰ(k) = ω_k - i Γ. The rates and source terms
are then given by,


    Λ^0_X ψ  = Im∫_0^∞ k^2   k /π ^2 ω_ψ [    
      Z^X_P ℰ_ψ+ k/( ℰ_ψ + ℰ^X_P )^2( n_f(ℰ_ψ) + n_f(ℰ^X_P) ) 
      
     +    
      Z^X *_P ℰ_ψ- k/(ℰ_ψ - ℰ^X *_P)^2( n_f(ℰ_ψ) - n_f(ℰ^X *_P) ) 
      
     +    
      Z^X_H ℰ_ψ- k/(ℰ_ψ + ℰ^X_H)^2( n_f(ℰ_ψ) + n_f(ℰ^X_H)  ) 
      
     +    
      Z^X *_H ℰ_ψ + k/(ℰ_ψ - ℰ^X *_H)^2( n_f(ℰ_ψ) - n_f(ℰ^X *_H) )  ],


    Λ^±_X ψ =  -Im∫_0^∞ k^2  k /T π^2 ω_ψ[
          Z^X_Pℰ_ψ +k/ℰ_ψ +ℰ^X_P(h_f(ℰ_ψ ) ∓  h_f(ℰ^X_P) )
      
     +    
      Z^X *_Pℰ_ψ-k/ℰ_ψ-ℰ^X *_P( h_f(ℰ_ψ) ∓ h_f(ℰ^X *_P) ) 
      
     +     Z^X_Hℰ_ψ -k/ℰ_ψ +ℰ^X_H( h_f(ℰ_ψ ) ∓ h_f(ℰ^X_H) ) 
      
     +    
      Z^X *_Hℰ_ψ+k/ℰ_ψ-ℰ^X *_H( h_f(ℰ_ψ) ∓ h_f(ℰ^X *_H) ) ],


    Λ^±_Y X =  Im∫_0^∞ k^2  k /T π^2 [   Z^X_P Z^Y_P/ℰ^X_P + ℰ^Y_P( h_f(  ℰ^X_P) ∓ h_f(  ℰ^Y_P) )
      
     +    Z^X_P Z^Y *_H/ℰ^X_P - ℰ^Y *_H( h_f(  ℰ^X_P) ∓ h_f(  ℰ^Y *_H) )    +  ( P ↔ H ) ],


    ℰ_ψ(k)    = ω_ψ(k) - i Γ_ψ = √(k^2 + M'^2_ψ) - i Γ_ψ
                            ,                                                                                                          
      n_f(x)    = 1/1 + e^x/T
               , 
    
                   h_f(x)    = - T n_f'(x) = e^x/T/(1 + e^x/T)^2  
                            .


The formulae ℰ^X_P,H(k) and Z^X_P,H(k) are derived in Ref. <cit.>. 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. <cit.>.


The Yukawa-loop self energy function I_F is taken from Ref. <cit.>,


    ℐ_F (m_1, m_2, m_ϕ) 
                   =
                   1/16 π^3 T (m_1^2 + m_2^2 - m_ϕ^2) ∫_m_1^∞ω_1 ∫_ω_ϕ^-^ω_ϕ^+ω_ϕ
           ×{
                   n_B (ω_ϕ) [1 - n_F(ω_1)]n_F(ω_1 - ω_ϕ) [θ(m_1 - m_2 - m_ϕ) - θ(m_ϕ - m_1 - m_2)]
                   
              + n_B(ω_ϕ) n_F(ω_1) [1 - n_F(ω_1 + ω_ϕ)] θ(m_2 - m_1 - m_ϕ)
                   },


    ω^±_ϕ =   1/2 m_1^2( ω_1 | m_ϕ^2 + m_1^2 - m_2^2 |
         ± √(                                         (ω_1^2 - m_1^2) (m_1^2 - (m_2 + m_ϕ)^2)(m_1^2 - (m_2 - m_ϕ)^2)	)) .


In the event where I_F = 0 due to the mass-thresholds we utilise the
approximation for the four-body rates introduced in Refs. <cit.>. For our
benchmarks this is only relevant for Y_Q_3 u_3 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 <ref>. In this benchmark the CPV source across the bubble
wall is,

    ∫ d z S_L ψ^CPV(z) = -8.60 · 10^-2 GeV^3




