X-ray Searches for Axions from Super Star Clusters

Axions may be produced in abundance inside stellar cores and then convert into observable X-rays in the Galactic magnetic fields. We focus on the Quintuplet and Westerlund 1 super star clusters, which host large numbers of hot, young stars including Wolf-Rayet stars; these stars produce axions efficiently through the axion-photon coupling. We use Galactic magnetic field models to calculate the expected X-ray flux locally from axions emitted from these clusters. We then combine the axion model predictions with archival NuSTAR data from 10 - 80 keV to search for evidence of axions. We find no significant evidence for axions and constrain the axion-photon coupling $g_{a\gamma\gamma} \lesssim 3.6 \times 10^{-12}$ GeV$^{-1}$ for masses $m_a \lesssim 5 \times 10^{-11}$ eV at 95% confidence.

Ultralight axions that couple weakly to ordinary matter are natural extensions to the Standard Model. For example, string compactifications often predict large numbers of such pseudo-scalar particles that interact with the Standard Model predominantly through dimension-five operators [1,2]. If an axion couples to quantum chromodynamics (QCD) then it may also solve the strong CP problem [3-6].
Axions may interact electromagnetically through the operator L = −g aγγ aF µνF µν /4, where a is the axion field, F is the electromagnetic field-strength tensor, with F its Hodge dual, and g aγγ is the dimensionful coupling constant. This operator allows both the production of axions in stellar plasmas through the Primakoff Process [7] and the conversion of axions to photons in the presence of static external magnetic fields. Strong constraints on g aγγ for low-mass axions come from the CAST experiment [8], which searches for axions produced in the Solar plasma that free stream to Earth and then convert to X-rays in the magnetic field of the CAST detector. CAST has excluded axion couplings g aγγ 6.6 × 10 −11 GeV −1 for axion masses m a 0.02 eV at 95% confidence [8]. Primakoff axion production also opens a new pathway by which stars may cool, and strong limits (g aγγ 6.6 × 10 −11 GeV −1 at 95% confidence for m a keV) are derived from observations of the horizontal branch (HB) star lifetime, which would be modified in the presence of axion cooling [9].
In this work, we produce some of the strongest constraints to-date on g aγγ for m a 10 −9 eV through Xray observations with the NuSTAR telescope of super star clusters (SSCs). The SSCs contain large numbers of hot, young, and massive stars, such as Wolf-Rayet (WR) stars. We show that these stars, and as a result the SSCs, are highly efficient at producing axions with energies ∼10-100 keV through the Primakoff process. These axions may then convert into photons in the Galactic magnetic fields, leading to signatures observable with spacebased X-ray telescopes such as NuSTAR. We analyze archival NuSTAR data from the Quintuplet SSC near the Galactic Center (GC) along with the nearby Westerlund 1 (Wd1) cluster and constrain g aγγ 3.6 × 10 −12 The locations of the stars are indicated in black, while the 90% energy containment region for emission associated with the SSC is indicated by the black circle, accounting for the NuSTAR PSF. RA0 and DEC0 denote the locations of the cluster center. We find no evidence for axion-induced emission from this SSC, which would follow the spatial counts template illustrated in the inset panel.
GeV −1 at 95% confidence for m a 5 × 10 −11 eV. In Fig. 1 we show the locations of the stars within the Quintuplet cluster that are considered in this work on top of the background-subtracted NuSTAR counts, from 10 -80 keV, with the point-spread function (PSF) of NuSTAR also indicated. In the Supplementary Material (SM) we show that observations of the Arches SSC yield similar but slightly weaker limits.
Our work builds upon significant previous efforts to use stars as laboratories to search for axions. Some of the strongest constraints on the axion-matter couplings, for example, come from examining how HB, white dwarf arXiv:2008.03305v1 [hep-ph] 7 Aug 2020 (WD), red giant, and neutron star (NS) cooling would be affected by an axion [9][10][11][12][13][14][15][16][17][18]. When the stars have large magnetic fields, as is the case for WDs and NSs, the axions can be converted to X-rays in the stellar magnetospheres [19][20][21][22]. Intriguingly, in [21,22] observations of the Magnificent Seven nearby isolated NSs found evidence for a hard X-ray excess consistent with the expected axion spectrum from nucleon bremsstrahlung. This work extends these efforts by allowing the axions to convert to X-rays not just in the stellar magnetic fields but also in the Galactic magnetic fields [23][24][25]. Axion production in SSCs.-SSCs are efficient sources of axions because they host large numbers of hot and massive stars, and the axion production rates in stars increase with temperature. A massive star spends the majority of its life on the main sequence (MS), where it burns hydrogen. Once core hydrogen is exhausted helium fusion begins with an initial core temperature of T 10 keV, although the temperature increases as fusion continues. At core helium ignition, the hydrogen envelope expands and the star evolves onto a super-or hyper-giant branch. Stars of higher masses burn through their material more quickly than those of low masses and therefore live shorter lives. Further evolution cools the envelope, and the star changes classification from a hot O supergiant to the cool red supergiant (RSG) phase through the blue supergiant (BSG), luminous blue variable (LBV), and yellow hypergiant (YHG) phases, and depending on its initial mass it may go supernova in any of the last three phases [26].
During helium burning, particularly massive stars may undergo considerable mass loss, especially through either rotation or binary interaction, which can begin to peel away the hydrogen envelope, revealing the hot layers underneath and reversing the cooling trend. Stars undergoing this process are known as WR stars. If the star has a small (<40% abundance) remaining hydrogen envelope, it is classified as a WNh star; at <5% hydrogen abundance it is classified as a WN star; otherwise, it is classified as WC or WO, which indicates the presence of >2% carbon, and oxygen, respectively, in the atmosphere.
Axions are produced through the photon coupling g aγγ in the high-mass stars in SSCs through the Primakoff process γ + (e − , Z) → a + (e − , Z). This process converts a stellar photon to an axion in the screened electromagnetic field of the nucleons and electrons. The massive stars are high-temperature and low-density and therefore form nonrelativistic nondegenerate plasmas. The Primakoff emission rate was calculated in [7,27] and is described in detail in the SM. Importantly, the axion luminosity scales with the axion-photon coupling as g 2 aγγ and with the plasma temperature between T 4 and T 7 , depending on the density. Given the temperature and density profiles from the stellar simulations, we may calculate the axion luminosities for each star.
To compute the axion luminosity in a given star, we use the stellar evolution code Modules for Experiments in Stellar Astrophysics (MESA) [28,29] to find, at any particular time in the stellar evolution, radial profiles of temperature, density, and composition. We use MESA to construct a set of models for each spectroscopic classification within each cluster. The census [30] for Quintuplet lists the spectroscopic classifications for each star in our sample. We adopt the spectroscopic classification definitions in [31,32]. We then obtain a set of models that describe the possible states for each spectroscopic classification in the cluster. The states are specified by an initial metallicity Z, an initial stellar mass, and an initial rotation velocity. The initial metallicity is taken to be constant for all stars. In the SM we show that the Quintuplet and Arches clusters, which are both near the GC, are likely to have initial metallicities in the range Z ∈ (0.018, 0.035), consistent with the conclusions of previous works with place the initial metallicities of these clusters near solar (recall that solar metallicity is Z ≈ 0.02) [33,34]. We vary Z over the range quoted above and find that that the most conservative limits are obtained for Z at the high end of this range and thus take that metallicity as our fiducial value. Note that higher metallicities may lead to the stars entering the WR classifications sooner.
Rotation may also cause certain massive stars to be classified as WR stars at younger ages. We model the initial rotation distribution as a Gaussian distribution with mean µ rot and standard deviation σ rot for non-negative rotation speeds [35,36]. Refs. [35,36] found µ rot ≈ 100 km/s and σ rot ≈ 140 km/s, but to assess systematic uncertainties we vary µ rot between 50 and 150 km/s [35].
We draw initial stellar velocities from the velocity distribution described above (from 0 to 500 km/s) and initial stellar masses from the Kroupa initial mass function [37] (from 15 to 200 M ). We use MESA to evolve the stars from pre-MS to near-supernova. At each time step we assign each stellar model a spectroscopic classification using the definitions in [31,32]. We then construct an ensemble of models for each spectroscopic classification by joining together the results of the different simulations that result in the same classification for stellar ages between 3.0 and 3.6 Myr (for Quintuplet) [30], which is the uncertainty range for the age of Quintuplet. Note that each simulation generally provides multiple representative models, taken at different time steps. In total we compute 10 5 models per stellar classification.
Quintuplet hosts 71 stars of masses 50M , with a substantial WR cohort [30]. In particular it has 14 WC + WN stars, and we find that these stars dominate the predicted axion flux. For example, at g aγγ = 10 −12 GeV −1 we compute that the total axion luminosity from the SSC (with Z = 0.035 and µ rot = 150 km/s) is 2.1 +0.7 −0.4 × 10 35 erg/s, with WC + WN stars contributing ∼70% of that flux. Note that the uncertainties arise from performing multiple (500) draws of the stars from our ensembles of representative models. In the 10 -80 keV energy range relevant for NuSTAR the total luminosity is 1.7 +0.4 −0.3 ×10 35 erg/s. We take Z = 0.035 and µ rot = 150 km/s because these choices lead to the most conservative limits. For example, taking the metallicity at the lower-end of our range (Z = 0.018) along with µ rot = 100 km/s the predicted 10 -80 keV flux increases by ∼60%. At fixed Z = 0.035 changing µ rot from 150 km/s to 100 km/s increases the total luminosity (over all energies) by ∼10%, though the luminosity in the 10 -80 keV range is virtually unaffected.
The Wd1 computations proceed similarly. Wd1 is measured from parallax to be a distance d ∈ (2.2, 4.8) kpc from the Sun [38], accounting for both statistical and systematic uncertainties [39]. Wd1 is estimated to have an age between 4.5 and 7.1 Myr from isochrone fitting, which we have broadened appropriately from [40] accounting for expanded distance uncertainties. In our fiducial analysis we simulate the stars in Wd1 for initial metallicity Z = 0.035 and µ rot = 150 km/s as this leads to the most conservative flux predictions, even though it is likely that the metallicity is closer to solar for Wd1 [41], in which cases the fluxes are larger by almost a factor of two (see the SM). We model 153 stars in Wd1 [40], but the axion flux is predominantly produced by the 8 WC and 14 WN stars. In total we find that the 10 -80 keV luminosity, for g aγγ = 10 −12 GeV, is 9.02 +1.2 −1.1 × 10 35 erg/s, which is ∼5 times larger than that from Quintuplet. Axion conversion in Galactic fields.-The axions produced within the SSCs may convert to X-rays in the Galactic magnetic fields. The axion Lagrangian term L = g aγγ aE · B, written in terms of electric and magnetic fields E and B, causes an incoming axion state to rotate into a polarized electromagnetic wave in the presence of an external magnetic field (see, e.g., [42]). When the conversion probability p a→γ is sufficiently less than unity, it may be approximated by [20]: where B i , for i = 1, 2, denote the two orthogonal projections of the magnetic field onto axes perpendicular to the direction of propagation. The integrals are over the line of sight, with the source located a distance d from Earth, and r = 0 denoting the location of the source. We have also defined ∆ a ≡ −m 2 a /(2E) and ∆ || (r) ≡ −ω pl (r) 2 /(2E), with E the axion energy and ω pl (r) the location-dependent plasma mass. The plasma mass may be related to the number density of free electrons n e by ω pl ≈ 3.7 × 10 −12 (n e /10 −2 cm −3 ) −1/2 eV.
To compute the energy-dependent conversion probabilities p a→γ for our targets we need to know the magnetic field profiles and electron density distributions along the lines of sight. For our fiducial analysis we use the regular components of the JF12 Galactic magnetic field model [43,44] and the YMW16 electron density model [45] (though in the SM we show that the ne2001 [46] model gives similar results). The JF12 model does not cover the inner kpc of the Galaxy, however, which may be the important region for studies of the Quintuplet cluster. Outside of the inner kpc the conversion probability for Quintuplet is dominated by the out-of-plane (X-field) component in the JF12 model. We conservatively assume that the magnitude of the vertical magnetic field within the inner kpc is the same as the value at 1 kpc (|B z | ≈ 3 µG). In our fiducial magnetic field model the conversion probability is p a→γ ≈ 2.4×10 −4 (7×10 −5 ) for g aγγ = 10 −12 GeV −1 for axions produced in the Quintuplet SSC with m a 10 −11 eV and E = 80 keV (E = 10 keV). Completely masking the inner kpc reduces these conversion probabilities to p a→γ ≈ 1.0 × 10 −4 (p a→γ ≈ 3.2 × 10 −5 ), for E = 80 keV (E = 10 keV). On the other hand, changing global magnetic field model to that presented in [47] (PTKN11), which has a larger in-plane component than the JF12 model but no out-of-plane component, leads to conversion probabilities at E = 80 and 10 keV of p a→γ ≈ 4.9 × 10 −4 and p a→γ ≈ 4.2 × 10 −5 , respectively, with the inner kpc masked. We note that since the X-ray flux scales with the axion coupling as g 4 aγγ , the differences between the magnetic field model variations discussed above correspond to ∼20% variations on the inferred values of g aγγ .
In reality the magnetic field is almost certainly larger than the assumed 3 µG within the inner kpc. Note that the local interstellar magnetic field, as measured directly by the Voyager missions [48], indirectly by the Interstellar Boundary Explorer [49], inferred from polarization measurements of nearby stars [50], and inferred from pulsar dispersion measure and the rotation measure data [51], has magnitude B ∼ 3 µG, and all evidence points to the field rising significantly in the inner kpc [52]. For example, Ref. [53] bounded the magnetic field within the inner 400 pc to be at least 50 µG, and more likely 100 µG (but less than ∼400 µG [54]), by studying non-thermal radio emission in the inner Galaxy. Localized features in the magnetic field in the inner kpc may also further enhance the conversion probability beyond what is accounted for here. For example, the line-of-sight to the Quintuplet cluster overlap with the GC radio arc non-thermal filament, which has a ∼3 mG vertical field over a narrow filament of cross-section ∼(10 pc) 2 (see, e.g., [55]). Accounting for the magnetic fields structures described above in the inner few hundred pc may enhance the conversion probabilities by over an order of magnitude relative to our fiducial scenario (see the SM).
When computing the conversion probabilities for Wd1 we need to account for the uncertain distance d to the SSC (with currently-allowable range given above). In the JF12 model we find the minimum p a→γ /d 2 (for m a 10 −11 eV) is obtained for d ≈ 2.6 kpc, which is thus the value we take for our fiducial distance in order to be conservative. At this distance the conversion probability is p a→γ ≈ 2.4×10 −6 (p a→γ ≈ 1.5×10 −6 ) for E = 10 keV (E = 80 keV), assuming g aγγ = 10 −12 GeV −1 and m a 10 −11 eV. We note that the conversion probabilities are over 10 times larger in the PTKN11 model (see the SM), since there is destructive interference (for d ≈ 2.6 kpc) in the JF12 model towards Wd1. We do not account for turbulent fields in this analysis; inclusion of these fields may further increase the conversion probabilities for Wd1, although we leave this modeling for future work. Data analysis.-We reduce and analyze 39 ks of archival NuSTAR data from Quintuplet with observation ID 40010005001. This observation was performed as part of the NuSTAR Hard X-ray Survey of the GC Region [56,57]. The NuSTAR data reduction was performed with the HEASoft software version 6.24 [58]. This process leads to a set of counts, exposure, and background maps for every energy bin and for each exposure (we use data from both Focal Plane Modules A and B). The astrometry of each exposure is calibrated independently using the precise location of the source 1E 1743.1-2843 [59], which is within the field of view. The background maps account for the cosmic Xray background, reflected solar X-rays, and instrumental backgrounds such as Compton-scattered gamma rays and detector and fluorescence emission lines [60]. We then stack and re-bin the data sets to construct pixelated images in each of the energy bins. We use 14 5-keV-wide energy bins between 10 and 80 keV. We label those images d i = {c p i }, where c p i stands for the observed counts in energy bin i and pixel p. The pixelation used in our analysis is illustrated in Fig. 1.
For the Wd1 analysis we reduced 138 ks of Focal Plane Module A and B data from observation IDs 80201050008, 80201050006, and 80201050002. This set of observations was performed to observe outburst activity of the Wd1 magnetar CXOU J164710.2-45521 [61], which we mask at 0.5 in our analysis. (The magnetar is around 1.5' away from the cluster center.) Note that in [61] hard X-ray emission was only detected with the NuSTAR data from 3 -8 keV from CXOU J164710.2-45521 -consistent with this, removing the magnetar mask does not affect our extracted spectrum for the SSC above 15 keV. We use the magnetar in order to perform astrometric calibration of each exposure independently. The Wd1 exposures suffer from ghost-ray contamination [62] from a nearby point source that is outside of the NuSTAR field of view at low energies (below ∼15 keV) [61]. The ghost-ray contamination affects our ability to model the background below 15 keV and so we remove the 10 -15 keV energy bin from our analysis.
In each energy bin we perform a Poissonian template fit over the pixelated data to constrain the number of counts that may arise from the template associated with axion emission from the SSC. To construct the signal template we use a spherically-symmetric approximation to the NuSTAR PSF [63] and we account for each of the stars in the SSC individually in terms of spatial location and expected flux, which generates a non-spherical and extended template. We label the set of signal templates by S p i . We search for emission associated with the signal templates by profiling over background emission. We use the set of background templates described above and constructed when reducing the data, which we label B p i . Given the set of signal and background templates we construct a Poissonian likelihood in each energy bin: We then construct the profile likelihood p i (d i |{S i }) by maximizing the log likelihood at each fixed S i over the nuisance parameter A B . Note that when constructing the profile likelihood we use the region of interest (ROI) where we mask pixels further than 2.0' from the SSC center. The 90% containment radius of NuSTAR is ∼1.74', independent of energy, as indicated in Fig. 1. We use a localized region around our source to minimize possible systematic biases from background mismodeling. However, as we show in the SM our final results are not strongly dependent on the choice of ROI. We also show in the SM that if we inject a synthetic axion signal into the real data and analyze the hybrid data, we correctly recover the simulated axion parameters.
The best-fit flux values and 1σ uncertainties extracted from the profile likelihood procedure are illustrated in Fig. 2 for the Quintuplet and Wd1 data sets. We  Figure 2. The spectra associated with the axion-induced templates from the Quintuplet and Wd1 SSCs constructed from the NuSTAR data analyses, with best-fit points and 1σ uncertainties indicated. In red we show the predicted spectra from an axion with ma 10 −11 eV and indicated gaγγ. Note that for Wd1 we do not analyze the 10 -15 keV energy bin because of ghost-ray contamination.
compare the spectral points to the axion model prediction to constrain the axion model. More precisely, we combine the profile likelihoods together from the individual energy bins to construct a joint likelihood that may be used to search for the presence of an axion signal: p(d|{m a , g aγγ }) where R i (m a , g aγγ ) denotes the predicted number of counts in the i th energy bin given an axion-induced X-ray spectrum with axion model parameters {m a , g aγγ }. The values R i (m a , g aγγ ) are computed using the forwardmodeling matrices constructed during the data reduction process.
Joint 95% limit Figure 3. The 95% upper limits (black) on gaγγ as a function of the axion mass from the Quintuplet and Wd1 data analyses. We compare the limits to the 1σ (green band) and 2σ (yellow band) expectations under the null hypothesis, along with the median expectations (dotted). The joint 95% upper limit, combining Quintuplet and Wd1, is also indicated (expected joint limit not shown). At low masses our limits may be surpassed by those from searches for X-ray spectral modulations from NGC 1275 [64], though we caution that those limits have been called into question recently [65].
In Fig. 3 we illustrate the 95% power-constrained [66] upper limits on g aγγ as a function of the axion mass m a found from our analyses. The joint limit (red in Fig. 3), combining the Quintuplet and Wd1 profile likelihoods, becomes g aγγ 3.6 × 10 −12 GeV −1 at low axion masses. At fixed m a the upper limits are constructed by analyzing the test statistic q(g aγγ |m a ) ≡ 2 ln p(d|{m a , g aγγ }) − 2 ln p(d|{m a ,ḡ aγγ }), whereḡ aγγ is the signal strength that maximizes the likelihood, allowing for the possibility of negative signal strengths as well. The 95% upper limit is given by the value g aγγ >ḡ aγγ such that q(g aγγ |m a ) ≈ 2.71 (see, e.g., [67]). The 1σ and 2σ expectations for the 95% upper limits under the null hypothesis, constructed from the Asimov procedure [67], are also shown in Fig. 3. The evidence in favor of the axion model is ∼0.3σ (0σ) local significance at low masses for Quintuplet (Wd1).
We compare our upper limits with those found from the CAST experiment [8], the non-observation of gammarays from SN1987a [68] (see also [69][70][71] along with [72], who recently questioned the validity of these limits), and the NGC 1275 X-ray spectral modulation search [64]. It was recently pointed out, however, that the limits in [64] are highly dependent on the intracluster magnetic field models and could be orders of magnitude weaker, when accounting for both regular and turbulent fields [65]. Discussion.-We present limits on the axion-photon coupling g aγγ from a search with NuSTAR hard X-ray data for axions emitted from the hot, young stars within SSCs and converting to X-rays in the Galactic magnetic fields. We find the strongest limits from analyses of data towards the Quintuplet and Wd1 clusters. Our limits represent some of the strongest and most robust limits todate on g aγγ for low-mass axions. We find no evidence for axions. Promising targets for future analyses could be nearby supergiant stars, such as Betelgeuse [23,73], or young NSs such as Cas A. There is likely a significant amount of archival X-ray data already that can be used to further search for the existence of axions, and future targeted observations may also prove fruitful. Acknowledgments. - This Supplementary Material contains additional results and explanations of our methods that clarify and support the results presented in the main Letter. First, we present additional details regarding the data analyses, simulations, and calculations performed in this work. We then show additional results beyond those presented in the main Letter. In the last section we provide results of an auxiliary analysis used to derive the metallicity range considered in this work. In this section we first provide additional details needed to reproduce our NuSTAR data reduction, before giving extended discussions of our MESA simulations, axion luminosity calculations, and conversion probability calculations.

A. Data Reduction and analysis
To perform the NuSTAR data reduction, we use the NuSTARDAS software included with HEASoft 6.24 [58]. We first reprocess the data with the NuSTARDAS task nupipeline, which outputs calibrated and screened events files. We use the strict filtering for the South Atlantic Anomaly. We then create counts maps for both focal plane modules (FPMs) of the full NuSTAR FOV with nuproducts in energy bins of width 5 keV from 5 − 80 keV. We additionally generate the ancillary response files (ARFs) and the redistribution matrix files (RMFs) for each FPM. We generate the corresponding exposure maps with nuexpomap, which produces exposure maps with units [s]. To obtain maps in exposure units [cm 2 s keV] that we can use to convert from counts to flux, we multiply in the mean effective areas in each bin with no PSF or vignetting correction.
Once the data is reduced, we apply the analysis procedure described in the main text to measure the spectrum associated with the signal template in each energy bin. However, to compare the signal-template spectrum to the axion model prediction, we need to know how to forward-model the predicted axion-induced flux, which is described in more detail later in the SM, through the instrument response. In particular, we pass the signal flux prediction through the detector response to obtain the expected signal counts that we can compare to the data: Here, t e is the exposure time corresponding to the exposure e in [s], while the signal is the expected intensity spectrum in [erg/cm 2 /s/keV]. We have now obtained the expected signal counts µ e S,i (θ S ) that may be integrated into the likelihood given in (2).

B. MESA Simulations
MESA is a one-dimensional stellar evolution code which solves the equations of stellar structure to simulate the stellar interior at any point in the evolution. In our fiducial analysis, we construct models at a metallicity Z = 0.035, initial stellar masses from 15 to 200 M , and initial surface rotations from 0 km/s to 500 km/s as indicated in the main text. We use the default inlist for high-mass stars provided with MESA. This inlist sets a number of parameters required for high-mass evolution, namely the use of Type 2 opacities. We additionally use the Dutch wind scheme [75] as in the high rotation module.
On this grid, we simulate each star from the pre-MS phase until the onset of neon burning around 1.2 × 10 9 K. At that point, the star only has a few years before undergoing supernova. Given that no supernova has been observed in the SSCs since the observations in 2012-2015, this end-point represents the most evolved possible state of stars in the SSCs at time of observation. The output is a set of radial profiles at many time steps along the stellar evolution. The profiles describe, for example, the temperature, density, and composition of the star. These profiles allow us to compute the axion spectrum at each time step by integrating the axion volume emissivity over the interior.
Here we show detailed results for a representative star of mass 85 M with initial surface rotation of 300 km/s. This star is a template star for the WC phase (and other WR phases) in the Quintuplet Cluster, which dominates the Quintuplet axion spectrum in the energies of interest. In the left panel of Fig. S1 Years before End of Run Figure S1. (Left) The HR diagram for the Quintuplet template star of mass 85 M and initial surface rotation of 300 km/s. The coloring indicates the year before the run was stopped, approximately a few years from supernova. We mark with black squares, in order of occurrence, when the star enters the WNh phase, when it is 3 Myr old, when its core undergoes helium ignition, when it enters the WN, WC, and WO phases, and finally when the run ends at 3.85 Myr. (Right) A logT-log ρ diagram for the template star with the same points of interest marked. We also show the relevant degeneracy zones, showing that the star is entirely in the nonrelativistic nondegenerate regime.
Eventually, the core runs out of hydrogen fuel and is forced to ignite helium to prevent core collapse (see Fig. S2 left). Because helium burns at higher temperatures, the star contracts the core to obtain the thermal energy required to ignite helium (see Fig. S3). At the same time, the radiation pressure in stellar winds cause heavy mass loss in the outer layers, which peels off the hydrogen envelope (see Fig. S4). When the surface is 40% hydrogen, the star enters the WNh phase; when it is 5% hydrogen, the star enters the WN phase. Further mass loss begins to peel off even the helium layers, and the star enters the WC and WO phases when its surface is 2% carbon and oxygen by abundance [32], respectively (see Fig. S2 right).  Fig. S1. With dashed-black vertical lines, we mark several points of interest: "WNh" indicates the time the star enters the WNh phase, "He ignition" when its core undergoes helium ignition, and "WN","WC", and "WO" indicate the beginning of the WN, WC, and WO phases, respectively. (Right) The same as in the left panel, but for surface abundances.

C. Axion Production in SSCs
In this section we overview how we use the output of the MESA simulations to compute axion luminosities and spectra.   Figure S4. The stellar mass (black) and radius (red) as a function of time from the simulation described in Fig. S1. The dashed-black vertical lines retain their meanings from Fig. S2.

The Axion Energy Spectrum
Here we focus on the calculation of the axion energy spectrum [erg/cm 2 /s/keV]. The axion production rate is [76] Γ p (E) = g 2 aγγ T κ 2 32π 1 + κ 2 4E 2 ln 1 + where κ 2 = 4πα T i Z 2 i n i gives the Debye screening scale, which is the finite reach of the Coulomb field in a plasma and cuts off the amplitude. To obtain the axion energy spectrum, this is to be convolved with the photon density, such that where we have defined the dimensionless parameter ξ 2 = κ 2 4T 2 . To obtain the axion emissivity for a whole star, we integrate over the profiles produced with MESA, and we show results for this calculation in Sec. I C 2. Finally, the axion-induced photon spectrum at Earth is given by with the conversion probability P a→γ computed in Sec. I D.

Results for Template Star
In this section, we show our expectation for the axion luminosity from the template star we detailed in Sec. I B, which is a 85 M star used in modelling the WN and WC stars in Quintuplet.
In the left panel of Fig. S5, we show the axion emissivity from the radial slices of the MESA profile, using the model at the start of the WC evolutionary stage. As expected, the stellar core is by far the most emissive due to its high temperature and density. We also show the temperature profile in the star. Note that the axion volume emissivity does not have the same profile shape as the temperature because the emissivity also depends on the density and composition which are highly nonuniform over the interior. In the right panel of Fig. S5, we show how the axion luminosity changes over the stellar lifetime. We see that before helium ignition, the axion luminosity is rather low, and the axion spectrum reaches its maximum around 10 keV, owing to the low core temperature-the star is still hydrogen burning at core temperatures well below 10 keV. During helium ignition, the luminosity increases quickly due to the sudden increase in temperature. During helium burning, the core temperature continues to increase; for this reason, more evolved stars will be more luminous in axions.

D. Magnetic field model and conversion probability
The conversion probabilities in this work are calculated using (1), which is valid under the assumptions outlined in the main Letter. To perform the integral we need to know (i) the free electron density along the line of sight to the target, and (ii) the orthogonal projections of the magnetic field along the line-of-sight. In this section we give further details behind the electron-density and magnetic-field profiles used in this Letter.
The Quintuplet and Arches SSCs are both ∼30 pc away from the GC and thus are expected to have approximately the same conversion probabilities for conversion on the ambient Galactic magnetic fields. It is possible, however, that local field configurations near the GC could enhanced the conversion probabilities for one or both of these sources. For example, the axions are expected to travel through or close to the GC radio arc, which has a strong magnetic field ∼mG over a cross-section ∼(10 pc) 2 [55]. Magnetic fields within the clusters themselves may also be important.
Our fiducial magnetic field model for Quintuplet and Arches is illustrated in the left panel of Fig. S6. In the right panel we show the magnetic field profiles relevant for the Wd1 observations. The components of the B-field along the two transverse directions are denoted by B 1 and B 2 . For the Quintuplet and Arches analyses, the propagation direction is very nearly aligned with −x (in Galactic coordinates), so we may take B 1 to point in theẑ direction, towards the north Galactic pole, and B 2 to point in the directionŷ (the approximate direction of the local rotation). Note that the targets are slightly offset from the origin of the Galactic coordinate system, so the actual basis vectors have small components in the other directions. As Wd1 is essentially within the plane of the disk, one of the transverse components points approximately in theẑ direction (B 1 ).
The dominant magnetic field towards the GC within our fiducial B-field model is the vertical direction (B 1 ), which is due to the out-of-plane X-shaped halo component in the JF12 model [43,44]. However, in the JF12 model that component is cut off within 1 kpc of the GC, due to the fact that in becomes difficult to model the B-field near the GC. The B-field is expected to continue rising near the GC -for example, in [53] it was claimed that the B-field  Figure S6. We denote the projections of the Galactic magnetic field onto the plane normal to the propagation direction by B1, B2. (Left) The transverse magnetic field components in our fiducial model (the JF12 model, black) and alternate model (PTKN11, orange) towards the Quintuplet and Arches clusters. Note that in our fiducial B-field model we extend the JF12 model to distances less than 1 kpc from the GC using the field values at 1 kpc. The true magnetic field values in the inner kpc almost certainly surpass those from this conservative model (see text for details). (Right) The two field components towards the Wd1 cluster, which is taken to be at a distance of 2.6 kpc from the Sun. The conversion probabilities towards Wd1 are much larger in the alternate model (PTKN11) than in our fiducial model (JF12), though we stress that random fields are not included and could play an important role in the conversion probabilities towards Wd1.
should be at least 50 µG (and likely 100 µG) within the inner 400 pc. However, to be conservative in our fiducial B-field model we simply extend the B-field to the GC by assuming it takes the value at 1 kpc (about 3 µG) at all distances less than 1 kpc from the GC. We stress that this field value is likely orders of magnitude less than the actual field strength, but this assumption serves to make our results more robust. The extended field model is illustrated in Fig. S6.
To understand the level of systematic uncertainty arising from the B-field models we also show in Fig. S6 the magnetic field profiles for the alternative ordered B-field model PTKN11 [47]. This model has no out-of-plane component, but the regular B-field within the disk is stronger than in the JF12 model. In the case of Quintuplet and Arches we find, as discussed below, that the PTKN11 model leads to similar but slightly enhanced conversion probabilities relative to the JF12 model. On the other hand, the conversion probabilities in the PTKN11 model towards Wd1 are significantly larger than in the JF12 model.
There is a clear discrepancy in Fig. S6 between the magnetic field values observed at the solar location, in both the JF12 model and the PTKN11 model, and the local magnetic field strength, which is ∼3 µG [49]. The reason is that the magnetic field profiles shown in Fig. S6 are only the regular components; additional random field components are expected. For example, in the JF12 model the average root-mean-square random field value at the solar location is ∼6.6 µG [43,44]. The random field components could play an important role in the axion-to-photon conversion probabilities, especially for the nearby source Wd1, but to accurately account for the random field components one needs to know the domains over which the random fields are coherent. It is expected that these domains are ∼100 pc [44], in which case the random fields may dominate the conversion probabilities, but since the result depends sensitively on the domain sizes, which remain uncertain, we conservatively neglect the random-field components from the analyses in this work (though this would be an interesting subject for future work).
To compute the conversion probabilities we also need the free-electron densities. We use the YMW16 model [45] as our fiducial model, but we also compare our results to those obtained with the older ne2001 model [46] to assess the possible effects of mismodeling the free-electron density. In the left panel of Fig. S7 we compare the free electron densities between the two models as a function of distance away from the Sun towards the GC, while in the right panel we show the free electron densities towards Wd1. The differences between these models result in modest differences between the computed conversion probabilities, as discussed below.
Combining the magnetic field models in Fig. S6 and the free-electron models in Fig. S7 we may compute the axionphoton conversion probabilities, for a given axion energy E. These conversion probabilities are presented in the left panels of Fig free-electron models and various magnetic field configurations.
In the top left panel our fiducial conversion probability model is shown in solid black. Changing to the ne2001 model would in fact slightly enhance the conversion probabilities at most energies, as shown in the dotted black, though the change is modest. Completely removing the B-field within 1 kpc of the GC leads only to a small reduction to the conversion probabilities, as indicated in red. Changing magnetic field models to that of [47] (PTKN11), while also removing the B-field within the inner kpc, leads to slightly enhanced conversion probabilities, as shown in orange (for both the YMW16 and ne2001 n e models). Note that the conversion probabilities exhibit clear constructive and destructive interference behavior in this case at low energies, related to the periodic nature of the disk-field component, though including the random field component it is expected that this behavior would be largely smoothed out.
As discussed previously the magnetic field is expected to be significantly larger closer in towards the GC than in our fiducial B-field model. As an illustration in blue we show the conversion probabilities computed, from the two different free-electron models, when we only include a B-field component of magnitude 50 µG pointing in theẑ direction within the inner 400 kpc (explicitly, in this case we do not include any other B-field model outside of the inner 400 kpc). The conversion probabilities are enhanced in this case by about an order of magnitude across most energies relative to in our fiducial model. The inner Galaxy also likely contains localized regions of even strong field strengths, such as non-thermal filaments with ∼mG ordered fields. As an illustration of the possible effects of such fields on the conversion probabilities, in Fig. S8 we show in grey the result we obtain for the conversion probability when we assume that the axions traverse the GC radio arc, which we model as a 10 kpc wide region with a vertical field strength of 3 mG and a free-electron density n e = 10 cm −3 [55,77]. Due to modeling uncertainties in the non-thermal filaments and the ambient halo field in the inner hundreds of pc, we do not include such magnetic-field components in our fiducial conversion probability model. However, we stress that in the future, with a better understanding of the Galactic field structure in the inner Galaxy, our results could be reinterpreted to give stronger constraints.
The Wd1 conversion probabilities change by over an order of magnitude going between the JF12 and PTKN11 models, as seen in the bottom left panel of Fig. S8, though it is possible that this difference would be smaller when random fields are properly included on top of the JF12 model (though again, we chose not to do this because of sensitivity to the random-field domain sizes).
The effects of the different conversion probabilities on the g aγγ limits may be seen in the top right panel for Quintuplet (Arches gives similar results, since the conversion probabilities are the same) and Wd1 in the bottom right panel of Fig. S8. Note that the observed fluxes scale linearly with p a→γ but scale like g 4 aγγ , so differences between conversion probability models result in modest differences to the g aγγ limits. Still, it is interesting to note that the Wd1 limits with the PTKN11 model are stronger than the fiducial Quintuplet limits, which emphasizes the importance of better understanding the B-field profile towards Wd1. For Quintuplet (and also Arches) we see that depending on the field structure in the inner ∼kpc, the limits may be slightly stronger and extend to slightly larger masses (because of field structure on smaller spatial scales) than in our fiducial B-field mode. The conversion probabilities for axions produced in the Quintuplet or Arches clusters for different modeling assumptions for the Galactic magnetic field and free-electron density. Our fiducial result is shown in solid black. Note that the plasma mass, induced by the free-electron density, becomes more important at lower axion energies and induces the lower-energy features. The dashed black curve shows the effect of changing from the YMW16 free-electron model to the ne2001 model. Removing the B-field within the inner kpc leads to the results in red, while only modeling a 50 µG field in the inner 400 pc leads to the results in blue. Changing to the PTKN11 model (and masking the inner kpc) gives the results in orange. We estimate that if the axions traverse the GC radio arc, located near the Quintuplet and Arches clusters, the conversion probabilities could be enhanced to the values shown in grey. (Bottom Left) As in the top left panel but for axions emitted from the Wd1 cluster. (Right Column) The effects of the different conversion probability models on the 95% upper limits on gaγγ for Quintuplet (top right) and Wd1 (bottom right). Note that Arches is similar to Quintuplet, since they are both assumed to have the same conversion-probability models.

II. EXTENDED DATA ANALYSIS RESULTS
In this section we present additional results from the data analyses summarized in the main Letter.

A. Quintuplet
In this subsection we give extended results for the Quintuplet data analysis. Our main focus is to establish the robustness of the flux spectra from the NuSTAR data analysis (shown in Fig. 2) that go into producing the limits on g aγγ shown in Fig. 3.

Data and templates
First we take a closer look at the stacked data and models that go into the Quintuplet data analysis. The stacked counts data in the vicinity of Quintuplet are shown in the left panel of Fig. S9. We show the counts summed from 10 -80 keV. Note that the circle in that figure indicates 2 , which the radius of our fiducial analysis ROI. 1 As in Fig. 1 we also indicate the locations of the individuals stars in Quintuplet that may contribute axion-induced X-ray flux. The middle panel shows the expected background flux from our background template. The template is generally uniform over the ROI, with small variations. On the other hand, the right panel shows the axion-induced signal counts template, normalized for g aγγ = 7 × 10 −12 GeV −1 , which is localized about the center of the SSC. Note that the signal template is generated by accounting for the PSF of NuSTAR in addition to the locations and predicted fluxes of the individual stars.

Axion Luminosity
We now show the axion luminosity and spectra that go into the right panel of Fig. S9. For each star in the cluster, we assign it a set of possible MESA models based on its spectral classification as described in the main text. In the upper left panel of Fig. S10, we show the mean expected axion luminosity, as a function of energy, of the Quintuplet cluster, assuming g aγγ = 10 −12 GeV −1 . The luminosity peaks around 40 keV, but the effective area of NuSTAR, also shown, rapidly drops above 10 keV. Due to the much higher effective area at low energies, most of the sensitivity is at lower energies. There is also considerable flux above 80 keV, although NuSTAR does not have sensitivity at these energies. In the upper right panel, we show the median contribution of each spectral classification in Quintuplet to this luminosity, summed over all stars with the given classification. For all energies of interest, the WC stars dominate the cluster luminosity. This is because WR stars have the hottest cores and there are 13 WC stars in Quintuplet (there is 1 WN star). In the bottom panel, we show the 10 -80 keV luminosity distribution for each spectral classification, along with the 1σ containment bands and the mean expectation. The distribution depends principally on whether or not core helium is ignited while the star is assigned a given classification. The O, BSG, and WNh stars all can be either hydrogen or helium burning, in which case they have 10 -80 keV luminosities of ∼ 10 31 or ∼ 10 33 erg/s, respectively-recall that the jump in temperature during helium ignition is a factor ∼ 3. The LBV phase is always core helium burning, and the star may go supernova in this phase. The same is true of the WR phases WN and WC, although the stars undergoing supernova in this phase are typically more massive.  Table I. The number of stars Nstar for each stellar class in the Quintuplet cluster, along with the predicted axion luminosities (all in erg/s). Note that Quintuplet is ∼30 pc away from the GC. Except in the last column, the axion luminosities are summed over all energies. All entries assume gaγγ = 10 −12 GeV −1 and are summed over all stars for the given stellar class.
The luminosities in Fig. S10 are computed for our fiducial choices of Z = 0.035 and µ rot = 150 km/s. To better understand the importance of these choices we show in Tab. I how the luminosities depend on the initial metallicity Z and mean rotation speed µ rot . Note that each entry in that table shows the luminosity summed over the stellar sub-types (with the number of stars indicated), and except in the two last columns the luminosities are summed over all stars. The uncertainties in the entries in Tab. I come from performing 500 draws from the representative models and account for the variance expected from star-to-star within a given classification. As discussed in the main text, the 10 -80 keV luminosity could be ∼70% larger than in our fiducial model, depending on the initial Z and µ rot .

Injecting an axion signal
As a first test of the robustness of the Quintuplet analysis we inject a synthetic axion signal into the real stacked data and then pass the hybrid real plus synthetic data through our analysis pipeline. Our goal from this test is to ensure that if a real axion signal were in the data with sufficiently high coupling to photons then we would be able to detect it. The results from this test are shown in Fig. S11.  Figure S11. (Left) We inject a synthetic axion signal into the Quintuplet NuSTAR data with axion coupling g inj aγγ . We then pass the hybrid synthetic plus real data through our analysis pipeline and show the best-fit coupling g rec aγγ , along with the recovered 1σ and 2σ uncertainties. (Middle) The discovery TS for the axion signal for the test illustrated in the left panel. The square root of the TS is approximately the discovery significance. (Right) The 95% upper limit recovered for the injected signal test. Importantly, the 95% upper limit is above the injected signal value, for all injected signal strengths, and the upper limit is consistent with the 68% and 95% expectations for the upper limit under the null hypothesis, which are indicated in green and gold, respectively.
The left panel of Fig. S11 shows the best-fit g rec.
aγγ as a function of the simulated g inj. aγγ used to produce the axioninduced counts that are added to the real NuSTAR stacked data. Importantly, as we increase the injected signal strength the recovered signal parameter converges towards the injected value, which is indicated by the dashed curve. Note that the band shows the 68% containment region for the recovered signal parameter from the analysis. As the injected signal strength increases, so to does the significance of the axion detection. This is illustrated in the middle panel, which shows the discovery TS as a function of the injected signal strength. Recall that the significance is approximately √ TS. Perhaps most importantly, we also verify that the 95% upper limit does not exclude the injected signal strength. In the right panel of Fig. S11 we show the 95% upper limit found from the analyses of the hybrid data sets at different g inj aγγ . Recall that all couplings above the g rec aγγ curve are excluded, implying that indeed we do not exclude the injected signal strength. Moreover, the 95% upper limit is consistent with the expectation for the limit under the signal hypothesis, as indicated by the shaded regions at 1σ (green) and 2σ (yellow) containment. Note that we do not show the lower 2σ containment region, since we power-constrain the limits. These regions were computed following the Asimov procedure [66].

Changing region size
As a systematic test of the data analysis we consider the sensitivity of the inferred spectrum associated with the axion model template to the ROI size. In our fiducial analysis, with spectrum shown in Fig. 2, we use an ROI size of r max = 2 . Here we consider changing the ROI size to r max = 1.5 and 2.5 . The resulting spectra are shown in Fig. S12. The spectrum does not appear to vary significantly when extracted using these alternate ROIs, indicating that significant sources of systematic uncertainty related to background mismodeling are likely not at play.

B. Westerlund 1
In this subsection we provide additional details and cross-checks of the Wd1 analysis.  Table II. As in Tab. I but for Wd1.

Data and templates
In Fig. S13 we show, in analogy with Fig. S9, the data, background, and signal maps summed from 15 -80 keV. We note that the background templates are summed using their best-fit normalizations from the fits to the null hypothesis of background-only emission. The signal template is noticeably extended in this case beyond a point-source template and is shown for g aγγ = 8 × 10 −12 GeV −1 and m a 10 −11 eV. The location of the magnetar CXOU J164710.2-45521 is indicated by the red star.

Axion Luminosity
We now show the axion luminosity and spectra that go into the right panel of Fig. S13. In the upper left panel of Fig. S14, we show the mean expected axion luminosity, as a function of energy, of the Wd1 cluster, assuming g aγγ = 10 −12 GeV −1 . In the upper right panel, we show the contribution of each spectral classification in Wd1 to this luminosity, summed over all stars with the given classification. For all energies of interest, the WN stars dominate the cluster luminosity, although the WC stars are important as well. As in Quintuplet, this is due to the fact that WR stars have the hottest cores, but in this case there are more WN stars than WC stars. In the bottom panel, we show the 10 -80 keV luminosity distribution for each spectral classification, along with the 1σ bands and the mean expectation. Again, the more evolved stars produce more axion flux, because their core temperatures increase with time. As in the case of Quintuplet, the O and BSG stars may be pre-or post-helium ignition. The LBV, YHG, and RSG stars are all post-helium ignition, although have generically cooler cores than the WR stars. The WNh stars are entirely helium burning.
In Tab. II we provide detailed luminosities for each of the stellar sub-types for different choices of initial Z and µ rot for Wd1, as we did in Tab. I. Note that we assume Z = 0.035 and µ rot = 150 km/s for our fiducial analysis, even though it is likely that the initial Z is closer to solar (in which case the luminosities would be enhanced, as seen in Tab. II).  Figure S13. As in Fig. S9, but for the Wd1 cluster NuSTAR analysis. The red star indicates the location of the magnetar CXOU J164710.2-45521, which is masked at 0.5'. Also shown is the background-subtracted count data, as in Fig. 1.

Systematics on the extracted spectrum
In analogy to the Quintuplet analysis we may profile over emission associated with the background template to measure the spectrum from 15 -80 keV associated with the axion-induced signal template shown in Fig. S13. That spectrum is reproduced in Fig. S15. For our default analysis we use the ROI with all pixels contained with r max = 2.0 of the cluster center, except for those in the magnetar mask, as indicated in Fig. S13. However, as a systematic test we also compute the spectrum associated with the axion-induced template for r max = 2.5 and 1.5 , as shown in Fig. S15. We measure a consistent spectrum across ROIs at these energies.

C. Arches
In this subsection we present results from the analysis of archival NuSTAR data for an axion-induced signal from the Arches cluster. The Arches cluster is at a similar location, ∼30 pc from the GC, as the Quintuplet cluster. Arches hosts even younger and more extreme (e.g., hotter and more massive) stars than the nearby Quintuplet cluster. Indeed, it is estimated that all ∼105 spectroscopically classified stars within Arches may become core-collapse supernovae within the next ∼10 Myr [78]. A priori, the Arches and Quintuplet clusters should have similar sensitivities to axions, though as we discuss below the axion prediction from Arches is less robust to uncertainties in the initial metallicity than the Quintuplet prediction.

Axion Luminosity
We now describe the axion luminosity and spectra for Arches. In the upper left panel of Fig. S16, we show the mean expected axion luminosity, as a function of energy, of the Arches cluster, assuming g aγγ = 10 −12 GeV −1 . The luminosity peaks at very low energies, although we could not analyze these energies due to contamination from the molecular cloud. As shown by the upper right panel, the Arches luminosity is dominated by the O stars, since the WNh stars are always hydrogen burning with our assumed metallicity of Z = 0.035 and there are many more O stars than WNh stars. In the bottom panel, we show the 10 -80 keV luminosity distribution for the O and WNh stars, along with the 1σ bands and the mean expectation.
However, unlike for the Quintuplet and Wd1 clusters we find that the Arches luminosity is a strong function of the initial metallicity Z, as illustrated in Tab. III. As seen in that table, changing the metallicity from Z = 0.035 to Z = 0.018 increases the flux by over an order of magnitude. This is because at the higher metallicity values the WNh stars are typically not in the He burning phase, while decreasing the initial metallicity slightly causes the WNh stars to enter the He burning phase. Note that at solar initial metallicity (Z = 0.02, and also taking µ rot = 100 km/s) we find that the 10-80 keV flux is 8.7 +9.4 −5.6 × 10 34 erg/s, comparable to but slightly larger than that found for Z = 0.018. Thus, it is possible that the sensitivity of the Arches observations is comparable to that from Quintuplet, but given the larger uncertainties related to the stellar modeling of the Arches stars the limit is, at present, less robust. We stress that the qualitative difference between Arches and Quintuplet that is responsible for this difference is that Quintuplet has a large cohort of WC and WN stars, which are robustly He burning, while Arches does not have any stars in these stellar classes.

Data analysis, results, and systematic tests
We reduce and analyze 370 ks of archival NuSTAR data from Arches. The Arches observations (IDs 40010005001, 40101001004, 40101001002, 40202001002, 40010003001) were performed as part of the same GC survey as the Quintuplet observations as well as for dedicated studies of the Arches cluster below 20 keV. Note that we discard data from the Focal Plane Module B instrument for observations 40101001004, 40101001002, 40202001002, and 40010003001 because of ghost-ray contamination. We perform astrometric calibration using the low-energy data on the Arches cluster itself, which is a bright point source above 3 keV. In the Arches analysis it is known that there is a nearby molecular cloud that emits in hard X-rays [79]. We follow [79] and model emission associated with this extended cloud as a 2D Gaussian centered at R.A.=17 h 45 m 50.62 s , Dec.=−28 • 49 47.17 with a FWHM of 72.4 . The hard X-ray spectrum associated with the molecular cloud has been observed to extend to approximately 40 keV [79]; indeed, we see that including the molecular cloud template, with a free normalization parameter, at energies below 40 keV affects the spectrum that we extract for the axion template, but it does not significantly affect the spectrum extraction above 40 keV. The non-thermal flux associated with the molecular cloud is expected to be well described by a power-law with spectral index Γ ≈ 1.6 and may arise from the collision of cosmic-ray ions generated within the star cluster with gas in the nearby molecular cloud [80]. With this spectral index the molecular cloud should be a sub-dominant source of flux above ∼20 keV, and we thus exclude the 10-20 keV energy range from the Arches analysis, though e.g. including the 15-20 keV bin results in nearly identical results (as does excluding the 20 -40 keV energy range).
The molecular cloud template is illustrated in the bottom left panel of Fig. S17. In that figure we also show the data, background templates, signal template, and background-subtracted counts, as in Fig. S9 for the Quintuplet analysis. Note that we profile over emission associated both the background template and with the halo template when constraining the flux in each energy bin associated with the signal template. As a systematic test of our signal extraction procedure we show in Fig. S18 (left panel) the spectrum extracted for axion emission from the Arches cluster both with and without the halo template. The two spectra diverge below ∼20 keV but give consistent results above this energy. Similarly, we find that the spectrum is relatively insensitive to the ROI size for energies above ∼20 keV, as shown in the right panel of Fig. S18, which is analogous to the Quintuplet Fig. S12.  Figure S18. (Left) The Arches spectrum measured with and without the halo template. Note that we use the spectrum with the halo template in our fiducial analysis, though the difference between the two results is relatively minor above ∼20 keV. (Right) As in Fig. S12 but for the Quintuplet analysis. Note that these spectra are computed while profiling over halo emission. Above ∼20 keV the different ROIs produce consistent results.
In Fig. S19 we show the 95% upper limit we obtain on g aγγ from the Arches analysis, using the conservative modeling with Z = 0.035 and µ rot = 150 km/s. We find no evidence for an axion-induced signal from this search. Note that, as in indicated in Fig. S18, we do not include data below 20 keV in this analysis. CAST Arches 95% limit median expected Figure S19. As in Fig. 3 but from the analysis towards the Arches SSC. No evidence for axions is found from this search.

III. INITIAL METALLICITY DETERMINATION FOR QUINTUPLET AND ARCHES
In our fiducial analysis we assumed the cluster metallicity was Z = 0.035, which we take as the highest allowed metallicity in the Quintuplet cluster. In this subsection we show how we arrived at this value. The cluster metallicity is an important parameter in that it affects the mass loss rates in the stellar winds, the lifetime of individual evolutionary stages, and the surface abundances. Here we use measurements of the nitrogren abundances of WNh stars in the Arches cluster to estimate the uncertainty on the cluster metallicities. The nitrogen abundance during the WNh phase reaches a maximum that depends only on the original CNO content, and as such is a direct tracer of stellar metallicity (and increases with increasing metallicity). Ref. [33] measured the nitrogen abundance in the WNh stars in the Arches cluster at present to be 0.0157 ± 0.0045. We run MESA simulations of the Arches WNh stars on a grid of metallicities from Z = 0.01 to Z = 0.04 and find this measurement implies that the Arches initial metallicity is between Z = 0.018 and Z = 0.035. The results are shown in Fig. S20, where we see that the nitrogen abundance during the WNh phase intersects with the measurement only for the initial metallicities in that range. Although there are no measurements of the Quintuplet WNh nitrogren abundance, note that a similar abundance was found in the nearby GC SSC of 0.0143 ± 0.0042 [81]. Given the similarity of these two measurements, we assume the same metallicity range for Quintuplet as computed for Arches.  [33].