Primordial Magnetogenesis in the Two-Higgs-doublet Model

$\gamma$-ray emission of blazars infer the presence of large-scale magnetic fields in the intergalactic medium, but their origin remains a mystery. Using recent data from MAGIC, H.E.S.S. and $\textit{Fermi}$-LAT, we investigate whether the large-scale magnetic fields in the intergalactic medium could have been generated by a first-order electroweak phase transition in the two-Higgs-doublet model (2HDM). We study two representative scenarios where we vary the initial conditions of the magnetic field and the plasma, assuming either a primordial magnetic field with maximal magnetic helicity or a primordial magnetic field with negligible magnetic helicity in a plasma with kinetic helicity. By considering a primordial magnetic field with maximal helicity and applying the conservative constraints derived from MAGIC and $\textit{Fermi}$-LAT data, we demonstrate that a first-order electroweak phase transition within the 2HDM may account for the observed intergalactic magnetic fields in the case of the strongest transitions. We show that this parameter space also predicts strong gravitational wave signals in the reach of space-based detectors such as LISA, providing a striking multi-messenger signal of the 2HDM.


I. INTRODUCTION
Observations of blazars indirectly infer the presence of coherent magnetic fields in the intergalactic medium [1][2][3].The origin of these intergalactic magnetic fields (IGMF) is a longstanding puzzle, with two main scenarios proposed for their genesis: astrophysical and cosmological.Weak initial magnetic fields stem from localized effects within astrophysical objects (i.e., such as the Biermann battery mechanism [4]) that are subsequently amplified through dynamo effects [5] are an example of astrophysical phenomena that can create such long-correlated magnetic fields.Such mechanisms face challenges when it comes to explaining the magnetic fields in cosmic voids with a large volume filling factor [6].This makes potential sources originating from early universe processes especially compelling.Primordial magnetogenesis has been addressed in the context of inflation [7][8][9], post-inflationary reheating [10], the electroweak (EW) phase transition [11,12] or the QCD phase transition [13,14].Notably, the idea of magnetogenesis during a first-order electroweak phase transition (FOEWPT) was first introduced in Ref. [11].In this context, magnetic fields are generated through EW sphaleron decays [15,16].Since the bubbles nucleated during the EW phase transition expand, collide and merge, they stir the primordial plasma at high Reynolds number, and the magnetic fields enter a regime of magnetohydrodynamic (MHD) turbulence [17][18][19][20][21][22][23].Moreover, baryon number violating processes can produce magnetic fields with non-vanishing net helicity due to changes of the Chern-Simons number [15,16,24].Magnetic helicity has a significant impact on the evolution of magnetic turbulence over time, particularly by instigating an inverse cascade of magnetic energy, and transferring the energy from smaller to larger length scales [16,25], * mariaolalla.olearomacho@phys.ens.frand, thus, forming coherent magnetic structures at scales much larger than the ones where the energy was initially injected.The net magnetic helicity could be significant when there is substantial violation of the C and CP symmetries [23,26], as required by the Sakharov conditions to explain the matter-antimatter asymmetry of the universe [27].Therefore, the observed diffuse gamma ray sky may also hold information about the cosmological helical magnetic field and CP violation in the early universe [28,29].Furthermore, the violation of these symmetries might also affect the motion of the plasma, resulting in non-zero helicity of the velocity field [23].The presence of kinetic helicity in the plasma can ultimately source inverse cascades even in the context of non-helical magnetic fields [30].In this letter, we will study these two representative scenarios where we vary the initial conditions of the plasma and the magnetic field generated during a FOEWPT, assuming either a primordial magnetic field with maximal magnetic helicity or a primordial magnetic field with negligible magnetic net helicity in a plasma with initial kinetic helicity.
Extended Higgs sectors have the potential to explain the observed baryon asymmetry of the universe through EW baryogenesis [31].The simplest constructions to achieve EW baryogenesis extends the Standard Model (SM) scalar sector by a second Higgs doublet [32], referred to as the Two-Higgs-Doublet Model (2HDM).The FOEWPT in the 2HDM has been extensively studied in the literature [33][34][35][36][37][38][39][40], with many of these analyses taking a multi-faceted approach to the topic, combining collider signatures of this scenario with prospects for detecting a stochastic primordial gravitational wave (GW) background at future space-based GW observatories, such as the Laser Interferometer Space Antenna (LISA) [41,42].This study contributes to this multi-messenger investigation of the EW phase transition in the 2HDM by presenting predictions for the generation of primordial magnetic fields during the EW phase transition, and comparing them to the most recent constraints on the strength of the IGMF presented by the Fermi -LAT and H.E.S.S. col-laborations [2] and the Fermi -LAT and MAGIC collaborations [1].
This paper is structured as follows: in Sec.II A, we provide a brief overview of the CP-conserving 2HDM, together with the experimental and theoretical constraints that shape its physically allowed parameter space.We also present a 2HDM scenario that allows for a strong FOEWPT.In Sec.II B, we summarize the basis for the description of strong FOEWPTs and the computation of the GW spectrum.In Sec.II C, we outline how to compute the magnetic field spectrum today, depending on the two qualitatively different scenarios for the decay of MHD turbulence.In Sec.III, we discuss the results, and we conclude in Sec.IV.

A. The two-Higgs-doublet model
In this section, we provide a concise overview of the CP-conserving 2HDM which consists of two Higgs doublets, Φ 1 and Φ 2 , each bearing a hypercharge of 1/2. 1e impose a softly broken Z 2 symmetry [44,45], Φ 1 → Φ 1 , Φ 2 → −Φ 2 , that prevents the existence of flavor changing neutral currents at tree-level when extended to the Yukawa sector.The tree-level scalar potential is given by where all parameters are real as a result of requiring CP-conservation.The fields vacuum expectation values (vevs) are After spontaneous symmetry breaking, the CP-conserving 2HDM gives rise to five physical mass eigenstates in the scalar sector: two CP-even neutral scalars h and H, one CP-odd neutral pseudoscalar A and a pair of charged states H ± .The rotation from the gauge basis to the mass basis involves two mixing matrices with the mixing angles α and β for the CP-even and the CP-odd/charged sector, respectively.
The state h is conventionally chosen as the lightest CP-even scalar and, for the remainder of the discussion, plays the role of the discovered Higgs boson [46,47] h 125 at m h = 125 GeV.Furthermore, we adopt the so-called alignment limit [48], where the lightest Higgs boson h with mass m h = 125 GeV has the SM-Higgs boson couplings at tree-level.As the two fields Φ 1 and Φ 2 transform differently under the Z 2 symmetry, they cannot be coupled both to the same SM fermions, which leads to four 2HDM configurations/types that avoid flavor changing neutral currents at tree-level, characterised by the Z 2 charge assignment of the fermion fields in the Yukawa sector.Therefore, to generically describe the 2HDM parameter space we must specify the Yukawa type, along with the following set of independent parameters: After considering both theoretical and experimental constraints and demanding a strongly FOEWPT, we will focus on a benchmark scenario with just three free parameters: m A , m H , and t β (see below).

Constraints
Employing the publicly available thdmTools code [39], we impose a set of constraints on the parameter space of the 2HDM.These constraints encompass both theoretical and experimental aspects, including: • Perturbative unitarity.

A benchmark scenario for a FOEWPT in the 2HDM
We explore a 2HDM parameter region that is especially compelling for the realization of a strong FOEWPT.In this scenario, the CP-odd scalar A and the charged scalars H ± are enforced to be mass-degenerate, m A = m H ± , a choice motivated by EW precision data.Moreover, we assume the alignment limit cos(β − α) = 0, which is favored by the measurements of the signal rates of h 125 [50,51].The Z 2 -breaking mass scale M 2 = m 2 12 /(s β c β ) is set equal to the mass of the heavy CPeven scalar H, i.e.M = m H .These conditions allow for large m A − m H mass splittings, driven by the requirement for large quartic couplings in the 2HDM potential, which enable a FOEWPT [33,34,38] while agreeing with the LHC data on the 125 GeV Higgs boson, the measurements of EW precision observables, and other theoretical constraints.
Furthermore, here we will concentrate on the Yukawa type II.Given that all four Yukawa types share identical couplings between the neutral scalars and the top quark, we expect only minor differences among these types concerning the phase transition dynamics.
Taking into account the above, the remaining free parameters are m H , m A = m H ± , and t β .

B. Thermal history analysis
To investigate the dynamics of the EW phase transition in the 2HDM, we will follow the formalism of the temperature-dependent effective potential (see e.g.Ref. [52] for a review).The full effective potential is given by The temperature-independent part consists of the treelevel scalar potential V tree , given by Eq. ( 1), and the oneloop Coleman-Weinberg potential V CW [53].We also include a set of UV-finite counterterms V CT which enforces the zero-temperature loop-corrected vacuum expectation values, scalar masses and mixing angles to be equal to their tree-level values [36].The temperature-dependent part comprises the one-loop thermal corrections V T to the scalar potential [54], and the resummation of the daisydiagrams V daisy in the Arnold-Espinosa scheme [55].
It is important to note that this approach carries significant theoretical uncertainties [56][57][58][59], and therefore the parameter space explored in the subsequent analysis should be considered only as indicative of a first-order phase transition.
In a first-order phase transition, the scalar fields involved in the transition undergo a discontinuous change from a high-temperature symmetric phase (the false vacuum) to a low-temperature broken phase (the true vacuum) as the universe cools down.In this scenario, bubbles of true vacuum nucleate and expand in the background of the false vacuum, eventually converting all space into the true vacuum.The false and true vacuums are characterized by the values of the scalar fields expectation value, which are zero and nonzero, respectively.The onset of the transition depends on the nucleation rate per unit time and volume [60][61][62][63] where S(T ) is the three-dimensional Euclidean action of the O(3) symmetric bounce solution.By requiring that on average one bubble is nucleated per horizon volume, we define the nucleation temperature [42] In the 2HDM, the temperature at which bubbles of the new phase start to form (nucleation) and the temperature at which they fill the whole space (percolation) are very similar [39].Therefore, we use the same symbol T * to refer to both temperatures, and we call it the transition temperature.We track the co-existing minima of the effective potential as a function of the temperature and compute the nucleation rate by using cosmoTransitions [64].
Apart from the transition temperature T * , the firstorder cosmological phase transition is characterized by three other quantities: (i) the strength at the transition temperature α (normalized to the energy density in the plasma at the time of nucleation), (ii) the inverse duration of the transition (in Hubble units) β/H, and (iii) the bubble wall velocity in the rest frame of the fluid (away from the bubble) v w .The strength of the transition is given by the ratio of the difference of the trace of the energy-momentum tensor between the two phases (false and true vacua) to the background radiation energy density, i.e. [41,42] where being the values of the potential in the false and true vacuum.The inverse duration of the transition in Hubble units β/H can be defined as Finally, the the fourth quantity that characterizes a cosmological first-order phase transition is the bubble wall velocity v w .Except for the case of ultrarelativistic bubbles [65,66], computing v w is a challenging task that requires solving a coupled system of Boltzmann and scalar field equations in a model-dependent approach (see e.g.Refs.[67][68][69][70][71][72][73][74][75][76].Recent results suggest that phase transition bubbles tend to expand with either v w ≈ c s (c s being the speed of sound of the plasma)2 or v w → 1 [72,77] (see also Ref. [70] for more discussion on bubble wall velocity estimates in BSM theories).In this work we fix v w = 0.6 as an optimistic case. 3

Electroweak baryogenesis
The baryon asymmetry of the universe can arise dynamically from EW baryogenesis [31].In this scenario, the FOEWPT provides the out-of-equilibrium condition, one of the requirements for successful baryogenesis according to the Sakharov conditions [27].To avoid the washout of the baryon asymmetry after the phase transition, the following criterion has been often used in the literature [79] where This condition also characterizes a strong FOEWPT.Besides this requirement, the SM needs additional sources of CP-violation to generate the matter-antimatter asymmetry.We will not perform a detailed study of EW baryogenesis here (see Ref. [40] for a recent study), but we will focus only on one of its essential elements, the FOEWPT, assuming that the properties of the transition are not significantly affected by the inclusion of CP-violation.

Gravitational wave relics
Cosmological phase transitions that are first-order may produce a stochastic GW background [17,18] that can be probed by future GW interferometers.For a FOEWPT, the GW spectrum peaks around mHz frequencies, which matches the optimal sensitivity range of the space-based interferometer LISA [80,81].In the 2HDM, the GW spectrum is mainly sourced by the plasma motions after the bubble collisions, in the form of sound waves and MHD turbulence, rather than by the bubble wall collisions themselves [35].We use numerical power-law fits to the results of hydrodynamical simulations of the thermal plasma, which modeled the GW stochastic background as a function of the four parameters defined above.The formulas for the GW spectral shapes, amplitudes and peak frequencies are taken from Ref. [38], which follows Refs.[41,82].The detectability of a stochastic GW signal at a GW observatory depends on the signal-to-noise ratio (SNR), which is given by where T is the duration of the experiment, h 2 Ω Sens is the nominal sensitivity of the detector, according to the mission requirements [83], and h 2 Ω GW = h 2 Ω sw +h 2 Ω turb is the spectral shape of the GW signal.We focus on the GW detectability with LISA, for which we assume an operation time of T = 7 years, and consider a GW signal to be detectable if SNR > 1.
C. Primordial magnetic fields

Magnetic field spectrum today
The decay of MHD turbulence is influenced significantly by initial conditions and physical processes which are presently uncertain.However, some predictions can be derived based on different assumptions [84].One of the factors that affects the decay of MHD turbulence is the magnetic helicity of the initial seed field [85], which is a measure of the twist and linkage of the magnetic field lines.The average helicity is defined as ⟨A • B⟩, where B = ∇ × A. In a highly conductive medium, the magnetic helicity is approximately conserved.This means that a maximally helical field has to increase its correlation length as it loses magnetic energy, which results in an inverse cascade of magnetic energy where the energy is transferred from smaller scales to larger scales, forming coherent magnetic structures at scales much larger than the ones where the energy was initially injected.Therefore, this phenomenon could have played a crucial role in the survival and evolution of primordial magnetic fields, which are correlated at very large length scales today.
For maximally helical magnetic fields, the MHD turbulence, B, decays as a power law in conformal time, t, with the following scaling relations for the magnetic energy and correlation length during the radiation-dominated epoch [86]: Numerical simulations have shown that non-helical magnetic fields, which have zero or negligible net helicity, can also undergo an inverse cascade of magnetic energy in the presence of a plasma with initial kinetic helicity [30]. 4n this case, the magnetic energy and correlation length, λ, scale as respectively, during the radiation-dominated epoch.These scaling laws are valid until the epoch of matterdomination, when the scale factor grows linearly with conformal time, i.e. a ∼ t.After recombination, the magnetic field redshifts like radiation, i.e.B ∼ a −2 .To express these two scenarios in a compact way, we define the parameters for the power-laws where the cases b = 0 and b = 1 correspond to the maximally helical and non-helical scenarios described above, respectively.
In the 2HDM scenario, where the expanding bubbles do not enter the run-away regime [35], the magnetic field generated by the bubble collisions can be safely ignored.Thus, the magnetic field energy density at percolation temperature can be estimated by [12,87] Here, ρ * = 3M 2 p H 2 * is the total energy density at the percolation temperature and K = κα/(1 + α) is the ratio of the kinetic energy density of the sound waves to the total energy density.κ is the fraction of the released vacuum energy that is converted into the kinetic energy density of the plasma [88] and it was obtained following Ref.[88].Additionally, ε turb denotes the fraction of sound wave kinetic energy density expended on magnetic field generation, with numerical simulations suggesting ε turb ≈ 0.1 [23,89,90].
The magnetic field spectrum today can be computed as [12] B 0 (λ) ≡ B(λ, t 0 ) = ( 16) for λ < λ 0 , which assumes a power-law spectrum for the magnetic field strength, with a spectral index of n = 2 at large scales.Here λ 0 denotes the field coherence scale redshifted to today [12] where the initial correlation length λ * is given by the bubble size at percolation The redshift factors are computed as where g * corresponds to the total number of relativistic degrees of freedom in entropy at the transition temperature T * .

Experimental constraints from blazar emissions
Active galactic nuclei jets that are approximately pointed in our direction are called blazars, and they are sources of TeV γ-rays.When these γ-rays collide with photons from the extragalactic background, they give rise to e + e − pairs.Subsequently, these pairs interact with photons from the cosmic microwave background (CMB), producing γ-rays with energies in the GeV range.This in turn changes the original spectrum of the blazars by reducing the number of TeV γ-rays and increasing the number of GeV γ-rays [91,92].In the presence of IGMF, the trajectories of the e + e − pairs can be deflected due to the Lorentz force and, if the field is strong enough, the final GeV photons are no longer directed towards the observer.Non-observation of these photons has been used to set lower bounds on the strength of the IGMF.
We considered two recent analyses of blazar observations that measure the minimum strength of the IGMF.
The first one, performed by the MAGIC γ-ray observatory, sets a lower bound of B > 1.8 • 10 −17 G for correlation lengths λ ≥ 0.2 Mpc [1].The second one, based on the data from the Fermi -LAT and H.E.S.S. collaborations, establishes a limit of B > 7.1 • 10 −16 G for coherence lengths of λ ≥ 1 Mpc, assuming a blazar duty cycle of t d = 10 yrs [2].The duty cycle duration is the main source of systematic uncertainty in the latter analysis, so in our conservative approach, we only use the bounds obtained for t d = 10 yrs.

III. RESULTS AND DISCUSSION
We considered two different cases for the evolution of the MHD turbulence during the cosmic expansion.We varied the initial conditions of the magnetic field and the plasma, assuming either a primordial magnetic field with maximal magnetic helicity (b = 0) or a primordial magnetic field with negligible magnetic helicity (b = 1) in a plasma with kinetic helicity.We expect some fraction of magnetic helicity to be generated after the phase transition from the decay of EW sphalerons [15].We can also anticipate that some helicity will be produced as the plasma vortices create twisted field lines, due to the interaction between magnetic field and plasma elements in an MHD system [12,93].It is plausible to consider that a realistic scenario may lie somewhere between the extreme case of maximal initial magnetic helicity (b = 0) and the case with negligible magnetic helicity and nonzero kinetic helicity (b = 1).Given a 2HDM benchmark, it is typically expected that the predictions for the value of B 0 and λ 0 would be larger for the case with b = 0.
In order to explore the connection between GWs at LISA and primordial magnetic fields, we show the relation between the peak value of the magnetic field B 0 and its coherence length λ 0 for our sample of points in Fig. 1.We varied the values of m H between 360 GeV ≤ m H ≤ 900 GeV and m A between 560 GeV ≤ m H ≤ 1050 GeV.We fixed the value of t β to 3, with the rest of the parameters set as described in Sec.II A 2. We imposed all the theoretical and experimental constraints listed in Sec.II A 1. We only considered the points that have ξ n ≳ 1, which are preferred by EW baryogenesis.The two sets of data points correspond to the two different cases for the evolution of the MHD turbulence (b = 0 and b = 1).The colorbar indicates the SNR in LISA for points that satisfy SNR > 10 −3 .For b = 0, the points with SNR < 10 −3 are indicated in blue, whereas for b = 1, we use gray.The solid coral line corresponds to the lower bound on the magnetic field set by Ref. [1], whereas the dashed line corresponds to that set by Ref. [2] for a duty cycle of t d = 10 yrs.Upper limits on the magnetic field strength derived from CMB and Big Bang Nucleosynthesis data exist for values of B 0 ≳ 1 nG [5].
Regardless of the choice for b, all benchmark points with SNR > 10 −3 lie above the most conservative bound set by MAGIC/Fermi -LAT.Moreover, data points with For b = 1 (non-helical magnetic fields in a plasma with kinetic helicity), the peak coherence length of the scan points ranges from 3.80 × 10 −9 Mpc ≤ λ 0 ≤ 3.12 × 10 −5 Mpc, and the magnetic field value lies within the interval 7.69×10 −7 nG ≤ B 0 ≤ 2.02×10 −4 nG.In this case, if we focus on the higher lower limit by Fermi -LAT and H.E.S.S., only the points with the strongest GW signals could account for the origin of the primordial magnetic fields.
For both scenarios (b = 0 and b = 1), this work exemplifies the need for continued improvement in the analysis of magnetic fields via blazar observations to exclude primordial seed fields sourced by the EW phase transition in the 2HDM.Conversely, GW interferometry can serve as an additional probe of a 2HDM origin of primordial magnetic fields over certain coherence length scales.
In Fig. 2, we show the (m H ,m A ) mass plane for the aforementioned settings for the rest of the parameters of the model.As mentioned above, only the points with ξ n ≳ 1 were considered, since they are favoured by EW baryogenesis.The largest value found for ξ n is 4.07.The SNR in LISA, assuming the experimental choices indicated under Eq.( 10) is depicted in the colorbar.
For b = 0, all parameter points in Fig. 2 produce primordial magnetic fields that could explain the observations made by both experiments.
Assuming b = 1, only the points inside the black dashed line exceed the less stringent lower bound on B 0 from MAGIC and Fermi -LAT observations.This region of parameter space also corresponds to the prediction of the strongest gravitational waves, which can be tested in future LHC runs in searches such as gg → A → ZH → ℓ + ℓ − t t [39].This is another experimental probe that could reveal whether primordial magnetic fields may have a 2HDM origin or not (for relatively low values of t β ).

IV. CONCLUSIONS
We explored the possibility that the large-scale magnetic fields in the intergalactic medium could have originated from a first-order electroweak phase transition in the two-Higgs-doublet model (2HDM).We considered two scenarios with different initial conditions for the magnetic field and the plasma, assuming either a primordial magnetic field with maximal magnetic helicity (defined by b = 0) or a primordial magnetic field with negligible magnetic helicity in a plasma with kinetic helicity (defined by b = 1).We calculated the maximum strength of the magnetic field spectrum B 0 for a representative parameter plane (see Fig. 2) that realizes a first-order electroweak phase transition in the 2HDM.We confronted these predictions with observational constraints on the intergalactic magnetic fields from MAGIC, H.E.S.S. and Fermi-LAT experiments.
We found that the strongest transitions, with a SNR for LISA larger than 10 −3 , can explain the γ-ray observations by Fermi-LAT/MAGIC, regardless of our choice of the initial conditions for the plasma and the magnetic field (see Fig. 1).
Fermi-LAT/MAGIC provide the most conservative bound on B 0 as a function of λ 0 , compared to the limit obtained by Fermi-LAT/H.E.S.S. Maximally helical magnetic fields, spanning scales of about 4.27 × 10 −7 Mpc ≤ λ 0 ≤ 2.94 × 10 −3 Mpc, are consistent with the bounds of both experiments.For non-helical magnetic fields, only a few points in our parameter space could exceed the stricter limit by Fermi-LAT/H.E.S.S.
However, all points that predict very strong gravitational wave signals in LISA (SNR ≥ 1) also predict a large peak magnetic field above the Fermi-LAT/H.E.S.S. lower bound, offering a multi-messenger test of the 2HDM origin of primordial magnetic fields.Additionally, the low-t β parameter space (see Fig. 2) can be explored in upcoming LHC runs.The process gg → A → ZH → ℓ + ℓ − t t [39] offers a collider-based approach to assess the 2HDM contribution to the genesis of intergalactic magnetic fields.
We concluded that a first-order electroweak phase transition within the 2HDM may explain (over certain coherence length scales) the intergalactic magnetic fields inferred from blazar observations, depending on the initial conditions of the magnetic field and the plasma.We also demonstrated that this scenario could be tested by future LHC runs and LISA observations, providing a tantalizing multi-messenger signal of the 2HDM.We also emphasized the need for improved analysis of magnetic fields via blazar observations to exclude or confirm such model predictions, i.e. primordial magnetic fields sourced by the electroweak phase transition in the 2HDM.

2 SNRFIG. 1 .
FIG. 1. Magnetic field strength as a function of its coherence length.The solid coral line corresponds to the lower bound on the magnetic field set by Ref. [1], while the dashed line corresponds to the one set by Ref. [2] for a duty cycle of t d = 10 yrs.The points indicate the values of the peak magnetic field strengths B0 and the corresponding coherence lengths λ0, for all the points of the scan in Fig. 2 (for b = 0 and b = 1).The colorbar indicates the SNR in LISA.Points for which SNR < 10 −3 are shown in blue (b = 0) and in gray (b = 1).

2 SNRFIG. 2 .
FIG. 2. (mH ,mA) mass plane for t β = 3.The colorbar indicates the SNR at LISA.The dashed black line comprises the points that exceed the lower bound on B0 from MAGIC and Fermi-LAT observations for the case of non-helical magnetic fields (b = 1).