# Magnetogenesis with gravitational waves and primordial black hole dark matter

Shyam Balaji, Malcolm Fairbairn, Maria Olalla Olea-Romacho

► To cite this version:

Shyam Balaji, Malcolm Fairbairn, Maria Olalla Olea-Romacho. Magnetogenesis with gravitational waves and primordial black hole dark matter. Physical Review D, 2024, 109 (7), pp.075048. 10.1103/PhysRevD.109.075048. hal-04467898

HAL Id: hal-04467898

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

Submitted on 29 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
Magnetogenesis with gravitational waves
and primordial black hole dark matter

Shyam Balaji<sup>1,*</sup> Malcolm Fairbairn,<sup>1,†</sup> and María Olalla Olea-Romacho<sup>2,‡</sup>

<sup>1</sup>Physics Department, King's College London, Strand, London WC2R 2LS, United Kingdom

<sup>2</sup>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 16 February 2024; accepted 11 April 2024; published 30 April 2024)

Strongly supercooled first-order phase transitions (FOPTs) can produce primordial black hole (PBH) dark matter (DM) along with observable gravitational waves (GWs) from bubble collisions. Such FOPTs may also produce coherent magnetic fields generated by bubble collisions and by turbulence in the primordial plasma. Here we find that the requirement for PBH DM can produce large primordial magnetic fields which subsequently yield intergalactic magnetic fields in the present Universe (with magnitude ≲20 pG across coherence length scales of ≈0.001–0.01 Mpc, assuming maximally helical magnetic fields) that easily exceed lower bounds from blazar observations. We follow a largely model-independent approach and highlight the possibility of producing DM and observable multimessenger magnetic fields and GW signals visible in next generation experiments.

DOI: 10.1103/PhysRevD.109.075048

I. INTRODUCTION

Dark matter (DM) remains one of the most important open questions in physics and cosmology. One popular possibility is for all of the DM to be composed of primordial black holes (PBHs) [1–4] if their masses are in the asteroid-mass range [5,6]. It has been demonstrated that PBHs can be copiously generated in the supercooling regime of a first-order phase transition (FOPT) when the energy density of the Universe is controlled by the latent heat of a phase transition [7,8]. The latent heat behaves as a cosmological constant which drives the expansion of the Universe until the transition ends and the energy is transformed into radiation when bubbles form and merge, subsequently reheating the plasma. Stochastically late-nucleated Hubble patches lead to horizon-scale overdensities that collapse into PBHs upon completion of the transition.

Primordial magnetogenesis has also been studied in
relation to FOPTs, such as, e.g., the electroweak (EW)
phase transition [9–11] or the QCD phase transition
[12,13]. Notably, the idea of magnetogenesis during a

first-order EW phase transition was first proposed in Ref. [9]. In this scenario, magnetic fields are created by EW sphaleron decays [14,15]. As the bubbles formed during the EW phase transition grow, collide, and merge, they cause the primordial plasma to move at high Reynolds number, and the magnetic fields enter a state of magnetohydrodynamic (MHD) turbulence [16–22].

Coherent intergalactic magnetic fields (IGMFs) are indirectly suggested by blazar observations [23–25]. The origin of these IGMFs is a long-standing mystery, with two main possibilities: astrophysical or cosmological. However, the uniformity and large-scale correlation of magnetic fields observed in cosmic voids are difficult to account for solely by astrophysical sources. This is because the mechanisms by which these fields could spread from galactic sources into the voids and maintain a near-uniform strength across vast distances are not well understood. However, an example of astrophysical sources are weak initial magnetic fields from local effects in astrophysical objects (e.g., the Biermann battery mechanism [26]) that are later enhanced by dynamo effects [27], which can create long-correlated magnetic fields. Cosmological explanations from FOPTs can also easily accommodate magnetic fields with very large coherence lengths. Hence, potential sources from early Universe cosmological processes such as FOPTs are very attractive.

In this paper, we will use model-independent results for the inverse duration of the FOPT duration and reheating temperature that produce PBHs as all DM, we then compute the corresponding IGMFs in the present Universe and compare them with blazar lower bounds.

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>shyam.balaji@kcl.ac.uk</footnote>

<footnote>malcolm.fairbairn@kcl.ac.uk</footnote>

<footnote>mariaolalla.olearomacho@phys.ens.fr</footnote>
Strong FOPTs also produce a stochastic background of gravitational waves (GWs) [16,17,28]. By enforcing the production of PBH DM in the transition, we also compute the corresponding GW signal from the bubble collisions, leading to observable signals in upcoming GW interferometers well above stochastic astrophysical foregrounds [29]. In this way, we highlight the relation between multimessenger cosmological signals of shared origin in two disparate channels, magnetic fields and GWs, while simultaneously providing an explanation for DM. We also briefly remark on a possible model realization via the classically conformal $U(1)_{B-L}$ model for completeness.

This paper is structured as follows: In Sec. II, we provide (i) a brief overview of PBH formation from FOPTs, (ii) we show how to compute the resulting primordial magnetic field and coherence length and outline IGMF bounds from blazars, and (iii) we outline how to calculate GW spectrum and verify observability in various experiments. In Sec. III, we discuss the main results of this work and finally we conclude in Sec. IV.

## II. PRIMORDIAL BLACK HOLES, MAGNETIC FIELDS, AND GRAVITATIONAL WAVES

### A. Primordial black hole formation

In a FOPT, the Universe can change from a metastable symmetric vacuum to a symmetry-breaking vacuum, through a process of bubble nucleation, growth, and merging. Supercooling occurs when the Universe remains in the metastable phase for an extended period of time, such that the expansion of the Universe is largely driven by the false vacuum energy rather than radiation. The vacuum energy inflates the Universe until it converts to radiation upon bubble coalescence, ultimately reheating the plasma to a temperature $T_{reh}$. Given the stochastic nature of the phase transition, PBHs can form in supercooled cosmological phase transitions. For the average background Hubble patch, nucleation occurs at the cosmic time $\tau_{nuc}$. Causal regions (labeled by i) that nucleate at a time $t_{nuc}^i$ later than $\tau_{nuc}$ remain vacuum dominated for a longer period. The false vacuum energy does not change with the expansion of the Universe, while the energy density of both the bubble walls and the radiation decreases rapidly. For this reason, late-nucleating patches quickly become overdense with respect to the surrounding background patches, and if the overdensity exceeds a certain threshold $\delta_c$, these regions may collapse into PBHs. We use the density contrast in radiation to determine the collapse condition as

$$ \delta(t) \equiv \frac{\rho_r^{\text{late}}(t; t_{\text{nuc}}^i) - \rho_r^{\text{bkg}}(t)}{\rho_r^{\text{bkg}}(t)} > \delta_c, \quad (1) $$

where $\rho_r^{\text{late}}(t; t_{nuc}^i)$ is the radiation density in a region where nucleation is delayed to $t_{nuc}^i$ and $\rho_r^{\text{bkg}}(t)$ is the radiation energy density of a background Hubble patch at cosmic

time $t$. The density contrast peaks shortly after the late patch percolation, as the energy density in the background patch begins diluting slightly earlier, while the late patch keeps a constant energy density due to the vacuum energy density.

The inverse timescale of the transition is given by

$$ \frac{\beta}{H} \equiv \frac{1}{H\Gamma_{\text{bub}}} \frac{d\Gamma_{\text{bub}}}{dt}, \quad (2) $$

where $\Gamma_{\text{bub}}$ is the bubble nucleation rate per unit volume and $H$ is the Hubble rate during the phase transition.<fnref n="1" /> The Hubble rate in the false vacuum reads

$$ H \equiv H_{\text{false}} = \sqrt{\frac{8\pi}{3M_{\text{Pl}}^2} \rho_{\text{vac}}}, \quad (3) $$

where $\rho_{vac}$ is the vacuum energy density just before the onset of the phase transition.

In order to quantify the temperature at the moment that the phase transition completes, and assuming that the Universe rapidly reheats to the reheating temperature at a rate much greater than the expansion rate of the Universe. We define the reheating temperature as

$$ T_{\text{reh}} = \left( \frac{30\rho_{\text{vac}}}{g_r(T_{\text{reh}})\pi^2} \right)^{1/4}, \quad (4) $$

where $g_r$ are the radiation degrees of freedom. This temperature should be understood as the temperature in the false vacuum patch, which does not contain nucleated bubbles or partially reheated plasma.

For a monochromatic distribution, the mass of these PBHs is given by the energy inside the sound horizon at the time of the collapse [30]

$$ M_{\text{PBH}}(T_{\text{reh}}) = M_{\odot} \left( \frac{20}{g_r(T_{\text{reh}})} \right)^{1/2} \left( \frac{140 \text{ MeV}}{T_{\text{reh}}} \right)^2, \quad (5) $$

whereas the PBH abundance normalized to the DM relic abundance reads [30]

$$ f_{\text{PBH}}(\beta/H, T_{\text{reh}}) = \frac{e^{-a_1(\beta/H)^{a_2}}(1+\delta_c)^{a_3(\beta/H)_*}}{2.2 \times 10^{-8}} \frac{T_{\text{reh}}}{0.14}. \quad (6) $$

Here, $a_1 = 0.56468$, $a_2 = 1.266$, and $a_3 = 0.6639$, and $\beta/H$ is evaluated at the percolation time, yielding

$$ (\beta/H)_* = \frac{a_2 W_0 \left[ \frac{a_3}{a_2} \left( \frac{19.5983 + \log T_{\text{reh}}}{a_1} \right)^{\frac{1}{a_2}} \log(1 + \delta_c) \right]}{a_3 \log(1 + \delta_c)}, \quad (7) $$

<footnote n="1">¹In the following, we will use a first-order Taylor expansion to approximate the bubble nucleation rate near $\tau_{nuc}$. In Ref. [29], it was shown that the approximate form yields accurate predictions around the nucleation time.</footnote>
where $W_0$ is the zeroth-order Lambert W function. Calculations, based on full general relativistic simulations, often suggest a range of 0.4–0.66 for $\delta_c$ in the context of overdensities from inflation reentering the Hubble horizon [31–39]. Deviations from these assumptions are expected in the phase transition scenario. However, in this study, we align with previous approaches [29,30] and choose $\delta_c = 0.45$ to enable comparison. The precise value of $\delta_c$ impacts $f_{PBH} = 1$, yet the strong magnetic field generated remain insensitive to this choice.

# 1. Experimental constraints for PBHs as the dark matter

Observationally, $f_{PBH} = 1$ is allowed for nearly monochromatic distributions of PBHs with masses falling within what is known as the asteroid-mass window [40]

$$10^{-16} M_{\odot} \lesssim M_{PBH} \lesssim 10^{-10} M_{\odot}. \qquad (8)$$

Here the lower bound is set by the particle flux from Hawking evaporation, which would affect the cosmic microwave background (CMB) [41] or be detected by the Voyager [42] probe. The upper bound is set by microlensing measurements [43]. According to Eq. (5), the range of reheating temperatures that gives rise to PBHs with masses within the asteroid-mass window is

$$10 \lesssim T_{\text{reh}} \lesssim 10^4 \text{ TeV}. \qquad (9)$$

To get $f_{PBH} = 1$ in this temperature range, we use Eq. (7).

# B. Primordial magnetic fields

## 1. Magnetic field spectrum today

The initial conditions of the magnetic field and the plasma, which are not well understood, have a significant impact on how MHD turbulence decays. For instance, the initial seed field magnetic helicity is one of the factors that influences the decay of MHD turbulence [44]. The magnetic helicity is a quantity that describes how twisted and linked the magnetic field lines are. It is given by $\langle \mathbf{A} \cdot \mathbf{B} \rangle$, where $\mathbf{B} = \nabla \times \mathbf{A}$, where $\mathbf{A}$ and $\mathbf{B}$ are the usual vector potential and magnetic field, respectively. The magnetic helicity is nearly conserved in a highly conductive medium. This implies that a field with maximum helicity has to increase its correlation length as it loses magnetic energy, which leads to an inverse cascade of magnetic energy where the energy moves from smaller scales to larger scales, creating coherent magnetic structures at scales much bigger than the ones where the energy was initially injected. Therefore, this process could have been very important for the survival and evolution of primordial magnetic fields, which would be correlated at very large length scales today. Helical magnetic fields originate from changes in the Chern-Simons number, a process associated with the violation of the combined baryon and lepton number.

This phenomenon typically occurs during the EW phase transition, in which sphaleron decays produce helical magnetic fields, as supported by numerical simulations and theory [14,15,45].

The MHD turbulence for maximally helical magnetic fields decays as a power law in conformal time $\eta$, with the following relations for the comoving magnetic field and correlation length during the radiation-dominated epoch [46]:

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

This relation is applicable for a highly conductive plasma in which magnetic helicity is conserved [10], illustrating that, at a specific scale, the magnetic field strength decreases with conformal time as magnetic energy is redistributed from smaller to larger scales, increasing the coherence length. It has been shown by numerical simulations that magnetic fields with zero or negligible net helicity can also experience an inverse transfer of magnetic energy when there is a plasma with initial kinetic helicity [47].<fnref n="2" /> In this scenario, comoving magnetic field and correlation length, $\lambda$ scale as

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

respectively, during the radiation-dominated epoch. These scaling laws are valid during the epoch of radiation domination, when the scale factor increases linearly with conformal time, i.e., $a \sim \eta$. After recombination, the magnetic field decreases like radiation, i.e., $B \sim a^{-2}$. To express these two scenarios in a concise way, we introduce the parameters

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

for the power laws

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

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

The magnetic field energy density at percolation can be estimated by [10,49]

$$\rho_{B,*} = 0.1 \frac{\kappa_{\text{col}} \alpha}{1 + \alpha} \rho_{*} \approx \frac{\pi^2}{3} T_{\text{reh}}^4 \approx 0.1 \rho_{\text{vac}}. \qquad (14)$$

<footnote n="2">²However, this has been challenged recently in Ref. [48], where a weaker inverse cascade of magnetic energy for non-helical fields was found, compared to previous studies in the literature.</footnote>
Here, $\rho_* = \frac{3M_{\text{pl}}^2}{8\pi} H_*^2 = \frac{g_r(T_{\text{reh}})\pi^2}{30} T_{\text{reh}}^4 \approx \rho_{\text{vac}}$ is the total energy density at the percolation temperature. Here $\kappa_{\text{col}}$ is the fraction of the released vacuum energy that is used in accelerating the bubbles. The efficiency for converting bulk fluid motion of the plasma into magnetic fields via MHD turbulence was assumed to be of 10% [22,50,51]. For supercooled phase transitions, we can safely assume $\alpha \to \infty$ and $\kappa_{\text{col}} = 1$, leading to the approximation in Eq. (14).

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

$$B_0(\lambda) = B(\lambda, t_0) \\ = \left(\frac{a_{\text{reh}}}{a_{\text{rec}}}\right)^{p_b/2} \left(\frac{a_{\text{reh}}}{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} \quad (15)$$

which assumes a power-law spectrum for the magnetic field strength, with a spectral index of $n = 2$ at large scales. For the magnetic field evolution above, we do not consider the slope of the original spectrum at length scales smaller than $\lambda_0$, as we generically expect turbulence to fully develop relatively quickly at these smaller length scales, driven by turbulent eddies at this coherence scale. However, for the subsequent initial spectrum (that evolves with time to $B_0$) that emerges from the original spectrum wherein the turbulent eddies have developed on all scales $\lambda < \lambda_0$, magnetic and kinetic energy are subsequently transferred to increasingly smaller scales via a direct energy cascade, until the energy is eventually lost as heat at the dissipation scale of the plasma. Technically, the power-law spectrum should cut off at the dissipation scale. However, we omit this from (15) since physics beyond the dissipation scale is irrelevant for our purposes.

We note $\lambda_0$ denotes the field coherence scale redshifted to today [10],

$$\lambda_0 = \lambda_B(t_0) = \left(\frac{a_{\text{rec}}}{a_{\text{reh}}}\right)^{q_b} \left(\frac{a_0}{a_{\text{reh}}}\right) \lambda_*, \quad (16)$$

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

$$\lambda_* = R_* = \frac{(8\pi)^{1/3}}{H_*} \left(\frac{\beta}{H}\right)^{-1}, \quad (17)$$

where, for the case of supercooled phase transitions, a bubble wall velocity of $v_w = 1$ was assumed.

Taking $v_w \to 1$ is justified in a strongly supercooled transition, because the pressure difference between the false vacuum (the metastable state) and the true vacuum (the stable state) becomes very large. This large pressure

difference results in a significant force on the bubble walls, causing them to accelerate rapidly. As the Universe expands and cools down, the potential energy difference between the two vacua increases, leading to a larger driving force for the bubble walls to expand. In the absence of significant friction or other retarding forces, the bubble walls can reach relativistic speeds as they move to convert the false vacuum into the true vacuum. The peak value of the magnetic field spectrum is denoted as $B_0 = B_0(\lambda_0)$. The redshift factors are computed as

$$\frac{a_{\text{reh}}}{a_0} = 8 \times 10^{-14} \left( \frac{100}{g_r(T_{\text{reh}})} \right)^{1/3} \left( \frac{1 \text{ GeV}}{T_{\text{reh}}} \right), \quad (18)$$

$$\frac{a_{\text{reh}}}{a_{\text{rec}}} = 8 \times 10^{-11} \left( \frac{100}{g_r(T_{\text{reh}})} \right)^{1/3} \left( \frac{1 \text{ GeV}}{T_{\text{reh}}} \right), \quad (19)$$

where $g_r(T_{\text{reh}})$ corresponds to the total number of relativistic degrees of freedom in entropy at the reheating temperature $T_{\text{reh}}$. For this work, we will assume $g_r(T_{\text{reh}}) = 106.75$ as in the case of the Standard Model (SM). However, this value is model dependent and would have to be adjusted to reflect the dynamics of the UV complete theory.

If the transition is driven by dark scalars, coupled through a portal to the SM Higgs field, a small fraction of the bubble wall energy can be stored in the Higgs field bubble. We take into account this suppression by defining an efficiency parameter $\kappa_h$ which indicates how much energy has been transferred to the SM Higgs field.

## 2. Experimental constraints from blazar emissions

Blazars are active galactic nuclei jets that are roughly pointed toward us, and they emit TeV $\gamma$ rays. These $\gamma$ rays can collide with photons from the extragalactic background, creating $e^+e^-$ pairs [53]. Then, these pairs can interact with photons from the CMB, producing $\gamma$ rays with energies in the GeV range. This in turn modifies the original spectrum of the blazars by decreasing the number of TeV $\gamma$ rays and increasing the number of GeV $\gamma$ rays [54,55]. In the presence of IGMF, the $e^+e^-$ pairs can be deflected by the Lorentz force and, if the field is strong enough, the final GeV photons are no longer directed toward the observer. The absence 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, done by the MAGIC $\gamma$-ray observatory, sets a lower bound of $B > 1.8 \times 10^{-17}$ G for correlation lengths $\lambda \ge 0.2$ Mpc [23]. 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}$ G for coherence lengths of $\lambda \ge 1$ Mpc, assuming a blazar duty cycle of $t_d = 10$ yr [24].
### C. Gravitational waves

We use the results from hybrid simulations [56] to estimate the anisotropic stress from a bubble collision, where the authors first simulated the collision in one spatial dimension, then used the result as a source at points at which walls collide in a three-dimensional lattice simulation with thin walls. This way, the effects of scalar gradients and gauge field production during the collision can be captured by the lower-dimensional simulation. The spectrum is given by

$$h^2\Omega_{gw}(f) = 5.10 \times 10^{-9} \left(\frac{100}{g_r(T_{reh})}\right)^{1/3} \left(\frac{10}{(\beta/H)_*}\right)^2 S_{hy}(f), \quad (20)$$

where the shape is

$$S_{hy}(f) = \frac{695}{\left[2.41 \left(\frac{f}{f_{hy}}\right)^{-0.557} + \left(\frac{f}{f_{hy}}\right)^{0.574}\right]^{4.20}} \quad (21)$$

and peak frequency

$$\tilde{f}_{hy} = 22 \text{ mHz} \left(\frac{g_r(T_{reh})}{100}\right)^{1/6} \left(\frac{(\beta/H)_*}{10}\right) \left(\frac{T_{reh}}{10^2 \text{ TeV}}\right). \quad (22)$$

The signal-to-noise ratio is given by [57–61]

$$\text{SNR} = \sqrt{t_{\text{obs}} \int \left( \frac{\Omega_{\text{gw}}^2}{\Omega_{\text{sens}}^2 + 2\Omega_{\text{gw}}\Omega_{\text{sens}} + 2\Omega_{\text{gw}}^2} \right) df}, \quad (23)$$

where $t_{\text{obs}}$ is the observation time and $\Omega_{\text{sens}}$ is the sensitivity of the interferometer. The predictions for the GWs produced by the supercooled FOPT that we consider yield compelling phenomenology in the observational window of Laser Interferometer Space Antenna (LISA), Einstein Telescope (ET), Cosmic Explorer (CE), and Deci-hertz Interferometer Gravitational Wave Observatory (DECIGO). At higher frequency, we have a band that would be observable with the ET [62], which would be sensitive to the range $10–10^3$ Hz. In between Square Kilometre Array [63] and ET frequencies, we expect LISA [64] and DECIGO [65–67] to be applicable. In order to compare theoretical predictions with experimental projections, we use the recently derived future sensitivity curves for these GW limits using the experimental design specification along with the peak integrated sensitivity curves computed in Ref. [68] and assume an observation time of 1 yr.

## III. RESULTS

Assuming a supercooled phase transition that saturates the DM content of the Universe as PBHs, we can compute the properties of the primordial magnetic fields created during the transition, such as their strength and coherence

<figure><img src="image_3.png" /><figcaption>FIG. 1. Magnetic field peak value (blue line) in the present Universe $B_0$ vs the coherence length $\lambda_0$ for a reheating temperature of $T_{reh} \in [10^4, 10^7]$ GeV, $b = 0$, and $f_{pbh} = 1$. The solid red line corresponds to the lower bound on the magnetic field set by Ref. [23], whereas the dashed red line corresponds to that set by Ref. [24] for a duty cycle of $t_d = 10$ yr.</figcaption></figure>

length. We can also evaluate the GW signal generated in this transition.

In Fig. 1, we show the peak magnetic field amplitude $B_0$, computed using (15), for a maximally helical magnetic field ($b=0$), as a function of the corresponding coherence length $\lambda_0$ at the peak magnetic field value, from (16), spanning over reheating temperatures of $[10^4, 10^7]$ GeV. The coherence length scales for combinations of reheating temperatures and $\beta/H$ values that produce $f_{pbh} = 1$ lie in a relatively narrow window from $6 \times 10^{-4}$ to $7 \times 10^{-3}$ Mpc. Across this range, the magnetic field amplitude varies from around 1 to 15 pG. In the same figure, we show the peak magnetic field for an efficiency factor of $\kappa_h = 10^{-9}$ (note $B_0$ scales like $\sqrt{\kappa_h}$) to illustrate the possible effect of bubble wall energy being stored in the SM Higgs field if the transition were driven by dark scalar through a portal to the Higgs field. For strongly supercooled transitions, $\kappa_h = 10^{-9}$ is an estimate of the possible suppression inspired in scale-invariant extensions of the SM [10] (see Sec. III A). The red solid line indicates the lower bound on $B_0$ from MAGIC and Fermi-LAT observations at the corresponding coherence length $\lambda_0$, whereas the dashed line corresponds to the one based on data from the Fermi-LAT and H.E.S.S. Collaborations. In the absence of suppression ($\kappa_h = 1$), the values of $B_0$ could explain the blazar observations, since the predicted values of $B_0$ easily exceed the lower bounds set by the experiments. Taking into account illustrative examples for the suppression featuring values of $\kappa_h$ as low as $10^{-9}$, we obtain values of $B_0$ that can exceed the lower bound set by H.E.S.S. and Fermi-LAT but are too weak to surpass the limits established by MAGIC and Fermi-LAT over most of the range of coherence lengths. However, the values for $B_0$ are still above the most conservative blazar bound set by H.E.S.S. and Fermi-LAT over coherence lengths of $\approx 3.4 \times 10^{-3}$ Mpc.
We note that, in the case of a nonhelical magnetic field ($b=1$), the peak magnetic field gets suppressed by more than 2 orders of magnitude and lies in the range $2 \times 10^{-3}-7 \times 10^{-2}$ pG over smaller coherence length scales of $\simeq 9 \times 10^{-7}-3 \times 10^{-5}$ Mpc. This is of significantly less experimental interest, so we focus our discussion on the more relevant, maximally helical case. The maximally helical case is also phenomenologically interesting since a cosmological origin is plausible for correlated magnetic fields over large scales. Such mechanisms frequently involve an inverse cascade process that requires an initial net magnetic helicity. This is because the correlation length of the fields generated in the early Universe cannot exceed the causal region at that time, so a way to transfer the magnetic energy to larger scales is needed. The inverse cascade, which occurs for helical fields, is a natural candidate for this. In the following, we choose to disregard the case of nonhelical magnetic fields ($b=1$), as the magnetic field amplitudes in this scenario are relatively suppressed.

We also plot the coherence length and magnetic field amplitude as a function of the PBH DM mass, as shown in Fig. 2. In the top and bottom panels, we show the peak magnetic field amplitude and coherence length $\lambda_0$ as a function of the PBH mass, respectively, computed using (5), which spans the aforementioned asteroid-mass range of $[10^{-16}, 10^{-10}] M_\odot$ for the reheating temperatures considered. PBHs as DM are correlated with coherent magnetic fields at length scales from around 0.001 to 0.01 Mpc.

<figure><img src="image_4.png" /><figcaption>FIG. 2. Magnetic field peak value in the present Universe $B_0$ vs the primordial black hole mass (top). Coherence length $\lambda_0$ vs the primordial black hole mass (bottom) is also shown. In both cases, we consider a reheating temperature of $T_{\text{reh}} \in [10^4, 10^7]$ GeV, $b=0$, and $f_{\text{pbh}} = 1$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIG. 3. Gravitational wave signal-to-noise ratio as a function of the magnetic field peak value in the present Universe. We show the signal-to-noise ratio for LISA (blue), ET (purple), DECIGO (gray), and CE (cyan). We consider a reheating temperature of $T_{\text{reh}} \in [10^4, 10^7]$ GeV, $b=0$, and $f_{\text{pbh}} = 1$.</figcaption></figure>

In Fig. 3, we show the signal-to-noise ratio from GWs produced by a supercooled phase transition against the peak magnetic field amplitude. We use the respective sensitivity curves for the LISA, ET, DECIGO, and CE high frequency GW experiments described in Sec. II C along with those predicted by the supercooled transition using (23) to calculate the signal-to-noise ratio. Future runs of experiments such as LIGO have a lower sensitivity than other experiments that cover the same subset of the frequency range such as ET, hence the prospects for the latter are more relevant. Remarkably, we see that the signal-to-noise ratios are around $\simeq 4000$ for DECIGO and span between around 1000 and 3000 for LISA. For ET, the signal-to-noise ratio is around 3000 for a magnetic field of 2 pG, but drops to unity at 4 pG. For CE, the signal-to-noise ratio is around 4000 for a magnetic field of 2 pG, but drops to unity around 6 pG. We therefore observe the potential for a tantalizing multimessenger signal yielding detectable GWs while producing a correlated cosmological magnetic field and all the requisite DM content of the Universe in the form of PBHs.

### A. A model example: The classically conformal $B-L$ model

To provide a concrete realization, we consider EW symmetry breaking in the classically scale-invariant $U(1)_{B-L}$ extension of the SM, where we introduce three right-handed neutrinos $N_i$, with $Q_{B-L} = -1$, and a complex scalar $\rho$, with $Q_{B-L} = -2$, which breaks the $B-L$ gauge symmetry. The radiative breaking of the $B-L$ symmetry induces a negative mass term for the EW Higgs field $H$ via a portal coupling $\lambda_{\rho H}$, thus eventually breaking EW symmetry. It is during this process that magnetic fields can be generated due to Higgs mean-field
gradients [55]. Subsequently, plasma effects restore the EW symmetry and, finally, the Universe cools down postphase transition reheating and the EW symmetry breaks in a similar fashion to the SM, i.e., through a crossover. The kinetic term for $\rho$ is given by $(D_\mu\rho)^*(D^\mu\rho)$, where $D_\mu$ includes the gauge coupling $g_{B-L}$, leading to a mass for the $B-L$ gauge boson upon $\rho$ acquiring a vacuum expectation value. The scalar potential is described by

$$V(H, \rho) = \lambda_{\rho}|\rho|^4 + \lambda_{\rho h}|\rho|^2|H|^2 + \lambda_h|H|^4, \quad (24)$$

and the Yukawa interactions include terms $y_{\nu ij} \overline{l}_{Li} \tilde{H} N_j + \frac{1}{2} y_{Ni} \rho \tilde{N} N^c + \text{H.c.}$, with $l_{Li}$ being the SM lepton doublets and $\tilde{H}$ defined as $i\sigma_2 H^*$. For further details on the model and its finite temperature treatment in cosmological phase transitions, see, e.g., Ref. [29].

Considering a $Z'$ with a mass of $m_{Z'} = 56.64$ TeV and $B-L$ gauge coupling of $g_{B-L} = 0.2832$: we get a nucleation temperature of $T_{\text{nuc}} = 2.03$ GeV, a reheating temperature of $T_{\text{reh}} = 5.03$ TeV, and the inverse duration of the transition $\beta/H = 6.81$. In this context, the predicted PBH abundance and PBH mass are

$$f_{\text{PBH}} = 1, \quad (25)$$

$$M_{\text{PBH}} = 5 \times 10^{-10} M_{\odot}. \quad (26)$$

At the same time, the peak strength of the generated magnetic fields is $B_0 = 1.9 \times 10^{-11}$ G at a peak coherence length of $\lambda_0 = 0.008$ Mpc. In this case we have estimated $\kappa_h$ as $\kappa_h = \Delta V_h / \Delta V_\phi$, where $\Delta V_\phi = V(v_\phi, 0) - V(0, 0)$ and $V_\phi = V(v_\phi, 0) - V(v_\phi, v_h)$. $\kappa_h$ can be understood as the ratio between the heights of the scalar potential in each of the two field directions. Including the effect of $\kappa_h \simeq 4.8 \times 10^{-9}$ yields a magnetic field strength of $B_0 = 1.3 \times 10^{-15}$ G at a coherence length scale $\lambda_0 = 0.008$ Mpc.

## IV. CONCLUSION

We investigate the intriguing connection between strongly supercooled FOPTs, PBH DM generation, and the simultaneous production of primordial magnetic fields with observable consequences. Assuming a FOPT that saturates the DM content of the Universe through PBHs, we explore the resulting properties of primordial magnetic fields, including their strength and coherence length.

We compute the Universe's peak magnetic field amplitude and its coherence length, considering reheating temperatures within the range of $10^4-10^7$ GeV and PBH DM masses in the range $[10^{-16}, 10^{-10}] M_\odot$. We observe a narrow window of coherence length scales, approximately 0.001-0.01 Mpc, producing magnetic field amplitudes ranging from around 1 to 15 pG (assuming maximally helical magnetic fields). This easily exceeds lower limits on intergalactic magnetic fields from blazar observations.

Finally, we compute the signal-to-noise ratio from GWs produced by a supercooled phase transition against the peak magnetic field amplitude. Calculations using sensitivity curves for various GW experiments demonstrate substantial signal-to-noise ratios, reaching around $\gtrsim 10^3$ for DECIGO, Laser Interferometer Space Antenna, Einstein Telescope, and Cosmic Explorer. This presents the exciting prospect of a multimessenger signal, featuring detectable GWs alongside the correlated generation of cosmological magnetic fields and the formation of PBHs as the dominant DM component in the Universe.

## ACKNOWLEDGMENTS

M. O. O. R. thanks I. Baldes, K. Kamada, and T. Vachaspati, while M. F. thanks V. Vaskonen for useful discussions. S. B. and M. F. are supported by the STFC under Grant No. ST/X000753/1. M. O. O. R. is supported by the European Union's Horizon 2020 research and innovation program under Grant Agreement No. 101002846, ERC CoG "CosmoChart."

[1] Y. B. Zel'dovich and I. D. Novikov, Sov. Astron. AJ (Engl. Transl.) **10**, 602 (1967).

[2] S. Hawking, Mon. Not. R. Astron. Soc. **152**, 75 (1971).

[3] B. J. Carr and S. W. Hawking, Mon. Not. R. Astron. Soc. **168**, 399 (1974).

[4] B. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama, Rep. Prog. Phys. **84**, 116902 (2021).

[5] S. Balaji, J. Silk, and Y.-P. Wu, J. Cosmol. Astropart. Phys. **06** (2022) 008.

[6] W. Qin, S. R. Geller, S. Balaji, E. McDonough, and D. I. Kaiser, Phys. Rev. D **108**, 043508 (2023).

[7] H. Kodama, M. Sasaki, and K. Sato, Prog. Theor. Phys. **68**, 1979 (1982).

[8] S. D. H. Hsu, Phys. Lett. B **251**, 343 (1990).

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

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

[11] M. O. Olea-Romacho, Phys. Rev. D **109**, 015023 (2024).

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

[13] A. G. Tevzadze, L. Kisslinger, A. Brandenburg, and T. Kahniashvili, Astrophys. J. **759**, 54 (2012).
[14] T. Vachaspati, Phys. Rev. Lett. **87**, 251302 (2001).

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

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

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

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

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

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

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

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

[23] V. A. Acciari et al. (MAGIC Collaboration), Astron. Astrophys. **670**, A145 (2023).

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

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

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

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

[28] R. Roshan and G. White, arXiv:2401.04388.

[29] I. Baldes and M. O. Olea-Romacho, J. High Energy Phys. **01** (2024) 133.

[30] Y. Gouttenoire and T. Volansky, arXiv:2305.04942.

[31] B. J. Carr, Astrophys. J. **201**, 1 (1975).

[32] M. Shibata and M. Sasaki, Phys. Rev. D **60**, 084002 (1999).

[33] I. Musco, J. C. Miller, and L. Rezzolla, Classical Quantum Gravity **22**, 1405 (2005).

[34] T. Harada, C.-M. Yoo, and K. Kohri, Phys. Rev. D **88**, 084051 (2013); **89**, 029903(E) (2014).

[35] I. Musco, Phys. Rev. D **100**, 123524 (2019).

[36] I. Musco, V. De Luca, G. Franciolini, and A. Riotto, Phys. Rev. D **103**, 063538 (2021).

[37] C. Germani and I. Musco, Phys. Rev. Lett. **122**, 141302 (2019).

[38] A. Escrivà, C. Germani, and R. K. Sheth, Phys. Rev. D **101**, 044022 (2020).

[39] I. D. Stamou, Phys. Rev. D **108**, 063515 (2023).

[40] B. Carr and F. Kuhnel, SciPost Phys. Lect. Notes **48**, 1 (2022).

[41] B. J. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama, Phys. Rev. D **81**, 104019 (2010).

[42] M. Boudaud and M. Cirelli, Phys. Rev. Lett. **122**, 041104 (2019).

[43] H. Niikura et al., Nat. Astron. **3**, 524 (2019).

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

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

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

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

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

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

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

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

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

[53] S. Balaji, Phys. Lett. B **845**, 138157 (2023).

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

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

[56] M. Lewicki and V. Vaskonen, Eur. Phys. J. C **81**, 437 (2021); **81**, 1077(E) (2021).

[57] B. Allen and J. D. Romano, Phys. Rev. D **59**, 102001 (1999).

[58] H. Kudoh, A. Taruya, T. Hiramatsu, and Y. Himemoto, Phys. Rev. D **73**, 064006 (2006).

[59] E. Thrane and J. D. Romano, Phys. Rev. D **88**, 124032 (2013).

[60] C. Caprini, D. G. Figueroa, R. Flauger, G. Nardini, M. Peloso, M. Pieroni, A. Ricciardone, and G. Tasinato, J. Cosmol. Astropart. Phys. **11** (2019) 017.

[61] D. Brzeminski, A. Hook, and G. Marques-Tavares, J. High Energy Phys. **11** (2022) 061.

[62] M. Maggiore et al., J. Cosmol. Astropart. Phys. 03 (2020) 050.

[63] S. Balaji, G. Domènech, and G. Franciolini, J. Cosmol. Astropart. Phys. **10** (2023) 041.

[64] E. Barausse et al., Gen. Relativ. Gravit. **52**, 81 (2020).

[65] K. Yagi and N. Seto, Phys. Rev. D **83**, 044011 (2011); **95**, 109901(E) (2017).

[66] S. Kawamura et al., Prog. Theor. Exp. Phys. **2021**, 05A105 (2021).

[67] S. Balaji, G. Domènech, and J. Silk, J. Cosmol. Astropart. Phys. **09** (2022) 016.

[68] K. Schmitz, J. High Energy Phys. **01** (2021) 097.