Interaction of inhomogeneous axions with magnetic fields in the early universe

We study the system of interacting axions and magnetic fields in the early universe after the quantum chromodynamics phase transition, when axions acquire masses. Both axions and magnetic fields are supposed to be spatially inhomogeneous. We derive the equations for the spatial spectra of these fields, which depend on conformal time. In case of the magnetic field, we deal with the spectra of the energy density and the magnetic helicity density. The evolution equations are obtained in the closed form within the mean field approximation. We choose the parameters of the system and the initial condition which correspond to realistic primordial magnetic fields and axions. The system of equations for the spectra is solved numerically. We compare the cases of inhomogeneous and homogeneous axions. The evolution of the magnetic field in these cases is different only within small time intervals. Generally, magnetic fields are driven mainly by the magnetic diffusion. We find that the magnetic field instability takes place for the amplified initial wavefunction of the homogeneous axion. This instability is suppressed if we account for the inhomogeneity of the axion.


Introduction
The existence of axions, which are pseudoscalar particles, was theoretically proposed in Ref. [1] to provide the conservation of the CP symmetry in the quantum chromodynamics (QCD). The conserved CP symmetry is required, in particular, because of the experimental nonobservation of an electric dipole moment of a neutron (see, e.g., Ref. [2]). There are numerous studies of the models for QCD axions which are summarized in Ref. [3].
Currently, there are active searches for axions in laboratory experiments. These particles can be revealed in their interactions with electromagnetic fields, fermion spins, and nuclear electric dipole moments. The difficulty of a direct axion observation is because of the smallness of corresponding coupling constants. Major present experimental studies of axions are reviewed in Ref. [4]. There are recent claims (see, e.g., Ref. [5]) that some experimental data can be interpreted as the axions observation. However, a more careful analysis to confirm such findings is required.
Besides laboratory studies, there are attempts to detect axions from astrophysical sources and cosmological media. One can try to capture axions produced in the Sun using existing

Magnetic field in the presence of inhomogeneous axions
The equations for the electromagnetic field F µν , interacting with a scalar axion field ϕ, in a curved spacetime with the coordinates x µ = (t, x) read whereF µν = 1 2 E µναβ F αβ is the dual counterpart of F µν , E µναβ = 1 √ −g ε µναβ is the covariant antisymmetric tensor, ε 0123 = +1, g = det(g µν ), g µν is the metric tensor, g aγ is the coupling constant, and J µ = (ρ, J/a) is the external current. We consider the Friedmann-Robertson-Walker (FRW) metric g µν = diag(1, −a 2 , −a 2 , −a 2 ),where a(t) is the scale factor. Using the conformal variables [19] E c = a 2 E, B c = a 2 B, ρ c = a 3 ρ, and J c = a 3 J, we rewrite Eq. (2.1) in the form, where the prime means the derivative with respect to the conformal time η defined by dt = adη.
Supposing that J c = σ c E c , where σ c = 10 2 is the conformal conductivity of ultrarelativistic plasma, and neglecting the displacement current E ′ c in Eq. (2.2), we get that the electric field has the form, Here we keep the terms linear in ∇ϕ. Setting ρ c = 0 in Eq. (2.4), we derive the equation for B c , where we use Eqs. (2.3) and (2.6). Equation (2.7) can be simplified if we assume that axions are inhomogeneous but isotropic. It means that, after the spatial averaging, one has that ∇ϕ = 0, whereas (∇ϕ) 2 = 0 and ∆ϕ = 0. Keeping the terms linear in g aγ in Eq. (2.7) and spatially averaging it, we get the modified induction equation, where we omit the averaging brackets for the sake of brevity. In Eq. (2.8), the α-dynamo parameter α(x), which is the coefficient in the magnetic instability term ∝ (∇ × B c ), is different from that in Refs [15,16]. It gets the contribution from the inhomogeneous axion field ∝ ∆ϕ. It is convenient to deal with a Fourier transform, which has the form, is real. Then we define the spectra of the magnetic energy density and the magnetic helicity, Note that the total magnetic helicity, has the standard form.
Taking the Fourier transform of Eq. (2.8) and forming the binary combinations of the fields, one gets the following differential equations: for the spectra ρ p (η) and h p (η) in Eq. (2.9). In Eq. (2.11), α k = ∫ d 3 xe ikx α(x) is the Fourier transform of the α-dynamo parameter in Eq. (2.8). We mention that the equations for the spectra in Eq. (2.9) cannot be written in a closed form, analogous to that, e.g., in Refs. [20,21], for a general coordinate dependent α-dynamo parameter.

Axions in the presence of an electromagnetic field
The evolution of an axion field under the influence of an electromagnetic field obeys the following equation in a curved spacetime: where m is the axion mass. We can rewrite Eq. (3.1) in the form, where we take the FRW metric and use the conformal variables.
In the Fourier representation, Eq. (3.2) takes the form, The integral in the right hand side of Eq. (3.3) differs from the derivative of the magnetic helicity, which is given in Eq. (2.10). Nevertheless, we assume that the electromagnetic field is inhomogeneous but isotropic. Thus all odd correlators, like k i , k i k j k n etc., vanish. In this case, we rewrite Eq. (3.4) as where h p (η) is the spectrum of the magnetic helicity density in Eq. (2.9) and is the scalar function.

Mean field approximation
Equations (2.11) and (3.3) are difficult to analyze in general case. Thus, we apply the mean field approximation. It implies the consideration of nontrivial spectra which are attributed, however, to the zero mode only. First, we assume in Eq.
In this approximation, we set q = 0 in Eq. (2.11) and consider the mean spectra of the magnetic energy and the magnetic helicity, , which were introduced, e.g., in Refs. [20,21]. They are related to the total spectra by the expressions, where ρ p and h p are given in Eq. (2.9). Using Eq. (2.11) in the considered approximation, we get that the mean spectra obey the equations, which formally coincide with those derived in Refs. [20,21]. However, Eq. (4.3) involves a nontrivial spectrum of the axion field in Eq. (4.1). Then, we consider the limit k → 0 in the integrand in Eq. (3.3). In this case, f which is similar to that used in Ref. [16]. However, we keep the term ∝ k 2 ϕ k in the left hand side of Eq. (4.4). Equations (4.4) and (4.3), together with Eq. (4.1), form the closed set of the evolution equations, which define the evolution of the magnetic and axion fields. We shall solve them numerically in Secs. 5 and 6.

Initial condition and the parameters of the system
We consider the evolution of the system after the QCD phase transition which happens at T QCD = 160 MeV. In this situation, the axion mass is taken to be constant, m = 10 −6 eV. We use the parameterization in which the scale factor a = 1/T is dimensional, where T is the primordial plasma temperature. The conformal time is defined up to a constant. We take GeV is the Planck mass, and g * = 17.25 is the number of the relativistic degrees of freedom at T QCD [24]. In this case η(T QCD ) = 0. We assume that the seed spectrum of the magnetic energy is Kolmogorov, ρ 0 (k) ∝ k n with n = −5/3. The seed spectrum of the magnetic helicity is h 0 (k) = 2qρ 0 (k)/k, where 0 ≤ q ≤ 1 is the parameter fixing the initial helicity. Note that k is the dimensionless conformal momentum. It spans the range k min < k < k max , where k min = T QCD /M Pl = 9.1 × 10 −20 is the reciprocal horizon size at T QCD . The maximal conformal momentum k max is a free parameter in our model. We take that k max < r −1 D = 0.1, where r D is the conformal Debye length, to guarantee the plasma electroneutrality.
There is not so much knowledge about the seed axion spectrum. We assume that ϕ There is the approximate relation between m and f a [9]: m = (5.7 ± 0.007) µeV 10 12 GeV fa . Thus, we take f a = 10 12 GeV since we have chosen m = 10 −6 eV. The suggested values of m and f a are in the sensitivity range of the operating and future lumped element detectors [8] (see also Refs. [22,23]). We take that 0 ≤ k < k max in the axion seed spectrum to be able to probe small momenta. Note that ϕ In this case, we can reproduce the zero mode approximation studied in Ref. [16]. We shall consider δ ∼ k max as a rule. Following Ref. [16], we estimate the initial first derivative ∂ η ϕ (0) k basing on the virial theorem. Computing the zero component of the energy-momentum tensor of an axion, we get that We introduce the new variables in Eqs. (4.3) and (4.4), where 0 ≤ τ ≤ τ max = 2k 2 maxMPl /T now , κ m < κ < 1 in the magnetic spectra, whereas 0 < κ < 1 in the axion spectrum, T now = 2.7 K is the current universe temperature, and κ m = k min /k max . Using Eq.
where µ = σcm 2k 2 max T QCD is the effective axion mass, K = σ c /2k max ,ã = 1+λτ , and λ = σcT QCD 2k 2 maxMPl . The initial condition for the new variables in Eq. (5.1) reads where B (0) c = B now /T 2 now = 3.6 × 10 −4 is the seed conformal magnetic field which corresponds to the present day magnetic field B now = 10 −9 G. Such a strength is within the current constrains on the extragalactic magnetic field obtained in Ref. [25] basing on the cosmic rays observations. However, lower bounds on cosmic magnetic fields are considered in Ref. [26]. The initial condition for the dimensionless axion wavefunction is where we take into account that [9] g aγ ≈ αem 2πfa and α em = 7.3 × 10 −3 is the fine structure constant.

Results
In this section, we show the results of the numerical solution of the system in Eq.   The system analogous to that in Eq. (5.2) was studied in Ref. [16] for a spatially homogeneous axion. Here, we are mainly focused on the influence of a nontrivial seed spectrum of an axion on the evolution of the system. That is why we compare the inhomogeneous axion MHD with the homogeneous one. As we mentioned in Sec. 5, the parameter δ ∼ k max . In Fig. 1, we compare the cases of δ = 10 −10 (the inhomogeneous axion) and δ = 0 (the homogeneous axion). The rest of the parameters and the initial condition are the same for both the inhomogeneous and homogeneous cases.
Unfortunately, we can proceed in numerical simulation only for T T QCD for inhomogeneous axions because of the very rapid oscillations of the axion wavefunction (see Fig. 2 below). We checked the dependence of the evolution of the system on the seed magnetic helicity q defined in Sec. 5. It turns out to be negligible. Thus, we set q = 1 which corresponds to the maximally helical fields.
The measurable quantities in the axion MHD are the conformal magnetic field, and the total axion energy in the cooling universe, E a ∝ ∫ d 3 xT 00 , where T 00 is the zeroth component of the axion energy-momentum tensor. We normalize E a on its initial value E (0) a . In case of a homogeneous axion, we deal with the axion energy density ρ a = (∂ τ Φ) 2 +ã 2 µ 2 Φ 2 , which is also normalized on ρ We show B c versus T in the universe cooling down from T QCD in Figs. 1(a) and 1(b) in inhomogeneous and homogeneous cases respectively. The magnetic field decays in both cases, mainly because of the diffusion terms, ∝ −κ 2 R and ∝ −κ 2 H, in Eq. (5.2).
The evolution time in Fig. 1(a) approximately corresponds to that in the inset in Fig. 1(b). We can see the ripple in the inset in Fig. 1(b) (homogeneous axion), whereas the line in Fig. 1(a) (inhomogeneous axion) is smooth. This ripple is not noticeable on great evolution times in Fig. 1(b). Comparing Fig. 1(a) and the inset in Fig. 1(b), we can see that the inhomogeneous case approximately corresponds to the mean value of the homogeneous axion wavefunction.
We depict E a in Fig. 1(c) and ρ a in Fig. 1(d) versus T . By the end of the evolution, both quantities are about 0.1 from their initial values. However, the decay in the inhomogeneous case in Fig. 1(c) happens much faster than in the homogeneous situation in Fig. 1(d). Such a behavior is owing to the additional term ∝ κ ′2 Φ in the integrand in the last line of Eq. (5.2). Therefore the damping of the inhomogeneous axion is greater effectively.
In Fig. 2, we show the evolution of the axion wavefunction integrated over the spectrum, which is again normalized on its initial value at T QCD . One can see that it oscillates vary rapidly. As mentioned above, this fact makes the solution of the system in Eq. (5.2) more complicated than for the homogeneous axion studied in Ref. [16]. Moreover, as one can see in the inset in Fig. 2, there are both rapid oscillations and a sort of beatings in the evolution ofφ. These beatings appear since we integrate different wavenumbers over the spectrum in Eq. (6.3). We take that k max = 10 −10 in our simulations. It means that the present length scale obeys l now > (k max T now ) −1 ≈ 10 −2 R ⊙ , where R ⊙ is the solar radius. Thus the perturbations of the axion field considered in our work, can result in the formation of hypothetical astronomical objects like Bose stars; cf. Ref. [12].

Magnetic field instability
We showed above that there are certain differences in the behavior of the measurable quantities in MHD in the presence of inhomogeneous and homogeneous axions. Nevertheless, these differences are smeared at low temperatures ∼ T now . The axion wavefunction contributes  Fig. 1 except for the initial axion wavefunction which is amplified by the factor Ξ = 5 × 10 7 in both panels now.
to the α-dynamo parameter in Eq. (4.1), which is known to result in the instability of the magnetic field. However, we did not observe this instability in the cases studied above. The behavior of B c was mainly affected by the magnetic diffusion.
It happens since we take that the initial axion wavefunction is ϕ 0 ∼ f a at T QCD . This initial condition is realistic (see, e.g., Ref. [15]). However, the impact of the axion on the magnetic field is small in this case. Now we consider the situation when the initial axion wavefunction is increased. Namely, we take that ϕ (0) k = 8π 3/2 Ξf a exp(−k 2 /δ 2 )/δ 3 , where the factor Ξ = 5 × 10 7 . We also consider the homogeneous axion with ϕ 0 = Ξf a , where Ξ has the same value as in the inhomogeneous case. The initial axion wavefunctions, which are greater than f a , were mentioned in Ref. [15] to potentially influence magnetic fields. We shall see soon that it is the case. The rest of the parameters and the initial condition are the same as in Fig. 1.
In Fig. 3, we show the evolution of the magnetic field defined in Eq. (6.1). We study the inhomogeneous and homogeneous cases in the same time interval. One can see in Fig. 3(b) that the magnetic field becomes unstable. The evolution of the magnetic field in the inhomogeneous case shown in Fig. 3(a) does not reveal the instability.
To illustrate the development of the instability of the magnetic field in the homogeneous case, we plot the spectra of the densities of the magnetic energy and the magnetic helicity in  Figure 4: The spectra of the densities of (a) the magnetic energy; and (b) the magnetic helicity in case of the homogeneous axion. The parameters of the system and the initial condition are the same as in Fig. 3(b). Solid lines are final spectra at the end of the system evolution. Dashed lines are the seed Kolmogorov spectra. Fig. 4. Unlike Eq. (5.1), here we use different new variables, R ′ (κ, τ ) = g 2 aγ T 2 QCD ρ(k, η)/k max and H ′ (κ, τ ) = g 2 aγ T 2 QCD h(k, η)/2. We studied the evolution of magnetic fields with Ξ < 5×10 7 and did not reveal the magnetic field instability, at least in the considered time interval.
In Fig. 4, by dashed lines, we show the seed Kolmogorov spectra which correspond to T QCD . We also depict the final spectra which are at the end of the evolution of the system. One can see that parts of the spectra at great momenta, i.e. at small scales, increase. One can say that it is equivalent to the direct energy cascade triggered by the homogeneous axion.
The α-dynamo parameter in Eq. (4.1) has two axion contributions in the integrand: (i) ∝ ϕ ′ k ; and (ii) ∝ k 2 ϕ k . Despite we take a greater axion wavefunction, the inhomogeneous axion oscillates more rapidly; cf. Fig. 2. Thus, the term (i) prevails over the term (ii). Since the term (i) is sign changing, the consideration of the inhomogeneous axion prevents the development of the magnetic field instability.

Conclusion
In the present work, we have studied the evolution of the system which consists of axions and magnetic fields. Both axions and magnetic fields are spatially inhomogeneous. The evolution of such a system was considered in the early universe cooling down from the QCD phase transition, with the universe expansion rate being accounted for exactly.
First, in Sec. 2, we have started with the general equations for an electromagnetic field, which interacts with axions, in an arbitrary curved spacetime. Then, taking the FRW metric, we have derived the equations for the conformal electric and magnetic fields. Basing on this result, we have obtained the general induction Eq. (2.7) for the magnetic field, which accounts for the spatial distribution of axions. This equation has been further modified to Eq. (2.8) in the isotropic medium. Then, we have obtained the general Eq. (2.11) for the spectra of the magnetic energy and the magnetic helicity.
To describe the evolution of axions in the presence of an electromagnetic field, in Sec. 3, we have considered the general Eq. (3.1) for ϕ in a curved spacetime. This equation has been rewritten in the Fourier space in Eq. (3.3). An arbitrary magnetic field contributes to the evolution of axions via the conformal time derivative of the magnetic helicity spectrum in Eq. (3.5).
In order to get the closed system of the evolution equations for a magnetic field and an axion, in Sec. 4, we have applied the mean field approximation which consists in the consideration of zero mode perturbations. In this situation, we could define the mean spectra, ρ(k, η) and h(k, η), in Eq. (4.2), as well as derive Eqs. (4.3) and (4.4). Formally, Eqs. (4.3) and (4.4) coincide with those obtained in Ref. [16]. However, they account for the nontrivial spectrum of the α-dynamo parameter.
In Sec. 5, we have chosen the initial conditions and have defined all the parameters in Eqs. (4.3) and (4.4). We have considered the evolution of the system below the QCD phase transition when the axion mass takes its vacuum value. It should be noted that there is not so much information about the seed spectrum of axions. We have adopted the initial Gaussian distribution for ϕ k which is characterized by the phenomenological parameter δ. One can reproduce, e.g., the zero mode case, discussed in Refs. [15,16], by setting δ → 0. For the magnetic field, we have chosen the seed Kolmogorov spectrum with realistic parameters. Finally, we have rewritten the evolution equations in the form convenient for numerical simulations; cf. Eq. (5.2).
The numerical solution of Eq. (5.2) has been present in Sec. 6. We have been mainly interested in the dependence of the results on the seed axion spectrum. For this purpose we have compared the cases of the inhomogeneous, with δ = 10 −10 , and the homogeneous, with δ = 0, axions. In both cases, the magnetic field evolution is different only in small time intervals. It is smooth for inhomogeneous axions and reveal a ripple for homogeneous axions; cf. Figs. 1(a) and 1(b). Globally, the evolution of magnetic fields is driven mainly by the diffusion terms in Eq. (5.2). The energy of the inhomogeneous axion decays faster than for the homogeneous one; cf. Figs. 1(c) and 1(d).
The possibility to get the magnetic field instability has been considered in Sec. 6.1. For this purpose we have taken the initial axion wavefunction amplified by the factor Ξ = 5 × 10 7 . The impact of axions on the magnetic field was mentioned in Ref. [15] to be possible if the initial axion wavefunction is greater than f a . We have obtained that the magnetic field can be unstable in the case of a homogeneous axion; cf. Fig. 3(b). This instability is washed out for the inhomogeneous axion (see Fig. 3(a)) owing to the very rapid sign alteration of the α-dynamo parameter.