Constraints on the coupling with photons of heavy axion-like-particles from Globular Clusters

We update the globular cluster bound on massive ($m_a$ up to a few 100 keV) axion-like particles (ALP) interacting with photons. The production of such particles in the stellar core is dominated by the Primakoff $\gamma + Ze\to Ze +a$ and by the photon coalescence process $\gamma+\gamma\to a$. The latter, which is predominant at high masses, was not included in previous estimations. Consequently, our result strengthens the previous constraints considerably, especially for $m_a \gtrsim 50$ keV. The combined constraints from Globular Cluster stars, SN 1987A, and beam-dump experiments leave a small triangularly shaped region open in the parameter space around $m_a \sim 0.5-1\,$ MeV and $g_{a\gamma} \sim 10^{-5}$ GeV$^{-1}$. This is informally known as the ALP"cosmological triangle"since it can be excluded only using standard cosmological arguments. As we shall discuss, however, there are viable cosmological models that are compatible with axion-like particles with parameters in such region. We also discuss possibilities to explore the cosmological triangle experimentally in upcoming accelerator experiments.


I. INTRODUCTION
Axion-like-particles (ALPs) with masses m a in the keV-MeV range emerge in different extension of the Standard Model, as Pseudo-Goldstone bosons of some broken global symmetry. The theoretical speculation about superheavy axion models began long ago (see Sec. 6.7 of Ref. [1] for a recent review), in an attempt to get rid of the strong astrophysical bounds on the axion coupling, which made it effectively invisible. In this context, superheavy means heavier than about 100 keV, so that the axion production in most stars (supernovae and neutron stars being an exception) is Boltzmann suppressed and the majority of the stellar axion bounds are relaxed. Nowadays, several mechanisms exist to increase the axion mass independently from its couplings, without spoiling the solution of the strong CP problem (a list of references can be found in [1]).
Besides QCD axions, heavy ALPs emerge in compactification scenarios of string theory [2][3][4], or in the context of "relaxion" models [5]. Heavy ALPs have also recently received considerable attention in the context of Dark Matter model-building. Indeed, they may act as mediators for the interactions between the Dark Sector and Standard Model (SM) allowing to reproduce the correct Dark Matter relic abundance via thermal freezeout [6]. ALPs with masses below the MeV scale can have a wide range of implications for cosmology and astrophysics (see [7] for a review), affecting for example Big Bang Nucleosynthesis (BBN), the Cosmic Microwave Background (CMB) [8][9][10] and the evolution of stars. Colliders and beam-dump experiments are also capable to explore this mass range, indeed reaching the m a ∼ O(GeV) frontier, which is not covered by any as-trophysical or cosmological considerations [7,11,12].
In this work we are interested in ALPs interacting exclusively with photons. Additional couplings with SM fields, particularly with electrons, may spoil some of our conclusions. For such ALPs, the collection of all the astrophysical and experimental constraints leaves a triangular area in the parameter space, for masses m a ∼ 0.5 − 1 MeV and couplings g aγ ∼ 10 −5 GeV −1 , open. Although the existence of ALPs with such parameters is in tension with standard cosmological arguments [8,9], the region of such masses and couplings passes the current experimental tests and all the known astrophysical arguments, and is also permitted in viable non-standard cosmological scenarios [10]. Because of that, this parameter area is sometimes dubbed as the ALP "cosmological triangle". As we shall discuss in Sec. V, this region is now the target of several direct investigations, as more and more experiments are reaching the sensitivity to probe those masss and couplings, and there is a chance that such area might be covered in the next decade or so. Redefining the boundaries of the cosmological triangle is, therefore, particularly timely and relevant to guide the experimental investigations.
In this work we revisit the globular cluster bound on heavy ALPs, which defines the low-mass boundary of the cosmological triangle. Globular Clusters (GC) are gravitationally bound systems of stars, typically harboring a few millions stars. Being among the oldest objects in the Milky Way, their population is made of lowmass stars (M < 1M ). Most of these stars belong to the so-called Main Sequence, which corresponds to the H burning evolutionary phase. However, there are two other well defined evolutionary phases, i.e., the Red Giant Branch (RGB) and the Horizontal Branch (HB). The first is made by cool giant stars, burning H in a thin shell surrounding a compact He-rich core. During the RGB phase the stellar luminosity increases and the core contracts, until the temperature rises enough to ignite He. Then, stars leave the RGB and enter the HB phase, during which they burn He, in the core, and H, in the shell.
The number of stars found in the different evolutionary phases depends linearly on the time spent by a star in each of them. For this reason, stellar counts provide a powerful tool to investigate the efficiency of the energy sources and sinks in stellar interiors, those that affect the stellar lifetime τ in a given stage of the stellar evolution. In this context, the GC R parameter, defined as the number ratio of horizontal branch to red giants branch stars, i.e.: is a powerful observable often used to investigate stellar physics. In particular, it has been also exploited to constrain the axion-photon coupling g aγ [13][14][15][16], at least for ALPs light enough, m a ∼ < 30 keV, that their production is not Boltzmann suppressed. At such low masses, the most relevant axion production mechanism induced by the photon coupling is the Primakoff process, γ + Ze → Ze + a, i.e. the conversion of a photon into an ALP in the electric field of nuclei and electrons in the stellar plasma (cf. Sec. II). This process is considerably more efficient in HB than in RGB stars, since in the latter case it is suppressed due to the larger screening scale and plasma frequency (see Sec. II). Therefore, the energy-loss caused by the production of ALPs with a sizable g aγ would imply a reduction of the HB lifetime and, in turn, a reduction of the R parameter. As it turns out, R has a substantial dependence, approximately linear, on the helium abundance of the cluster and, if ALPs are also included, a quadratic dependence on the axionphoton coupling. On the other hand, the R parameter is only marginally affected by a variation of the cluster age and metallicity. Thus, once the He abundance is known from direct or indirect measurements, bounds (or hints) on the axion-photon coupling can be obtained from the comparison of the R parameter measured in Globular clusters with the theoretical expectations obtained by varying g aγ [15][16][17][18]. An accurate application of this method, based on photometric data for 39 GCs, was discussed in [15] by some of us, who found an upper bound g aγ < 0.66 × 10 −10 GeV −1 at 95 % confidence level, a value more recently experimentally confirmed by the CAST collaboration [19]. The goal of the present work is to extend the GC bound on g aγ to higher axion masses. Given the typical temperature T ∼ 10 keV in the stellar core of a HB star, one expects the thermal production of particles to be Boltzmann suppressed for m a > ∼ 30 keV, relaxing the bound on g aγ . A quantitative analysis was carried out in [9], where a bound was derived from the requirement that the axion energy emitted per unit time and mass, ε a , averaged over a typical HB core, satisfies the requirement ε a ∼ < 10 erg g −1 s −1 [20]. However, that analysis neglected the contribution of the photon coalescence process γγ → a (cf. Sec. II), to the ALP production in stars. At low masses, this process is subdominant and it is forbidden for m a < 2 ω pl , where ω pl is the plasma frequency at the position where the process takes place. Hence, the inclusion of the photon coalescence does not affect the bound obtained in Ref. [15], which remains valid for light axions (m a ∼ < 10 keV). As we shall show in Sec. II, however, for masses m a > ∼ 50 keV, the coalescence production dominates and becomes several times larger than the Primakoff at masses > ∼ 100 keV. Thus, not surprisingly, the inclusion of this process strengthens the results considerably at masses of 100 keV and above.
The plan of our work is the following. In Sec. II, we revise the axion emissivity via the Primakoff process and the photon coalescence. In Sec. III, we discuss our procedure and present our bound on g aγ for massive ALPs. Then, we show the complementary of our bound with other constraints. In Sec. IV, with that from SN 1987A (in the trapping regime), and in Sec. V with the experimental bounds from beam-dump searches. Finally, in Sec. VI we summarize our results and we conclude.

II. AXION EMISSIVITY
The ALP-two photon vertex is described by the Lagrangian term where g aγ is the ALP-photon coupling constant (which has dimension of an inverse energy), F the electromagnetic field andF its dual. The primary production mechanisms for ALPs interacting with photons in the core of a HB star are: • the Primakoff process γ + Ze → Ze + a, where a thermal photon in the stellar core converts into an axion in the Coulomb fields of nuclei and electrons; • the photon coalescence process γγ → a, where two photons in a medium of sufficiently high density annihilate producing an axion.
As we shall see, the former dominates at low ALP masses (m a ∼ < 50 keV) while at large mass the photon coalescence takes over.
There is a vast literature on the axion Primakoff production rate. The interested reader may consult Ref. [9,13] for a detailed discussion. Here, we provide only a brief review and present some results applicable in the typical plasma conditions relevant for this work. In general, the axion emission rate (energy per mass per time) via the Primakoff process is given by the expression where the factor 2 comes from the photon degrees of freedom, ρ is the local density, f (E) = (e E/T − 1) −1 is the Bose-Einstein distribution, and Γ γ→a is the photon-axion transition rate, In the last expression, E and p = E 2 − m 2 a are, respectively, the ALP energy and momentum. The photon obeys the dispersion relation k = ω 2 − ω 2 pl where k is the photon momentum, ω its energy, and ω 2 pl 4παn e /m e is the plasma frequency (or effective "photon mass"). In a photon-axion transition the energy is conserved because we ignore recoil effects. Therefore, we use ω = E. Finally κ is the Debye-Hückel screening scale with ρ the mass density, m u = 1.66 × 10 −24 g the atomic mass unit, Y e the number of electrons per baryon, and Y j the number per unit mass of the ions with nuclear charge Z j . At low density, the electrons are non-degenerate and the ion correlation can be neglected, so that the Debye-Hückel theory provides a valid description of the plasma screening. This condition is certainly fulfilled in the core of a HB star. Note also, that in the center of a HB star, T ∼ 8.6 keV, ρ ∼ 10 4 g cm −3 , and ω pl ∼ 3 keV. Thus, the plasma frequency is considerably smaller than the thermal energy. Nevertheless, to achieve a higher accuracy our numerical Primakoff emission rate includes also the effects induced by a finite plasma frequency (a detailed description of the adopted emission rate can be found the Appendix in Ref. [22], which we have generalized at finite ALP mass). The axion coalescence process, γγ → a, has a kinematic threshold, vanishing for m a ≤ 2ω pl [23]. As we shall see, above this threshold the production rate is a steep function of the mass and dominates over the Primakoff at m a > ∼ 50 keV. In order to calculate the axion coalescence rate in a thermal medium, it is convenient to approximate the Bose-Einstein photon distribution with a Maxwell-Boltzmann, f (E) → e −E/T , for the photon occupation number [23]. This approximation is justified since we are interested only in axion masses (and thus axion energies) of the order of the temperature or larger (for m a ∼ < T the coalescence process is practically negligible). Therefore, the production rate per unit volume of ALPs of energy between E and E + dE is [23] The temperature and density profiles within the Herich core of a typical HB stellar model are shown in Fig. 1. The model has been evolved starting from the pre-main sequence up to the end of the core He burning phase. For the initial structure (t = 0 model) we have adopted a mass M=0.82 M and, as usual, a homogeneus composition, namely: Y = 0.25 and Z = 0.001. After ∼ 13 Gyr the central He burning begins (zero age HB). At that time the stellar mass is m ∼ 0.72 M , while the mass of the He-rich core is m ∼ 0.5 M . Fig. 1 is a snapshot of the stellar core taken when the central mass fraction of He reduces down to X He ∼ 0.6. The corresponding Primakoff and photon coalescence emission rates are compared in Fig. 2. The quantity reported in the vertical axis is the ratio of the energy-loss rate, in units of erg g −1 s −1 , and the square of the axion-photon coupling, g 10 ≡ g aγ /10 −10 GeV −1 . The Primakoff and photon coalescence emission rates have been computed for two different values of the axion mass, namely: m a = 30 keV and m a = 80 keV. In the case of m a = 30 keV, the Primakoff energy-loss rate (in the center of the star) is a factor of ∼ 3 larger than the photon coalescence rate. Conversely, for m a = 80 keV the photon coalescence dominates and the contribution of the Primakoff is effectively negligible. Fig. 3 shows how the ALP luminosity, depends on the ALP mass m a . The integration is extended from the center (r = 0) to the stellar surface. According to our expectations, the coalescence process is sub-leading for m a ∼ < 50 keV, but it dominates at higher masses. Note that the stellar luminosity, as due to the photon emission, is L ∼ 2 × 10 35 erg s −1 , which is about 2 orders of magnitude larger than the axion luminosity at g 10 ∼ 1.

III. GLOBULAR CLUSTER BOUND
In order to derive a bound on g aγ for massive ALPs, we have computed several evolutionary sequences of stellar models, from the pre-main-sequence to the end of the core He burning. The models have been computed by means of FuNS (Full Network Stellar evolution), an hydrostatic 1D stellar evolution code [22]. In general, the inclusion of the axion energy-loss in stellar model computations leads to a reduction of the R parameter, defined in Eq. (1). On the other hand, the larger the initial He abundance the larger the estimated R. In practice, the upper bound on the axion-photon coupling is obtained when the largest possible value of the He abundance is assumed. Analyzing the He abundance measured in molecular clouds with metallicity in the same range of those of galactic GCs, in Ref. [15] it was estimated a conservative upper limit for the He abundance, specifically Y = 0.26. Adopting this value of Y , it was shown that the R parameter obtained from photometric observations of 39 GCs, R = 1.39 ± 0.03, implies the stringent upper bound g aγ = 0.66 × 10 −10 GeV −1 (95 % C.L.). However, this bound is only valid for light axions.
Since ALPs interacting only with photons are not efficiently produced in the core of RGB stars, and hence affect minimally the RGB lifetime (τ RGB in Eq. (1)), the variation of R due to an axion production is essentially a consequence of the reduction of the HB lifetime (τ HB in Eq. (1)). We have computed τ HB for a GC benchmark. Specifically, we used: age 13 Gyr, metallicity Z = 0.001, and Y = 0.26, corresponding to the conservative upper limit for the GC He abundance reported in [15]. In the standard case, when no exotic energy-loss process is included, we found τ HB = 8.84 × 10 7 yr. The addition of light axions with with g aγ = 0.66 × 10 −10 GeV −1 reduces the HB lifetime down to τ HB = 7.69 × 10 7 yr. Requiring the HB lifetime to be within these values guarantees that the predicted R parameter is consistent, within 2 σ, with the observed one.
The argument was generalized to massive ALPs by searching for the ALP-photon coupling that reduces the HB lifetime down to τ HB = 7.69 × 10 7 yr at each fixed ALP mass. The result is shown in the exclusion plot reported in Fig. 4. The area shaded in light red, delimited by a red line, represents the ALP parameter region excluded at 95 % C.L. according to the criterion discussed above. 1 It is evident how the bound loses its strength for masses above ∼ 30 keV, because of the Boltzmann suppression of the axion emissivity, reaching g aγ ∼ < 10 −5 GeV −1 for m a ∼ 500 keV. Comparing with the gray dashed line, which shows the bound obtained assuming only the Primakoff process, we see that for masses of a few 10 keV, the photon coalescence is no longer negligible. Indeed, we find that including this channel the bound is strengthened by a factor > ∼ 4 at m a > ∼ 100 keV, and by over an order of magnitude for m a > ∼ 200 keV. For reference, in the figure we are also showing, in 1 Notice that our derivation of the bound fails at very large couplings and masses. The reason is that we are assuming that ALPs are free streaming out of the star and do not decay into photons before leaving it. These assumptions are not satisfied for gaγ > ∼ 10 −3 GeV −1 and ma ∼ a few 100 keV [9]. This is outside the region shown in our Fig. 4. However, the reader should be warned that the HB boundary shown in Fig. 5 loses validity at those high couplings (very top of the figure). In any case, this parameter region is already excluded by direct experimental searches (see Sec. V). light green, the region excluded by SN 1987A in the regime of ALPs trapped in the SN core (see Sec. IV), and in blue the parameters excluded by direct searches at beam dump experiments (see Sec. V). Interestingly, the combination of all the astrophysical and experimental bounds leave a small triangular area, roughly at m a ∼ 0.5−1 MeV and g aγ ∼ 10 −5 GeV −1 , unconstrained. This is the ALP cosmological triangle. Standard cosmological arguments, particularly concerning BBN and the allowed effective number of relativistic species, N eff , can be used to exclude this area [9,10]. Nevertheless, in non-standard cosmological scenarios, e.g. in low-reheating models, the cosmological bounds can be relaxed all the way to the GC bound calculated in this work [10]. Thus, the cosmological triangle is still a viable region of the ALP parameter space, open to experimental and phenomenological investigations.

IV. SN 1987A BOUND FROM AXION TRAPPING
For the sake of completeness, in this section we present briefly our derivation of the SN 1987A constraint on heavy ALPs presented in Fig. 4 and 5. A detailed study of this constraint, based on state-of-the-art SN models [24], is currently ongoing and will be the topic of a forthcoming work by some of us [25]. Here, we just present a succinct discussion of the SN argument to constraint the ALP-photon couplings at the bottom edge of the cosmological triangle.
Heavy ALPs can be copiously produced in a supernova (SN) core via Primakoff and coalescence processes. Due to the higher core temperature, T ∼ O(30) MeV, SNe can be used to probe ALP masses considerably larger than those probed by GCs (see, e.g. [7,26,27]). For couplings of interest in this work, g aγ ∼ O(10 −5 ) GeV −1 , ALPs would be trapped in the SN, having a mean-free path smaller than the size of the SN core (R ∼ 10 km) [7,26]. In this case, ALPs may contribute significantly to the energy transport in the star, modifying the SN evolution. Since SN 1987A neutrino data are in a reasonable agreement with core-collapse SN models without the emission of exotic species, one should require that ALPs interact more strongly than the particles which provide the standard mode of energy transfer, i.e. neutrinos.
When ALPs interact strongly enough to be trapped in the SN core, they are emitted from an axion-sphere, a spherical shell whose radius r a is fixed by the optical depth being about unity. More specifically, we calculated r a imposing that the optical depth where κ a is the axion opacity, satisfies the condition τ a (r a ) 2/3. This is analogous to the neutrino last scattering surface, i.e. the "neutrino-sphere", with radius r ν . Trapped ALPs have a black-body emission with a luminosity L a ∝ r 2 a T 4 (r a ). In order to obtain the bound on g aγ one should impose [20,21] L a ∼ < L ν .
We are concerned mostly with a time posterior to 0.5-1 s, where the outer core has settled and the shock has begun to escape. Specifically, in our numerical calculation we refer to the SN model used in [28], for a representative post-bounce time t pb = 1 s.
We calculated the ALP opacity following the prescriptions in [29] (see [27] for an alternative approach). For masses m a ∼ < a few MeV, the dominant contribution to the axion opacity is due to the inverse Primakoff process, a + Ze → γ + Ze, where λ a→γ is the mean free-path, and β E = (1 − m 2 a /E 2 ) 1/2 . The inverse Primakoff rate is Γ a→γ = 2Γ γ→a , with Γ γ→a given in Eq. (4).
From κ a→γ one can calculate the mean ALP Rosseland opacity [20] where is the ALP thermal spectrum.
We derived our bound on axion coupling from the luminosity condition in Eq. (10), taking the axion-sphere radius that satisfies Eq. (9). As shown in Fig. 4, for m a < 10 MeV, the luminosity condition excludes the values of the photon-axion coupling g aγ ∼ < 8 × 10 −6 GeV −1 , in agreement with previous results [7,26].
Note that the SN 1987A bound should not be considered at the same level of confidence as the GCs one, since it is not based on a self-consistent SN simulations. Performing such a simulation, which should include also the trapped ALPs, would be a challenging task (see, e.g., [30] for a recent investigation in the context of dark photons), and demand a separated investigation.

V. DIRECT EXPERIMENTAL TESTS OF THE COSMOLOGICAL TRIANGLE
As discussed above, the ALP region at masses of a few MeV is the target of numerous investigations. In this section, we briefly comment on the existing experimental limits near the cosmological triangle and future prospects to test that region directly in experiments. Fig. 5 shows an overview plot of the status of the search for heavy ALPs with our updated bound as discussed in this work in red, labeled "HB". Colored regions are excluded at 95% C.L. Other limits are compiled from references [7,12] and detailed therein. The experimental limits which are "nose-like-shaped" (E137, CHARM, nuCal, E141) are from beam-dump setups, in which the ALP needs to live long enough to reach the detection volume (boundary at "large" couplings and masses). However, it should not be so long-lived that it can excape from it (boundary at "small" couplings and masses).
The most efficient experiment to "touch" the cosmic triangle was E137, shown as a blue-shaded region in Fig. 4. This bound is based on data published by the experiment E137 [31] and its revisit in [7]: around 2×10 20 electrons were dumped into an aluminum target, potentially yielding to Primakoff-production of ALPs. However, no excess of expected photon signals was observed at a distance of ∼ 200 m, leading to an exclusion limit.
The small-coupling-limit of E137 relevant for us in this context is largely determined by how long-lived ALPs can be while still being detected by the experiment. The limit estimated [7] for this reason seems robust as late ALP decays will suffer little from their non-negligible probability of showering in air. We thus show this limit in Fig. 5. Roughly spoken, a long baseline together with a relatively soft ALP spectrum (compared to proton dumps whose lower limits are at much larger couplings [12]), made E137 an ideal fixed-target in probing the cosmological triangle at the top section of its parameter space.
As for the possibilities to probe the remaining region at m a ∼ 1 MeV, Ref. [7] details on prospects to significantly probe the cosmic triangle at Belle-II at a statistics of 50 ab −1 . Sensitivity is also expected at "active" beam dumps such as LDMX-type set-ups, that can infer the presence of ALPs through a "missing-momentum signature" [32]. The running experiment PADME, at Frascati, does not currently have the potential to reach the cosmological triangle [33,34] but could be potentially sensitive to this area after a luminosity upgrade.
It is worth stressing that far more experimental options to probe this triangle exist if the axion-coupling is not limited strictly to direct photon couplings [27]. However, this possibility is outside the assumptions made in our work.

VI. DISCUSSIONS AND CONCLUSIONS
In this work we have extended the GCs bound on the ALP-photon coupling to masses m a > ∼ 10 keV, in the region of the parameter space where the Boltzmann suppression of the axion emission rate can no longer be neglected. Our analysis improves on the previous work by including the coalescence process, γ + γ → a, which is the dominating axion production mechanism at masses above ∼ 50 keV. The bound is shown in Fig. 4 (red line), where we also compare it to the bound obtained ignoring the coalescence process (dashed gray line). The inclusion of the coalescence reduced the allowed value of the axion photon coupling by a factor of ∼ 4 at masses ∼ 100 keV, and by over an order of magnitude at m a > ∼ 200 keV. Our bound reduces slightly the region at masses ∼ 0.5 to 1 MeV and couplings g aγ ∼ 10 −5 GeV −1 (the cosmological triangle) cornered on the other sides by SN 1987A and direct searches (see Fig. 5). Though excluded by standard cosmological arguments, this is a viable region in non-standard cosmological scenarios, e.g. in low-reheating models, which relax substantially the cosmological bounds [10]. Thus, the cosmological triangle is an area of great experimental interest, as shown in our Fig. 5. Indeed, a number of theoretical models permit ALPs (and even QCD axions) with parameters in this region, as discussed in Sec. I, making this a possible target area for future experimental investigations. Interestingly, a detection of an axion signal in this region would have dramatic cosmological consequences, requiring non-standard cosmological scenarios. This intriguing possibility confirms once more the nice complementarity between astrophysical, cosmological arguments and direct searches in order to corner or luckily discover axionlike-particles.  [25] and the experimental limits, compiled from Refs. [7,12], are also shown. Prospects to experimentally probe the viable region are commented on in the text.