Stability of axion dark matter-photon conversion

It is known that a coherently oscillating axion field is a candidate of the dark matter. In the presence of the oscillating axion, the photon can be resonantly produced through the parametric amplification. In the universe, there also exist cosmological magnetic fields which are coherent electromagnetic fields. In the presence of magnetic fields, an axion can be converted into a photon, and vice versa . Thus, it is interesting to investigate what happens for the axion-photon system in the presence of both the axion dark matter and the magnetic fields. This system can be regarded as a coupled system of the axion and the photon whose equations contain the Mathieu type terms. We find that the instability condition is changed in the presence of magnetic fields in contrast to the conventional Mathieu equation. The positions of bifurcation points between stable and unstable are shifted and new instability bands appear. This is because the resonantly amplified axion can be converted to photon, and vice versa .


I. INTRODUCTION
The cosmological dark matter problem has been studied in the context of beyond the standard model of particle physics. One of such dark matter candidates is an axion which has been originally proposed as a solution for the strong CP problem [1][2][3][4]. This original axion is called a QCD axion. String theory predicts axionlike particles (ALPs) with a broad mass range [5,6]. Throughout this paper, we will simply use a word "axion," for both QCD axion and ALPs. The axion has feeble interaction with standard model particles and could be produced in the early universe by nonthermal mechanism. This is the reason why the axion can be the dark matter [7][8][9]. In particular, it is known that ultralight axions called fuzzy dark matter [10] can resolve the issues in ΛCDM, e.g., the core-cusp problem and the missing satellite problem. We can treat axion dark matter as a classical field. Then, the axion is coherently oscillating with a frequency determined by the mass. There are various experiments to search for the axion dark matter [11,12].
It is known that, in the presence of the axion dark matter, the propagation of photons is governed by the Mathieu equation [13]. The properties of the Mathieu equation are well studied in mathematics [14][15][16]. It is known that the system becomes unstable for specific parameter regions.
In the universe, on top of the axion dark matter which is a coherent axion field, there exist cosmological magnetic fields which is a coherent electromagnetic field. Remarkably, in the presence of magnetic field, there occurs the axion-photon conversion [17,18]. The axion-photon conversion has been investigated in the context of astrophysics. Indeed, the axion-photon conversion can explain the fact that high energy photons can reach the Earth through intergalactic magnetic fields without disappearing. On the other hand, in the CAST experiment [19], strong magnetic field is applied to the detector in order to detect axions produced in the sun by converting axions into photons. The fact that no signal of axions has been detected until now has given constraints on the mass of the axion and the coupling constant between an axion and two photons.
As we explained in the above, it is natural to consider the axion dark matter and the magnetic fields at the same time. Hence, in this paper, we investigate what happens for the axion-photon system in the presence of both the axion dark matter and the magnetic fields. More precisely, we study the stability of such system in terms of both numerical and analytical methods. Although there are related papers which discuss behavior of axion dark matter and photon with and without magnetic field [20][21][22][23][24][25][26][27][28][29][30][31], to the best of our knowledge, no one did the stability analysis focusing on the axion-photon conversion.
This paper is organized as follows. In Sec. II, we introduce basic equations of axion electrodynamics. We derive basic equations by separating a background and perturbed quantities. Then, we show numerical results in Sec. III. They show stability of the solutions for the basic equations. In Sec. IV, we give an analytical derivation of the numerical results. In particular, we will show you how to determine the boundaries between stable and unstable region in the parameter space. In Sec. V, we will discuss an interpretation of our numerical results and a possible application. The final Sec. VI is devoted to the conclusion.

II. AXION ELECTRODYNAMICS
In this section, we introduce basic equations of axion electrodynamics. Then, we consider an oscillating axion field and a static uniform magnetic field as a background. Given the background, we derive perturbative equations for describing propagation of axions and photons. With these equations, we can study mixing between axions and the photons and the stability of the system.

A. Basic equations of axion electrodynamics
We consider the following system: where a is an axion field with mass m a , and g aγγ is a coupling constant. The field strength F μν of the electromagnetic field A μ ð⃗ x; tÞ and its dualF μν are given by Using the potential A μ ð⃗ x; tÞ ¼ ½−ϕð⃗ x; tÞ; ⃗ Að⃗ x; tÞ, the electric and magnetic fields are defined by Here, we have chosen the Lorenz gauge: ∇ · ⃗ Að⃗ x; tÞ þ ∂ t ϕð⃗ x; tÞ ¼ 0: Equations (5) and (6) are basic equations of axion electrodynamics [32]. For the full analysis, we need to resort to lattice calculations. Here, we use the perturbative analysis.

B. Background equations
Now, we assume both the axion dark matter and the magnetic fields as a background. The background magnetic field is static and uniform, We introduce here coordinate basis so that the propagation is in the direction ⃗ e z ¼ ½0; 0; 1, one of the rests is parallel to the magnetic field: ⃗ e k ¼ ½1; 0; 0, and the other is ⃗ e ⊥ ¼ ½0; 1; 0. The background equation for axion is and for photon Note that we have chosen the radiation gauge This is because the source term of the scalar potential ϕðtÞ equation vanish Solving the Eq. (10), we see that the electric field is induced by the axion oscillation: Substituting (13) into (9), we can geẗ where we replace time variable t with τ ≡ m a t, and express a derivative with respect to τ by dot. Here we also introduced new dimensionless parameters β and Ω β as ð15Þ Note that conservatively we have a constraint g aγγ ≤ 10 −11 GeV −1 . Recall the relation 1G ¼ 1.95 × 10 −2 eV 2 , for the cosmological magnetic fields ∼nG and the axion mass m a > 10 −22 eV, we can neglect the effect of magnetic fields β ≪ 1. However, for more strong magnetic fields, we need to consider the effect of β.
In the end, there are uniform static magnetic field, oscillating axion field and oscillating electric field in the background: a 0 ðτÞ ¼ā cos ðΩ β τÞ; ð17Þ E 0k ðτÞ ¼ m a βā cos ðΩ β τÞ: We introduce energy density ρ as follows: then background energy density ρ BG is given by We determine axion amplitudeā by the energy density ρ, Thus, we found the following expressions The first order equations of (5) and (6) are given by Although the time translational symmetry is broken by the time dependent coherent oscillation of the axion field, the system has the spatial translation invariance. Hence, it is useful to use Fourier transformation where α denotes k or ⊥. Using this transformation, we can write the equations as follows: Now, we need to substitute the background solutions (22) and (23) into Eqs. (36)-(38). Then, we get following equations: where we introduced dimensionless parameters κ and ϵ as follows: From now on, for simplicity, we use an approximation neglecting higher order terms in β and ϵ, namely, we take into account up to the first order in β and ϵ. Thus, we obtain δÄ k ðk; τÞ þ κ 2 δA k ðk; τÞ ¼ −βδ _ aðk; τÞ − iϵ sinðτÞδA ⊥ ðk; τÞ; ð44Þ δÄ ⊥ ðk; τÞ þ κ 2 δA ⊥ ðk; τÞ ¼ iϵ sinðτÞδA k ðk; τÞ: We can see that when ϵ ¼ 0, Eqs. (43)-(45) describe the axion-photon conversion [17,18]. For β ¼ 0, they describe photon propagation in the presence of only axion dark matter [13]. Taking the circular polarization basis we see original Eqs. (39)-(41) are rewritten as follows: δÄ L=R ðk; τÞ þ ½κ 2 AE ϵ sinðτÞδA L=R ðk; τÞ Under the approximation we are considering, Eqs. (43)-(45) are rewritten as follows: Here, we should mention the previous work [27]. They investigated the similar system, but they neglected the parametric resonance. In this paper, we consider Mathieu type terms and focus on the resonance instability.

III. STABILITY ANALYSIS-INCE-STRUTT CHART
In this section, we numerically investigate the behavior of solutions for the basic equations Eqs. (43)-(45). First, we give a short review of the Mathieu equation for comparison. Next, we show numerical results for axion dark matter-photon conversion which have both similarities and differences with the Mathieu equation's. In the next Sec. IV, we will provide analytical derivation of the numerical results.

A. Without background magnetic field
A photon propagating in the axion dark matter obeys following equations [13].
This can be obtained by putting β ¼ 0 in Eq. (50). The equation (51) represents harmonic oscillator whose frequency also oscillates, and this type of equation is called the Mathieu equation [14][15][16]. The solutions can be stable or unstable, depending on dimensionless parameters, κ and ϵ. The Floquet theorem [33] divide the (κ − ϵ) plane into two regions ( Fig. 1), stable and unstable, and this chart is called Ince-Strutt chart [34,35]. Please refer the reader to [16] for the Floquet theorem and the Ince-Strutt chart.
The bifurcation points on the κ axis appear at and the boundaries between the stable and unstable region are called transition curves. On the transition curves, the equation (51) has periodic solutions. Here, we introduce dimension less parameter χ: whereκ is given by Eq. (52). For nonzero ϵ ≠ 0, a wave number κ deviates fromκ in order for the solution of Eq. (51) to still have a period T ¼ 2π=κ, and the deviation is represented by χ.
For example, on the transition curves originated at κ ¼ 1=2, there is a periodic solution with T ¼ 4π. For small jϵj, the transition curves are approximately given by In the case ofκ ¼ 1, the transition curves are given by On these curves, (51) has a periodic solution with T ¼ 2π.
In the following two sections, we will investigate how these results are changed when the background magnetic field is taken into account by solving Eqs. (43)-(45) numerically. At the same time, we show some transition curves, and its analytical derivation is given in Sec. IV.

B. Shift of bifurcation points
The Fig. 2 shows that bifurcation points of transition curves appear again aroundκ ¼ n=2ðn more precise, bifurcation points are shifted even on the κ axis (ϵ ¼ 0) due to the background magnetic field.
As can be seen from Fig. 2(a), the starting point of transition curvesκ ∼ 1=2 is shifted by magnetic field as In the case ofκ ∼ 1 [ Fig. 2(b)], the unstable region splits into two regions.
The first two curves are exactly the same as the conventional one (55). On the other hand, the other two curves are shifted by magnetic field β. The region which intervene between χ 1;− and χ 1;þ represents the instability of parallel photon component which does interact with axion through magnetic field. We verified that transition curves (56)-(57) and (58)-(61) are still valid for full equations (39)-(41). However, in the case of (39)-(41), a new bifurcation point appears between κ ¼ 1 and κ ¼ 1 þ β 2 =2 due to higher order contributions.

C. New bifurcation points
It seems that the axion dark matter-photon conversion has other bifurcation points. From our numerical calculations, we empirically found the condition for the bifurcation points ffiffiffiffiffiffiffiffiffiffiffiffi ffi 1 þκ 2 Solving this with respect toκ, we get the following relation: In the case of n ¼ 2, we depicted the unstable region in Fig. 3. We will see that transition curves can be derived in an analytical way in the next section as follows: We checked that the Eq. (64) is still valid for full equations (39)-(41). However, in the case of (39)-(41), the width of unstable band is more broader due to more higher order contributions. Moreover, a new bifurcation point appears on the left side of κ ¼ 3=4 − 5β 2 =32. In this paper, we shall restrict ourselves to the instability of the leading order equations (43)-(45).

IV. ANALYTIC EXPRESSIONS OF TRANSITION CURVES
As we have seen in the previous section, there are differences between the conventional Mathieu equation and axion dark matter-photon conversion. However, the results are obtained numerically. In this section, we would like to give an analytical support to our findings. We show how the boundaries between stable and unstable regions are determined by treating parameters χ, β, ϵ as small quantities.

5: ð68Þ
Let us put an ansatz ⃗ xðτÞ ¼ e λτ ½The superpositions of many photon's overtones:; and substitute it into the Eq. (67). In order for (67) to have a nontrivial solution, i.e., ⃗ x ≠ 0, the determinant of the coefficient matrix obtained in this way must vanish. Here, we introduce growth rate λ ∈ C under the condition, jχj ∼ jλj ≪ 1. It is the real part ℜ½λ that determines the stability of the solutions to the Eq. (67), and the imaginary part ℑ½λ detune the frequency of the solutions. The criteria of the stable and unstable is given as follows: In the case of ℜ½λ > 0, the growth rate after one period T is given by roughly e Tℜ½λ , ⃗ xðτ þ TÞ ≃ e Tℜ½λ ⃗ xðτÞ: From the explicit calculations, it turns out that the determinant of the coefficient matrix depends only on q ≡ λ 2 . Therefore, the criteria (69) can be replaced by following ones: all other condition∶ unstable: ð71Þ Note that small parameters χ, β, ϵ, λ may have different relative magnitude relationship depending on bifurcation points. Before moving on to the concrete analysis, we define the 3 × 3 matrices which compose coefficient matrix: ð72Þ where I 3 denotes the identity matrix, and n is non-negative integer, n ¼ 0; 1; 2; …. The matrices K, B, E have been already defined in (68).
det½R 1=2 ðq; χ; β; ϵÞ ¼ 0: From numerical results, we see the hierarchy of the order λ ∼ χ ∼ ϵ ∼ β 2 . The leading order contribution to the determinant is given by det½R 1=2 ðq; χ; β; ϵÞ leading The next leading order is not relevant here, however, you will soon see that it should be taken into account when you consider the new bifurcation points. In the present case, we have the following next order contribution: det½R 1=2 ðq; χ; β; ϵÞ next leading Evaluating the determinant at the leading order det½R 1=2 ðq; χ; β; ϵÞ leading ¼ 0; we get four λ: Therefore, the range where the criterion for stability (69) is broken is as follows: and the transition curves are derived from the condition ℜ½λ ¼ 0, Here, we would like to comment on the growth rate (80). If you consider the situation where there is only coherent oscillating axion in the background like [13,30,31], then the growth rate does not depend on axion mass m a . In fact, if you choose β ¼ 0 in (80), then you can confirm this fact. However, in our case β ≠ 0, note that growth rate m a ℜ½λ become to depend on axion mass due to the presence of background magnetic field β.
det½R 3=4 ðq; χ; β; ϵÞ ¼ 0: From numerical results, we see the hierarchy of the order λ ∼ χ ∼ ϵ 2 ∼ β 2 . Evaluating the determinant at the leading order det½R 3=4 ðq; χ; β; ϵÞ leading ¼ 0; we get a cubic equation for q ≡ λ 2 , We obtained solutions as follows: Equation (96) has three negative real solution q 1 , q 2 , q 3 , so λ must be pure imaginary. Hence, no instability occurs. Unlike κ ¼ 1=2 and 1, all the relations among the parameters derived from the condition q ¼ 0, do not give transition curves on the (κ − ϵ) plane. Now, let us go into the cubic equation with respect to q ¼ λ 2 (96) in more detail. We refer to the discriminant of the cubic equation (96) as D leading . Equation (96) says D leading ≥ 0 as long as the determinant of coefficient matrix is evaluated at leading order. In the case of D leading > 0, even if higher order contributions are considered, it still remain D higher > 0. However, if D leading ¼ 0, the sign of discriminant D higher can be minus due to higher order contributions. This means that two q out of three are complex, and that the solution (90) is always unstable.
Before considering higher order contributions, we derive the condition that the equation (96) has multiple roots, i.e., D leading ¼ 0. There are three possibilities: It is expected that instability aroundκ ¼ 3=4 is caused by the coupling of the axion and the photon (k) through the magnetic field, so q 1 ¼ q 2 may be meaningful. Solving the equation we obtain the relation among parameters, In order to clarify the origin of instability, we need to proceed to the next order. Please refer the reader to the Appendix for concrete expressions. The determinant at next leading order det½R 3=4 ðq; χ; β; ϵÞ next leading is also cubic equation. Up to the next leading order, we have det½R 3=4 ðq; χ; β; ϵÞ up to next leading ≡ det½R 3=4 ðq; χ; β; ϵÞ leading þ det½R 3=4 ðq; χ; β; ϵÞ next leading ≡ C 3 ðχ; β; ϵÞq 3 þ C 2 ðχ; β; ϵÞq 2 þ C 1 ðχ; β; ϵÞq where we labeled the coefficients of each orders of q as C 3 , C 2 , C 1 , C 0 . The discriminant of cubic equation (104)  Let us find a correction term X to the relation among parameters, On the curve that satisfies this relationship (106), we expect that the cubic equation (104) will have multiple roots, ½q þ g 1 ðβ; ϵÞ½q þ g 2 ðβ; ϵÞ 2 ¼ 0: In the Eq. (106), a higher order contribution X is incorporated into the leading order relation (102). The correction X is chosen so that the leading order of D up to next leading vanish: Thus, we can get correction terms, Up to the next leading order, the particular relationships among parameters are given by Remarkably, the original curve (102) splits into two curves (110). In the region which intervene between (110), the all order of discriminant D up to next leading is negative, i.e., D up to next leading < 0. Hence, in the region between the two curves (110), the cubic equation (104) can be rewritten as follows: ½q þ g 1 ðβ; ϵÞ½q þ h 1 ðβ; ϵÞ þ ih 2 ðβ; ϵÞ × ½q þ h 1 ðβ; ϵÞ − ih 2 ðβ; ϵÞ ¼ 0; ð111Þ and the criterion for stability (71) is not satisfied. Thus, we have found that two curves (110) is nothing but the transition curves forκ ¼ 3=4.

V. DISCUSSION
In the case of the conventional Mathieu equation, bifurcation points are located atκ ¼ n=2. This bifurcation point can be also rewritten as follows: where m a is axion mass, and k is the wave number. Diagrammatically, Eq. (112) for n ¼ 1 can be interpreted as in Fig. 4. Namely, the parametric resonance is nothing but a coherent decay of axions into photons. As described in the Sec. III C and Sec. IV C, in the situation where axion and magnetic field coexist in the background, a new bifurcation point ffiffiffiffiffiffiffiffiffiffiffiffi ffi 1 þκ 2 p þκ ¼ n arises. This bifurcation point can be also rewritten as follows: The case for n ¼ 2 is illustrated in Fig. 5. In this case, an axion and a photon are generated through the coherent decay of axions and photons in the background. Finally, let us consider what system needs to be arranged to give rise to instability seen in Sec. III C and Sec. IV C.
Appropriate numerical values depend on the wavelength L, The dimensionless parameter κ determine a relation between the axion mass and the wavelength of electromagnetic waves as follows: The parameter β characterize the strength of magnetic fields. In the case of a neutron star, and trying to detect the signal with microwaves (L ¼ 10 2 cm), the value of β is given by The axion dark matter-photon conversion could be effective. On the other hand, ϵ has a extremely small value in the case of (115) and (116). If you assume ultralight axions, ϵ has a suitable value: The de Broglie wavelength L dB where spatial variation of axion can be neglected is as follows: where we took a typical velocity in the galaxy v ∼ 10 −3 .
Since the coherence length is sufficiently long, we can expect parametric amplification of electromagnetic waves with the wavelength 10 18 cm ∼ 1 pc. Here, we would like to comment on the coherence of the axion dark matter. We assumed coherence of the axion dark matter, and use the values 0.3 GeV=cm 3 and v ∼ 10 −3 as a rough parameter estimate. However, for the more precise analysis, we need to compare the bandwidth of instability with velocity dispersion as discussed in [31]. Now, the question is whether both the conversion and the resonance can be important at the same time. In the case of radio waves, three parameters included in basic equations (43)-(45) are given as follows: : ð121Þ A strong magnetic field B 0 ∼ 10 7 G can be realized with a white dwarf. However, energy density ρ needs 10 8 times as much as the average density of dark matter near the solar system. Next, in the case of an ultralight axion, three parameters included in basic equations (43)-(45) are given as follows: ð124Þ It might be difficult to find the astrophysical situation with the strength of magnetic fields, B 0 ¼ 10 −1 G and to detect electromagnetic waves L ¼ 10 18 cm ∼ 1 pc. Devising a smart way, we may be able to realize the situation where both the axion dark matter-photon conversion and the resonance are relevant in the laboratory.

VI. CONCLUSION
We studied the stability of axion dark matter-photon conversion numerically and analytically. Since the axion field is coupled with the electromagnetic field, axions can be converted into photons and vice versa. On the other hand, axion is one of the candidates for dark matter, and photon propagating in axion dark matter obeys the Mathieu equations. Therefore, it is important to understand the behavior of the system where the axion dark matter and the magnetic field coexist. First, we derived basic equations describing axion dark matter-photon conversion. Then, we found the instability band by numerical calculations. Remarkably, we found the bands different from those in the conventional Mathieu equation. Moreover, we found the shift of bifurcation point due to the magnetic fields. More importantly, we confirmed numerical findings by using the analytical method. In the course of the analysis, we found that the condition for the new boundary curves between the stability and the instability requires different method from that of the conventional instability condition. Finally, we gave graphical interpretation to the instability condition, and comment on a possible physical application.