Astrophysical axions, and using MESA simulate the lives of stars

The story of axion pumps

Back in 1991, Wick and his student Kar Lee were studying the potential for axions (a pseudoscalar particle hypothesized to explain, among other things, the lack of CP violation in the strong interaction) to be produced by nuclei at the hearts of stars [Ha91]. As pseudoscalars, axions should compete with magnetic photons in the nuclear magnetic dipole (M1) transition. However, nuclear M1 transitions, which involve a $\pm1$ change in either the total angular momentum magnitude $J$ or projection $M_J$ quantum numbers, are generally difficult to excite inside stars: Due to fusion’s affinity for $\alpha$-particle chains, the most common stellar isotopes ($^4_2$He, $^{12}_6$C, $^{16}_8$O, $^{20}_{10}$Ne, …) have even numbers of protons and neutrons. Because the strong interaction likes to group like-nucleons in $\upharpoonleft\downharpoonright$ spin pairs, it takes a lot of energy to flip a single nucleon’s spin between $\pm1/2$ in these even-even nuclei — usually several MeV, much higher than the $1\sim100$ keV ($10^7\sim10^9$ K) seen in stellar interiors.

Luckily, there are a few stable isotopes inside stars with an odd nucleon left unpaired to mediate a lower-energy M1 transition: Two that Haxton & Lee found were $^{57}_{26}$Fe, with an odd neutron and a 14.4 keV transition, and $^{23}_{11}$Na, with an odd proton and a 440 keV transition. Their potential to steadily emit axions in place of photons at stellar fusion temperatures earned these nuclei the nickname of nuclear axion pumps. In contrast to nuclear bremsstrahlung, which would produce astrophysical axions over a continuous energy range, nuclear axion pumps would emit them at narrow line energies. This enables searches for photons converted from axions (in the presence of either an astrophysical or artificial magnetic field) to fine-tune their sensitivity to very specific frequencies.

As it was well known that trace amounts of relic $^{57}_{26}$Fe from supernova explosions always find themselves sunken into the cores of stars, including our Sun — current models estimate a solar abundance of $\sim 3\times10^{-5}\,M_\odot$ — projects like the CERN Axion Solar Telescope (CAST) soon began using these iron axion pumps to place limits on the axion-nucleon couplings [An09]. Less attention was paid to the sodium axion pump: Although also present in trace amounts in the Sun, $^{23}_{11}$Na’s 440 keV M1 state is much harder to excite in the $1.5\times 10^7$ K ($1.3$ keV) solar core than $^{57}_{26}$Fe’s 14.4 keV.

But Wick noticed another source of heated $^{23}_{11}$Na: Massive stars burning $^{12}_6$C + $^{12}_6$C above $5\times10^8$ K (40 keV) produce excited $^{24}_{12}$Mg nuclei, which occasionally decays into $^{23}_{11}$Na instead of $^{20}_{10}$Ne + $^4_2$He. Compared to axion fluxes from relic solar abundances, the greater temperature and active $^{23}_{11}$Na production could make up for the greater distances and shorter lifespans of these massive stars. However, prior to the 2000s, the exact stellar mass threshold for carbon ignition was unclear: For stars born with zero-age main-sequence mass $\sim8\, M_\odot$, carbon is expected to ignite off-center, perhaps multiple times, before the flame propagates to the entire core — a difficult behavior to simulate [Ta13] Neither was the amount of $^{23}_{11}$Na produced certain — various codes point to $0.02\sim0.03\,M_\odot$ in $10\sim15\,M_\odot$ progenitors [La76].

Fortunately, increasing computing power and more sophisticated stellar evolution codes have since improved our understanding of $^{23}_{11}$Na production. Off-center carbon ignitions now reliable start for stars with initial masses between $7\sim9.5\,M_\odot$ [Li24]. In addition, after incorporating many more intermediate-mass isotopes and nuclear reactions, including weak interactions, these newer codes predict a slightly greater $^{23}_{11}$Na abundance — between $0.05\sim0.06\,M_\odot$ during carbon burning [Zh19].

This enables us to perform a few order-of-magnitude estimates. From nuclear physics, the expected axion emission rate for a single carbon-burning star is [Ha25]

\[ \tag{1} \Gamma_a \sim \prn{\f{g_{aNN}^\text{eff}}{10^{-9}}}^2 \cdot \f{M_{23\text{Na}}}{0.05\,M_\odot} \cdot \prn{\f{T}{5\times 10^8\text{ K}}}^9 \cdot 10^{44}\text{ axions/s}, \]

where $g_{aNN}^\text{eff} = g_{app} + 0.0638\,g_{ann}$ is a combination of the axion-proton and axion-neutron couplings, currently constrained to $\lesssim10^{-9}$ for axion masses below several $\mu$ev; As expected, it is more sensitive to the odd proton in $^{23}_{11}$Na. With a galactic magnetic field strength of a few $\mu$G, and an axion-photon coupling currently constrained to $\lesssim10^{-11}$ GeV$^{-1}$, the axion-photon conversion probability is approximately $0.01\%$ for axion masses below the neV scale [Ca24]. As a result, a single carbon-burning star at galactic distances $\sim8$ kpc should, at current levels of constraints for the various axion couplings, leave a sharp 440 keV gamma-ray line on Earth with a flux of about $0.1$ photons/cm$^2\cdot$s.

Now, according to intermediate-mass evolution studies done by Ken‘ichi Nomoto’s group, the threshold between whether a star becomes an O-Ne white dwarf or undergoes electron-capture supernova is an initial mass somewhere between $8\sim8.5\,M_\odot$. Applying the standard Salpeter initial mass function, this means $75\sim85\%$ of carbon-burning stars to eventually go supernova [Li24]. In the Milky Way, we have historically observed a supernova every $50\sim100$ years. As carbon burning typically lasts $\sim$10 thousand years, this means there should be a couple hundred sodium axion pumps active in our galaxy at any moment, boosting the all-sky flux to several photons/cm$^2\cdot$s at current coupling constraints. This is well above the sensitivity of upcoming wide-area gamma-ray space telescopes such as the Compton Spectrometer and Imager (COSI) experiment. We are therefore motivated to investigate this possibility with more rigor.

Getting started with MESA

Modules for Experiments in Stellar Astrophysics (MESA) is a suite of open-source Fortran libraries that, when chained together, can perform pretty darn good 1D star simulations. As a pre-/co-requisite for understanding its outputs, I recommend thoroughly reading a modern stellar astrophysics textbook such as Dina Prialnik’s An Introduction to the Theory of Stellar Structure and Evolution [Pr09], or taking an upper-division stellar astrophysics course. I would say the rough level of necessary proficiency is

• A detailed, qualitative understanding of the evolutionary stages of 1, 5, 10, 15, and 20 $M_\odot$ stars, and
• How to plot their evolution on $(\log T,\log L)$ and $(\log T,\log \rho)$ planes.

(I have to admit, I didn’t know about the $(\log T,\log\rho)$ plane until halfway into this research project, to my own detriment.) Beyond these basic skills, depending on how willing you are to blindly follow instructions, and the specificities of your interests, it may be helpful to also know about the effects of different initial metalicities, wind and mass loss models, nuclear reaction networks, pulsation mechanisms, etc, but to each their own.

To install MESA, I found it fairly straightforward to follow their official documentation. I ran into a few missing packages on Linux (Ubuntu 22.04.05 LTS) which were quickly fixed by searching on StackExchange. To become familiarized with MESA’s interface, I worked through every section in the “Using Mesa” tutorials. These tutorials are great for getting a hang of the basics.

After the tutorials, I recommend working through MESA’s test suite, which contains example configurations for simulating a wide collection of common stellar phenomena, from white dwarfs and supernova, to starspots pulsating variables. I believe running through some projects in the test suite is an excellent complement to any stellar astrophysics course. They really give you a feeling that the stars are evolving right before your own eyes, especially if you set pgstar, the built-in plotter, to the highest refresh rate.

As for data analysis, MESA lets us specify which parts of its simulation outputs we want saved to disc. Broadly speaking, there are three kinds of data; From smallest to largest typical file sizes, they are:

• Data series with one datapoint per time-slice, such as total mass, surface luminosity, effective temperature, etc; Configured in the project’s history_columns.list file and logged as the history.data file.

• Pre-made plots by pgstar, usually configured in the project’s inlist_pgstar file, which you can optionally display live and/or save to a directory.

• Data series with one datapoint per layer in the 1D star interior, per time-slice; Includes enclosed mass, local chemical abundances, nuclear burning power, local temperature, etc; Configured in the project’s profile_columns.list file and logged as profile$N$.data files, where $N$ is the index of the time-slice.

We use MESA’s recommended third-party Python library, mesa_reader, to read the data series files into Python, then process them with your garden-variety numpy, scipy and matplotlib.