# Magnetogenesis with gravitational waves and primordial black hole dark
matter

Shyam Balaji<EMAIL_ADDRESS>Physics Department, King’s College
London, Strand, London, WC2R 2LS, United Kingdom  Malcolm Fairbairn
<EMAIL_ADDRESS>Physics Department, King’s College London, Strand,
London, WC2R 2LS, United Kingdom  María Olalla Olea-Romacho
<EMAIL_ADDRESS>Laboratoire de physique de l’École normale
supérieure, ENS, Université PSL, CNRS, Sorbonne Université, Université Paris
Cité, F-75005 Paris, France

###### Abstract

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
$\lesssim 20$ pG across coherence length scales of $\simeq 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 multi-messenger
magnetic fields and GW signals visible in next generation experiments.

##  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 comprised of
primordial black holes (PBHs) Zel’dovich and Novikov (1967); Hawking (1971);
Carr and Hawking (1974); Carr  _et al._ (2021) if their masses are in the
asteroid-mass range Balaji  _et al._ (2022a); Qin  _et al._ (2023). 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 Kodama  _et
al._ (1982); Hsu (1990). 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 Vachaspati (1991); Ellis  _et al._
(2019); Olea-Romacho (2024) or the QCD phase transition Sigl  _et al._ (1997);
Tevzadze  _et al._ (2012). Notably, the idea of magnetogenesis during a first-
order electroweak phase transition was first proposed in Ref. Vachaspati
(1991). In this scenario, magnetic fields are created by EW sphaleron decays
Vachaspati (2001); Copi  _et al._ (2008). 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 Witten (1984); Hogan (1986); Kamionkowski
_et al._ (1994); Brandenburg  _et al._ (1996); Christensson  _et al._ (2001);
Kahniashvili  _et al._ (2010); Brandenburg  _et al._ (2017).

Coherent intergalactic magnetic fields (IGMFs) are indirectly suggested by
blazar observations Acciari  _et al._ (2023); Aharonian  _et al._ (2023);
Neronov and Vovk (2010). The origin of these IGMFs is a long-standing mystery,
with two main possibilities: astrophysical or cosmological. An example of
astrophysical sources is weak initial magnetic fields from local effects in
astrophysical objects (e.g. the Biermann battery mechanism Biermann and
Schlüter (1951)) that are later enhanced by dynamo effects Alves Batista and
Saveliev (2021), which can create such 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. Strong FOPTs also produce a stochastic background of
GWs Witten (1984); Hogan (1986); Roshan and White (2024). By enforcing the
production of PBH DM in the transition, we also compute the corresponding
gravitational wave (GW) signal from the bubble collisions, leading to
observable signals in upcoming GW interferometers well above stochastic
astrophysical foregrounds Baldes and Olea-Romacho (2024). In this way, we
highlight the relation between mulitmessenger cosmological signals of shared
origin in two disparate channels, magnetic fields and GWs, whilst
simultaneously providing an explanation for DM. We also briefly remark on a
possible model realisation 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

###  II.1 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_{\rm 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_{\rm nuc}$. Causal regions
(labelled by $i$) that nucleate at a time $t_{\rm nuc}^{i}$ later than
$\tau_{\rm 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}^{\rm
late}(t;t_{\textrm{{nuc}}}^{i})-\rho_{r}^{\rm bkg}(t)}{\rho_{r}^{\rm
bkg}(t)}>\delta_{c},$ (1)

where $\rho_{r}^{\rm late}(t;t_{\textrm{{nuc}}}^{i})$ is the radiation density
in a region where nucleation is delayed to $t_{\textrm{{nuc}}}^{i}$ and
$\rho_{r}^{\rm 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_{\rm bub}}\frac{d\Gamma_{\rm
bub}}{dt},$ (2)

where $\Gamma_{\rm bub}$ is the bubble nucleation rate per unit volume and $H$
is the Hubble rate during the phase transition.111In the following, we will
use a first-order Taylor expansion to approximate the bubble nucleation rate
near $\tau_{\rm nuc}$. In Ref. Baldes and Olea-Romacho (2024), it was shown
that the approximate form yields accurate predictions around the nucleation
time. The Hubble rate in the false vacuum reads,

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

where $\rho_{\rm 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_{\rm reh}=\left(\frac{30\rho_{\rm
vac}}{g_{r}(T_{\textrm{reh}})\pi^{2}}\right)^{1/4},$ (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 Gouttenoire and
Volansky (2023)

$M_{\rm
PBH}\left(T_{\textrm{reh}}\right)=M_{\odot}\left(\frac{20}{g_{r}(T_{\textrm{reh}})}\right)^{1/2}\left(\frac{140\,\textrm{MeV}}{T_{\textrm{reh}}}\right)^{2},$
(5)

whereas the PBH abundance normalised to the DM relic abundance reads
Gouttenoire and Volansky (2023)

$\displaystyle
f_{\textrm{PBH}}\left(\beta/H,T_{\textrm{reh}}\right)=\frac{e^{-a_{1}(\beta/H)_{\ast}^{a_{2}}\left(1+\delta_{c}\right)^{a_{3}(\beta/H)_{\ast}}}}{2.2\times
10^{-8}}\frac{T_{\textrm{reh}}}{0.14}.$ (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

$\displaystyle(\beta/H)_{*}=\frac{a_{2}W_{0}\left[\frac{a_{3}}{a_{2}}\left(\frac{19.5983+\log
T_{\textrm{reh}}}{a_{1}}\right)^{\frac{1}{a_{2}}}\log(1+\delta_{c})\right]}{a_{3}\log(1+\delta_{c})},$
(7)

where $W_{0}$ is the zeroth order Lambert $W$ function. Calculations, based on
full general relativistic simulations, often suggest a range of 0.4 to 0.66
for $\delta_{c}$ in the context of overdensities from inflation re-entering
the Hubble horizon Carr (1975); Shibata and Sasaki (1999); Musco  _et al._
(2005); Harada  _et al._ (2013); Musco (2019); Musco  _et al._ (2021); Germani
and Musco (2019); Escrivà  _et al._ (2020); Stamou (2023). Deviations from
these assumptions are expected in the phase transition scenario. However, in
this study, we align with previous approaches Gouttenoire and Volansky (2023);
Baldes and Olea-Romacho (2024), and choose $\delta_{c}=0.45$ to enable
comparison. The precise value of $\delta_{c}$ impacts $f_{\rm PBH}=1$, yet the
strong magnetic field generated remain insensitive to this choice.

####  II.1.1 Experimental constraints for PBHs as the dark matter

Observationally, $f_{\rm PBH}=1$ is allowed for nearly monochromatic
distributions of PBHs with masses falling within what is known as the
asteroid-mass window Carr and Kuhnel (2022)

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

Here the lower bound is set by the particle flux from Hawking evaporation,
which would affect the CMB Carr  _et al._ (2010) or be detected by the Voyager
Boudaud and Cirelli (2019) probe. The upper bound is set by microlensing
measurements Niikura  _et al._ (2019). According to Eq. (5), the range of
reheating temperatures that gives rise to PBHs with masses within the
asteroid-mass window is

$10\mathinner{\mathrm{TeV}}\lesssim T_{\rm reh}\lesssim
10^{4}\mathinner{\mathrm{TeV}}.$ (9)

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

###  II.2 Primordial magnetic fields

####  II.2.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 Banerjee and Jedamzik (2004). 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.

The MHD turbulence for maximally helical magnetic fields decays as a power law
in conformal time, $\eta$, with the following relations for the magnetic
energy and correlation length during the radiation-dominated epoch Biskamp and
Müller (1999)

$B\sim\eta^{-1/3}\,\,\mathrm{and}\,\,\lambda\sim\eta^{2/3}.$ (10)

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 Brandenburg  _et
al._ (2019).222However, this has been challenged recently in Ref. Armua  _et
al._ (2023), where a weaker inverse cascade of magnetic energy for non-helical
fields was found, compared to previous studies in the literature. In this
scenario, the magnetic energy and correlation length, $\lambda$, scale as

$B\sim\eta^{-1/2}\,\,\,\mathrm{and}\,\,\,\lambda\sim\eta^{1/2}\,,$ (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}\,\,\,\,\mathrm{and}\,\,\,\,p_{b}=\frac{2}{b+3}(b+1)\,,$
(12)

for the power-laws

$B\sim\eta^{-p_{b}/2}\,\,\,\,\mathrm{and}\,\,\,\,\lambda\sim\eta^{q_{b}}\,,$
(13)

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

The magnetic field energy density at percolation can be estimated by Ellis
_et al._ (2019); Roper Pol  _et al._ (2023)

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

Here, $\rho_{*}=\frac{3M_{\rm
Pl}^{2}}{8\pi}H_{*}^{2}=\frac{g_{r}(T_{\textrm{reh}})\pi^{2}}{30}T_{\rm
reh}^{4}\approx\rho_{\rm vac}$ is the total energy density at the percolation
temperature. Here $\kappa_{\rm 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\%$ Kahniashvili  _et al._ (2011); Durrer and Neronov
(2013); Brandenburg  _et al._ (2017). For supercooled phase transitions, we
can safely assume $\alpha\rightarrow\infty$ and $\kappa_{\rm col}=1$, leading
to the approximation in Eq. (14).

The magnetic field spectrum today can be computed as Ellis  _et al._ (2019)

$\displaystyle B_{0}(\lambda)\equiv B(\lambda,t_{0})=$ (15)
$\displaystyle\left(\frac{a_{\rm reh}}{a_{\rm
rec}}\right)^{p_{b}/2}\left(\frac{a_{\rm
reh}}{a_{0}}\right)^{2}\sqrt{\frac{17}{10}\,\rho_{B,*}}\,\begin{cases}\left(\frac{\lambda}{\lambda_{0}}\right)^{-5/2}&\mbox{
for }~{}~{}\lambda\geq\lambda_{0}\\\
\left(\frac{\lambda}{\lambda_{0}}\right)^{1/3}&\mbox{ for
}~{}~{}\lambda<\lambda_{0}\,,\end{cases}$

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 Ellis  _et al._ (2019)

$\lambda_{0}\equiv\lambda_{B}(t_{0})=\left(\frac{a_{\rm rec}}{a_{\rm
reh}}\right)^{q_{b}}\left(\frac{a_{0}}{a_{\rm reh}}\right)\lambda_{*},\,$ (16)

where the initial correlation length $\lambda_{*}$ is given by the bubble size
at percolation $R_{*}$ Caprini  _et al._ (2020)

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

where, for the case of supercooled phase transitions, a bubble wall velocity
of $v_{w}=1$ was assumed. The peak value of the magnetic field spectrum is
denoted as $B_{0}=B_{0}(\lambda_{0})$. The redshift factors are computed as

$\displaystyle\frac{a_{\rm reh}}{a_{0}}=8\times
10^{-14}\left(\frac{100}{g_{r}(T_{\rm
reh})}\right)^{1/3}\left(\frac{1\mathinner{\mathrm{GeV}}}{T_{\rm
reh}}\right)\,,$ (18) $\displaystyle\frac{a_{\rm reh}}{a_{\rm rec}}=8\times
10^{-11}\left(\frac{100}{g_{r}(T_{\rm
reh})}\right)^{1/3}\left(\frac{1\mathinner{\mathrm{GeV}}}{T_{\rm
reh}}\right)\,,$ (19)

where $g_{r}(T_{\rm reh})$ corresponds to the total number of relativistic
degrees of freedom in entropy at the reheating temperature $T_{\rm reh}$. For
this work, we will assume $g_{r}(T_{\rm 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.

####  II.2.2 Experimental constraints from blazar emissions

Blazars are active galactic nuclei jets that are roughly pointed towards us,
and they emit TeV $\gamma$-rays. These $\gamma$-rays can collide with photons
from the extragalactic background, creating $e^{+}e^{-}$ pairs Balaji (2023).
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 Vachaspati (2017, 2021). 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 towards 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}\,\mathrm{G}$ for
correlation lengths $\lambda\geq 0.2\,\mathrm{Mpc}$ Acciari  _et al._ (2023).
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}\,\mathrm{G}$ for
coherence lengths of $\lambda\geq 1\,\mathrm{Mpc}$, assuming a blazar duty
cycle of $t_{d}=10\,\mathrm{yrs}$ Aharonian  _et al._ (2023).

###  II.3 Gravitational waves

We use the results from hybrid simulations Lewicki and Vaskonen (2021) 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_{\rm gw}(f)=5.10\times
10^{-9}\left(\frac{100}{g_{r}(T_{\textrm{reh}})}\right)^{1/3}\left(\frac{10}{(\beta/H)_{*}}\right)^{2}S_{\rm
hy}(f),$ (20)

where the shape is

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

and peak frequency

$\tilde{f}_{\rm
hy}=22\;\mathrm{mHz}\,\left(\frac{g_{r}(T_{\textrm{reh}})}{100}\right)^{1/6}\left(\frac{(\beta/H)_{*}}{10}\right)\,\left(\frac{T_{\rm
reh}}{10^{2}\;\mathrm{TeV}}\right).$ (22)

The signal-to-noise ratio is given by Allen and Romano (1999); Kudoh  _et al._
(2006); Thrane and Romano (2013); Caprini  _et al._ (2019); Brzeminski  _et
al._ (2022)

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

where $t_{\rm obs}$ is the observation time and $\Omega_{\rm sens}$ 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 DECIGO (Deci-hertz Interferometer
Gravitational wave Observatory). At higher frequency, we have a band that
would be observable with the ET Maggiore  _et al._ (2020), which would be
sensitive to the range $10$–$10^{3}$ Hz. In between Square Kilometre Array
Balaji  _et al._ (2023) and ET frequencies we expect LISA Barausse  _et al._
(2020) and DECIGO Yagi and Seto (2011); Kawamura  _et al._ (2021); Balaji  _et
al._ (2022b) 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. Schmitz (2021)
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
length. We can also evaluate the GW signal generated in this transition.

Figure 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_{\textrm{reh}}\in[10^{4},10^{7}]$ GeV, $b=0$ and $f_{\textrm{pbh}}=1$. The
solid red line corresponds to the lower bound on the magnetic field set by
Ref. Acciari  _et al._ (2023), whereas the dashed red line corresponds to that
set by Ref. Aharonian  _et al._ (2023) for a duty cycle of
$t_{d}=10\,\mathrm{yrs}$.

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}]\mathinner{\mathrm{GeV}}$. The coherence length scales for
combinations of reheating temperatures and $\beta/H$ values that produce
$f_{\textrm{pbh}}=1$ lie in a relatively narrow window of $6\times 10^{-4}$ to
$7\times 10^{-3}\,\rm Mpc$. Across this range, the magnetic field amplitude
varies from around 1-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 Ellis  _et al._ (2019) (see
Sec. III.1). 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 $\simeq 3.4\times 10^{-3}$ Mpc.

We note that in the case of a non-helical magnetic field ($b=1$), the peak
magnetic field gets suppressed by more than two 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 non-helical magnetic fields
($b=1$), as the magnetic field amplitudes in this scenario are relatively
suppressed.

Figure 2: Magnetic field peak value in the present universe $B_{0}$ vs the
primordial black hole mass (top panel). Coherence length $\lambda_{0}$ vs the
primordial black hole mass (bottom panel) is also shown. In both cases, we
consider a reheating temperature of $T_{\textrm{reh}}\in[10^{4},10^{7}]$ GeV,
$b=0$ and $f_{\textrm{pbh}}=1$.

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\,\rm Mpc$.

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.3 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-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
multi-messenger signal yielding detectable GWs whilst producing a correlated
cosmological magnetic field and all the requisite DM content of the universe
in the form of PBHs.

Figure 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_{\textrm{reh}}\in[10^{4},10^{7}]$
GeV, $b=0$ and $f_{\textrm{pbh}}=1$.

###  III.1 A model example: the classically conformal $B-L$ model

To provide a concrete realization, we consider electroweak symmetry breaking
in the classically scale-invariant $U(1)_{B-L}$ extension of the SM. The
radiative breaking of the $B-L$ symmetry induces a negative mass term for the
EW Higgs field via a portal coupling, thus eventually breaking EW symmetry. It
is during this process that magnetic fields can be generated. Subsequently,
plasma effects restore the EW symmetry and, finally, the universe cools down
post-reheating and the EW symmetry breaks in a similar fashion to the SM. For
details on the model and its finite temperature treatment in cosmological
phase transitions, see e.g. Ref. Baldes and Olea-Romacho (2024). Considering a
$Z^{\prime}$ with a mass of $m_{Z^{\prime}}=56.64\mathinner{\mathrm{TeV}}$ and
$B-L$ gauge coupling of $g_{B-L}=0.2832$: we get a nucleation temperature of
$T_{\rm nuc}=2.03\mathinner{\mathrm{GeV}}$, a reheating temperature of $T_{\rm
reh}=5.03\mathinner{\mathrm{TeV}}$, and the inverse duration of the transition
$\beta/H=6.81$. In this context, the predicted PBH abundance and PBH mass are

$\displaystyle f_{\rm PBH}=1,$ (24) $\displaystyle M_{\rm PBH}=5\times
10^{-10}M_{\odot}.$ (25)

At the same time, the peak strength of the generated magnetic fields is
$B_{0}=1.9\times 10^{-11}\,\textrm{G}$ at a peak coherence length of
$\lambda_{0}=0.008\,\textrm{Mpc}$. In this case we have estimated $\kappa_{h}$
as $\kappa_{h}=\Delta V_{h}/\Delta V_{\varphi}$, where $\Delta
V_{\varphi}=V(v_{\varphi},0)-V(0,0)$ and
$V_{\varphi}=V(v_{\varphi},0)-V(v_{\varphi},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}\,\rm{G}$ at a coherence length scale
$\lambda_{0}=0.008\,\textrm{Mpc}$.

##  IV Conclusion

We investigate the intriguing connection between strongly supercooled first-
order phase transitions (FOPTs), primordial black hole (PBH) dark matter (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 to 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 Deci-
hertz Interferometer Gravitational wave Observatory, Laser Interferometer
Space Antenna, Einstein Telescope and Cosmic Explorer. This presents the
exciting prospect of a multi-messenger 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

MOOR thanks I. Baldes, K. Kamada and T. Vachaspati while MF thanks V. Vaskonen
for useful discussions. SB and MF are supported by the STFC under grant
ST/X000753/1. MOOR is supported by the European Union’s Horizon 2020 research
and innovation programme under grant agreement No 101002846, ERC CoG
“CosmoChart”.

## References

  * Zel’dovich and Novikov (1967) Y. B. Zel’dovich and I. D. Novikov, Soviet Astron. AJ (Engl. Transl. ), 10, 602 (1967).
  * Hawking (1971) S. Hawking, Mon. Not. Roy. Astron. Soc. 152, 75 (1971).
  * Carr and Hawking (1974) B. J. Carr and S. W. Hawking, Mon. Not. Roy. Astron. Soc. 168, 399 (1974).
  * Carr  _et al._ (2021) B. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama, Rept. Prog. Phys. 84, 116902 (2021), arXiv:2002.12778 [astro-ph.CO] .
  * Balaji  _et al._ (2022a) S. Balaji, J. Silk, and Y.-P. Wu, JCAP 06, 008 (2022a), arXiv:2202.00700 [astro-ph.CO] .
  * Qin  _et al._ (2023) W. Qin, S. R. Geller, S. Balaji, E. McDonough, and D. I. Kaiser, Phys. Rev. D 108, 043508 (2023), arXiv:2303.02168 [astro-ph.CO] .
  * Kodama  _et al._ (1982) H. Kodama, M. Sasaki, and K. Sato, Prog. Theor. Phys. 68, 1979 (1982).
  * Hsu (1990) S. D. H. Hsu, Phys. Lett. B 251, 343 (1990).
  * Vachaspati (1991) T. Vachaspati, Phys. Lett. B 265, 258 (1991).
  * Ellis  _et al._ (2019) J. Ellis, M. Fairbairn, M. Lewicki, V. Vaskonen, and A. Wickens, JCAP 09, 019 (2019), arXiv:1907.04315 [astro-ph.CO] .
  * Olea-Romacho (2024) M. O. Olea-Romacho, Phys. Rev. D 109, 015023 (2024), arXiv:2310.19948 [hep-ph] .
  * Sigl  _et al._ (1997) G. Sigl, A. V. Olinto, and K. Jedamzik, Phys. Rev. D 55, 4582 (1997), arXiv:astro-ph/9610201 .
  * Tevzadze  _et al._ (2012) A. G. Tevzadze, L. Kisslinger, A. Brandenburg, and T. Kahniashvili, Astrophys. J. 759, 54 (2012), arXiv:1207.0751 [astro-ph.CO] .
  * Vachaspati (2001) T. Vachaspati, Phys. Rev. Lett. 87, 251302 (2001), arXiv:astro-ph/0101261 .
  * Copi  _et al._ (2008) C. J. Copi, F. Ferrer, T. Vachaspati, and A. Achucarro, Phys. Rev. Lett. 101, 171302 (2008), arXiv:0801.3653 [astro-ph] .
  * Witten (1984) E. Witten, Phys. Rev. D 30, 272 (1984).
  * Hogan (1986) C. J. Hogan, Mon. Not. Roy. Astron. Soc. 218, 629 (1986).
  * Kamionkowski  _et al._ (1994) M. Kamionkowski, A. Kosowsky, and M. S. Turner, Phys. Rev. D 49, 2837 (1994), arXiv:astro-ph/9310044 .
  * Brandenburg  _et al._ (1996) A. Brandenburg, K. Enqvist, and P. Olesen, Phys. Rev. D 54, 1291 (1996), arXiv:astro-ph/9602031 .
  * Christensson  _et al._ (2001) M. Christensson, M. Hindmarsh, and A. Brandenburg, Phys. Rev. E 64, 056405 (2001), arXiv:astro-ph/0011321 .
  * Kahniashvili  _et al._ (2010) T. Kahniashvili, A. Brandenburg, A. G. Tevzadze, and B. Ratra, Phys. Rev. D 81, 123002 (2010), arXiv:1004.3084 [astro-ph.CO] .
  * Brandenburg  _et al._ (2017) A. Brandenburg, T. Kahniashvili, S. Mandal, A. Roper Pol, A. G. Tevzadze, and T. Vachaspati, Phys. Rev. D 96, 123528 (2017), arXiv:1711.03804 [astro-ph.CO] .
  * Acciari  _et al._ (2023) V. A. Acciari _et al._ (MAGIC), Astron. Astrophys. 670, A145 (2023), arXiv:2210.03321 [astro-ph.HE] .
  * Aharonian  _et al._ (2023) F. Aharonian _et al._ (H.E.S.S., Fermi-LAT), Astrophys. J. Lett. 950, L16 (2023), arXiv:2306.05132 [astro-ph.HE] .
  * Neronov and Vovk (2010) A. Neronov and I. Vovk, Science 328, 73 (2010), arXiv:1006.3504 [astro-ph.HE] .
  * Biermann and Schlüter (1951) L. Biermann and A. Schlüter, Phys. Rev. 82, 863 (1951).
  * Alves Batista and Saveliev (2021) R. Alves Batista and A. Saveliev, Universe 7, 223 (2021), arXiv:2105.12020 [astro-ph.HE] .
  * Roshan and White (2024) R. Roshan and G. White, (2024), arXiv:2401.04388 [hep-ph] .
  * Baldes and Olea-Romacho (2024) I. Baldes and M. O. Olea-Romacho, JHEP 01, 133 (2024), arXiv:2307.11639 [hep-ph] .
  * Gouttenoire and Volansky (2023) Y. Gouttenoire and T. Volansky, (2023), arXiv:2305.04942 [hep-ph] .
  * Carr (1975) B. J. Carr, Astrophys. J. 201, 1 (1975).
  * Shibata and Sasaki (1999) M. Shibata and M. Sasaki, Phys. Rev. D 60, 084002 (1999), arXiv:gr-qc/9905064 .
  * Musco  _et al._ (2005) I. Musco, J. C. Miller, and L. Rezzolla, Class. Quant. Grav. 22, 1405 (2005), arXiv:gr-qc/0412063 .
  * Harada  _et al._ (2013) T. Harada, C.-M. Yoo, and K. Kohri, Phys. Rev. D 88, 084051 (2013), [Erratum: Phys.Rev.D 89, 029903 (2014)], arXiv:1309.4201 [astro-ph.CO] .
  * Musco (2019) I. Musco, Phys. Rev. D 100, 123524 (2019), arXiv:1809.02127 [gr-qc] .
  * Musco  _et al._ (2021) I. Musco, V. De Luca, G. Franciolini, and A. Riotto, Phys. Rev. D 103, 063538 (2021), arXiv:2011.03014 [astro-ph.CO] .
  * Germani and Musco (2019) C. Germani and I. Musco, Phys. Rev. Lett. 122, 141302 (2019), arXiv:1805.04087 [astro-ph.CO] .
  * Escrivà  _et al._ (2020) A. Escrivà, C. Germani, and R. K. Sheth, Phys. Rev. D 101, 044022 (2020), arXiv:1907.13311 [gr-qc] .
  * Stamou (2023) I. D. Stamou, Phys. Rev. D 108, 063515 (2023), arXiv:2306.02758 [astro-ph.CO] .
  * Carr and Kuhnel (2022) B. Carr and F. Kuhnel, SciPost Phys. Lect. Notes 48, 1 (2022), arXiv:2110.02821 [astro-ph.CO] .
  * Carr  _et al._ (2010) B. J. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama, Phys. Rev. D 81, 104019 (2010), arXiv:0912.5297 [astro-ph.CO] .
  * Boudaud and Cirelli (2019) M. Boudaud and M. Cirelli, Phys. Rev. Lett. 122, 041104 (2019), arXiv:1807.03075 [astro-ph.HE] .
  * Niikura  _et al._ (2019) H. Niikura _et al._ , Nature Astron. 3, 524 (2019), arXiv:1701.02151 [astro-ph.CO] .
  * Banerjee and Jedamzik (2004) R. Banerjee and K. Jedamzik, Phys. Rev. D 70, 123003 (2004), arXiv:astro-ph/0410032 .
  * Biskamp and Müller (1999) D. Biskamp and W.-C. Müller, Phys. Rev. Lett. 83, 2195 (1999).
  * Brandenburg  _et al._ (2019) A. Brandenburg, T. Kahniashvili, S. Mandal, A. Roper Pol, A. G. Tevzadze, and T. Vachaspati, Phys. Rev. Fluids. 4, 024608 (2019), arXiv:1710.01628 [physics.flu-dyn] .
  * Armua  _et al._ (2023) A. Armua, A. Berera, and J. C. Figueroa, Phys. Rev. E 107, 055206 (2023), arXiv:2212.02418 [physics.plasm-ph] .
  * Roper Pol  _et al._ (2023) A. Roper Pol, A. Neronov, C. Caprini, T. Boyer, and D. Semikoz, (2023), arXiv:2307.10744 [astro-ph.CO] .
  * Kahniashvili  _et al._ (2011) T. Kahniashvili, A. G. Tevzadze, and B. Ratra, Astrophys. J. 726, 78 (2011), arXiv:0907.0197 [astro-ph.CO] .
  * Durrer and Neronov (2013) R. Durrer and A. Neronov, Astron. Astrophys. Rev. 21, 62 (2013), arXiv:1303.7121 [astro-ph.CO] .
  * Caprini  _et al._ (2020) C. Caprini _et al._ , JCAP 03, 024 (2020), arXiv:1910.13125 [astro-ph.CO] .
  * Balaji (2023) S. Balaji, Phys. Lett. B 845, 138157 (2023), arXiv:2211.03807 [astro-ph.HE] .
  * Vachaspati (2017) T. Vachaspati, Phys. Rev. D 95, 063505 (2017), arXiv:1606.06186 [astro-ph.CO] .
  * Vachaspati (2021) T. Vachaspati, Rept. Prog. Phys. 84, 074901 (2021), arXiv:2010.10525 [astro-ph.CO] .
  * Lewicki and Vaskonen (2021) M. Lewicki and V. Vaskonen, Eur. Phys. J. C 81, 437 (2021), [Erratum: Eur.Phys.J.C 81, 1077 (2021)], arXiv:2012.07826 [astro-ph.CO] .
  * Allen and Romano (1999) B. Allen and J. D. Romano, Phys. Rev. D 59, 102001 (1999), arXiv:gr-qc/9710117 .
  * Kudoh  _et al._ (2006) H. Kudoh, A. Taruya, T. Hiramatsu, and Y. Himemoto, Phys. Rev. D 73, 064006 (2006), arXiv:gr-qc/0511145 .
  * Thrane and Romano (2013) E. Thrane and J. D. Romano, Phys. Rev. D 88, 124032 (2013), arXiv:1310.5300 [astro-ph.IM] .
  * Caprini  _et al._ (2019) C. Caprini, D. G. Figueroa, R. Flauger, G. Nardini, M. Peloso, M. Pieroni, A. Ricciardone, and G. Tasinato, JCAP 11, 017 (2019), arXiv:1906.09244 [astro-ph.CO] .
  * Brzeminski  _et al._ (2022) D. Brzeminski, A. Hook, and G. Marques-Tavares, JHEP 11, 061 (2022), arXiv:2203.13842 [hep-ph] .
  * Maggiore  _et al._ (2020) M. Maggiore _et al._ , JCAP 03, 050 (2020), arXiv:1912.02622 [astro-ph.CO] .
  * Balaji  _et al._ (2023) S. Balaji, G. Domènech, and G. Franciolini, JCAP 10, 041 (2023), arXiv:2307.08552 [gr-qc] .
  * Barausse  _et al._ (2020) E. Barausse _et al._ , Gen. Rel. Grav. 52, 81 (2020), arXiv:2001.09793 [gr-qc] .
  * Yagi and Seto (2011) K. Yagi and N. Seto, Phys. Rev. D 83, 044011 (2011), [Erratum: Phys.Rev.D 95, 109901 (2017)], arXiv:1101.3940 [astro-ph.CO] .
  * Kawamura  _et al._ (2021) S. Kawamura _et al._ , PTEP 2021, 05A105 (2021), arXiv:2006.13545 [gr-qc] .
  * Balaji  _et al._ (2022b) S. Balaji, G. Domenech, and J. Silk, JCAP 09, 016 (2022b), arXiv:2205.01696 [astro-ph.CO] .
  * Schmitz (2021) K. Schmitz, JHEP 01, 097 (2021), arXiv:2002.04615 [hep-ph] .