# Evolution of axions in the presence of primordial magnetic fields

Maxim Dvornikov<EMAIL_ADDRESS>V. B. Semikoz<EMAIL_ADDRESS>Pushkov
Institute of Terrestrial Magnetism, Ionosphere and Radiowave Propagation
(IZMIRAN), 108840 Moscow, Troitsk, Russia

###### Abstract

We study the evolution of axions interacting with primordial magnetic fields
(PMFs) starting just from the QCD phase transition in the expanding universe.
This interaction is owing to the Primakoff effect. Adopting the zero mode
approximation for axions, we derive the system of equations for axions and
magnetic fields, where the expansion of the universe and the spectra of
magnetic fields are accounted for exactly. We find that the contribution of
the Primakoff effect to the dynamics of axions and magnetic fields is rather
weak. It confirms some previous estimates leading to analogous conclusions,
when accounting here for the Hubble expansion both for an uniform axion field
and non-uniform PMFs using Fourier spectra for their energy and helicity
densities. We solve the corresponding system of the evolution equations and
find that the axion zero mode, when evolving during radiation era, has its
amplitude at the level sufficient for that axion to be a good candidate for
the cold dark matter.

##  I Introduction

As the solution to the CP problem in quantum chromodynamics (QCD), Peccei and
Quinn suggested a mechanism that naturally gives rise to a Nambu-Goldstone
boson, the so-called axion PQ ; Weinberg:1977ma ; Wilczek . In addition, the
axion also provides a possible candidate for the cold dark matter (CDM) of the
Universe Rin12 . There are active searches for axions Gra15 . A progress in
the experimental studies of these particles is reported in Ref. Apr20
recently.

The simultaneous presence of primordial magnetic fields (PMFs) and an axion
field in a hot universe plasma near the QCD phase transition (QCDPT) allows to
study their mutual influence during the radiation era. We adopt that in the
universe cooling $T_{\mathrm{PQ}}\to T_{\mathrm{QCD}}$, where
$T_{\mathrm{PQ}}\simeq f_{a}=(10^{10}\div 10^{12})\,{\rm GeV}$ is the
temperature at the Peccei-Quinn phase transition,
$T_{\mathrm{QCD}}\simeq(100\div 200)\,{\rm MeV}$ is the temperature at QCDPT,
and $f_{a}$ is the Peccei-Quinn parameter, a small axion mass $m_{a}(T)$
increases and gets after mixing with $\pi^{0}$ mesons its “cold” (fixed) value
$m_{a}(0)=m_{a}$ just near the initial time in our problem
$t_{0}=t_{\mathrm{QCD}}$,

$m_{a}=6\times 10^{-6}\,{\rm eV}\left(\frac{10^{12}\,{\rm
GeV}}{f_{a}}\right).$ (1)

In its turn, PMFs originate from the hypermagnetic fields $B_{\mathrm{Y}}$
which exist before the electroweak phase transition (EWPT),
$T>T_{\mathrm{EW}}\simeq(100\div 160)\,{\rm GeV}$. Such PMF, being a seed for
observable galactic magnetic fields $B\sim 10^{-6}\,{\rm G}$, can be strong
enough at $T_{\mathrm{QCD}}$. Returning the present bound on the extragalactic
PMF, $B\leq B_{\mathrm{now}}=10^{-9}\,{\rm G}$, back to early universe times,
one gets the PMF strength
$B_{\mathrm{QCD}}=B_{\mathrm{now}}(1+z_{\mathrm{QCD}})^{2}\simeq 1.6\times
10^{14}\,{\rm G}$ for the red shift $z_{\mathrm{QCD}}=4\times 10^{11}$ at
$T_{\mathrm{QCD}}\sim 100\,{\rm MeV}$, or $B_{\mathrm{EW}}=1.6\times
10^{20}\,{\rm G}$ at $T_{\mathrm{EW}}\sim 100\,{\rm GeV}$.

On one hand, such a strong field can not influence the Universe expansion
since the PMF energy density $\rho_{\mathrm{B}}=B^{2}/2$ is much less than the
ultrarelativistic matter energy density $\rho=g^{*}\pi^{2}T^{4}/30$,
$\rho_{\mathrm{B}}\ll\rho$. Here $g^{*}$ is the effective number of degrees of
freedom (Gorbunov:2011zz, , p. 95). On the other hand, it is intriguing to
check whether such a field can be enlarged by the axion instability term
$\sim\dot{\varphi}\nabla\times{\bf B}$ entering the induction equation
Long:2015cza ,

$\frac{\partial{\bf B}}{\partial t}=\frac{\nabla^{2}{\bf
B}}{\sigma_{\mathrm{cond}}}+\frac{g_{a\gamma}}{\sigma_{\mathrm{cond}}}\frac{{\rm
d}\varphi}{{\rm d}t}\nabla\times{\bf B},$ (2)

or, vice versa, the PMF can influence evolution of that uniform axion field,
$\nabla\varphi=0$,

$\ddot{\varphi}+m_{a}^{2}\varphi=\frac{g_{a\gamma}}{\sigma_{\mathrm{cond}}}(\nabla\times{\bf
B})\cdot{\bf B}.$ (3)

Here $g_{a\gamma}\approx\alpha_{\mathrm{em}}/2\pi f_{a}$ is the coupling
constant for axion interaction with electromagnetic fields,
$\alpha_{\mathrm{em}}\approx(137)^{-1}$ is the fine structure constant,
$\sigma_{\mathrm{cond}}\simeq 100T$ is the electric conductivity in the hot
universe plasma.

The problem of the mutual influence of axions and PMFs was discussed in Ref.
Long:2015cza . However, in Ref. Long:2015cza , the Hubble expansion was
neglected for both a uniform axion gas and primordial magnetic fields. We
study the similar problem accounting for both that expansion when using
conformal variables and an inhomogeneity of PMFs through their Kolmogorov
spectra for the densities of the magnetic energy and the magnetic helicity.

Our work is organized as follows. In Sec. II, we present the complete system
of the evolution equations for the homogeneous axion field $\varphi(\eta)$ and
the spectra of the conformal PMF energy density $\rho_{c}(k_{c},\eta)$ and the
conformal PMF helicity density $h_{c}(k_{c},\eta)$. Here $\eta$ is the
conformal time and $k_{c}$ is the conformal momentum. Details for derivation
of such a system are given in Appendix A. In Sec. II.2, we formulate the
initial conditions for this system of evolution equations. In Sec. III, we
present the results of the numerical calculations illustrated by plots both
for the maximum extragalactic field $B_{\mathrm{now}}=10^{-9}\,{\rm G}$ and
without PMFs at all. In Sec. IV, we discuss our results.

##  II Axions in the early universe plasma with PMFs

For a homogeneous axion field $\nabla\varphi=0$, or for the zero axion mode f1
, the axion evolution equation in the presence of PMF and in the Friedmann-
Robertson-Walker (FRW) metric,
$\mathrm{d}s^{2}=a^{2}(t)(\mathrm{d}\eta^{2}-\mathrm{d}{\bf x}^{2})$, where
$a$ is the scale factor and $\mathrm{d}\eta=dt/a(t)$, takes the following form
in the variables $x^{\mu}=(\eta,{\bf x})$ (see Eqs. (41) and (42) in Appendix
A):

$\frac{\mathrm{d}^{2}\varphi}{\mathrm{d}\eta^{2}}+\frac{2}{a}\frac{\mathrm{d}a}{\mathrm{d}\eta}\frac{\mathrm{d}\varphi}{\mathrm{d}\eta}+m_{a}^{2}a^{2}\varphi=-\frac{g_{a\gamma}}{2a^{2}}\frac{\mathrm{d}h_{c}(\eta)}{\mathrm{d}\eta}.$
(4)

Here $h_{c}(\eta)=a^{3}h(t)=V^{-1}\int\mathrm{d}^{3}x{\bf A}_{c}(\eta,{\bf
x})\cdot{\bf B}_{c}(\eta,{\bf x})$ is the conformal PMF helicity density,
${\bf A}_{c}(\eta,{\bf x})=a{\bf A}$ is the vector potential, ${\bf
B}_{c}(\eta,{\bf x})=\nabla\times{\bf A}_{c}(\eta,{\bf x})=a^{2}{\bf B}$ is
the comoving magnetic field in conformal variables f2 .

We can calculate the derivative $\mathrm{d}h_{c}(\eta)/\mathrm{d}\eta$ for the
magnetic helicity density $h_{c}(\eta)=\int\mathrm{d}k_{c}h(\eta,k_{c})$ given
by its isotropic Fourier spectrum $h_{c}(k,\eta)=k^{2}_{c}{\bf
A}_{c}(k_{c},\eta){\bf B}_{c}^{*}(k_{c},\eta)/(2\pi^{2}V)$ as
$\partial_{\eta}h_{c}(\eta)=\int\mathrm{d}k_{c}\partial_{\eta}h_{c}(\eta,k_{c})$,
where $k_{c}=ak=\text{const}$ is the comoving momentum that does not depend on
time. We use the evolution equations for PMF spectra for this purpose (see
Appendix A),

$\displaystyle\frac{\partial h_{c}(k_{c},\eta)}{\partial\eta}$
$\displaystyle=-\frac{2k_{c}^{2}}{\sigma_{c}}h_{c}(k_{c},\eta)+\frac{4g_{a\gamma}}{\sigma_{c}}\frac{\mathrm{d}\varphi(\eta)}{\mathrm{d}\eta}\rho_{c}(k_{c},\eta),$
(5) $\displaystyle\frac{\partial\rho_{c}(k_{c},\eta)}{\partial\eta}$
$\displaystyle=-\frac{2k_{c}^{2}}{\sigma_{c}}\rho_{c}(k_{c},\eta)+k_{c}^{2}\frac{g_{a\gamma}}{\sigma_{c}}\frac{\mathrm{d}\varphi(\eta)}{\mathrm{d}\eta}h_{c}(k_{c},\eta),$
(6)

where the magnetic energy density
$\rho_{c}(\eta)=\int\mathrm{d}k_{c}\rho_{c}(\eta,k_{c})=B^{2}_{c}(\eta)/2$ is
given by the isotropic Fourier spectrum $\rho_{c}(k_{c},\eta)=k_{c}^{2}{\bf
B}_{c}(k_{c},\eta){\bf B}_{c}^{*}(k_{c},\eta)/(4\pi^{2}V)$, $k_{c}=|{\bf
k}_{c}|$ is the absolute value of the comoving Fourier momentum, and
$\sigma_{c}=a\sigma_{\mathrm{cond}}$ is given by the electric conductivity
$\sigma_{\mathrm{cond}}\simeq 100T$ in the hot universe plasma.

Substituting the derivative
$\mathrm{d}h_{c}(\eta)/\mathrm{d}\eta=\int_{k_{\mathrm{min}}}^{k_{\mathrm{max}}}\mathrm{d}k_{c}\partial_{\eta}h_{c}(\eta,k_{c})$,
with the integrand given by Eq. (5), we can rewrite the axion evolution Eq.
(4) as

$\displaystyle\frac{\mathrm{d}^{2}\varphi(\eta)}{\mathrm{d}\eta^{2}}+$
$\displaystyle\frac{2}{a(\eta)}\frac{\mathrm{d}\varphi(\eta)}{\mathrm{d}\eta}\left(\frac{\mathrm{d}a(\eta)}{\mathrm{d}\eta}+\frac{g_{a\gamma}^{2}}{a(\eta)\sigma_{c}}\int_{k_{\mathrm{min}}}^{k_{\mathrm{max}}}\rho_{c}(k_{c},\eta)\mathrm{d}k_{c}\right)$
$\displaystyle+m_{a}^{2}a^{2}(\eta)\varphi(\eta)=\frac{g_{a\gamma}}{a^{2}(\eta)\sigma_{c}}\int_{k_{\mathrm{min}}}^{k_{\mathrm{max}}}k_{c}^{2}h_{c}(k_{c},\eta)\mathrm{d}k_{c},$
(7)

that is the generalization of the axion wave Eq. (40a) in Ref. Long:2015cza
accounting for the Hubble expansion and the Fourier spectra of non-uniform
electromagnetic fields.

The system of equations that describes an evolution of the axion field in PMF,
given by Eq. (7), is completed by the differential Eqs. (5) and (6) for the
spectra $h_{c}(k_{c},\eta)$ and $\rho_{c}(k_{c},\eta)$. Here, in the causal
scenario, the momentum limits $k_{\mathrm{min}}=k_{1}a(t)$ and
$k_{\mathrm{max}}=k_{2}a(t)$ obey the inequality
$k_{2}>k_{1}=l_{\mathrm{H}}^{-1}(T_{0})=\text{const}$ where the horizon size
$l_{\mathrm{H}}(T_{0})$ is chosen at the QCD phase transition time
$t_{0}=(T_{0}/{\rm MeV})^{-2}\,\mathrm{s}$ when the axion mass $m_{a}(T)$
reaches its maximum (“cold”) value, given in Eq. (1), during cooling of the
Universe, $T\to T_{0}$. We consider such a time $t_{0}\simeq
10^{-4}\,\mathrm{s}$ as the initial moment in our simulations.
Correspondingly, for the conformal time during radiation epoch,
$\eta(t)=\int_{t_{0}}^{t}\mathrm{d}t/a(t)+\eta(t_{0})$, for which
$a(t)=a(t_{0})\sqrt{t/t_{0}}$ is the scale factor in FRW metric, we put the
initial value $\eta(t_{0})=0$ resulting in

$\eta(t)=\frac{2t_{0}}{a(t_{0})}\left[\sqrt{\frac{t}{t_{0}}}-1\right],$ (8)

where the initial scale factor
$a(t_{0})=T_{\mathrm{now}}/T_{0}=(1+z_{\mathrm{QCD}})^{-1}$ is given by by the
red shift at the QCD phase transition, $z_{\mathrm{QCD}}\simeq 4\times
10^{11}$, since $T_{\mathrm{now}}=2.725\,\mathrm{K}$ and
$T_{0}=T_{\mathrm{QCD}}=10^{12}\,\mathrm{K}\simeq 100\,{\rm MeV}$.

We assume that PMF survives against the ohmic diffusion in the early Universe
in a wide region of Fourier momenta $k_{1}<k_{2}$. Obviously, such a field
survives for the minimal and fixed value
$k_{1}=l_{\mathrm{H}}^{-1}=(2/3)\times 10^{-20}\,{\rm GeV}$ given by the
horizon size $l_{\mathrm{H}}$, or maximum PMF scale
$\Lambda^{(\mathrm{max})}_{\mathrm{B}}=l_{\mathrm{H}}$ at the initial time
$t_{0}$. Additionally, we demand that the expansion time
$t=t_{0}(T_{0}/T)^{2}$ be faster than the diffusion time
$t_{\mathrm{diff}}=\sigma_{\mathrm{cond}}(\Lambda^{(\mathrm{min})}_{\mathrm{B}})^{2}$,
or $t<t_{\mathrm{diff}}$ even for the minimal PMF scale
$\Lambda^{(\mathrm{min})}_{\mathrm{B}}=k_{2}^{-1}$. For the initial time
$t_{0}$ this means

$t_{0}\leq\frac{\sigma_{c}}{a(t_{0})k_{2}^{2}},~{}~{}~{}~{}T\leq T_{0},$ (9)

where the conductivity $\sigma_{\mathrm{cond}}=\sigma_{c}/a(t_{0})$ was
substituted.

The system of the evolution Eqs. (6)-(7), rewritten for the three
dimensionless functions
$\mathcal{H}(\kappa,\tau)=g_{a\gamma}^{2}h_{c}(k_{c},\eta)/2$,
$\mathcal{R}(\kappa,\tau)=(g_{a\gamma}^{2}/k_{\mathrm{max}})\rho_{c}(k_{c},\eta)$,
$\Phi(\tau)=2(k_{\mathrm{max}}g_{a\gamma}/\sigma_{c})\varphi(\eta)$ with
change of the conformal time $\eta$ to the dimensionless one,
$\tau=2k_{\mathrm{max}}^{2}\eta/\sigma_{c}$, reads,

$\displaystyle\frac{{\rm d}\mathcal{H}}{{\rm d}\tau}=$
$\displaystyle-\kappa^{2}\mathcal{H}(\kappa,\tau)+\frac{{\rm
d}\Phi(\tau)}{{\rm d}\tau}\mathcal{R}(\kappa,\tau)$ (10)
$\displaystyle\frac{{\rm d}\mathcal{R}}{{\rm d}\tau}=$
$\displaystyle-\kappa^{2}\mathcal{R}(\kappa,\tau)+\kappa^{2}\frac{{\rm
d}\Phi(\tau)}{{\rm d}\tau}\mathcal{H}(\kappa,\tau)$ (11)
$\displaystyle\frac{{\rm d}^{2}\Phi}{{\rm d}\tau^{2}}$
$\displaystyle+\frac{2}{a(\tau)}\frac{{\rm d}\Phi}{{\rm
d}\tau}\left(\frac{{\rm d}a(\tau)}{{\rm
d}\tau}+\frac{1}{2a(\tau)}\int_{\kappa_{m}}^{1}\mathcal{R}(\kappa,\tau)d\kappa\right)$
$\displaystyle+a^{2}\mu^{2}\Phi(\tau)=\frac{1}{a^{2}(\tau)}\int_{\kappa_{m}}^{1}\kappa^{2}\mathcal{H}(\kappa,\tau)d\kappa.$
(12)

Here, in Eq. (12), $\mu=m_{a}\sigma_{c}/2k_{\mathrm{max}}^{2}$ is the
dimensionless parameter given by the axion mass in Eq. (1). The lower limit
$\kappa_{m}=k_{\mathrm{min}}/k_{\mathrm{max}}=k_{1}/k_{2}$ for the
dimensionless Fourier momentum $\kappa=k_{c}/k_{\mathrm{max}}=k/k_{2}$ obeys
the constraints

$2.6\times 10^{-11}\leq\kappa_{m}\ll 1,$ (13)

where we account for Eq. (9).

Note that we have three free parameters for PMF in our problem:

  * •

a wide width of the Fourier spectrum for PMF given by a small varying
parameter $\kappa_{m}$ in Eq. (13), or by varying
$k_{\mathrm{max}}=k_{\mathrm{min}}/\kappa_{m}$ for the fixed
$k_{\mathrm{min}}=a(t_{0})k_{1}$;

  * •

the helicity parameter $q$, $0\leq|q|\leq 1$, entering the initial magnetic
helicity spectrum (see Eq. (25) below), where $q=0$ corresponds to the non-
helical initial PMF, whereas $|q|=1$ to the maximum helical PMF;

  * •

an undetermined cosmological magnetic field $B_{\mathrm{now}}$ at the present
red shift $z=0$ originated in our scenario from the initial PMF
$B_{0}=B_{\mathrm{now}}(1+z_{\mathrm{QCD}})^{2}$ at $z_{\mathrm{QCD}}=4\times
10^{11}$ in the early Universe.

These three parameters fully describe necessary characteristics of PMF: the
spatial scale, its structure (the topology) and the strength.

Using the definition of the conformal time in Eq. (8), the scale in the
radiation era,
$a(t)=a(t_{0})\sqrt{t/t_{0}}=a(t_{0})[1+a(t_{0})\eta(t)/2t_{0}]$, can be
rewritten in the new variables as $a(\tau)=a(t_{0})[1+K\tau]\equiv
a(t_{0})\tilde{a}(\tau)$, with the factor

$\tilde{a}(\tau)=1+K\tau,\quad K=3.7\times 10^{20}\kappa_{m}^{2},$ (14)

where $\tau\geq 0$ f3 .

###  II.1 Axion energy density

The important characteristic of axions is the axion energy density
$\rho_{a}(t)=T_{00}$, where $T_{00}$ is the component of the energy-momentum
tensor for a scalar field,

$T_{\mu\nu}=\partial_{\mu}\varphi\partial_{\nu}\varphi-\frac{g_{\mu\nu}}{2}\left(g^{\lambda\rho}\partial_{\lambda}\varphi\partial_{\rho}\varphi-
m_{a}^{2}\varphi^{2}\right),$ (15)

In case of a spatially homogeneous field, basing on Eq. (15), one gets that

$\displaystyle\rho_{a}(t)=$
$\displaystyle\frac{1}{2}[\dot{\varphi}^{2}+m_{a}^{2}\varphi^{2}]=\frac{1}{2}\left(\frac{k_{\mathrm{max}}}{g_{a\gamma}a}\right)^{2}$
$\displaystyle\times\left[\left(\frac{{\rm d}\Phi}{{\rm
d}\tau}\right)^{2}+\mu^{2}a^{2}\Phi^{2}(\tau)\right],$ (16)

where $k_{\mathrm{max}}=k_{2}a=k_{1}a/\kappa_{m}$ is given by the parameter
$\kappa_{m}$ in Eq. (13) and the minimal momentum $k_{1}$ in our causal
scenario at $t_{0}=t_{\mathrm{QCD}}$,
$k_{1}=l_{\mathrm{H}}(t_{0})^{-1}=(2/3)\times 10^{-20}\,\text{GeV}$. This
results in the dimensional factor ahead the brackets in Eq. (II.1),

$C_{m}=\frac{1}{2}\left(\frac{k_{\mathrm{max}}}{g_{a\gamma}a}\right)^{2}=\frac{2}{9}\left(\frac{10^{-10}}{\kappa_{m}^{2}}\right)\,\text{GeV}^{4},$
(17)

that depends on the PMF spectrum width $\kappa_{m}$.

###  II.2 Initial conditions

The initial axion zero mode amplitude $\varphi(t_{0})\simeq f_{a}$,
$\bar{\theta}(t_{0})\simeq 1$, see, e.g., Ref. Long:2015cza , corresponds to

$\Phi(\tau=0)=\frac{2k_{\mathrm{max}}g_{a\gamma}}{\sigma_{c}}\varphi(t_{0})=\frac{1.5\times
10^{-24}}{\kappa_{m}},$ (18)

where $\kappa_{m}$ is the free parameter constrained in Eq. (13).

The initial derivative $\mathrm{d}\Phi/\mathrm{d}\tau|_{\tau=0}$, necessary
for solving Eq. (12), is derived using the virial theorem
$\langle\dot{\theta}^{2}\rangle=m_{a}^{2}\langle\theta^{2}\rangle$ resulting
from the axion energy density in Eq. (II.1) around its potential minimum
Duffy:2009ig ,

$\rho=\frac{f_{a}^{2}}{2}[\dot{\theta}^{2}+m_{a}^{2}(t)\theta^{2}],\quad\theta(t)=\frac{\varphi}{f_{a}}.$
(19)

We put roughly $({\rm d\varphi}/{\rm d}t)_{t=t_{0}}=m_{a}\varphi(t_{0})$
relying on that virial theorem. Then, one gets

$\left.\frac{{\rm d}\Phi}{{\rm
d}\tau}\right|_{\tau=0}=\left(\frac{g_{a\gamma}a(t_{0})}{k_{\mathrm{max}}}\right)m_{a}f_{a}=10^{3}\kappa_{m},$
(20)

where we substituted the “cold” axion mass $m_{a}$ in Eq. (1) at the moment of
the QCD phase transition, $m_{a}(T)\to m_{a}$ at $T\to T_{\mathrm{QCD}}$.

Note that the mass term, $(\mu a^{2})\Phi$, in Eq. (12)

$(\mu a)^{2}=2.85\times 10^{30}\kappa_{m}^{2}(1+K\tau)^{2}.$ (21)

is rather great. We shall see in Sec. III that this term is dominant.

The initial PMF spectra $\mathcal{R}(\kappa,\tau=0)$ and
$\mathcal{H}(\kappa,\tau=0)$,

$\displaystyle\mathcal{R}(\kappa,\tau=0)=\frac{g_{a\gamma}^{2}}{k_{\mathrm{max}}}\rho_{c}(k_{c},\eta=0),$
(22)
$\displaystyle\mathcal{H}(\kappa,\tau=0)=\frac{g_{a\gamma}^{2}}{2}h_{c}(k_{c},\eta=0),$
(23)

are given by the PMF spectra in the known form f4 ,

$\displaystyle\rho_{c}(k_{c},\eta=0)=$ $\displaystyle
Ck_{c}^{n_{\mathrm{B}}},$ $\displaystyle C=$
$\displaystyle\frac{(n_{\mathrm{B}}+1)B_{\mathrm{now}}^{2}}{2k_{\mathrm{max}}^{n_{\mathrm{B}}+1}(1-\kappa_{m}^{n_{\mathrm{B}}+1})},$
(24) $\displaystyle h_{c}(k_{c},\eta=0)=$
$\displaystyle\frac{2q}{k_{c}}\rho_{c}(k_{c},\eta=0),\quad 0\leq|q|\leq 1.$
(25)

The constant $C$ for the spectrum in Eq. (II.2) is given by normalization on
the initial magnetic energy density
$B_{c}^{2}(\eta=0)/2=\int\mathrm{d}k_{c}\rho_{c}(k_{c},\eta=0)$ for the
comoving PMF strength $B_{c}(\eta=0)=B_{\mathrm{now}}$. Thus, substituting
Eqs. (II.2) and (25) into Eqs. (22) and (23) correspondingly we obtain the
initial PMF spectra

$\displaystyle\mathcal{R}(\kappa,\tau=0)=$
$\displaystyle\frac{C_{0}(n_{\mathrm{B}}+1)}{2(1-\kappa_{m}^{n_{\mathrm{B}}+1})}\kappa^{n_{\mathrm{B}}},$
(26) $\displaystyle\mathcal{H}(\kappa,\tau=0)=$ $\displaystyle
q\frac{C_{0}(n_{\mathrm{B}}+1)}{2(1-\kappa_{m}^{n_{\mathrm{B}}+1})}\kappa^{n_{\mathrm{B}}-1},$
$\displaystyle C_{0}=$
$\displaystyle\left(\frac{g_{a\gamma}B_{\mathrm{now}}}{k_{\mathrm{max}}}\right)^{2}.$
(27)

We consider below the Kolmogorov’s spectrum with $n_{\mathrm{B}}=-5/3$.

##  III Results

In this section, we present the numerical solution of Eqs. (10)-(12) with the
initial conditions fixed in Sec. II.2.

We rely below in Fig. 1, on a maximum value for PMF at present time
$B_{\mathrm{now}}=10^{-9}\,{\rm G}$ that obeys the bounds

  * •

$B_{1\,\mathrm{Mpc}}<4.4\times 10^{-9}\,{\rm G}$ found in Ref. Plank  from CMB
anisotropies (effect on CMB polarization induced by Faraday rotation) for zero
magnetic helicity ($q=0$);

  * •

$B_{1\,\mathrm{Mpc}}<5.6\times 10^{-9}\,{\rm G}$ for a maximally helical
magnetic field ($q=1$).

Notice that, at the initial time $t_{0}=t_{\mathrm{QCD}}$ chosen in our
problem for the causal scenario, the maximum spatial PMF scale
$\Lambda_{\mathrm{B}}^{(\mathrm{max})}=l_{\mathrm{H}}(t_{0})\simeq 3\times
10^{6}\,\mathrm{cm}$ occurs too small after its following growth till present
time. It grows upto the scale $\Lambda_{\mathrm{B}}(z=0)\simeq 1\,\mathrm{pc}$
only since the horizon expands during the cooling much faster,
$l_{\mathrm{H}}\sim T^{-2}$, than the correlation length
$\Lambda_{\mathrm{B}}\sim T^{-1}$. However, the inverse cascade in the
relativistic MHD, accounting for non-linear terms in Navier-Stokes equation
(we do not touch here the full MHD approach), can rearrange the Fourier MHD
spectra in such a way that $\Lambda_{\mathrm{B}}$ as a measure of the
coherence length of the magnetic field, can increase in the coherence by the
five orders of magnitude Brandenburg:1996fc . This fact allows us to get
closer to the present scale $\Lambda_{\mathrm{B}}\sim 1\,{\rm Mpc}$ for PMF
bounded in Ref. Plank . The other problem is that, around the time of the
recombination, the photon diffusion becomes very large and so-called the Silk
mechanism could destroy the PMF characteristics. This danger was refuted in
Ref. Brandenburg:1996sa , where nonlinear effects were shown to prevent most
likely this problem from happening.

In Fig. 1, the normalized axion field
$\Phi(\tau)/\Phi(0)=\varphi(t)/\varphi(t_{0})=\bar{\theta}(t)$, starting from
$\bar{\theta}(t_{0})=\varphi(t_{0})/f_{a}=1$, grows forty times just at the
beginning (see inset in Fig. 1), and then reduces soon to the same value
$\bar{\theta}(t)\sim 1$ acceptable for a future axion contribution to CDM.
This result does not depend on the presence of PMF, see in Fig.1 where we put
$g_{a\gamma}=0$ and $q=0$ excluding interaction of axions with the PMF. In
both cases, this evolution of $\bar{\theta}(t)$-amplitudes happens during
lepton era in the hot universe plasma.

Figure 1: (a) The evolution of the axion zero mode amplitude
$\Phi(\tau)/\Phi(\tau=0)\equiv\varphi(t)/f_{a}=\bar{\theta}(t)$ for the
Peccei-Quinn parameter $f_{\mathrm{PQ}}\simeq f_{a}=10^{12}\,{\rm GeV}$ and
the Fourier spectrum width $\kappa_{m}=10^{-10}$ in Eq. (13), choosing the
maximal $B_{\mathrm{now}}=10^{-9}\,{\rm G}$. The seed magnetic field is
supposed to be nonhelical, i.e. $q=0$ in Eq. (25). (b) The evolution of the
axion zero mode amplitude
$\Phi(\tau)/\Phi(\tau=0)\equiv\varphi(t)/f_{a}=\bar{\theta}(t)$ for the
Peccei-Quinn parameter $f_{\mathrm{PQ}}\simeq f_{a}=10^{12}\,{\rm GeV}$
without PMF contribution, $g_{a\gamma}=0$ and $q=0$.

One can explain the appearance of a great amplitude $\Phi(\tau)/\Phi(0)\sim
40$ just after the start at $t\geq t_{0}=10^{-4}\,\mathrm{s}$ using initial
conditions in Eqs. (18)-(21). For a slow change of the scale factor $a(\tau)$
and neglecting the PMF terms, we get oscillations of the axion field,
$\Phi(\tau)=\mathcal{A}\cos\omega\tau+\mathcal{B}\sin\omega\tau$ obeying the
simplified equation in Eq. (12), $\Phi_{\tau\tau}+\omega^{2}\Phi(\tau)=0$,
$\omega=\mu a$, $\tau\ll 1$. Since $\mathcal{A}=\Phi(0)$ and
$\mathcal{B}=\mathrm{d}\Phi/\mathrm{d}\tau|_{\tau=0}/\omega$, one finds from
the initial conditions the corresponding ratio $\Phi(\tau)/\Phi(0)$ for the
amplitude of $\Phi(\tau)$ at $\omega\tau=\pi/2$,

$\displaystyle\frac{\Phi(\tau)}{\Phi(0)}=$
$\displaystyle\frac{\mathcal{B}}{\mathcal{A}},\quad\text{or}\quad\frac{\Phi(\tau)}{\Phi(0)}=\frac{(\mathrm{d}\Phi/\mathrm{d}\tau)_{\tau=0}}{\Phi(0)(\mu
a)}$ $\displaystyle=\frac{10^{3}\kappa_{m}}{1.5\times
10^{-24}(\sqrt{2.85}\times 10^{15})}\simeq 40,$ (28)

where we substitute the values of the parameters from Eqs. (18)-(21) and put
$\kappa_{m}=10^{-10}$.

In what follows the phase $\bar{\theta}=\varphi(t)/f_{a}=\Phi(\tau)/\Phi(0)$
oscillates and diminishes with the universe expansion as $\sim t^{-3/4}$ (see
Eq. (45) in Ref. Sikivie:2006ni ) for a zero mode evolution neglecting
electromagnetism for axions, or as $\sim 1/a(\eta)\sqrt{a(\eta)}$ during
radiation era, $a\sim\sqrt{t}$ (see Ref. (Gorbunov:2011zz, , p. 40)). Of
course, the oscillation frequency $\omega(\tau)=\mu a(\tau)$ rises with the
expansion that blackens oscillations in Fig. 1 for times $t\gg t_{0}$. Compare
the insets in that plot where oscillations with the large amplitude
$\Phi/\Phi(0)\sim 40$ are distinguishable at $t\sim t_{0}$.

The more smooth evolution curves are seen in Fig. 2 for the normalized axion
energy density $\rho_{a}(t)=C_{m}\tilde{\rho}(\tau)$ given by Eqs. (II.1) and
(17). The ratio of $\rho_{a}(t)/\rho(t_{0})$ corresponds to the ratio of
dimensionless functions $\tilde{\rho}(\tau)$,

$\tilde{\rho}(\tau)=\left[\left(\frac{{\rm d}\Phi}{{\rm
d}\tau}\right)^{2}+\mu^{2}a^{2}\Phi^{2}(\tau)\right],$ (29)

and $\tilde{\rho}(0)$ plotted in Fig. 2,

$\frac{\tilde{\rho}(\tau)}{\tilde{\rho}(0)}=\frac{({\rm d}\Phi/{\rm
d}\tau)^{2}+\mu^{2}a^{2}\Phi^{2}(\tau)}{({\rm d}\Phi/{\rm
d}\tau)_{\tau=0}^{2}+\mu^{2}a^{2}(0)\Phi^{2}(0)}.$ (30)

For $\kappa_{m}=10^{-10}$ the huge dimensional factor $C_{m}=(2/9)\times
10^{10}~{}{\rm GeV}^{4}$ in Eq. (II.1) is compensated in the product
$\rho_{a}(t_{0})=C_{m}\tilde{\rho}(\tau=0)$ by the small value
$\tilde{\rho}(\tau=0)=10^{-14}$ calculated numerically from the system of
evolution Eqs. (10)-(12), or
$\rho_{a}(t_{0})=C_{m}\tilde{\rho}(\tau=0)=(2/9)\times 10^{-4}~{}{\rm
GeV}^{4}$ is the initial axion energy density. On the other hand, the ultra-
relativistic matter energy density $\rho=\pi^{2}T^{4}g^{*}/90$ in Friedmann’s
law is given by the number of relativistic degrees of freedom $g^{*}=10.75$
and temperature $T_{\mathrm{QCD}}=0.1~{}{\rm GeV}$ at the same moment $t_{0}$,
$\rho=1.18\times 10^{-4}~{}{\rm GeV}^{4}$. Thus, we get at the start a
remarkable ratio of energy densities for the axion as a candidate to CDM,
$\rho_{a}(t_{0})/\rho=0.19$. Notice that this ratio does not depend on a width
of the PMF Fourier spectrum $\kappa_{m}$ in Eq. (13) since that parameter
shrinks in the initial axion energy density
$\rho_{a}(t_{0})=C_{m}\tilde{\rho}(0)$ where $C_{m}\sim\kappa_{m}^{-2}$ in Eq.
(17) and $\tilde{\rho}(0)\sim(d\Phi/d\tau)_{\tau=0}^{2}\sim\kappa_{m}^{2}$ in
Eq. (20).

Figure 2: (a) The evolution of the axion energy density
$\tilde{\rho}(\tau)/\tilde{\rho}(\tau=0)$ for the Peccei-Quinn parameter
$f_{\mathrm{PQ}}\simeq f_{a}=10^{12}\,{\rm GeV}$ and the Fourier spectrum
width $\kappa_{m}=10^{-10}$ in Eq. (13), choosing the maximal
$B_{\mathrm{now}}=10^{-9}\,{\rm G}$. The seed magnetic field is supposed to be
nonhelical, i.e. $q=0$ in Eq. (25). (b) The evolution of the axion energy
amplitude $\tilde{\rho}(\tau)/\tilde{\rho}(\tau=0)$ for the Peccei-Quinn
parameter $f_{\mathrm{PQ}}\simeq f_{a}=10^{12}\,{\rm GeV}$ without PMF
contribution, $g_{a\gamma}=0$ and $q=0$.

Then the axion energy density $\rho_{a}\sim\tilde{\rho}(\tau)$ reduces as seen
in Fig. 2 as well as the matter density $\rho$ decreases due to the universe
cooling, $\rho\sim T^{4}$, therefore supporting approximately the same ratio
$\rho_{a}(t)/\rho(t)\simeq 0.2$.

We are not guaranteed to observe in future the maximum PMF
$B_{\mathrm{now}}=10^{-9}\,{\rm G}$ while some predictions with a lower bound
for $B_{1\,{\rm Mpc}}(z=0)>10^{-16}\,{\rm G}$ exist Neronov:1900zz . They
arise, in particular, from $\gamma$-ray observation in the Fermi experiment
Fermi .

Recently, in Ref. Jedamzik , it was suggested to look for PMF at a level of
the comoving strength $B_{\mathrm{now}}=Ba^{2}=\text{const}\simeq
10^{-11}\,{\rm G}$. In the last photon scattering (for which CMB instruments
are sensitive), such a field corresponds to the PMF strength
$B(z_{r})=B_{\mathrm{now}}(1+z_{r})^{2}\simeq 10^{-5}\,{\rm G}$, whereas, in
our problem, to $B_{0}=B_{\mathrm{now}}(1+z_{\mathrm{QCD}})^{2}\simeq
1.6\times 10^{12}\,{\rm G}$.

##  IV Discussion

In the present work, we have studied the evolution of a zero axion mode,
$\nabla\varphi=0$, in a hot plasma of the early Universe in the presence of
PMF that do not exceed $B_{\mathrm{now}}=10^{-9}\,{\rm G}$ at the present
time. We conclude that such a mode can contribute to the axion CDM since its
amplitude tends to $\bar{\theta}(t)=\varphi(t)/f_{a}\to 1$ evolving during
radiation era, see Fig. 1. Simultaneously the corresponding axion energy
density $\rho_{a}(t)$ reduces in such a way that its part in the total matter
energy density $\rho$ tends to the CDM ratio $\rho_{a}(t)/\rho\simeq 0.2$, see
comments on that after Eq. (30) concerning the axion energy evolution in Fig.
2.

Notice that our PMF amplitudes
$\sqrt{\int_{\kappa_{m}}^{1}\mathrm{d}\kappa{\mathcal{R}}(\kappa,\tau)}\sim
Ba^{2}=\text{const}$ correspond to the comoving strength for a given
$B_{\mathrm{now}}=Ba^{2}=\text{const}$ and $\kappa_{m}$. Therefore, the ratio

$\sqrt{\frac{\int_{\kappa_{m}}^{1}\mathrm{d}\kappa{\mathcal{R}}(\kappa,\tau)}{\int_{\kappa_{m}}^{1}\mathrm{d}\kappa{\mathcal{R}}(\kappa,\tau=0)}}\approx
1,$ (31)

does not depend on the conformal time $\eta\sim\tau$, being slightly violated
by PMF interaction with axions. The similar conservation concerns the magnetic
helicity density,
$\int_{\kappa_{m}}^{1}\tilde{\mathcal{H}}(\kappa,\tau)\mathrm{d}\kappa\sim
ha^{3}=\text{const}$, calculated from the same self-consistent system of Eqs.
(10)-(12) for each given $B_{\mathrm{now}}$ and $\kappa_{m}$. Hence, the
magnetic helicity density for the case $q=0$,

$\int_{\kappa_{m}}^{1}{\mathcal{H}}(\kappa,\tau)\mathrm{d}\kappa\approx\int_{\kappa_{m}}^{1}{\mathcal{H}}(\kappa,\tau=0)\mathrm{d}\kappa=0,$
(32)

does not depend on time as well. The similar independence on time occurs in
the case of a maximally helical PMF with $q=1$. For the PMF amplitude one
obtains precisely the result in Eq. (31), while for the magnetic helicity
density in Eq. (32) changes to

$\int_{\kappa_{m}}^{1}{\mathcal{H}}(\kappa,\tau)\mathrm{d}\kappa\approx\int_{\kappa_{m}}^{1}{\mathcal{H}}(\kappa,\tau=0)\mathrm{d}\kappa=\text{const},$
(33)

where the constant depends on the initial conditions.

The conservation of the magnetic helicity density in Eq. (33) is not
surprising in the absence of the axion interaction with PMF’s through
Primakoff’s effect. Indeed, the comoving mean energy density
$\varepsilon_{\mathrm{M}}=\langle{\bf
B}^{2}\rangle/2=\smallint\mathrm{d}k_{c}\rho_{c}(k_{c},\eta)$ decays with time
such that $\varepsilon_{\mathrm{M}}\propto\eta^{-2/3}$ while the correlation
length grows like $\xi_{\mathrm{M}}\propto\eta^{2/3}$ as it follows from the
dimensionless analysis based on the Kolmogorov-type approach for the magnetic
turbulence, see, e.g., Refs. Kahniashvili:2012uj ; Brandenburg:2017neh . This
results in conservation of the magnetic helicity density,
$h_{c}(\eta)=\smallint\mathrm{d}k_{c}h_{c}(k_{c},\eta)=\langle{\bf
B}^{2}\rangle\xi_{\mathrm{M}}=\text{const}$. If we switch additionally on the
interaction of axions with PMF’s, the violation of that conservation law is
negligible because of a very small coupling constant
$g_{a\gamma}=\alpha_{\mathrm{em}}/2\pi f_{a}$; cf. Eqs (32) and (33).

Equations (31)-(33) were directly checked by our numerical simulations. Thus,
we confirm the conclusion of Ref. Long:2015cza  that the axion zero mode
($\nabla\varphi=0$) practically does not influence any PMF characteristics.
The main discrepancies of our work and Ref. Long:2015cza  consist in the exact
accounting for the universe expansion in the dynamics of axions and
electromagnetic fields. We have also taken into account the nontrivial spectra
of the magnetic helicity and the magnetic energy densities.

On the first glance, there is a similarity of the axion electromagnetism with
the chiral magnetic effect (CME) Boyarsky:2011uy ; Tashiro:2012mf  that was a
motivation for us to apply the Primakoff’s mechanism and study the mutual
influence of an axion zero mode and PMFs in the early Universe.

However, it is not surprising that the CME is significantly more efficient for
the amplification of PMFs than in the case of the axion interaction with
magnetic fields. In the induction Eq. (2), the helicity parameter
$\alpha(t)=(g_{a\gamma}\mathrm{d}\varphi/\mathrm{d}t)/\sigma_{\mathrm{cond}}$,
that enters there the PMF instability term $\alpha(t)(\nabla\times{\bf B})$,
occurs to be much less than the corresponding CME parameter
$\alpha_{\mathrm{CME}}(t)=2\alpha_{\mathrm{em}}\mu_{5}(t)/\pi\sigma_{\mathrm{cond}}$,
where $\mu_{5}=(\mu_{\mathrm{R}}-\mu_{\mathrm{L}})/2$ and $\mu_{\mathrm{R,L}}$
are the chemical potentials for right and left charged leptons. Indeed, the
CME parameter, $\alpha_{\mathrm{CME}}$, is given by a great pseudoscalar
$\mu_{5}/T\sim 10^{-5}$ at the temperature $T_{\mathrm{QCD}}\sim 10^{8}\,{\rm
eV}$, or $\mu_{5}\sim 10^{3}\,{\rm eV}$ (see Fig. 1 in Ref. Boyarsky:2011uy ).
While in the case of axions the pseudoscalar $\alpha$, when accounting for the
initial derivative $\mathrm{d}\varphi/\mathrm{d}t|_{t=t_{0}}\simeq
m_{a}\varphi(t_{0})=m_{a}f_{a}$ [see below Eq. (18)], is estimated as

$\alpha=\left(\frac{\alpha_{\mathrm{em}}}{2\pi}\right)\frac{m_{a}}{\sigma_{\mathrm{cond}}}\ll\alpha_{\mathrm{CME}}=\left(\frac{2\alpha_{\mathrm{em}}}{\pi}\right)\frac{\mu_{5}}{\sigma_{\mathrm{cond}}},$
(34)

where, obviously, the axion mass $m_{a}$ in Eq. (1) is quite small,
$m_{a}\ll\mu_{5}$.

Thus, we conclude that the axion zero mode, $\nabla\varphi=0$, and PMFs evolve
independently each of other. Secondly, for reasonable PMF characteristics, the
axion zero mode that appears in the case when the inflation occurs after the
Peccei-Quinn phase transition, $T_{\mathrm{PQ}}>T_{\mathrm{R}}$,
$T_{\mathrm{PQ}}\simeq 10^{12}\,{\rm GeV}$, can be a candidate to the CDM. If
there is no inflation after the PQ phase transition,
$T_{\mathrm{PQ}}<T_{\mathrm{R}}$, the axion field is spatially varying,
$\nabla\varphi\neq 0$. Axion strings can produce such nonzero momentum modes
of the axion field which could contribute to the CDM Sikivie:2006ni ;
Davidson:2014hfa . The interaction of these modes with PMF is beyond scope of
the present work.

###### Acknowledgements.

We are grateful to V. A. Berezin for remarks concerning description of massive
scalar fields in the general relativity with conformal variables.

##  Appendix A Complete set of evolution equations for axion zero mode and PMF
spectra

In this appendix, we derive the equations for the spectra of the magnetic
energy and helicity, as well as for the axion field.

The action for the axion field $\varphi$ with the mass $m_{a}$, interacting
with the electromagnetic field $A^{\mu}$, in curved spacetime is

$\displaystyle S=$
$\displaystyle\int\mathrm{d}^{4}x\sqrt{-g}\bigg{[}\frac{1}{2}\left(g^{\mu\nu}\partial_{\mu}\varphi\partial_{\nu}\varphi-
m_{a}^{2}\varphi^{2}\right)-\frac{1}{4}F_{\mu\nu}F_{\lambda\rho}g^{\mu\lambda}g^{\nu\rho}$
$\displaystyle-\frac{g_{a\gamma}\varphi}{4}F_{\mu\nu}\tilde{F}{}_{\lambda\rho}g^{\mu\lambda}g^{\nu\rho}-A_{\mu}J_{\nu}g^{\mu\nu}\bigg{]},$
(35)

where $g=\det(g_{\mu\nu})$, $g_{\mu\nu}$ is the metric tensor,
$F_{\mu\nu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}$ is the
electromagnetic stress tensor, $g_{a\gamma}$ is the coupling constant,
$\tilde{F}^{\mu\nu}=\tfrac{1}{2}E^{\mu\nu\alpha\beta}F_{\alpha\beta}$ is the
dual counterpart of $F_{\mu\nu}$,
$E^{\mu\nu\alpha\beta}=\tfrac{1}{\sqrt{-g}}\varepsilon^{\mu\nu\alpha\beta}$ is
the covariant antisymmetric tensor, $\varepsilon^{0123}=+1$, and $J^{\mu}$ is
the external current.

Using Eq. (35), one gets the equations for $\varphi$ and $A^{\mu}$ in the
form,

$\displaystyle\frac{1}{\sqrt{-g}}\partial_{\mu}\left(\sqrt{-g}\partial^{\mu}\varphi\right)+m_{a}^{2}\varphi+\frac{g_{a\gamma}}{4}F_{\mu\nu}\tilde{F}^{\mu\nu}$
$\displaystyle=0,$ (36)
$\displaystyle\frac{1}{\sqrt{-g}}\partial_{\nu}(\sqrt{-g}F^{\mu\nu})+\partial_{\nu}\varphi\tilde{F}^{\mu\nu}+J^{\mu}$
$\displaystyle=0.$ (37)

In Eq. (37), we take into account that $\nabla_{\nu}\tilde{F}^{\mu\nu}=0$,
where $\nabla_{\mu}$ is the covariant derivative. This fact can be proven
basing on the Maxwell equation,
$\nabla_{\mu}F_{\alpha\beta}+\nabla_{\alpha}F_{\beta\mu}+\nabla_{\beta}F_{\mu\alpha}=0$,
and the definition of $\tilde{F}^{\mu\nu}$.

Now we choose the FRW metric, which has the form,

$g_{\mu\nu}=\text{diag}(1,-a^{2},-a^{2},-a^{2}),$ (38)

where $a=a(t)$ is the scale factor. In this case, $\sqrt{-g}=a^{3}$. We also
take the zero mode approximation for the axion field, $\partial_{0}\varphi\neq
0$ and $\partial_{i}\varphi=0$. The four-current $J^{\mu}$ in Eq. (37) can be
represented in the form, $J^{\mu}=(\rho,\mathbf{J}/a)$, where $\rho$ is the
charge density and $\mathbf{J}$ is the electric three-current.

Following Ref. Brandenburg:1996fc , we introduce the conformal variables,

$\mathbf{E}_{c}=a^{2}\mathbf{E},\quad\mathbf{B}_{c}=a^{2}\mathbf{B},\quad\rho_{c}=a^{3}\rho,\quad\mathbf{J}_{c}=a^{3}\mathbf{J}.$
(39)

Then, Eq. (37) takes the form,

$\displaystyle(\nabla\cdot\mathbf{E}_{c})$
$\displaystyle=\rho_{c},\quad\frac{\partial\mathbf{E}_{c}}{\partial\eta}+\mathbf{J}_{c}+g_{a\gamma}\mathbf{B}_{c}\frac{\mathrm{d}\varphi}{\mathrm{d}\eta}=(\nabla\times\mathbf{B}_{c}),$
$\displaystyle(\nabla\cdot\mathbf{B}_{c})$
$\displaystyle=0,\quad\frac{\partial\mathbf{B}_{c}}{\partial\eta}=-(\nabla\times\mathbf{E}_{c}),$
(40)

where $\eta=\smallint\tfrac{\mathrm{d}t}{a}$ is the conformal time. Equation
(40) should be supplied with the relation between $\mathbf{E}_{c}$ and
$\mathbf{J}_{c}$: $\mathbf{J}_{c}=\sigma_{c}\mathbf{E}_{c}$, where
$\sigma_{c}=a\sigma=\text{const}$ and $\sigma$ is the conductivity.

After the volume integration, $V^{-1}\int\mathrm{d}^{3}x(...)$, Eq. (36) for
the axion zero mode can be rewritten in the form,

$\frac{\mathrm{d}^{2}\varphi}{\mathrm{d}\eta^{2}}+\frac{2}{a}\frac{\mathrm{d}a}{\mathrm{d}\eta}\frac{\mathrm{d}\varphi}{\mathrm{d}\eta}+m_{a}^{2}a^{2}\varphi=\frac{g_{a\gamma}}{Va^{2}}\int\mathrm{d}^{3}x(\mathbf{E}_{c}(\eta,{\bf
x})\cdot\mathbf{B}_{c}(\eta,{\bf x})),$ (41)

where we substituted $F_{\mu\nu}\tilde{F}^{\mu\nu}=-4({\bf E}{\bf B})$
multiplying Eq. (36) by $a^{2}$, and assumed that $\varphi$ does not change
under the conformal transformations. The right hand side in Eq. (41) enters
our master Eq. (4) as,

$-\frac{g_{a\gamma}}{2a^{2}}\frac{\mathrm{d}h_{c}(\eta)}{\mathrm{d}\eta}=\frac{g_{a\gamma}}{Va^{2}}\int\mathrm{d}^{3}x(\mathbf{E}_{c}(\eta,{\bf
x})\cdot\mathbf{B}_{c}(\eta,{\bf x})).$ (42)

Such a change follows from the definition of the conformal PMF helicity
density $h_{c}=a^{3}h=V^{-1}\int\mathrm{d}^{3}x{\bf A}_{c}(\eta,{\bf
x})\cdot{\bf B}_{c}(\eta,{\bf x})$ obeying Maxwell Eq. (40). Here
$h_{c}(\eta)=\int\mathrm{d}k_{c}h_{c}(k_{c},\eta)$ is given by a spectrum of
the PMF helicity.

We introduce the conformal spectra of the PMF energy and the PMF helicity
density,

$\displaystyle\rho_{c}(k_{c},\eta)=$
$\displaystyle\frac{k_{c}^{2}}{4\pi^{2}V}\mathbf{B}_{c}(k_{c},\eta)\cdot\mathbf{B}_{c}^{*}(k_{c},\eta),$
$\displaystyle h_{c}(k_{c},\eta)=$
$\displaystyle\frac{k_{c}^{2}}{2\pi^{2}V}\mathbf{A}_{c}(k_{c},\eta)\cdot\mathbf{B}_{c}^{*}(k_{c},\eta),$
(43)

where $k_{c}=ka$ is the conformal momentum, $\mathbf{A}_{c}(k_{c},\eta)$ and
$\mathbf{B}_{c}(k_{c},\eta)$ are the Fourier components of the vector
potential and the magnetic field, and $V$ is the normalization volume.

Basing on Eqs. (40) and (41), as well as using the results of Ref.
Dvornikov:2013bca , we get that,

$\displaystyle\frac{\partial h_{c}(k_{c},\eta)}{\partial\eta}$
$\displaystyle=-\frac{2k_{c}^{2}}{\sigma_{c}}h_{c}(k_{c},\eta)+\frac{4g_{a\gamma}}{\sigma_{c}}\frac{\mathrm{d}\varphi(\eta)}{\mathrm{d}\eta}\rho_{c}(k_{c},\eta),$
$\displaystyle\frac{\partial\rho_{c}(k_{c},\eta)}{\partial\eta}$
$\displaystyle=-\frac{2k_{c}^{2}}{\sigma_{c}}\rho_{c}(k_{c},\eta)+k_{c}^{2}\frac{g_{a\gamma}}{\sigma_{c}}\frac{\mathrm{d}\varphi(\eta)}{\mathrm{d}\eta}h_{c}(k_{c},\eta),$
(44)

and

$\displaystyle\frac{\mathrm{d}^{2}\varphi(\eta)}{\mathrm{d}\eta^{2}}+$
$\displaystyle\frac{2}{a(\eta)}\frac{\mathrm{d}\varphi(\eta)}{\mathrm{d}\eta}\left(\frac{\mathrm{d}a(\eta)}{\mathrm{d}\eta}+\frac{g_{a\gamma}^{2}}{a(\eta)\sigma_{c}}\int_{k_{\mathrm{min}}}^{k_{\mathrm{max}}}\rho_{c}(k_{c},\eta)\mathrm{d}k_{c}\right)$
$\displaystyle+m_{a}^{2}a^{2}(\eta)\varphi(\eta)=\frac{g_{a\gamma}}{a^{2}(\eta)\sigma_{c}}\int_{k_{\mathrm{min}}}^{k_{\mathrm{max}}}k_{c}^{2}h_{c}(k_{c},\eta)\mathrm{d}k_{c},$
(45)

where $k_{\mathrm{min},\mathrm{max}}$ are the ranges of the $k_{c}$ variation.

Now we define the dimensionless variables,

$\kappa=\frac{k_{c}}{k_{\mathrm{max}}},\quad\tau=\frac{2k_{\mathrm{max}}^{2}}{\sigma_{c}}\eta,$
(46)

and the dimensionless functions of these variables

$\displaystyle\mathcal{H}(\kappa,\tau)=$
$\displaystyle\left(\frac{g_{a\gamma}^{2}}{2}\right)h_{c}(k_{c},\eta),$
$\displaystyle\mathcal{R}(\kappa,\tau)=$
$\displaystyle\left(\frac{g_{a\gamma}^{2}}{k_{\mathrm{max}}}\right)\rho_{c}(k_{c},\eta),$
$\displaystyle\Phi(\tau)=$
$\displaystyle\left(\frac{2k_{\mathrm{max}}g_{a\gamma}}{\sigma_{c}}\right)\varphi(\eta).$
(47)

Equations (44) and (45) can be rewritten in the form,

$\displaystyle\frac{\partial\mathcal{H}(\tau,\kappa)}{\partial\tau}=$
$\displaystyle-\kappa^{2}\mathcal{H}(\tau,\kappa)+\frac{\mathrm{d}\Phi(\tau)}{\mathrm{d}\tau}\mathcal{R}(\tau,\kappa),$
$\displaystyle\frac{\partial\mathcal{R}(\tau,\kappa)}{\partial\tau}=$
$\displaystyle-\kappa^{2}\mathcal{R}(\tau,\kappa)+\kappa^{2}\frac{\mathrm{d}\Phi(\tau)}{\mathrm{d}\tau}\mathcal{H}(\tau,\kappa),$
$\displaystyle\frac{\mathrm{d}^{2}\Phi(\tau)}{\mathrm{d}\tau^{2}}+$
$\displaystyle\frac{2}{a(\tau)}\frac{\mathrm{d}\Phi(\tau)}{\mathrm{d}\tau}\left(\frac{\mathrm{d}a(\tau)}{\mathrm{d}\tau}+\frac{1}{2a(\tau)}\int_{\kappa_{m}}^{1}\mathcal{R}(\tau,\kappa)\mathrm{d}\kappa\right)$
$\displaystyle=-\mu^{2}a^{2}(\tau)\Phi(\tau)+\frac{1}{a^{2}(\tau)}\int_{\kappa_{m}}^{1}\kappa^{2}\mathcal{H}(\tau,\kappa)\mathrm{d}\kappa,$
(48)

where $\mu=m_{a}\sigma_{c}/2k_{\mathrm{max}}^{2}$ and
$\kappa_{m}=k_{\mathrm{min}}/k_{\mathrm{max}}$.

## References

  * (1) R. D. Peccei and H. R. Quinn, CP Conservation in the Presence of Instantons, Phys. Rev. Lett. 38, 1440 (1977).
  * (2) S. Weinberg, A New Light Boson?, Phys. Rev. Lett. 40, 223 (1978).
  * (3) F. Wilczek, Problem of Strong $P$ and $T$ Invariance in the Presence of Instantons, Phys. Rev. Lett. 40, 279 (1978).
  * (4) A. Ringwald, Exploring the role of axions and other WISPs in the dark universe, Phys. Dark Univ. 1, 116–135 (2012) [arXiv:1210.5081].
  * (5) P. W. Graham, I. G. Irastorza, S. K. Lamoreaux, A. Lindner, and K. A. van Bibber, Experimental Searches for the Axion and Axion-Like Particles, Annu. Rev. Nucl. Part. Sci. 65, 485–514 (2015) [arXiv:1602.00039].
  * (6) E. Aprile et al. (XENON Collaboration), Excess electronic recoil events in XENON1T, Phys. Rev. D 102, 072004 (2020) [arXiv:2006.09721].
  * (7) D. S. Gorbunov and V. A. Rubakov, Introduction to the theory of the early Universe: Hot Big Bang theory (World Scientific, Singapore, 2011).
  * (8) A. Long and T. Vachaspati, Implications of a Primordial Magnetic Field for Magnetic Monopoles, Axions, and Dirac Neutrinos, Phys. Rev. D 91, 103522 (2015) [arXiv:1504.03319].
  * (9) If the inflation occurs with the reheating temperature, $T_{\mathrm{R}}$, smaller than $T_{\mathrm{PQ}}\sim f_{a}$, $T_{\mathrm{R}}<T_{\mathrm{PQ}}$, the axion field is homogenized over enormous distances. The subsequent evolution of this zero momentum mode is relatively simple in comparison with the multimode case, in which $\nabla\varphi\neq 0$ at $T_{\mathrm{R}}>T_{\mathrm{PQ}}$. In this case, strings, domains as topological defects in the axion field are present Sikivie:2006ni .
  * (10) P. Sikivie, Axion cosmology, Lect. Notes Phys. 741, 19–50 (2008) [astro-ph/0610440].
  * (11) Accounting for the scale factor $a=(1+z)^{-1}$, we vary the undetermined PMF strength $B_{c}=a^{2}B=B_{\mathrm{now}}$, where $B_{\mathrm{now}}(z=0)\leq 10^{-9}\,{\rm G}$ for extragalactic fields Plank . For the initial time at the QCD phase transition the PMF occurs much stronger, $B_{c}(t_{0})=B_{\mathrm{now}}(1+z_{\mathrm{QCD}})^{2}$.
  * (12) For example, substituting the parameter $\kappa_{m}=10^{-10}$ obeying inequalities in Eq. (13) one gets $\tilde{a}=1+3.7\tau$.
  * (13) P. A. R. Ade et al. (Plank Collaboration), Plank 2015 results. XIX. Constraints on primordial magnetic field, Astron. Astrophys. 594, A19 (2016) [arXiv:1502.01594].
  * (14) L. Duffy and K. Bibber, Axions as Dark Matter Particles, New J. Phys. 11, 105008 (2009) [arXiv:0904.3346].
  * (15) The equivalence to physical Fourier variables for the corresponding PMF spectrum $\rho_{\mathrm{B}}(k,t_{0})=C_{1}k^{n_{\mathrm{B}}}$ comes from its normalization on $B^{2}(t_{0})/2=B_{\mathrm{now}}^{2}/2a^{4}(t_{0})$ and the relation $\rho_{c}(\eta=0)=a^{4}(t_{0})\rho_{\mathrm{B}}(t_{0})$ for the conformal PMF energy density $\rho_{c}(\eta=0)$ and the standard one, $\rho_{\mathrm{B}}(t_{0})=\int_{k_{1}}^{k_{2}}\mathrm{d}k\rho_{\mathrm{B}}(k,t_{0})$, where $k$ are physical momenta, $k=k_{c}/a$.
  * (16) A. Brandenburg, K. Enqvist, and P. Olesen, Large-scale magnetic fields from hydromagnetic turbulence in the very early universe, Phys. Rev. D 54, 1291 (1996) [astro-ph/9602031].
  * (17) A. Brandenburg, K. Enqvist, and P. Olesen, The Effect of Silk damping on primordial magnetic fields, Phys. Lett. B 392, 395 (1997) [hep-ph/9608422].
  * (18) A. Neronov and I. Vovk, Evidence for strong extragalactic magnetic fields from Fermi observations of TeV blazars, Science, 328, 73 (2010) [arXiv:1006.3504].
  * (19) W. B. Atwood et al., The large area telescope on the fermi gamma-ray space telescope mission, Astrophys. J. 697, 1071–1102 (2009).
  * (20) K. Jedamzik and L. Pogosian, Reliving the Hubble tension with primordial magnetic fields, Phys. Rev. Lett. 125, 181302 (2020) [arXiv:2004.09487].
  * (21) S. Davidson, Axions: Bose Einstein Condensate or Classical Field?, Astropart. Phys. 65 101 (2015) [arXiv:1405.1139].
  * (22) M. Dvornikov and V. B. Semikoz, Instability of magnetic fields in electroweak plasma driven by neutrino asymmetries, J. Cosmol. Astropart. Phys. 05 (2014) 002 [arXiv:1311.5267].
  * (23) T. Kahniashvili, A. G. Tevzadze, A. Brandenburg, and A. Neronov, Evolution of Primordial Magnetic Fields from Phase Transitions, Phys. Rev. D 87, 083007 (2013) [arXiv:1212.0596].
  * (24) A. Brandenburg, T. Kahniashvili, S. Mandal, A. R. Pol, A. G. Tevladze, and T. Vachaspati, Evolution of hydromagnetic turbulence from the electroweak phase transition, Phys. Rev. D 96, 123528 (2017) [arXiv:1711.03804].
  * (25) A. Boyarsky, J. Fröhlich, and O. Ruchayskiy, Self-consistent evolution of magnetic fields and chiral asymmetry in the early Universe, Phys. Rev. Lett. 108, 031301 (2012) [arXiv:1109.3350].
  * (26) H. Tashiro, T. Vachaspati, and A. Vilenkin, Chiral Effects and Cosmic Magnetic Fields, Phys. Rev. D 86, 105033 (2012) [arXiv:1206.5549].