X-ray Signatures of Axion Conversion in Magnetic White Dwarf Stars

TL;DR

This paper proposes that axions produced in white dwarf stars can convert into X-rays in their magnetic fields, and uses observations to set new constraints on axion properties, enhancing the potential to test axion-related stellar cooling explanations.

Contribution

It introduces a novel method to detect axions via X-ray signatures from magnetic white dwarfs and provides the strongest current constraints on axion couplings using existing X-ray data.

Findings

Set the strongest constraints on axion-electron and axion-photon couplings.

Showed that dedicated X-ray observations could significantly improve sensitivity.

Proposed a new observational approach to test axion-like particles in stellar environments.

Abstract

White dwarf (WD) stars may radiate keV-energy axions produced in their stellar cores. This has been extensively studied as an extra channel by which WDs may cool, with some analyses even suggesting that axions can help explain the observed WD luminosity function. We show that the radiated axions may convert into X-rays in the strong magnetic fields surrounding the WDs, leading to observable X-ray signatures. We use Suzaku observations of the WD RE J0317-853 to set the strongest constraints to date on the combination of the axion-electron ($g_{aee}$) times axion-photon ($g_{a\gamma\gamma}$) couplings, and we show that dedicated observations of magnetic WDs by telescopes such as Chandra, XMM-Newton, and NuSTAR could increase the sensitivity to $|g_{aee} g_{a\gamma\gamma}|$ by over an order of magnitude, allowing for a definitive test of the axion-like-particle explanation of the stellar cooling anomalies.

Tables

Table 1. Table 1: MWD stars that make good candidates for measurement of their secondary, axion-induced X 𝑋 X -ray flux. The columns correspond to the star’s mass in solar masses, radius in solar radii, luminosity in solar luminosities, effective temperature in Kelvin, magnetic field strength in mega-Gauss, distance from Earth in parsecs, and predicted X 𝑋 X -ray flux from 2 − 10 ​ keV 2 10 keV 2-10\,\mathrm{keV} in erg / cm 2 / s erg superscript cm 2 s \mathrm{erg}/\mathrm{cm}^{2}/\mathrm{s} , calculated assuming m a = 10 − 9 ​ eV subscript 𝑚 𝑎 superscript 10 9 eV m_{a}=10^{-9}\,\mathrm{eV} and g a ​ γ ​ γ ​ g a ​ e ​ e = 10 − 24 ​ GeV − 1 subscript 𝑔 𝑎 𝛾 𝛾 subscript 𝑔 𝑎 𝑒 𝑒 superscript 10 24 superscript GeV 1 g_{a\gamma\gamma}g_{aee}=10^{-24}\ \mathrm{GeV}^{-1} . The parameters were obtained by merging the catalogs in Refs.  Kleinman  et al. ( 2013 ); Ferrario  et al. ( 2015 ); Brown  et al. ( 2018 ) . We infer the mass and radius of WD 2010+310 as discussed in the text.

	$M_{\text{WD}}\left[M_{\odot }\right]$	$R_{\text{WD}}\left[R_{\odot }\right]$	$L_{\gamma }\left[L_{\odot }\right]$	$T_{\mathrm{eff}}\left[K\right]$	$B\left[\mathrm{MG}\right]$	$d_{\text{WD}}\left[\mathrm{pc}\right]$	$F_{2-10}\left[\mathrm{erg}/{\mathrm{cm}}^2/s\right]$
RE J0317-853	$1.32$	$0.00405$	$0.0120$	$30000$	$200$	$29.54$	$6.8\times {10}^{-14}$
WD 2010+310	$1^{\ast }$	${0.00643}^{\ast }$	$0.00566$	$19750$	$520$	$30.77$	$4.4\times {10}^{-14}$
WD 0041-102 (Feige 7)	$1.05$	$0.00756$	$0.00635$	$18750$	$35$	$31.09$	$3.0\times {10}^{-14}$
WD 1031+234	$0.937$	$0.00872$	$0.0109$	$20000$	$200$	$64.09$	$2.3\times {10}^{-14}$
WD 1533-057	$0.717$	$0.0114$	$0.0121$	$18000$	$31$	$68.96$	$1.3\times {10}^{-14}$
WD 1017+367	$0.730$	$0.0111$	$0.0082$	$16500$	$65$	$79.24$	$7.1\times {10}^{-15}$
WD 1043-050	$1.02$	$0.00787$	$0.00388$	$16250$	$820$	$83.33$	$5.4\times {10}^{-15}$
WD 1211-171	$1.06$	$0.00754$	$0.00992$	$21000$	$50$	$92.61$	$5.4\times {10}^{-15}$
SDSS 131508.97+093713.87	$0.848$	$0.00968$	$0.01347$	$20000$	$14$	$101.7$	$3.5\times {10}^{-15}$
WD 1743-520	$1.13$	$0.00681$	$0.00184$	$14500$	$36$	$38.93$	$2.9\times {10}^{-15}$

Table 2. Table 2: The X 𝑋 X -ray flux of RE J0317-853 is calculated for different values of M WD subscript 𝑀 WD M_{\text{WD}} , R WD subscript 𝑅 WD R_{\text{WD}} , and T eff subscript 𝑇 eff T_{\mathrm{eff}} , which were determined by  Külebi  et al. ( 2010 ) . The luminosities, L γ subscript 𝐿 𝛾 L_{\gamma} , are inferred from the Stefan-Boltzmann law. We calculate F 2 − 10 subscript 𝐹 2 10 F_{2-10} for m a = 10 − 9 ​ eV subscript 𝑚 𝑎 superscript 10 9 eV m_{a}=10^{-9}\,\mathrm{eV} and g a ​ γ ​ γ ​ g a ​ e ​ e = 10 − 24 ​ GeV − 1 subscript 𝑔 𝑎 𝛾 𝛾 subscript 𝑔 𝑎 𝑒 𝑒 superscript 10 24 superscript GeV 1 g_{a\gamma\gamma}g_{aee}=10^{-24}\ \mathrm{GeV}^{-1} .

	$M_{\text{WD}}\left[M_{\odot }\right]$	$R_{\text{WD}}\left[0.01R_{\odot }\right]$	$L_{\gamma }\left[L_{\odot }\right]$	$T_{\mathrm{eff}}\left[K\right]$	$B\left[\mathrm{MG}\right]$	$d_{\text{WD}}\left[\mathrm{pc}\right]$	$F_{2-10}\left[\mathrm{erg}/{\mathrm{cm}}^2/\sec \right]$
CO-low-T	$1.32\pm 0.020$	$0.405\pm 0.011$	$0.0120$	$30000$	$200$	$29.54$	$6.8\times {10}^{-14}$
CO-high-T	$>1.46$	$0.299\pm 0.008$	$0.0503$	$50000$	$200$	$29.54$	$2.4\times {10}^{-13}$
ONe-low-T	$1.28\pm 0.015$	$0.416\pm 0.011$	$0.0126$	$30000$	$200$	$29.54$	$7.2\times {10}^{-14}$
ONe-high-T	$1.38\pm 0.020$	$0.293\pm 0.008$	$0.0483$	$50000$	$200$	$29.54$	$2.2\times {10}^{-13}$

Equations

$$\begin{array}{rl}\frac{L_a}{L_{\odot }}& \approx 1.6\times 10^{-4}{\left(\frac{g_{aee}}{10^{-13}}\right)}^2\left(\frac{M_{\text{WD}}}{1M_{\odot }}\right){\left(\frac{T_c}{10^7\text{K}}\right)}^4,\end{array}$$

$$\begin{array}{r}\left[i{\partial }_r+E+\left(\begin{array}{cc}{\Delta }_{\parallel }& {\Delta }_B\\ {\Delta }_B& {\Delta }_a\end{array}\right)\right]\left(\begin{array}{c}{A}_{\parallel }\\ a\end{array}\right)=0,\end{array}$$

$$\begin{split}p_{a\to\gamma}=\Bigl{|}\int_{R_{\text{WD}}}^{\infty}\!\mathrm{d}r^{\prime}\Delta_{B}(r^{\prime})\,e^{i\Delta_{a}r^{\prime}-i\int_{R_{\text{WD}}}^{r^{\prime}}\!\mathrm{d}r^{\prime\prime}\Delta_{\parallel}(r^{\prime\prime})}\Bigr{|}^{2}\,.\end{split}$$

$$\begin{array}{r}\frac{d{F}_{{\gamma }_a}}{dE}(E)=\frac{d{L}_a}{dE}(E)\times p_{a\to \gamma }(E)\times \frac{1}{4\pi d_{\text{WD}}^2},\end{array}$$

$$\begin{split}T_{c}\simeq\bigl{(}3\times 10^{6}\,\,{\rm K}\bigr{)}\left(\frac{L_{\gamma}}{10^{-4}L_{\odot}}\right)^{0.4}\,.\end{split}$$

$$\mathcal{L}$$

$$\displaystyle\qquad-A_{\mu}j^{\mu}-\frac{1}{4}g_{a\gamma\gamma}\,a\,F_{\mu\nu}\tilde{F}^{\mu\nu}+\frac{\alpha_{\mathrm{em}}^{2}}{90\,m_{e}^{4}}\Bigl{[}\bigl{(}F_{\mu\nu}F^{\mu\nu}\bigr{)}^{2}+\frac{7}{4}\bigl{(}F_{\mu\nu}\tilde{F}^{\mu\nu}\bigr{)}^{2}\Bigr{]}\,,$$

$$\ddot{a}-{\nabla }^{2}a+m_a^{2}a$$

$$-\nabla \cdot \dot{\mathbf{A}}$$

$$\displaystyle\qquad+\frac{16\alpha_{\mathrm{em}}^{2}}{45\,m_{e}^{4}}{\bm{\nabla}}\cdot\Bigl{[}\frac{1}{2}\bigl{(}|\dot{{\bm{A}}}|^{2}-|{\bm{\nabla}}\times{\bm{A}}|^{2}\bigr{)}\dot{{\bm{A}}}+\frac{7}{4}\dot{{\bm{A}}}\cdot({\bm{\nabla}}\times{\bm{A}})({\bm{\nabla}}\times{\bm{A}})\Bigr{]}$$

$$\ddot{\mathbf{A}}-{\nabla }^2\mathbf{A}+\nabla (\nabla \cdot \mathbf{A})+m_A^2\mathbf{A}$$

$$\displaystyle\qquad-\frac{16\alpha_{\mathrm{em}}^{2}}{45\,m_{e}^{4}}\partial_{t}\Bigl{[}\frac{1}{2}\bigl{(}|\dot{{\bm{A}}}|^{2}-|{\bm{\nabla}}\times{\bm{A}}|^{2}\bigr{)}\dot{{\bm{A}}}+\frac{7}{4}\dot{{\bm{A}}}\cdot({\bm{\nabla}}\times{\bm{A}})({\bm{\nabla}}\times{\bm{A}})\Bigr{]}$$

$$\displaystyle\qquad-\frac{16\alpha_{\mathrm{em}}^{2}}{45\,m_{e}^{4}}{\bm{\nabla}}\times\Bigl{[}\frac{1}{2}\bigl{(}|\dot{{\bm{A}}}|^{2}-|{\bm{\nabla}}\times{\bm{A}}|^{2}\bigr{)}({\bm{\nabla}}\times{\bm{A}})-\frac{7}{4}\dot{{\bm{A}}}\cdot({\bm{\nabla}}\times{\bm{A}})\dot{{\bm{A}}}\Bigr{]}\,.$$

$$\rho (\mathbf{x},t)=\bar{\rho }(\mathbf{x}),\mathbf{j}(\mathbf{x},t)=\bar{\mathbf{j}}(\mathbf{x}),a(\mathbf{x},t)=\bar{a}(\mathbf{x})+\delta a(\mathbf{x},t),\text{and}\mathbf{A}(\mathbf{x},t)=\bar{\mathbf{A}}(\mathbf{x})+\delta \mathbf{A}(\mathbf{x},t),$$

$$\delta a(\mathbf{x},t)$$

$$\text{with}a(z)=\tilde{a}(z)e^{-i\omega t+ikz}$$

$$\delta \mathbf{A}(\mathbf{x},t)$$

$$\text{with}{A}_a(z)={\tilde{A}}_a(z)e^{-i\omega t+ikz}.$$

$$\left[i{\partial }_z+k+\left(\begin{array}{ccc}{\Delta }_x& {\Delta }_{xy}& {\Delta }_{ax}\\ {\Delta }_{xy}& {\Delta }_y& {\Delta }_{ay}\\ {\Delta }_{ax}& {\Delta }_{ay}& {\Delta }_a\end{array}\right)\right]\left(\begin{array}{c}{A}_x(z)\\ A_y(z)\\ a(z)\end{array}\right)=0,$$

$${\Delta }_a$$

$${\Delta }_x$$

$${\Delta }_y$$

$${\Delta }_{ax}$$

$${\Delta }_{ay}$$

$${\Delta }_{xy}$$

$$\begin{array}{r}{p}_{a\to \gamma }={\sin }^{2}2\theta {\sin }^2({\Delta }_{\mathrm{osc}}{R}_{\text{WD}}/2),\end{array}$$

$$\begin{array}{rl}\tan 2\theta & =\frac{2{\Delta }_B}{{\Delta }_{\parallel }-{\Delta }_a},{\Delta }_{\mathrm{osc}}=\frac{{\Delta }_{\parallel }-{\Delta }_a}{\cos 2\theta }.\end{array}$$

$$\begin{split}p_{a\to\gamma}&\approx\frac{(\Delta_{B,0}R_{\text{WD}})^{2}}{(\Delta_{\parallel,0}R_{\text{WD}})^{\frac{4}{5}}}\Bigl{|}\frac{\Gamma(\frac{2}{5})-\Gamma(\frac{2}{5},-\frac{i}{5}\,\Delta_{\parallel,0}R_{\text{WD}})\bigr{]}}{5^{\frac{3}{5}}}\Bigr{|}^{2}\,,\end{split}$$

$$\displaystyle\bar{\bm{B}}=\frac{|{\bm{m}}|}{4\pi}\,\frac{1}{r^{3}}\,\Bigl{(}3(\hat{\bm{m}}\cdot\hat{\bm{x}})\,\hat{\bm{x}}-\hat{\bm{m}}\Bigr{)}\qquad\text{for}\qquad r>R_{\text{WD}}\,,$$

$$\begin{array}{r}{\mathcal{L}}_{a,\Lambda }\supset -\frac{1}{4}{g}_{aWW}^{\Lambda }aW\tilde{W}-\frac{1}{4}{g}_{aBB}^{\Lambda }aB\tilde{B}+\frac{C_e^{\Lambda }}{2}\frac{{\partial }_{\mu }a}{f_a}\bar{e}{\gamma }^{\mu }{\gamma }_{5}e,\end{array}$$

$$\begin{array}{r}{C}_e^{\mu }=C_e^{\Lambda }+\frac{3}{8{\pi }^2}{\alpha }_{\mathrm{em}}^2\left(\frac{3}{8}\frac{C_W^{\Lambda }}{s_w^4}+\frac{5}{8}\frac{C_B^{\Lambda }}{c_w^4}\right)\log \frac{{\Lambda }^2}{{\mu }^2},\end{array}$$

$$\begin{array}{r}{C}_e^{\mu }\approx C_e^{\Lambda }+(5.2\times 10^{-4})C_W^{\Lambda }+(7.9\times 10^{-5})C_B^{\Lambda },\end{array}$$

$$\begin{array}{r}|C_e^{\mu }|\sim (5\times 10^{-4})|C_{\gamma }|\end{array}$$

$$\begin{array}{r}(2\times 10^{-4})|g_{a\gamma \gamma }|\mathrm{GeV}<|g_{aee}|<2.8\times 10^{-13}.\end{array}$$

Full text

††preprint: LCTP-19-05

The quantum chromodynamics (QCD) axion, originally proposed to solve the strong CP problem Peccei and Quinn (1977a, b); Weinberg (1978); Wilczek (1978), is a well-motivated extension of the Standard Model of particle physics. The QCD axion is a light pseudoscalar particle that couples to the QCD operator $G\tilde{G}$, with $G$ the QCD field strength. Additionally, the axion has dimension-5 couplings to electromagnetism and to matter. Studies of string compactifications show that, in addition to the QCD axion, there may exist a number of additional light pseudoscalar particles, with couplings to electromagnetism and matter but not to QCD Svrcek and Witten (2006); Arvanitaki et al. (2010). These pseudoscalars are called axion-like particles (ALPs), though throughout this work we will refer to all such particles as axions. In this work we present a novel method, using $X$-ray observations of magnetic white dwarf (WD) stars (MWDs), to probe the existence of axions.

WDs have long been used as probes of axions by studying the possibility of energy loss from axion emission Raffelt (1986). Axions are emitted by axion bremsstrahlung in electron-nucleon scattering. The extra energy loss would modify WD cooling and thus change the luminosity function of WDs. Comparisons to the observed luminosity function have been used to set stringent constraints on the axion-electron coupling Isern et al. (2008, 2009); Miller Bertolami et al. (2014). Moreover, it has been suggested that the observed WD luminosity function actually prefers the existence of an axion Isern et al. (2008, 2009), a claim further supported by period-drift measurements of WDs undergoing pulsations Isern et al. (2010).

In this work we propose to use $X$-ray observations of MWDs to detect the small fraction of emitted axions that convert to $X$-ray photons outside of the MWD in the strong surrounding magnetic field, as illustrated in Fig. 1. The conversion of axions to photons utilizes the axion-photon coupling. The proposed method uses the following key properties of isolated MWDs. First, the surrounding magnetic fields can be quite high, \sim$$10^{9} G Kepler et al. (2013); Ferrario et al. (2015), enhancing the axion-photon conversion probability. Second, while the core temperature of WDs is typically in the $X$-ray band, $T_{c}\sim 10^{7}\,\mathrm{K}$, the effective surface temperature is significantly lower, $T_{\rm eff}\sim 10^{4}\,\mathrm{K}$ (see, e.g., Ref. Chabrier et al. (2000)). Therefore an isolated MWD should not produce $X$-rays in the absence of axions. $X$-ray energy axions may escape the core and then convert into real $X$-ray photons in the magnetic field surrounding the WD, leading to a nearly thermal $X$-ray flux at the temperature of the core. In this work, we show that $X$-ray telescope observations of MWDs have the potential to probe a wide range of axion masses and couplings. However, as we will show, the $X$-ray observations are only sensitive to low-mass ALPs and not to the QCD axion, as the axion-photon conversion probability is suppressed for large axion masses.

A similar approach to axion detection was previously suggested for neutron stars (NSs) Morris (1986). In that case, nuclear bremsstrahlung in the NS generates an outgoing flux of axions, which may convert into $X$-rays in the strong magnetic field surrounding the NS. However, it was pointed out in Ref. Raffelt and Stodolsky (1988) that vacuum birefringence effects from quantum electrodynamics (QED), in the presence of strong magnetic fields, stymie the axion-photon conversion process. The result is that the axion-induced $X$-ray flux from NSs is expected to be negligible for a KSVZ Dine et al. (1981); Zhitnitsky (1980); Srednicki (1985) or DFSZ Kim (1979); Shifman et al. (1980) QCD axion, but the flux may still be significant for ALPs Fortin and Sinha (2018, 2019). We will discuss these QED effects for MWDs below.

In many ways our proposal is reminiscent of the CAST experiment Arik et al. (2014, 2015); Anastassopoulos et al. (2017), which looks for axions produced within the Sun. These keV-energy axions travel to Earth where they are converted into $X$-ray photons by the strong magnetic field of the CAST experiment and detected with $X$-ray optics. The key conceptual difference between our MWD proposal and CAST, in addition to using MWDs as the axion source instead of the Sun, is that the conversion to $X$-rays takes place not in the lab but rather in the magnetosphere surrounding the MWD itself. A similar approach has in fact been suggested for the Sun (see, e.g., Ref. Zioutas et al. (2009)), whereby one looks at spectral and morphological distortions to the solar $X$-ray spectrum from axion-photon conversion in the solar magnetic field, but this is complicated by the fact that the Sun is already a strong $X$-ray source.

We note that a host of additional astrophysical probes of axions and axion dark matter (DM) have been proposed. These include radio signatures of axion DM conversion in NSs Pshirkov and Popov (2009); Huang et al. (2018); Hook et al. (2018); Safdi et al. (2018) and radio signatures of axion decay Caputo et al. (2018), supernova cooling Raffelt (2008); Fischer et al. (2016), energy loss in horizontal branch stars in globular clusters Ayala et al. (2014), and photon-axion oscillations leading to increased transparency of TeV gamma rays Simet et al. (2008); Sánchez-Conde et al. (2009); De Angelis et al. (2011); Ajello et al. (2016). Axion DM is also the subject of significant laboratory efforts at present Battaglieri et al. (2017); Shokair et al. (2014); Du et al. (2018); Brubaker et al. (2017a); Al Kenany et al. (2017); Brubaker et al. (2017b); Caldwell et al. (2017); Kahn et al. (2016); Foster et al. (2017); Ouellet et al. (2018); Chaudhuri et al. (2015); Silva-Feaver et al. (2016); Budker et al. (2014); Bähre et al. (2013); Bogorad et al. (2019).

${\bm{X}}$-ray flux calculation — It is useful to recall how an axion interacts with matter (see Ref. Tanabashi et al. (2018) for a review). The QCD axion couples to gluons through the operator ${\mathcal{L}}\supset-a\alpha_{s}G\tilde{G}/(8\pi f_{a})$, where $a$ is the axion field, $\alpha_{s}$ is the strong fine structure constant, $f_{a}$ is the axion decay constant, and $G$ is the QCD field strength tensor. The axion may also couple to electromagnetism through the operator ${\mathcal{L}}\supset-g_{a\gamma\gamma}aF\tilde{F}/4$, where $g_{a\gamma\gamma}$ is the axion-photon coupling with units of $\mathrm{GeV}^{-1}$, and $F$ is the QED field strength tensor. Finally the axion can interact with electrons through the operator ${\mathcal{L}}\supset g_{aee}/(2m_{e})\bar{e}\gamma^{\mu}\gamma_{5}e\partial_{\mu}a$, where $g_{aee}$ is the dimensionless axion-electron coupling, $m_{e}$ is the electron’s mass, and $e$ is the electron field. It is customary to write $g_{a\gamma\gamma}=C_{\gamma}\alpha_{\mathrm{em}}/(2\pi f_{a})$ and $g_{aee}=C_{e}m_{e}/f_{a}$, where $\alpha_{\mathrm{em}}$ is the electromagnetic fine structure constant, and $C_{\gamma}$ and $C_{e}$ are dimensionless parameters.

The coupling of axions to matter allows axions to be emitted from WDs. WDs have dense cores that are supported by electron degeneracy pressure and consist predominantly of a hot plasma of carbon, oxygen, and electrons. The temperature of this isothermal plasma, which makes up the bulk of the matter of the WD, can be $T_{c}\approx 10^{7}$ K or higher, depending on the WD. Axions are emitted most efficiently in WDs from bremsstrahlung off of electrons in electron-nuclei scattering, through the diagram illustrated in Fig. 1. The luminosity radiated from a WD in axions, denoted by $L_{a}$, is calculated in Refs. Raffelt (1986); Nakagawa et al. (1988, 1987) (see also Ref. Raffelt (1990)), and for a carbon-oxygen WD with mass density ratio $0.5$ and total mass $M_{\text{WD}}$, we find

[TABLE]

where $L_{\odot}$ and $M_{\odot}$ are the luminosity and mass of the Sun, respectively. The energy spectrum of the axion emission is found to be thermal at temperature $T_{c}$: $dL_{a}/dE\propto E^{3}/(e^{E/T_{c}}-1)$ Nakagawa et al. (1988).

The emitted axions may be converted into $X$-ray photons in the magnetic field surrounding the MWD. Under the approximation where we assume all axions travel along radial trajectories originating from the MWD center, we may derive a simple analytic expression for the conversion probability. Note that in the Supplemental Material (SM) we derive and numerically solve the relevant equations of motion for the more general trajectories and field configurations appropriate for a realistic MWD with finite extent, but the results are similar to those found by assuming radial trajectories. The axion-photon mixing equations in the presence of an external magnetic field can be reduced to a system of first-order differential equations using a WKB approximation. This approximation assumes that the scale of variation in the magnetic field is much larger than the axion’s de Broglie wavelength. By working in the Weyl gauge, $A^{0}=0$, and focusing on axion trajectories along which the angle with respect to the magnetic field does not change, these equations take the form Raffelt and Stodolsky (1988)

[TABLE]

where $E$ is the axion’s energy. Here $A_{\parallel}(r)$ denotes the vector potential component in the plane normal to the direction of propagation and parallel to the external magnetic field, as a function of the radial coordinate $r$, while $a(r)$ is the axion field. The probability for an axion to convert into a photon $p_{a\to\gamma}$ is determined by solving Eq. (2) and comparing the magnitude squared of an initial pure axion state with the asymptotic solution for the electromagnetic vector potential. In the background magnetic field, the axion-photon interaction induces a mixing, which is parametrized by $\Delta_{B}(r)=(g_{a\gamma\gamma}/2)B(r)\sin\Theta$, where $B(r)$ is the strength of the magnetic field at radius $r$ and $\Theta$ is the angle between the radial propagation direction and the magnetic field, which is $r$-independent. The term $\Delta_{a}=-m_{a}^{2}/(2E)$ incorporates the axion mass and is responsible for the slightly different momenta between the axion and photon states. Strong-field QED effects in vacuum give rise to the term $\Delta_{\parallel}(r)=(7/2)E\xi(r)\sin^{2}\Theta$, with $\xi(r)=(\alpha_{\mathrm{em}}/45\pi)[B(r)/B_{\rm crit}]^{2}$ with $B_{\rm crit}=m_{e}^{2}/e\approx 4.41\times 10^{13}\ \mathrm{G}$ the critical field strength. In general $\Delta_{\parallel}$ also contains a term related to the photon’s effective plasma mass $\omega_{\rm pl}$, $\Delta_{\parallel}=-\omega_{\mathrm{pl}}^{2}/(2E)$, but this term is subdominant to the QED one for the systems that we consider.

Since the $B$ field strength decreases as we move farther away from the surface of the MWD, it is important to solve the equation of motion Eq. (2) including the changing magnetic field profile with $r$. In the weak-mixing limit we can solve Eq. (2) using the formalism of time-dependent perturbation theory, and the axion-photon conversion probability is found to be Raffelt and Stodolsky (1988); Fortin and Sinha (2019)

[TABLE]

The integral starts at the star’s surface, $r=R_{\text{WD}}$, since it is assumed that $X$-ray photons produced inside the MWD cannot escape. To first approximation we may model the magnetic field outside of the MWD as a magnetic dipole, such that $B(r)=(R_{\text{WD}}/r)^{3}B_{0}$ with $B_{0}$ the field at the surface and $r>R_{\text{WD}}$. In the SM we consider more general magnetic field configurations, but the results are largely the same. In general, we evaluate Eq. (3) numerically. Note that for typical MWD parameters and asymptotically small $m_{a}$, the conversion probabilities are of order $10^{-4}\times(g_{a\gamma\gamma}/10^{-11}\ \mathrm{GeV}^{-1})^{2}$.

Having calculated both the spectrum of axion radiation, $dL_{a}/dE$, and the axion-to-photon conversion probability, $p_{a\to\gamma}$, the flux of axion-induced $X$-ray photons at Earth is calculated as

[TABLE]

where $d_{\text{WD}}$ is the distance between Earth and the MWD. Note that for a typical MWD $d_{\rm WD}\approx 30$ pc away, combining Eq. (4) with the conversion probability estimate and the axion luminosity Eq. (1) leads to axion-induced $X$-ray fluxes of order $F_{\gamma a}\sim 10^{-15}\times(g_{aee}g_{a\gamma\gamma}/10^{-24}\,\,{\rm GeV}^{-1})^{2}$ erg$/$cm2$/$s at low axion masses, though in practice we compute $dF_{\gamma_{a}}/dE$ precisely using the above formalism given the properties of the MWD under consideration. In practice we compute the integrated flux $F_{2-10}$, defined as the integral of Eq. (4) from 2 to 10 keV, to compare to data.

The only parameter that appears in Eq.(4) that cannot be measured directly for a given MWD is the core temperature, $T_{c}$, which affects the axion luminosity through Eq.(1). The core temperature is not directly observable, since the thin WD atmosphere is largely opaque to radiation. The effective temperature $T_{\rm eff}$ of the atmosphere, which determines the observed luminosity $L_{\gamma}$, is much smaller than the temperature of the isothermal core $T_{c}$. Understanding the relation between $L_{\gamma}$ and $T_{c}$ requires detailed modeling of the layers connecting the atmosphere to the degenerate core (see e.g., Refs. Hansen (1999); Salaris et al. (2000); Chabrier et al. (2000); Renedo et al. (2010)). We use the result of models presented in Ref. Chabrier et al. (2000), which finds that for luminosities above $L_{\gamma}\sim 10^{-3.8}$ $L_{\odot}$ and below \sim$$10^{-1.5} $L_{\odot}$ the luminosity-core temperature relation is well approximated by

[TABLE]

As we discuss further in the SM, the WD models in Refs. Hansen (1999); Salaris et al. (2000); Chabrier et al. (2000); Renedo et al. (2010) agree to within $O(10\%)$ uncertainty.

It is important to note that while the effective temperatures $T_{\rm eff}$ are often at the level of $10^{4}$ K, the emission from the WD does not follow a pure thermal distribution at this temperature since higher frequencies probe deeper within the WD atmosphere Shipman (1976). Still, for the hot WDs we consider, the expected thermal hard $X$-ray flux (e.g., from 2 to 10 keV) is negligible compared to the axion-induced fluxes that will be probed.

${\bm{X}}$-ray data analysis and projections — Although there are over $70000$ known WD stars within approximately $100\,\mathrm{pc}$ of Earth Brown et al. (2018), the total number of known MWDs with well-measured properties is only around $200$ Ferrario et al. (2015). In the SM we list several MWDs that are expected to have the largest axion-induced $X$-ray flux, while in this section we focus on the most promising candidate.

The magnetic WD RE J0317-853 Barstow et al. (1995) is especially hot, has a strong magnetic field, and is relatively nearby, making it an excellent target for $X$-ray searches for axions. Optical and UV observations suggest an effective temperature between $T_{\rm eff}=30000$ K and $60000$ K, while at the same time being incredibly massive, $M_{\text{WD}}\sim 1.3$ $M_{\odot}$, and compact, $R_{\text{WD}}\sim 0.003-0.004$ $R_{\odot}$ Barstow et al. (1995); Külebi et al. (2010). Note that the best-fit temperature presented in Ref. Barstow et al. (1995) is $T_{\rm eff}=4.93_{-0.12}^{+0.22}\times 10^{4}$ K, though a wider range was considered in Ref. Külebi et al. (2010). The MWD rotates with a period of \sim$$725 s, and over this period the surface magnetic field strength is seen to oscillate. Reference Burleigh et al. (1999) used time-resolved far-UV spectroscopic data from the Hubble Space Telescope to study the variations in the locations and magnitudes of Lyman $\alpha$ transition to measure the magnetic field profile. They found that the observable surface magnetic field varies between 200 and 800 MG over the WD rotation period, which is well fit by an offset dipole model Burleigh et al. (1999). They also find a best-fit temperature $T_{\rm eff}\sim 40000$ K.

In the SM we model in detail the magnetic field structure of this MWD and assess the sensitivity of our results to uncertainties in the value of $T_{\mathrm{eff}}$, for example. Here, however, we simply calculate the predicted $X$-ray flux using Eq. (3) with the most conservative choice of parameters for the magnetic field strength and effective temperature of RE J0317-853. We use a central dipole field with surface-field value $B_{0}=200\ \mathrm{MG}$ and $\sin\Theta=1$, corresponding to the lowest surface-field value observed in Ref. Burleigh et al. (1999) anywhere on the WD surface over its period. We also take the fiducial values $T_{\rm eff}=30000$ K, $M_{\text{WD}}=1.32$ $M_{\odot}$, and $R_{\text{WD}}=0.00405$ $R_{\odot}$ Külebi et al. (2010), which let us infer $L_{\gamma}=0.0120$ $L_{\odot}$ from the Stefan-Boltzmann law and $T_{c}=2.0\times 10^{7}$ K from Eq. (5). We note that higher $T_{\rm eff}$ lead to higher flux rates, and our fiducial value for $T_{\rm eff}$ is the lowest end of the range discussed in Refs. Barstow et al. (1995); Külebi et al. (2010). Parallax measurements from Gaia-DR2 Brown et al. (2018) place RE J0317-853 at a distance of $d_{\text{WD}}=29.54\pm 0.04$ pc from Earth.

The axion-induced $X$-ray spectrum from RE J0317-853 should peak around $E\sim 3\,T_{c}\sim 5\,\text{keV}$. Observations with Suzaku, using approximately $60\,\text{ks}$ of exposure time, detected no astrophysical $X$-ray emission from this MWD and set a flux limit in the $2$-$10\,\text{keV}$ range of $F_{2-10}<1.7\times 10^{-13}\,\text{erg}/\text{cm}^{2}/\text{s}$ at 95% confidence Harayama et al. (2013). We note that the limit in Ref. Harayama et al. (2013) required background subtraction and modeling; the limit itself is dominated by systematic uncertainties in modeling the cosmic $X$-ray background and the non-$X$-ray background. We also caution that the limit in Ref. Harayama et al. (2013) is formally only valid for an energy spectrum that resembles the above backgrounds, given the energy-dependence of the Suzaku effective area, though we have checked that this only induces a $\sim$10% difference in the predicted counts and thus can be ignored for our purposes. We translate the flux limit into a 95% constraint on the axion coupling constants $|g_{a\gamma\gamma}\,g_{aee}|$ using the fiducial values for RE J0317-853 and the formalism for axion-photon conversion developed above. Our results are presented in Fig. 2, which shows our constraint on the axion parameter space; the region above the blue curve is excluded at $95\%$ confidence for the fiducial stellar parameters ($T_{\mathrm{eff}}$, $B_{0}$, etc) given above. Note that the nontrivial structure in the limit at high $m_{a}$ is due to transitioning across the regime where $\Delta_{a}\ll\Delta_{\parallel}$ and then again to the regime where $\Delta_{a}\gg R_{\text{WD}}$.

We compare our result to the previous best limits on this coupling combination. By searching for axions produced inside of the Sun though the axion-electron coupling, the CAST experiment obtains a direct limit on $|g_{a\gamma\gamma}\,g_{aee}|$ Barth et al. (2013). Our constraints are approximately 2 orders of magnitude stronger than those from CAST at low axion masses. Several experiments derive upper limits on $|g_{a\gamma\gamma}|$ alone. CAST provides an upper limit on $|g_{a\gamma\gamma}|$ from axions produced in the Sun through the Primakoff process Anastassopoulos et al. (2017). Similarly an upper bound on $|g_{a\gamma\gamma}|$ is derived from the nonobservation of a gamma ray flux coincident with SN 1987A Payez et al. (2015), which is expected to be seen if axions are produced during the core-collapse supernova and subsequently convert into $\gamma$-ray photons in the intergalactic magnetic field. By combining these limits on $g_{a\gamma\gamma}$ with the limit $|g_{aee}|<2.8\times 10^{-13}$, which arises from modeling the WD luminosity function Miller Bertolami et al. (2014), we obtain a second $95\%$ confidence upper limit, which is stronger than CAST’s direct limit on $|g_{a\gamma\gamma}\,g_{aee}|$.

In Fig. 2 we also show the best-fit value for $g_{a\gamma\gamma}\times g_{aee}$ (“cooling hints”) found in Ref. Giannotti et al. (2017) from a global fit to the available stellar cooling data, which slightly favors nonzero axion-electron and axion-photon couplings. Unfortunately the global fits say little about the axion mass, since the stellar cooling probes are not sensitive to the mass of the particle causing the cooling. From our work we are able to say that the axion mass at these best-fit couplings needs to be above \sim$$2\times 10^{-5} eV. The coupling-mass relations for the DFSZ Kim (1979); Shifman et al. (1980) and KSVZ Dine et al. (1981); Zhitnitsky (1980); Srednicki (1985) QCD axion models are also shown in Fig. 2.

The Suzaku observations of MWD RE J0317-853 from Harayama et al. (2013) are not ideal for searching for the proposed axion-induced $X$-ray signal. This is because Suzaku has a rather poor point-source sensitivity in the $\sim$2-10 keV band compared to other telescopes like XMM-Newton or Chandra. For instance, we estimate that a $\sim$400 ks observation with Chandra would yield a 95% confidence flux sensitivity in the $2$-$10$ keV band at the level of $6\times 10^{-16}\,\text{erg}/\text{cm}^{2}/\text{s}$  Chandra X-Ray Observatory (2019), which is over 2 orders of magnitude better than the sensitivity achieved in Harayama et al. (2013) with Suzaku. (Specifically, this flux sensitivity estimate includes $X$-ray and particle backgrounds as estimated in Chandra X-Ray Observatory (2019) and is for the ACIS-I instrument in the timed exposure mode with no grating and with CCD I3, along with Very Faint telemetry.) This projected sensitivity appears in Fig. 2 as the red dashed curve. Note that this sensitivity projection assumes that the MWD does not produce $X$-ray emission in the energy range of interest at the flux levels that would be probed by XMM-Newton and Chandra. This assumption is justified because the thermal emission is exponentially suppressed in the energy range 2-10 keV and well below the levels that would be probed by future observations. Possible nonthermal emission mechanisms include synchrotron and curvature radiation in the strong magnetic field surrounding the MWD. However, these processes are suppressed by the rather large spin period of the MWD, which leads to a small accelerating potential compared to e.g. the fast-spinning cataclysmic variable (CV) MWD AE Aquarii, for which pulsed nonthermal emission has been possibly observed Terada et al. (2008); Harayama et al. (2013).

Discussion — In this Letter we used Suzaku observations of the nearby MWD RE J0317-853 to set the strongest limits to date on $|g_{a\gamma\gamma}\,g_{aee}|$ for axion masses $\lesssim$ meV. Dedicated observations, with existing telescopes such as XMM-Newton, Chandra, or NuSTAR of this MWD and others have the potential to improve the sensitivity by over an order of magnitude and perhaps more. If the previously observed stellar cooling hints are due to low-mass axion-like particles, $X$-ray observations of MWDs should detect excess hard $X$-ray flux.

If a hard $X$-ray signal is seen from a MWD, the first question should be if the MWD is accreting, which would be the case for a CV star. CVs, which often emit hard $X$-ray spectra are distinguishable from e.g. their variability and emission line spectra. It is also possible that some MWDs emit hard $X$-rays due to chromospheric activity, though this has yet to be observed Brinkworth et al. (2004). Our predicted axion signal has the unique feature of having an order unity modulation fraction over the MWD period for most alignment angle configurations, as discussed more in the SM, which could help differentiate it from other backgrounds.

The signal discussed in this work requires the axion to interact with both electromagnetism and electrons. While both interactions are expected in a generic ALP theory, from an effective field theory point of view, this fact makes it hard to directly compare our sensitivity to e.g. the CAST limits on $g_{a\gamma\gamma}$ alone. However, in the SM we show that even if one assumes that the axion-electron operator is generated through the renormalization group alone due to $W$ and $Z$-boson loops, the constraints from this work may be translated to constraints on $g_{a\gamma\gamma}$ that are comparable to those from CAST.

In addition to the axion-electron coupling, a generic axion is also expected to interact derivatively with quarks. These interactions cause NSs to radiate $\sim$keV axions, as in the WD case described in this work. NSs are also promising targets for axion-induced hard $X$-ray signals, and this will be explored in future work Buschmann et al. (2021).

I Additional MWD candidates for $X$-ray observation

In the main text we have focused our analysis on the MWD star RE J0317-853, since it is expected to have a particularly strong $X$-ray flux, and because $X$-ray data is already available from Suzaku. However, there are over $200$ MWD stars with well-measured field strengths, temperatures, and distances. For each of these stars, we calculate the expected axion-induced $X$-ray flux in the $2-10\ \mathrm{keV}$ energy window, denoted by $F_{2-10}$, assuming $m_{a}=10^{-9}\ \mathrm{eV}$ and $g_{a\gamma\gamma}g_{aee}=10^{-24}\ \mathrm{GeV}^{-1}$. The flux is insensitive to the axion mass in the limit $m_{a}\to 0$, and more generally the flux has an overall scaling with the couplings, $F_{2-10}\propto(g_{a\gamma\gamma}g_{aee})^{2}$.

Our results are summarized in Table 1, which shows the ten MWDs with the largest predicted $X$-ray flux, $F_{2-10}$. We constructed a full list by merging the SDSS DR7 magnetic WD catalog Kleinman et al. (2013), a review MWD catalog Ferrario et al. (2015), and Gaia DR2 WD catalog Brown et al. (2018). The former two provide the magnetic field strengths and temperatures of the WDs, while the latter provides distances, luminosities, masses, and radii, if known. Since the mass and radius of WD 2010+310 are not known, we take its mass to be $M_{\text{WD}}=1\,M_{\odot}$, and we infer its radius from the Stefan-Boltzmann law.

II The probability for axion-photon conversion in a general magnetic field background

In this section we present a more general formalism to calculate the axion-photon conversion probability, and in the following section we apply this calculation to study the $X$-ray emission from RE J0317-853 in more detail. Interactions between the axion and photon are described by the following Lagrangian:

[TABLE]

where $m_{a}$ is the axion’s mass, $m_{A}$ is the photon’s effective mass, $j^{\mu}$ is the electromagnetic current density, $g_{a\gamma\gamma}$ is the axion-photon coupling, $\alpha_{\mathrm{em}}$ is the electromagnetic fine structure constant, and $m_{e}$ is the electron’s mass. An environment with a nonzero electron density, $n_{e}$, gives rise to an effective mass, $m_{A}=\omega_{\mathrm{pl}}=\sqrt{4\pi\alpha_{\mathrm{em}}n_{e}/m_{e}}$, and although this term will be negligible for our calculation, we retain it here for the sake of generality. The last term in (S1), known as the Euler-Heisenberg term, describes the photon’s self-interaction that arises at energies below the electron’s mass, where the electron can be integrated out of the theory.

One can derive the field equations for $a({\bm{x}},t)$ and $A_{\mu}({\bm{x}},t)$ by applying the variational principle to (S1). It is customary to work in the Weyl gauge, $A_{0}=A^{0}=0$ Raffelt and Stodolsky (1988). Upon writing $A_{i}=-A^{i}=-{\bm{A}}_{i}$ and $j^{\mu}=(\rho,{\bm{j}})$, the field equations are expressed as

[TABLE]

Here $\partial_{0}=\partial_{t}$ is denoted by a dot, and $\partial_{i}={\bm{\nabla}}_{i}$.

To study axion-photon conversion, we are interested in the field equations for linearized perturbations around static background fields. This motivates us to write

[TABLE]

where $\bar{a}({\bm{x}})$ and $\bar{\bm{B}}({\bm{x}})={\bm{\nabla}}\times\bar{\bm{A}}$ are the background axion and magnetic fields. In particular, we are interested in systems for which $\bar{a}=0$ and $\bar{\bm{B}}$ can be approximated as a magnetic dipole.

We assume that gradients in the background field, $\bar{\bm{B}}$, are small compared to the wavenumber of the axions, $k=|{\bm{k}}|$, which lets us apply the WKB approximation. Consider the trajectory ${\bm{x}}(z)={\bm{x}}_{0}+z\,\hat{\bm{k}}$ where ${\bm{x}}_{0}$ is a point inside the WD, $z\geq 0$ parametrizes the distance, and the unit vector $\hat{\bm{k}}={\bm{k}}/k$ is the trajectory’s orientation. If gradients in the directions transverse to $\hat{{\bm{k}}}$ are small, then we are motivated to adopt the plane-wave Ansatz

[TABLE]

For relativistic particles, it is an excellent approximation to take the natural frequency as $\omega\approx k$. We choose a set of orthonormal basis vectors, $\{\hat{\bm{e}}_{x},\,\hat{\bm{e}}_{y},\,\hat{\bm{e}}_{z}\}$, such that $\hat{\bm{e}}_{z}=\hat{\bm{k}}$ and therefore $\hat{\bm{e}}_{x}\cdot\hat{\bm{k}}=\hat{\bm{e}}_{y}\cdot\hat{\bm{k}}=0$. We can also write $B_{x}(z)=\bar{\bm{B}}\cdot{\bm{e}}_{x}$, $B_{y}(z)=\bar{\bm{B}}\cdot{\bm{e}}_{y}$, $B_{T}(z)=[B_{x}^{2}+B_{y}^{2}]^{1/2}$, $B_{L}(z)=\bar{\bm{B}}\cdot\hat{\bm{k}}$, and $B(z)=|\bar{\bm{B}}|=[B_{T}^{2}+B_{L}^{2}]^{1/2}$. In the WKB approximation, we assume that the background changes slowly along the axion’s trajectory. This lets us drop derivatives of the background field and also second derivatives of the perturbations. With these simplifications, the Weyl-gauge field equations (S2) become

[TABLE]

where we have defined

[TABLE]

In deriving (S5), we have dropped terms that are nonlinear in the perturbations. The wavefunction $\tilde{A}_{z}$ is not dynamical, meaning that its field equation is algebraic rather than differential, and we can remove it from the system of equations. We regain (2) in the main text by taking $z=r$, $\omega=k$, $m_{A}=0$, and by focusing on isotropic field configurations for which $B_{x}=0$, $B_{y}(r)=B(r)\sin\Theta$, and $B_{z}(r)=B(r)\cos\Theta$.

The axion-to-photon conversion probability is calculated by solving (S5). We assume that photons produced inside of the star, $r<R_{\text{WD}}$, will be unable to escape as $X$-ray emission, and therefore we are only interested in conversions that occur outside of the star. Along a given trajectory, let $z=z_{0}$ denote the surface of the star. We solve (S5) along with the initial condition $a(z_{0})=a_{0}$, $A_{x}(z_{0})=0$, and $A_{y}(z_{0})=0$, and then the axion-to-photon conversion probability is given by p_{a\to\gamma}=\bigl{(}|A_{x}(\infty)|^{2}+|A_{y}(\infty)|^{2}\bigr{)}/a_{0}^{2}. Since (S5) is linear, the initial condition $a_{0}$ cancels when calculating $p_{a\to\gamma}$.

III Analytic approximations to the conversion probability

In order to develop our intuition, let us consider axion-photon conversion in a constant magnetic field, $B(r)=B$, which extends for a distance $\Delta r=R_{\text{WD}}$ outside of the MWD, with $R_{\text{WD}}$ the star’s radius. Then (2) may be solved exactly, and the axion-photon conversion probability is found to be

[TABLE]

where the mixing angle, $\theta$, and the oscillation length scale, $\Delta_{\mathrm{osc}}^{-1}$, can be written as

[TABLE]

Note that $\Delta_{\parallel}-\Delta_{a}$ is always positive.

We will generally be interested in the weak-mixing regime where $\theta\ll 1$, and for the sake of illustration let us also focus on the regime where the axion mass is low, such that $|\Delta_{a}|\ll\Delta_{\parallel}$. With these assumptions, (S7) becomes $p_{a\to\gamma}\approx 4(\Delta_{B}^{2}/\Delta_{\parallel}^{2})\,\sin^{2}(\Delta_{\parallel}R_{\text{WD}}/2)$. If the magnetic field is weak, such that $\Delta_{\parallel}R_{\text{WD}}\ll 1$, then the conversion probability becomes $p_{a\to\gamma}\approx\Delta_{B}^{2}R_{\text{WD}}^{2}\propto g_{a\gamma\gamma}^{2}B^{2}$. Alternatively, if the magnetic field is very strong, such that $1\ll\Delta_{\parallel}R_{\text{WD}}$, then we have instead $p_{a\to\gamma}\propto g_{a\gamma\gamma}^{2}/(BE)^{2}$. Here we see that the conversion probability is suppressed in the strong-magnetic field regime, which is a consequence of the QED birefringence effects, as anticipated in Raffelt and Stodolsky (1988). The axion-photon conversion probability is maximal for $\Delta_{\mathrm{osc}}R_{\text{WD}}\sim 1$. In practice, for frequencies $E\sim 1$ (10) keV, this transition occurs for magnetic fields around \sim$$10^{7} (\sim$$5\times 10^{6}) G. The dipole magnetic fields around the MWDs fall with distance from the WD surface, which leads to the result that if the magnetic field at the surface of the WD is high, the conversion probability is suppressed until the field drops to where $\Delta_{\mathrm{osc}}R_{\text{WD}}\sim 1$.

As we increase the axion mass, we enter the regime where $\Delta_{\parallel}\ll|\Delta_{a}|$ and $p_{a\to\gamma}\approx 4(\Delta_{B}^{2}/\Delta_{a}^{2})\,\sin^{2}(\Delta_{a}R_{\text{WD}}/2)$. For asymptotically large axion masses we have also $|\Delta_{a}|R_{\text{WD}}\gg 1$, in which case $p_{a\to\gamma}\propto g_{a\gamma\gamma}^{2}B^{2}E^{2}/m_{a}^{4}$. Thus, contrary to the search for low-mass axions, searches for high-mass axions benefit from the largest magnetic fields and largest frequencies available.

With the dipole field configuration we instead need to compute the integral given in (3). While in general this must be done numerically, we note that an analytic approximation is available in the regime of small axion mass, such that $|\Delta_{a}|\ll\Delta_{\parallel}$, where we find

[TABLE]

with $\Gamma(z)$ the gamma function, $\Gamma(a,z)$ is the incomplete gamma function, $\Delta_{B,0}=\Delta_{B}(R_{\text{WD}})$, and $\Delta_{\parallel,0}=\Delta_{\parallel}(R_{\text{WD}})$.

To help illustrate the different regimes for the conversion probability, in Fig. S1 we illustrate the conversion probability computed using (S9) as a function of the transverse magnetic field $B_{T,0}$ at $r=r_{\text{WD}}$ for an axion with an arbitrary small mass, a reference $g_{a\gamma\gamma}=10^{-11}\ \mathrm{GeV}^{-1}$, and a MWD with $R_{\text{WD}}=0.00405$ $R_{\odot}$. For definiteness we show the conversion probability for $E=10\ \mathrm{keV}$. At small values of $B_{T,0}$ the conversion probability rises like $B_{T,0}^{2}$. However, the Euler-Heisenberg term becomes important for $B_{T,0}\gtrsim 10$ MG, causing the conversion probability to scale less strongly with $B_{T,0}$ beyond this point.

IV A detailed analysis of the axion-induced $X$-ray flux from RE J0317-853

Our analysis of the $X$-ray flux from RE J0317-853 in the main text makes several simplifying assumptions. In this section we revisit and refine these assumptions.

IV.1 Magnetic field structure

In the main text we calculated the probability for axions to convert into $X$-ray photons in a background magnetic field by assuming that the background field is isotropic along the axion’s trajectory and decreases with propagation distance like $B\sim r^{-3}$. This is a reasonable assumption for a dipole field configuration if axions propagate on radial trajectories, as if they were produced at the star’s center. However, in general the axion’s trajectory will not originate at the origin and, moreover, the background field will deviate from a simple dipole configuration. Here we assess the effect of relaxing these assumptions.

Let us first understand how our results depend on the axion’s trajectory through the magnetosphere. Consider a dipolar magnetic field centered on the center of the star. Such a field can be written as

[TABLE]

where ${\bm{m}}=|{\bm{m}}|\,\hat{\bm{m}}$ is the magnetic dipole moment, $|{\bm{m}}|/4\pi=B_{0}R_{\text{WD}}^{3}/2$ is the relationship with the polar field $B_{0}$, ${\bm{x}}=r\,\hat{\bm{x}}$ is the spatial coordinate, and $r=|{\bm{x}}|$ is the distance from the center of the MWD. Note that here $B_{0}$ is the value of the magnetic field at the star’s surface in the direction of the magnetic pole, whereas in the main text below (3) we used $B_{0}$ to denote the value at the surface in any arbitrary direction. Note that the orientation of $\bar{\bm{B}}$ along a radial trajectory is invariant, which allows us to calculate the axion-photon conversion probability, $p_{a\to\gamma}$, using (3), in this case. The result is presented in the left panel of Fig. S2 as the blue curve, where we show the conversion probability as a function of the angle between the magnetic pole and the viewing angle, $\hat{\bm{m}}\cdot\hat{\bm{x}}=\cos\theta$. As one expects, the probability vanishes at $\theta=0$ and $\pi$ where the transverse magnetic field vanishes, and the probability peaks at $\theta=\pi/2$ where the transverse magnetic field, $B_{T}=B\,\sin\theta$, is maximal. It is a rough approximation to suppose that axions only propagate on radial trajectories, and more accurately we should allow the axion’s trajectory to originate anywhere inside the MWD. In general the field direction is not invariant along such a trajectory, and we cannot use (3) to calculate $p_{a\to\gamma}$, but rather we must solve (S5) directly. Doing so leads to the red-dashed curve in the left panel of Fig. S2, where we have sampled $3,\!000$ points in the star’s interior, uniformly distributed in ${\bm{x}}$ with $|{\bm{x}}|<R_{\text{WD}}$. The effects of this averaging are to remove the oscillatory behavior and enhance the probability at $\theta=0,\pi$. However, the rough approximation is reliable up to an order unity correction factor.

Let us finally address the specific magnetic field structure in our MWD candidate, RE J0317-853. Using phase-resolved far-ultraviolet spectroscopy data, Burleigh et al. (1999) constructed several models for the field configuration of RE J0317-853. It was found that the observations are well described by a magnetic dipole, with a polar field strength of $B_{0}=363\,\mathrm{MG}$, which is offset from the star’s center along the magnetic axis by $19\%$ of $R_{\text{WD}}$, which is misaligned with the star’s rotation axis by an angle of $20^{o}$, and which is viewed from an angle of $51^{o}$ to the rotation axis. We have calculated the axion-photon conversion probability for this field configuration by solving (S5) and averaging over $5000$ trajectories. Our results appear in the right panel of Fig. S2. Since the polar viewing angle is known, we vary the azimuthal angle, $\phi$, instead. Here we see that the axion-to-photon conversion probability experiences an $O(50\%)$ variation as the star revolves, which corresponds to a period of \sim$$725 s for RE J0317-853, and this translates into a corresponding modulation of the $X$-ray flux. Note that across all masses we find that when integrating out in radius $r$ to compute the conversion probability the probability has researched its asymptotic values by $r\sim 10\times R_{\rm WD}$.

Finally let us compare our limit on $|g_{a\gamma\gamma}\,g_{aee}|$, calculated using two different models for the background magnetic field. In the main text we used a configuration with $B(r)=(r/R_{\text{WD}})^{-3}\,B_{0}$ for $B_{0}=200\,\mathrm{MG}$, and we obtained the limit appearing in Fig. 2. If we perform the same calculation using the offset dipole configuration described above, then we obtain instead the limit appearing in Fig. S3 as the brown curve. Note that the limit derived in the main text is reliable up to an $O(1)$ factor, but it is also a conservative estimate. At small values of the axion mass, we find that the two limits are comparable. This result is consistent with Fig. S2 where we see $p_{a\to\gamma}\simeq 3.0\times 10^{-4}$ at $\theta=\pi/2$ in the left panel while $\langle p_{a\to\gamma}\rangle\approx 2.9\times 10^{-4}$ in the right panel after averaging over $\phi$. We note, though, that improved sensitivity may be obtained by using an analysis procedure that incorporates the expected non-trivial light curve, in the case of the displaced dipole model.

IV.2 Effective temperature uncertainties

The MWD RE J0317-853 has garnered significant interest, due in part to its high effective temperature, which takes a value in the range from $30,\!000$ to $55,\!000$ K Barstow et al. (1995). This uncertainty in $T_{\mathrm{eff}}$ translates into an uncertainty in our predicted $X$-ray flux, $F_{2-10}$, and therefore an uncertainty in our constraint on the axion’s parameter space. In order to quantify this uncertainty we compare our fiducial model from the main text with four benchmark models, given by Table 6 of Külebi et al. (2010). The results are summarized in Tab. 2, which shows that the predicted $X$-ray flux can vary by a factor of $\sim 4$ depending on the specific stellar model parameters used to describe RE J0317-853. Since the flux depends on the axion’s couplings through $F_{2-10}\propto(g_{a\gamma\gamma}g_{aee})^{2}$, the effect on our limit is a factor of $\sim 2$. The flux uncertainty follows primarily from the uncertainty in $T_{\mathrm{eff}}$. In Fig. S4 we show how these uncertainties translate into the uncertainty in the 95% limit from Suzaku observations and (projected) Chandra observations. The bands encompass the range of limits obtained by cycling through the parameters presented in Tab. 2.

IV.3 Core temperature model uncertainties

Since the WD’s core temperature is not directly observable, it is necessary to have a model for the WD interior in order to relate $T_{c}$ to an observable quantity like the luminosity, $L_{\gamma}$. Moreover, since the axion emission depends sensitively on the WD core temperature, with the axion luminosity going as $L_{a}\sim T_{c}^{4}$ in (1), it is therefore important to quantitatively assess how our prediction for the axion-induced X-ray signal depends on this modeling. In this section we specifically investigate how the modeling uncertainties affect our limits on the axion parameter space for RE J0317-853.

The physics of WD cooling is a balance between the energy stored in the hot core and the transparency of the envelope. The transport of thermal energy also depends on various physical inputs including the thermal conductivity of the degenerate matter, neutrino emission rates, and chemical diffusion. The cooling of old hydrogen-rich DA WDs has been studied extensively, and we summarize the relevant results of several prominent studies Hansen (1999); Salaris et al. (2000); Chabrier et al. (2000); Renedo et al. (2010) in Fig. S5. In particular we are interested in the predicted relation between the core temperature and surface luminosity. The left panel of Fig. S5 shows that all four models agree very well with our fiducial formula (5) for $10^{-4}\lesssim L_{\gamma}/L_{\odot}\lesssim 10^{-1}$.

At the fiducial luminosity of RE J0317-853, $L_{\gamma}=0.012\ L_{\odot}$, the predictions for the core temperature vary from $T_{c}\simeq 1.8$ to $2.0\times 10^{7}\ \mathrm{K}$, and thus we can associate the WD cooling model with an $O(10\%)$ uncertainty on $T_{c}$. Since the axion luminosity goes as $L_{a}\sim T_{c}^{4}$, we infer that the WD cooling model leads to an $O(40\%)$ uncertainty on our axion-induced X-ray flux signal and an $O(20\%)$ uncertainty on our axion parameter space limits for RE J0317-853 since $|g_{a\gamma\gamma}\times g_{aee}|\sim T_{c}^{2}$. This is illustrated in the right panel of Fig. S5.

IV.4 Spectrum

Fig. S6 shows the predicted spectra of axion and photon emission from our candidate MWD star, RE J0317-853. The shape of the axion spectrum (blue curve) is very well approximated by a blackbody at temperature $T_{c}\simeq 2\times 10^{7}\,\mathrm{K}\simeq 1.7\,\mathrm{keV}$ whereas the amplitude of the spectrum is set by the magnitude of the axion-electron coupling according to (1). To draw the blue curve we take $g_{aee}=10^{-13}$. The spectrum of secondary, axion-induced photons (red curve) tracks the thermal spectral shape up to an additional energy dependence coming from the axion-photon conversion probability (3). To draw the red curve we take $m_{a}=10^{-9}\,\mathrm{eV}$ and $g_{a\gamma\gamma}=10^{-11}\,\mathrm{GeV}^{-1}$, such that the product $|g_{a\gamma\gamma}\,g_{aee}|=10^{-24}$ is marginally consistent with the Suzaku limit. We assume that all axions are emitted isotropically and homogeneously from throughout the interior of the WD core, and the axion-photon conversion probability is calculated using (S5).

V The running of the axion-electron coupling and sensitivity to $g_{a\gamma\gamma}$

In the main text we have shown the limit on $|g_{a\gamma\gamma}\,g_{aee}|$ that follows from $X$-ray observations of MWD stars. Here we would like to translate this into a limit on simply $g_{a\gamma\gamma}$. For any value of $g_{a\gamma\gamma}$ there is a “reasonable range” of values for $g_{aee}$: it is bounded from above by direct observation, and it is bounded from below because an axion-electron coupling can be induced radiatively from an axion-photon coupling. In this section we estimate the lower bound.

Consider an axion-like particle, which does not couple to QCD. The effective theory describing the interactions of this axion with the Standard Model electroweak gauge fields is generated at a high scale, which we denote by $\Lambda\sim f_{a}\sim 10^{10}$ GeV, and the corresponding Lagrangian can be written as

[TABLE]

where $W\equiv W^{a}_{\mu\nu}$ is the $SU(2)_{L}$ field strength tensor, $B\equiv B_{\mu\nu}$ is the $U(1)_{Y}$ field strength tensor, $g_{aWW}^{\Lambda}=C_{W}^{\Lambda}\alpha_{W}/(2\pi f_{a})$, and $g_{aBB}^{\Lambda}=C_{B}^{\Lambda}\alpha_{B}/(2\pi f_{a})$. The fine structure constants can be written as $\alpha_{W}=\alpha_{\mathrm{em}}/s_{w}^{2}$ and $\alpha_{B}=\alpha_{\mathrm{em}}/c_{w}^{2}$ where the electromagnetic fine structure constant is $\alpha_{\mathrm{em}}\simeq 1/137$ and the weak mixing angle’s sine and cosine are $s_{w}=\sin\theta_{w}\simeq 0.48$ and $c_{w}=\cos\theta_{w}\simeq 0.88$. The parameters $C_{W}^{\Lambda}$, $C_{B}^{\Lambda}$, and $C_{e}^{\Lambda}$ are dimensionless constants that parameterize the UV theory. We are interested in how these parameters evolve, under the renormalization group, to scales $\mu\ll\Lambda$, since the physical processes we consider are at much lower energies. We assume that the scale $\mu$ is larger than the electroweak scale $\Lambda_{\rm EW}\sim 100$ GeV so that the Lagrangian (S11) is the correct description of the axion-gauge-boson interactions. Below this scale we should instead map to the axion-photon coupling, and we will return to this point shortly.

The one-loop diagrams that contribute to the running of $C_{e}$ were computed in Srednicki (1985); Bauer et al. (2017) and consist of triangle diagrams with electron and axion final states connected by a loop of $SU(2)_{L}$ or $U(1)_{Y}$ gauge bosons. Evaluating the divergent part of these diagrams and applying the renormalization procedure gives the beta function for the axion-electron coupling. The running coupling at the scale $\mu$, denoted by $C_{e}^{\mu}$, is given by Bauer et al. (2017)

[TABLE]

which illustrates how the axion-electron coupling is induced radiatively from the axion-gauge couplings in (S11). Evaluating (S12) at the electroweak scale, $\mu=\Lambda_{\rm EW}\simeq 100\ \mathrm{GeV}$ gives

[TABLE]

for typical values of $\Lambda\sim 10^{10}\ \mathrm{GeV}$.

Below the scale of electroweak symmetry breaking we are interested in the axion-photon coupling, $\mathcal{L}\supset-g_{a\gamma\gamma}aF\tilde{F}/4$ with $g_{a\gamma\gamma}=C_{\gamma}\alpha_{\mathrm{em}}/(2\pi f_{a})$. To leading order we have $C_{\gamma}=C_{W}^{\Lambda}+C_{B}^{\Lambda}$, which illustrates how the axion-photon interaction arises at low energies from the axion-$W_{\mu}^{a}$ and/or axion-$B_{\mu}$ interactions in the UV. Thus, barring any accidental cancellations between the $C_{W}$ and $C_{B}$ terms, we anticipate a relation between the axion-electron and axion-photon couplings, which is

[TABLE]

in the IR for theories where the axion-electron coupling is not present in the UV. Running the axion-electron coupling from the weak scale down to the electron mass scale, further enhances $C_{e}$ by a factor of $\sim$1-2. Additionally, the case of the QCD axion is slightly more complicated since $C_{\gamma}$ receives an IR contribution from pion mixing in that case. Therefore we simply use (S14) for the following estimates.

In Fig. 2 we expressed our sensitivity in terms of $|g_{a\gamma\gamma}\,g_{aee}|$, and now by using the expression for $C_{e}^{\mu}$ from (S14), we can map the sensitivity onto $g_{a\gamma\gamma}$ directly. A reasonable range of values for the axion-electron coupling is given by

[TABLE]

The lower limit follows from (S14) since $g_{aee}=C_{e}m_{e}/f_{a}$ and $g_{a\gamma\gamma}=\alpha_{\mathrm{em}}C_{\gamma}/2\pi f_{a}$. The upper limit is an empirical $95\%$ confidence constraint from modeling the WD luminosity function Miller Bertolami et al. (2014). Using this range of values for $g_{aee}$, our sensitivity curves from Fig. 2 translate into bands, which are shown in Fig. S7. We also show the CAST experiment’s limit on the axion-photon coupling Anastassopoulos et al. (2017), and the predictions of the KSVZ and DFSZ models of the QCD axion. For $m_{a}\lesssim 10^{-5}\ \mathrm{eV}$ our limit is comparable to the CAST 95% CL exclusion, even under the “pessimistic” assumption that $g_{aee}$ is “minimal” as in (S14), which leads to the weakest limit on $g_{a\gamma\gamma}$. On the other hand, dedicated Chandra observations would lead to a conservative upper limit on $g_{a\gamma\gamma}$ that is significantly stronger than the CAST bound at low masses. We also show an upper bound on $|g_{a\gamma\gamma}|$ from gamma-ray flux associated with SN1987a Payez et al. (2015).

Our results also constrain the axion explanations of the very-high-energy gamma-ray transparency anomalies previously observed with Cherenkov telescope data Horns and Meyer (2012); Meyer et al. (2013). The gamma-ray observations indicate that the Universe is more transparent than previously thought to high-energy gamma-rays, and one explanation is that gamma-rays could oscillate into axions in astrophysical magnetic fields. The allowed parameter space from Meyer et al. (2013) to fit these anomalies is indicated in shaded gray in Fig. S7, though we note that stringent constraints from e.g. the Fermi Large Area Telescope Ajello et al. (2016) and H.E.S.S. Abramowski et al. (2013) also exist below \sim$${\rm few}\times 10^{-8} eV. Furthermore, the transparency anomalies are subject to uncertainties from assumptions about galactic magnetic field models. Future dedicated observations of MWDs by e.g. XMM-Newton or Chandra would be able to probe much of the motivated parameter space to explain these anomalies.

For low-mass axions with $m_{a}\lesssim 10^{-6}\ \mathrm{eV}$, the upper limits on axion-matter couplings become insensitive to the axion’s mass, $m_{a}$. We summarize these limits in Fig. S8, where we also compare our limit with previous limits on the axion’s coupling to photons, $g_{a\gamma\gamma}$ and the axion’s coupling to electrons, $g_{aee}$.