# Primordial Magnetogenesis in the Two-Higgs-doublet Model

María Olalla Olea-Romacho

► To cite this version:

María Olalla Olea-Romacho. Primordial Magnetogenesis in the Two-Higgs-doublet Model. Physical Review D, 2024, 109 (1), pp.015023. 10.1103/PhysRevD.109.015023. hal-04291515

HAL Id: hal-04291515

https://hal.science/hal-04291515v1

Submitted on 28 Aug 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Primordial magnetogenesis in the two-Higgs-doublet model

María Olalla Olea-Romacho*

Laboratoire de Physique de l'École Normale Supérieure, ENS, Université PSL,
CNRS, Sorbonne Université, Université Paris Cité, F-75005 Paris, France

<figure><img src="image_2.png" /></figure>

(Received 14 November 2023; accepted 3 January 2024; published 22 January 2024)

γ-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 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 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 multimessenger signal of the 2HDM.

DOI: 10.1103/PhysRevD.109.015023

# I. INTRODUCTION

Observations of blazars indirectly infer the presence of coherent magnetic fields in the intergalactic medium [1–3]. The origin of these intergalactic magnetic fields (IGMFs) is a long-standing 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] and 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–9], postinflationary 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–23]. Moreover, baryon number-violating processes can produce magnetic fields with nonvanishing 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], 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 γ-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 nonzero helicity of the velocity field [23]. The presence of kinetic helicity in the plasma can ultimately source inverse cascades even in the context of nonhelical magnetic fields [30]. In this paper, 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

Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI. Funded by SCOAP³.

<footnote>* mariaolalla.olearomacho@phys.ens.fr</footnote>
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–40], with many of these analyses taking a multifaceted 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 multimessenger 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. Collaborations [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.

## II. THEORETICAL BACKGROUND

### 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, $\Phi_1$ and $\Phi_2$, each bearing a hypercharge of $1/2$.<fnref n="1" /> We impose a softly broken $Z_2$ symmetry [44,45], $\Phi_1 \to \Phi_1, \Phi_2 \to -\Phi_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

$$V_{\text{tree}}(\Phi_1, \Phi_2) = m_{12}^2 \Phi_1^\dagger \Phi_1 + m_{22}^2 \Phi_2^\dagger \Phi_2 - m_{12}^2 (\Phi_1^\dagger \Phi_2 + \Phi_2^\dagger \Phi_1) + \frac{1}{2} \lambda_1 (\Phi_1^\dagger \Phi_1)^2 + \frac{1}{2} \lambda_2 (\Phi_2^\dagger \Phi_2)^2 + \lambda_3 (\Phi_1^\dagger \Phi_1) (\Phi_2^\dagger \Phi_2) \\ + \lambda_4 (\Phi_1^\dagger \Phi_2) (\Phi_2^\dagger \Phi_1) + \frac{1}{2} \lambda_5 ((\Phi_1^\dagger \Phi_2)^2 + (\Phi_2^\dagger \Phi_1)^2), \quad (1)$$

where all parameters are real as a result of requiring CP conservation. The field vacuum expectation values are $\langle\Phi_1\rangle = v_1$ and $\langle\Phi_2\rangle = v_2$, with $v_1^2 + v_2^2 = v^2 \approx 246$ GeV and $v_2/v_1 \equiv \tan\beta \equiv t_\beta$. 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^\pm$. The rotation from the gauge basis to the mass basis involves two mixing matrices with the mixing angles $\alpha$ and $\beta$ 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],

$$\cos(\beta - \alpha) = 0, \quad (2)$$

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 $\Phi_1$ and $\Phi_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, characterized 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:

$$t_\beta, m_{12}^2, v, \cos(\beta - \alpha), m_h, m_H, m_A, m_{H^\pm}. \quad (3)$$

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_\beta$ (see below).

### 1. 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

(i) perturbative unitarity,

(ii) vacuum stability, and

(iii) direct searches for beyond the Standard Model (BSM) scalars using HiggsTools [49].

<footnote n="1">¹See, e.g., Ref. [43] for a complete review of the two-Higgs-doublet model.</footnote>
## 2. 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(\beta - \alpha) = 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_{12}^2/(s_\beta c_\beta)$ is set equal to the mass of the heavy CP-even 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

$$V_{\text{eff}}(\phi_i, T) = V_{\text{tree}}(\phi_i) + V_{\text{CW}}(\phi_i) + V_{\text{CT}}(\phi_i) + V_{\text{T}}(\phi_i, T) + V_{\text{daisy}}(\phi_i, T). \quad (4)$$

The temperature-independent part consists of the tree-level scalar potential $V_{\text{tree}}$, given by Eq. (1), and the one-loop Coleman-Weinberg potential $V_{\text{CW}}$ [53]. We also include a set of UV-finite counterterms $V_{\text{CT}}$ that enforce 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 daisy diagrams $V_{\text{daisy}}$ in the Arnold-Espinosa scheme [55].

It is important to note that this approach carries significant theoretical uncertainties [56–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 field expectation values, which are zero and nonzero, respectively. The onset of the transition depends on the nucleation rate per unit time and volume [60–63],

$$\Gamma(T) = A(T)e^{-S(T)/T}, \quad (5)$$

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]

$$S(T_n)/T_n \sim 140. \quad (6)$$

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 coexisting 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 first-order cosmological phase transition is characterized by three other quantities: (i) the strength at the transition temperature $\alpha$ (normalized to the energy density in the plasma at the time of nucleation), (ii) the inverse duration of the transition (in Hubble units) $\beta/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],

$$\alpha = \frac{1}{\rho_R} \left( \Delta V(T_*) - \left. \left( \frac{T}{4} \frac{\partial \Delta V(T)}{\partial T} \right) \right|_{T_*} \right), \quad (7)$$

where $\Delta V(T) = V_f - V_t$, with $V_f \equiv V_{\text{eff}}(\phi_f)$ and $V_t \equiv V_{\text{eff}}(\phi_t)$ being the values of the potential in the false and true vacuum. The inverse duration of the transition in Hubble units $\beta/H$ can be defined as

$$\frac{\beta}{H} = T_* \left( \frac{d}{dT} \frac{S(T)}{T} \right) \Big|_{T_*} . \quad (8)$$

Finally, 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–76]). Recent results suggest that phase transition bubbles tend to
expand with either $v_w \approx c_s$ ($c_s$ being the speed of sound of the plasma)$^2$ or $v_w \to 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.<fnref n="3" />

## 1. 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]:

$$ \xi_n \equiv v(T_n)/T_n \gtrsim 1, \qquad (9) $$

where $v(T_n) = \sqrt{v_1(T_n)^2 + v_2(T_n)^2}$. This condition also characterizes a strong FOEWPT. In addition to 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.

## 2. 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 millihertz 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

$$ \text{SNR} = \sqrt{\mathcal{T} \int_{-\infty}^{+\infty} df \left[ \frac{h^2 \Omega_{\text{GW}}(f)}{h^2 \Omega_{\text{Sens}}(f)} \right]^2}, \qquad (10) $$

where $\mathcal{T}$ is the duration of the experiment, $h^2\Omega_{\text{Sens}}$ is the nominal sensitivity of the detector, according to the mission requirements [83], and $h^2\Omega_{\text{GW}} = h^2\Omega_{\text{sw}} + h^2\Omega_{\text{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 $\mathcal{T} = 7$ yr and consider a GW signal to be detectable if SNR > 1.

## C. Primordial magnetic fields

### 1. Magnetic field spectrum today

The decay of MHD turbulence is influenced significantly by initial conditions and physical processes that 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 $\langle \mathbf{A} \cdot \mathbf{B} \rangle$, where $\mathbf{B} = \nabla \times \mathbf{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]:

$$ B \sim t^{-1/3} \quad \text{and} \quad \lambda \sim t^{2/3}. \qquad (11) $$

Numerical simulations have shown that nonhelical 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].<fnref n="4" /> In this case, the magnetic energy and correlation length $\lambda$ scale as

$$ B \sim t^{-1/2} \quad \text{and} \quad \lambda \sim t^{1/2}, \qquad (12) $$

respectively, during the radiation-dominated epoch. These scaling laws are valid until the epoch of matter domination, when the scale factor grows linearly with conformal time, i.e., $a \sim t$. After recombination, the magnetic field redshifts

<footnote n="2">²For a relativistic perfect fluid, $c_s = 1/\sqrt{3} \approx 0.577$.</footnote>

<footnote n="3">³See Ref. [78] for estimates of the bubble wall velocity in a similar model.</footnote>

<footnote n="4">⁴This has been contested recently in Ref. [84], where a weaker inverse transfer of magnetic energy for nonhelical fields was found, compared to previous studies in the literature.</footnote>
like radiation, i.e., $B \sim a^{-2}$. To express these two scenarios in a compact way, we define the parameters

$$q_b = \frac{2}{b+3} \quad \text{and} \quad p_b = \frac{2}{b+3}(b+1), \qquad (13)$$

for the power laws

$$B \sim t^{-p_b/2} \quad \text{and} \quad \lambda \sim t^{q_b}, \qquad (14)$$

where the cases $b=0$ and $b=1$ correspond to the maximally helical and nonhelical scenarios described above, respectively.

In the 2HDM scenario, where the expanding bubbles do not enter the runaway 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]

$$\rho_{B,*} = \epsilon_{\text{turb}} K \rho_{*} = 0.1 \frac{\kappa \alpha}{1 + \alpha} \rho_{*}. \qquad (15)$$

Here, $\rho_*=3M_p^2H_*^2$ is the total energy density at the percolation temperature and $K=\kappa\alpha/(1+\alpha)$ is the ratio of the kinetic energy density of the sound waves to the total energy density. $\kappa$ 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, $\epsilon_{\text{turb}}$ denotes the fraction of sound wave kinetic energy density expended on magnetic field generation, with numerical simulations suggesting $\epsilon_{\text{turb}} \approx 0.1$ [23,89,90].

The magnetic field spectrum today can be computed as [12]

$$B_0(\lambda) \equiv B(\lambda, t_0) = \left(\frac{a_*}{a_{\text{rec}}}\right)^{p_b/2} \left(\frac{a_*}{a_0}\right)^2 \sqrt{\frac{17}{10}} \rho_{B,*} \begin{cases} \left(\frac{\lambda}{\lambda_0}\right)^{-5/2} & \text{for } \lambda \ge \lambda_0 \\ \left(\frac{\lambda}{\lambda_0}\right)^{1/3} & \text{for } \lambda < \lambda_0, \end{cases} \qquad (16)$$

which assumes a power-law spectrum for the magnetic field strength, with a spectral index of $n=2$ at large scales. Here $\lambda_0$ denotes the field coherence scale redshifted to today [12],

$$\lambda_0 \equiv \lambda_B(t_0) = \left(\frac{a_{\text{rec}}}{a_*}\right)^{q_b} \left(\frac{a_0}{a_*}\right)^{-1} \lambda_*, \qquad (17)$$

where the initial correlation length $\lambda_*$ is given by the bubble size at percolation $R_*$,

$$\lambda_* = R_* = \frac{(8\pi)^{1/3}}{H_*} \left(\frac{\beta}{H}\right)^{-1} v_w. \qquad (18)$$

The redshift factors are computed as

$$\frac{a_*}{a_0} = 8 \times 10^{-14} \left(\frac{100}{g_*}\right)^{1/3} \left(\frac{1 \text{ GeV}}{T_*}\right), \qquad (19)$$

$$\frac{a_*}{a_{\text{rec}}} = 8 \times 10^{-11} \left(\frac{100}{g_*}\right)^{1/3} \left(\frac{1 \text{ GeV}}{T_*}\right), \qquad (20)$$

where $g_*$ corresponds to the total number of relativistic degrees of freedom in entropy at the transition temperature $T_*$.

## 2. 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 $\gamma$ rays. When these $\gamma$ 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 $\gamma$ rays with energies in the GeV range. This, in turn, changes the original spectrum of the blazars by reducing the number of TeV $\gamma$ rays and increasing the number of GeV $\gamma$ 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 toward the observer. Nonobservation 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 $\gamma$-ray observatory, sets a lower bound of $B > 1.8 \times 10^{-17} \text{ G}$ for correlation lengths $\lambda \ge 0.2 \text{ 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 \times 10^{-16} \text{ G}$ for coherence lengths of $\lambda \ge 1 \text{ Mpc}$, assuming a blazar duty cycle of $t_d = 10 \text{ yr}$ [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 \text{ yr}$.

## III. RESULTS AND DISCUSSION

We considered two scenarios 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.

In Fig. 1, we show the peak value of the magnetic field spectrum $B_0 = B_0(\lambda_0)$ computed with Eq. (16) together with its corresponding coherence length $\lambda_0$ for the two different types of turbulence decay for a 2HDM benchmark point with $m_H = 383.13 \text{ GeV}$, $m_A = 632.67 \text{ GeV}$, and $t_\beta = 3$, which corresponds to the point with the maximum
<figure><img src="image_3.png" /></figure>

FIG. 1. Magnetic field strength as a function of its coherence length. The peak value of the magnetic field spectrum $B_0$ is plotted together with its corresponding coherence length $\lambda_0$ for two scenarios of turbulence decay: the scenario with $b=0$ (magenta star) and $b=1$ (blue star). The 2HDM benchmark point is characterized by $m_H = 383.13$ GeV, $m_A = 632.67$ GeV, and $t_\beta = 3$. The transition temperature is $T_* = 59$ GeV. 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$ yr.

<figure><img src="image_4.png" /></figure>

FIG. 2. ($m_H, m_A$) mass plane for $t_\beta = 3$. The color bar indicates the peak value of the magnetic field spectrum $B_0$ (assuming $b=0$). The dashed black line comprises the points that exceed the lower bound on $B_0$ from MAGIC and Fermi-LAT observations at the corresponding coherence length.

value for $\xi_n = 4.07$ for the mass plane shown in Fig. 2 (see discussion below). The rest of the parameters were set as described in Sec. II A 2. The first-order phase transition for this point is characterized by a nucleation temperature of

$T_n = 59$ GeV, a strength $\alpha = 0.12$, and an inverse time-scale of $\beta/H = 48$. We find that the scenario with no initial magnetic helicity produces a magnetic field spectrum that is much weaker than the lower limit on the IGMF derived from the blazar observations for a typical strongly FOEWPT in the 2HDM (such as the benchmark scenario shown in Fig. 1 with $\xi_n = 4.07$). 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]. Therefore, we will consider the optimistic hypothesis that the scenario with maximal initial magnetic helicity ($b=0$) is more realistic for the initial conditions of the magnetic field, and we will focus on this case in the subsequent analysis. The benchmark point in Fig. 1 shows that the peak value of the magnetic field $B_0$ for maximally helical fields is above the blazar bound from the MAGIC/Fermi-LAT and Fermi-LAT/H.E.S.S. Collaborations for a coherence length scale $\lambda_0$ around $2.9 \times 10^{-3}$ Mpc.

In Fig. 2, we show the ($m_H, m_A$) mass plane for a fixed value of $t_\beta = 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 $\xi_n \gtrsim 1$, which are preferred by EW baryogenesis. The largest value found for $\xi_n$ is 4.07, which also corresponds to the benchmark point analyzed in Fig. 1. Assuming $b=0$, we calculated the peak value of the magnetic field $B_0$, which is depicted in the color bar. The black dashed line indicates the points that exceed the lower bound on $B_0$ from MAGIC and Fermi-LAT observations at the corresponding coherence length $\lambda_0$.

The peak magnetic field strengths vary within the range of approximately $3 \times 10^{-17}$ to about $10^{-14}$ G across the entire dataset and from $5 \times 10^{-16}$ to roughly $10^{-14}$ G for the data points enclosed within the black dashed contour. This parameter region can be tested in future LHC runs in searches such as $gg \to A \to ZH \to \ell^+\ell^- t\bar{t}$ [39], which could reveal whether primordial magnetic fields could have a 2HDM origin or not (for relatively low values of $t_\beta$).

In order to explore the connection between GWs at LISA and primordial magnetic fields in the 2HDM, we show the relation between the peak value of the magnetic field and its coherence length for our sample of points in Fig. 3. The color bar indicates the SNR in LISA for points that satisfy SNR $> 10^{-3}$. Upper limits on the magnetic field strength derived from CMB and big bang nucleosynthesis data exist for values of $B_0 \gtrsim 10^{-9}$ G [5]. We find that the points with the highest SNR (SNR $> 1$), which might be detected by LISA, also have the highest peak values of the magnetic field $B_0$, and they all exceed the lower bound from blazar observations by MAGIC and Fermi-LAT for coherence length scales $\lambda_0$ of about $10^{-4} \le \lambda_0 \le 2.9 \times 10^{-3}$ Mpc.
<figure><img src="image_5.png" /></figure>

FIG. 3. 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$ yr. The points indicate the values of the peak magnetic field strengths $B_0$ and the corresponding coherence lengths $\lambda_0$, for all the points of the scan in Fig. 2 (assuming $b=0$). The color bar indicates the SNR in LISA. Points for which SNR $< 10^{-3}$ are shown in gray. The blue line depicts the limits derived in Ref. [94], where an upper limit on the magnetic field is derived when the observed $\gamma$ rays are considered as a combination of those emitted by the blazars and those stemming from cosmic-ray interactions.

However, if we consider the larger lower bound by Fermi-LAT and H.E.S.S., only a few points in our scan could account for the purported origin of primordial magnetic fields. Hence, this work exemplifies the need for continued improvement in the analysis of magnetic fields via blazar observations to exclude such model predictions, i.e., primordial seed fields sourced by the EW phase transition in the 2HDM. Since the MAGIC/Fermi-LAT constraint is the more conservative of the two, we conclude that GW interferometry can serve as an additional probe of a 2HDM origin of primordial magnetic fields over certain coherence length scales.

## IV. CONCLUSIONS

In this work, we explored the possibility that the large-scale magnetic fields observed in the intergalactic medium could have originated from a first-order electroweak phase transition in the two-Higgs-doublet model. 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 or a primordial magnetic field with negligible magnetic helicity in a plasma with kinetic helicity. We confronted these two

scenarios with observational constraints on the intergalactic magnetic fields from MAGIC, H.E.S.S., and Fermi-LAT experiments. We found that the strongest transitions, in the case of a maximally helical primordial magnetic field, could account for the magnetic fields spanning scales of about $10^{-4} \le \lambda_0 \le 2.9 \times 10^{-3}$ Mpc, which lie above the most conservative blazar bound. However, only a few points in our parameter space could exceed the more stringent bound by Fermi-LAT and H.E.S.S. We also showed that all the points, whose peak magnetic field values lie above the Fermi-LAT and MAGIC lower bound, also predict strong GW signals that could be detected by LISA, offering a multimessenger test of the 2HDM origin of primordial magnetic fields.

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 found that the maximum magnetic field strengths range from about $3 \times 10^{-17}$ to $10^{-14}$ G for all the data points and from $5 \times 10^{-16}$ to $10^{-14}$ G for the data points inside the black dashed contour, which marks the region where $B_0$ exceeds the lower limit from Fermi-LAT and MAGIC observations. Moreover, the low $t_\beta$ region can be explored in future LHC runs by searching for $gg \to A \to ZH \to \ell^+\ell^- t\bar{t}$ [39], which provides a collider probe of a possible primordial origin for the intergalactic magnetic fields within the 2HDM.

We also investigated the GW signals generated by the electroweak phase transition and found that the points with SNR $> 1$ also had the highest peak values of the magnetic field, and they all exceeded the lower bound from blazar observations by MAGIC and Fermi-LAT. However, if we considered the stronger bound by Fermi-LAT and H.E.S.S., only a few points in our scan could account for the cosmological origin of primordial magnetic fields in this work.

We concluded that a first-order electroweak phase transition within the 2HDM may explain the observed intergalactic magnetic fields over certain coherence length scales ($10^{-4} \le \lambda_0 \le 2.9 \times 10^{-3}$ Mpc), 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 multimessenger 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.

## ACKNOWLEDGMENTS

I am thankful to Shyam Balaji, Iason Baldes, Thomas Biekötter, Alberto Roper Pol, and Ramkishor Sharma for useful discussions. My work is supported by the European Union's Horizon 2020 research and innovation program under Grant Agreement No. 101002846, ERC CoG "CosmoChart."
[1] V. A. Acciari et al. (MAGIC Collaboration), Astron. Astrophys. **670**, A145 (2023).

[2] F. Aharonian et al. (H.E.S.S., Fermi-LAT Collaboration), Astrophys. J. Lett. **950**, L16 (2023).

[3] A. Neronov and I. Vovk, Science **328**, 73 (2010).

[4] L. Biermann and A. Schlüter, Phys. Rev. **82**, 863 (1951).

[5] R. Alves Batista and A. Saveliev, Universe **7**, 223 (2021).

[6] K. Dolag, M. Kachelriess, S. Ostapchenko, and R. Tomas, Astrophys. J. Lett. **727**, L4 (2011).

[7] M. S. Turner and L. M. Widrow, Phys. Rev. D **37**, 2743 (1988).

[8] B. Ratra, Astrophys. J. Lett. **391**, L1 (1992).

[9] J. Martin and J. Yokoyama, J. Cosmol. Astropart. Phys. 01 (2008) 025.

[10] T. Kobayashi, J. Cosmol. Astropart. Phys. 05 (2014) 040.

[11] T. Vachaspati, Phys. Lett. B **265**, 258 (1991).

[12] J. Ellis, M. Fairbairn, M. Lewicki, V. Vaskonen, and A. Wickens, J. Cosmol. Astropart. Phys. 09 (2019) 019.

[13] G. Sigl, A. V. Olinto, and K. Jedamzik, Phys. Rev. D **55**, 4582 (1997).

[14] A. G. Tevzadze, L. Kisslinger, A. Brandenburg, and T. Kahniashvili, Astrophys. J. **759**, 54 (2012).

[15] T. Vachaspati, Phys. Rev. Lett. **87**, 251302 (2001).

[16] C. J. Copi, F. Ferrer, T. Vachaspati, and A. Achucarro, Phys. Rev. Lett. **101**, 171302 (2008).

[17] E. Witten, Phys. Rev. D **30**, 272 (1984).

[18] C. J. Hogan, Mon. Not. R. Astron. Soc. **218**, 629 (1986).

[19] M. Kamionkowski, A. Kosowsky, and M. S. Turner, Phys. Rev. D **49**, 2837 (1994).

[20] A. Brandenburg, K. Enqvist, and P. Olesen, Phys. Rev. D **54**, 1291 (1996).

[21] M. Christensson, M. Hindmarsh, and A. Brandenburg, Phys. Rev. E **64**, 056405 (2001).

[22] T. Kahniashvili, A. Brandenburg, A. G. Tevzadze, and B. Ratra, Phys. Rev. D **81**, 123002 (2010).

[23] A. Brandenburg, T. Kahniashvili, S. Mandal, A. Roper Pol, A. G. Tevzadze, and T. Vachaspati, Phys. Rev. D **96**, 123528 (2017).

[24] Y.-Z. Chu, J. B. Dent, and T. Vachaspati, Phys. Rev. D **83**, 123530 (2011).

[25] T. Kahniashvili, A. G. Tevzadze, A. Brandenburg, and A. Neronov, Phys. Rev. D **87**, 083007 (2013).

[26] M. M. Forbes and A. R. Zhitnitsky, Phys. Rev. Lett. **85**, 5268 (2000).

[27] A. D. Sakharov, Pis'ma Zh. Eksp. Teor. Fiz. **5**, 32 (1967).

[28] H. Tashiro, W. Chen, F. Ferrer, and T. Vachaspati, Mon. Not. R. Astron. Soc. **445**, L41 (2014).

[29] W. Chen, B. D. Chowdhury, F. Ferrer, H. Tashiro, and T. Vachaspati, Mon. Not. R. Astron. Soc. **450**, 3371 (2015).

[30] A. Brandenburg, T. Kahniashvili, S. Mandal, A. Roper Pol, A. G. Tevzadze, and T. Vachaspati, Phys. Rev. Fluids **4**, 024608 (2019).

[31] V. A. Kuzmin, V. A. Rubakov, and M. E. Shaposhnikov, Phys. Lett. B **155**, 36 (1985).

[32] T. D. Lee, Phys. Rev. D **8**, 1226 (1973).

[33] G. C. Dorsch, S. J. Huber, and J. M. No, J. High Energy Phys. 10 (2013) 029.

[34] G. C. Dorsch, S. J. Huber, K. Mimasu, and J. M. No, Phys. Rev. Lett. **113**, 211802 (2014).

[35] G. C. Dorsch, S. J. Huber, T. Konstandin, and J. M. No, J. Cosmol. Astropart. Phys. 05 (2017) 052.

[36] P. Basler, M. Krause, M. Muhleitner, J. Wittbrodt, and A. Wlotzka, J. High Energy Phys. 02 (2017) 121.

[37] D. Gonçalves, A. Kaladharan, and Y. Wu, Phys. Rev. D **105**, 095041 (2022).

[38] T. Biekötter, S. Heinemeyer, J. M. No, M. O. Olea-Romacho, and G. Weiglein, J. Cosmol. Astropart. Phys. 03 (2023) 031.

[39] T. Biekötter, S. Heinemeyer, J. M. No, K. Radchenko, M. O. O. Romacho, and G. Weiglein, arXiv:2309.17431.

[40] D. Gonçalves, A. Kaladharan, and Y. Wu, Phys. Rev. D **108**, 075010 (2023).

[41] C. Caprini et al., J. Cosmol. Astropart. Phys. 03 (2020) 024.

[42] P. Auclair et al., arXiv:2204.05434.

[43] G. C. Branco, P. M. Ferreira, L. Lavoura, M. N. Rebelo, M. Sher, and J. P. Silva, Phys. Rep. **516**, 1 (2012).

[44] P. Fayet, Nucl. Phys. **B78**, 14 (1974).

[45] P. Fayet, Nucl. Phys. **B90**, 104 (1975).

[46] G. Aad et al. (ATLAS Collaboration), Phys. Lett. B **716**, 1 (2012).

[47] S. Chatrchyan et al. (CMS Collaboration), Phys. Lett. B **716**, 30 (2012).

[48] J. F. Gunion and H. E. Haber, Phys. Rev. D **67**, 075019 (2003).

[49] H. Bahl, T. Biekötter, S. Heinemeyer, C. Li, S. Paasch, G. Weiglein, and J. Wittbrodt, Comput. Phys. Commun. **291**, 108803 (2023).

[50] G. Aad et al. (The ATLAS Collaboration) Nature (London) **607**, 52 (2022); **612**, E24 (2022).

[51] A. Tumasyan et al. (CMS Collaboration), Nature (London) **607**, 60 (2022).

[52] A. Masiero, and G. Senjanovic, and A. Smirnov, in Proceedings, Summer School in High-energy physics and cosmology (World Scientific, Singapore, 1999), pp. 187–259.

[53] S. R. Coleman and E. J. Weinberg, Phys. Rev. D **7**, 1888 (1973).

[54] L. Dolan and R. Jackiw, Phys. Rev. D **9**, 3320 (1974).

[55] P. B. Arnold and O. Espinosa, Phys. Rev. D **47**, 3546 (1993); **50**, 6662(E) (1994).

[56] D. Curtin, P. Meade, and H. Ramani, Eur. Phys. J. C **78**, 787 (2018).

[57] D. Croon, O. Gould, P. Schicho, T. V. I. Tenkanen, and G. White, J. High Energy Phys. 04 (2021) 055.

[58] O. Gould and T. V. I. Tenkanen, J. High Energy Phys. 06 (2021) 069.

[59] O. Gould and T. V. I. Tenkanen, arXiv:2309.01672.

[60] S. R. Coleman, Phys. Rev. D **15**, 2929 (1977); **16**, 1248(E) (1977).

[61] Curtis G. Callan Jr. and S. R. Coleman, Phys. Rev. D **16**, 1762 (1977).

[62] A. D. Linde, Phys. Lett. B **100**, 37 (1981).

[63] A. D. Linde, Nucl. Phys. B**216**, 421 (1983); B**223**, 544(E) (1983).

[64] C. L. Wainwright, Comput. Phys. Commun. **183**, 2006 (2012).

[65] D. Bodeker and G. D. Moore, J. Cosmol. Astropart. Phys. 05 (2009) 009.

[66] D. Bodeker and G. D. Moore, J. Cosmol. Astropart. Phys. 05 (2017) 025.
[67] G. D. Moore and T. Prokopec, Phys. Rev. D **52**, 7182 (1995).

[68] T. Konstandin, G. Nardini, and I. Rues, J. Cosmol. Astropart. Phys. 09 (2014) 028.

[69] J. Kozaczuk, J. High Energy Phys. 10 (2015) 135.

[70] G. C. Dorsch, S. J. Huber, and T. Konstandin, J. Cosmol. Astropart. Phys. 12 (2018) 034.

[71] M. Barroso Mancha, T. Prokopec, and B. Swiezewska, J. High Energy Phys. 01 (2021) 070.

[72] B. Laurent and J. M. Cline, Phys. Rev. D **102**, 063516 (2020).

[73] G. C. Dorsch, S. J. Huber, and T. Konstandin, J. Cosmol. Astropart. Phys. 04 (2022) 010.

[74] W.-Y. Ai, B. Garbrecht, and C. Tamarit, J. Cosmol. Astropart. Phys. 03 (2022) 015.

[75] M. Lewicki, M. Merchand, and M. Zych, J. High Energy Phys. 02 (2022) 017.

[76] B. Laurent and J. M. Cline, Phys. Rev. D **106**, 023501 (2022).

[77] W.-Y. Ai, B. Laurent, and J. van de Vis, J. Cosmol. Astropart. Phys. 07 (2023) 002.

[78] S. Jiang, F. P. Huang, and X. Wang, Phys. Rev. D **107**, 095005 (2023).

[79] S. Dimopoulos and L. Susskind, Phys. Rev. D **18**, 4500 (1978).

[80] P. Amaro-Seoane et al. (LISA Collaboration), arXiv:1702.00786.

[81] P. Auclair et al. (LISA Cosmology Working Group Collaboration), Living Rev. Relativity **26**, 5 (2023).

[82] P. Auclair, C. Caprini, D. Cutting, M. Hindmarsh, K. Rummukainen, D. A. Steer, and D. J. Weir, J. Cosmol. Astropart. Phys. 09 (2022) 029.

[83] Science requirements document, https://www.cosmos.esa.int/web/lisa/lisa-documents, Technical Report No. ESA-L3-EST-SCI-RS-001, 2018.

[84] A. Armua, A. Berera, and J. C. Figueroa, Phys. Rev. E **107**, 055206 (2023).

[85] R. Banerjee and K. Jedamzik, Phys. Rev. D **70**, 123003 (2004).

[86] D. Biskamp and W.-C. Müller, Phys. Rev. Lett. **83**, 2195 (1999).

[87] A. Roper Pol, A. Neronov, C. Caprini, T. Boyer, and D. Semikoz, arXiv:2307.10744.

[88] J. R. Espinosa, T. Konstandin, J. M. No, and G. Servant, J. Cosmol. Astropart. Phys. 06 (2010) 028.

[89] T. Kahniashvili, A. G. Tevzadze, and B. Ratra, Astrophys. J. **726**, 78 (2011).

[90] R. Durrer and A. Neronov, Astron. Astrophys. Rev. **21**, 62 (2013).

[91] T. Vachaspati, Phys. Rev. D **95**, 063505 (2017).

[92] T. Vachaspati, Rep. Prog. Phys. **84**, 074901 (2021).

[93] D. Biskamp, Plasma Phys. Controlled Fusion **45**, 1827 (2003).

[94] W. Essey, S. Ando, and A. Kusenko, Astropart. Phys. **35**, 135 (2011).