Photon directional profile from stimulated decay of axion clouds with arbitrary momentum distributions

We model clusters of axions with spherically symmetric spacial but arbitrary momentum distributions and study the directional profile of photos produced in their evolution through spontaneous and stimulated decay of axions via the process a -->gamma + gamma. Several specific examples are presented.


INTRODUCTION
The strong interactions conserve CP invariance to a high degree of accuracy, but the Lagrangian for QCD L = 1 4 F 2 + θFF + m Fψ ψ appears to have two allowed sources of CP violation. One is the instant term θFF , and the other is the fermion mass matrix m F . The full contribution can be writtenθ = θ + arg[det(m F )] and compared with the best experimental limit from the neutron electric dipole moment [1] isθ ≤ 10 −10 radians. To avoid such a fine-tuning, Peccei and Quinn [2,3] proposed a solution to this so called strong CP problem whereθ is promoted to a field a, the axions, who's potential when minimized relaxes the field (and affectivelyθ) to zero. While this axion has a KeV mass [4,5] and interacts too strongly to be allowed by experiment, other "invisible axion" variants are lighter and allowed [6].
Invisible axions are weakly coupled light pseudo scalars that are a perfect dark matter candidate. Current experiments limit the axion to the mass range 10 −3 to 10 −5 eV. They are nonrelativistic if produced at the electroweak (EW) phase transition, and are unlike other particles in this mass range like neutrinos which were produced and thermalized in the Big Bang. The lower limit of the axion mass range is due to the fact that if they were lighter they would be so abundant that they would over close the Universe. Details can be found in [7][8][9]. Early reviews of axions and the strong CP problem and axions in cosmology and astrophysics are [10][11][12][13].
The physics of the axion is similar to the π 0 in that they are both pseudo scalars and can decay to two photons through the one loop triangle anomaly diagram related to breaking the axial U(1) global symmetry of massless QCD. The coupling of the axion field a(x) to photons is where α is the fine structure constant, K is an O(1) model dependent constant and F a is the axion decay constant which can be related to the pion mass and decay constant through m a F a ≈ m π F π where m a is the axion mass. The axion lifetime scales with the π 0 lifetime τ a ≈ m π m a 5 τ π 0 ∼ 10 8 eV m a

5
Gyr so axions are very long lived compared to the age of the Universe. Since axions are long lived and copiously produced at rest in the EW phase transition, they are a candidate for cold dark matter and can be the source of early universe density perturbation that grow in a way consistent with observations of the cosmic microwave background (CMB), which is unlike the results for free streaming hot dark matter perturbation that do not grow at small length scales.
Axions in vacuum can undergo spontaneous decay a → γγ, and if their density is high enough, stimulated decay can cause lasing to commence and rapidly deplete the axion number density, while the photon number grows exponentially. Sufficiently dense regions of axions can release enough monochromatic electromagnetic energy for possible detection [14,15].
The case of lasing axions for a spherical symmetric distribution, both in momentum and coordinate space has been considered in [14][15][16]. However, resent studies of axions as dark matter suggests they may be in several other phase space configurations. E.g., they may fall along caustics in galaxies [17,18], or they may be produced by superradiance, a Penrose type of process, around black holes [19] and fill a hydrogen-like orbit with quantum numbers (n, l, m) = (2, 1, 1). (Since axions are bosons there is no fundamental limit on their numbers in this state.) For these reason we have decided to begin an exploration of the general case of lasing for arbitrary axion distributions in terms of a spherical harmonic expansion.
The axion mass range can be expanded somewhat by fine tuning, but this begins to defeat the purpose of introducing the axions in the first place to avoid fine tuningθ. The range can also be loosened by considering axion-like particles not necessarily designed to solve the strong CP problem. The analysis we present here can be applied to any of these cases as long as the axion or axion-like particles have a dominant decay mode to two photons.
In this paper we focus on arbitrary momentum distributions while keeping the coordinate space distribution spherically symmetric. In particular we write a general expansion of momentum modes in spherical harmonics. Many of the algebraic details are relagated to the appendices. Once we have derived the general form we then consider a few specific examples of physical interest. In future work we plan to also relax the coordinate space spherical symmetry requirement.

SETUP AND PREPARATION
To fix notation and set the stage for our analysis we first review the setup of the sphericrally symmetric phase space axion cluster model [14,16]. Axions of mass density ρ a are produced nonthermally during the QCD phase transition with a cosmological density parameter today [7][8][9] of where ρ c is the critical density. Hence axions are a dominant dark matter (DM) candidate if m a ∼ 10 −5 eV.
Our purpose here is to give detailed calculations of the stimulated emission rate for axion clusters. There are many possible sources of initial photons of frequency in the ∼ 10 −5 eV range, therefore spontaneous decays are not required to start the lasing process, although spontaneous decay of axions is a lower bound on axion cluster luminosity. We neglect all normal matter in the cluster, hence we neglect the attenuation of photons by nonaxionic matter, as well as any other affects due to the local environment.
For any species of particles with occupation number f (p, r, t), the particle number density is n(r, t) = d 3 p 8π 3 f (p, r, t), and the total particle number of these particles in volume V is Angular momentum conservation requires that a pair of photons emitted by decay of a spin zero particle (scalar or pseudoscalar) have the same helicity. The rate of change of the photon number density of helicity λ = ±1 within an axion cluster, due to the process a ↔ γγ can be written as a Lorentz invariant phase space integral where f a = f a ( p), f iλ = f iλ ( k i ) and M = M(a → γ(+)γ(+)) = M(a → γ(−)γ(−)). In detail, the phase space integral is where p is the axion momentum and k 1 and k 2 are the momenta of the photons.
Upon defining k = k 0 2 and using the above results, we arrive at the rate equation for the photon occupation number In the case of spherical symmetry we set f ( q) = f (| q|) = f (q) and the intergral simplifies to give In [14,16] a specific model with spherical symmetry was choosen where the initial axions were contained in a ball of radius R, with a maximum momentum value of p max ≈ m a β in the nonrelativistic case. I.e., the initial axion occupation number was chosen to be Where C a is a constant that can be written in terms of axion mass and initial number density. Details of this model were worked out in [14,16] and will be recovered below when we take the spherical symmetric limit of the general case.
The direction of axion momentum can be characterized by the infintestimal solid momentum angle Ω p = (θ p , φ p ). f a (p, r, Ω p , t) and n a (r, Ω p , t) are the axion occupation number and axion number density, respectively. A integral of f a (p, r, Ω p , t) over p or (p, Ω p ) gives number density n a (r, Ω p , t) or n a (r, t).
The photons produced in the decay of axions are initially contained in the ball of radius R, and in a momentum spherical shell of inner and outer radius k − = maγ 2 (1 − β) and k + = maγ 2 (1 + β), respectively. f λ (k, r, θ, t) and n λ (r, θ, t) are the photon occupation number and photon number density, respectively, of helicity λ = ±1. We assume that the number density of each helicity state is the same, so the total photon number density n γ can be written as n γ (r, θ, t) = n + (r, θ, t) + n − (r, θ, t) = 2n λ (r, θ, t).
Since the momentum distribution of axions is assumed to be spherically symmetric, there is an equal chance for an axion to be moving in any direction. Likewise the lasing process is equally likely to commence in any direction. In [14,16] this model was used to give lasing bounds and find the stimulated decay rate equations for axions and photons. Numerical results were also given for the evolution of such an axion cluster.
Here we want to investigate how the momenta of photons are distributed if a general non-spherically symmetric distribution of the axion momenta is specified. We do this by generalizing the model of [14,16] and proceeding without making this assumption of spherical symmetry in momentum space. We proceed by expanding in spherical harmonics and then integrating over all angles of photon momentum.

CUPATION NUMBER
Once we relax spherical symmetry we must take angular dependencies into account.
There are several sets of angles to consider, those for both axion and photon position and momentum, hence a potential total of eight angles on which our results may depend. Let us begin with the axions.
No spherical symmetry in momentum space means f a ( p) = f a (p, θ p , φ p ) is a general function of angles θ p and φ p . Since f a ( p) is the occupation number of axions, it has to be integrable.

Harmonic Expansion
We can write f a ( p) as the square of a square integrable function formed in a complex spherical harmonics expansion, where the Y m l s are complex spherical harmonics, and the C(l, l ′ , l ′′ |m, m ′ m ′′ )s are Clebsch-Gorden coefficients.
By regrouping and renaming coefficients, f a ( p) can be written directly as a real spherical harmonic expansion, Similarly, photon occupation numbers can also be written as The following calculations do not put any restrictions on the coefficients a lm 's or b lm 's, neither do the calculations require these coefficients to be positive or negative. But in real physical world, only positive occupation numbers are allowed. In short, these coefficients can be of any value, but we should scrutinize the resulting occupation numbers so that they are positive and have real world meaning. For example, occupation number f a ( p) = a 10 (p, t)Y 10 (Ω p ) can not describe any real scenario since it is negative in half of phase space, but occupation number f a ( p) = a 00 (p, t)Y 00 (Ω p )+a 10 (p, t)Y 10 (Ω p ) may be allowed since it can be positive everywhere by adjusting the coefficients a 00 and a 10 . Nevertheless the following calculations are applicable to both of these occupation numbers, whether they have real physical meaning or not.
Integration over k 1 , Ω k 1 We start our analysis from the general rate equation given in [16], which describes the evolution process of spontaneous and stimulated decay of axions to photons, and the back reaction of photons to axions. 2k where Γ a = 1/τ a is the axion decay width. Writing out the differential d 3 k 1 and the δ- canceling common factors and substituting from the harmonic expansions, we find We can do the integration over θ k 1 and φ k 1 most efficiently by changing from Ω k 1 to Ω p− k .
Upon use the fact that for photons k 0 1 = k 1 and k 0 = k along with the identity Integration over φ p : preparation Substituting d 3 p = p 2 dp dΩ p = (p 0 ) 2 − m 2 a p 0 dp 0 dΩ p and various relations from Appendix A into equation (5), we find Canceling common factors and writing the equation in a form that is suitable for doing the two angular integrations, we have The momentum space angular integration requires some algebra which can be found in Appendix C. Adopting θ p1 notation from eq. eq. (22) in Appendix C, eq. (6) become Evolution of individual components The δ-function in (4) imposes conservation of 4-momentum, since the integration has been carried out, k and p 0 in (23) of Appendix C are independent. Hence, the k in the denominator of (23) can be taken outside of the integral, Upon writing out the Ω k dependence of RHS of the equation above explicitly, and identifying the coefficients of Y lm (Ω k ), we obtain differential an equations for components b lm (k, t), Replacing p 0 with k 1 + k, also using the following substitution,  Axion and photon occupation numbers depend on position. This fact was not included in previous discussion, since it is not involved in the dynamical process. Here we will assume that spherical symmetry is present in position space and axions sit inside a ball of radius R, where the occupation number components can be written Hence, the axion number density is where the components n a lm (t) and a lm (t) are related through n a lm (t) = Photo occupation number and number density can be treated in a similar fashion. We and n λ (t, r, The components n λ lm (t) and b lm (t) are related through In Appendix C, we replaced all the occupation number components in (8) with the corresponding number density components and integrated over the momentum space, which gives us the evolution equations for the individual components of number density of photon, normal axion, and sterile axion.
The K 0 l ′ are defined to be constant coefficients describing spontaneous and half of stimulated decay. K 01 l ′ l ′′ m ′ are constant coefficients describing the other half of stimulated decay. The B 1 l ′ are constant coefficients describing back reaction of photons, which when necessary, can be split into N 1 l ′ and S 1 l ′ the part of back reactions that produce normal axions and sterile axions respectively, and hence B 1 The choice Y 00 corresponds to the momentum distribution of axion with no preferred direction. The number density of axion for such a momentum distribution is Thus, for any l = 0, m = 0, the number density component is zero, n a lm (t) = 0 (lm = 00) .
This means that the evolution equations (14) for these components reduce to There are two cases regarding the solution of the equation above.
Since the photon field is fixed in the Y 00 momentum distribution, the exact form of number density of normal axion, photon and sterile axion can be obtained by solving the evolution equations (13), (14) and (15). In this first case we find dn as These equations are the same as (34 ′ ), (37 ′ ) and (38 ′ ) of [16], when we set The second case is This solution to eq. (16,17,18) means that the effect of half of stimulated decay exactly cancels the effect of photon annihilation back into normal axion. This also expresses the l = 0, m = 0 component of axion number density as a sum of components of photon number density, Evolution equation (14) for l = 0, m = 0 component then becomes Evolution equations (13) for any l = 0, m = 0 components of photon number density in this second case is Photon annihilation back into sterile axion and surface loss are the only mechanisms that contribute to n γ lm (t)(l = 0, m = 0) modes. There is no source of decay providing photons to l = 0, m = 0 number density components. Therefore n γ lm (t)(l = 0, m = 0) modes are expected to vanish quickly. The majority of photons would be in n γ 00 (t) component.
In both cases, either n γ lm (t) = 0 or 3 if the axions are locked in Y 00 momentum state, then so would be the photons (or at least predominantly in the second case).

Y 20 momentum distribution
In this configuration the direction of momentum of axions is preferentially parallel to the polar axis. It describes the movement of axions between northern and southern hemispheres n a (t, r, The evolution equation (14) for n a 20 (t) becomes dn a At the same time, except for n a 00 (t) and n a 20 (t), evolution equation (14) becomes = 0 (lm = 00, 20) .
This shows that there are only two nonzero components for photon number density, n γ 00 (t) and n γ 20 (t).
Since the n a 00 (t) component of axion number density is 0, evolution equation (14) reduces to This can be used to simplify the photon evolution equation (13) for component n γ 00 (t), The n γ 00 (t) component is expected to vanish quickly since the only contributions to this mode are due to back reaction to sterile axions and surface losses. Meanwhile, the photon component n γ 20 (t) evolves according to equation (13), If axions were locked in a Y 20 momentum state, the majority of photons would be expected to be in this state, which means more photons would travel in the direction that is to some extent parallel to the polar axis. However, this is not an acceptable physical number density for particles because Y 20 is negative in some regions, but in the appropriate combinations with Y 00 , the total number density can be positive.
This configuration has the direction of momentum of axions preferentially parallel to the equatorial plane, and is a good approximation to the case of superradiant axions near a Kerr black hole [19]. It describes axions rotating around a polar axis. The axion number density has two nonzero components, n a 00 and n a 20 .
n a (t, r, To maintain the sin 2 θ distribution, there is a relation between these two nonzero components, Evolution equation (14) for n a 20 (t) becomes dn a while for components other than n a 00 (t) and n a 20 (t), evolution equation (14) becomes = 0 (lm = 00, 20) .
This shows that there are only two nonzero components for photon number density, n γ 00 (t) and n γ 20 (t).
In Appendix D, we find that the quantity n γ 00 (t)+ √ 5n γ 20 (t) evolves based on the following equation, Therefore n γ 00 (t) + √ 5n γ 20 (t) = 0 is a possible solution to this differential equation although it is not the uniquely solution, which indicates that axions of sin 2 θ momentem distribution can still generate photons of sin 2 θ momentem distribution. However, since only photon annihilation back into sterile axion and surface loss contribute to the change of the quantity n γ 00 (t) + √ 5n γ 20 (t), we expect, So approximately, decay from sin 2 θ momentum distribution of axions results in a similar directional profile of photons. Here we collect some useful kinematic relations needed in the text.
Choose the z-axis to align with photon momentum k, then where cos θ p0 is the root of the equation g(cos θ p ) − (p 0 − k) = 0. Clearly The derivative at cos θ p0 can also be calculated, Substituting the root and the derivative, and δ[| p − k| − (p 0 − k)] is converted to a δ function with respect to θ p , Then the common way to obtain various trigonometric quantities about p is through relations for the angles given in the appendix Similarly, trigonometric quantities about p − k can be obtained via following relations, Thus φ p− k and φ p are equal, Then the third φ p integration from equation (6) can be carried out,

APPENDIX C: MOMENTUM SPACE ANGULAR INTEGRATION
Integration over φ p : Y lm (Ω p ) and Y lm (Ω p− k ) The φ p integration of the first term in braces in eq. (6) can be calculated directly using eq. (1). Due to the factor e imφ contained in Y lm , the integration is 0 unless m = 0.
As stated before, k is chosen as the z-axis.
Again substituting the form of f a from eq. (1), we see the φ p integration in the second term of eq. (6) has to be performed on Y l ′ m ′ (Ω p )Y lm (Ω p− k ), The calculation is easier if we convert from real back to the complex spherical harmonics. But the relation between real and complex spherical harmonics depends on the signs of m ′ and m. Using the relations among spherical harmonics found in Appendix B facilitates the integration and we find, The third integral is similar to the first.
Collecting terms after all φ p integrations eq. (6) becomes Apart from writing down step function Θ(R − r), integrations are carried over (8) term by term in the following calculation, The integration on the RHS of (24) is over k and k 1 . The following calculations omit writing down common factor maΓa (2π) 3 . The first term on the RHS of (24) is where K 0 l ′ are defined to be constant coefficients describing spontaneous and half of stimulated decay, and are related to associated Legendre polynomials P 0 l ′ .
The second term on the RHS of (8) can be treated similarly, where K 01 l ′ l ′′ m ′ are constant coefficients describing the other half of stimulated decay and are related to associated Legendre polynomials P m ′ l ′ , P m ′ l ′′ .
The third term on the RHS of (24) is simplified by using newly defined coefficients K 0 l ′ . dk m 2 a 4k The last term on the RHS of (24) is dk m 2 a 4k where B 1 l ′ are constant coefficients describing back reaction of photons and are related to associated Legendre polynomials P 0 l ′ .
However, back reaction can produce sterile axions with energy higher than m a γ which need to be excluded from axion number counting since we focus on nonrelativistic axions. Dividing the integral inteval [ m 2 a 4k , k + ] into two parts gives two sets of constant coefficients N 1 l ′ and S 1 l ′ describing back reactions produce normal axions and sterile axions, respectively, where This means that the quantity n γ 00 (t) + √ 5n γ 20 (t) evolves based on the following equation,