# Primordial black holes induced stochastic axion-photon oscillations in
primordial magnetic field

Hai-Jun Li

###### Abstract

Primordial black holes (PBHs) can be produced in the very early Universe due
to the large density fluctuations. The cosmic background of axion-like
particles (ALPs) could be non-thermally generated by PBHs. In this paper, we
investigate the ALPs emitted by ultra-light PBHs with the mass range $10\,{\rm
g}\lesssim M_{\rm PBH}\lesssim 10^{9}\,\rm g$, in which PBHs would have
completely evaporated before the start of Big Bang Nucleosynthesis (BBN) and
can therefore not be directly constrained. In this case, the minimal scenario
that ALPs could interact only with photons is supposed. We study the
stochastic oscillations between the ALPs and photons in the cosmic magnetic
field in detail. The primordial magnetic field (PMF) can be modelled as the
stochastic background field model with the completely non-homogeneous
component of the cosmic plasma. Using the latest stringent limits on PMF, we
show the numerical results of ALP-photon oscillation probability distributions
with the homogeneous and stochastic magnetic field scenarios. The PBH-induced
stochastic ALP-photon oscillations in the PMF may have the effects on some
further phenomena, such as the cosmic microwave background (CMB), the cosmic
X-ray background (CXB), and the extragalactic gamma-ray background (EGB).

##  1 Introduction

Primordial black holes (PBHs) have recently gathered significant interest as
the candidates of dark matter (DM) [1]. They are a kind of black holes that
may be produced due to the large density fluctuations ($\rm i.e.$, the
radiation pressure could not resist the gravitational collapse in the
overdense spherically symmetric regions surrounded by the compensating void)
in the very early Universe [2, 3, 4], see $\rm e.g.$ refs. [5, 6, 7] for
recent reviews.

PBHs could exist over a wide range of mass scales from the Planck mass to
$\sim\mathcal{O}(10^{10})$ solar masses [8]. They emit the energetic particles
through Hawking evaporation. PBHs larger than $M_{\rm
PBH}\sim\mathcal{O}(10^{15})\,\rm g$ may have not entirely evaporated by now,
which could account for a significant fraction of the cold DM [9]. There are
several constraints on the abundance of non-evaporated PBHs, such as the
evaporation constraints [10, 11], the lensing constraints [12, 13], the
dynamical constraints [14, 15], the cosmic structure constraints [16, 17], and
the gravitational wave constraints [18, 19], etc. On the contrary, those PBHs
with the mass range less than $M_{\rm PBH}\sim\mathcal{O}(10^{15})\,\rm g$ may
have completely evaporated at present, while it could also lead the observable
signals on the cosmological consequences. There are also various constraints
on the evaporating PBHs, such as the Big Bang Nucleosynthesis (BBN)
constraints [20, 10], the cosmic microwave background (CMB) constraints [21,
22, 23], the PBH explosion constraints [24, 25], and the extragalactic
(Galactic) gamma-ray background (EGB) and cosmic-ray constraints [26, 27, 28,
29, 30, 31], etc.

Additionally, the cosmic axion(-like) background (C$a$B) could be non-
thermally generated by PBHs [32]. One of the leading hypotheses is that
ultralight spin-0 bosons such as axion or axion-like particle (ALP) are
potential DM candidates if they are non-thermally produced in the early
Universe by misalignment mechanism [33, 34, 35]. In refs. [36, 37, 38], they
recently investigated the phenomenological consequences of axion and ALP as
the DM triggered by PBHs. Axion is originally proposed with the spontaneous
breaking of the global $U(1)$ PQ symmetry to solve the strong CP problem in
quantum chromodynamics (QCD) [39, 40, 41, 42], which is also defined as the
QCD axion. While ALP is predicted by the extensions of Standard Model (SM) in
the high energy, such as the string theory [43, 44], without having to solve
the strong CP problem. See $\rm e.g.$ refs. [45, 46, 47] for recent reviews of
QCD axion and ALP. ALP is the ultralight pseudo-Nambu-Goldstone boson (pNGB)
$a$ with the two-photon vertex

$\displaystyle\mathscr{L}\supset-\frac{1}{4}g_{a\gamma}aF_{\mu\nu}\tilde{F}^{\mu\nu}\,,$
(1.1)

where $g_{a\gamma}$ is the coupling constant of ALP-photon. The ALP-photon
oscillation in the astrophysical magnetic fields would lead to several
detectable effects, such as the TeV transparency of the Universe [48, 49, 50,
51], the anomalous stellar cooling problem [52, 53], the galaxy cluster center
soft X-ray excess [54, 55], and the XENON1T low energy electronic recoil event
excess [56, 57], etc.111The continuously updated ALP-photon limits can be
found in [58].

Magnetic field is ubiquitous in the Universe, which would appear wherever
plasmas and currents can be found. The primordial magnetic field (PMF) might
have been generated with the high temperature environment of the Big Bang
[59], which also has the effects on various cosmological phenomena, such as
the inflation [60, 61], the cosmological phase transitions [62, 63, 64], the
BBN [65, 66], the CMB temperature and polarization anisotropies [67, 68], the
large scale structure (LSS) formation [69, 70], and the cosmic gravity wave
background [71, 72, 73], see $\rm e.g.$ refs. [74, 75, 59] for recent reviews.
The most widely used PMF model is the stochastic background field model with
the non-helical and helical components. The latest upper limits on PMF are
derived by the Planck collaboration using the Planck 2015 data [76], with the
magnetic field strength $B\sim\mathcal{O}(1)\,\rm nG$ for the scale
$\lambda\sim\mathcal{O}(1)\,\rm Mpc$. We also note that the very recent
stringent limits on PMF with $B\sim 8.9-47\,\rm pG$ [77] derived from the CMB
radiation.

In this work, we are more concerned about the phenomenon of ALP-photon
coupling induced by PBHs in the PMF. We note that the impacts of ALP-photon
oscillations triggered by PBHs on cosmic X-ray background (CXB) and X-ray
excess in galaxy clusters are recently investigated in ref. [78]. Similar to
refs. [79, 78], we focus our attention on the ALPs emitted by PBHs with the
mass range [80]

$\displaystyle\mathcal{O}(10)\,{\rm g}\lesssim M_{\rm
PBH}\lesssim\mathcal{O}(10^{9})\,\rm g\,,$ (1.2)

in which PBHs would have completely evaporated before the start of BBN. In
this case, the minimal scenario that ALPs could interact only with photons is
supposed. Differently from [78], we just study the stochastic oscillations
between ALPs and photons induced by PBHs in the cosmic magnetic field, but not
intend to investigate its phenomenological consequences in this context. The
PMF model and ALP-photon oscillation in the magnetic field are discussed in
detail. We consider the PMF model with the stochastic background field, which
can be modelled as the completely non-homogeneous component of the cosmic
plasma. We show the ALP-photon oscillation probability in different cases for
comparisons. In addition to the CXB, the ALP-photon oscillations induced by
PBHs may also have the effects on CMB and EGB, etc.

The rest of this paper is structured as follows. In section 2, we introduce
the PBH formation and evaporation, including the Schwarzschild black holes and
Kerr black holes, and a brief description of ALPs emitted by PBHs. In section
3, we describe the PMF model and the latest limits on PMF. Section 4 is the
ALP-photon oscillations in the PMF, including the general ALP-photon
oscillation in the magnetic field and the numerical results. Our numerical
results of the ALP-photon oscillation probability in the homogeneous and
stochastic magnetic fields are also presented in this section. Finally, we
comment on our results and conclude in section 5.

##  2 Primordial black hole formation and evaporation

PBHs could be generated due to large density perturbations in the very early
Universe [2]. When these density fluctuations reenter the horizon, they can
collapse to form a black hole if they exceed a threshold [3, 4]. In this
paper, we consider the formation of PBHs in the radiation dominated epoch. The
PBH initial mass at the formation time $t_{\rm in}$ can be described by [10]

$\displaystyle M_{\rm PBH}=\frac{4\pi}{3}\frac{\gamma\rho_{R}(t_{\rm
in})}{H^{3}(t_{\rm in})}\,,$ (2.1)

with the gravitational collapse factor $\gamma\sim\mathcal{O}(1)$ in the
radiation dominated epoch, the SM radiation energy density $\rho_{R}$, and the
Hubble expansion rate $H$. Using the Friedmann equation

$\displaystyle H=\sqrt{\frac{\rho_{R}}{3M_{\rm Pl}^{2}}}\,,$ (2.2)

we can derive

$\displaystyle M_{\rm PBH}=\frac{4\pi M_{\rm Pl}^{2}}{H(t_{\rm in})}\,,$ (2.3)

where $M_{\rm Pl}$ ($\equiv 1/\sqrt{8\pi G}$) is the reduced Planck mass.

After formation, PBH evaporates by releasing all possible degrees of freedom
in nature, which is also called Hawking evaporation. The Hawking evaporation
time of PBH can be described by [81]

$\displaystyle t_{\rm ev}=\frac{160}{\pi g_{*}}\frac{M_{\rm PBH}^{3}}{M_{\rm
Pl}^{4}}\,,$ (2.4)

where $g_{*}\sim\mathcal{O}(100)$ is the number of relativistic degrees of
freedom contributing to the SM energy density. Considering that PBHs have
completely evaporated before the start of BBN, we have

$\displaystyle H(t_{\rm ev})>H(t_{\rm BBN})\,.$ (2.5)

With $H(t_{\rm ev})=1/(2t_{\rm ev})$, $H(t_{\rm BBN})=\sqrt{\rho_{\rm
BBN}/(3M_{\rm Pl}^{2})}$, and $\rho_{\rm BBN}^{1/4}\sim 1\,\rm MeV$, we can
derive $M_{\rm PBH}\lesssim 10^{9}\,\rm g\,$. On the other hand, the PBHs
should be formed after the end of inflation, then we have

$\displaystyle H(t_{\rm in})<H(t_{\rm inf})\,.$ (2.6)

With $\rho_{\rm inf}^{1/4}\lesssim 10^{16}\,\rm GeV$, which can be indicated
from the current upper limit on tensor-to-scalar ratio in the single-field
slow-roll models of inflation [82], then we can derive the condition $M_{\rm
PBH}\gtrsim 10\,\rm g$. Therefore, the PBH mass range we considered in this
work is [80]

$\displaystyle 10\,{\rm g}\lesssim M_{\rm PBH}\lesssim 10^{9}\,\rm g\,.$ (2.7)

###  2.1 Schwarzschild black holes

For the simplest black hole scenario Schwarzschild black holes, they are
spherically symmetric and can only be described by the mass. The horizon
temperature is defined by [83]

$\displaystyle T_{\rm PBH}=\frac{M_{\rm Pl}^{2}}{M_{\rm PBH}}\simeq
10^{13}\left(\frac{M_{\rm PBH}}{1\,{\rm g}}\right)^{-1}\rm GeV\,.$ (2.8)

After formation, they emit the energetic particles through Hawking
evaporation. The emission rate of particles with type $i$ and spin $s$ from a
black hole per units of energy $E$ and time $t$ is given by [81]

$\displaystyle\frac{{\rm d}N_{i}}{{\rm d}E{\rm d}t}=\frac{1}{2\pi}\sum_{\rm
dof}\frac{\Gamma_{s_{i}}(M_{\rm PBH},E)}{\exp\left({E/T_{\rm
PBH}}\right)-(-1)^{2s_{i}}}\,,$ (2.9)

where the sum is over the particle total multiplicity, and
$\Gamma_{s_{i}}(M_{\rm PBH},E)$ is the dimensionless absorption coefficient,
called greybody factor

$\displaystyle\Gamma_{s_{i}}(M_{\rm PBH},E)=\frac{\sigma_{s_{i}}(M_{\rm
PBH},E)E^{2}}{\pi}\,.$ (2.10)

The greybody factor represents the spectral distortion of the blackbody,
corresponding to the probability that spherical wave of the particles induced
by the thermal fluctuation of the vacuum at the black hole horizon escape the
gravitational well. The absorption cross-section $\sigma_{s_{i}}(M_{\rm
PBH},E)$ would approach to the geometric optics limit $\sigma_{g}=27\pi
G^{2}M_{\rm PBH}^{2}$ in the high energy $E\gg T_{\rm PBH}$. While it is a
complicated function of $E$, $M_{\rm PBH}$, and spin $s$ in the low energy
region.

With eq. (2.9), we can derive the mass loss rate of an evaporating PBH through
Hawking evaporation [84]

$\displaystyle\frac{{\rm d}M_{\rm PBH}}{{\rm d}t}=-\frac{f(M_{\rm
PBH})}{M_{\rm PBH}^{2}}\,,$ (2.11)

where $f(M_{\rm PBH})$ is the Page factor [84, 85]

$\displaystyle\begin{aligned} f(M_{\rm PBH})&=-M_{\rm PBH}^{2}\frac{{\rm
d}M_{\rm PBH}}{{\rm d}t}\\\ &=M_{\rm
PBH}^{2}\int_{0}^{+\infty}\frac{E}{2\pi}\sum_{i}\sum_{\rm
dof}\frac{\Gamma_{s_{i}}(M_{\rm PBH},E)}{\exp\left({E/T_{\rm
PBH}}\right)-(-1)^{2s_{i}}}{\rm d}E\,,\end{aligned}$ (2.12)

which represents the number of relativistic degrees of freedom that can be
emitted in the evaporation process as a function of the instantaneous black
hole mass.

We show the instantaneous spectra of spin-0 particle emitted by a
Schwarzschild black hole as a function of energy in figure 1 with the public
code BlackHawk [86, 87]. Three typical values of PBH mass are selected for
comparisons. The black, blue, and red lines represent the Schwarzschild PBH
masses $M_{\rm PBH}=10^{5}\,\rm g$, $10^{7}\,\rm g$, and $10^{9}\,\rm g$,
respectively.

(a) Schwarzschild black holes. (b) Kerr black holes with $a_{*}=0.50$. (c)
Kerr black holes with $a_{*}=0.75$. (d) Kerr black holes with $a_{*}=0.99$.

Figure 1: The instantaneous spectra of the spin-0 particle emitted by a black
hole as a function of energy. The different panels correspond to the
Schwarzschild black holes and Kerr black holes (with $a_{*}=0.5$, $0.75$, and
$0.99$). The black, blue, and red lines correspond to the PBH masses $M_{\rm
PBH}=10^{5}\,\rm g$, $10^{7}\,\rm g$, and $10^{9}\,\rm g$, respectively.

###  2.2 Kerr black holes

Kerr black holes is an extension of the Schwarzschild black holes with the
non-zero angular momentum, which could acquire the spin through accretion [88]
or mergers [89]. Here the horizon temperature can be described by

$\displaystyle T_{\rm PBH}=\frac{1}{4\pi M_{\rm
PBH}}\frac{\sqrt{1-a_{*}}}{1+\sqrt{1-a_{*}}}\,,$ (2.13)

where $a_{*}$ is the reduced black hole spin parameter (the spin $a\equiv
J/M_{\rm PBH}\in\left[0,M_{\rm PBH}\right]$)

$\displaystyle a_{*}\equiv\frac{J}{M_{\rm PBH}^{2}}\in\left[0,1\right]\,,$
(2.14)

with the total angular momentum of the black hole $J$. Then the emission rate
of the particles becomes

$\displaystyle\frac{{\rm d}N_{i}}{{\rm d}E{\rm d}t}=\frac{1}{2\pi}\sum_{\rm
dof}\frac{\Gamma_{s_{i}}^{lm}(M_{\rm
PBH},E,a_{*})}{\exp\left({E^{\prime}/T_{\rm PBH}}\right)-(-1)^{2s_{i}}}\,,$
(2.15)

where

$\displaystyle E^{\prime}=E-\frac{m}{2M_{\rm
PBH}}\frac{a_{*}}{1+\sqrt{1-a_{*}^{2}}}\,,$ (2.16)

with the particles total $l$ and axial $m$ angular momentum numbers,
$m\in\\{-l,...,l\\}$.

Using eq. (2.15), the mass loss rate of these PBHs can be described by

$\displaystyle\frac{{\rm d}M_{\rm PBH}}{{\rm d}t}=-\frac{f(M_{\rm
PBH},a_{*})}{M_{\rm PBH}^{2}}\,,$ (2.17)

with the Page factor [84]

$\displaystyle\begin{aligned} f(M_{\rm PBH},a_{*})&=-M_{\rm PBH}^{2}\frac{{\rm
d}M_{\rm PBH}}{{\rm d}t}\\\ &=M_{\rm
PBH}^{2}\int_{0}^{+\infty}\frac{E}{2\pi}\sum_{i}\sum_{\rm
dof}\frac{\Gamma_{s_{i}}^{lm}(M_{\rm
PBH},E,a_{*})}{\exp\left({E^{\prime}/T_{\rm PBH}}\right)-(-1)^{2s_{i}}}{\rm
d}E\,.\end{aligned}$ (2.18)

The evolution of $a_{*}$ is given by [90]

$\displaystyle\begin{aligned} \frac{{\rm d}a_{*}}{{\rm d}t}&=\frac{{\rm
d}(J/M_{\rm PBH}^{2})}{{\rm d}t}=\frac{1}{M_{\rm PBH}^{2}}\frac{{\rm d}J}{{\rm
d}t}-2\frac{J}{M_{\rm PBH}^{3}}\frac{{\rm d}M_{\rm PBH}}{{\rm d}t}\\\
&=a_{*}\frac{2f(M_{\rm PBH},a_{*})-g(M_{\rm PBH},a_{*})}{M_{\rm
PBH}^{3}}\,,\end{aligned}$ (2.19)

with

$\displaystyle\begin{aligned} g(M_{\rm PBH},a_{*})&=-\frac{M_{\rm
PBH}}{a_{*}}\frac{{\rm d}J}{{\rm d}t}\\\ &=-\frac{M_{\rm
PBH}}{a_{*}}\int_{0}^{+\infty}\frac{m}{2\pi}\sum_{i}\sum_{\rm
dof}\frac{\Gamma_{s_{i}}^{lm}(M_{\rm
PBH},E,a_{*})}{\exp\left({E^{\prime}/T_{\rm PBH}}\right)-(-1)^{2s_{i}}}{\rm
d}E\,.\end{aligned}$ (2.20)

For comparison, we also show the instantaneous spectra of spin-0 particle
emitted by a Kerr black hole in figure 1 with the black hole spins
$a_{*}=0.50$, $0.75$, and $0.99$. In these panels, the black, blue, and red
lines represent the Kerr PBH masses $M_{\rm PBH}=10^{5}\,\rm g$, $10^{7}\,\rm
g$, and $10^{9}\,\rm g$, respectively. As we can see, the spectral
distributions in the high energy regions would shift to the right for the
black holes with large spin values, which is caused by the competition between
the power-law behaviour of the greybody factor and the exponential cut-off of
the Boltzmann factor. The Boltzmann factor is not greatly affected by the
black hole spin, while the greybody factor is enhanced as the spin increases.
Note that the nonphysical “oscillations” of the instantaneous spectra are the
consequences of the sum over many angular modes for the Hawking radiation.
Since each one of these angular modes corresponds to one “peak” in the
spectrum, the sum over all angular modes would generate a fake oscillatory
distribution. Therefore, the apparent “oscillations” have the nonphysical
significance, even if they increase with the black hole spin.

###  2.3 ALPs emitted by PBHs

In this subsection, we briefly introduce the ALPs emitted by PBHs, see ref.
[78] for more details. The ALP number density obeys the Boltzmann equation

$\displaystyle\left[\frac{\partial}{\partial t}-H(t)k\frac{\partial}{\partial
k}+2H(t)\right]\frac{{\rm d}n_{\rm ALP}(k,t)}{{\rm d}k}=n_{\rm
PBH}(t)\phi_{0}[k,M_{\rm PBH}(t)]\,,$ (2.21)

with

$\displaystyle n_{\rm PBH}(t)=\frac{n_{\rm PBH}(t_{\rm in})}{a^{3}(t)}\,,$
(2.22)

where $\phi_{0}$ is the emission rate of ALP, $k$ is the modulus of ALP
momentum, $n_{\rm PBH}$ and $n_{\rm ALP}$ are the number densities of PBH and
ALP, respectively. The solution of eq. (2.21) is

$\displaystyle\frac{{\rm d}n_{\rm ALP}(k^{\prime},t_{\rm ev})}{{\rm
d}k}=n_{\rm PBH}(t_{\rm ev})\int_{t_{\rm in}}^{t_{\rm
ev}}\frac{a^{\prime}}{a(t)}\phi_{0}[k^{\prime}\frac{a^{\prime}}{a(t)},M_{\rm
PBH}(t)]{\rm d}t\,,$ (2.23)

where $a^{\prime}$ corresponds to the PBH evaporation time $t_{\rm ev}$.
Considering the very light ALP, we have $k\simeq E$ with the ALP energy $E$.
The ALP energy density per unit of energy can be described by

$\displaystyle\frac{{\rm d}\rho_{\rm ALP}}{{\rm d}E}=E\frac{{\rm d}n_{\rm
ALP}(E,t)}{{\rm d}k}\,.$ (2.24)

Here we suppose the minimal scenario that ALPs could interact only with
photons. Considering the ALP to photon oscillations in the magnetic field, the
photon flux then can be described by

$\displaystyle\frac{{\rm d}n_{\gamma}(E,t)}{{\rm d}E}=P_{a\gamma}\frac{{\rm
d}n_{\rm ALP}(E,t)}{{\rm d}E}\,,$ (2.25)

where $P_{a\gamma}$ is the ALP-photon oscillation probability. The detailed
ALP-photon oscillations in the cosmic magnetic field can be found in section
4.

##  3 Primordial magnetic field

In this section, we describe the PMF model and the latest upper limits on PMF,
which is intended to discuss the ALP-photon oscillations in the PMF in next
section. One simplest PMF is the homogeneous magnetic field model [91, 92],
which can not be included in the homogeneous and isotropic expanding Universe
filled with the cosmic plasma. While the most widely accepted PMF model is the
stochastic background field model, which can be modelled as the completely
non-homogeneous component of the cosmic plasma [93, 94, 76]. In this case, the
contributions of energy density and anisotropic stress at the homogeneous
level can be neglected, the energy momentum tensor (EMT) components are on the
same basis as the cosmic perturbations.

###  3.1 Magnetic field model

The magnetic field fluctuation for the stochastic background field with the
non-helical and helical components can be described by the two-point
correlation functions in the Fourier space [94, 76]

$\displaystyle\langle
B_{i}(\textbf{k})B^{*}_{j}(\textbf{k}^{\prime})\rangle_{1}=\frac{\left(2\pi\right)^{3}}{2}\delta^{(3)}(\textbf{k}-\textbf{k}^{\prime})P_{ij}P_{B}(k)\,,$
(3.1) $\displaystyle\langle
B_{i}(\textbf{k})B^{*}_{j}(\textbf{k}^{\prime})\rangle_{2}=\frac{\left(2\pi\right)^{3}}{2}\delta^{(3)}(\textbf{k}-\textbf{k}^{\prime})\left[P_{ij}P_{B}(k)+i\epsilon_{ijl}\hat{k}_{l}P_{H}(k)\right]\,,$
(3.2)

where $P_{ij}=\delta_{ij}-\hat{k}_{i}\hat{k}_{j}$ is the projection operator,
$\hat{k}_{i}$ is the cartesian component of the normalized wave vector,
$\epsilon_{ijl}$ is the total anti-symmetric tensor, $P_{B}$ and $P_{H}$ are
the non-helical and helical components of the magnetic field power spectra

$\displaystyle P_{B}(k)=A_{B}k^{n_{B}}\,,$ (3.3) $\displaystyle
P_{H}(k)=A_{H}k^{n_{H}}\,,$ (3.4)

with the amplitudes $A_{B}$ and $A_{H}$, and the spectral indexes $n_{B}$ and
$n_{H}$.

For the non-helical component of PMF with the cut-off scale $k_{D}$, the
magnetic field fluctuation on the comoving scale of $\lambda=1\,\rm Mpc$ can
be described by

$\displaystyle B_{\lambda}^{2}=\frac{1}{2\pi^{2}}\int_{0}^{\infty}k^{2}{\rm
e}^{-k^{2}\lambda^{2}}P_{B}(k){\rm
d}k=\frac{A_{B}}{4\pi^{2}\lambda^{n_{B}+3}}\Gamma\left(\frac{n_{B}+3}{2}\right)\,.$
(3.5)

The cut-off scale is [95]

$\displaystyle k_{D}=\left(5.5\times
10^{4}\right)^{\frac{1}{n_{B}+5}}\left(\frac{B_{\lambda}}{\rm
nG}\right)^{-\frac{2}{n_{B}+5}}\left(\frac{2\pi}{\lambda/\rm
Mpc}\right)^{\frac{n_{B}+3}{n_{B}+5}}h^{\frac{1}{n_{B}+5}}\left(\frac{\Omega_{b}h^{2}}{0.022}\right)^{\frac{1}{n_{B}+5}}\rm
Mpc^{-1}\,,~{}~{}$ (3.6)

with the reduced Hubble constant $h$, and the baryon density parameter
$\Omega_{b}$. The magnetic field fluctuation of the helical component on the
comoving scale of $\lambda$ is given by

$\displaystyle
B_{\lambda}^{2}=\frac{\lambda}{2\pi^{2}}\int_{0}^{\infty}k^{3}{\rm
e}^{-k^{2}\lambda^{2}}|P_{H}(k)|{\rm
d}k=\frac{|A_{H}|}{4\pi^{2}\lambda^{n_{H}+3}}\Gamma\left(\frac{n_{H}+4}{2}\right)\,,$
(3.7)

with the helical spectral index $n_{H}>-4$.

The magnetic EMT can be regarded as the description of perturbations that
carry the energy density and anisotropic stress, and generate the Lorentz
force on the charged particles. For the non-helical component of PMF, the
relevant Fourier components of the EMT can be described by [76]

$\displaystyle\begin{aligned}
\Big{|}\rho_{B}(k)\Big{|}^{2}&=\frac{1}{1024\pi^{5}}\int
P_{B}(p)P_{B}(|\textbf{k}-\textbf{p}|)\left(1+\mu^{2}\right){\rm
d}^{3}p\,,\end{aligned}$ (3.8) $\displaystyle\begin{aligned}
\Big{|}L_{B}^{(\rm scalar)}(k)\Big{|}^{2}&=\frac{1}{128\pi^{2}a^{8}}\int
P_{B}(p)P_{B}(|\textbf{k}-\textbf{p}|)\left[1+\mu^{2}+4\gamma\beta\left(\gamma\beta-\mu\right)\right]{\rm
d}^{3}p\,,\end{aligned}$ (3.9) $\displaystyle\begin{aligned} \Big{|}\Pi^{(\rm
vector)}(k)\Big{|}^{2}&=\frac{1}{512\pi^{5}}\int
P_{B}(p)P_{B}(|\textbf{k}-\textbf{p}|)\left[\left(1-\gamma^{2}\right)\left(1+\beta^{2}\right)-\gamma\beta\left(\gamma\beta-\mu\right)\right]{\rm
d}^{3}p\,,\end{aligned}$ (3.10) $\displaystyle\Big{|}\Pi^{(\rm
tensor)}(k)\Big{|}^{2}=\frac{1}{512\pi^{5}}\int
P_{B}(p)P_{B}(|\textbf{k}-\textbf{p}|)\left(1+2\gamma^{2}+\gamma^{2}\beta^{2}\right){\rm
d}^{3}p\,,$ (3.11)

with
$\mu=\hat{\textbf{p}}\cdot(\textbf{k}-\textbf{p})/|\textbf{k}-\textbf{p}|$,
$\beta=\hat{\textbf{k}}\cdot(\textbf{k}-\textbf{p})/|\textbf{k}-\textbf{p}|$,
and $\gamma=\hat{\textbf{k}}\cdot\hat{\textbf{p}}$. While for the helical
component of PMF, the relevant components of the EMT are

$\displaystyle\begin{aligned}
\Big{|}\rho_{B}(k)\Big{|}^{2}&=\frac{1}{1024\pi^{5}}\int
P_{B}(p)P_{B}(|\textbf{k}-\textbf{p}|)\left(1+\mu^{2}\right){\rm d}^{3}p\\\
&-\frac{1}{512\pi^{5}}\int\mu P_{H}(p)P_{H}(|\textbf{k}-\textbf{p}|){\rm
d}^{3}p\,,\end{aligned}$ (3.12) $\displaystyle\begin{aligned}
\Big{|}L_{B}^{(\rm scalar)}(k)\Big{|}^{2}&=\frac{1}{128\pi^{2}a^{8}}\int
P_{B}(p)P_{B}(|\textbf{k}-\textbf{p}|)\left[1+\mu^{2}+4\gamma\beta\left(\gamma\beta-\mu\right)\right]{\rm
d}^{3}p\\\ &+\frac{1}{64\pi^{2}a^{8}}\int
pP_{H}(p)P_{H}(|\textbf{k}-\textbf{p}|)\left(\mu-2\gamma\beta\right){\rm
d}^{3}p\,,\end{aligned}$ (3.13) $\displaystyle\begin{aligned} \Big{|}\Pi^{(\rm
vector)}(k)\Big{|}^{2}&=\frac{1}{512\pi^{5}}\int
P_{B}(p)P_{B}(|\textbf{k}-\textbf{p}|)\left[\left(1-\gamma^{2}\right)\left(1+\beta^{2}\right)-\gamma\beta\left(\gamma\beta-\mu\right)\right]{\rm
d}^{3}p~{}~{}~{}\\\ &+\frac{1}{512\pi^{5}}\int
P_{H}(p)P_{H}(|\textbf{k}-\textbf{p}|)\left(\mu-\gamma\beta\right){\rm
d}^{3}p\,,\end{aligned}$ (3.14) $\displaystyle\begin{aligned} \Big{|}\Pi^{(\rm
tensor)}(k)\Big{|}^{2}&=\frac{1}{512\pi^{5}}\int
P_{B}(p)P_{B}(|\textbf{k}-\textbf{p}|)\left(1+2\gamma^{2}+\gamma^{2}\beta^{2}\right){\rm
d}^{3}p\\\ &+\frac{1}{128\pi^{5}}\int\gamma\beta
P_{H}(p)P_{H}(|\textbf{k}-\textbf{p}|){\rm d}^{3}p\,.\end{aligned}$ (3.15)

###  3.2 PMF limits

Here we list the latest PMF limits set by the Planck collaboration in 2015 and
by Jedamzik & Saveliev in 2019, which are defined as Planck 2015 and JS 2019,
respectively.

  * •

Planck 2015 [76]: The 95% confidence level comoving field amplitude limits of
the non-helical and helical components over the $1\,\rm Mpc$ region are
$B_{1\,\rm Mpc}<4.4\,\rm nG$ and $5.6\,\rm nG$, respectively. These results
are come from the CMB angular power spectra analysis by using the Planck
likelihood.

  * •

JS 2019 [77]: The 95% confidence level limits of the scale-invariant ($n=0$,
where $n$ is the spectral index) and a violet Batchelor spectrum ($n=5$) PMFs
are $B<47\,\rm pG$ and $8.9\,\rm pG$, respectively. These limits are derived
from the CMB radiation based on the non-observation of the magnetic field
induced clumping of baryonic density fluctuations during recombination epoch.

Since the limits of JS 2019 are more than one order of magnitude more
stringent than any prior PMF limits, we will use these results to discuss the
ALP-photon oscillations in the PMF in next section.

##  4 ALP-photon oscillations in primordial magnetic field

In this section, we discuss the ALP-photon oscillations in the PMF. We first
describe the general process of ALP-photon oscillation in the magnetic field.
Then we show the numerical results of the oscillation probability. Finally, we
briefly discuss the other effects on ALP-photon oscillations.

###  4.1 General ALP-photon oscillation

The ALP-photon interaction in the external magnetic field can be described by
the following effective Lagrangian [96]

$\displaystyle\mathscr{L}=-\frac{1}{4}F^{\mu\nu}F_{\mu\nu}+\frac{1}{2}\left(\partial_{\mu}a\partial^{\mu}a-m_{a}^{2}a^{2}\right)-\frac{1}{4}g_{a\gamma}aF_{\mu\nu}\tilde{F}^{\mu\nu}\,,$
(4.1)

with
$-\frac{1}{4}g_{a\gamma}aF_{\mu\nu}\tilde{F}^{\mu\nu}=g_{a\gamma}a\textbf{E}\cdot\textbf{B}$,
where $a$ is the ALP field, $m_{a}$ is the ALP mass, $g_{a\gamma}$ is the ALP-
photon coupling constant, $F_{\mu\nu}$ is the electromagnetic field tensor and
its dual tensor
$\tilde{F}^{\mu\nu}=\frac{1}{2}\epsilon^{\mu\nu\rho\sigma}F_{\rho\sigma}$, and
E, B are the local electric and magnetic field vectors.

Considering the ALP-photon system propagates along the direction $x_{3}$, the
system $\psi(x_{3})$ can be described by [97]

$\displaystyle\psi(x_{3})=\left(\begin{array}[]{c}A_{1}(x_{3})\\\
A_{2}(x_{3})\\\ a(x_{3})\end{array}\right)\,,$ (4.5)

where $A_{1}(x_{3})$ and $A_{2}(x_{3})$ correspond to the linear polarization
amplitudes of the photons in the perpendicular directions ($x_{1}\otimes
x_{2}$)

$\displaystyle|A_{1}\rangle=\left(\begin{array}[]{c}1\\\ 0\\\
0\end{array}\right)\,,\qquad|A_{2}\rangle=\left(\begin{array}[]{c}0\\\ 1\\\
0\end{array}\right)\,,\qquad|a\rangle=\left(\begin{array}[]{c}0\\\ 0\\\
1\end{array}\right)\,,$ (4.15)

and $a$ is the ALP. The equation of motion (EOM) for the ALP-photon system in
the magnetic field with the beam energy $E$ can be described by [96]

$\displaystyle\left(i\frac{\rm d^{2}}{{\rm
d}x_{3}^{2}}+E^{2}+2E\mathcal{M}(E,x_{3})\right)\psi(x_{3})=0\,.$ (4.16)

Considering the short-wavelength approximation, we have

$\displaystyle\begin{aligned} \left(\frac{\rm d^{2}}{{\rm
d}x_{3}^{2}}+E^{2}\right)\psi(x_{3})&=\left(i\frac{\rm d}{{\rm
d}x_{3}}+E\right)\left(-i\frac{\rm d}{{\rm d}x_{3}}+E\right)\psi(x_{3})\\\
&=2E\left(i\frac{\rm d}{{\rm d}x_{3}}+E\right)\psi(x_{3})\,.\end{aligned}$
(4.17)

Then we can rewrite the second-order eq. (4.16) to first-order as

$\displaystyle\left(i\frac{\rm d}{{\rm
d}x_{3}}+E+\mathcal{M}(E,x_{3})\right)\psi(x_{3})=0\,,$ (4.18)

where $\mathcal{M}(E,x_{3})$ is the mixing matrix

$\displaystyle\mathcal{M}(E,x_{3})=\left(\begin{array}[]{ccc}\Delta_{11}(E,x_{3}){}&\Delta_{12}(E,x_{3}){}&\Delta_{a\gamma,1}(x_{3})\\\
\Delta_{21}(E,x_{3}){}&\Delta_{22}(E,x_{3}){}&\Delta_{a\gamma,2}(x_{3})\\\
\Delta_{a\gamma,1}(x_{3}){}&\Delta_{a\gamma,2}(x_{3}){}&\Delta_{aa}(E)\end{array}\right)\,,$
(4.22)

with the terms

$\displaystyle\Delta_{11}(E,x_{3})$ $\displaystyle=$ $\displaystyle\Delta_{\rm
pl}(E,x_{3})+2\Delta_{\rm QED}(E,x_{3})+\Delta_{\rm CMB}(E)\,,$ (4.23)
$\displaystyle\Delta_{22}(E,x_{3})$ $\displaystyle=$ $\displaystyle\Delta_{\rm
pl}(E,x_{3})+\frac{7}{2}\Delta_{\rm QED}(E,x_{3})+\Delta_{\rm CMB}(E)\,,$
(4.24) $\displaystyle\Delta_{a\gamma,1}(x_{3})$ $\displaystyle=$
$\displaystyle\frac{1}{2}g_{a\gamma}B_{T,1}(x_{3})=0\,,$ (4.25)
$\displaystyle\Delta_{a\gamma,2}(x_{3})$ $\displaystyle=$
$\displaystyle\frac{1}{2}g_{a\gamma}B_{T,2}(x_{3})=\frac{1}{2}g_{a\gamma}B_{T}(x_{3})\,,$
(4.26) $\displaystyle\Delta_{aa}(E)$ $\displaystyle=$
$\displaystyle-\frac{m_{a}^{2}}{2E}\,,$ (4.27)

where $\Delta_{12}(E,x_{3})$ and $\Delta_{21}(E,x_{3})$ represent the Faraday
rotation terms and can be neglected. Here we assume that the transversal
magnetic field $B_{T}$ is aligned along $x_{2}$. We can rewrite the mixing
matrix $\mathcal{M}(E,x_{3})$ as

$\displaystyle\mathcal{M}(E,x_{3})=\left(\begin{array}[]{ccc}\Delta_{11}(E,x_{3}){}&0&0\\\
0&\Delta_{22}(E,x_{3}){}&\Delta_{a\gamma}(x_{3})\\\
0&\Delta_{a\gamma}(x_{3})&\Delta_{aa}(E)\end{array}\right)\,,$ (4.31)

with the terms

$\displaystyle\Delta_{\rm pl}(E,x_{3})$ $\displaystyle=$
$\displaystyle-\frac{\omega_{\rm pl}^{2}(x_{3})}{2E}\simeq-1.08\times
10^{-1}\left(\frac{n_{e}}{\rm cm^{-3}}\right)\left(\frac{E}{\rm
GeV}\right)^{-1}{\rm Mpc}^{-1}\,,$ (4.32)
$\displaystyle\Delta_{a\gamma}(x_{3})$ $\displaystyle=$
$\displaystyle\frac{1}{2}g_{a\gamma}B_{T}(x_{3})\simeq 1.52\times
10^{-2}\left(\frac{g_{a\gamma}}{10^{-11}\,\rm{GeV}^{-1}}\right)\left(\frac{B_{T}(x_{3})}{1\,\rm
nG}\right){\rm Mpc}^{-1}\,,$ (4.33) $\displaystyle\Delta_{aa}(E)$
$\displaystyle=$ $\displaystyle-\frac{m_{a}^{2}}{2E}\simeq-0.78\times
10^{2}\left(\frac{m_{a}}{10^{-9}\,\rm GeV}\right)^{2}\left(\frac{E}{\rm
GeV}\right)^{-1}{\rm Mpc}^{-1}\,,$ (4.34) $\displaystyle\Delta_{\rm
QED}(E,x_{3})$ $\displaystyle=$ $\displaystyle\frac{\alpha
E}{45\pi}\left(\frac{B_{T}(x_{3})}{B_{\rm cr}}\right)^{2}\simeq 4.10\times
10^{-12}\left(\frac{E}{\rm GeV}\right)\left(\frac{B_{T}(x_{3})}{1\,\rm
nG}\right)^{2}{\rm Mpc}^{-1}\,,$ (4.35) $\displaystyle\Delta_{\rm CMB}(E)$
$\displaystyle=$ $\displaystyle\rho_{\rm CMB}E\simeq 0.80\times
10^{-4}\left(\frac{E}{\rm GeV}\right){\rm Mpc}^{-1}\,.$ (4.36)

The term $\Delta_{\rm pl}(E,x_{3})$ represents the plasma effect when ALP-
photon system propagates in the plasma with the plasma frequency $\omega_{\rm
pl}=\sqrt{4\pi\alpha n_{e}/m_{e}}$, where $\alpha$ is the fine-structure
constant, $n_{e}$ is the free electron number density, and $m_{e}$ is the
electron mass. The term $\Delta_{\rm QED}(E,x_{3})$ represents the QED vacuum
polarization effect from the following Heisenberg-Euler-Weisskopf effective
Lagrangian [98]

$\displaystyle\mathscr{L}=\frac{2\alpha^{2}}{45m_{e}^{2}}\left[\left(\textbf{E}^{2}-\textbf{B}^{2}\right)^{2}+7\left(\textbf{E}\cdot\textbf{B}\right)^{2}\right]\,,$
(4.37)

with the critical magnetic field $B_{\rm cr}=m^{2}_{e}/|e|\simeq 4.41\times
10^{13}\rm\,G$. The term $\Delta_{\rm CMB}(E,x_{3})$ represents the CMB photon
dispersion effect with [99]

$\displaystyle\rho_{\rm CMB}\simeq 0.511\times 10^{-42}\,.$ (4.38)

Now considering the ALP-photon system propagates in a homogeneous magnetic
field $B_{T}$, the ALP-photon oscillation probability can be simply described
by

$\displaystyle\mathcal{P}_{a\gamma}(E,x_{3})=\left(\frac{g_{a\gamma}B_{T}L_{\rm
osc}(E)}{2\pi}\right)^{2}\sin^{2}\left(\frac{\pi x_{3}}{L_{\rm
osc}(E)}\right)\,,$ (4.39)

with the oscillation length

$\displaystyle\begin{aligned} L_{\rm
osc}(E)&=2\pi\left[\left(\Delta_{22}(E)-\Delta_{aa}(E)\right)^{2}+4\Delta_{a\gamma}^{2}\right]^{-\frac{1}{2}}\\\
&=2\pi\left[\left[\frac{|m_{a}^{2}-\omega_{\rm
pl}^{2}|}{2E}+E\left(\frac{7\alpha}{90\pi}\left(\frac{B_{T}}{B_{\rm
cr}}\right)^{2}+\rho_{\rm
CMB}\right)\right]^{2}+g_{a\gamma}^{2}B_{T}^{2}\right]^{-\frac{1}{2}}\,.~{}\end{aligned}$
(4.40)

###  4.2 Numerical results of ALP-photon oscillation probability

Using the PMF limits JS 2019 given in section 3.2 and the ALP-photon
oscillation process in the magnetic field discussed in section 4.1, we can
derive the oscillation probability for the ALP-photon system propagates in the
PMF. In this subsection, we show the numerical results of ALP-photon
oscillation probability with the homogeneous and stochastic magnetic field
scenarios, which is simulated with the domain-like structure.

In this case, the ALP-photon system in eq. (4.5) can be described by the
following density matrix

$\displaystyle\rho(x_{3})=\left(\begin{array}[]{c}A_{1}(x_{3})\\\
A_{2}(x_{3})\\\
a(x_{3})\end{array}\right)\otimes\left(\begin{array}[]{c}A_{1}(x_{3}),A_{2}(x_{3}),a(x_{3})\end{array}\right)^{*}\,,$
(4.45)

which satisfies the Von Neumann-like commutator equation

$\displaystyle i\frac{{\rm d}\rho(x_{3})}{{\rm
d}x_{3}}=\rho(x_{3})\mathcal{M}^{\dagger}(E,x_{3})-\mathcal{M}(E,x_{3})\rho(x_{3})\,.$
(4.46)

The solution of eq. (4.46) is

$\displaystyle\rho(x_{3})=\mathcal{T}(E,x_{3})\rho(0)\mathcal{T}^{\dagger}(E,x_{3})\,,$
(4.47)

with the initial density matrix $\rho(0)$

$\displaystyle\rho(0)=\frac{1}{2}\left(\begin{array}[]{ccc}1{}&0{}&0\\\
0{}&1{}&0\\\ 0{}&0{}&0\end{array}\right)\,,$ (4.51)

and the transport matrix $\mathcal{T}(E,x_{3})$

$\displaystyle\mathcal{T}(E,x_{3})=\prod^{N}_{i=1}\mathcal{T}(E_{i},x_{3,i})\,.$
(4.52)

Then the ALP-photon oscillation probability $\mathcal{P}_{a\gamma}$ for some
propagation distance can be described by

$\displaystyle\mathcal{P}_{a\gamma}=1-{\rm
Tr}\left[\left(\rho_{11}+\rho_{22}\right)\mathcal{T}(E,x_{3})\rho(0)\mathcal{T}^{\dagger}(E,x_{3})\right]\,,$
(4.53)

with

$\displaystyle{\rho}_{11}=\left(\begin{array}[]{ccc}1{}&0{}&0\\\ 0{}&0{}&0\\\
0{}&0{}&0\end{array}\right)\,,\qquad{\rho}_{22}=\left(\begin{array}[]{ccc}0{}&0{}&0\\\
0{}&1{}&0\\\ 0{}&0{}&0\end{array}\right)\,.$ (4.60)

(a) $E=10^{-3}\,\rm GeV$, $B_{T}=47\,\rm pG$. (b) $E=10^{-3}\,\rm GeV$,
$B_{T}=5\,\rm nG$. (c) $E=1\,\rm GeV$, $B_{T}=47\,\rm pG$. (d) $E=1\,\rm GeV$,
$B_{T}=5\,\rm nG$. (e) $E=10^{3}\,\rm GeV$, $B_{T}=47\,\rm pG$. (f)
$E=10^{3}\,\rm GeV$, $B_{T}=5\,\rm nG$.

Figure 2:  The distributions of ALP-photon oscillation probability in the
$m_{a}-g_{a\gamma}$ plane for several typical energy values. The top, middle,
and bottom panels correspond to $E=10^{-3}\,\rm GeV$, $1\,\rm GeV$, and
$10^{3}\,\rm GeV$, respectively. The left and right panels represent
$B_{T}=47\,\rm pG$ and $5\,\rm nG$, respectively.

(a) $E=10^{4}\,\rm GeV$, $B_{T}=47\,\rm pG$. (b) $E=10^{4}\,\rm GeV$,
$B_{T}=5\,\rm nG$. (c) $E=10^{5}\,\rm GeV$, $B_{T}=47\,\rm pG$. (d)
$E=10^{5}\,\rm GeV$, $B_{T}=5\,\rm nG$. (e) $E=10^{6}\,\rm GeV$,
$B_{T}=47\,\rm pG$. (f) $E=10^{6}\,\rm GeV$, $B_{T}=5\,\rm nG$.

Figure 3:  Same as figure 2 but for different energy values. The top, middle,
and bottom panels correspond to $E=10^{4}\,\rm GeV$, $10^{5}\,\rm GeV$, and
$10^{6}\,\rm GeV$, respectively. The left and right panels represent
$B_{T}=47\,\rm pG$ and $5\,\rm nG$, respectively.

Here we show the numerical results of the ALP-photon oscillation probability
for the homogeneous magnetic field with the PMF limits JS 2019. In this case,
the magnetic field strength and direction in each domain are same, which
corresponds to the scale-invariant limit $B<47\,\rm pG$ with the spectral
index $n=0$. In figure 2, we plot the oscillation probability distributions in
the ALP parameter ($m_{a}$, $g_{a\gamma}$) space for several typical values of
energy. The top, middle, and bottom panels correspond to $E=10^{-3}\,\rm GeV$,
$1\,\rm GeV$, and $10^{3}\,\rm GeV$, respectively. For comparison, the
oscillation probability distributions with the PMF limits Planck 2015 ($B\sim
5\,\rm nG$) are also shown. The left and right panels correspond to the
magnetic field strength $B_{T}=47\,\rm pG$ and $B_{T}=5\,\rm nG$,
respectively. The other typical values of energy ($E=10^{4}\,\rm GeV$,
$10^{5}\,\rm GeV$, and $10^{6}\,\rm GeV$) for the oscillation probability
distributions can be found in figure 3.

In these panels, the distributions of ALP-photon oscillation probability at
$0.01\lesssim\mathcal{P}_{a\gamma}\lesssim 0.50$ in the $m_{a}-g_{a\gamma}$
plane are shown, while the probability at $\mathcal{P}_{a\gamma}\lesssim 0.01$
could not be completely displayed. The oscillation probability distribution
changes significantly with the energy, and could be dramatically suppressed in
the energy above $\sim 10^{4}\,\rm GeV$. We also find the significant
difference between the strength $B_{T}=47\,\rm pG$ and $5\,\rm nG$ with the
homogeneous magnetic field scenario. Here we do not show the specific
oscillation probability, which can be found in next subsection with the
stochastic magnetic field scenario.

Figure 4:  Effective magnetic field strength distributions for the stochastic
magnetic field scenario with 100 domains for $1\,\rm Mpc$. The blue and red
lines correspond to the magnetic field strength in each domain with the two
stochastic magnetic field configurations. The dot-dashed black lines represent
$|B_{T}|=8.9\,\rm pG$.

####  4.2.1 Stochastic magnetic field

For the stochastic magnetic field scenario, if the transversal magnetic field
$B_{T}$ is not aligned along the direction of $x_{2}$ and forms the angle
$\phi$, the mixing matrix $\mathcal{M}(E,x_{3})$ in eq. (4.22) should be
converted as

$\displaystyle\mathcal{M}(E,x_{3})=\mathcal{V}^{\dagger}(\phi)\mathcal{M}^{\prime}(E,x_{3})\mathcal{V}(\phi)\,,$
(4.61)

with

$\displaystyle\mathcal{V}(\phi)=\left(\begin{array}[]{ccc}\cos\phi{}&-\sin\phi{}&0\\\
\sin\phi{}&\cos\phi{}&0\\\ 0&0&1\\\ \end{array}\right)\,.$ (4.65)

Then we can rewrite the mixing matrix $\mathcal{M}(E,x_{3},\phi)$ in eq.
(4.31) as

$\displaystyle\mathcal{M}(E,x_{3},\phi)=\left(\begin{array}[]{ccc}\Delta_{11}(E,x_{3})&0&\Delta_{a\gamma}(x_{3})\sin\phi\\\
0&\Delta_{22}(E,x_{3})&\Delta_{a\gamma}(x_{3})\cos\phi\\\
\Delta_{a\gamma}(x_{3})\sin\phi{}&\Delta_{a\gamma}(x_{3})\cos\phi{}&\Delta_{aa}(E)\end{array}\right)\,.$
(4.69)

(a) $m_{a}=10^{-10}\,\rm eV$, $g_{a\gamma}=10^{-9}\,\rm GeV^{-1}$. (b)
$m_{a}=10^{-8}\,\rm eV$, $g_{a\gamma}=10^{-9}\,\rm GeV^{-1}$. (c)
$m_{a}=10^{-10}\,\rm eV$, $g_{a\gamma}=10^{-10}\,\rm GeV^{-1}$. (d)
$m_{a}=10^{-8}\,\rm eV$, $g_{a\gamma}=10^{-10}\,\rm GeV^{-1}$.

Figure 5:  ALP-photon oscillation probability in the stochastic magnetic field
as a function of energy with the limit $B<8.9\,\rm pG$. Several typical ALP
parameter sets [$m_{a}=10^{-10}\,\rm eV$, $g_{a\gamma}=10^{-9}\,\rm GeV^{-1}$]
(a), [$m_{a}=10^{-8}\,\rm eV$, $g_{a\gamma}=10^{-9}\,\rm GeV^{-1}$] (b),
[$m_{a}=10^{-10}\,\rm eV$, $g_{a\gamma}=10^{-10}\,\rm GeV^{-1}$] (c), and
[$m_{a}=10^{-8}\,\rm eV$, $g_{a\gamma}=10^{-10}\,\rm GeV^{-1}$] (d) are
selected for comparisons. The brown and orange bounds represent the $1\sigma$
and $2\sigma$ contours with 1000 random magnetic field configurations,
respectively. The blue lines represent the oscillation probability with the
random one configuration.

Here we show the numerical results for the stochastic magnetic field scenario.
In this case, the magnetic field strength in every domain is same, while the
direction is random from the domain $i$ to domain $i+1$, which corresponds to
the PMF limit $B<8.9\,\rm pG$ with $n=5$. In figure 4, we plot the effective
magnetic field strength distributions with this PMF limit. The stochastic
magnetic field is simulated with 100 domains for $1\,\rm Mpc$. The blue and
red lines represent the effective magnetic field strength in each domain with
the two stochastic magnetic field configurations. The dot-dashed black lines
represent the magnetic field strength $|B_{T}|=8.9\,\rm pG$.

In figure 5, we show the ALP-photon oscillation probability as a function of
energy for the stochastic magnetic field scenario. We show the several typical
ALP parameter sets [$m_{a}=10^{-10}\,\rm eV$, $g_{a\gamma}=10^{-9}\,\rm
GeV^{-1}$], [$m_{a}=10^{-8}\,\rm eV$, $g_{a\gamma}=10^{-9}\,\rm GeV^{-1}$],
[$m_{a}=10^{-10}\,\rm eV$, $g_{a\gamma}=10^{-10}\,\rm GeV^{-1}$], and
[$m_{a}=10^{-8}\,\rm eV$, $g_{a\gamma}=10^{-10}\,\rm GeV^{-1}$] for
comparisons. Here the magnetic field is simulated with 1000 random magnetic
field configurations with the JS 2019 limit $B<8.9\,\rm pG$. The brown and
orange bounds represent the oscillation probability distributions with the
$1\sigma$ and $2\sigma$ contours, respectively. The blue lines represent the
oscillation probability with the random one magnetic field configuration.

As we can see, the ALP-photon oscillation probability for these typical ALP
parameter sets show the different distributions. For the same ALP mass
$m_{a}$, the oscillation probability can be dramatically enhanced by the
coupling $g_{a\gamma}$. Compared with $g_{a\gamma}$, the mass $m_{a}$ can
affect the oscillation probability in the low energy regions, while the
probability would be stable for the high energies, which can also be seen from
figure 3. Note that the stochastic magnetic field is simulated with the PMF
limit $B<8.9\,\rm pG$, the effect of magnetic field strength on ALP-photon
oscillation probability can be found from the homogeneous magnetic field
scenario as show in figures 2 and 3.

###  4.3 Other effects on oscillations

In order to obtain the current oscillation probability, the ALP-photon
oscillation in the Universe evolution process may be affected by some other
effects, such as the redshift and the photon absorption at the recombination
[79, 78]. Since we are more concerned about the phenomena of ALP-photon
oscillations in the PMF, we just present these effects in this subsection and
do not intend to show them in our numerical results.

Considering the expansion of the Universe, we have the notations $B_{T}\to
B_{T}\left(1+z\right)^{2}$, $E\to E\left(1+z\right)$, $l^{\prime}\to
l^{\prime}\left(1+z\right)^{-1}$, $n_{e}\to n_{e}\left(1+z\right)^{3}$,
$\Delta_{\rm pl}\to\Delta_{\rm pl}\left(1+z\right)^{2}$,
$\Delta_{a\gamma}\to\Delta_{a\gamma}\left(1+z\right)^{2}$,
$\Delta_{aa}\to\Delta_{aa}\left(1+z\right)^{-1}$, $\Delta_{\rm
QED}\to\Delta_{\rm QED}\left(1+z\right)^{5}$, and $\Delta_{\rm
CMB}\to\Delta_{\rm CMB}\left(1+z\right)^{5}$, where $z$ is the redshift, and
$l^{\prime}$ is the length of the magnetic field domain.

Additionally, the photon absorption by hydrogen and helium atoms formed at the
recombination should also be considered. We can rewrite the mixing matrix
$\mathcal{M}(E,x_{3})$ in eq. (4.31) as

$\displaystyle\mathcal{M}(E,x_{3})=\left(\begin{array}[]{ccc}\Delta_{11}(E,x_{3})-i\frac{\Gamma(E)}{2}{}&0&0\\\
0&\Delta_{22}(E,x_{3})-i\frac{\Gamma(E)}{2}{}&\Delta_{a\gamma}(x_{3})\\\
0&\Delta_{a\gamma}(x_{3})&\Delta_{aa}(E)\end{array}\right)\,.$ (4.73)

The term $\Gamma(E)$ represents the photon absorption rate,
$\Gamma(E)=\sigma_{\rm H}(E)n_{\rm H}+\sigma_{\rm He}(E)n_{\rm He}+\sigma_{\rm
KN}(E)n_{e}\,,$ with the photoelectric cross-sections of Hydrogen $\sigma_{\rm
H}$ and Helium $\sigma_{\rm He}$, and the Klein-Nishina cross-section for the
Compton effect on electron $\sigma_{\rm KN}$. The number densities of $\rm H$,
$\rm H_{e}$, and $e$ are $n_{\rm H}$, $n_{\rm He}$, and $n_{e}$, respectively.

Considering the above effects on ALP-photon oscillations, one could derive the
present-day oscillation probability $P^{\prime}_{a\gamma}$ and the present-day
photon flux

$\displaystyle\frac{{\rm d}n^{\prime}_{\gamma}(E,t_{0})}{{\rm
d}E}=P^{\prime}_{a\gamma}\frac{{\rm d}n^{\prime}_{\rm ALP}(E,t_{0})}{{\rm
d}E}\,.$ (4.74)

In this case, the ALP-photon oscillations induced by PBHs may have the effects
on some further phenomena for the photons in the different energy regions,
such as the CMB, the CXB [78], and the EGB.

##  5 Conclusions

In this paper, we have presented the stochastic ALP-photon oscillations
induced by PBHs in the PMF with the latest magnetic field limits. First we
introduce the PBH formation and evaporation, including the Schwarzschild black
holes and Kerr black holes. We investigate the ALPs emitted by ultra-light
PBHs with the mass range $10\,{\rm g}\lesssim M_{\rm PBH}\lesssim 10^{9}\,\rm
g$, in which PBHs would have completely evaporated before the start of BBN and
can therefore not be directly constrained. We show the instantaneous spectra
of the spin-0 particle emitted by Schwarzschild and Kerr PBHs with the public
code BlackHawk. The minimal scenario that ALPs could interact only with
photons is supposed. Then we study the stochastic oscillations between ALPs
and photons in the cosmic magnetic field.

We discuss the ALP-photon oscillations in the PMF in detail, focusing on the
distributions of oscillation probability with the homogeneous and stochastic
magnetic field scenarios. For the PMF, we consider the stochastic background
field model as the completely non-homogeneous component of the cosmic plasma.
The stringent PMF limits are taken from JS 2019, which are shown as the scale-
invariant ($n=0$, $B<47\,\rm pG$) and a violet Batchelor spectrum ($n=5$,
$B<8.9\,\rm pG$) PMFs, correspond to the magnetic field strength of the
homogeneous and stochastic magnetic field scenarios, respectively. Above
limits are derived from the CMB radiation based on the non-observation of the
magnetic field induced clumping of baryonic density fluctuations during
recombination epoch. For the homogeneous magnetic field, we show the ALP-
photon oscillation probability distributions in the ALP parameter ($m_{a}$,
$g_{a\gamma}$) space for several typical values of energy with the magnetic
field strength $B_{T}=47\,\rm pG$ (JS 2019) and $\sim 5\,\rm nG$ (Planck
2015). The oscillation probability distribution changes significantly with the
energy and magnetic field strength, and can be dramatically suppressed in the
energy above $\sim 10^{4}\,\rm GeV$. For the stochastic magnetic field, the
magnetic field strength in every domain is supposed to be same, while the
direction is random from the domain $i$ to $i+1$. We show the effective
magnetic field strength distributions with the JS 2019 PMF limit $B<8.9\,\rm
pG$, which is simulated with 100 domains for $1\,\rm Mpc$. We also show the
numerical results of oscillation probability for the stochastic magnetic
field, which is simulated with 1000 random magnetic field configurations. The
oscillation probability can be dramatically enhanced by the coupling
$g_{a\gamma}$. On the other hand, the oscillation probability can be affected
by the ALP mass $m_{a}$ in the low energy regions and becomes stable for the
high energies. Finally, we briefly introduce the effects of redshift and
photon absorption at the recombination on ALP-photon oscillations in the
Universe evolution process.

In addition, the stochastic ALP-photon oscillations induced by PBHs may have
the effects on some further phenomena, such as the CMB, the CXB, and the EGB,
which may be detected as the signals of ALP-photon oscillations in the next
generation ALP/photon measurements, $\rm e.g.$, ADMX [100, 101], IAXO [102,
103], ALPS-II [104, 105], LiteBIRD [106, 107], CMB-S4 [108, 109], THESEUS
[110, 111], eXTP [112, 113], LHAASO [114, 115], and CTA [116, 117]. The
stochastic ALP-photon oscillations in the PMF discussed in this paper may be
of great importance to these topics.

## Acknowledgments

The author would like to thank the anonymous referee for valuable comments and
suggestions improving the manuscript. We also thank Jérémy Auffinger, Junsong
Cang, and Wei Chao for helpful discussions. This work is supported by the
National Natural Science Foundation (NNSF) of China (Grants No. 11775025 and
No. 12175027).

## References

  * [1] S. Bird, I. Cholis, J. B. Muñoz, Y. Ali-Haïmoud, M. Kamionkowski, E. D. Kovetz et al., _Did LIGO detect dark matter?_ , _Phys. Rev. Lett._ 116 (2016) 201301, [1603.00464].
  * [2] S. Hawking, _Gravitationally collapsed objects of very low mass_ , _Mon. Not. Roy. Astron. Soc._ 152 (1971) 75.
  * [3] B. J. Carr and S. W. Hawking, _Black holes in the early Universe_ , _Mon. Not. Roy. Astron. Soc._ 168 (1974) 399–415.
  * [4] B. J. Carr, _The Primordial black hole mass spectrum_ , _Astrophys. J._ 201 (1975) 1–19.
  * [5] B. Carr, K. Kohri, Y. Sendouda and J. Yokoyama, _Constraints on primordial black holes_ , _Rept. Prog. Phys._ 84 (2021) 116902, [2002.12778].
  * [6] B. Carr and F. Kuhnel, _Primordial Black Holes as Dark Matter: Recent Developments_ , _Ann. Rev. Nucl. Part. Sci._ 70 (2020) 355–394, [2006.02838].
  * [7] B. Carr and F. Kuhnel, _Primordial black holes as dark matter candidates_ , _SciPost Phys. Lect. Notes_ 48 (2022) 1, [2110.02821].
  * [8] B. Carr, F. Kuhnel and M. Sandstad, _Primordial Black Holes as Dark Matter_ , _Phys. Rev. D_ 94 (2016) 083504, [1607.06077].
  * [9] G. F. Chapline, _Cosmological effects of primordial black holes_ , _Nature_ 253 (1975) 251–252.
  * [10] B. J. Carr, K. Kohri, Y. Sendouda and J. Yokoyama, _New cosmological constraints on primordial black holes_ , _Phys. Rev. D_ 81 (2010) 104019, [0912.5297].
  * [11] A. W. Strong, I. V. Moskalenko and O. Reimer, _A new determination of the extragalactic diffuse gamma-ray background from egret data_ , _Astrophys. J._ 613 (2004) 956–961, [astro-ph/0405441].
  * [12] A. Barnacka, J. F. Glicenstein and R. Moderski, _New constraints on primordial black holes abundance from femtolensing of gamma-ray bursts_ , _Phys. Rev. D_ 86 (2012) 043001, [1204.2056].
  * [13] H. Niikura, M. Takada, S. Yokoyama, T. Sumi and S. Masaki, _Constraints on Earth-mass primordial black holes from OGLE 5-year microlensing events_ , _Phys. Rev. D_ 99 (2019) 083503, [1901.07120].
  * [14] B. J. Carr and M. Sakellariadou, _Dynamical constraints on dark compact objects_ , _Astrophys. J._ 516 (1999) 195–220.
  * [15] S. M. Koushiappas and A. Loeb, _Dynamics of Dwarf Galaxies Disfavor Stellar-Mass Black Holes as Dark Matter_ , _Phys. Rev. Lett._ 119 (2017) 041102, [1704.01668].
  * [16] B. Carr and J. Silk, _Primordial Black Holes as Generators of Cosmic Structures_ , _Mon. Not. Roy. Astron. Soc._ 478 (2018) 3756–3775, [1801.00672].
  * [17] N. Afshordi, P. McDonald and D. N. Spergel, _Primordial black holes as dark matter: The Power spectrum and evaporation of early structures_ , _Astrophys. J. Lett._ 594 (2003) L71–L74, [astro-ph/0302035].
  * [18] M. Raidal, V. Vaskonen and H. Veermäe, _Gravitational Waves from Primordial Black Hole Mergers_ , _JCAP_ 09 (2017) 037, [1707.01480].
  * [19] D. Hooper, G. Krnjaic, J. March-Russell, S. D. McDermott and R. Petrossian-Byrne, _Hot Gravitons and Gravitational Waves From Kerr Black Holes in the Early Universe_ , 2004.00618.
  * [20] K. Kohri and J. Yokoyama, _Primordial black holes and primordial nucleosynthesis. 1. Effects of hadron injection from low mass holes_ , _Phys. Rev. D_ 61 (2000) 023501, [astro-ph/9908160].
  * [21] M. Ricotti, J. P. Ostriker and K. J. Mack, _Effect of Primordial Black Holes on the Cosmic Microwave Background and Cosmological Parameter Estimates_ , _Astrophys. J._ 680 (2008) 829, [0709.0524].
  * [22] Y. Ali-Haïmoud and M. Kamionkowski, _Cosmic microwave background limits on accreting primordial black holes_ , _Phys. Rev. D_ 95 (2017) 043534, [1612.05644].
  * [23] J. Cang, Y. Gao and Y. Ma, _Prospects of Future CMB Anisotropy Probes for Primordial Black Holes_ , _JCAP_ 05 (2021) 051, [2011.12244].
  * [24] K. Maki, T. Mitsui and S. Orito, _Local flux of low-energy anti-protons from evaporating primordial black holes_ , _Phys. Rev. Lett._ 76 (1996) 3474–3477, [astro-ph/9601025].
  * [25] E. T. Linton et al., _A new search for primordial black hole evaporations using the Whipple gamma-ray telescope_ , _JCAP_ 01 (2006) 013.
  * [26] D. N. Page and S. W. Hawking, _Gamma rays from primordial black holes_ , _Astrophys. J._ 206 (1976) 1–7.
  * [27] B. J. Carr, K. Kohri, Y. Sendouda and J. Yokoyama, _Constraints on primordial black holes from the Galactic gamma-ray background_ , _Phys. Rev. D_ 94 (2016) 044029, [1604.05349].
  * [28] J. H. MacGibbon and B. J. Carr, _Cosmic rays from primordial black holes_ , _Astrophys. J._ 371 (1991) 447–469.
  * [29] C. Bambi, A. D. Dolgov and A. A. Petrov, _Primordial black holes and the observed Galactic 511-keV line_ , _Phys. Lett. B_ 670 (2008) 174–178, [0801.2786]. [Erratum: Phys.Lett.B 681, 504–504 (2009)].
  * [30] C. Lunardini and Y. F. Perez-Gonzalez, _Dirac and Majorana neutrino signatures of primordial black holes_ , _JCAP_ 08 (2020) 014, [1910.07864].
  * [31] J. Cang, Y. Gao and Y.-Z. Ma, _21-cm constraints on spinning primordial black holes_ , _JCAP_ 03 (2022) 012, [2108.13256].
  * [32] J. A. Dror, H. Murayama and N. L. Rodd, _Cosmic axion background_ , _Phys. Rev. D_ 103 (2021) 115004, [2101.09287].
  * [33] J. Preskill, M. B. Wise and F. Wilczek, _Cosmology of the Invisible Axion_ , _Phys. Lett. B_ 120 (1983) 127–132.
  * [34] P. Sikivie, _Dark matter axions_ , _Int. J. Mod. Phys. A_ 25 (2010) 554–563, [0909.0949].
  * [35] D. J. E. Marsh, _Axion Cosmology_ , _Phys. Rept._ 643 (2016) 1–79, [1510.07633].
  * [36] N. Bernal, F. Hajkarim and Y. Xu, _Axion Dark Matter in the Time of Primordial Black Holes_ , _Phys. Rev. D_ 104 (2021) 075007, [2107.13575].
  * [37] N. Bernal, Y. F. Perez-Gonzalez, Y. Xu and O. Zapata, _ALP dark matter in a primordial black hole dominated universe_ , _Phys. Rev. D_ 104 (2021) 123536, [2110.04312].
  * [38] G. Choi and E. D. Schiappacasse, _PBH assisted search for QCD axion dark matter_ , _JCAP_ 09 (2022) 072, [2205.02255].
  * [39] R. Peccei and H. R. Quinn, _Constraints Imposed by CP Conservation in the Presence of Instantons_ , _Phys. Rev. D_ 16 (1977) 1791–1797.
  * [40] R. Peccei and H. R. Quinn, _CP Conservation in the Presence of Instantons_ , _Phys. Rev. Lett._ 38 (1977) 1440–1443.
  * [41] S. Weinberg, _A New Light Boson?_ , _Phys. Rev. Lett._ 40 (1978) 223–226.
  * [42] F. Wilczek, _Problem of Strong $P$ and $T$ Invariance in the Presence of Instantons_, _Phys. Rev. Lett._ 40 (1978) 279–282.
  * [43] P. Svrcek and E. Witten, _Axions In String Theory_ , _JHEP_ 06 (2006) 051, [hep-th/0605206].
  * [44] A. Arvanitaki, S. Dimopoulos, S. Dubovsky, N. Kaloper and J. March-Russell, _String Axiverse_ , _Phys. Rev. D_ 81 (2010) 123530, [0905.4720].
  * [45] L. Di Luzio, M. Giannotti, E. Nardi and L. Visinelli, _The landscape of QCD axion models_ , _Phys. Rept._ 870 (2020) 1–117, [2003.01100].
  * [46] K. Choi, S. H. Im and C. Sub Shin, _Recent Progress in the Physics of Axions and Axion-Like Particles_ , _Ann. Rev. Nucl. Part. Sci._ 71 (2021) 225–252, [2012.05029].
  * [47] G. Galanti and M. Roncadelli, _Axion-like Particles Implications for High-Energy Astrophysics_ , _Universe_ 8 (2022) 253, [2205.00940].
  * [48] D. Hooper and P. D. Serpico, _Detecting Axion-Like Particles With Gamma Ray Telescopes_ , _Phys. Rev. Lett._ 99 (2007) 231102, [0706.3203].
  * [49] H.-J. Li, J.-G. Guo, X.-J. Bi, S.-J. Lin and P.-F. Yin, _Limits on axion-like particles from Mrk 421 with 4.5-year period observations by ARGO-YBJ and Fermi-LAT_ , _Phys. Rev. D_ 103 (2021) 083003, [2008.09464].
  * [50] H.-J. Li, _Relevance of VHE blazar spectra models with axion-like particles_ , _JCAP_ 02 (2022) 025, [2112.14145].
  * [51] H.-J. Li, _Probing photon-ALP oscillations from the flat spectrum radio quasar 4C+21.35_ , _Phys. Lett. B_ 829 (2022) 137047, [2203.08573].
  * [52] M. Giannotti, I. Irastorza, J. Redondo and A. Ringwald, _Cool WISPs for stellar cooling excesses_ , _JCAP_ 05 (2016) 057, [1512.08108].
  * [53] M. Giannotti, I. G. Irastorza, J. Redondo, A. Ringwald and K. Saikawa, _Stellar Recipes for Axion Hunters_ , _JCAP_ 10 (2017) 010, [1708.02111].
  * [54] J. P. Conlon and M. C. D. Marsh, _Excess Astrophysical Photons from a 0.1–1 keV Cosmic Axion Background_ , _Phys. Rev. Lett._ 111 (2013) 151301, [1305.3603].
  * [55] S. Angus, J. P. Conlon, M. C. D. Marsh, A. J. Powell and L. T. Witkowski, _Soft X-ray Excess in the Coma Cluster from a Cosmic Axion Background_ , _JCAP_ 09 (2014) 026, [1312.3947].
  * [56] XENON collaboration, E. Aprile et al., _Excess electronic recoil events in XENON1T_ , _Phys. Rev. D_ 102 (2020) 072004, [2006.09721].
  * [57] L. Di Luzio, M. Fedele, M. Giannotti, F. Mescia and E. Nardi, _Solar axions cannot explain the XENON1T excess_ , _Phys. Rev. Lett._ 125 (2020) 131804, [2006.12487].
  * [58] C. O’HARE, _cajohare/axionlimits: Axionlimits_ , July, 2020.  10.5281/zenodo.3932430.
  * [59] K. Subramanian, _The origin, evolution and signatures of primordial magnetic fields_ , _Rept. Prog. Phys._ 79 (2016) 076901, [1504.02311].
  * [60] M. S. Turner and L. M. Widrow, _Inflation Produced, Large Scale Magnetic Fields_ , _Phys. Rev. D_ 37 (1988) 2743.
  * [61] B. Ratra, _Cosmological ’seed’ magnetic field from inflation_ , _Astrophys. J. Lett._ 391 (1992) L1–L4.
  * [62] T. Vachaspati, _Magnetic fields from cosmological phase transitions_ , _Phys. Lett. B_ 265 (1991) 258–261.
  * [63] D. Grasso and A. Riotto, _On the nature of the magnetic fields generated during the electroweak phase transition_ , _Phys. Lett. B_ 418 (1998) 258–265, [hep-ph/9707265].
  * [64] A. Brandenburg, T. Kahniashvili, S. Mandal, A. Roper Pol, A. G. Tevzadze and T. Vachaspati, _Evolution of hydromagnetic turbulence from the electroweak phase transition_ , _Phys. Rev. D_ 96 (2017) 123528, [1711.03804].
  * [65] D. Grasso and H. R. Rubinstein, _Limits on possible magnetic fields at nucleosynthesis time_ , _Astropart. Phys._ 3 (1995) 95–102, [astro-ph/9409010].
  * [66] T. Kahniashvili, A. G. Tevzadze, S. K. Sethi, K. Pandey and B. Ratra, _Primordial magnetic field limits from cosmological data_ , _Phys. Rev. D_ 82 (2010) 083005, [1009.2094].
  * [67] K. Jedamzik, V. Katalinic and A. V. Olinto, _A Limit on primordial small scale magnetic fields from CMB distortions_ , _Phys. Rev. Lett._ 85 (2000) 700–703, [astro-ph/9911100].
  * [68] K. E. Kunze and E. Komatsu, _Constraining primordial magnetic fields with distortions of the black-body spectrum of the cosmic microwave background: pre- and post-decoupling contributions_ , _JCAP_ 01 (2014) 009, [1309.7994].
  * [69] E.-j. Kim, A. Olinto and R. Rosner, _Generation of density perturbations by primordial magnetic fields_ , _Astrophys. J._ 468 (1996) 28, [astro-ph/9412070].
  * [70] E. Battaner, E. Florido and J. Jimenez-Vicente, _Magnetic fields and large scale structure in a radiation dominated universe_ , _Astron. Astrophys._ 326 (1997) 13–22, [astro-ph/9602097].
  * [71] C. Caprini and R. Durrer, _Gravitational wave production: A Strong constraint on primordial magnetic fields_ , _Phys. Rev. D_ 65 (2001) 023517, [astro-ph/0106244].
  * [72] C. Caprini and R. Durrer, _Gravitational waves from stochastic relativistic sources: Primordial turbulence and magnetic fields_ , _Phys. Rev. D_ 74 (2006) 063521, [astro-ph/0603476].
  * [73] A. Roper Pol, S. Mandal, A. Brandenburg, T. Kahniashvili and A. Kosowsky, _Numerical simulations of gravitational waves from early-universe turbulence_ , _Phys. Rev. D_ 102 (2020) 083512, [1903.08585].
  * [74] D. G. Yamazaki, T. Kajino, G. J. Mathew and K. Ichiki, _The Search for a Primordial Magnetic Field_ , _Phys. Rept._ 517 (2012) 141–167, [1204.3669].
  * [75] R. Durrer and A. Neronov, _Cosmological Magnetic Fields: Their Generation, Evolution and Observation_ , _Astron. Astrophys. Rev._ 21 (2013) 62, [1303.7121].
  * [76] Planck collaboration, P. A. R. Ade et al., _Planck 2015 results. XIX. Constraints on primordial magnetic fields_ , _Astron. Astrophys._ 594 (2016) A19, [1502.01594].
  * [77] K. Jedamzik and A. Saveliev, _Stringent Limit on Primordial Magnetic Fields from the Cosmic Microwave Background Radiation_ , _Phys. Rev. Lett._ 123 (2019) 021301, [1804.06115].
  * [78] F. Schiavone, D. Montanino, A. Mirizzi and F. Capozzi, _Axion-like particles from primordial black holes shining through the Universe_ , _JCAP_ 08 (2021) 063, [2107.03420].
  * [79] C. Evoli, M. Leo, A. Mirizzi and D. Montanino, _Reionization during the dark ages from a cosmic axion background_ , _JCAP_ 05 (2016) 006, [1602.08433].
  * [80] T. Papanikolaou, V. Vennin and D. Langlois, _Gravitational waves from a universe filled with primordial black holes_ , _JCAP_ 03 (2021) 053, [2010.11573].
  * [81] S. W. Hawking, _Black hole explosions_ , _Nature_ 248 (1974) 30–31.
  * [82] Planck collaboration, Y. Akrami et al., _Planck 2018 results. X. Constraints on inflation_ , _Astron. Astrophys._ 641 (2020) A10, [1807.06211].
  * [83] S. W. Hawking, _Particle Creation by Black Holes_ , _Commun. Math. Phys._ 43 (1975) 199–220. [Erratum: Commun.Math.Phys. 46, 206 (1976)].
  * [84] D. N. Page, _Particle Emission Rates from a Black Hole. 2. Massless Particles from a Rotating Hole_ , _Phys. Rev. D_ 14 (1976) 3260–3273.
  * [85] J. H. MacGibbon and B. R. Webber, _Quark and gluon jet emission from primordial black holes: The instantaneous spectra_ , _Phys. Rev. D_ 41 (1990) 3052–3079.
  * [86] A. Arbey and J. Auffinger, _BlackHawk: A public code for calculating the Hawking evaporation spectra of any black hole distribution_ , _Eur. Phys. J. C_ 79 (2019) 693, [1905.04268].
  * [87] A. Arbey and J. Auffinger, _Physics Beyond the Standard Model with BlackHawk v2.0_ , _Eur. Phys. J. C_ 81 (2021) 910, [2108.02737].
  * [88] D. N. Page and K. S. Thorne, _Disk-Accretion onto a Black Hole. Time-Averaged Structure of Accretion Disk_ , _Astrophys. J._ 191 (1974) 499–506.
  * [89] A. Buonanno, L. E. Kidder and L. Lehner, _Estimating the final spin of a binary black hole coalescence_ , _Phys. Rev. D_ 77 (2008) 026004, [0709.3839].
  * [90] D. N. Page, _Particle Emission Rates from a Black Hole: Massless Particles from an Uncharged, Nonrotating Hole_ , _Phys. Rev. D_ 13 (1976) 198–206.
  * [91] J. D. Barrow, P. G. Ferreira and J. Silk, _Constraints on a primordial magnetic field_ , _Phys. Rev. Lett._ 78 (1997) 3610–3613, [astro-ph/9701063].
  * [92] T. Kahniashvili, G. Lavrelashvili and B. Ratra, _CMB Temperature Anisotropy from Broken Spatial Isotropy due to an Homogeneous Cosmological Magnetic Field_ , _Phys. Rev. D_ 78 (2008) 063012, [0807.4239].
  * [93] A. Lewis, _CMB anisotropies from primordial inhomogeneous magnetic fields_ , _Phys. Rev. D_ 70 (2004) 043011, [astro-ph/0406096].
  * [94] D. Paoletti, F. Finelli and F. Paci, _The full contribution of a stochastic background of magnetic fields to CMB anisotropies_ , _Mon. Not. Roy. Astron. Soc._ 396 (2009) 523–534, [0811.0230].
  * [95] A. Mack, T. Kahniashvili and A. Kosowsky, _Microwave background signatures of a primordial stochastic magnetic field_ , _Phys. Rev. D_ 65 (2002) 123004, [astro-ph/0105504].
  * [96] G. Raffelt and L. Stodolsky, _Mixing of the Photon with Low Mass Particles_ , _Phys. Rev. D_ 37 (1988) 1237.
  * [97] A. De Angelis, G. Galanti and M. Roncadelli, _Relevance of axion-like particles for very-high-energy astrophysics_ , _Phys. Rev. D_ 84 (2011) 105030, [1106.1132]. [Erratum: Phys.Rev.D 87, 109903 (2013)].
  * [98] J. S. Schwinger, _On gauge invariance and vacuum polarization_ , _Phys. Rev._ 82 (1951) 664–679.
  * [99] A. Dobrynina, A. Kartavtsev and G. Raffelt, _Photon-photon dispersion of TeV gamma rays and its role for photon-ALP conversion_ , _Phys. Rev. D_ 91 (2015) 083003, [1412.4777]. [Erratum: Phys.Rev.D 95, 109905 (2017)].
  * [100] ADMX collaboration, N. Du et al., _A Search for Invisible Axion Dark Matter with the Axion Dark Matter Experiment_ , _Phys. Rev. Lett._ 120 (2018) 151301, [1804.05750].
  * [101] ADMX collaboration, T. Braine et al., _Extended Search for the Invisible Axion with the Axion Dark Matter Experiment_ , _Phys. Rev. Lett._ 124 (2020) 101303, [1910.08638].
  * [102] E. Armengaud et al., _Conceptual Design of the International Axion Observatory (IAXO)_ , _JINST_ 9 (2014) T05002, [1401.3233].
  * [103] IAXO collaboration, E. Armengaud et al., _Physics potential of the International Axion Observatory (IAXO)_ , _JCAP_ 06 (2019) 047, [1904.09155].
  * [104] R. Bähre et al., _Any light particle search II —Technical Design Report_ , _JINST_ 8 (2013) T09001, [1302.5647].
  * [105] ALPS collaboration, A. Spector, _ALPS II technical overview and status report_ , in _12th Patras Workshop on Axions, WIMPs and WISPs_ , pp. 133–136, 2017.  1611.05863.  DOI.
  * [106] T. Matsumura et al., _Mission design of LiteBIRD_ , _J. Low Temp. Phys._ 176 (2014) 733, [1311.2847].
  * [107] M. Hazumi et al., _LiteBIRD: A Satellite for the Studies of B-Mode Polarization and Inflation from Cosmic Background Radiation Detection_ , _J. Low Temp. Phys._ 194 (2019) 443–452.
  * [108] CMB-S4 collaboration, K. N. Abazajian et al., _CMB-S4 Science Book, First Edition_ , 1610.02743.
  * [109] K. Abazajian et al., _CMB-S4 Science Case, Reference Design, and Project Plan_ , 1907.04473.
  * [110] THESEUS collaboration, L. Amati et al., _The THESEUS space mission concept: science case, design and expected performances_ , _Adv. Space Res._ 62 (2018) 191–244, [1710.04638].
  * [111] THESEUS collaboration, G. Stratta et al., _THESEUS: a key space mission concept for Multi-Messenger Astrophysics_ , _Adv. Space Res._ 62 (2018) 662–682, [1712.08153].
  * [112] eXTP collaboration, S. N. Zhang et al., _eXTP – enhanced X-ray Timing and Polarimetry Mission_ , _Proc. SPIE Int. Soc. Opt. Eng._ 9905 (2016) 99051Q, [1607.08823].
  * [113] eXTP collaboration, S.-N. Zhang et al., _The enhanced X-ray Timing and Polarimetry mission—eXTP_ , _Sci. China Phys. Mech. Astron._ 62 (2019) 29502, [1812.04020].
  * [114] LHAASO collaboration, Z. Cao, _A future project at Tibet: The large high altitude air shower observatory (LHAASO)_ , _Chin. Phys. C_ 34 (2010) 249–252.
  * [115] LHAASO collaboration, A. Addazi et al., _The Large High Altitude Air Shower Observatory (LHAASO) Science Book (2021 Edition)_ , _Chin. Phys. C_ 46 (2022) 035001–035007, [1905.02773].
  * [116] CTA Consortium collaboration, M. Actis et al., _Design concepts for the Cherenkov Telescope Array CTA: An advanced facility for ground-based high-energy gamma-ray astronomy_ , _Exper. Astron._ 32 (2011) 193–316, [1008.3703].
  * [117] CTA Consortium collaboration, B. S. Acharya et al., _Introducing the CTA concept_ , _Astropart. Phys._ 43 (2013) 3–18.