Evolution of axions in the presence of primordial magnetic fields

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 [1][2][3]. In addition, the axion also provides a possible candidate for the cold dark matter (CDM) of the Universe [4]. There are active searches for axions [5]. A progress in the experimental studies of these particles is reported in Ref. [6] 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 PQ → T QCD , where T PQ ≃ f a = (10 10 ÷ 10 12 ) GeV is the temperature at the Peccei-Quinn phase transition, T QCD ≃ (100 ÷ 200) 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 π 0 mesons its "cold" (fixed) value m a (0) = m a just near the initial time in our problem t 0 = t QCD , m a = 6 × 10 −6 eV 10 12 GeV f a . (1.1) In its turn, PMFs originate from the hypermagnetic fields B Y which exist before the electroweak phase transition (EWPT), T > T EW ≃ (100÷160) GeV. Such PMF, being a seed for observable galactic magnetic fields B ∼ 10 −6 G, can be strong enough at T QCD . Returning the present bound on the extragalactic PMF, B ≤ B now = 10 −9 G, back to early universe times, one gets the PMF strength B QCD = B now (1 + z QCD ) 2 ≃ 1.6 × 10 14 G for the red shift z QCD = 4 × 10 11 at T QCD ∼ 100 MeV, or B EW = 1.6 × 10 20 G at T EW ∼ 100 GeV. * maxdvo@izmiran.ru † semikoz@yandex.ru On one hand, such a strong field can not influence the Universe expansion since the PMF energy density ρ B = B 2 /2 is much less than the ultrarelativistic matter energy density ρ = g * π 2 T 4 /30, ρ B ≪ ρ. Here g * is the effective number of degrees of freedom [7, p. 95]. On the other hand, it is intriguing to check whether such a field can be enlarged by the axion instability term ∼φ∇×B entering the induction equation [8], 2) or, vice versa, the PMF can influence evolution of that uniform axion field, ∇ϕ = 0, Here g aγ ≈ α em /2πf a is the coupling constant for axion interaction with electromagnetic fields, α em ≈ (137) −1 is the fine structure constant, σ cond ≃ 100T is the electric conductivity in the hot universe plasma. The problem of the mutual influence of axions and PMFs was discussed in Ref. [8]. However, in Ref. [8], 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 ϕ(η) and the spectra of the conformal PMF energy density ρ c (k c , η) and the conformal PMF helicity density h c (k c , η). Here η 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 B, 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 now = 10 −9 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 ∇ϕ = 0, or for the zero axion mode [9], the axion evolution equation in the presence of PMF and in the Friedmann-Robertson-Walker (FRW) metric, ds 2 = a 2 (t)(dη 2 − dx 2 ), where a is the scale factor and dη = dt/a(t), takes the following form in the variables x µ = (η, x) (see Eqs. (A7) and (A8) in Appendix A): is the comoving magnetic field in conformal variables [11]. We can calculate the derivative dh c (η)/dη for the magnetic helicity density h c (η) = dk c h(η, k c ) given by its isotropic Fourier spec- where k c = ak = 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), where the magnetic energy density ρ c (η) = dk c ρ c (η, k c ) = B 2 c (η)/2 is given by the isotropic Fourier spectrum ρ c (k c , η) = k 2 c B c (k c , η)B * c (k c , η)/(4π 2 V ), k c = |k c | is the absolute value of the comoving Fourier momentum, and σ c = aσ cond is given by the electric conductivity σ cond ≃ 100T in the hot universe plasma.
Substituting the derivative dh c (η)/dη = kmax kmin dk c ∂ η h c (η, k c ), with the integrand given by Eq. (2.2), we can rewrite the axion evolution Eq. (2.1) as that is the generalization of the axion wave Eq. (40a) in Ref. [8] 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. (2.4), is completed by the differential Eqs. (2.2) and (2.3) for the spectra h c (k c , η) and ρ c (k c , η). Here, in the causal scenario, the momentum limits k min = k 1 a(t) and k max = k 2 a(t) obey the inequality k 2 > k 1 = l −1 H (T 0 ) = const where the horizon size l H (T 0 ) is chosen at the QCD phase transition time t 0 = (T 0 /MeV) −2 s when the axion mass m a (T ) reaches its maximum ("cold") value, given in Eq. (1.1), during cooling of the Universe, T → T 0 . We consider such a time t 0 ≃ 10 −4 s as the initial moment in our simulations. Correspondingly, for the conformal time during radiation epoch, η(t) = t t0 dt/a(t) + η(t 0 ), for which a(t) = a(t 0 ) t/t 0 is the scale factor in FRW metric, we put the initial value η(t 0 ) = 0 resulting in where the initial scale factor a(t 0 ) = T now /T 0 = (1 + z QCD ) −1 is given by by the red shift at the QCD phase transition, z QCD ≃ 4 × 10 11 , since T now = 2.725 K and T 0 = T QCD = 10 12 K ≃ 100 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 −1 H = (2/3) × 10 −20 GeV given by the horizon size l H , or maximum PMF scale Λ where the conductivity σ cond = σ c /a(t 0 ) was substituted.
The system of the evolution Eqs. (2.3)-(2.4), rewritten for the three dimensionless functions H(κ, τ ) = Here, in Eq. (2.9), µ = m a σ c /2k 2 max is the dimensionless parameter given by the axion mass in Eq. (1.1). The lower limit κ m = k min /k max = k 1 /k 2 for the dimensionless Fourier momentum κ = k c /k max = k/k 2 obeys the constraints where we account for Eq. (2.6). 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 κ m in Eq. (2.10), or by varying k max = k min /κ m for the fixed k min = a(t 0 )k 1 ; • the helicity parameter q, 0 ≤ |q| ≤ 1, entering the initial magnetic helicity spectrum (see Eq. (2.22) 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 now at the present red shift z = 0 originated in our scenario from the initial PMF B 0 = B now (1 + z QCD ) 2 at z QCD = 4 × 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. (2.5), the scale in the radiation era, where τ ≥ 0 [12].

A. Axion energy density
The important characteristic of axions is the axion energy density ρ a (t) = T 00 , where T 00 is the component of the energy-momentum tensor for a scalar field, In case of a spatially homogeneous field, basing on Eq. (2.12), one gets that where k max = k 2 a = k 1 a/κ m is given by the parameter κ m in Eq. (2.10) and the minimal momentum k 1 in our causal scenario at t 0 = t QCD , GeV. This results in the dimensional factor ahead the brackets in Eq. (2.13), that depends on the PMF spectrum width κ m .

B. Initial conditions
The initial axion zero mode amplitude ϕ(t 0 ) ≃ f a , θ(t 0 ) ≃ 1, see, e.g., Ref. [8], corresponds to where κ m is the free parameter constrained in Eq. (2.10). The initial derivative dΦ/dτ | τ =0 , necessary for solving Eq. (2.9), is derived using the virial theorem θ2 = m 2 a θ 2 resulting from the axion energy density in Eq. (2.13) around its potential minimum [14], We put roughly (dϕ/dt) t=t0 = m a ϕ(t 0 ) relying on that virial theorem. Then, one gets where we substituted the "cold" axion mass m a in Eq. (1.1) at the moment of the QCD phase transition, Note that the mass term, (µa 2 )Φ, in Eq. (2.9) is rather great. We shall see in Sec. III that this term is dominant.
The initial PMF spectra R(κ, τ = 0) and H(κ, τ = 0), h c (k c , η = 0), (2.20) are given by the PMF spectra in the known form [15], The constant C for the spectrum in Eq. (2.21) is given by normalization on the initial magnetic energy density (2.24) We consider below the Kolmogorov's spectrum with n B = −5/3.

III. RESULTS
In this section, we present the numerical solution of Eqs. (2.7)-(2.9) with the initial conditions fixed in Sec. II B.
Notice that, at the initial time t 0 = t QCD chosen in our problem for the causal scenario, the maximum spatial PMF scale Λ (max) B = l H (t 0 ) ≃ 3 × 10 6 cm occurs too small after its following growth till present time. It grows upto the scale Λ B (z = 0) ≃ 1 pc only since the horizon expands during the cooling much faster, l H ∼ T −2 , than the correlation length Λ B ∼ 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 Λ B as a measure of the coherence length of the magnetic field, can increase in the coherence by the five orders of magnitude [16]. This fact allows us to get closer to the present scale Λ B ∼ 1 Mpc for PMF bounded in Ref. [13]. 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. [17], where nonlinear effects were shown to prevent most likely this problem from happening.
In Fig. 1(a), the normalized axion field Φ(τ )/Φ(0) = ϕ(t)/ϕ(t 0 ) =θ(t), starting fromθ(t 0 ) = ϕ(t 0 )/f a = 1, grows forty times just at the beginning (see inset in Fig.  1(a)), and then reduces soon to the same valueθ(t) ∼ 1 acceptable for a future axion contribution to CDM. This result does not depend on the presence of PMF, see in Fig.1(b) where we put g aγ = 0 and q = 0 excluding interaction of axions with the PMF. In both cases, this evolution ofθ(t)-amplitudes happens during lepton era in the hot universe plasma.

(3.3)
For κ m = 10 −10 the huge dimensional factor C m = (2/9) × 10 10 GeV 4 in Eq. (2.13) is compensated in the product ρ a (t 0 ) = C mρ (τ = 0) by the small valueρ(τ = 0) = 10 −14 calculated numerically from the system of evolution Eqs. (2.7)-(2.9), or ρ a (t 0 ) = C mρ (τ = 0) = (2/9) × 10 −4 GeV 4 is the initial axion energy density. On the other hand, the ultra-relativistic matter energy density ρ = π 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 QCD = 0.1 GeV at the same moment t 0 , ρ = 1.18 × 10 −4 GeV 4 . Thus, we get at the start a remarkable ratio of energy densities for the axion as a candidate to CDM, ρ a (t 0 )/ρ = 0.19. Notice that this ratio does not depend on a width of the PMF Fourier spectrum κ m in Eq. (2.10) since that parameter shrinks in the initial axion energy density ρ a (t 0 Then the axion energy density ρ a ∼ρ(τ ) reduces as seen in Fig. 2 as well as the matter density ρ decreases due to the universe cooling, ρ ∼ T 4 , therefore supporting approximately the same ratio ρ a (t)/ρ(t) ≃ 0.2.
We are not guaranteed to observe in future the maximum PMF B now = 10 −9 G while some predictions with a lower bound for B 1 Mpc (z = 0) > 10 −16 G exist [18]. They arise, in particular, from γ-ray observation in the Fermi experiment [19]. Recently, in Ref. [20], it was suggested to look for PMF at a level of the comoving strength B now = Ba 2 = const ≃ 10 −11 G. In the last photon scattering (for which CMB instruments are sensitive), such a field corresponds to the PMF strength B(z r ) = B now (1 + z r ) 2 ≃ 10 −5 G, whereas, in our problem, to B 0 = B now (1 + z QCD ) 2 ≃ 1.6 × 10 12 G.

IV. DISCUSSION
In the present work, we have studied the evolution of a zero axion mode, ∇ϕ = 0, in a hot plasma of the early Universe in the presence of PMF that do not exceed B now = 10 −9 G at the present time. We conclude that such a mode can contribute to the axion CDM since its amplitude tends toθ(t) = ϕ(t)/f a → 1 evolving during radiation era, see Fig. 1. Simultaneously the corresponding axion energy density ρ a (t) reduces in such a way that its part in the total matter energy density ρ tends to the CDM ratio ρ a (t)/ρ ≃ 0.2, see comments on that after Eq. Notice that our PMF amplitudes 1 κm dκR(κ, τ ) ∼ Ba 2 = const correspond to the comoving strength for a given B now = Ba 2 = const and κ m . Therefore, the ratio does not depend on the conformal time η ∼ τ , being slightly violated by PMF interaction with axions. The similar conservation concerns the magnetic helicity density, 1 κmH (κ, τ )dκ ∼ ha 3 = const, calculated from the same self-consistent system of Eqs. (2.7)-(2.9) for each given B now and κ m . Hence, the magnetic helicity density for the case q = 0, where the constant depends on the initial conditions. The conservation of the magnetic helicity density in Eq. (4.3) is not surprising in the absence of the axion interaction with PMF's through Primakoff's effect. Indeed, the comoving mean energy density ε M = B 2 /2 = ∫ dk c ρ c (k c , η) decays with time such that ε M ∝ η −2/3 while the correlation length grows like ξ M ∝ η 2/3 as it follows from the dimensionless analysis based on the Kolmogorov-type approach for the magnetic turbulence, see, e.g., Refs. [23,24]. This results in conservation of the magnetic helicity density, h c (η) = ∫ dk c h c (k c , η) = B 2 ξ M = 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γ = α em /2πf a ; cf. Eqs (4.2) and (4.3).
Equations (4.1)-(4.3) were directly checked by our numerical simulations. Thus, we confirm the conclusion of Ref. [8] that the axion zero mode (∇ϕ = 0) practically does not influence any PMF characteristics. The main discrepancies of our work and Ref. [8] 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) [25,26] 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. (1.2), the helicity parameter α(t) = (g aγ dϕ/dt)/σ cond , that enters there the PMF instability term α(t)(∇ × B), occurs to be much less than the corresponding CME parameter α CME (t) = 2α em µ 5 (t)/πσ cond , where µ 5 = (µ R − µ L )/2 and µ R,L are the chemical potentials for right and left charged leptons. Indeed, the CME parameter, α CME , is given by a great pseudoscalar µ 5 /T ∼ 10 −5 at the temperature T QCD ∼ 10 8 eV, or µ 5 ∼ 10 3 eV (see Fig. 1 in Ref. [25]). While in the case of axions the pseudoscalar α, when accounting for the initial derivative dϕ/dt| t=t0 ≃ m a ϕ(t 0 ) = m a f a [see below Eq. (2.15)], is estimated as where, obviously, the axion mass m a in Eq. (1.1) is quite small, m a ≪ µ 5 . Thus, we conclude that the axion zero mode, ∇ϕ = 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 PQ > T R , T PQ ≃ 10 12 GeV, can be a candidate to the CDM. If there is no inflation after the PQ phase transition, T PQ < T R , the axion field is spatially varying, ∇ϕ = 0. Axion strings can produce such nonzero momentum modes of the axion field which could contribute to the CDM [10,21]. The interaction of these modes with PMF is beyond scope of the present work.

ACKNOWLEDGMENTS
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 ϕ with the mass m a , interacting with the electromagnetic field A µ , in curved spacetime is where g = det(g µν ), g µν is the metric tensor, F µν = ∂ µ A ν − ∂ ν A µ is the electromagnetic stress tensor, g aγ is the coupling constant,F µν = 1 2 E µναβ F αβ is the dual counterpart of F µν , E µναβ = 1 √ −g ε µναβ is the covariant antisymmetric tensor, ε 0123 = +1, and J µ is the external current.
Following Ref. [16], we introduce the conformal variables, E c = a 2 E, B c = a 2 B, ρ c = a 3 ρ, J c = a 3 J. (A5) Then, Eq. (A3) takes the form, where η = ∫ dt a is the conformal time. Equation (A6) should be supplied with the relation between E c and J c : J c = σ c E c , where σ c = aσ = const and σ is the conductivity.
After the volume integration, V −1 d 3 x(...), Eq. (A2) for the axion zero mode can be rewritten in the form, Here h c (η) = dk c h c (k c , η) is given by a spectrum of the PMF helicity. We introduce the conformal spectra of the PMF energy and the PMF helicity density,