Introduction.— The Standard Model (SM) of particle physics has provided a remarkably successful framework for describing the fundamental particles and their interactions via three of the four known fundamental forces. Despite its predictive power, the SM leaves several key questions unresolved, notably the strong charge-parity (CP) problem in quantum chromodynamics (QCD) and the composition of dark matter. Intriguingly, a single theoretical extension proposed over four decades ago by Peccei and Quinn [1] offers a potential solution to both. This framework predicts the existence of a light, pseudo-scalar boson—the QCD axion [2]–[7]. Beyond this minimal realization, a broader class of particles known as axion-like particles (ALPs) [8] arises in many extensions of the SM. Like QCD axions, ALPs are pseudo-scalar bosons that couple to photons, but unlike QCD axions, they are not constrained by a specific relationship between their mass \(m_a\) and the Peccei-Quinn symmetry-breaking scale \(f_a\). Although ALPs do not solve the strong CP problem, they remain compelling candidates for dark matter. For the purposes of this Letter, we will refer to both QCD axions and ALPs generically as “axions.”

Numerous laboratory and astronomical experiments have been carried out to detect axions  [9]–[39], with additional efforts either ongoing or proposed [31], [40]–[42]. These investigations predominantly focus on the interaction of axions with photons [43], described by the effective Lagrangian term \(\mathcal{L}_{a\gamma} = g_{a\gamma}\, {\boldsymbol{E}}\cdot {\boldsymbol{B}}\, a\),

where \(a\) represents the axion field, \({\boldsymbol{E}}\) and \(\boldsymbol{B}\) denote electric and magnetic fields, respectively, and \(g_{a\gamma}\) is the coupling strength. Such an interaction enables the Primakoff effect, where photons convert into axions (and vice versa) in the presence of an electromagnetic field.

The Sun acts as a bright source of axions, which are generated in its core via the Primakoff process, wherein thermal photons are converted into axions through interactions with the electric fields of electrons and ions in the plasma. These axions can then reconvert into photons within the Sun’s atmospheric magnetic fields, yielding a distinct X-ray flux with spectral and spatial characteristics that allow differentiation from typical solar X-ray emissions and background signals in X-ray telescopes. Advances in solar X-ray observational capabilities, such as higher resolution and sensitivity, have improved prospects for identifying solar-origin axions. Initial efforts using satellite instruments for axion detection include the Yohkoh/SXT  [44] mission data from 1991-2001, utilized in [45] to establish early constraints on the axion-photon coupling. Later missions, including RHESSI [46], [47] and Hinode/XRT [48], applied similar methodologies, gradually refining observational techniques and constraints on axion interactions. However, due to the high background rates and optimization of these instruments for bright solar X-rays, dedicated solar missions are not optimized for weak-signal dark matter searches.

In this letter, we present the first analysis of solar observation X-ray data from NASA’s Nuclear Spectroscopic Telescope Array (NuSTAR) [49], specifically aimed at searching for signatures of solar axion conversion in the solar atmosphere. We find no evidence for an axion-induced signal above the expected background, allowing us to place stringent constraints on the axion-photon coupling over a significant region of parameter space. In particular, we constrain the axion-photon coupling \(g_{a\gamma}\) to less than \(7.3\times 10^{-12}\) GeV\(^{-1}\) at the \(95\%\) confidence level (CL) for axion masses \(m_a \lesssim 4\times 10^{-7}\) eV, substantially improving upon previous limits set by laboratory-based experiments [9]–[39], exceeding the results of planned experiments [31], [40]–[42] and competing with current astrophysical limits  [50], [50]–[68]. Fig. 1 shows the outcome of the present study in comparison to those from previous and forthcoming ground-based axion experiments.

Axion Conversion in Solar Atmospheric Magnetic Field— Accurate theoretical models for the solar axion flux are well-established, and our analysis is grounded in the results presented in [69], [70]. Since both, the generation of axions from X-rays in the solar core and their reconversion into photons in the solar atmosphere’s magnetic field occur via the Primakoff process, the resulting X-ray signal detected by satellite missions scales with \(g_{a\gamma}^{4}\). Axions generated in the core are highly relativistic and follow approximately radial trajectories to the solar surface. The probability of their conversion into photons is then provided by [71], \[\label{eq:probMainnew} P_{a\rightarrow\gamma}(E,h,m_{a}) = \frac{1}{4}g_{a\gamma}^2\Big|\int_0^h dh' B_\perp(h')\;e^{i\int_0^{h'}dh''q(h'',m_{a})}\;e^{-\frac{1}{2}\int_{h'}^h dh''\Gamma(h'')} \Big|^2 ,\enspace\tag{1}\] where \(h\) is the altitude above the visible surface of the Sun, \(B_\perp\) is the component of the magnetic field perpendicular to the direction of propagation of the photon, and \(q(h,m_{a}) = (\omega_{p}^2(h) - m_a^2)/2E\) is the momentum exchanged between the photon in the medium and the axion, both carrying energy \(E\). The plasma frequency of the medium is \(\omega_{p}(h) = \sqrt{e^2 n_e(h)/m_e}\), where \(e\) and \(m_e\) are the electric charge and mass of the electron, respectively, while \(n_e(h)\) is the height-dependent number density of free electrons plus that of bound electrons with ionization energy much smaller than the X-ray energy. Notably, in natural units, the plasma frequency corresponds to an effective photon mass. The term containing \(\Gamma\) represents the absorption between altitudes \(h'\) and \(h\) of the arising X-rays in the solar atmosphere [71], [72].

In this work, we characterize the solar magnetic field by means of advanced and realistic Magneto-Hydrodynamics (MHD) simulations that accurately represent the quiet-Sun atmosphere that corresponds to the conditions for the NuSTAR solar observations used in this letter. For altitudes below 400 km, we make use of magneto-convection simulations with small-scale dynamo activity [73] developed with the MURaM radiative MHD code [74], [75]. This three-dimensional (3D) model applies periodic boundary conditions in the horizontal plane orthogonal to the solar radius vector at the observed solar disk location. It accurately represents the quiet inter-network regions of the solar photosphere, which cover most of the solar disk across the solar activity cycle. The model’s magnetic field, tangled at scales beneath current telescope resolution, shows a mean field strength of 170 Gauss at the visible surface, matching observed depolarization effects in photospheric spectral line polarization [76], [77]. The blue curve in the top panel of Fig. 2 illustrates the altitude variation of \(\langle \lvert B_{\perp} \rvert \rangle\), calculated as the spatial average of \(\lvert B_{\perp} \rvert\) across all points in the plane perpendicular to the solar radius vector.

For modeling the coronal magnetic field, we utilize the publicly available Predictive Science Inc. (PSI) MHD simulation corresponding to the solar eclipse of July 2, 2019 [78], and apply it to the quiet phase of solar activity of our observations. See Supplemental Material (SM) for further details [79]. The PSI simulation provides a 3D representation of the coronal magnetic field up to an altitude of \(h=29 R_\odot\). To obtain \(\langle|B_\perp|\rangle\), we extract the perpendicular component of the magnetic field relative to the Earth-Sun line of sight, averaged over a disk of radius \(0.1R_\odot\) centered on our signal region, as represented by the orange line in the top panel of Fig. 2. We perform a power-law interpolation between the photospheric and coronal magnetic field models (dashed line in the top panel of Fig. 2), which is further validated by the Potential-Field-Source-Surface (PFSS) model for the specific day of our observation [80], denoted by the green line in the top panel of Fig. 2.

To accurately calculate the axion-photon conversion probability in the solar magnetic atmosphere, it is essential to model the free electron and hydrogen density. For this purpose, we utilize the model from [81], which is representative of the quiet Sun. This model builds on the VAL-C chromospheric model [82] and the coronal model [83]. Our focus is specifically on the quiet Sun, excluding the added coronal density contributions from the streamer belt. In addition to free electrons and hydrogen, we account for helium (He) in the plasma density, assuming it follows the density profile of hydrogen (H) with a relative abundance \(n_{\rm{H}} / n_{\rm{He}} = 0.06\) [84]–[86]. Contributions from heavier elements are not included in the computation of \(q\), as their low abundances render their impact on the plasma frequency negligible. In the chromosphere, where the ionization fraction is low, the plasma frequency is primarily determined by the density of bound electrons in atoms, mainly hydrogen, as shown in the bottom panel of Fig. 2. In the corona, with temperatures on the order of \(10^{6}\) K, hydrogen and helium are fully ionized. Given the higher mass of nuclei, their contribution to the plasma frequency is minimal compared to that of free electrons.

The absorption coefficient \(\Gamma\) is calculated as the sum over all species \(i\) such that \(\Gamma_i = n_i\sigma_i\), where \(n_i\) represents the number density of the species \(i\) and \(\sigma_i\) is the total photon cross-section for that species. We include elements up to atomic number \(Z=30\) in our calculation, using elemental abundances from the CHIANTI database  [87], [88], specifically based on the Schmelz extended model of coronal abundances derived in [84]–[86]. For each element, we obtain the total photon cross-section from the NIST XCOM database [89], which accounts for X-ray attenuation due to coherent (Rayleigh) and incoherent (Compton) scattering, as well as photoelectric absorption.

NuSTAR observations.— On February 21, 2020, NuSTAR was directed toward the center of the solar disk, capturing 23,814 seconds of live exposure data with flight module \(A\) and 24,983 seconds with flight module \(B\), to support both a solar axion search and the analysis of X-ray bright points on the solar photosphere [90]. The minimum of the 11-year solar activity cycle provided a unique opportunity to study the emission from the quiet Sun, and this campaign offers the best data for axion studies, since it involved a long dwell at the disk center during a period of very low solar activity, providing an ideal environment by minimizing interference from bright active regions or solar flares within the \(10'\times10'\) field of view [49]. To identify the predicted axion spectrum, counts were accumulated from pixels within \(0.1\) solar radii (\(1.6'\)) of the disk center, while background counts were gathered from an annulus between \(0.15\) and \(0.30\) solar radii.

Although the solar axion surface luminosity profile [91] indicates that an optimal source region for our analysis, given NuSTAR’s energy threshold, would be at approximately \(0.15\,R_{\odot}\), we adopt a more conservative selection radius of \(0.1\,R_{\odot}\) and account for remaining expected axion signal in the background region as discussed later on in the Analysis Section. This choice satisfies two key criteria: it maximizes the axion flux originating from the solar core while minimizing contamination from bright solar features present during the observation period.

For the analysis presented here, NuSTAR conducted a total of nine observations (orbits). Modules \(A\) and \(B\) have been analyzed independently, since time exposure, effective area and detector response are inherent to each particular module. To account for the motion of the solar center in the field of view (FoV), each observation has been divided into four segments, each approximately 700 seconds in duration. The top panels of Fig. 3 present the total counts in the signal region as a function of energy (red dots), along with the background (cyan line) derived from the outer annular region for module \(A\) (left) and module \(B\) (right). The background of the respective modules was adjusted by normalizing to the source area, smoothing via a running polynomial fit, and further correcting for the known gradient of cosmic background X-rays across the detector chips. For clarity, statistical uncertainties of the plotted background were omitted the detector background was measured over an area 5.3 times larger than the signal region, enhancing the accuracy of background estimation. The bottom panel shows the background subtracted spectrum (gray dots) including error propagation for modules \(A\) (left) and \(B\) (right) and we report, for illustrative purposes, the expected axion flux (green line) corresponding to the 95% CL of \(g_{a\gamma}\) for an axion with mass of \(10^{-7}\) eV. Detailed procedures for the NuSTAR solar data processing using NASA-HEASARC HEASoft 6.34, NuSTAR Data Analysis Software (NuSTARDAS v.2.1.4) and calibration database (CALDB v.20240325) are discussed in SM [79]. XSPEC tool [92] version 12.13.1 was used to read the spectra and export the data into ASCII format.

Data Analysis.— The expected photon count resulting from axion conversion within a given NuSTAR energy bin \(i\), originating from region \(j=s,b\) where \(s\) denotes the signal region (\(r<0.1R_\odot\)) and \(b\) the background region (\(0.15R_\odot<r<0.3R_\odot\))is computed separately for each module \(k=A,B\) and observation segment as:

\[N_{\gamma,i}^{j,k}(h,m_{a}) = \Delta t^{k} \int dE\, A_{\rm{eff}}^{j,k}(E)\epsilon_{RMF,i}^{k}\,(E)\frac{dN_a^j}{dA\, dt\, dE}(E)\,P_{a\to\gamma}(E,h,m_{a})\enspace,\]

where \(\Delta t^{k}\) is the observation exposure and \(A_{\rm{eff}}^{j,k}\) the effective area of a given module \(k\). The energy response function of each telescope in an energy bin \(i\) is denoted by \(\epsilon_{RMF,i}^{{k}}\)  [93]. The term \(\frac{dN_a^j}{dA\, dt\, dE}(E)\) corresponds to the differential axion flux originating from the solar core, while \(P_{a\to\gamma}(E,h,m_{a})\) represents the axion-photon conversion probability within the solar atmosphere following Eq. 1 . We calculate the total number of photons from axion conversion expected in each module by summing over all the relative segments and account for leftover signal in \(0.15-0.3R_\odot\). The likelihood analysis was performed using an 80 eV bin width. However, for improved visual clarity, Fig. 3 is presented with a bin width of 200 eV.

The total number of X-ray photons expected in the signal region for a given module and energy bin, is given by \[\lambda^{k}_{i} = N_{\gamma,i}^{s,k} (g_{10}^{4})+ z^{s,k}_{i}\enspace, \label{eq:nevents}\tag{2}\]

where \(z^{s,k}_{i}\) represents the background of bin \(i\) in the signal region and \(g_{10}\equiv g_{a\gamma}/10^{-10}{\,\rm GeV}^{-1}\). To estimate \(z^{s,k}_{i}\), we use the photon counts from the background region. Assuming that the background X-ray photon density is uniform across the solar disk, we can express this as \[z^{s,k}_{i} = (t^{k}_{i}- N_{\gamma,i}^{b,k}) \frac{A_\odot^s}{A_\odot^b} ,\]

where \(t^{k}_{i}\) denotes the total number of photons detected in the bin \(i\) of the background region, \(N_{\gamma,i}^{b,k}\) is the faint contribution from axions within the background annulus. Furthermore, \(A_\odot^s\) and \(A_\odot^b\) represent the areas of the signal and background regions, respectively.

We assume independent Poisson statistics for each energy bin and calculate the likelihood for each module as \[\label{eq:PoisLike} \mathcal{L}^{k} \propto \prod_i \frac{e^{-\lambda_i^{k}} \lambda_i^{n_i^{k}}}{n_i^k!}\enspace.\tag{3}\] In this analysis, \(n_i^k\) denotes the number of photons detected in the \(i\)-bin of the energy spectrum of a given module. To derive our limit, we focus on the energy range from \(4\) to \(11\) keV and apply the prior \(\Theta(g_{10}^4)\), where \(\Theta\) is the Heaviside step function. We then integrate the Bayesian posterior probability density function (PDF) from zero to the point that includes 95% of the total PDF area. To obtain the combined bound shown in Figure 1, we multiply the likelihoods of each module and follow the procedure described above applied to the total likelihood. The two modules yield similar bounds, the one from module \(A\) being slightly more stringent.

Our study considers various sources of systematic error, including the model uncertainties of the solar atmospheric magnetic field, which are the dominant contributors to the systematics, as well as uncertainties in the solar axion flux and the background of the NuSTAR detector. Overall, the combined effect of these uncertainties in our bound is calculated to be \(_{-21.2\%}^{+27.6\%}\), which is further discussed in SM [79].

Results and Discussion.— Using the statistical analysis outlined above, the NuSTAR solar axion observation data provides a constraint of \(g_{a\gamma}\lesssim 7.3\times 10^{-12}\) GeV\(^{-1}\) at 95% CL on the axion-photon coupling strength \(g_{a\gamma}\) for axion masses \(m_a \lesssim 4\times 10^{-7}\) eV. This limit is depicted in blue in the broader \(m_{a}-g_{a\gamma}\) parameter space shown in Fig. 1. While NuSTAR’s limits remain nearly constant for \(m_a \lesssim 4\times 10^{-7}\) eV, they weaken for higher axion masses due to the loss of coherence in the conversion probability, as described by Eq. 1 , when the axion mass exceeds the plasma frequency at lower solar atmospheric heights. At low axion masses, the \(\rm{NuSTAR}\) limit represents a significant improvement over the leading CAST constraint [10] (light red) and falls between the projected limits of future experiments, such as \(\rm{BabyIAXO}\) [31] (dashed-dotted, with data collection scheduled to begin in 2028) and IAXO [41] (dashed, expected in the mid 2030s). For comparison, we also report the projected sensitivity of the ALPS-II laser propagation experiment at DESY [40] (dotted line) in Fig.1. We note that the NuSTAR analysis presented here explores the \(m_{a}-g_{a\gamma}\) parameter space with unprecedented sensitivity over a broad range of axion masses. The region explored in this project includes areas not accessible to any other current or near future axion experiments and probes large regions that could address the transparency hints [94].

Unlike haloscope experiments [11]–[39], helioscopes, including our study, do not assume axions to be dark matter, as the axion flux is inherently produced by the Sun. To maintain consistency, the luminosity associated with the axion flux must remain a small fraction of the total solar luminosity, to avoid significant energy loss due to escaping axions [95]. The measured solar neutrino flux and helioseismology observations place a constraint of \(g_{10} \lesssim 4.1\) [96], which is well above the bounds reported here. In this regime, the theoretical precision of the axion flux is at the few-percent level [97], as outlined above, and therefore does not pose a significant systematic concern.

Acknowledgments.— We gratefully acknowledge the support of the NuSTAR operations, software, and calibration teams in the execution and analysis of these observations. Our thanks also extend to Georg G. Raffelt, Jiří Štěpan and T. O’Shea for insightful discussions. ET, MR, and MT are supported by the project “Theoretical Astroparticle Physics (TAsP)" funded by INFN and the”Grant for Internationalization" from the University of Torino. MT further acknowledges the research grant “Addressing Systematic Uncertainties in Searches for Dark Matter No. 2022F2843L" funded by MIUR. ET expresses gratitude to the University of Zaragoza for its hospitality during the initial phases of this work. This publication is based on work supported by COST Action COSMIC WISPers CA21106.

The code to generate the results of this study is openly available at [98].

1 Modeling of the Axion Flux and Solar Atmosphere.↩︎

1.1 Solar axion flux↩︎

In the core of the Sun, black-body photons can convert to axions in the presence of the electric fields generated by charged particles in the high-temperature plasma. This process follows the interaction \(\gamma+(e^{-}, Ze)\rightarrow a + (e^{-}, Ze)\) [43]. According to current models of stellar nucleosynthesis [99], nuclear fusion reactions within the Sun’s core serve as the primary energy source, with the proton-proton chain being a key mechanism. The core, characterized by a temperature of \(T_{c}\sim 1.3\,\rm{keV}\) and a density of \(\rho_{c}\sim 1.5 \times 10^{5}\,\rm{kg/m^{3}}\), occupies about 20% of the solar radius and is the only region where significant heat production occurs. Axion emission from the solar core is highly temperature-dependent, with an average energy of \(E_a\sim 2.7 \times T_c\), aligning with the temperature dependence of thermal X-ray production. Given the temperature gradient within the solar interior, the axion spectrum exhibits radial variation, as illustrated in Fig. 4. This model allows calculation of the total differential axion flux from the core, incorporating the standard solar model and the Sun’s luminosity \(L_{\odot}=3.85\times 10^{33}\,\rm{erg\,s^{-1}}\) [69], [100].

The differential solar X-ray flux resulting from axion-photon conversion in the Sun’s atmosphere depends on two key factors: the axion flux generated within the solar interior, and the probability of axions converting into photons upon passing through the solar atmosphere, denoted as \((P_{a\rightarrow\gamma})\), and can be written as: \[\label{eq:flux} \frac{dN_{\gamma}}{dE\,dA\,dt\,d\Omega}=\frac{dN_{a}}{dE\,dA\,dt\,d\Omega}\,P_{a\rightarrow\gamma}\enspace,\tag{4}\] with \(dE,\) \(dA\), \(dt\), and \(d\Omega\) being the elements of energy, area, time, and solid angle, respectively. As shown in Fig. 4, the resulting axion flux peaks at energies \(\sim 4{\, \rm keV}\) in the solar core. These features are imprinted on the spectral and morphological properties of the resulting X-ray emission in Eq. 4 . The axion flux from Primakoff conversion in the Sun’s core [100] is shown on the top panel of Fig. 4, where \(A_\odot\) represents the transverse area of the solar disk expressed in units of solar radius squared \(R_\odot^2\). The bottom panel of Fig. 4 shows the total axion flux and the axion flux arising from our signal and background regions.

1.2 Solar magnetic field↩︎

To accurately model axion conversion in the solar atmosphere, we require detailed knowledge of the magnetic field structure.

As explained in the main manuscript, the magnetic field in the quiet solar photosphere is obtained from [73]. Among the various magneto-convection 3D models, we selected the one with a mean field strength of 170 G at the model’s visible surface. This selection was made because it is the only model consistent with observational constraints on scattering polarization signals in the Sr I 460.7 nm line [78], see [77].

The coronal magnetic field can be described by magnetohydrodynamic (MHD) models based on photospheric magnetic field observations. We used the Predictive Science Inc. (PSI) model for the magnetic structure during the July 2, 2019, total solar eclipse, as detailed in [78]. Although this model represents an epoch several months before the NuSTAR observations on February 21, 2020, it provides a relevant baseline, despite some solar evolution over that period (see the decrease in solar activity evident in the two snapshots from the Hinode X-Ray Telescope shown in Fig.5). During the eclipse, an isolated active region was Earth-facing and represents the primary activity captured in the three-dimensional PSI model under consideration.

From the 3D PSI simulation we define our reference model for the coronal \(B_{\perp}\) as follows. We average \(B_{\perp}\) over the signal region, i.e. the disk within a radius of \(0.1\,R_{\odot}\). We repeat this procedure 120 times, corresponding to a complete solar rotation, with the PSI model rotated by an azimuthal angle of \(3^\circ\) in each iteration. At each altitude \(h\), the median of the resulting 120 values, \(\langle |B_\perp|\rangle\), is computed and corresponds to our baseline model used for the calculation of the axion-photon limit. In addition, we compute the lower 5, 10 and 15 percentiles of the distribution. These magnetic field profiles are shown in Fig. 2.

To validate this procedure, we compare these results with the prediction of the PFSS model [80]. The PFSS model constructs a global magnetic field by assuming a potential field, meaning that it does not incorporate currents threading through the corona as they physically would. Instead, it approximates these effects by setting a boundary condition where the magnetic field becomes purely radial at a distance of \(2.5\,\rm{R_\odot}\). This “source surface" creates artificial surface currents that adjust the internal field to reasonably reflect the influence of photospheric magnetic sources (such as active regions and sunspots) and permits a radial field component, allowing for a representation of the solar wind. These publicly accessible models, although not as advanced as the PSI model, allow us to examine the evolution of the coronal magnetic field between the PSI model date and the NuSTAR observation period. To this end, we compute \(|B_{\perp}|\) using the PFSS model specific for the time of the NuSTAR observation. As shown in Fig. 2 it align closely with the fiducial (median) PSI model described before, strengthening confidence in our results. Instead, the more conservative magnetic field profiles corresponding to the lower 5, 10 and 15 percentiles of the PSI model distribution lie significantly below the PFSS model.

The PSI simulation lacks sufficient resolution to resolve the magnetic field below \(h\sim 0.1R_\odot\), while the model for the photosphere from [73] is reliable only for altitudes \(h<400\) km. This leaves a gap in the precise understanding of the magnetic field between these two regions. To bridge this, we applied a simple power-law interpolation, depicted by the dotted line in the upper part of Fig. 2. This interpolation estimates a perpendicular magnetic field strength of approximately 1.4 G at an altitude of 10,000 km. We notice that this intermediate magnetic field is relevant only for axion masses \(\gtrsim 10^{-4}\) eV – for which the axion-photon conversion process loose coherence in these regions, see Fig. 7.

1.3 Expected X-ray signal↩︎

To calculate the axion-photon conversion probability, we utilize the components outlined in previous sections and perform a numerical integration of \(P_{a\rightarrow\gamma}\) (see Eq. 2 of the main manuscript) up to an altitude of \(h=29 R_{\odot}\), the upper boundary of our model. To illustrate the contributions from various layers within the solar atmosphere, Fig. 7 presents the conversion probability as a function of altitude above the photosphere, assuming conversion halts at altitude \(h\) and the resulting photons are subsequently detected by NuSTAR. This probability is shown for several axion masses at an energy of \(E=4\) keV, with a coupling constant \(g_{a\gamma}=10^{-10}\) GeV\(^{-1}\).

At altitudes \(h < 10^{-4}R_\odot\), the conversion probability remains low due to a lack of coherence (\(qh\gg 1\)), resulting in suppressed X-ray emission. The probability begins to increase as the plasma density decreases, enabling more effective conversion. When the plasma density drops sufficiently and the axion mass is negligible compared to the plasma frequency, the phase \(\varphi = \int dh\, q\) becomes approximately constant, effectively factoring out. Simultaneously, reduced absorption enhances the conversion process. When the plasma frequency reaches the axion mass, i.e., at resonance, \(q= (\omega_{p}^2 - m_a^2)/2E\) changes sign and becomes dominated by the axion mass. Beyond this point, the complex phase oscillation \(\varphi\) becomes significant, particularly for higher axion masses that lead to shorter oscillation wavelengths. This effect, coupled with the decreasing magnetic field, flattens the growth of conversion probability. Since the plasma frequency remains above \(10^{-9}\) eV within the altitude range of our model (see Fig. 2 of the main text), resonance is not met for \(m_a\lesssim 10^{-9}\,{\rm eV}\), and conversion probability flattens as the coronal magnetic field diminishes. For lower axion masses, the photon production rate, \(dP/dh\), peaks in the mid-to-upper corona. Photon production at even greater heights up to 1 AU is conservatively neglected.

Fig. 8 presents the expected photon flux from our signal region for different axion masses. For the largest axion mass shown, \(m_a = 4\times 10^{-3}\) eV, the effect of photoelectric absorption is clearly apparent. Specifically, increased absorption is observed for X-ray photons with energies above the K-edge of iron at approximately \(\sim 7.1\) keV. However, for the axion masses where the coherence is completely fulfilled, \(m_{a}\lesssim 5\times 10^{-4}\) eV , the presence of K-edges becomes negligible.

Fig. 9 shows the expected photon flux from our signal region \(R<0.1 R_\odot\). We can see that, as expected, the signal peaks at around 4 keV. Fig. 7 shows that for masses larger than \(10^{-6}\) eV, the conversion probability stops growing at lower and lower altitude in the solar atmosphere, leading to a reduced photon flux. This explains the drop in the photon production for \(m_a \gtrsim 10^{-6}\) eV seen in Figure 9.

2 NuSTAR’s Solar Observation Campaigns↩︎

Since 2014, NuSTAR has conducted several observations of the Sun, with an overview of these observations and quick-look plots available in [103]. These observations have primarily focused on weakly flaring active regions, often located away from the disk center and therefore, not ideal for axion searches. Even when NuSTAR targets quieter regions of the solar disk, stray X-rays—referred to as “ghost rays"—can still be detected hindering low-event studies. These rays result from a single reflection in the two-bounce optics and can originate from active regions outside the \(10'\times10'\) field of view [104]. As a result, for axion searches, the entire solar disk must be relatively quiet. The minimum of the 11-year solar activity cycle provided a unique opportunity to study the emission from the quiet Sun, and NuSTAR conducted several campaigns targeting this period from 2018 to 2020. Among these, the 21-Feb-2020 campaign offers the best data for axion studies: This campaign involved a long dwell at the disk center during a period of very low solar activity, with 9 orbits of observations yielding approximately 23.8 ks and 24.9 ks of data from NuSTAR’s telescope modules \(A\) and \(B\), respectively. Despite the quiet conditions, faint X-ray features, known as X-ray bright points, are present in the solar atmosphere. A detailed analysis of these features is presented in [90], where their spectra are found to be consistent with optically-thin thermal emission sources. The effective temperature of these emissions is typically \(\leq\)​3 MK, with occasional brief brightening reaching closer to 4 MK. Such sources produce X-ray spectra that steeply decrease with photon energy, making them nearly undetectable above 4 keV amidst the background noise. While other NuSTAR observations from this solar minimum exist, they were part of mosaics covering the entire solar disk, contributing only about 100 extra seconds of data.

With the next solar minimum expected in 2030, the current dataset represents the best X-ray satellite observation for solar axion analysis in the near future.

2.1 NuSTAR data processing for solar observations↩︎

NuSTAR collected quiet Sun data over nine spacecraft orbits between 04:51 GMT and 22:41 GMT on February 21, with each orbit providing about one hour of observation. The \(10'\) square field of view (FoV) covered \(\sim18\%\) of the solar disk. While the satellite’s orbit was selected to minimize the impact of passage through the South Atlantic Anomaly (SAA) [49], most observations were affected by observational gaps or dead times due to such passages. The SAA is a geomagnetic anomaly where the weakened magnetic field over the southern Atlantic Ocean exposes spacecraft in low-Earth orbit, like \(\rm{NuSTAR}\), to Earth’s inner radiation belt, extending down to the upper atmosphere [105]. As the Sun drifted slowly across the FOV during the observations, we adjusted the coordinates of the solar center in 10 to 15-minute intervals (segments) to minimize the effects of this drift. Separate data files were extracted for each orbit segment and the details of each orbit can be found in Table 1.

The data were processed using NuSTAR Data Analysis Software (NuSTARDAS) version 2.1.4 to generate scientific products. NuSTARDAS is integrated within the NASA-HEASARC HEASoft software framework (version 6.34), which was downloaded along with the CALDB calibration database (version 20240325) [106]. XSPEC tool [92] version 12.13.1 was used to read the spectra and export the data into ASCII format. The data processing involved three main steps: calibration, screening, and product extraction. All 32 grades of NuSTAR data were collected without filtering, with each grade corresponding to a distinct pattern of energy deposition across the detector pixels [106]. This resulted in level 0 data, which consists of raw telemetry packets that were then formatted into Level 1 FITS files prior to calibration. We processed the Level 1 data using the software module nupipeline to generate Level 1a calibrated files and Level 2 calibrated, cleaned files.

The calibration process includes several steps, such as nuhotpix (for identifying hot pixels), nucalcpha (for assigning energy and grade to each event), and nucoord (for converting raw coordinates into detector and sky coordinates). Screening, also referred to as cleaning, involves additional processes like nucalcsaa (for calculating SAA passages), nulivetime (for correcting dead times in the event files), and nusplitsc (for splitting the cleaned event files into separate files based on the spacecraft’s switching between cameras used for attitude determination).

The initial processing was followed by selection of the source and background regions. The source region is a circle centered at the solar disk center with radius \(0.1 R_{\odot}\) (1.6’), and the background region is an annulus from \(0.15\) to \(0.30\) solar radius. Wedges were occasionally excluded from the background annulus where examination of the energy range from 3.5-4.0 keV showed an X-ray bright point. Generally, no obvious excess was visible above 4.0 keV, but we nevertheless removed these regions in case their aggregate contained a slight excess not clearly visible in any one image. This method avoids bias by not using the data that actually contribute to the analysis (\(>\)​4.0 keV) to decide what areas to excise. The editing of the annulus region was done separately for each 10–15 min sub-interval orbit. Fig. 10 shows the signal and background regions superimposed on data for all the time intervals of one orbit.

The first part of the entire nuproducts product extraction sequence includes the following functions: nuproducts extracting refined PHA files and light curves by slicing out our chosen time intervals and the shapes of our supplied region files; ‘numkarf’, creating ancillary response files (ARFs) which detail the telescope modules’ responses as a function of position and energy; and ‘numkrmf’, creating redistribution matrix files (RMFs) which detail the detectors’ response as a function of photon energy. We executed nuproducts for each source and background region by supplying directories, IDs, event files, region files, time interval files, and setting automatic background scaling and extraction to ‘no’. The second part of the nuproducts sequence is called ‘nubackscale’, where we extracted our own background scaling variables separately. These variables define what proportion our source and background regions are of the entire FOV without hot/bad pixels and detector gaps. We executed ‘nubackscale’ by supplying the PHA files of both the source and background from the previous step as well as other reference files pertaining to the observation ID. This added a  keyword to new PHA files, which we used in follow-on steps.

“Stray light" (cosmic diffuse X-rays that enter the detector without encountering the optics) is our dominant source of background and has a distinctive gradient pattern across the detector plane, based on variable shadowing by the spacecraft structure and optics [107]. Using the nuskybgd package, we determined the effect of this pattern on each of our source and background annulii regions. The XSPEC software tool (version 12.13.1) was used as well to obtain spectra.

3 Sensitivity of the results to model uncertainties. Systematic studies↩︎

The sources of systematic uncertainty in our bounds include the solar magnetic field, the solar axion flux, and the NuSTAR background. To assess how these factors influence our results, we first explore how the bounds vary with changes in the profile shape and normalization of the coronal magnetic field model. Our photospheric model is already a conservative estimate, as it represents the inter-network regions of the quiet Sun, excluding contributions from network and active regions. To estimate the uncertainty arising from the shape of the coronal magnetic field profile, we recompute our bound in the limit of vanishing axion mass using 120 profiles, each obtained by rotating the model by \(3^\circ\) increments (see section on the solar magnetic field) and to isolate the dependence on the profile shape, we rescale each profile so that the average magnetic field strength \(\langle|B_\perp|\rangle\) at an altitude of \(h=0.1R_\odot\) matches that of our fiducial model. This also aligns with the PFSS model for the day of observation. We then recompute the bound for each of the 120 profiles and report the systematic uncertainty in Table 2, expressed as the percentage variation in the bound.

The uncertainty arising from the normalization of the coronal magnetic field is estimated from the variance of the PFSS model over the perpendicular disk of radius \(0.1R_\odot\) at an altitude of \(h=0.1R_\odot\). In our reference model, the value of \(\langle|B_\perp|\rangle\) at this altitude is 0.174 G. From the PFSS model, we estimate the \(1\sigma\) range to be \(0.155~\mathrm{G}<|B_\perp|<0.190~\mathrm{G}\). The variance at this altitude is larger than that at higher altitudes (\(h>0.1R_\odot\)), and \(h=0.1R_\odot\) is the lowest altitude at which we use the coronal model from PSI. To estimate the corresponding uncertainty, we rescale our fiducial magnetic field by an overall factor so that its value at \(h=0.1R_\odot\) is either 0.155 G or 0.190 G. The resulting change in the bound is reported in Table 2 as a percentage variation.

Another source of uncertainty is the solar axion flux, that has been recently revisited in [97], estimating a flux uncertainty of 6% due to the solar model. This rescaling of the axion flux has a minor effect on our bounds, as shown in Table 2. Additionally, we considered the potential impact of the Cosmic X-ray Background (CXB) as a possible source of background variability. The dominant background component at the energies where we expect an axion signal is the aperture Cosmic X-Ray Background (aCXB), which could present spatial variations across the NuSTAR detectors, potentially leading to an imperfect subtraction of this component in our analysis. Based on previous systematic studies of the NuSTAR instrument [108], we account for background uncertainties of \(\pm 1.5\%\). The final line of Table 2 shows how our bound for \(m_a\to 0\) is affected by this uncertainty in the photon counts within the annular region if the background estimation is adjusted by this factor.

4 Comparison to astrophysical bounds↩︎

In addition to the experimental limits presented in the main body, there are several other competitive, yet more model-dependent, astrophysical bounds on the axion-to-photon coupling \(g_{a\gamma}\) that have been reported in the literature. These limits are shown in Fig. 11.

Studies in Refs. [50]–[56] focus on detecting the imprint of axion-photon mixing in the energy spectra of extra-galactic gamma-ray and X-ray sources. Other notable targets for constraining axion-photon conversion include super-giant stars [57], stars in external galaxies [58], magnetic white dwarfs [59], [60], pulsars [61], super star clusters [62], and the supernova SN1987A [63], [64]. However, although these systems provide excellent opportunities for axion searches, they come with the hurdle of estimating the size of uncertainties for the intergalactic magnetic field and plasma density models, which significantly affect the derived bounds.

Other astrophysical bounds shown in Fig. 11, are derived from neutron stars [65], [66] or the explosion of hypothetical axion stars [67] and, are based on the assumption that all dark matter consists of axions. In contrast, the limits obtained in this work, as well as those discussed earlier, do not rely on such an assumption.

Lastly, we also present in Fig. 11 limits from globular clusters, derived from the R parameter [113], which measures the ratio of horizontal branch stars to red giant stars. These bounds are comparable to those obtained from the \({\rm R}_2\) parameter [114], which measures the ratio of stellar populations on the asymptotic giant branch to those on the horizontal branch.

As illustrated in Fig. 11, many of the methods discussed in this section show promising potential for axion detection. However, the associated limits often come with larger uncertainties than those derived in this study. As mentioned earlier, these uncertainties arise primarily from the complexity of modeling axion production and conversion in the far away astrophysical objects and intergalactic environments studied. In this context, the search presented in this paper offers a distinct advantage: The Sun provides an environment that is studied in great detail, allowing for a robust assessment of the model uncertainties, and leading to more reliable bounds on axions.

References↩︎