Gravitationally bound axions and how one can discover them
GGravitationally bound axions and how one can discover them
Xunyu Liang ∗ and Ariel Zhitnitsky † Department of Physics & Astronomy, University of British Columbia, Vancouver, Canada
As recently advocated in [1], there is a fundamentally new mechanism for the axion productionin the Sun and Earth. However, the role of very slow axions in previous studies were neglectedbecause of its negligible contribution to the total axion production by this new mechanism. In thepresent work we specifically focus on analysis of the non-relativistic axions which will be trappedby the Sun and Earth due to the gravitational forces. The corresponding emission rate of these lowenergy axions (below the escape velocity) is very tiny. However, these axions will be accumulatedby the Sun and Earth during their life-times, i.e. 4.5 billion of years, which greatly enhances thediscovery potential. The computations are based on the so-called Axion Quark Nugget (AQN) DarkMatter Model. This model was originally invented as a natural explanation of the observed ratioΩ dark ∼ Ω visible when the DM and visible matter densities assume the same order of magnitudevalues, irrespectively to the axion mass m a or initial misalignment angle θ . This model, withoutadjustment of any parameters, gives a very reasonable intensity of the extreme UV (EUV) radiationfrom the solar corona as a result of the AQN annihilation events with the solar material. This extraenergy released in corona represents a resolution, within AQN framework, a long standing puzzleknown in the literature as the “solar corona heating mystery”. The same annihilation events alsoproduce the axions. The flux of these axions is unambiguously fixed in this model and expressed interms of the EUV luminosity from solar corona. We make few comments on the potential discoveryof these gravitationally bound axions. ∗ [email protected] † [email protected] a r X i v : . [ h e p - ph ] J a n I. INTRODUCTION
The Peccei-Quinn (PQ) mechanism accompanied by the axions remains the most compelling resolution of the strongCP problem, see original papers [2–4] and recent reviews [5–13] on the subject. We refer to the review articles for thediscussions and analysis on the recent activities in the field of the axion searches by a numerous number of differentgroups using very different instruments.For the purposes of the present work it is sufficient to mention that the conventional dark matter galactic axions areproduced due to the misalignment mechanism [14] when the cosmological field θ ( t ) oscillates and emits cold axionsbefore it settles down at its final destination θ final = 0. Another mechanism is due to the decay of the topologicalobjects [15–20]. There is a number of uncertainties and remaining discrepancies in the corresponding estimates. Weshall not comment on these subtleties by referring to the original papers [15–20]. It is important that in both casesthe produced axions are non-relativistic particles with typical v axion /c ∼ − , and their contribution to the darkmatter density scales as Ω axion ∼ m − / a . This scaling unambiguously implies that the axion mass must be fine-tuned m a (cid:39) − eV to saturate the DM density today, see footnote 1, while larger axion mass will contribute very little toΩ DM . The cavity type experiments have a potential to discover these non-relativistic axions.Axions can be also produced as a result of the Primakoff effect in a stellar plasma at high temperature [21]. Theseaxions are ultra-relativistic as the typical average energy of the axions emitted by the Sun is (cid:104) E (cid:105) = 4 . v AQNaxion (cid:39) . c . These features should be contrasted withconventional galactic non-relativistic axions v axion ∼ − c and solar ultra-relativistic axions with typical energies (cid:104) E (cid:105) = 4 . v (cid:28) c . Thispart of the spectrum represents very tiny portion of the produced axions. Therefore, it had been ignored in theoriginal studies [1]. However, the upgraded CAST instrument will be highly sensitive to the low energy part of thespectrum. Therefore, it is highly desirable to develop a new computational technique which would allow to carry outthe computations in the region of small velocities v (cid:28) c .Furthermore, the low-energy axions produced in the Sun might be trapped by strong gravitational force such that v ≤ v trapped will be trapped by the Sun indefinitely. The v trapped is numerically the same as the free fall velocity, v trapped c = (cid:115) GM (cid:12) R (cid:12) (cid:39) · − . (1)While the portion of these low energy axions is tiny as we shall estimate below, these trapped axions may play animportant role in physics as they will be accumulated around the Sun during entire life time of the solar system,i.e. around 4.5 billion years. The effects related to the trapped axions are not new, and discussed previously in theliterature [23]. The goal here is to present some numerical estimates for our specific AQN model when the axionswhich are produced as a result of the annihilation events in the solar atmosphere and will be indefinitely bounded tothe Sun.Therefore, the main goal of the present studies is to develop a new technique to generalize the results of ref. [1]to perform the self-consistent computations of the axion spectrum in the regime when the axion velocities are small v (cid:28) c . The corresponding generalization of the results [1] requires to abandon the so-called “thin wall approximation”and develop some new technical tools which proper describe the regime with v (cid:28) c .We should emphasize that the present work is a natural continuation of the previous studies [1]. Therefore, to avoidrepetition we refer the readers to that paper with detail discussions of the AQN framework itself, its motivation, itsconsequences and predictions. The only comment we would like to make here is as follows. The AQN framework isconsistent with all known astrophysical, cosmological, satellite and ground based constraints. In fact, in a number ofcases some observables become very close to present day constraints. Furthermore, in few cases the predictions of themodel may explain a number of the long standing mysteries as overviewed in [1]. According to the most recent computations presented in ref.[20], the axion contribution to Ω DM as a result of decay of the topologicalobjects can saturate the observed DM density today if the axion mass is in the range m a = (2 . ± . − eV, while the earlierestimates suggest that the saturation occurs at a larger axion mass. One should also emphasize that the computations [15–20] havebeen performed with assumption that PQ symmetry was broken after inflation. The paper is organized as follows. In Section II we develop a new technique which allows to generalize thesecomputations for low energy portion of the spectrum when v (cid:28) c . We use the corresponding results in Section IIIto discuss the physics of the trapped axions and we highlight the basic ideas how to discover them. We conclude inSection IV with few thoughts on the future developments. II. AXIONS FROM AQNS: INTENSITY AND THE SPECTRUM
The AQN model was invented long ago [24] (though a specific formation mechanism of the nuggets was developedin much more recent papers [25–27]) as a natural explanation of the observed ratio Ω dark ∼ Ω visible . In context of thepresent work the argument supporting the AQN model goes as follows. It has been known for quite some time thatthe total intensity of the observed EUV and soft x-ray radiation (averaged over time) can be estimated as follows, L (cid:12) (from Corona) ∼ · GeVs ∼ · ergs , (2)which represents (since 1939) the renowned “the solar corona heating puzzle”. The observations (2) imply that thecorona has the temperature T (cid:39) K which is 100 times hotter than the surface temperature of the Sun, andconventional astrophysical sources fail to explain the EUV and soft x ray radiation from corona.It turns out that if one estimates the extra energy being produced within the AQN dark matter scenario oneobtains the total extra energy ∼ erg / s which precisely reproduces (2) for the observed EUV and soft x-rayintensities [28]. The full scale Monte Carlo simulations [29] support this estimate. One should add that the estimate ∼ erg / s for extra energy is derived exclusively in terms of known dark matter density ρ DM ∼ . − anddark matter velocity v DM ∼ − c surrounding the sun without adjusting any parameters of the model. We interpretthis “numerical coincidence” as an additional hint supporting the AQN model. Our original remark relevant forthe present work is that if one accepts the explanation that the solar corona heating puzzle is resolved within AQNscenario than the axion flux will be unambiguously fixed in terms of the EUV observed luminosity (2) as the axionfield represents the crucial element in the AQN construction.We start our presentation with subsection II A where we highlight the basic results from ref. [1] by providing aself-contained text for the convenience of the readers. In next subsections II B and II C we explain the computationalframework and present the results of the computations, referring to Appendix A for the technical details. A. Intensity
The axions play a key role in construction of the AQNs as they provide an additional pressure to stabilize thenuggets. The corresponding axion contribution into the total nugget’s energy density has been computed in [27].Depending on parameters the axion’s contribution to the nugget’s mass represents about 1/3 of the total mass. Itcan be translated in terms of the axion luminosity from the Sun as follows [1]: L (cid:12) (axion) (cid:39) . · · ergs . (3)The corresponding axion flux measured on Earth can be computed as follows [1]Φ(solar axions) ∼ L (cid:12) (axion) π (cid:104) E a (cid:105) D (cid:12) ∼ . · s (cid:18) − eV m a (cid:19) , D (cid:12) (cid:39) · km , (4)where we assume that the axion’s energy when the antinuggets get annihilated is slightly relativistic E a (cid:39) . m a , butnever becomes very relativistic. The corresponding energy flux is [1] m a Φ(solar axions) ∼ · eVcm s . (5)These estimates should be compared with conventional cold dark matter galactic axion contribution assuming theaxions saturate the observed DM density: m a Φ(galactic axions) ∼ ρ DM · v DM (cid:39) . v DM (cid:39) eV cm s . (6)Similar estimates can be also carried out for Earth. In this case as explained in [1] the observations of the E & M showers due to the nuggets entering the Earth’s atmosphere (before hitting the Earth’s surface) require very large areadetectors. The nuggets will continue to radiate E & M energy in the deep underground. However, this radiation byobvious reasons cannot be recovered and observed. At the same time the observation of the axions (which have beenproduced as a result of the annihilation events in the very deep underground) is possible, and in fact very promising.Indeed, the corresponding axion flux can be estimated as follows [1] m a Φ(Earth axions) ∼ · (cid:18) ∆ BB (cid:19) eVcm s , (7)where ∆ B/B is the portion of the AQNs which get annihilated in the Earth’s interior. Interestingly, the axion flux(7) which is generated due to the AQN annihilation events is much larger than the flux (5) generated due to the AQNannihilation events in the solar corona and measured on Earth. At the same time, the axion flux (7) is the same orderof magnitude as the conventional cold dark matter galactic axion contribution (6). This is because the parameter∆
B/B ∼ v ∼ . c , in contrast with conventionalgalactic isotropic axions with v ∼ − c . Therefore, these two distinct contributions can be easily discriminated. B. Spectral properties. General Comments
The basic idea of the computation of the spectrum is as follows. Consider an AQN loosing its mass due to theannihilation with surrounding material, while that the axion portion to the energy remains the same, as it is notlinked to the annihilation processes. One should comment here that the axion domain wall in the equilibrium doesnot emit any axions as a result of pure kinematical constraint: the domain wall axions are off-shell axions in theequilibrium. The time dependent perturbation obviously changes this equilibrium configuration. In other words, theconfiguration becomes unstable because the total energy of the system is no longer at its minimum. To retrieve itsground state, an AQN will therefore intend to lower its domain wall mass by radiating the axions. To summarize:the emission of axions is an inevitable consequence during the annihilation of antinuggets in simply for the reason tomaintain the AQN stability.Now, we want to identify a precise mechanism which produces the on-shell freely propagating axions emitted by theaxion domain wall. In this section we overview the basic idea of the computational technique to be used. To addressthis question, we consider the general form of a domain wall solution: φ ( R ) = φ w ( R ) + χ (8)where R is the radius of the AQN, φ w is the classical solution of the domain wall, while χ describes the excitations dueto the time dependent perturbation. We should note that, φ w is clearly off-shell classical solution, while χ describesthe on-shell propagating axions. Thus, whenever the domain wall is excited, namely χ (cid:54) = 0, freely propagating axionswill be produced and emitted by the excitation modes.Few comments are in order before we proceed in subsection II C with description of the technical details andcorresponding results. First, if the domain wall can be considered to be infinitely large in xy direction such thatthe profile function depends on a single variable z the computations can be carried out easily as the classical profilefunction φ w ( z ) is known exactly. This is precisely the procedure which has been adopted in previous paper [1]. Thecorresponding technique is justified when a typical size L x ∼ L y (cid:29) m − a along x, y is much larger than the width ofthe domain wall of order m − a . If the wave length of the emitted axion is small, i.e. λ a ∼ m − a the axions cannotcarry any information about the finite size of the system and the approximation is marginally justified ( λ a standsfor the de Brogile wavelength of the emitted axion). This is precisely the approximation, the so- called “thin wallapproximation” which has been adopted in computations [1]. This approximate treatment is marginally justifiedfor relativistic axions with v ∼ c , and we expect that accounting for the finite size of the system cannot drasticallychange the results in the relativistic domain v ∼ c . This will be explicitly confirmed below by present computationsaccounting for finite size of the system.Secondly, it is quiet obvious that the “thin wall approximation” is badly broken for non-relativistic axions with v (cid:28) c when λ a (cid:29) m − a and a new technique must be developed for proper analysis. The basic idea of computation accountingfor finite size of the system R goes as follows. Suppose an AQN is traveling in vacuum where no annihilation eventtakes place, we expect the solution stays in its ground state φ ( R ) = φ w ( R ) which corresponds to the minimumenergy state. Since there is no excitation (i.e. χ = 0), no free axion can be produced. However, the scenariodrastically changes when some baryon charge from the AQN get annihilated. Due to these annihilation processes,the AQN starts to loose a small amount of its mass, and consequentially its size shrinks from R to a slightly smallerradius R new = R − ∆ R . Note that its quantum state φ ( R ) = φ w ( R ) is no longer the ground state, becausea lower energy state φ w ( R new ) becomes available. Then, we may write the current state of the domain wall as φ ( R ) = φ w ( R new ) + φ (cid:48) w ( R new )∆ R , so the domain wall now has a nonzero exciting mode χ = φ (cid:48) w ( R new )∆ R and freeaxions can be produced during oscillations of the domain wall. To reiterate: the annihilation of antinuggets withsurrounding matter forces the domain wall to oscillate. These domain wall oscillations generate excitation modeswhich ultimately lead to radiation of the propagating axions.Our last comment deals with terminology and notations. The results for the spectrum obtained using the “thin-wallapproximation” is coined as 1D spectrum. As we mentioned above it is marginally justified when λ a ∼ m − a , and itadmits mathematically exact treatment which was previously presented in [1]. In the present work we mostly dealswith 3D computations when a finite size of the system plays a key role, which is always the case for λ a (cid:29) m − a . Thepotential pitfall is that some technical simplifications are required to treat the 3D case. Consequentially, the obtainedresults might be sensitive to these technical simplifications. In order to characterize the sensitivity to our technicalsimplifications we introduce a tunable parameter δ which varies from 0 to 1, so δ will serve as a probe to test thesensitivity with respect to numerical simplifications. As we shall see below, the obtained results are not very sensitiveto the choice of δ . Therefore we conclude that our 3D results are robust and reliable.In what follows we will express the normalized spectrum as a function of the speed of emitted axion v a /c definedin the nugget’s frame, defined as follows ρ ( v a ) ≡ totaxions ddv a Φ axions ( v a ) , (cid:90) dv a ρ ( v a ) = 1 , (9)where Φ totaxions is the axion flux inserted to eq. (9) for normalization purposes. It assumes the magnitude Φ(solar axions)given by eq. (5) for the solar axions, and the value Φ(Earth axions) given by (7) for the axions emitted from theEarth’s core. C. Spectral properties. Results
We follow the procedure described above in subsection II B and present the axion field in time dependent backgroundas follows φ ( t, r ) = φ w,δ ( r − R ) + χ ( t, r ) (10)where φ w,δ ( r − R ) satisfies the classical equation of motion while χ ( t, r ) describes the time-dependent excitations. Asexact solution accounting for the finite size of the nugget is not known we parameterize different simplified solutionsby parameter δ . We consider parameter δ as a probe as explained above.The next step is to expand the action S [ φ ] by keeping the quadratic terms only, S [ φ ] = S [ φ w,δ ] + (cid:90) dt (cid:90) d x (cid:20)
12 ˙ χ − χL [ δ ] χ (cid:21) + O ( χ ) . (11)where L [ δ ] is the second order linear differential operator which depends on classical solution φ w,δ ( r − R ), and theparameter δ introduced here is a result of approximation to the true solution φ w ( r − R ), see Appendix A for thetechnical details. The next step, as usual, is to expand the fluctuations χ in terms of complete basis and compute thecoefficients a plm in this expansion. The result for the total radiated energy E rad is given by eq.(A24) from AppendixA. It can be presented it in the following form E rad = (cid:90) d x χ (cid:20) − ∂ ∂t + L [ δ ] (cid:21) χ = (cid:88) lm (cid:90) d p E a | a plm | = (cid:88) lm (cid:90) ∞ m a dE a · πp E a | a plm | , (12)where the coefficients a plm can be explicitly computed and are given by (A22). The expression for the radiated energy(12) allows us to compute the desired spectrum ρ ( v a ) defined by (9). The results of the computations are presentedon Fig. 1a with three different choices of parameter: δ = 0, 0.5, and 1 for physically realistic conditions, see AppendixA for the details. The low energy portion of the spectrum with 0 ≤ v/c ≤ .
01 is shown on Fig.1b.Few comments are in order. Parameter δ in our treatment of the problem was introduced as a probe to testour computational scheme which requires to compute all the modes in the background of the classical solutionparameterized by parameter δ . While the classical solution itself can be computed numerically, we need some analyticalform to proceed with computations of the modes. Parameter δ is precisely introduced in order to parametrize thisanalytical expression entering the differential operator L [ δ ]. As mentioned above, the parameter δ roughly variesfrom 0 to 1 in physically realistic circumstances. With the purpose of the test we performed the computations fordifferent values of δ shown in Figs. 1, where we also included the “unphysical value” for parameter δ = 8 exclusivelyfor illustrative purposes.One can explicitly see that the results for the spectrum are not very sensitive to parameter δ . As we discussbelow, the crucial factor ξ to be introduced in next section and which enters all final formulae is also not sensitive toparameter δ . To reiterate: the basic qualitative results are not very sensitive to choice of parameter 0 ≤ δ ≤ v a /c (cid:38) . v a /c ≤ .
01 which is the subject of the present work. Indeed, the 3D spectrumin this domain behaves as ρ ( v a ) ∼ v a as shown in Fig. 1b, in contrast with linear dependence in simplified treatmentin ref. [1]. This difference in behaviour at small v a /c (cid:28) ∼ d k in 3D case for λ a (cid:29) m − a . / c ρ ( v ) i n t en s i t y (a) spectrum for 0 ≤ v/c ≤ / c ρ ( v ) i n t en s i t y ( x - ) (b) zoom in portion of the spectrum with 0 ≤ v/c ≤ . FIG. 1: ρ ( v a , δ ) vs v a /c . Different values of δ are chosen respectively: 0 (blue), 0.5 (orange), 1 (green), and 8 (black). III. GRAVITATIONALLY TRAPPED AXIONS
In the previous section we computed the portion of the axions which have sufficiently low velocities (below escapevelocity) such that they will be orbiting the Sun as long as it exists. This portion of the non-relativistic axions isextremely tiny as shown on Fig. 1b. Nevertheless, the effect could be drastically enhanced as we discuss below dueto accumulation of these axions during entire life of the solar system, i.e. for ∼ v trapped ,defined as v trapped (cid:12) c = (cid:115) GM (cid:12) c R (cid:12) (cid:39) · − , v trapped ⊕ c = (cid:115) GM ⊕ c R ⊕ (cid:39) . · − , (13)such that all axions with v ≤ v trapped (cid:12) will be trapped by the Sun and the axions with v ≤ v trapped ⊕ will be trappedby the Earth. The effect of the trapped axions is not new, and discussed previously in the literature [23]. The goalhere is to present some numerical estimates for our specific AQN model when the axions which are produced as aresult of the annihilation events can be trapped in the solar atmosphere. These estimates will play a key role in ourdiscussions on the discovery potential of these axions. A. Solar corona background. Non-resonance case.
According to Fig. 1b these highly non-relativistic axions represent a very tiny portion of the produced axions.The energy which is accumulated in the solar atmosphere per unit time as a result of trapping these axions can beestimated as follows dE (cid:12) dt (trapped axions) (cid:39) . · · ξ · ergs (cid:39) · (cid:18) ξ − (cid:19) · ergs (14)where we used the expression (3) for the rate of the energy transfer to the axions. We also introduced the suppressionfactor ξ to account for the small fraction of the trapped axions with v ≤ v trapped . For numerical estimates in formula(14) we use suppression factor ξ ∼ − computed in previous section and presented on Fig. 1b.The axions (14) could not leave the system during entire life time of the Sun, i.e. 4.5 billion years (cid:39) s. Therefore,the total energy accumulated by the Sun and related to AQN annihilation events radiating the slow velocity axionscan be estimated as follows E (cid:12) (trapped axions) (cid:39) · (cid:18) ξ − (cid:19) · ergs · s (cid:39) (cid:18) ξ − (cid:19) erg . (15)This energy can be expressed in terms of extra solar mass ∆ M (cid:12) accumulated by the Sun and represented by thetrapped axions ∆ M (cid:12) (trapped axions) (cid:39) (cid:18) ξ − (cid:19) kg , (16)which of course represents a very tiny fraction of the solar mass M (cid:12) (cid:39) · kg.The energy (15) corresponds to the following total number of the axions accumulated by the Sun during its life-time: N axions (cid:12) ∼ E (cid:12) (trapped axions) m a c (cid:39) · (cid:18) ξ − (cid:19) · (cid:18) − eV m a (cid:19) . (17)If we assume that the majority of these axions are localized within 2 solar radius R (cid:12) , we arrive to the followingestimate for the average axion energy density inside this volume ρ axions (cid:12) ∼ E (cid:12) (trapped axions) π (2 R (cid:12) ) ∼ . · (cid:18) ξ − (cid:19) GeVcm , (18)which is 3 orders of magnitude larger than the present average dark matter density today ρ DM (cid:39) . GeVcm . One shouldcomment here that this enhancement of the DM density in the vicinity of the Sun obviously not in contradictionwith most precise observational upper limits on solar system (SS) -bound DM, which is normally expressed as ρ SS < · , see e.g.[30]. It is also interesting to note that some authors [31, 32] previously argued that the DM inthe SS might be greatly enhanced (on the level of 10 ) as a result of capturing of DM particles from the Galactichalo due to the 3 body interaction (the Sun, a planet and DM particle). Other authors [33, 34] estimated that theeffect of capturing is small. We refer to these original papers for the discussions and details. The only comment wewould like to make here is that the effect estimated in eq. (18) is fundamentally distinct in nature in comparison withpreviously discussed effect [31–34]. The novel effect which is the subject of this work is entirely rooted to the AQNmodel when the nuggets get disintegrated when enter the solar atmosphere. The corresponding annihilation eventsproduce the low velocities axions with v ≤ v trapped (cid:12) . These axions which behave as DM particles surrounding the Sunhave no relation to the effect discussed in [31–34].Now we want to estimate the number density n axions (cid:12) of these axions assuming, as before, that the majority of theaxions are localized within 2 solar radius R (cid:12) . n axions (cid:12) = N axions (cid:12) π (2 R (cid:12) ) (cid:39) . · (cid:18) ξ − (cid:19) · (cid:18) − eV m a (cid:19) . (19)Can these axions be observed? These axions cannot decay as the axion life time τ ( a → γ ) is very long. However,these axions can be converted to photons in the background of external magnetic field. The corresponding probabilityof this conversion is determined by the formula [35, 36]: P a → γ = (cid:88) q = q ± (cid:18) g aγ B q (cid:19) sin (cid:18) qL (cid:19) , q ± = ± ω − (cid:112) ω − m a (20) To demonstrate the insensitivity to parameter δ , we note that ξ shows very moderate changes between (0 . − × − when δ variesbetween 0 and 1. For the “unphysical value” δ = 8 the parameter ξ (cid:39) . × − , see Appendix A. where L is a typical distance where the magnetic field B is present. For non-relativistic axions one can approximate q ± (cid:39) ± ω . Furthermore, for our present analysis we assume that typical B ∼
300 G in the solar atmosphere, while L is very large such that sin (cid:16) qL (cid:17) can be approximated as . Therefore, probability of the conversion can beapproximated as follows P a → γ (cid:39) (cid:18) g aγ B m a (cid:19) where g aγ m a (cid:39) α π ( m π f π ) · (cid:18) EN −
23 4 + z z (cid:19) z √ z , (21)where z = m u /m d (cid:39) .
56 and parameter
E/N = 0 for KSVZ model, and
E/N = 8 / E/N = 0 to arrive to the following estimate P a → γ (cid:39) (cid:18) g aγ B m a (cid:19) ∼ − (cid:18) B G (cid:19) . (22)The number of the produced photons (as a result of the conversion from the axions) per unit volume with the frequency ω = m a is estimated as follows dN ( a → γ ) dV (cid:39) n axions (cid:12) · P a → γ (cid:39) − (cid:18) B G (cid:19) (cid:18) ξ − (cid:19) (cid:18) − eV m a (cid:19) (23)where n axions (cid:12) is estimated in (19). These converted photons obviously can leave the system. The total number ofphotons leaving the system through area ∼ π (2 R (cid:12) ) per unit time is given by d Φ( a → γ ) dt = dN ( a → γ ) dV (cid:2) π (2 R (cid:12) ) (cid:3) c (cid:39) (cid:18) B G (cid:19) (cid:18) ξ − (cid:19) (cid:18) − eV m a (cid:19) . (24)These photons are very monochromatic with ω = m a with accuracy of order 10 − . Potentially, it gives us some chanceto observe them on Earth. The corresponding count of photons dF ( a → γ ) arriving from the Sun with monochromaticfrequency ω = m a (due to the axion-photon conversion) is estimated as dF ( a → γ ) dA · dt ∼ d Φ( a → γ ) /dt πD (cid:12) ∼ − (cid:18) B G (cid:19) (cid:18) ξ − (cid:19) (cid:18) − eV m a (cid:19) · s . (25)This count, of course, is extremely low. However, these estimates were based on rate (22) corresponding a → γ conversion in vacuum. As it is known since [35] the rate could be drastically enhanced if the system is placed in amedia with non-vanishing plasma frequency ω p exactly matching the axion mass, i.e. ω p = m a , which represents thetopic for the next subsection. B. Solar Corona background. Resonance conversion in solar plasma
We start with numerical estimation for the plasma frequency ω p in the solar corona where the most axions arereleased as a result of the AQN’s annihilation events, ω p ≡ (cid:114) παnm e (cid:39) . · − · (cid:16) n cm − (cid:17) eV . (26)The numerical similarity between ω p and the expected value for the axion mass m a from allowed window m a ∈ (10 − − − ) eV represents the basic motivation for analysis in this subsection. In other words, our goal here isto study possible observational consequences of the resonance case when the condition ω p = m a could occur in thecorona, which is explicit manifestation of the so-called “level-crossing effect” as formulated in ref.[35].If the condition ω p = m a is fulfilled the corresponding resonance a → γ conversion in media is determined byformula [35]: P a → γ = sin (∆ M L ) , ∆ M = B M sin θ, M ≡ g − aγ , cos θ ≡ ˆ (cid:126) B · ˆ (cid:126)k (27) A rough estimate in the following subsection suggests that Lm a ∼ if the resonance condition is satisfied, see Eq. (33). For moregeneral non-resonant case a typical length scale is even larger within the classical axion windows, 10 − eV < m a < − eV becausethere is no requirement for the variation of the plasma frequency on scale L to be small. where we adopted the notations for ∆ M from [35] and expressed the fundamental PQ mass scale M from [35] in termsof the original definition for g aγ . Of course we do not expect that this condition can be exactly satisfied in realityin nature. Furthermore, the oscillation length l deg = π/ ∆ M is very long, much longer than the size of the systemsuch that P a → γ never becomes of order one effect. However, our goal here is different, and we present formula (27)exclusively for illustrative purposes to illuminate the role of the distance scale where the conversion occurs. With thispurpose we expand the resonance expression (27) assuming that ∆ M L (cid:28) P a → γ (cid:39) (∆ M L ) (cid:39) (cid:18) g aγ B m a (cid:19) · (cid:18) m a L (cid:19) , (28)where we consider special case θ = π/ ∼ m − a where this conversiontakes place. Indeed, the first brackets in (28) identically coincides with formula (22) describing the conversion innon-resonance case. Precisely this first term describes a huge suppression factor.For our present studies it is important to emphasize that the same formula (28) also explicitly shows that thissuppressed conversion (22) can be greatly enhanced with the second factor ∼ ( m a L ) if one can increase the coherencelength L by maintaining ω p = m a . If the coherence can be maintained on much larger scale than m − a such that( m a L ) (cid:29) P a → γ will be strongly enhanced in comparison with (22) by this large factor( m a L ) (cid:29) l deg = π/ ∆ M . However, some enhancement in comparison with (22)still can be achieved.The same conclusion also follows from the following expression which was derived using the perturbation theory bytreating the inhomogeneities of the magnetic field and plasma density as small perturbations [35] P a → γ = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:90) L dz ∆ M ( z ) · exp (cid:18) i ∆ a z − i (cid:90) z dz (cid:48) ∆ || ( z (cid:48) ) (cid:19)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , ∆ a = − m a ω , ∆ || = − ω p ω , (29)where we neglected ∆ vac || ∼ B which is numerically very small for relatively weak typical solar magnetic field.The conversion rate given by eq. (29) generates the enhancement proportional to the large length L if the phasesmaintain the coherence and the cancellations between different phases do not occur due to the fast fluctuations. Therequirement that the coherence is maintained up to the scale L is determined by the following condition (cid:32) ∆ a L − (cid:90) L dz (cid:48) ∆ || ( z (cid:48) ) (cid:33) (cid:46) π. (30)If this condition is fulfilled then the conversion rate given by eq. (29) reduces to our previous expression (28) withenhancement factor ∼ L , i.e. P a → γ ∼ (∆ M L ) ∼ L . (31)This supports our previous conclusion that the enhancement factor ( m a L ) is a result of constructive interference.The corresponding length scale L is determined by condition (30).Now we want to address the following question: What is the typical length scale L where the condition (30) can besatisfied in solar atmosphere? We limit our analysis with the trapped axions which have non-relativistic velocities with ω (cid:39) m a as discussed in previous section III A. These axions are distributed in the entire solar atmosphere. Therefore,there is always an extended region in corona or chromosphere where the electron density n is such that the plasmafrequency (26) equals the axion mass, i.e. m a = ω p . The only question remains to be answered is what are the typicallength scales where the average value (cid:104) n (cid:105) for the electron density varies .To estimate the corresponding scale L we notice that a typical variation of the density (and the plasma frequency) n by factor ∼
10 occurs when the altitude changes by ∼ km. Assuming a linear extrapolation (excluding very fastchanges in the transition region) one should expect that the variation of the density δn/n ∼ l ∼ km. This estimate implies that the relative variation (mismatch) of the density on the coherence scale Local fluctuations of the density always occur as a result of different types of waves, including the sound waves, in plasma. However, itis expected that these oscillations do not change the integral entering (30). In other words, we are interested in steady and sustainedvariation of the average density (cid:104) n (cid:105) with latitude and altitude, rather than numerous conventional fluctuations which always occur inhot plasma but do not modify the average magnitude of the integral (30). L of the axion/photon oscillation must not exceed λ/L to be consistent with (30). In other words, the scale L wherecoherence (30) can be maintained must satisfy the following condition L ∼ l λL ⇒ L ∼ (cid:112) λl ∼ . · (cid:115) − eV m a cm . (32)Precisely at this coherence scale L the mismatch in plasma frequency is sufficiently small as ( δn/n ) L ∼ λ/L . At thesame time ( δn/n ) l ∼ ( λ/L ) · ( l /L ) ∼ ∼ l where coherence, of course,cannot be maintained. One should emphasize that this very rough estimate assumes a linear extrapolation. Thisassumption may or may not be justified in reality in solar atmosphere. One should emphasize here that any variationof the magnetic field entering (29) do not modify our estimate for the coherence length (32). This is because theestimate (32) is sensitive to the phase variation (rather than to the amplitude changes ∼ ∆ M ( z )) determined by asteady and systematic variation of the plasma frequency ω p ∼ √ n in corona.If one literally accepts the estimate (32) the corresponding enhancement factor can be approximated as follows,( Lm a ) ∼ · (cid:16) m a − eV (cid:17) (cid:29) . (33)Our previous (non-resonance) estimate (25) should be multiplied by the enhancement factor (33) for case if theresonance conditions can be satisfied and the linear extrapolation is justified. The rate of conversion with this factorbecomes dF ( a → γ ) dA · dt (cid:12)(cid:12)(cid:12) resonance ≈ − (cid:18) B G (cid:19) (cid:18) ξ − (cid:19) · s . (34)The corresponding energy flux can be estimated as m a dF ( a → γ ) dA · dt (cid:12)(cid:12)(cid:12) resonance ≈ − (cid:18) B G (cid:19) (cid:18) ξ − (cid:19) (cid:16) m a − eV (cid:17) W cm , (35)where we expressed the intensity in conventional ( W/ cm ) units using the relations 1 W = 10 (erg / s) and 1 eV =1 . · − erg.While the count (34), (35) is still very low, some hope is that this is an unique monochromatic line. Furthermore,the intensity of this line must be correlated with the EUV emission from corona. In addition, during the flaresthe magnetic field B might be very large in the solar system which provides some enhancement factor and possiblecorrelations with the flares. Finally, this monochromatic line can be, in principle, discriminated from the backgroundnoise as it should appear only along the line-of-sight in the direction of the Sun.It is very instructive to compare the intensity (35) with corresponding conventional energy flux from the Sun in thefrequency band ω ≈ m a . To proceed with the estimates we recall that the total solar intensity (integrated over allfrequency bands) at Earth surface is about 0 . · W cm , which of course many orders of magnitude higher than the rate(35). However, we should compare (35) not with the total intensity from the Sun measured on Earth, but rather withthe solar energy flux from the low energy frequency band with ω ≤ m a . The corresponding estimates can be easilyperformed as the corresponding spectral properties are determined by the Reyleigh- Jeans formula for low energy tailof the black-body (BB) radiation, i.e. dE ω = Tπ ω dω, E BBtot = π T , E ω ≤ m a E BBtot ∼ π (cid:16) m a T (cid:17) . (36)The intensity of the conventional BB solar radiation in the low energy frequency band with ω ≤ m a is estimated asfollows (cid:18) . · W cm (cid:19) · π (cid:16) m a T (cid:17) ∼ . · − (cid:16) m a − eV (cid:17) W cm , (37)which is still much higher than the energy flux due to the axion conversion (35) within conventional axion masswindow. Only for ultra light axions with m a ≤ − eV the background (37) becomes below the signal (35). Thecorresponding case of ultra light axion is not part of this work and shall not be further elaborated.Therefore, one should completely remove the emission from the photosphere with large background (37) for analyzingof the axion conversion (35) from corona. Such removing occurs naturally during the solar eclipses. In practice,astronomers in the past have developed a number of technical tools which allowed to study a weak emission fromcorona by removing a much stronger radiation from photosphere. A high resolution instrument would be very beneficial1to study the monochromatic line (35) with the width ∆ ω , which is determined by the escape velocity v trapped (cid:12) from(13), i.e. ∆ ω ∼ m a v trapped (cid:12) ∼ − m a .To recapitulate: the expected signal from the solar corona is very low as our estimate (35) suggests. It can beonly studied if the background radiation (37) from photosphere can be removed from analysis and the instrumentalresolution is sufficiently high to study the highly monochromatic emission (with the width ∆ ω/ω ∼ − ) from corona.Needless to say that a strong magnetic field in a detector is not required for the observation of these photons on Earthbecause the axion-photon conversion occurs in the solar atmosphere rather than on Earth. In a sense we use entireSun as a one big helioscope where the trapped axions have been accumulated during 4.5 billion years and where theaxions can be converted to photons in the entire solar atmosphere.In next subsection we consider much more optimistic case when the axions are trapped by the Earth. As we shallsee below, the density of such axions could be the same order of magnitude as the galactic axion density, which isthe conventional normalization point for the most presently operational (or under construction, or in stage of design)axion search experiments. C. Axions from the Earth’s Underground
According to ref. [1] the axion flux due to the AQN annihilation events in the very deep underground is given by(7). We integrate this rate over entire surface to arrive dE ⊕ dt (trapped axions) ∼ ξ ⊕ · (cid:18) ∆ BB (cid:19) · πR ⊕ eVs ∼ · (cid:18) ξ ⊕ − (cid:19) · (cid:18) ∆ BB (cid:19) eVs (38)where we used the expression (3) for the rate of the energy transfer to the axions. We also introduced the suppressionfactor ξ ⊕ to account for the small fraction of the trapped axions with v ≤ v trapped . For numerical estimates in formula(14) we use suppression factor ξ ⊕ ∼ ξ (cid:12) · ( v ⊕ /v (cid:12) ) ∼ − computed in previous section and given by (13). This isof course very tiny rate even when ∆ B/B ∼ (cid:39) s.Therefore, the total energy accumulated by the Earth and related to AQN annihilation events radiating the slowvelocity axions can be estimated as follows E ⊕ (trapped axions) (cid:39) · (cid:18) ξ ⊕ − (cid:19) (cid:18) ∆ BB (cid:19) · eVs · s (cid:39) (cid:18) ξ ⊕ − (cid:19) (cid:18) ∆ BB (cid:19) eV . (39)This energy can be expressed in terms of extra Earth’s mass ∆ M ⊕ accumulated by the Earth and represented by thetrapped axions ∆ M ⊕ (trapped axions) (cid:39) . (cid:18) ξ ⊕ − (cid:19) (cid:18) ∆ BB (cid:19) kg , (40)which of course represents a very tiny fraction of the Earth mass M ⊕ (cid:39) . · kg.The energy (39) corresponds to the following total number of the axions accumulated by the Earth during itslife-time: N axions ⊕ ∼ E ⊕ (trapped axions) m a c (cid:39) · (cid:18) ξ ⊕ − (cid:19) (cid:18) ∆ BB (cid:19) · (cid:18) − eV m a (cid:19) . (41)If we assume that the majority of these axions are localized within radius R ⊕ , we arrive to the following estimate forthe average axion energy density inside this volume ρ axions ⊕ ∼ E ⊕ (trapped axions) πR ⊕ ∼ . (cid:18) ξ ⊕ − (cid:19) (cid:18) ∆ BB (cid:19) GeVcm . (42)which is amazingly close to the average dark matter density today ρ DM (cid:39) . GeVcm . The eq. (42) should be viewedas the order of magnitude estimate at the very best. The main uncertainty here is that the trapped axions are notdistributed uniformly, as assumed in (42). Instead, they are obviously distributed in a highly nontrivial way determinedby the position of the nugget when emission occurs (in deep underground) and the direction of the velocity at themoment of emission. Though the estimate (42) is rough, it is also very promising as it suggests that the galactic axiondensity and the axion density produced by the AQN mechanism could be the same order of magnitude.2One may wonder if the terrestrial geomagnetic field can be used as the axion converter, and if the resonanceconversion may occur on Earth, similar to our discussions in previous Sect. III B devoted to the resonance conversionin solar corona. Unfortunately, for the conventional axion mass window m a ∈ (10 − − − ) eV the resonanceconditions cannot be satisfied. Indeed, while the Earth’s ionosphere is highly ionized, the corresponding electrondensity is very low: n ∼ cm − for the so-called F-layer which extends from the altitude 150 km for few hundredkilometers. The corresponding plasma frequency in ionosphere ω p ≡ (cid:114) παnm e (cid:39) . · − · (cid:16) n cm − (cid:17) eV (43)is well below the typical axion mass. This estimate suggests that the resonance case cannot be realized for theconventional axion mass window. Therefore, one should use non-resonance formula for conversion: m a dF ⊕ ( a → γ ) dA · dt (cid:39) ρ axions ⊕ P a → γ c (cid:39) − (cid:18) ξ ⊕ − (cid:19) (cid:18) ∆ BB (cid:19) (cid:18) B . G (cid:19) Wcm . (44)This estimate indicates that the corresponding rate is too low to be observed if one uses the terrestrial geomagneticfield as the axion converter. The crucial suppression factors here are the small terrestrial geomagnetic magnetic field B ∼ . ξ ⊕ ∼ − .Similar arguments also apply to other planets such as Jupiter, which has no resonant enhancement nor sufficientlystrong magnetic field comparable to the solar sunspots. Therefore, while the AQNs obviously get annihilated in theJupiter’s underground producing additional internal heat, the corresponding axion emission would be even smallerthan from the Sun (35) due to a number of additional suppression factors such as smaller mass (and therefore, smallerimpact parameter leading to a smaller AQN flux hitting Jupiter), larger distance from Earth, smaller escape velocity(leading to a smaller parameter ξ ) in comparison with the solar ξ (cid:12) , etc.Therefore we return to our main and most promising estimate (42) which indicates that the density of the boundaxions could be the same order of magnitude as the galactic DM axions, and therefore the conventional instrumentsoriginally designed for the galactic axion searches can be also used to study the trapped axions. The distinct featureof the AQN trapped axions is very large wave length λ a = ( m a v a ) − as the typical trapped velocity v a ≤ v trapped ⊕ is much smaller than a typical galactic DM velocity ∼ − c according to (13). This unique feature of the trappedaxions might be the “smoking gun” leading to their discovery. IV. CONCLUSION AND FUTURE DIRECTIONS
This work represents a natural generalization of the previous studies [1] to properly account for the production ofthe low energy axions when the AQNs get annihilated in the Sun or Earth and emit axions with v ≤ v trapped (cid:12) in thesolar corona or v ≤ v trapped ⊕ in the deep Earth’s underground. This portion of the non-relativistic axions is extremelytiny as shown on Fig. 1b. However, the effect is drastically enhanced as argued in Section III due to accumulation ofthese axions during entire life of the solar system, i.e. for ∼ ρ axions ⊕ . We think this estimate has a hugediscovery potential because ρ axions ⊕ is relatively large and comparable with the average galactic dark matter densitytoday ρ DM (cid:39) . GeVcm . What is more important is that the spectral features of the trapped axions are very distinctfrom conventional galactic axions because the typical trapped velocity v a ≤ v trapped ⊕ is much smaller than a typicalgalactic DM velocity ∼ − c according to (13). Therefore, the typical wave length λ a = ( m a v a ) − of these axions ismuch longer in comparison with galactic axions. This unique feature makes the trapped axions are very distinct fromconventional galactic axions. These axions obviously can be easily discriminated from anything else. The discoveryof such axions would be a “smoking gun” for the entire AQN proposal unifying the DM and baryogenesis (separationof charges) problems.This new mechanism of the axion production is entirely based on the unorthodox AQN dark matter model. Whywe think that this new AQN framework (and accompanying the axion emission) should be taken seriously? We referto [1] for overview of this DM model. Nevertheless, we want to make few comments here suggesting that the AQNframework should be indeed taken seriously. We are thankful to anonymous Referee who suggested to produce such estimates. dark ∼ Ω visible between visible and dark matter densities. In context of the present work the most important featureof this model is that it may potentially resolve the old renowned puzzle (since 1939) known in the community underthe name “the Solar Corona Mystery”. In particular, this model, without adjusting any parameters, generates theobserved EUV luminosity (2) as reviewed in [1].Furthermore, the AQN resolution of the solar corona puzzle also resolves another mystery [37] where it was claimedthat a number of highly unusual phenomena observed in solar atmosphere might be related to the gravitational lensingof “invisible” streaming matter towards the Sun which is correlated with positions of the planets. This is really aweird correlation because one should not expect any connections between the flare occurrences, the intensity of theEUV radiation, and the position of the planets. At the same time, such “weird” correlations naturally occur withinAQN framework. This is because the dark matter AQNs, being the “invisible streaming matter” (in terminology ofref. [37]) can play the role of the triggers sparking the large flares [38]. Therefore, the observation of the correlationbetween the EUV intensity, the frequency of the flares and positions of the planets can be considered as an additionalsupporting argument of the dark matter explanation of the observed EUV irradiation (2), because both effects areoriginated from the same dark matter AQNs.Last, but not least. The AQN model offers a very natural resolution of the so-called “Primordial Lithium Puzzle”as recently argued in [39]. This problem has been with us for at least two decades, and conventional astrophysicaland nuclear physics proposals could not resolve this longstanding mystery. In the AQN framework this puzzle isautomatically and naturally resolved without adjusting any parameters as shown in [39]. This resolution representsyet another, though indirect, support for this new AQN framework.All these arguments obviously represent indirect support for the AQN paradigm. The discovery of the trappedaxions with energy density (42) and with drastically distinct spectral features (in comparison with conventionalgalactic axions) would be the direct support for this model as it is hard to imagine any other model which couldproduce the axions with v a ≤ v trapped ⊕ with sharp cutoff in density for v a > v trapped ⊕ . We conclude this work on thisoptimistic note. ACKNOWLEDGEMENTS
One of us (AZ) is thankful to Konstantin Zioutas for enormous number of questions which motivated these studies.We are thankful to Ludo van Waerbeke for discussions on feasibility and perspective to drastically improve the orderof magnitude estimates (18), (42) using numerical Monte Carlo simulations accounting for all possible axion orbits.This research was supported in part by the Natural Sciences and Engineering Research Council of Canada.
Appendix A: Technical details. Axion emission from the domain wall. 3D case
In this Appendix we want to study the spectral properties of the axion’s emission as a result of time-dependentperturbations of the axion domain wall. We want to focus on the axion portion of the axion DW, which also includesother fields such as π, η (cid:48) , see [40]. It also contains a phase describing the baryon charge distribution on the surfaceof the nugget as discussed in [25]. Exact features of the profile functions for all these fields are not important for ourpurposes. Therefore, one can simplify our computations by considering the following effective Lagrangian with twodegenerate vacuum states . S [ φ ] = (cid:90) d x (cid:34)
12 ( ∂ µ φ ) − g (cid:18) φ − π f a (cid:19) (cid:35) , (A1)where g = √ π m a f a , and we set the effective axion angle as φ/f a ≡ θ + arg det M + π/ π/ ∂ ∂r φ ( r ) + 2 r ∂∂r φ ( r ) = g φ ( r ) · (cid:20) φ ( r ) − π f a (cid:21) , φ ( R ) = 0 (A2) In our previous studies [25–27] we always discussed the so-called N = 1 domain walls. It implies that the vacuum is unique and theDW solution interpolates between one and the same physical vacuum. This interpolation always occurs as a result of variation of theaxion field together with another fields, such as π or η (cid:48) as discussed in [40]. These additional fields do not generate much changes inthe domain wall tension, nor they affect our analysis of the axion production, which is the subject of the present work. Therefore, weignore these fields to simplify notations and qualitative analysis in this work. R defines the boundary which separates two distinct physical vacua and it coincides with the radius of theAQN in equilibrium. While the exact domain wall solution to Eq. (A2) is hard to solve, the approximate solutiongives φ w,R ( r ) (cid:39) π f a · R eff r tanh (cid:20) m a ( r − R ) (cid:21) , R (cid:46) r ≤ R trans π f a · tanh (cid:20) m a ( r − R + δ R ) (cid:21) , r > R trans (A3)where R eff and δ R are functions of a tunable parameter R trans δ R (cid:39) R ( R trans − R ) ≡ m a δ, (A4a) R eff (cid:39) R trans tanh[ R m a R ]tanh[ m a R trans ] (A4b)within domain R < R trans (cid:46) √ m a R · m − a . One can explicitly check this approximate solution (A3) is continuousand first order differentiable. Also, it is precisely the exact solution in the near-field limit r ∼ R and the far-fieldlimit r (cid:29) m − a . Hence, the only unknown part the solution is the “transition” regime between these two limits, wherewe introduce a tunable parameter R trans to account for this type of error source. We will keep this parameter in thefollowing calculations, so it serves as a probe to test whether the final result is sensitive to our crude approach in thetransition regime. As we will see, the final result is not sensitive to the tuning of R trans .Lastly, instead of using R trans directly, it is more convenient to define a simple parameter δ ≡ m a δ R which roughlyvaries from 0 to 1. As we will see, δ is the only parameter entering the final result.We are now ready to compute the excitations χ ( t, z ) in the time dependent background. These excitations will beeventually identified as the axions emitted by the axions DW. To achieve this task we expand φ ( t, z ) = φ w ( z − R ) + χ ( t, z ), which gives S [ φ ] = S [ φ w ] + (cid:90) dt (cid:90) d x (cid:20)
12 ˙ χ − χL χ (cid:21) + O ( χ ) . (A5)where L is a linear differential operator of the second order, L χ = − r ∂ ( rχ ) ∂r − r (cid:20) θ ∂∂θ (cid:18) sin θ ∂χ∂θ (cid:19) + 1sin θ ∂ χ∂φ (cid:21) + (cid:2) g φ χ + g ( φ − v ) χ (cid:3)(cid:12)(cid:12)(cid:12)(cid:12) φ = φ w,R = − r ∂ ( rχ ) ∂r − r (cid:20) θ ∂∂θ (cid:18) sin θ ∂χ∂θ (cid:19) + 1sin θ ∂ χ∂φ (cid:21) + 12 m a v (3 φ w,R − v ) χ. (A6)The corresponding equation of motion is therefore ∂ ∂t χ = − L χ. (A7)To look for the initial conditions, we now want to describe the emission of axions in one cycle of oscillation. Asmentioned in Sec. II B, annihilation of baryon charge results in oscillations of domain wall. Assuming the oscillationis approximately adiabatic, it is sufficient to only analyze the first half of an oscillation – say, the “contraction period”–where the domain wall shrinks from R to a slightly smaller size R − ∆ R . We assumed the rest half of the cycle,the “expansion period”, is just the time-reversed and produces an equivalent contribution. We may write down suchinitial conditions as φ (0 , r ) = φ w,R ( r ) (A8a) φ ( 12 t osc , r ) = φ w,R − ∆ R ( r ) + (excitations) (A8b) Note that interaction between axion and other fields such as π and η (cid:48) becomes strong within r (cid:46) R , see [40]. Hence, we should set acutoff range at r (cid:46) R where Eq. (A2) is no longer valid. t osc denotes the period of one full oscillation. The excitation modes in condition (A8b) is unknown and dependson the conversion rate from excitation modes to freely propagating axions. In terms of χ , the initial conditions (A8)imply χ (0 , r ) = 0 (A9a) χ (0 , r ) = η ( θ, ϕ ) ∂ R [ φ w,R ( r )]∆ R + O (∆ R ) (A9b)where we introduce a free parameter η ( θ, ϕ ) which may be interpreted as the “amplitude of efficiency” of the conversionrate from excitations to free axions, so η must vary between 0 to 1. However, η here may be also interpreted as acorrection term like δ in the approximate solution (A3) within the transition regime R (cid:28) r (cid:46) m − a , so η can begreater than 1 in general. Nonetheless, we will expect η ∼ η can be expanded by partial waves η ( θ, ϕ ) = ∞ (cid:88) l =0 l (cid:88) m = − l η lm Y lm ( θ, ϕ ) ,η lm = (cid:90) π dϕ (cid:90) π dθ sin θ Y ∗ lm ( θ, ϕ ) η ( θ, ϕ ) . (A10)If we assume a good spherical symmetry preserves during the most period of the annihilation process of AQN, then η will be the dominant contribution and η be the next order correction.To solve for the excitation mode, it is convenient to write χ in terms of some normalized basis. The expansion forfree wave is conventionally χ ( t, r, θ, ϕ ) = ∞ (cid:88) l =0 l (cid:88) m = − l (cid:90) d p a plm ( t ) χ plm ( r, θ, ϕ ) , χ plm ( r, θ, ϕ ) = 1 √ π E a j l ( pr ) Y lm ( θ, ϕ ) (A11)where j l ( x ) is the spherical Bessel function, and we have implicitly used two orthogonalities (cid:90) ∞ dr r j l ( pr ) j l ( qr ) = π p δ ( p − q ) , (A12a) (cid:90) π dϕ (cid:90) π dθ sin θ Y ∗ lm ( θ, ϕ ) Y l (cid:48) m (cid:48) ( θ, ϕ ) = δ ll (cid:48) δ mm (cid:48) . (A12b)Note that L is diagonal in basis of χ plm (cid:90) d xχ ∗ ql (cid:48) m (cid:48) ( r, θ, ϕ ) L χ plm ( r, θ, ϕ ) = 18 πE a δ ( p − q ) + m a π E a K ( l ) pq (cid:90) dr r j l ( qr ) j l ( pr )= δ ( p − q )8 πE a p ( p + K ( l ) p,q m a ) , (A13)where K ( l ) p,q is a coefficient defined as K ( l ) p,q ≡ lim L →∞ (cid:82) L dr r j l ( pr ) j l ( qr ) (cid:34) (cid:18) v φ w,R ( r ) (cid:19) − (cid:35)(cid:82) L dr r j l ( pr ) j l ( qr ) (A14)for simplicity of calculation. In B we can show K ( l ) p,q δ ( p − q ) = δ ( p − q ). Then Eq. (A7) is simplified to d dt a plm ( t ) = − E a ( p ) a plm ( t ) , E a ( p ) ≡ (cid:112) p + m a , (A15)which clearly has solution a plm ( t ) = b plm sin E a t, (A16)6following the initial condition (A9a), where b plm is an time-independent coefficient to be determined. To find b plm ,we should impose the second initial condition (A9b) which implies b plm = π f a m a ∆ R η lm sin( E a t osc ) (cid:114) E a π (cid:40)(cid:90) R trans dr r · R eff r sech (cid:104) m a r − R ) (cid:105) j l ( pr ) − (cid:90) R trans dr r · sech (cid:104) m a r − R + δ R ) (cid:105) j l ( pr )+ (cid:90) ∞ dr r · sech (cid:104) m a r − R + δ R ) (cid:105) j l ( pr ) (cid:27) . (A17)Note that only the last term in the curly bracket is dominant because R trans (cid:28) m − a largely suppresses the first twoterms. Thus, we conclude b plm (cid:39) π f a m a ∆ R η lm sin( E a t osc ) (cid:114) E a π (cid:26)(cid:90) ∞ dr r · sech (cid:104) m a r + δ R ) (cid:105) j l ( pr ) + O ( R l +2trans ) (cid:27) , (A18)where we have also drop R in the hyperbolic secant function because it is of order R trans . This integral can beevaluated precisely if we expand the hyperbolic secant assech ( 12 x ) = e − x ∞ (cid:88) n =0 ( − n ( n + 1) 12 n ( e − x − n = ∞ (cid:88) n =0 n (cid:88) k =0 n ( n + 1)! k !( n − k )! ( − k e − ( k +1) x (A19)and use the fact (cid:90) ∞ dρ ρ · e − ( k +1) ρ j l ( pρ ) = √ π l +1 p l ( k + 1) l +3 Γ( l + 3) · f (cid:18)
12 ( l + 3) ,
12 ( l + 4) , l + 32 ; − p ( k + 1) (cid:19) f ( a, b, c ; z ) ≡ c ) F ( a, b, c, z ) , (A20)where F ( a, b, c ; z ) is the Gauss hypergeometric function, and f ( a, b, c ; z ) is defined to be the regularized version of F ( a, b, c, z ) in a conventional way, see Refs. [41, 42] and recent article [43]. As discussed in Sec. III, we are especiallyinterested in the non-relativistic domain, in this limit we have f (cid:18)
12 ( l + 3) ,
12 ( l + 4) , l + 32 ; − p ( k + 1) (cid:19) (cid:39) l + ) (cid:20) − ( l + 3)( l + 4)4( k + 1) Γ( l + )Γ( l + ) p + O ( p ) (cid:21) (A21)Combing Eqs. (A16), (A18), (A20), and (A21), we conclude a plm ( t ) = η lm f a ∆ R e − δ l +3 m a (cid:112) πE a sin( E a t )sin( E a t osc ) Γ( l + 3)Γ( l + ) (cid:18) pm a (cid:19) l H l ( p, δ ) (cid:39) η lm f a ∆ R e − δ l +3 m a (cid:112) πE a sin( E a t )sin( E a t osc ) Γ( l + 3)Γ( l + ) (cid:18) pm a (cid:19) l H l (0 , δ ) (cid:2) O ( p/m a ) (cid:3) (A22)where we define H l ( p, δ ) to be the summation series H l ( p, δ ) ≡ ∞ (cid:88) n =0 n (cid:88) k =0 e − kδ n ( n + 1)! k !( n − k )! ( − k ( k + 1) l +3 Γ( l + 32 ) f (cid:18)
12 ( l + 3) ,
12 ( l + 4) , l + 32 ; − ( p/m a ) ( k + 1) (cid:19) . (A23) More specifically, due to the fact R (cid:28) m − a , we have the hierarchy R < R trans (cid:46) √ m a R · m − a (cid:28) m − a . E rad of the domain wall is obviously E rad = (cid:90) d x χ (cid:20) − ∂ ∂t + L (cid:21) χ = (cid:88) lm (cid:90) d p E a | a plm | = (cid:88) lm (cid:90) ∞ m a dE a · πp E a | a plm | . (A24)More generally, assuming now the flux is produced within a “cavity of radiation” V rad , the density of radiation energy(per unit volume) is therefore E rad /V rad . Then the net flux Φ rad going through the boundary of the cavity is clearly1 S rad ddE a Φ rad = pE a ddE a (cid:18) E rad V rad (cid:19) = (cid:88) lm πp V rad | a plm | . (A25)Let R rad ≡ V rad S rad defines the effective size of cavity of radiation, we obtain ddE a Φ rad = (cid:88) lm η lm R rad f a ∆ R m a π e − δ l +5 (cid:20) Γ( l + 3)Γ( l + ) (cid:21) (cid:20) sin( E a t )sin( E a t osc ) (cid:21) E a (cid:18) pm a (cid:19) l +2 | H l ( p, δ ) | (cid:39) (cid:88) lm η lm R rad f a ∆ R m a π e − δ l +5 (cid:20) Γ( l + 3)Γ( l + ) (cid:21) (cid:20) sin( E a t )sin( E a t osc ) (cid:21) E a (cid:18) pm a (cid:19) l +2 | H l (0 , δ ) | + O ( p/m a ) l +4 . (A26)A few comments should be made regarding to the magnitude of R rad . First, V rad is defined as the cavity whereradiation happens, so R rad (cid:39) ∆ R in 1D case where thin-wall approximation is assumed. However, in 3D R rad mayextend to order of R or even m − a . More generally, it is reasonable to conjecture R rad can depend on the angularmomentum l . It is clear that to compute or even estimate the order of R rad is very difficult. Thus, we should notbother the details of R rad , but rather treat it as a tunable normalization parameter and maybe absorb it into η lm ifapplicable.We also express the spectra as a function of flux velocity ddv a Φ rad = (cid:88) lm η lm R rad f a ∆ R m a π e − δ l +5 (cid:20) Γ( l + 3)Γ( l + ) (cid:21) (cid:20) sin( E a t )sin( E a t osc ) (cid:21) E a (cid:18) pm a (cid:19) l +3 | H l ( p, δ ) | (cid:39) (cid:88) lm η lm R rad f a ∆ R m a π e − δ l +5 (cid:20) Γ( l + 3)Γ( l + ) (cid:21) (cid:20) sin( E a t )sin( E a t osc ) (cid:21) E a (cid:18) pm a (cid:19) l +3 | H l (0 , δ ) | + O ( p/m a ) l +5 . (A27)If the spherical symmetry is well preserved during most period of the annihilation of AQNs, then η is the dominantterm and Eq. A27 can be simplified considerably. We plot this (normalized) result in Figs. 1. One can see the finalresult is not very sensitive to the parameter δ . Such spectrum indicates an average energy (cid:104) E a (cid:105) (cid:39) . m a and anaverage velocity (cid:104) v a (cid:105) (cid:39) . c . Comparing to the 1D case where (cid:104) E a (cid:105) (cid:39) . m a and (cid:104) v a (cid:105) (cid:39) . c [1], we conclude thatthe general features of the spectrum in the relativistic regime for v a ≥ . c is qualitatively consistent between 1D and3D cases as anticipated in the original work [1]. In particular, the difference between these two cases is about 20%for average velocity (cid:104) v a (cid:105) , and about 14% for average energy (cid:104) E a (cid:105) . However, the spectra in the non-relativistic regime v a (cid:28) c are dramatically different, see Fig. 2.Lastly, it is instructive to compare our approximate analytical solution (A3) to the exact numerical solution. Todo this, we plot the corresponding solutions in Fig. 3 with R chosen to be 0.01 m − a for demonstration purpose. Wefind the exact numerical solution in general has a steeper growth, and approaches the outer vacuum expectation valuefaster comparing to approximate solutions with 0 ≤ δ (cid:46)
1. Beyond δ (cid:46) δ = 8,but we should also note such solution badly violates the continuity and first-order differentiability at the transitionzone r ∼ R trans . We also plot δ = 8 in the flux spectrum Figs. 1, which again gives qualitatively consistent answer.As a final remark, one should not consider the spectrum with δ = 8 any better than other values of δ for tworeasons. First, solution with δ = 8 is not physical for its bad violation of continuity and smoothness. Second, evenan exact numerical solution gives no better quantitative prediction because the effective Lagrangian (A1) is only aphenomenological model for qualitative analysis, see footnote 6.8 / c ρ ( v ) FIG. 2: Axion flux spectrum: 1D versus 3D case. Here 1D case (gray dotted) computed in [1] is compared with the3D case (blue solid, δ = 0). m a r0.00.20.40.60.81.0 ( f a ) ( r ) = 0= 1= 8Exact FIG. 3: Approximate and exact numerical solutions to Eq. (A2). Here we choose m a R = 0 . Appendix B: About K ( l ) p,q This appendix is devoted to prove K ( l ) p,q δ ( p − q ) = δ ( p − q ). Before we proceed to the proof, it is convenient to definetwo operators d l = ddr + l + 1 r , d † l = − ddr + l + 1 r . (B1)9And we will need some useful identities: d l d † l [ rj l ( r )] = rj l ( r ) . (B2a) d † l [ rj l ( r )] = rj l +1 ( r ) (B2b) d † l d l = d l +1 d † l +1 (B2c) (cid:90) ∞ [ d l A ( r )] · B ( r ) = (cid:90) ∞ A ( r ) · [ d † l B ( r )] (if A, B = 0 at r = 0 , ∞ ) (B2d)Now we are able to prove by induction. First, the proof in case of l = 0 is quite trivial: K (0) p,q = lim L →∞ (cid:82) L dr sin( pr ) sin( qr ) (cid:104) (cid:0) v φ w,R ( r ) (cid:1) − (cid:105)(cid:82) L dr sin( pr ) sin( qr )= lim L →∞ (cid:82) L dr [cos( p − q ) r − cos( p + q ) r ] (cid:104) (cid:0) v φ w,R ( r ) (cid:1) − (cid:105)(cid:82) L dr [cos( p − q ) r − cos( p + q ) r ] . (B3)Up to this point, we note that K (0) p,q is always finite for any positive p, q >
0. Now, if we multiply both sides by δ ( p − q ) K (0) p,q δ ( p − q ) = δ ( p − q ) lim L →∞ (cid:82) L dr [1 − cos( p + q ) r ] (cid:104) (cid:0) v φ w,R ( r ) (cid:1) − (cid:105)(cid:82) L dr [1 − cos( p + q ) r ]= δ ( p − q ) , (B4)where we know the integral is quickly dominant by the the term (cid:82) L dr ·
1, so that the fraction in the limit L → ∞ gives trivial result.Now, we want to show that K ( l ) p,q δ ( p − q ) = K ( l +1) p,q δ ( p − q ) for all l = 0 , , ... . First, let us see that K ( l ) p,q = lim L →∞ (cid:82) L dr · d l d † l [ prj l ( pr )] · d l d † l [ qrj l ( qr )] · (cid:104) (cid:0) v φ w,R ( r ) (cid:1) − (cid:105)(cid:82) L dr · d l d † l [ prj l ( pr )] · d l d † l [ qrj l ( qr )]= lim L →∞ (cid:82) L dr · prj l +1 ( pr ) · d † l d l [ qrj l +1 ( qr )] · (cid:104) (cid:0) v φ w,R ( r ) (cid:1) − (cid:105)(cid:82) L dr · prj l +1 ( pr ) · d † l d l [ qrj l +1 ( qr ) −− lim L →∞ (cid:82) L dr · prj l +1 ( pr ) · qrj l +1 ( qr ) · ddr (cid:104) (cid:0) v φ w,R ( r ) (cid:1) − (cid:105)(cid:82) L dr · prj l +1 ( pr ) · d † l d l [ qrj l +1 ( qr )]= K ( l +1) p,q − lim L →∞ (cid:82) L dr · r j l +1 ( pr ) · j l +1 ( qr ) · ddr (cid:104) (cid:0) v φ w,R ( r ) (cid:1) − (cid:105)(cid:82) L dr · r j l +1 ( pr ) · j l +1 ( qr ) , (B5)where we have applied Eqs. (B2a) to (B2d) in the intermediate steps. Again, the integral is finite for any p, q > p = q , the second term must vanish. It is because the numerator is obviously finite,while the denominator tends to infinity in the large L limit as indicated by Eq. A12a. Thus, we conclude K ( l ) p,q δ ( p − q ) = δ ( p − q ) (B6)for all l = 0 , , , ... as expected. [1] H. Fischer, X. Liang, Y. Semertzidis, A. Zhitnitsky and K. Zioutas, Phys. Rev. D , no. 4, 043013 (2018) [arXiv:1805.05184[hep-ph]]. [2] R. D. Peccei and H. R. Quinn, Phys. Rev. D , 1791 (1977);S. Weinberg, Phys. Rev. Lett. , 223 (1978);F. Wilczek, Phys. Rev. Lett. , 279 (1978).[3] J.E. Kim, Phys. Rev. Lett. (1979) 103;M.A. Shifman, A.I. Vainshtein, and V.I. Zakharov, Nucl. Phys. B166 (1980) 493(KSVZ-axion).[4] M. Dine, W. Fischler, and M. Srednicki, Phys. Lett.
B104 (1981) 199;A.R. Zhitnitsky, Yad.Fiz. (1980) 497; Sov. J. Nucl. Phys. (1980) 260 (DFSZ-axion).[5] K. van Bibber and L. J. Rosenberg, Phys. Today , 30 (2006);[6] S. J. Asztalos, L. J. Rosenberg, K. van Bibber, P. Sikivie, K. Zioutas, Ann. Rev. Nucl. Part. Sci. , 293-326 (2006).[7] Pierre Sikivie, Lect. Notes Phys. , 19 (2008) arXiv:0610440v2 [astro-ph].[8] G. G. Raffelt, Lect. Notes Phys. , 51 (2008) [hep-ph/0611350].[9] P. Sikivie, Int. J. Mod. Phys. A , 554 (2010) [arXiv:0909.0949 [hep-ph]].[10] L. J. Rosenberg, Proc. Nat. Acad. Sci. (2015),[11] P. W. Graham, I. G. Irastorza, S. K. Lamoreaux, A. Lindner and K. A. van Bibber, Ann. Rev. Nucl. Part. Sci. , 485(2015) [arXiv:1602.00039 [hep-ex]].[12] D. J. E. Marsh, Phys. Rept. , 1 (2016) [arXiv:1510.07633 [astro-ph.CO]].[13] A. Ringwald, PoS NOW , 081 (2016) [arXiv:1612.08933 [hep-ph]].[14] J. Preskill, M. B. Wise, and F. Wilczek, Phys.Lett. B120 , 127 (1983);L. Abbott and P. Sikivie, Phys.Lett. B , 133 (1983);M. Dine and W. Fischler, Phys.Lett. B , 137 (1983).[15] S. Chang, C. Hagmann and P. Sikivie, Phys. Rev. D , 023505 (1999) [hep-ph/9807374].[16] T. Hiramatsu, M. Kawasaki, K. Saikawa and T. Sekiguchi, Phys. Rev. D , 105020 (2012) Erratum: [Phys. Rev. D ,089902 (2012)] [arXiv:1202.5851 [hep-ph]].[17] M. Kawasaki, K. Saikawa and T. Sekiguchi, Phys. Rev. D , no. 6, 065014 (2015) [arXiv:1412.0789 [hep-ph]].[18] L. Fleury and G. D. Moore, JCAP , 004 (2016) [arXiv:1509.00026 [hep-ph]].[19] M. Gorghetto, E. Hardy and G. Villadoro, JHEP , 151 (2018) [arXiv:1806.04677 [hep-ph]].[20] V. B. Klaer and G. D. Moore, JCAP , no. 11, 049 (2017) [arXiv:1708.07521 [hep-ph]].[21] P. Sikivie, Phys. Rev. Lett. , 1415 (1983) Erratum: [Phys. Rev. Lett. , 695 (1984)].[22] S. Andriamonje et al. [CAST Collaboration], JCAP , 010 (2007) [hep-ex/0702006].[23] L. DiLella, K. Zioutas, Astroparticle Physics, , 145 (2003). [arXiv:astro-ph/0207073v1 ].[24] A. R. Zhitnitsky, JCAP , 010 (2003) [hep-ph/0202161].[25] X. Liang and A. Zhitnitsky, Phys. Rev. D , 083502 (2016) [arXiv:1606.00435 [hep-ph]].[26] S. Ge, X. Liang and A. Zhitnitsky, Phys. Rev. D , no. 6, 063514 (2017) [arXiv:1702.04354 [hep-ph]].[27] S. Ge, X. Liang and A. Zhitnitsky, Phys. Rev. D , no. 4, 043008 (2018) [arXiv:1711.06271 [hep-ph]].[28] A. Zhitnitsky, JCAP , no. 10, 050 (2017) [arXiv:1707.03400 [astro-ph.SR]].[29] N. Raza, L. van Waerbeke and A. Zhitnitsky, Phys. Rev. D , no. 10, 103527 (2018) [arXiv:1805.01897 [astro-ph.SR]].[30] S. L. Adler, Phys. Rev. D , 023505 (2009) doi:10.1103/PhysRevD.79.023505 [arXiv:0805.2895 [astro-ph]].[31] X. Xu and E. R. Siegel, arXiv:0806.3767 [astro-ph].[32] I. B. Khriplovich and D. L. Shepelyansky, Int. J. Mod. Phys. D , 1237 (1988).[36] A. N. Ioannisian, N. Kazarian, A. J. Millar and G. G. Raffelt, JCAP , no. 09, 005 (2017) [arXiv:1707.00701 [hep-ph]].[37] S. Bertolucci, K.Zioutas, S. Hofmann, M. Maroudas, Phys. Dark Univ. , 13 (2017) [arXiv:1602.03666 [astro-ph.CO][38] A. Zhitnitsky, Phys. Dark Univ. , 1 (2018) [arXiv:1801.01509 [astro-ph.SR]].[39] V. V. Flambaum and A. R. Zhitnitsky, Phys. Rev. D , xxx (2019). arXiv:1811.01965 [hep-ph].[40] M. M. Forbes and A. R. Zhitnitsky, JHEP , 013 (2001) [hep-ph/0008315].[41] M. Abramowitz, Chap. 15 Hypergeometric Functions, Handbook of Mathematical Functions, edited by M. Abramowitzand I.A. Stegun, National Bureau of Standards, Applied Mathematics Series - 55 (1972).[42] M. Abramowitz, Chap. 6 Gamma function and related functions, Handbook of Mathematical Functions, edited by M.Abramowitz and I.A. Stegun, National Bureau of Standards, Applied Mathematics Series - 55 (1972).[43] N. Michel, M.V. Stoitsov, “Fast computation of the Gauss hypergeometric function with all its parameters complex withapplication to the PschlTellerGinocchio potential wave functions,” Computer Physics Communications178(7)