# Constraining small-scale primordial magnetic fields from the abundance of
primordial black holes

Ashu Kushwaha<EMAIL_ADDRESS>Department of Physics, Indian Institute of
Science C. V. Raman Road, Bangalore 560012, India  Department of Physics,
Indian Institute of Technology Bombay, Mumbai 400076, India  Teruaki Suyama
Department of Physics, Tokyo Institute of Technology, 2-12-1 Ookayama, Meguro-
ku, Tokyo 152-8551, Japan

###### Abstract

The presence of magnetic fields in the early universe affects the cosmological
processes, leading to the distinct signature, which allows constraining their
properties and the genesis mechanisms. In this study, we revisit the method to
constrain the amplitude of the magnetic fields on small scales in the
radiation-dominated era from the abundance of primordial black holes.
Constraints in the previous work were based on the fact that the density
perturbations sourced by stronger magnetic fields become large enough to
gravitationally collapse to form PBHs. However, we demonstrate that this
picture is incomplete because magnetic fields also increase the threshold
value of the density contrast required for PBH formation. The increase in
threshold density contrast is more pronounced on smaller scales, and in
extreme cases, it might even prevent PBH production despite the presence of
significant magnetic field. Taking into account the relevant physical effects
on the magnetized overdense region, we establish an upper-limit on the
amplitude of comoving magnetic fields, approximately $0.12-0.14{\rm\mu G}$ at
a scale of $10^{17}{\rm Mpc}^{-1}$. Additionally, we compare our constraints
with various small-scale probes.

##  I Introduction

Observational probes have confirmed the presence of magnetic fields in the
universe at all length scales. They have been observed in the galaxies and
galaxy clusters with the typical strength of order micro-Gauss with a
coherence length over kpc to Mpc Kronberg (1994); Beck  _et al._ (1996);
Widrow (2002); Vallée (2004); Gaensler  _et al._ (2004); Giovannini (2004);
Barrow  _et al._ (2007); Kandus  _et al._ (2011); Durrer and Neronov (2013);
Subramanian (2016). The origin of these magnetic fields is still unknown and
they could have cosmological or astrophysical origin. The FERMI and High
Energy Stereoscopic System (HESS) measurements of $\gamma$-rays from blazars
provide a lower bound of the order of $10^{-16}\rm{G}$ in intergalactic voids
Neronov and Vovk (2010). In the literature, various proposals argue that this
lower bound suggests that these magnetic fields could be of primordial origin,
which are generated in the early universe and later amplified due to dynamo
mechanism Turner and Widrow (1988); Ratra (1992); Vachaspati (1991); Dolgov
(1993); Martin and Yokoyama (2008); Mukohyama (2016); Kushwaha and
Shankaranarayanan (2020); Bamba  _et al._ (2021); Giovannini (2021); Tripathy
_et al._ (2021); Nandi (2021). In most of these scenarios, either the magnetic
field is too weak to be amplified via dynamo action or the coherence length is
too small. In contrast to these proposals, in a recent study, it is shown that
the primordial magnetic fields (PMFs) on large scales are ruled out from
baryon isocurvature constraints Kamada  _et al._ (2021).

The magnetic fields in the early universe affect the cosmological processes,
leading to distinct signatures that constrain the properties of the PMFs. For
example, magnetic field energy density at the Big Bang Nucleosynthesis (BBN)
epoch might increase ${\rm He}$ production beyond the observational bounds,
which allows to constrain its comoving amplitude $<1.5{\rm\mu G}$ independent
of the coherence scale Kawasaki and Kusakabe (2012). Observations of the
temperature and polarization anisotropies of CMB provide one of the most
promising ways to detect or constrain the properties of PMFs. A homogeneous
PMFs due to its preferred spatial direction would lead to anisotropic
expansion around that direction, hence leading to a quadrupole anisotropy in
the CMB, which puts an upper-limit on the amplitude of order nano-Gauss (on
Mpc scales) at the present epoch Durrer and Neronov (2013); Subramanian
(2016). Furthermore, before recombination, PMFs induce small-scale baryonic
density fluctuations, leading to an inhomogeneous recombination process that
affects the CMB anisotropies. Using magnetohydrodynamic (MHD) effects, authors
in Ref. Jedamzik and Saveliev (2019) put the most stringent upper-limit of
$\lesssim 47{\rm pG}$ (present day) for scale-invariant PMFs power spectra.

While most of the above-mentioned studies focus on explaining or constraining
the properties of large-scale magnetic fields, small-scale PMFs have not
received much attention. Given that the early universe was magnetized and the
most favorable scenarios to explain the origin of the PMFs should occur before
the recombination epoch, we can ask if there is a way to constrain the
properties of PMFs on small scales without being affected by sub-Horizon scale
processes such as MHD or dynamo action. By doing that, we can not only gain a
better understanding of the nature of _seed_ magnetic field, but also
constrain the early-time magnetogenesis scenarios. Small-scale magnetic fields
in the radiation-dominated (RD) era can lead to some interesting consequences;
for example, they can create matter-antimatter asymmetry Kamada and Long
(2016); Fujita and Kamada (2016); Kushwaha and Shankaranarayanan (2021), CMB
distortion Kunze and Komatsu (2014), magnetic reheating Saga  _et al._
(2018a), gravitational wave background (GWB) Saga  _et al._ (2018b) and might
also affect the primordial black hole (PBH) formation Saga  _et al._ (2020).
To constrain the upper-limit of the amplitude of small scale PMFs, we need to
investigate the observable signatures of PMFs on small-scale cosmological
processes, for example GWB and PBHs formation.

In this work, we present a systematic study to derive the upper-limit on the
amplitude of PMFs by taking into account the physical effects affecting the
formation of PBHs from magnetized overdense regions. This program can be
achieved by using the fact that the magnetic fields surpassing a critical
threshold source the primordial curvature perturbation to a degree that
results in the over-production of PBHs, exceeding the observational upper
limit on their abundance. Upper limits on PMFs based on this consideration has
been already given in the previous work Saga  _et al._ (2020). In the present
paper, we point out a new effect which was over-looked in the previous study:
change of the threshold value of the PBH formation by the pressure of the
magnetic fields. This effect acts as weakening the upper limits on the PMFs.
Yet, as it will turn out, our constraints are severer than those provided in
Ref. Saga  _et al._ (2020) because of another effect (i.e., temperature
dependence of the relativistic degrees of freedom) which is not included in
the previous study. To simplify the calculations and to make the direct
comparison with the previous work, we consider the same power spectrum for the
PMFs as the one considered in Saga  _et al._ (2020). We also compare the
constraints for both Gaussian and non-Gaussian type probability distribution
functions of the density contrast.

Before closing this section, we stress that our constraints are fairly generic
in the sense that they are independent of the generation mechanism.
Furthermore, our constraints will also provide a great insight into
understanding the magnetogenesis scenarios and might rule out certain class of
models.

This work is organized in the following way: in section II, we discuss the
evolution of curvature perturbations sourced by PMFs-induced anisotropic
stress and calculate the density contrast. Section III discusses the effects
of magnetic fields on the overdense region, then in section IV, we compute
abundance of PBHs and V discusses our constraints on the PMFs and comparison
with other observational probes. We work with the natural units where
$\hbar=c=k_{B}=1$ and the reduced Planck mass $M_{\rm Pl}=2.43\times
10^{18}{\rm GeV}$. For the magnetic field, we follow the electromagnetic CGS
units where $1\rm{G}=6.92\times 10^{-20}{\rm GeV^{2}}$.

##  II Anisotropic stress and density contrast sourced by PMFs

###  II.1 Curvature perturbation from PMFs

In this section, we discuss the evolution of primordial curvature perturbation
sourced by the anisotropic stress of the PMFs. We assume that the
inhomogeneous magnetic fields are generated before the radiation-dominated era
(for example, during inflation), with the energy density of the magnetic field
being the first-order perturbation to the Friedmann–Lemaître–Robertson–Walker
(FLRW) spacetime Mack  _et al._ (2002); Suyama and Yokoyama (2012). Since the
primordial plasma, which is produced by the decay of the inflaton after
inflation, is highly conducting, magnetic fields are frozen-in and decay as
$1/a^{2}(\eta)$. This allows to separate the spatial and temporal dependence
of the magnetic fields on sufficiently large-scales as
$\textbf{B}(\textbf{x},\eta)=\textbf{B}(\textbf{x})/a^{2}(\eta)$, where $\eta$
is the conformal time, and we set the scale factor (in cosmic time $t$) at
present epoch to be unity, $a(t_{0})=1$. The Fourier transform of the magnetic
field is defined as:

$\textbf{B}(\textbf{k})=\int
d^{3}x\textbf{B}(\textbf{x})e^{-i\textbf{k}\cdot\textbf{x}}.$ (1)

Assuming $\textbf{B}(\textbf{k})$ is a Gaussian random field, the statistical
properties of the field $\textbf{B}(\textbf{k})$ are determined by the two-
point correlation function given by Mack  _et al._ (2002); Shaw and Lewis
(2010); Durrer and Neronov (2013):

$\displaystyle\langle
B_{i}(\textbf{k})B_{j}(\textbf{k}^{\prime})\rangle=(2\pi)^{3}\delta^{3}(\textbf{k}-\textbf{k}^{\prime})\frac{P_{ij}(\hat{k})}{2}P_{B}(k)$
(2)

where $\hat{k}_{i}=k_{i}/k$ are unit wavenumber components,
$\delta^{3}(\textbf{k}-\textbf{k}^{\prime})$ is 3D Dirac-delta function.
$P_{ij}(\hat{k})=\delta_{ij}-\hat{k}_{i}\hat{k}_{j}$ is projection tensor that
comes from the divergence condition $\nabla\cdot\textbf{B}=0$ and $P_{B}(k)$
is the magnetic power spectrum.

Now, let us consider the spatial part of the energy-momentum tensor due to the
magnetic fields $T_{ij}$. Following Shaw and Lewis (2010), we express $T_{ij}$
as $T_{ij}=p_{\gamma}(g_{ij}\Delta_{B}+\Pi_{ij}^{B})$, where $\Delta_{B}$ is
the perturbations in energy density, $p_{\gamma}$ is the pressure of photons
and $\Pi_{ij}^{B}$, which sources the curvature perturbation through the mode-
mode coupling, is the anisotropic stress due to magnetic fields which is
traceless i.e., ${\Pi^{B}}^{i}_{i}=0$. The scalar part of the anisotropic
stress in momentum space is given by Shaw and Lewis (2010):

$\displaystyle\Pi_{B}(\mathbf{k},T)$
$\displaystyle=-\frac{3}{2}\left(\hat{k}^{i}\hat{k}^{j}-\frac{\delta^{ij}}{3}\right)\Pi^{ij}_{B}$
$\displaystyle=\left(\hat{k}^{i}\hat{k}^{j}-\frac{\delta^{ij}}{3}\right)\frac{9}{8\pi\rho_{\gamma}(T)a^{4}(T)}\int\frac{d^{3}\textbf{k}^{\prime}}{(2\pi)^{3}}B_{i}(\mathbf{k}^{\prime})B_{j}(\mathbf{k}-\textbf{k}^{\prime})$
(3)

and the spatial indices are raised and lowered by the Kronecker delta
$\delta_{ij}$ and $\rho_{\gamma}(T)$ is the photon energy density which varies
with temperature of the Universe. Using the conservation of entropy i.e.,
$g_{*s}(T)a^{3}(T)T^{3}=\text{constant}$, we obtain

$\displaystyle\rho_{\gamma}(T)=\rho_{\gamma,0}\left(\frac{a(T_{0})}{a(T)}\right)^{4}\left(\frac{g_{*s}(T_{0})}{g_{*s}(T)}\right)^{4/3}$
(4)

where $\rho_{\gamma,0}$ is the photon energy density at present epoch
($T_{0}$) and $g_{*s}$ is a number of relativistic degree of freedom. Thus,
using Eq. (4) in Eq. (II.1) gives the scalar part of the anisotropic stress as

$\displaystyle\Pi_{B}(\textbf{k},T)=\left(\hat{k}^{i}\hat{k}^{j}-\frac{\delta^{ij}}{3}\right)\frac{9}{8\pi\rho_{\gamma,0}}\left(\frac{g_{*s}(T)}{g_{*s}(T_{0})}\right)^{4/3}\int\frac{d^{3}\textbf{k}^{\prime}}{(2\pi)^{3}}B_{i}(\mathbf{k}^{\prime})B_{j}(\mathbf{k}-\textbf{k}^{\prime}).$
(5)

Note that the above expression differs from Ref. Saga  _et al._ (2020) wherein
the dependence on relativistic degrees of freedom is missing. The two-point
correlation function of the scalar part of the anisotropic stress can be
obtained as

$\displaystyle\langle\Pi_{B}(\textbf{k},T)\Pi_{B}^{*}(\textbf{p},T)\rangle$
$\displaystyle=\delta^{3}(\textbf{k}-\textbf{p})\left(\frac{9}{16\pi\rho_{\gamma,0}^{2}}\right)\left(\frac{g_{*s}(T)}{g_{*s}(T_{0})}\right)^{8/3}\int
dk^{\prime}{k^{\prime}}^{2}P_{B}(k^{\prime})\int_{-1}^{1}d\gamma
P_{B}(|\textbf{k}-\mathbf{k}^{\prime}|)$
$\displaystyle{}\qquad\qquad\times\left[\frac{k^{2}(1+\gamma^{2})+kk^{\prime}(2\gamma-6\gamma^{3})+{k^{\prime}}^{2}(5+9\gamma^{4}-12\gamma^{2})}{4(k^{2}+{k^{\prime}}^{2}-2\gamma
kk^{\prime})}\right]$ (6)

where $\gamma=\hat{k}\cdot\hat{k}^{\prime}$.

The anisotropic stress due to primordial magnetic fields acts as a source for
the curvature perturbation. In particular, it induces two kinds of
contributions to the primordial curvature perturbations Shaw and Lewis (2010):
(i) _passive modes_ , which are induced before neutrino decoupling on both
super-horizon scales and sub-horizon scales, and (ii) _compensated modes_ ,
which are induced only on sub-horizon scales. In this work, we will only focus
on passive modes because only super-horizon modes are relevant for the
formation of PBHs in the radiation-dominated era. Note that after neutrino
decoupling time $\eta_{\nu}$, the anisotropic stress due to free streaming
neutrinos compensates for primordial magnetic fields. Therefore, the
significant growth of the PMFs sourced curvature perturbation
$\zeta_{B}(\textbf{k},\eta)$ occur between the time of magnetic field
generation $\eta_{B}$ and neutrino decoupling time $\eta_{\nu}$ on super-
horizon scales. During this growing regime, the PMF-induced primordial
curvature perturbation on super-horizon scales ($k\eta\ll 1$) on comoving
slices is given by Shaw and Lewis (2010); Durrer and Neronov (2013);
Subramanian (2016)

$\displaystyle\zeta_{B}(\textbf{k},\eta)=-\frac{1}{3}\xi(\eta)R_{\gamma}(\eta)\Pi_{B}(\textbf{k},\eta)$
(7)

where

$\displaystyle\xi(\eta)=\begin{cases}\log\left(\frac{\eta}{\eta_{B}}\right)+\frac{\eta_{B}}{2\eta}-\frac{1}{2}&\qquad(\eta_{B}<\eta<\eta_{\nu})\\\
\log\left(\frac{\eta_{\nu}}{\eta_{B}}\right)+\left(\frac{5}{8R_{\nu}}-1\right)&\qquad(\eta>\eta_{\nu})\end{cases}$
(8)

and $R_{\gamma}(\eta)=\frac{\rho_{\gamma}}{\rho_{\rm tot}}$ refers to the
energy fraction of the photon component in the total radiation. As we can see
the passive mode depends logarithmically on the magnetic field generation
epoch $\eta_{B}$ which in our case is determined by the scale of inflation.
The dimensionless power-spectrum $\mathcal{P}_{\zeta_{B}}(k,\eta)$ is related
to the two-point correlation function of comoving curvature perturbation as

$\displaystyle\langle\zeta_{B}(\textbf{k},\eta)\zeta_{B}^{*}(\textbf{p},\eta)\rangle=(2\pi)^{3}\delta^{3}(\textbf{k}-\textbf{p})\frac{2\pi^{2}}{k^{3}}\mathcal{P}_{\zeta_{B}}(k,\eta).$
(9)

Using Eq.(II.1) and Eq.(7), we can calculate the two-point correlation
function of the comoving curvature perturbation as

$\displaystyle\langle\zeta_{B}(\textbf{k},\eta)\zeta_{B}^{*}(\textbf{p},\eta)\rangle$
$\displaystyle=\frac{1}{9}\xi^{2}(\eta)R_{\gamma}^{2}(\eta)\langle\Pi_{B}(\textbf{k},\eta)\Pi_{B}^{*}(\textbf{p},\eta)\rangle$
$\displaystyle=\delta^{3}(\textbf{k}-\textbf{p})\frac{1}{9}\xi^{2}(\eta)R_{\gamma}^{2}(\eta)\left(\frac{9}{16\pi\rho_{\gamma,0}^{2}}\right)\left(\frac{g_{*s}(T)}{g_{*s}(T_{0})}\right)^{8/3}\int
dk^{\prime}{k^{\prime}}^{2}P_{B}(k^{\prime})\int_{-1}^{1}d\gamma
P_{B}(|\textbf{k}-\mathbf{k}^{\prime}|)$
$\displaystyle{}\qquad\qquad\times\left[\frac{k^{2}(1+\gamma^{2})+kk^{\prime}(2\gamma-6\gamma^{3})+{k^{\prime}}^{2}(5+9\gamma^{4}-12\gamma^{2})}{4(k^{2}+{k^{\prime}}^{2}-2\gamma
kk^{\prime})}\right]$ (10)

By comparing Eq.(9) and Eq.(II.1), we obtain the dimensionless power-spectrum
for the passive comoving curvature perturbation mode as

$\displaystyle\frac{2\pi^{2}}{k^{3}}\mathcal{P}_{\zeta_{B}}(k,\eta)$
$\displaystyle=\frac{1}{(2\pi)^{3}}\xi^{2}(\eta)R_{\gamma}^{2}\left(\frac{1}{16\pi\rho_{\gamma,0}^{2}}\right)\left(\frac{g_{*s}(T)}{g_{*s}(T_{0})}\right)^{8/3}\int
dk^{\prime}{k^{\prime}}^{2}P_{B}(k^{\prime})\int_{-1}^{1}d\gamma
P_{B}(|\textbf{k}-\mathbf{k}^{\prime}|)$
$\displaystyle{}\qquad\qquad\times\left[\frac{k^{2}(1+\gamma^{2})+kk^{\prime}(2\gamma-6\gamma^{3})+{k^{\prime}}^{2}(5+9\gamma^{4}-12\gamma^{2})}{4(k^{2}+{k^{\prime}}^{2}-2\gamma
kk^{\prime})}\right].$ (11)

Apart from the factor accounting for the change of the relativistic degrees of
freedom, this result recovers the one given in Saga  _et al._ (2020). Eq.
(II.1) will be used for estimating the variance of the magnetized overdense
region, which determines the properties of the PBH formation.

###  II.2 Density contrast and its variance

The curvature perturbation $\zeta_{B}$ is not a quantity that directly
represents the local (=horizon-scale) spatial curvature. For the purpose of
computing the PBH abundance, the density contrast on the comoving slice
$\delta_{B}$ is more suitable quantity than $\zeta_{B}$ Young  _et al._
(2014).

In the radiation-dominated epoch, the relation between comoving curvature
perturbation $\zeta_{B}$ and density contrast $\delta_{B}(\eta,k)$ is given by
Green  _et al._ (2004):

$\displaystyle\delta_{B}(\textbf{k},\eta)=\frac{4}{9}\left(\frac{k}{aH}\right)^{2}\zeta_{B}(\textbf{k},\eta)=-\frac{4}{27}\left(\frac{k}{aH}\right)^{2}\left[\log\left(\frac{\eta}{\eta_{B}}\right)+\frac{\eta_{B}}{2\eta}-\frac{1}{2}\right]R_{\gamma}(\eta)\Pi_{B}(\textbf{k},\eta)$
(12)

where we have used Eq. (7) in obtaining the last equality. The smoothed
density perturbation over the comoving radius $R$ is given by

$\displaystyle\delta_{B}(\textbf{x},\eta;R)=\int
d^{3}x^{\prime}W_{R}(|\textbf{x}-\textbf{x}^{\prime}|)\delta_{B}(\textbf{x}^{\prime},\eta)=\int\frac{d^{3}k}{(2\pi)^{3}}W(kR)\delta_{B}(\textbf{k},\eta)e^{i\textbf{k}\cdot\textbf{x}}$
(13)

where $W(kR)=e^{-k^{2}R^{2}/2}$ is the Window function. In calculating the PBH
abundance, the smoothing scale $R$ is usually taken to be the comoving horizon
scale at the PBH formation Sasaki  _et al._ (2018); Carr and Kuhnel (2020).
The variance is related to the power spectrum of the density perturbation by
the relation Nakama and Suyama (2016):

$\displaystyle\sigma^{2}_{B}(R,\eta)\equiv\langle\delta_{B}(\textbf{x},\eta;R)^{2}\rangle-\langle\delta_{B}(\textbf{x},\eta;R)\rangle^{2}=\int\frac{dk}{k}W^{2}(kR)\mathcal{P}_{\delta_{B}}(k,\eta),$
(14)

where $\mathcal{P}_{\delta_{B}}(k,\eta)$ is the dimensionless power spectrum
of $\delta_{B}$. The fact that magnetic fields are divergenceless, and from
the Eqns.(2, 5), the mean value of the curvature perturbation induced by PMFs
vanishes in the real space $\langle\zeta_{B}(\textbf{x})\rangle=0$. Thus, from
Eq.(14), we see that the mean density perturbation is zero,
$\langle\delta_{B}(\textbf{x},\eta;R)\rangle=0$. Using (12), the variance of
$\delta_{B}$ can be obtained as

$\displaystyle\sigma^{2}_{B}(R,\eta)$
$\displaystyle=\frac{16}{81}\int\frac{dk}{k}W^{2}(kR)(k\eta)^{4}\mathcal{P}_{\zeta_{B}}(k,\eta)$
(15)

We would like to mention that up to this point, the analysis is generic in the
sense that we have not used any specific form for the magnetic power spectrum.
However, to study the effect of the magnetic fields on the formation of PBHs,
we need to consider the magnetic field power spectrum. We will discuss this in
section IV.

##  III Effects of PMF on overdense region

In this section, we discuss the effect of PMFs on the overdense region and
study the consequences which are particularly important for PBH formation.

###  III.1 Effect of the magnetic fields on the Jeans length

To determine the PBH abundance from the given variance of the density contrast
(i.e., Eq. (15)), knowledge of the threshold of the PBH formation is needed.
In the case of the standard scenario where PBH forms purely from perturbation
of radiation, the threshold has been investigated in detail both analytically
and numerically Carr (1975); Shibata and Sasaki (1999); Musco  _et al._
(2005); Harada  _et al._ (2013); Nakama  _et al._ (2014); Musco (2019); Young
(2019); Escrivà (2020). The simplest version, which was derived in Carr (1975)
based on the Newtonian gravity, gives the threshold value $\delta_{\rm
th}=c_{s}^{2}(=\frac{1}{3})$ and it still serves as a physically reasonable
estimate of threshold at least qualitatively. In what follows, we provide a
physical argument that presence of the magnetic field changes the threshold
and derive the correction to the threshold coming from the magnetic field
based on the Jeans criterion in Newtonian dynamics. General relativistic
effect will change our result, but we believe that our correction captures the
essential features of the effect caused by the magnetic field and is valid
within a factor of ${\cal O}(1)$ uncertainty.

According to Jeans criterion, the perturbations (amplitude) would grow
(against the radiation pressure) if the size of the perturbed region is
greater than the Jeans length $\lambda_{J}$. If the uniform magnetic field
pervades the medium (which is assumed to be highly conducting), the dispersion
relation for the perturbations changes depending on the angle between the
direction of propagation of perturbations and the magnetic field vector B. For
example, the dispersion relation for the modes propagating perpendicular to
the direction the magnetic field is given by Chandrasekhar (1961)

$\displaystyle\omega^{2}=\frac{B^{2}k^{2}}{4\pi\rho}+c_{s}^{2}k^{2}-4\pi
G\rho$ (16)

where $B$ is the magnetic field strength, $c_{s}$ is the sound speed, $\rho$
is the energy density of the Universe. Gravitational instability occurs for
modes having negative $\omega^{2}$. Using the modified dispersion relation in
Eq.(16), we can obtain the modified Jeans length referred to as _magnetic
Jeans length_ $\lambda_{B}$, which for the modes traveling perpendicular to
the direction of the magnetic field is given by Chandrasekhar (1961);
Strittmatter (1966)

$\displaystyle\lambda_{B}\equiv\lambda_{J}\,\sqrt{1+\frac{2\rho_{B}}{c_{s}^{2}\rho}}=\lambda_{J}\,\sqrt{1+\frac{2r_{B}}{c_{s}^{2}}}.$
(17)

In what follows, we adopt this equation as the representative Jeans length in
the presence of the magnetic field. This is in sharp contrast to Ref. Saga
_et al._ (2020) in which the effect of the magnetic field on the Jeans length
was neglected. To the best of our knowledge, this paper is the first to
incorporate the effect of the magnetic field as $\lambda_{B}$ defined above in
the context of PBH formation. Note that the additional term in the above
equation is similar to the effect of rotation on the Jeans length discussed in
He and Suyama (2019). As we can see, the magnetic field increases the Jeans
length. For instance, for $r_{B}=0.1$, we have
$\left.\lambda_{B}\right|_{r_{B}=0.1}=1.265\lambda_{J}$. Although the increase
is only by $26\%$, this does not necessarily mean that the impact of the
magnetic field on the PBH abundance is mild as well since the abundance is
typically exponentially sensitive to the Jeans length. We will later discuss
this effect in detail in section IV.

###  III.2 Estimation of threshold density contrast in magnetized environment

In the previous subsection, we have shown that the presence of a uniform
magnetic field affects the Jeans criterion by increasing the Jeans length (17)
by a factor that depends on the magnetic field strength. To estimate the
effect of the magnetic field on the threshold density contrast, we first
briefly review how the threshold is derived in the standard case (i.e., no
magnetic field) Carr (1975). To this end, let us consider the overdense region
with the evolution given by Friedmann equation for the closed Universe Sasaki
_et al._ (2018); Carr and Kuhnel (2020)

$\displaystyle H^{2}=\frac{8\pi G}{3}\rho-\frac{1}{a^{2}}$ (18)

where $\rho\equiv\bar{\rho}(1+\delta)$ with $\bar{\rho}$ being the unperturbed
background energy density of the Universe and
$\delta\equiv(\rho-\bar{\rho})/\rho=\delta\rho/\rho$ is the density contrast.
Because of the curvature term, the overdense region turns from the expansion
phase to the contraction phase; we refer to this turnaround as the maximum
expansion and add the subscript “max” to any quantities evaluated at that
time. The time when the proper size of the overdense region becomes equal to
the Hubble horizon (horizon re-entry) is denoted by “hc” which refers to
horizon-crossing. The proper size of the overdense region at this maximum
expansion time is $\left(a_{\rm max}/a_{\rm hc}\right)H_{\rm hc}^{-1}$ and the
Jeans length evaluated at the maximum expansion time is given by:

$\displaystyle\lambda_{J}=\frac{c_{s}}{\sqrt{G\rho_{\rm
max}}}=\frac{c_{s}}{\sqrt{G\bar{\rho}_{\rm max}(1+\delta_{\rm max})}}\simeq
c_{s}a_{\rm max}.$ (19)

From Eq.(18), we have $a_{\rm max}=\sqrt{\frac{3}{8\pi G\rho_{\rm max}}}$.
Imposing the Jeans criterion at maximum expansion time gives the condition for
the formation of PBHs

$\displaystyle\frac{a_{\rm max}}{a_{\rm hc}}\frac{1}{H_{\rm
hc}}>\lambda_{J}\simeq c_{s}a_{\rm max}.$ (20)

The relation between $a_{\rm max}$ and $a_{\rm hc}$ in terms of density
contrast at horizon crossing is given by

$\displaystyle\frac{a_{\rm max}}{a_{\rm hc}}=\sqrt{\frac{1+\delta_{\rm
hc}}{\delta_{\rm hc}}},$ (21)

and constancy of $\rho(1+\delta)a^{4}$ gives

$\displaystyle H_{\rm hc}a_{\rm max}=\frac{\sqrt{1+\delta_{\rm
hc}}}{\delta_{\rm hc}}.$ (22)

Plugging these relations in Jeans criterion (20), we obtain the threshold
density contrast Carr (1975)

$\displaystyle\delta_{\rm hc}>\delta_{\rm th}=c_{s}^{2}.$ (23)

Now, let us include the presence of the magnetic fields in the above picture,
and investigate how the presence of uniform magnetic fields in the Universe
changes the threshold density contrast (23). As we have discussed in the
previous subsection, the presence of the uniform magnetic fields modifies the
Jeans length which is given by Eq.(17). Note that even in the presence of
uniform magnetic fields, the original argument that the size of the overdense
region at the maximum expansion time needs to be greater than the Jeans length
to form a black hole should still hold. This means the condition for the PBH
formation in the presence of the magnetic field is obtained by replacing the
original Jeans length with $\lambda_{B}$. Following the calculations as done
previously but with $\lambda_{B}$ on RHS of Eq.(20), we obtain

$\displaystyle\frac{a_{\rm max}}{a_{\rm hc}}>c_{s}(a_{\rm max}H_{\rm
hc})\sqrt{1+\frac{2\rho_{B}}{c_{s}^{2}\rho}}$ (24)

where we used $\lambda_{J}=c_{s}a_{\rm max}$. Now, using Eq.(21) and Eq.(22)
in Eq.(24), we can calculate the threshold density contrast for the formation
of PBHs in the presence of uniform magnetic field as

$\displaystyle\delta^{B}_{\rm
hc}>c_{s}^{2}\left(1+\frac{2r_{B}}{c_{s}^{2}}\right)\quad\implies\quad\delta_{\rm
th}^{B}=c_{s}^{2}+2r_{B}.$ (25)

From Eq. (25), we can see that a sufficiently strong magnetic field
significantly enhances the threshold density contrast than the non-magnetized
case. This can be understood as follows: a sufficiently strong magnetic field
would increase the total pressure due to an extra contribution from magnetic
pressure, which would act against the gravitational attraction. Therefore,
higher density contrast is required to form PBHs from the gravitational
collapse of the magnetized overdense region.

Before closing this section, we would like to mention that the threshold (25)
has been derived analytically based on the simple Newtonian picture and it is
merely approximate. According to general relativistic simulations (without the
magnetic field) Musco  _et al._ (2005), the threshold is not only different
from $c_{s}^{2}$ but also varies in the range $0.4-0.5$ depending on the
density profile of the overdense region. We expect the contribution of the
magnetic field, which is given as $2r_{B}$ in (25), also possesses such ${\cal
O}(1)$ variations in the realistic cases. Yet, since there are no simulations
of PBH formation in the presence of the magnetic field, to keep consistency of
our analysis, we use the threshold (25) in the following analysis.

##  IV Computation of PBH abundance

In this section, we provide a formalism to compute the PBH abundance in the
context of our scenario.

###  IV.1 Spiky power spectrum of PMF

In what follows, for simplicity, we consider a Delta-function type power
spectrum given by

$\displaystyle P_{B}(k)=\frac{2\pi^{2}}{k^{3}}B_{p}^{2}k_{p}\delta(k-k_{p})$
(26)

where amplification is peaked at $k=k_{p}$ and $B_{p}$ is the comoving
magnetic field at the scale $k=k_{p}$. Physically, $B_{p}$ represents the
amplitude of the virtual magnetic field at present time obtained by simply
evolving the primordial magnetic field in an adiabatic manner throughout the
whole redshift range. In reality, magneto-hydrodynamics in general completely
changes the evolution of the magnetic field once the scale $k_{p}$ enters the
eddy turnover scale, and $B_{p}$ does not represent the actual magnitude of
the present-day magnetic field. Nevertheless, in order to make a direct
comparison with the previous work Saga  _et al._ (2020), we use this quantity
to place the upper limit on the primordial magnetic field.

Magnetic power spectrum with a sharp spike is realized in the models where
magnetic fields are generated due to the parametric resonance mechanism Byrnes
_et al._ (2012); Patel  _et al._ (2020); Sasaki  _et al._ (2023). In these
mechanisms, the large amplification of EM modes happens only for the specific
modes lying in the resonance band determined by peak wavenumber $k_{p}$ and
the model’s parameter space. For the Delta-function type magnetic field power-
spectrum (26), the dimensionless power-spectrum for the passive comoving
curvature perturbation (II.1) becomes

$\displaystyle\mathcal{P}_{\zeta_{B}}(k,\eta)$
$\displaystyle=\frac{1}{(2\pi)^{3}}\xi^{2}(\eta)R_{\gamma}^{2}\left(\frac{1}{32\pi\rho_{\gamma,0}^{2}}\right)\left(\frac{g_{*s}(T)}{g_{*s}(T_{0})}\right)^{8/3}$
(27)
$\displaystyle\qquad\times\frac{\pi^{2}k^{2}B_{p}^{4}\left(k^{4}-16k^{2}k_{p}^{2}+80k_{p}^{4}\right)\theta(2k_{p}-k)}{16k_{p}^{6}}$

Using the above Eq.(27) in Eq. (15), we can obtain the variance as

$\displaystyle\sigma_{B}^{2}(R,\eta)$
$\displaystyle=\frac{B_{p}^{4}}{5184\pi^{2}\rho_{\gamma,0}^{2}}\left(\frac{g_{*s}(T)}{g_{*s}(T_{0})}\right)^{8/3}\xi(\eta)^{2}R_{\gamma}^{2}\left(\frac{\eta}{R}\right)^{4}\left(\frac{e^{-4(k_{p}R)^{2}}}{(k_{p}R)^{6}}\right)$
$\displaystyle\times\left(e^{4k_{p}^{2}R^{2}}(3-12k_{p}^{2}R^{2}+20k_{p}^{4}R^{4})+4(k_{p}R)^{4}-16(k_{p}R)^{6}-64(k_{p}R)^{8}-3\right).$
(28)

At the time of the horizon reentry at which $\eta=R$, the variance becomes

$\displaystyle\sigma_{B}^{2}(R,\eta=R)\simeq 1.5\times
10^{-3}{\left(\frac{B_{p}}{1~{}{\rm\mu
G}}\right)}^{4}\left(\frac{g_{*s}(\eta)}{g_{*s}(\eta_{0})}\right)^{8/3}(\xi(\eta)R_{\gamma}(\eta))^{2}.$
(29)

Or in terms of $r_{B}$, this relation becomes

$\sigma_{B}^{2}(R,\eta=R)\simeq\frac{r_{B}^{2}}{8.5}{\xi(\eta)}^{2}.$ (30)

From the above equation, we see that due to the function
$\xi(\eta)=\log\left(\frac{\eta}{\eta_{B}}\right)+\frac{\eta_{B}}{2\eta}-\frac{1}{2}$,
the variance $\sigma_{B}^{2}(R,\eta)$ depends on the epoch of magnetic field
generation $\eta_{B}$. In what follows, we take $T_{B}=10^{14}{\rm GeV}$. We
would like to emphasize that the variance $\sigma_{B}^{2}(R,\eta)$ also
depends significantly on the variation of relativistic degrees of freedom. To
our knowledge, this effect is not discussed in the literature in the
estimation of $\sigma_{B}^{2}(R,\eta)$.

Figure 1: The behaviour of $\xi R_{\gamma}$ (left) and the variance
$\sigma_{B}^{2}/B_{\rm nG}^{4}$ (right) with respect to $k$.

The left panel of 1 shows $\xi(\eta)R_{\gamma}(\eta)$ as a function of
$k_{p}$, where the evaluation time is taken to be the horizon crossing
($k=aH$). The right panel shows the variance of the density contrast smoothed
over the horizon scale $R\sim(aH)^{-1}=\eta$ evaluated at the horizon
crossing. This behaviour is not seen in Ref. Saga  _et al._ (2020) as they
consider the constant value of the quantity $\xi R_{\gamma}$. However, in this
work, we incorporated the variation of relativistic degrees of freedom with
temperature by using the data provided in Ref. Saikawa and Shirai (2018)

###  IV.2 PBH abundance

To investigate the abundance of PBHs, we define the parameter $\beta$,
representing the mass fraction (energy density fraction) of PBHs at the time
of their formation Inomata  _et al._ (2018); Sasaki  _et al._ (2018):

$\displaystyle\left.\beta=\frac{\rho_{\rm PBH}}{\rho_{\rm tot}}\right|_{\rm
form}.$ (31)

Denoting $P(\delta)$ by the probability distribution function of
$\delta_{B}(R,R)$, the fraction $\beta$ can be written as

$\beta=\gamma\int_{\delta>\delta^{B}_{\rm th}}P(\delta)d\delta,$ (32)

where $\gamma\sim 0.2$ is the numerical factor which appears due to the effect
that the PBH mass is not equal to the horizon mass at the time of the horizon
re-entry 111 In the literature, often the upper limit of integration in Eq.
(36) is taken as $\delta_{\rm max}=\infty$ (see Ref.Carr and Kuhnel (2020);
Bagui  _et al._ (2023)). However, both give approximately the same result but
for $\delta_{\rm max}=1$ case, we need to restrict $r_{B}$ to avoid
$\delta_{\rm th}^{B}$ becoming greater than unity which leads to the negative
value of $\beta$ (unphysical).. The lower limit of integration is non-trivial
because $\delta_{\rm th}^{B}$ itself depends on $\delta_{B}$, which means the
lower limit is obtained as a solution of the algebraic equation

$\delta_{B}=c_{s}^{2}+2r_{B}(\delta_{B}).$ (33)

It is hard to precisely evaluate the second term on the right-hand side. Here
we make an approximate identification $\delta_{B}\simeq\sigma_{B}$ and use the
relation (30) to finally relate $\delta_{B}$ with $r_{B}$, by which the above
equation reduces to a simple linear equation (notice that $\sigma_{B}$ in the
present context refers to the one defined for an over-dense region that
collapses to PBH contrary to the other parts in this paper where $\sigma_{B}$
represents the quantity averaged over the cosmological scale.). Solving this
equation, we obtain (see appendix A)

$\delta_{B}=\frac{c_{s}^{2}}{\left(1-\frac{5.8}{\xi(\eta)}\right)}.$ (34)

We adopt the right-hand side of this equation as the lower limit of the
integration (32).

In Ref. Carr  _et al._ (2021), constraints on the abundance of PBHs over
different mass ranges are provided. Applying those constraints to the fraction
(32) in our scenario, we can place upper limit on the amplitude of the
primordial magnetic fields $B_{p}$ on the scale $k_{p}$ by making a
correspondence between $k_{p}$ and the PBH mass $M$ given as

$\displaystyle k_{p}\simeq 7.5\times
10^{5}\gamma^{1/2}\left(\frac{g_{*}}{10.75}\right)^{-1/6}\left(\frac{M}{30M_{\odot}}\right)^{-1/2}\rm{Mpc}^{-1}.$
(35)

For the Gaussian distribution function of $\delta_{B}(R,R)$, the fraction is
given by Sasaki  _et al._ (2018):

$\displaystyle\beta(M)$ $\displaystyle=\gamma\int_{\delta^{B}_{\rm
th}}^{1}d\delta~{}~{}\frac{1}{\sqrt{2\pi}\sigma_{B}}e^{-\frac{\delta^{2}}{2\sigma^{2}_{B}}}\simeq\frac{\gamma}{\sqrt{2\pi}}\left(\frac{\sigma_{B}}{\delta_{\rm
th}^{B}}\right)\text{exp}\left(-\frac{{\delta_{\rm
th}^{B}}^{2}}{2\sigma_{B}^{2}}\right).$ (36)

In reality, however, the density contrast is highly non-Gaussian since it
depends quadratically on the magnetic field which is Gaussian. Therefore, to
correctly estimate PBH abundance, it is important to consider the non-Gaussian
probability distribution function (PDF). In this work, we use the non-Gaussian
distribution function studied in Ref. Saga  _et al._ (2020), which is given
by:

$\displaystyle P_{\rm
NG}(x)=C\exp\left[-1.5\sqrt{1+\frac{x^{2}}{b^{2}}}~{}~{}\right]$ (37)

where $C$ and $b$ are free parameters which can be fixed by using the
following conditions

$\displaystyle\int_{-\infty}^{\infty}P_{\rm
NG}(x)dx=1\qquad\int_{-\infty}^{\infty}x^{2}P_{\rm NG}(x)dx=\sigma_{B}^{2},$
(38)

where $\sigma_{B}^{2}$ is given by Eq. (29), which incorporates the effects of
variation of relativistic degrees of freedom. Using the above conditions gives

$\displaystyle C=2.13\sigma_{B}^{-1}~{},\qquad b=0.84\sigma_{B}.$ (39)

In the next section, we show and compare the constraints on the amplitude of
the PMFs for the non-Gaussian distribution function with the Gaussian
distribution and also compare our constraints with Ref. Saga  _et al._ (2020).

##  V Constraints on the strength of PMFs from the abundance of PBHs

Figure 2: The constraints on $B_{p}$ for Gaussian (blue) and non-Gaussian
(black) PDFs, dashed curves are for $\delta_{\rm th}$. Constraints from Ref.
Saga  _et al._ (2020) are shown in gray curves, where solid is for non-
Gaussian and dotted is for Gaussian PDFs.

Having discussed the effects of magnetic fields on the PBH abundance for both
types of distribution functions, our aim is now to constrain the amplitude of
PMFs in the RD era. To do that, we present the constraints in terms of the
comoving magnetic field amplitude $B_{p}$.

As we mentioned in the previous section, we employ the result provided in Ref.
Carr  _et al._ (2021) as the observational upper limits on the PBH fraction
$\beta$ 222 In Ref. Carr  _et al._ (2021), the upper limits are given in terms
of the rescaled fraction $\beta^{\prime}$ defined by
$\displaystyle\beta^{\prime}(M)\equiv\gamma^{1/2}\left(\frac{g_{*}}{106.75}\right)^{-1/4}\left(\frac{h}{0.67}\right)^{-2}\beta(M).$
(40) . 2 shows the constraints on the amplitude of PMFs for both the Gaussian
probability distribution function and the non-Gaussian one given by Eq. (37).
We find that inclusion of the effect of magnetic fields in the estimation of
threshold density contrast ($\delta_{\rm th}^{B}$) affects the constraints for
both Gaussian and non-Gaussian PDFs in the same manner. This is expected
because the net effect of the magnetic field is simply to shift the threshold
from $c_{s}^{2}$ to the one given by Eq. (34) irrespective of the shape of the
probability distribution function.

From 2, for the non-Gaussian distribution function (37), we estimate
$B_{p}\lesssim 0.12-0.14{\rm\mu G}$ on scale $k_{p}\sim 10^{17}{\rm Mpc}^{-1}$
and $B_{p}\lesssim 0.07{\rm\mu G}$ on scale $k_{p}\sim 10^{8}-10^{11}{\rm
Mpc}^{-1}$. 2 also compares our constraints (black and blue curves) with Saga
et al. Saga  _et al._ (2020) (shown in gray curves). We would like to mention
that our constraints incorporate new factors which are not included in the
previous work Ref. Saga  _et al._ (2020): (i) effect of the magnetic field in
the threshold density contrast and (ii) variation of relativistic degrees of
freedom.

Having derived the constraints from PBHs abundance, let us briefly discuss the
necessary relations to incorporate other small-scale observational probes to
constrain the PMFs. For the constraints from GWB and magnetic reheating, we
follow Refs. Saga  _et al._ (2018a, b). However, as we have discussed in the
previous section, concerning the variation of relativistic degrees of freedom
with respect to temperature changes, the correct relation can be obtained by
using Eq.(4). As we have shown, including the variation of the relativistic
degrees of freedom significantly affects the estimation of PBHs abundance.
Since the stochastic GW background is also affected by the relativistic
degrees of freedom, we must modify the relevant equation obtained in Refs.
Saga  _et al._ (2018a, b) for the precise estimation of the upper limit on the
amplitude of PMFs. Within our framework, the energy density of GWB at the
present epoch (corresponding to Ref.Saga  _et al._ (2018b)) is given as

$\displaystyle\Omega_{\rm
GW}(\eta_{0},k)=\frac{R_{\gamma}(\eta)^{2}}{512\pi^{2}}\left(\frac{B^{2}}{\rho_{\gamma,0}}\right)^{2}\left(\frac{g_{*s}(T)}{g_{*s}(T_{0})}\right)^{8/3}\left(\frac{k}{H_{0}}\right)^{2}a_{\rm
eq}^{2}h_{T}^{2}(k,\eta_{\rm eq})\left(\frac{k}{k_{\rm
p}}\right)^{2}\left(1+\frac{k^{2}}{4k_{p}^{2}}\right)^{2}\Theta_{H}\left(1-\frac{k}{2k_{p}}\right)$
(41)

where $h_{T}(k,\eta_{\rm eq})$ is the transfer function of GWs given in Ref.
Saga  _et al._ (2018b). Note that Eq. (41) tells us that given the upper
bounds on $\Omega_{\rm GW}$, we can put an upper limit on the amplitude of
PMFs.

Figure 3: Upper limits on the primordial magnetic field $B_{p}$.

Following Ref. Saga  _et al._ (2018b), we focus on the direct detection
measurements of GWB from pulsar timing arrays (PTAs) and SKA, space-based
observatories like LISA, and ground-based observatories such as LIGO, SKA, and
LISA. Note that the PTA and LIGO bounds are obtained from the observed
results; however, SKA and LISA are expected upper bounds in the future.

Another interesting effect of the PMFs on small scales is the magnetic
reheating, recently proposed in Ref. Saga  _et al._ (2018a). This mechanism is
based on the fact that the dissipation of PMFs changes the baryon-to-photon
number ratio and can constrain the PMFs on small scales. Again, for this case
as well, within our framework, the injected energy fraction to the CMB energy
by decaying PMFs should also be modified to

$\displaystyle\frac{\Delta\rho_{\gamma}}{\rho_{\gamma}}=\frac{B^{2}}{8\pi\rho_{\gamma,0}}\left(\frac{g_{*s}(T)}{g_{*s}(T_{0})}\right)^{4/3}\left[e^{-\frac{2k^{2}}{k_{D}^{2}(z_{i})}}-e^{-\frac{2k^{2}}{k_{D}^{2}(z_{\mu})}}\right]$
(42)

where $z_{\mu}=2\times 10^{6}$ and $z_{\nu}=10^{9}$ are the redshifts
corresponding to the CMB distortion era and neutrino decoupling era Saga  _et
al._ (2018a). The damping scale $k_{D}$ is related to the redshift as
$k_{D}(z)=7.44\times 10^{-6}(1+z)^{3/2}\rm{Mpc^{-1}}$. 3 shows the comparison
of the constraints on the amplitude of the PMFs on small scales and comparison
with various observations such as GW, magnetic reheating, CMB distortion, and
BBN. We would like to mention that the constraints in 3 give a more precise
estimation of the upper bounds on the amplitude of PMFs than presented in
Ref.Saga  _et al._ (2020). This is because one should consider two key issues:
increase in the threshold density contrast due to PMFs and the variation of
relativistic degrees of freedom with temperature.

##  VI Conclusion and Discussion

There are early-universe models which predict production of primordial
magnetic fields on small scales. Those magnetic fields will dissipate in the
early times and do not exist in the present universe, which makes probing such
fields non-trivial. In this work, we have systematically studied the procedure
to constrain the amplitude of PMFs in the RD era from the abundance of PBHs.
Given that anisotropic stress due to the PMFs sources the primordial curvature
perturbations, for sufficiently strong PMFs the density perturbations might
become large enough to gravitationally collapse to form PBHs in the RD era.
However, we have shown that this picture is not completely correct because the
magnetic field increases the threshold value of the density contrast required
for the formation of PBHs. We computed this correction term to the threshold
density contrast (due to magnetic field) and revised the constraints by using
the new threshold. Furthermore, if we naively extrapolate the new threshold to
very small scales (for $k_{p}\gtrsim 10^{18}~{}{\rm Mpc}^{-1}$), the threshold
density contrast becomes greater than one, which might stop the formation of
PBHs. Since our estimation is based on the simple Newtonian analysis, more
sophisticated works based on general relativity are needed to reach a robust
conclusion.

To make the estimates more precise as compared to the previous studies, we
have also considered the variation of relativistic degrees of freedom with
respect to the formation epoch and shown that it affects the variance by
$\mathcal{O}(2)$ on smaller scales. Thus, by taking into account the effect of
magnetic fields on the threshold density contrast, variation of relativistic
degrees of freedom, and using the most updated constraints on PBHs formation
parameter, we estimated an upper-limit $0.12-0.14{\rm\mu G}$ on scale
$10^{17}{\rm Mpc}^{-1}$ and $0.07{\rm\mu G}$ on scale $10^{8}-10^{11}{\rm
Mpc}^{-1}$ for non-Gaussian distribution function. We also compared our
constraints with those of other small-scale probes such as GWB, CMB
distortion, magnetic reheating, and BBN. These constraints also help to
support or rule-out certain classes of early-time magnetogenesis models.

###### Acknowledgements.

AK thanks Avijit Chowdhury and Joseph P Johnson for discussions. AK would like
to thank IIT Bombay for the financial support. AK would also like to
acknowledge the financial support from the Indo-French Centre for the
Promotion of Advanced Research (CEFIPRA) grant no. 6704-4 under the
Collaborative Scientific Research Programme. This work was supported by the
JSPS KAKENHI Grant Number JP23K03411 (TS).

##  Appendix A Derivation of the parameter $r_{B}$

In this appendix, we calculate the correction to the threshold density
contrast due to magnetic field i.e., the quantity $r_{B}$. Let us consider the
definition of $r_{B}$ which incorporates the effect of variation of
relativistic DOF:

$\displaystyle r_{B}\equiv\frac{\rho_{B}}{\rho_{\rm
tot}}=\frac{\rho_{B}}{\rho_{\gamma}}\frac{\rho_{\gamma}}{\rho_{\rm
tot}}=\frac{B_{f}^{2}}{8\pi\rho_{\gamma,0}}\left(\frac{a(T)}{a(T_{0})}\right)^{4}\left(\frac{g_{*s}(T)}{g_{*s}(T_{0})}\right)^{4/3}R_{\gamma}(\eta)=\frac{B_{0}^{2}}{8\pi\rho_{\gamma,0}}\left(\frac{g_{*s}(T)}{g_{*s}(T_{0})}\right)^{4/3}R_{\gamma}(\eta).$
(43)

Figure 4: The behaviour of $\delta_{th}^{B}$.

Using the relation between density perturbation and curvature perturbation

$\displaystyle\delta_{B}(\textbf{k},\eta)=\frac{4}{9}\left(\frac{k}{aH}\right)^{2}\zeta_{B}(\textbf{k},\eta),$
(44)

with Eq.(14), we can obtain the approximate relation
$\delta_{B}(R,\eta)\approx\sqrt{\langle\delta_{B}(\textbf{x},\eta;R)^{2}\rangle}=\sigma_{B}(R,\eta)$.
Now, using this relation with Eq. (IV.1), we obtain

$\displaystyle r_{B}=2.9\frac{\delta_{B}}{\xi(\eta)}.$ (45)

Note that the lower limit of integration for the PBH formation parameter
$\beta$ itself depends on the integration variable $\delta_{B}$ through
$r_{B}$. In the integration, only $\delta_{B}$ satisfying the condition
$\delta_{B}>c_{s}^{2}+2r_{B}$ contributes to the fraction $\beta$. Therefore,
using Eq.(45) and Eq.(33), we obtain

$\displaystyle\delta_{B}>\frac{c_{s}^{2}}{\left(1-\frac{5.8}{\xi(\eta)}\right)}.$
(46)

This is an important relation where RHS gives the threshold density contrast
$\delta_{\rm th}^{B}$ required for the formation of PBHs due to magnetized
overdense region. To see the effect of the correction term due to magnetic
field, we plot Eq. (46) in 4, which shows that for $k\gtrsim 10^{18}{\rm
Mpc}^{-1}$ the threshold becomes greater than unity.

## References

  * Kronberg (1994) P. P. Kronberg, Rept. Prog. Phys. 57, 325 (1994).
  * Beck  _et al._ (1996) R. Beck, A. Brandenburg, D. Moss, A. Shukurov, and D. Sokoloff, Annual Review of Astronomy and Astrophysics 34, 155 (1996).
  * Widrow (2002) L. M. Widrow, Rev. Mod. Phys. 74, 775 (2002).
  * Vallée (2004) J. P. Vallée, New Astron. Rev. 48, 763 (2004).
  * Gaensler  _et al._ (2004) B. M. Gaensler, R. Beck, and L. Feretti, New Astronomy Reviews 48, 1003 (2004).
  * Giovannini (2004) M. Giovannini, Int. J. Mod. Phys. D 13, 391 (2004), arXiv:astro-ph/0312614 .
  * Barrow  _et al._ (2007) J. D. Barrow, R. Maartens, and C. G. Tsagas, Phys. Rept. 449, 131 (2007), arXiv:astro-ph/0611537 [astro-ph] .
  * Kandus  _et al._ (2011) A. Kandus, K. E. Kunze, and C. G. Tsagas, Phys. Rept. 505, 1 (2011), arXiv:1007.3891 [astro-ph.CO] .
  * Durrer and Neronov (2013) R. Durrer and A. Neronov, Astron. Astrophys. Rev. 21, 62 (2013), arXiv:1303.7121 [astro-ph.CO] .
  * Subramanian (2016) K. Subramanian, Rept. Prog. Phys. 79, 076901 (2016), arXiv:1504.02311 [astro-ph.CO] .
  * Neronov and Vovk (2010) A. Neronov and I. Vovk, Science 328, 73 (2010).
  * Turner and Widrow (1988) M. S. Turner and L. M. Widrow, Phys. Rev. D37, 2743 (1988).
  * Ratra (1992) B. Ratra, Astrophys. J. Lett. 391, L1 (1992).
  * Vachaspati (1991) T. Vachaspati, Phys. Lett. B 265, 258 (1991).
  * Dolgov (1993) A. D. Dolgov, Phys. Rev. D 48, 2499 (1993).
  * Martin and Yokoyama (2008) J. Martin and J. Yokoyama, JCAP 01, 025 (2008), arXiv:0711.4307 [astro-ph] .
  * Mukohyama (2016) S. Mukohyama, Phys. Rev. D 94, 121302 (2016).
  * Kushwaha and Shankaranarayanan (2020) A. Kushwaha and S. Shankaranarayanan, Phys. Rev. D 101, 065008 (2020), arXiv:1912.01393 [hep-th] .
  * Bamba  _et al._ (2021) K. Bamba, E. Elizalde, S. D. Odintsov, and T. Paul, JCAP 04, 009 (2021), arXiv:2012.12742 [gr-qc] .
  * Giovannini (2021) M. Giovannini, JCAP 11, 058 (2021), arXiv:2110.02632 [hep-th] .
  * Tripathy  _et al._ (2021) S. Tripathy, D. Chowdhury, R. K. Jain, and L. Sriramkumar, (2021), arXiv:2111.01478 [astro-ph.CO] .
  * Nandi (2021) D. Nandi, JCAP 08, 039 (2021), arXiv:2103.03159 [astro-ph.CO] .
  * Kamada  _et al._ (2021) K. Kamada, F. Uchida, and J. Yokoyama, JCAP 04, 034 (2021), arXiv:2012.14435 [astro-ph.CO] .
  * Kawasaki and Kusakabe (2012) M. Kawasaki and M. Kusakabe, Phys. Rev. D 86, 063003 (2012).
  * Jedamzik and Saveliev (2019) K. Jedamzik and A. Saveliev, Phys. Rev. Lett. 123, 021301 (2019).
  * Kamada and Long (2016) K. Kamada and A. J. Long, Phys. Rev. D 94, 063501 (2016), arXiv:1606.08891 [astro-ph.CO] .
  * Fujita and Kamada (2016) T. Fujita and K. Kamada, Phys. Rev. D 93, 083520 (2016), arXiv:1602.02109 [hep-ph] .
  * Kushwaha and Shankaranarayanan (2021) A. Kushwaha and S. Shankaranarayanan, Phys. Rev. D 104, 063502 (2021), arXiv:2103.05339 [hep-ph] .
  * Kunze and Komatsu (2014) K. E. Kunze and E. Komatsu, JCAP 01, 009 (2014), arXiv:1309.7994 [astro-ph.CO] .
  * Saga  _et al._ (2018a) S. Saga, H. Tashiro, and S. Yokoyama, Mon. Not. Roy. Astron. Soc. 474, L52 (2018a), arXiv:1708.08225 [astro-ph.CO] .
  * Saga  _et al._ (2018b) S. Saga, H. Tashiro, and S. Yokoyama, Phys. Rev. D 98, 083518 (2018b).
  * Saga  _et al._ (2020) S. Saga, H. Tashiro, and S. Yokoyama, JCAP 05, 039 (2020), arXiv:2002.01286 [astro-ph.CO] .
  * Mack  _et al._ (2002) A. Mack, T. Kahniashvili, and A. Kosowsky, Phys. Rev. D 65, 123004 (2002), arXiv:astro-ph/0105504 .
  * Suyama and Yokoyama (2012) T. Suyama and J. Yokoyama, Phys. Rev. D 86, 023512 (2012), arXiv:1204.3976 [astro-ph.CO] .
  * Shaw and Lewis (2010) J. R. Shaw and A. Lewis, Phys. Rev. D 81, 043517 (2010), arXiv:0911.2714 [astro-ph.CO] .
  * Young  _et al._ (2014) S. Young, C. T. Byrnes, and M. Sasaki, JCAP 07, 045 (2014), arXiv:1405.7023 [gr-qc] .
  * Green  _et al._ (2004) A. M. Green, A. R. Liddle, K. A. Malik, and M. Sasaki, Phys. Rev. D 70, 041502 (2004), arXiv:astro-ph/0403181 .
  * Sasaki  _et al._ (2018) M. Sasaki, T. Suyama, T. Tanaka, and S. Yokoyama, Class. Quant. Grav. 35, 063001 (2018), arXiv:1801.05235 [astro-ph.CO] .
  * Carr and Kuhnel (2020) B. Carr and F. Kuhnel,  (2020), 10.1146/annurev-nucl-050520-125911, arXiv:2006.02838 [astro-ph.CO] .
  * Nakama and Suyama (2016) T. Nakama and T. Suyama, Phys. Rev. D 94, 043507 (2016), arXiv:1605.04482 [gr-qc] .
  * 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] .
  * Nakama  _et al._ (2014) T. Nakama, T. Harada, A. G. Polnarev, and J. Yokoyama, JCAP 01, 037 (2014), arXiv:1310.3007 [gr-qc] .
  * Musco (2019) I. Musco, Phys. Rev. D 100, 123524 (2019), arXiv:1809.02127 [gr-qc] .
  * Young (2019) S. Young, Int. J. Mod. Phys. D 29, 2030002 (2019), arXiv:1905.01230 [astro-ph.CO] .
  * Escrivà (2020) A. Escrivà, Phys. Dark Univ. 27, 100466 (2020), arXiv:1907.13065 [gr-qc] .
  * Chandrasekhar (1961) S. Chandrasekhar,  _Hydrodynamic and hydromagnetic stability_ (1961).
  * Strittmatter (1966) P. A. Strittmatter, MNRAS 131, 491 (1966).
  * He and Suyama (2019) M. He and T. Suyama, Phys. Rev. D 100, 063520 (2019), arXiv:1906.10987 [astro-ph.CO] .
  * Byrnes  _et al._ (2012) C. T. Byrnes, L. Hollenstein, R. K. Jain, and F. R. Urban, JCAP 03, 009 (2012), arXiv:1111.2030 [astro-ph.CO] .
  * Patel  _et al._ (2020) T. Patel, H. Tashiro, and Y. Urakawa, JCAP 01, 043 (2020), arXiv:1909.00288 [astro-ph.CO] .
  * Sasaki  _et al._ (2023) M. Sasaki, V. Vardanyan, and V. Yingcharoenrat, Phys. Rev. D 107, 083517 (2023), arXiv:2210.07050 [astro-ph.CO] .
  * Saikawa and Shirai (2018) K. Saikawa and S. Shirai, JCAP 05, 035 (2018), arXiv:1803.01038 [hep-ph] .
  * Inomata  _et al._ (2018) K. Inomata, M. Kawasaki, K. Mukaida, and T. T. Yanagida, Phys. Rev. D 97, 043514 (2018), arXiv:1711.06129 [astro-ph.CO] .
  * Bagui  _et al._ (2023) E. Bagui _et al._ (LISA Cosmology Working Group), (2023), arXiv:2310.19857 [astro-ph.CO] .
  * 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] .