# Cosmological constraints on decaying axion-like particles: a global analysis

Csaba Balázs, Sanjay Bloor, Tomás E. Gonzalo, Will Handley, Sebastian Hoof, Felix Kahlhoefer, Marie Lecroq, David J.E. Marsh, Janina J. Renk, Pat Scott, et al.

► To cite this version:

Csaba Balázs, Sanjay Bloor, Tomás E. Gonzalo, Will Handley, Sebastian Hoof, et al. Cosmological constraints on decaying axion-like particles: a global analysis. JCAP, 2022, 12, pp.027. 10.1088/1475-7516/2022/12/027. hal-03693170

HAL Id: hal-03693170

https://hal.science/hal-03693170

Submitted on 17 Mar 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
PAPER • OPEN ACCESS

# Cosmological constraints on decaying axion-like particles: a global analysis

To cite this article: Csaba Balázs et al JCAP12(2022)027

View the [article online](#) for updates and enhancements.

You may also like

- [GRB 221009A: A Light Dark Matter Burst or an Extremely Bright Inverse Compton Component?](#)

M. M. González, D. Avila Rojas, A. Pratts et al.

- [Constraints on axion-like particles with different magnetic field models from the PKS 2155-304 energy spectrum](#)

Jia Bu and Ya-Ping Li

- [Physics beyond colliders at CERN: beyond the Standard Model working group report](#)

J Beacham, C Burrage, D Curtin et al.
# Cosmological constraints on decaying axion-like particles: a global analysis

Csaba Balázs,ª Sanjay Bloor,ᵇ Tomás E. Gonzalo,ᶜ,ᵈ

Will Handley,ᵉ, f Sebastian Hoof,ᵉ, g Felix Kahlhoefer,ᶜ, d

Marie Lecroq,ª,h David J.E. Marsh,ᵢ Janina J. Renk,ᵇ, j, k

Pat Scottᵇ, k and Patrick Stöckerᶜ, l

Received June 10, 2022

Revised October 14, 2022

Accepted November 15, 2022

Published December 19, 2022

<figure><img src="image_2.png" /></figure>

<footnote>ªSchool of Physics and Astronomy, Monash University,
Melbourne, VIC 3800, Australia</footnote>

<footnote>ᵇDepartment of Physics, Imperial College London, Blackett Laboratory,
Prince Consort Road, London, SW7 2AZ, U.K.</footnote>

<footnote>ᶜInstitute for Theoretical Particle Physics and Cosmology (TTK), RWTH Aachen University,
D-52056 Aachen, Germany</footnote>

<footnote>ᵈInstitute for Theoretical Particle Physics (TTP), Karlsruhe Institute of Technology (KIT),
76128 Karlsruhe, Germany</footnote>

<footnote>ᵉCavendish Laboratory, University of Cambridge,
JJ Thomson Avenue, Cambridge, CB3 0HE, U.K.</footnote>

<footnote>fKavli Institute for Cosmology, University of Cambridge,
Madingley Road, Cambridge, CB3 0HA, U.K.</footnote>

<footnote>gInstitut für Astrophysik und Geophysik, Georg-August-Universität Göttingen,
Friedrich-Hund-Platz 1, 37077 Göttingen, Germany</footnote>

<footnote>hInstitut d'Astrophysique de Paris,
98 bis boulevard Arago, 75014 Paris, France</footnote>

<footnote>iTheoretical Particle Physics and Cosmology Group, Department of Physics,
King's College London, Strand, London, WC2R 2LS, U.K.</footnote>

<footnote>jOskar Klein Centre for Cosmoparticle Physics, AlbaNova University Centre,
SE-10691 Stockholm, Sweden</footnote>

<footnote>kSchool of Mathematics and Physics, The University of Queensland,
St. Lucia, Brisbane, QLD 4072, Australia</footnote>

<footnote>lInstitute for Astroparticle Physics (IAP), Karlsruhe Institute of Technology (KIT),
Hermann-von-Helmholtz-Platz 1, D-76344 Eggenstein-Leopoldshafen, Germany
E-mail: tomas.gonzalo@kit.edu, hoof@kit.edu, kahlhoefer@kit.edu,
stoecker@physik.rwth-aachen.de</footnote>
**Abstract.** Axion-like particles (ALPs) decaying into photons are known to affect a wide range of astrophysical and cosmological observables. In this study we focus on ALPs with masses in the keV–MeV range and lifetimes between $10^4$ and $10^{13}$ seconds, corresponding to decays between the end of Big Bang Nucleosynthesis and the formation of the Cosmic Microwave Background (CMB). Using the CosmoBit module of the global fitting framework GAMBIT, we combine state-of-the-art calculations of the irreducible ALP freeze-in abundance, primordial element abundances (including photodisintegration through ALP decays), CMB spectral distortions and anisotropies, and constraints from supernovae and stellar cooling. This approach makes it possible for the first time to perform a global analysis of the ALP parameter space while varying the parameters of ΛCDM as well as several nuisance parameters. We find a lower bound on the ALP mass of around $m_a > 300$ keV, which can only be evaded if ALPs are stable on cosmological timescales. Future observations of CMB spectral distortions with a PIXIE-like mission are expected to improve this bound by two orders of magnitude.

**Keywords:** particle physics - cosmology connection, axions

**ArXiv ePrint:** 2205.13549
# Contents

<table><tr><td>1 Introduction</td><td>1</td></tr><tr><td>2 Production and decay of axion-like particles coupled to photons</td><td>4</td></tr><tr><td>3 Constraints and likelihoods</td><td>7</td></tr><tr><td>3.1 CMB power spectra</td><td>7</td></tr><tr><td>3.1.1 Modifications of N<sub>eff</sub> and η<sub>b</sub></td><td>8</td></tr><tr><td>3.2 Primordial nucleosynthesis</td><td>10</td></tr><tr><td>3.2.1 Likelihood implementation and uncertainties</td><td>11</td></tr><tr><td>3.2.2 The <sup>7</sup>Li problem</td><td>12</td></tr><tr><td>3.3 Spectral distortions</td><td>13</td></tr><tr><td>3.3.1 ALP-induced spectral distortions</td><td>13</td></tr><tr><td>3.3.2 Likelihood implementation</td><td>15</td></tr><tr><td>3.4 ALP decays after SN1987A</td><td>15</td></tr><tr><td>3.5 Stellar evolution and cluster counts</td><td>16</td></tr><tr><td>3.6 Other likelihoods</td><td>17</td></tr><tr><td>4 Results</td><td>17</td></tr><tr><td>4.1 Frequentist scans</td><td>17</td></tr><tr><td>4.2 Bayesian scans</td><td>23</td></tr><tr><td>5 Conclusions</td><td>25</td></tr><tr><td>A GAMBIT implementation</td><td>28</td></tr><tr><td>A.1 New models</td><td>28</td></tr><tr><td>A.2 New capabilities provided in CosmoBit and DarkBit</td><td>29</td></tr><tr><td>A.3 New and updated backend interfaces</td><td>31</td></tr><tr><td>B Details of implementation of entropy injection between BBN and recombination</td><td>32</td></tr><tr><td>B.1 Numerical implementation</td><td>33</td></tr><tr><td>C Details of the SN1987A limits from ALP decay</td><td>33</td></tr><tr><td>C.1 Reference frame transformations</td><td>34</td></tr><tr><td>C.2 Monte Carlo procedure</td><td>35</td></tr><tr><td>C.3 Results</td><td>35</td></tr><tr><td>C.4 Technical details</td><td>36</td></tr></table>

# 1 Introduction

Cosmological observables such as the Cosmic Microwave Background (CMB) have now been measured with such precision that they not only constrain the properties of particles in the primordial thermal plasma, but also even probe out-of-equilibrium species with sub-dominant cosmological abundances. For example, unstable particles that decay around the time of recombination leave an observable imprint on the CMB temperature anisotropies, even if
they change the energy density of the electron-photon plasma only at the level of one part per billion [1–3]. Earlier decays are in turn tightly constrained by observations of CMB spectral distortions [4, 5] and the primordial element abundances produced in Big Bang Nucleosynthesis (BBN) [6–9]. This raises the exciting possibility of using cosmology to search for the effects of decaying exotic particles across many different cosmological eras, ranging from one second to billions of years [10].

A particularly well-motivated class of exotic particles that may be probed in this way are axion-like particles (ALPs) [11, 12], which generically arise in theories with spontaneously broken global chiral symmetries and theories with compactified extra dimensions [13]. In analogy with the QCD axion proposed to solve the strong CP problem [14, 15], ALPs would have an approximate shift symmetry that explains both the smallness of the ALP mass and the weakness of its couplings in terms of the large scale of symmetry breaking. Although ALPs may in general have a wide range of different couplings [16, 17], the most phenomenologically interesting one is the coupling to photons. Indeed, there is a vast literature on how to constrain ALP-photon couplings using cosmological [18–21], astrophysical [22–24] and laboratory data [25, 26]. For more details, we refer the interested reader to reviews on QCD axion and ALP cosmology [27] and experimental search strategies [28], as well as a continuously updated online repository for limits and sensitivity forecasts [29].

Most studies of cosmological ALP constraints make one of two assumptions on the ALP abundance. Either ALPs are assumed to be produced by the realignment mechanism [22, 30–34] in such a way that their present-day abundance is consistent with the observed abundance of dark matter (DM) in the universe, or they are assumed to be in thermal equilibrium with the bath of Standard Model (SM) particles at some high temperature but then decay, depleting their abundance. Both of these scenarios turn out to be very tightly constrained, such that ALPs are often believed to be excluded unless they are either very light (with a mass at or below the keV-scale) or very heavy (with a mass at or above the GeV-scale) [20]. It has only recently been pointed out that a low reheating temperature can open up otherwise ruled-out parameter space for intermediate-mass ALPs, if realignment production is inefficient and their couplings are small enough to avoid attaining thermal equilibrium [35]. In this case the dominant contribution to ALP production arises from the freeze-in mechanism at temperatures close to the reheating temperature, and the resulting abundance may be small enough to evade cosmological constraints arising from their subsequent decays [20, 36].

In the present work we perform a detailed analysis of ALPs in the keV–GeV mass range with tiny couplings to SM photons, such that they never enter into thermal equilibrium in the early universe. While these ALPs may in principle be sufficiently long-lived to survive until the late universe, we focus on ALPs that decay between BBN and recombination. In this case the ALP decays affect a wide range of observables, such as

• the primordial element abundances through photodisintegration of the nuclei formed during BBN,

• the black-body spectrum of the CMB through the injection of electromagnetic energy into the plasma,

• the CMB anisotropies through modifications of the recombination history,

• the ratio of the photon and neutrino temperatures (encoded in the relativistic number of effective degrees of freedom) through a relative increase of the former,
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Summary plot of our results in the $m_a-g_{a\gamma}$ plane. We show a part of the cosmological ALP region commonly excluded by Horizontal Branch (HB) stars from cluster counts [40], supernova SN1987A (see ref. [41] and section 3.4), supernova energy deposition [42], ionisation fraction constraints [20], BBN and $N_{\text{eff}}$ [35], spectral distortions [20], comparison to the extragalactic background light (EBL) [20], and X-ray searches [20, 43–45] (we used data provided by ref. [29] for some of the contours). For nonthermal ALPs considered in this work, the grey-shaded region around the best-fitting point (star) is however not excluded.</figcaption></figure>

* auxiliary astrophysical observables, such as axion emission and decays from supernova SN1987A or stellar cooling.

We make use of the recently released CosmoBit module [37] of the GAMBIT framework [38, 39] to perform a combined analysis of all these observables. Exploring the combined likelihood function has the advantage of automatically incorporating any complementary information contained in the data. Similar analyses, such as the ones presented in refs. [20, 21], investigated the individual likelihoods instead. Compared to ref. [20] we include additional constraints from SN1987A and a more rigorous treatment of spectral distortions and other CMB constraints. In contrast to our work, ref. [21] focuses on models with lifetimes $\tau_a \lesssim 10^7$ s and not all possible types of spectral distortions. In addition, our analysis improves upon the results of ref. [35] by performing a more accurate calculation of the ALP abundance as a function of the reheating temperature and by including state-of-the-art likelihoods for all cosmological data sets — in particular CMB spectral distortions and anisotropies.

A visual summary of our findings in the mass-coupling plane is shown in figure 1. We identify large regions of viable parameter space where ALPs would constitute a sizeable fraction of DM in the very early universe, but decay away before recombination. Considering regions of low reheating temperature $T_R \sim 5$ MeV, which is still consistent with BBN, opens up parameter regions that are usually considered to be excluded by cosmology. Moreover, by combining the individual observations into a global likelihood function we find that ALP decays can in fact improve the agreement between BBN and the CMB compared to $\Lambda$CDM.
Intriguingly, the corresponding parameter regions may be fully probed by future missions that measure CMB spectral distortions.

The remainder of this work is structured as follows. In section 2 we discuss the freeze-in production and subsequent decay of ALPs with tiny couplings to photons. The various ways in which these couplings can be constrained through observations and the corresponding likelihood functions are presented in section 3. Section 4 contains the setup of our parameter scans and the results that we obtain. These results, along with a GAMBIT input file to reproduce them, and a plotting script, can be downloaded from Zenodo [46].

## 2 Production and decay of axion-like particles coupled to photons

In this work we will focus on ALPs $a$ with an effective coupling to photons given by

$$ \mathcal{L} = \frac{g_{a\gamma}}{4} a F^{\mu\nu} \tilde{F}_{\mu\nu}, \qquad (2.1) $$

where $F^{\mu\nu}$ is the electromagnetic (EM) field strength tensor, $\tilde{F}^{\mu\nu} = \frac{1}{2}\epsilon^{\mu\nu\rho\sigma}F_{\rho\sigma}$ its dual, and the ALP-photon coupling $g_{a\gamma}$ has mass dimension $-1$.

**The ALP lifetime.** Axion-like particles are unstable and decay into two photons with a lifetime of $\tau_a = \Gamma_a^{-1}$, where

$$ \Gamma_a = \frac{g_{a\gamma}^2 m_a^3}{64\pi} = \frac{1}{1.32 \times 10^8 \text{s}} \left( \frac{g_{a\gamma}}{10^{-12} \text{GeV}^{-1}} \right)^2 \left( \frac{m_a}{10 \text{MeV}} \right)^3. \qquad (2.2) $$

We can broadly distinguish three different lifetime regimes:

1. ALPs decay before neutrino decoupling ($\tau_a \lesssim 1$ s). In this case the decay photons will fully thermalise and not affect any cosmological observables.

2. ALPs decay between neutrino decoupling and the present day ($1 \text{s} \lesssim \tau_a \lesssim 10^{17} \text{s}$). In this case ALPs may significantly affect various cosmological observables.

3. ALPs survive until the present day ($\tau_a \gtrsim 10^{17} \text{s}$). In this case ALPs contribute to the DM abundance; constraints in the mass range of interest arise from searches for mono-energetic X-ray or $\gamma$-ray photons.

The second case is particularly interesting because non-trivial constraints on the ALP parameter space can be expected. Indeed, the previous studies mentioned in section 1 showed that there are very tight constraints on EM decays of DM particles during or after recombination, or during BBN. Less attention has been given to the case of ALPs decaying between BBN and recombination, i.e. $10^4 \text{s} < \tau_a < 10^{13} \text{s}$, which will be the focus of this work.

**The ALP abundance.** In addition to $m_a$ and $g_{a\gamma}$, the two model parameters introduced above, the ALP abundance can be introduced as a third parameter. The abundance will be decisive for the subsequent analysis. Since in our set-up there are no processes that can produce ALPs in the temperature range that we will be interested in,<fnref n="1" /> the comoving number density $a^3n_a$ only changes through ALP decays. We can therefore define the would-be

<footnote n="1">¹A possible exception would be the inverse decay $\gamma\gamma \to a$. This process however turns out to be irrelevant in the viable parameter space of the model.</footnote>
present-day abundance in the absence of decays $n_{a,0}^{\tau_a \to \infty}$ and the corresponding would-be dark matter fraction

$$ \xi = \frac{m_a n_{a,0}^{\tau_a \to \infty}}{\rho_{\text{DM},0}}, \qquad (2.3) $$

where $\rho_{\text{DM},0}$ is the present-day DM density. Because the comoving DM density is constant, the ALP fraction is temperature-independent in the absence of decays. When including ALP decays, however, the actual present-day abundance of ALPs becomes completely negligible for the lifetimes considered in this work. Hence, $\xi$ is not directly constrained by the observed DM abundance, and it is possible to have $\xi > 1$.

To study the cosmological effects of ALPs at the keV–MeV scale, ref. [20] assumes that ALPs enter thermal equilibrium in the early universe and subsequently undergo thermal freeze-out. In this case, which has also been considered in ref. [21], $\xi$ can be calculated in terms of the other model parameters and, typically, $\xi \gg 1$.

Here we instead consider the case where ALPs’ interactions are sufficiently weak that they never thermalise. This assumption implies an upper bound on the reheating temperature $T_R$: because ALPs couple to photons through an effective interaction, the production rate of ALPs for $T \gg m_a$ scales as $\Gamma_{\text{prod}} \sim \alpha T^3 g_{a\gamma}^2$, where $\alpha \approx 1/137$ is the EM fine-structure constant. Thermalisation occurs if $\Gamma_{\text{prod}} > H$, where $H = \dot{R}/R$ is the Hubble rate in terms of the scale factor of the universe $R$.<fnref n="2" /> During radiation domination the Hubble rate depends only on the total energy density in radiation $\rho_{\text{rad}}$:

$$ H^2(T) \approx \frac{\rho_{\text{rad}}}{3\bar{m}_{\text{P}}^2} = \frac{\pi^2}{90\bar{m}_{\text{P}}^2} g_*(T) T^4 \approx 0.11 g_*(T) \frac{T^4}{\bar{m}_{\text{P}}^2}, \qquad (2.4) $$

where $g_*$ denotes the number of effective relativistic degrees of freedom and $\bar{m}_{\text{P}} \equiv M_{\text{P}}/\sqrt{8\pi} \approx 2.435 \times 10^{18}$ GeV is the reduced Planck mass. Hence, to first approximation, no thermalisation occurs if $g_{a\gamma}$ satisfies

$$ g_{a\gamma} \lesssim 10^{-7} \text{GeV}^{-1} \left( \frac{5 \text{MeV}}{T_R} \right)^{\frac{1}{2}}. \qquad (2.5) $$

If ALPs never enter into thermal equilibrium, their abundance depends on the details of how inflation and reheating proceed, and on the entire subsequent cosmological evolution. In particular, there may be a relevant contribution to the ALP abundance from the vacuum realignment mechanism, which depends on the initial misalignment angle of the ALP field, or from the decay of cosmic strings. For all practical purposes one can therefore treat the ALP fraction $\xi$ as a free parameter.

When doing so, it is however essential to realise that $\xi$ cannot be made arbitrarily small. The reason is that there is an irreducible<fnref n="3" /> contribution to the ALP abundance from the freeze-in mechanism, i.e. the direct production of nonthermal ALPs from the thermal bath of SM particles via the Primakoff effect [20, 36]: $\gamma\ell^\pm \to a\ell^\pm$, where $\ell^\pm$ denotes the charged

<footnote n="2">²We denote the scale factor by $R$ instead of $a$ to avoid confusion with the axion field.</footnote>

<footnote n="3">³We emphasise that this argument implicitly assumes that $g_{a\gamma}$ is constant during the epochs of interest. If the ALP is generated by breaking a U(1) symmetry, the associated radial modes have masses $m_r \gtrsim 1/g_{a\gamma}$, and are thus very heavy. In the case of low $T_R$ we therefore expect that these modes (and hence $g_{a\gamma}$) do not evolve. Nevertheless, if $g_{a\gamma}$ is allowed to vary, the ALP abundance may be suppressed further. The investigation of such scenarios requires a more detailed description of reheating, involving e.g. an early phase of matter domination, which is beyond the scope of this work.</footnote>
leptons.<fnref n="4" /> Let us therefore briefly review the freeze-in production of ALPs and calculate the corresponding contribution to ξ.

The freeze-in production is dominated by temperatures close to $T_R$, so that the final abundance of ALPs will be proportional to $T_R$ [36, 47]. The resulting bound on $T_R$ is substantially stronger than the one obtained from the requirement of nonthermalisation via eq. (2.5). Conversely, as $T_R$ cannot be arbitrarily low, we obtain an irreducible contribution to the abundance of ALPs in the early universe.

Let us focus on temperatures below the muon mass $m_μ$ for simplicity and consider the case that the ALP mass is negligible: $m_a \ll T \ll m_μ \sim 100$ MeV. The Boltzmann equation governing the evolution of the ALP number density $n_a$ is then given by [36]

$$ \frac{dn_a}{dt} + 3Hn_a = \frac{\alpha \zeta(3) g_{a\gamma}^2 T^6}{12\pi^2} \left[ 2 \ln \left( \frac{T}{m_\gamma} \right) + 0.8194 \right], \quad (2.6) $$

where $\zeta(3) \approx 1.202$ is Apéry's constant and $m_γ^2 = 4\pi\alpha T^2/9$ is the plasmon mass for a fully ionised universe.<fnref n="5" /> Defining $x = T_R/T$ and $Y_a = n_a/s$ with $s = \frac{2\pi^2}{45} g_*^s(T)T^3$ and $g_*^s = g_* = 10.75$ in the MeV range, this becomes

$$ \frac{dY_a}{dx} = 7.82 \times 10^{-5} \frac{g_{a\gamma}^2 \bar{m}_P T_R}{x^2}. \quad (2.7) $$

This equation can be readily integrated from $x = 1$ to $x \to \infty$, to yield

$$ \xi = 7.82 \times 10^{-5} \frac{m_a g_{a\gamma}^2 \bar{m}_P T_R s_0}{\rho_{\text{DM},0}} = 0.022 \left( \frac{m_a}{1 \text{ MeV}} \right) \left( \frac{T_R}{5 \text{ MeV}} \right) \left( \frac{g_{a\gamma}}{10^{-10} \text{ GeV}^{-1}} \right)^2, \quad (2.8) $$

where $s_0 = 2891.2 \text{ cm}^{-3}$ denotes the present-day entropy density. In particular, we find that $\xi \propto T_R$. The requirement that BBN proceeds as usual implies the lower bound $T_R \gtrsim 5$ MeV [48, 49], which in turn leads to a lower bound on $\xi$ as a function of $m_a$ and $g_{a\gamma}$.

In the derivation of eq. (2.8) we made a few simplifications that are not necessarily valid. In particular, we treated both ALPs and electrons as massless. However, when $m_a > T_R$, ALP production via the Primakoff process becomes exponentially suppressed, and the lower bound on $\xi$ is significantly relaxed. This effect can be approximately captured by multiplying eq. (2.8) by a Boltzmann suppression factor $e^{-m_a/T_R}$.

To more accurately include the effect of finite masses, we also implemented the ALP model in micrOMEGAs v5.2.11. Figure 2 provides a comparison between the ALP fraction computed by multiplying eq. (2.8) by $e^{-m_a/T_R}$, and that computed by micrOMEGAs. Both approaches include quantum statistics and the plasma mass of the photon. Eq. (2.8) includes a more careful treatment of the contribution of soft momenta from ref. [36], whereas micrOMEGAs includes the sub-leading contribution from $e^+e^- \to a\gamma$. Despite these differences, eq. (2.8) and the results from micrOMEGAs agree very well for $m_e$, $m_a \ll T_R$. At higher masses, the curves are qualitatively similar, but diverge. In order to accurately capture the suppression of ALP production for large masses, we will always use the result from micrOMEGAs in what follows. Since these calculations are computationally expensive, we perform them in advance and interpolate the tabulated results.

<footnote n="4">⁴Analogous Primakoff processes involving charged hadrons give a negligible contribution for the temperature range that we consider.</footnote>

<footnote n="5">⁵We emphasise that ALPs are typically produced with relativistic energy, such that it is not directly possible to infer the evolution of $\rho_a$ from eq. (2.6). Only at sufficiently late times, when $T \ll m_a$, can we make use of the relation $\rho_a = m_a n_a$.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** The irreducible freeze-in contribution to the would-be ALP dark matter fraction $\xi$ in the absence of decays as a function of $m_a$ (left) and $T_R$ (right) for a reference value of $g_{a\gamma} = 10^{-11}$ GeV$^{-1}$. We compare the numerical result from micrOMEGAs v5.2.11 (solid lines) with the approximate result from eq. (2.8) multiplied by a naïve Boltzmann suppression (dashed lines).</figcaption></figure>

To conclude this section, we note that the relative contribution of ALPs to the total energy density of the universe, $\rho_a(T)/\rho_{tot}(T)$, will be largest shortly before they start decaying, i.e. for $T \approx T_{dec}$ with $H(T_{dec}) = 1/\tau_a$. Assuming that these decays happen after neutrino decoupling but before matter-radiation equality, we can estimate

$$ \frac{\rho_a(T_{\text{dec}})}{\rho_{\text{tot}}(T_{\text{dec}})} = \xi \frac{\rho_{\text{DM}}(T_{\text{dec}})}{\rho_{\text{rad}}(T_{\text{dec}})} = \xi \frac{\Omega_{\text{DM}}}{\Omega_{\text{rad}}} \frac{T_{\text{dec}}}{T_0} \sim 10^{-4} \xi \left( \frac{\tau_a}{10^4 \text{s}} \right)^{1/2}. \quad (2.9) $$

We will see below that in the allowed regions of parameter space this ratio is always several orders of magnitude below unity, such that ALPs never dominate the expansion history of the universe. Nevertheless, their decays can have a number of profound effects, which we will discuss in the following section.

# 3 Constraints and likelihoods

In this section we discuss the various constraints on the ALP model and the corresponding likelihood functions, which we implement within the GAMBIT software framework.

## 3.1 CMB power spectra

We vary all parameters of the cosmological concordance model ΛCDM: the present-day Hubble expansion rate $H_0$; the present-day baryon abundance $\omega_{b,0} = \rho_{b,0}h^2/\rho_{\text{crit},0}$, where $H_0 = 100h \text{ km s}^{-1} \text{ Mpc}^{-1}$ and $\rho_{\text{crit},0} = 3\bar{m}_P^2 H_0^2$; the present-day dark matter abundance $\omega_{\text{DM},0}$; the redshift of reionisation $z_{\text{reio}}$;<fnref n="6" /> the amplitude $A_s$; and the spectral index $n_s$ of the

<footnote n="6">⁶Unlike the parametrisation used in the initial CosmoBit release [37], $\tau_{\text{reio}}$ is replaced with $z_{\text{reio}}$, as the latter is the quantity that is internally used within the calculations of CLASS and $\tau_{\text{reio}}$ is a derived quantity. The parametrisation used here avoids numerical instabilities which can occur when $\tau_{\text{reio}}$ is used as primary parameter in scenarios of strong exotic energy injection.</footnote>
initial power spectrum. We fix the temperature of the CMB<fnref n="7" /> to $T_{\text{CMB}} = 2.725\,48$ K, the sum of the neutrino masses to $\sum m_{\nu} = 0.06$ eV and the pivot scale $k_0 = 0.05$ Mpc$^{-1}$.

For a given set of parameters we use CLASS [50] to calculate the CMB power spectra. We then compare these predictions to data using the Planck 2018 TT,TE,EE+lowE+lensing likelihoods (referred to as the “baseline” likelihood in ref. [51], see also ref. [37]), which are publicly available through plc [52]. For high multipoles we use the ‘lite’ version of the likelihood, which requires only one additional nuisance parameter, $A_{\text{Planck}}$. As done by the Planck collaboration, we apply a Gaussian prior on this parameter [52]. This is done through an additional Gaussian likelihood on this parameter, in order to allow for both, Frequentist and Bayesian analyses. The details of the interface between GAMBIT, CLASS and plc were presented in ref. [37].

Crucially, we include the effects of ALP decays on the CMB power spectra. For $\tau_a > 10^{12}$ s, ALPs are still decaying when the electrons and protons in the primordial plasma combine to form neutral hydrogen. In this case the energy injected through these decays affects the free electron fraction and thereby delays recombination. The resulting constraints on decaying dark matter subcomponents were studied in e.g. ref. [3]. The interfaces to the CLASS extension ExoCLASS, as well as to all modern versions of CLASS (v3.0 onwards) in GAMBIT provide a simple way to calculate these constraints for a wide range of models, based on transfer functions calculated with DarkAges [53]. The transfer functions used in this code, however, are designed to derive the impact of energy injection at redshifts around recombination and later, but not at the redshifts prior to recombination, at which spectral distortions primarily occur. At these redshifts, the universe is opaque and all injected energy will be transferred to heat up the plasma, such that it suffices to encode all impact of energy injection in a single effective efficiency factor $f_{\text{eff}}$. As we assume the ALP to decay into photons exclusively, this efficiency factor is assumed to be $f_{\text{eff}} = 1$. To include the constraints from energy injection in the present work we use the interface between GAMBIT and of CLASS, which was first developed in ref. [37], further refined in ref. [54], and which we updated to work with the latest version of CLASS (v3.1).

For $\tau_a < 10^{12}$ s the dominant effect of ALP decays is a change of the photon energy density relative to those of neutrinos and baryons, which will be discussed in more detail in the following.

### 3.1.1 Modifications of $N_{\text{eff}}$ and $\eta_b$

In the absence of any beyond-the-SM physics, the only particle species that are still relativistic after BBN are photons and neutrinos. Hence the radiation energy density is given by

$$ \rho_{\text{rad}}(T_{\gamma}, T_{\nu}) = \frac{\pi^2}{15} \left( T_{\gamma}^4 + 3 \times \frac{7}{8} T_{\nu}^4 \right), \qquad (3.1) $$

where the first term accounts for the contribution from photons and the second term accounts for the contribution from neutrinos. As the two species are thermally decoupled, their respective temperatures $T_{\gamma}$ and $T_{\nu}$ evolve independently. However, because all SM processes that could modify either of the two heat baths are inefficient after neutrino decoupling, the respective temperatures are only subject to redshifting. Their ratio is therefore fixed at neutrino decoupling. In the instantaneous-decoupling approximation, the ratio’s value

<footnote n="7">⁷We note that the tiny uncertainty of the CMB temperature, $\sigma(T_{\text{CMB}}) = 5 \times 10^{-4}$ K, does become relevant in the context of CMB spectral distortions. This uncertainty is automatically marginalised over when we calculate the spectral distortion likelihood (see below).</footnote>
can be derived by assuming conservation of entropy in both heat baths separately, such that $T_ν^{inst} = (4/11)^{1/3} T_γ$. Non-instantaneous neutrino decoupling leads to a slightly larger temperature, and it is thus convenient to define the effective number of neutrino species via [55–58]

$$N_{\text{eff}} = \frac{3}{(\frac{4}{11})^{4/3}} \left(\frac{T_ν}{T_γ}\right)^4 = 3.0440 \pm 0.0002, \qquad (3.2)$$

such that eq. (3.1) can be rewritten as

$$\rho_{\text{rad}}(T_{\gamma}) = \frac{\pi^2}{15} \left[ 1 + \frac{7}{8} N_{\text{eff}} \left( \frac{4}{11} \right)^{\frac{4}{3}} \right] T_{\gamma}^4. \qquad (3.3)$$

For ALPs with a lifetime $\tau_a > 10^4$ s, the value of $N_{\text{eff}}$ during BBN will be equal to its SM prediction, given that the contribution of ALPs produced via the freeze-in mechanism to the total energy density during radiation domination is completely negligible. Subsequent ALP decays into photons will however increase $T_γ$ relative to $T_ν$ and thereby decrease the value of $N_{\text{eff}}$ inferred from the CMB anisotropies. Furthermore, a modification of the photon temperature between BBN and recombination changes the baryon-to-photon ratio $\eta_b = n_b/n_γ$. Including ALP decays, the evolution of the photon temperature is given by [7]

$$\frac{dT_γ}{dt} = \frac{15 m_a n_a(t) T_γ^{-3}}{4\pi^2 \tau_a} - H(t) T_γ, \qquad (3.4)$$

where $n_a(t)$ is the ALP number density, given by

$$n_a(t) = n_a(t_{\text{ini}}) \exp \left( -3 \int_{t_{\text{ini}}}^{t} dt' H(t') \right) e^{-t/\tau_a}, \qquad (3.5)$$

where $t_{\text{ini}}$ defines the initial time, which we take to be $10^3$ s, i.e. sufficiently below the lower bound $\tau_a \ge 10^4$ s. The first term of eq. (3.4) accounts for the temperature modification through the injection of entropy, whereas the second term accounts for the adiabatic cooling of the photon heat bath as the universe expands.

As $H(t)$ depends on both temperatures $T_γ$ and $T_ν$, it is important to also track the evolution of the neutrino temperature $T_ν$ as function of time. As ALPs do not interact with neutrinos, the temperature evolution is given by adiabatic cooling. Solving the corresponding differential equation yields

$$T_{\nu}(t) = T_{\nu}(t_{\text{ini}}) \exp \left( - \int_{t_{\text{ini}}}^{t} dt' H(t') \right). \qquad (3.6)$$

We solve the evolution of the photon temperature in eq. (3.4) numerically up to a final time $t_{\text{fin}} = 2 \times 10^{12}$ s, corresponding approximately to the time of matter-radiation equality.<fnref n="8" /> In appendix B we provide details of our numerical implementation.

As the contribution from entropy injection is strictly positive, photons do not cool as fast as the neutrinos, and $N_{\text{eff}}$ decreases between BBN and recombination. For $t_{\text{ini}} \le t \le t_{\text{fin}}$ the value of $N_{\text{eff}}$ is

$$N_{\text{eff}}(t) = N_{\text{eff, BBN}} \left( \frac{T_{\nu}(t)}{T_{\gamma}(t)} \right)^4 \left( \frac{11}{4} \right)^{4/3}, \qquad (3.7)$$

<footnote n="8">⁸Note that our calculation assumes that $\tau_a < t_{\text{fin}}$. For larger ALP lifetimes, the CMB will be directly sensitive to the exotic energy injection, and the resulting constraints are much stronger than those obtained from the modification of $N_{\text{eff}}$.</footnote>
where $N_{\text{eff,BBN}} \equiv N_{\text{eff}}(t_{\text{ini}})$ is the value of $N_{\text{eff}}$ at the end of BBN, given by eq. (3.2). Since CMB observations are in agreement with the SM prediction, we can exclude parameter points that predict too large deviations from this value. To calculate the resulting constraints, we pass the $\Lambda$CDM parameters and the calculated value of $N_{\text{eff}}(t_{\text{fin}})$ to CLASS.

Furthermore, a change in $T_\gamma$ also implies a change in the comoving photon number density and hence in the baryon-to-photon ratio. As $\eta_b(t_{\text{fin}}) \approx \eta_{b,0}$ is tightly constrained, we invert the calculation to obtain the baryon-to-photon ratio for earlier times $t_{\text{ini}} \le t \le t_{\text{fin}}$:

$$ \eta(t) = \eta(t_{\text{fin}}) \left( \frac{T_{\gamma}(t_{\text{fin}})}{T_{\gamma}(t)} \right)^3 \exp \left( 3 \int_{t}^{t_{\text{fin}}} dt' H(t') \right). \quad (3.8) $$

The baryon-to-photon ratio during BBN is then simply given by $\eta_{b,\text{BBN}} = \eta_b(t_{\text{ini}})$.

Note that the discussion above can be easily generalised to the case of additional contributions from dark radiation during BBN, viz. $N_{\text{eff,BBN}} = 3.044 + \Delta N_{\text{eff}}$. In this case the energy density in dark radiation scales exactly as the neutrino energy density, and the ratio $N_{\text{eff,BBN}}/N_{\text{eff,CMB}}$ is independent of the initial value of $N_{\text{eff,BBN}}$.

## 3.2 Primordial nucleosynthesis

Axion-like particles with $\tau_a > 10^4$ s decay after BBN, and the synthesis of primordial elements thus proceeds as in standard cosmology. Still, BBN constraints do play an important role for this model, as the decays of heavy ALPs with masses above $m_a \gtrsim 1$ MeV may contribute to the photodisintegration of some of the light elements [9]. The evolution of the primordial element abundances is thus affected by the injected EM spectrum through ALP decays and can be computed by solving the Boltzmann equation [59]

$$ \frac{dY_N}{dt} = \sum_j \int_0^\infty dE f_\gamma (Y_j \sigma_{j\gamma \to N} - Y_N \sigma_{N\gamma \to j}), \quad (3.9) $$

where $N \in \{n,p,D,^3\text{H}, ^3\text{He}, ^4\text{He}, ^6\text{Li}, ^7\text{Li}, ^7\text{Be}\}$ are the primordial elements, $Y_N = n_N/n_b$ their abundances, $\sigma_p$ the cross sections for the various disintegration processes of the light elements, and $f_\gamma$ the phase space distribution of nonthermal photons produced through ALP decays. Photons with energy $E_\gamma \gtrsim m_e^2/(22T)$, where $m_e$ denotes the electron mass, disappear efficiently via $e^+e^-$ creation [60]. Hence, disintegration processes only become efficient for temperatures $T \lesssim m_e^2/(22E_{\text{th}})$, where $E_{\text{th}}$ denotes the threshold energy for each nuclear reaction. As typical threshold energies lie in the MeV range, photodisintegration becomes relevant only for $T \ll 10$ keV, well after the end of BBN. The various disintegration processes and their corresponding threshold energies are given in table 1.

In scenarios where these depletion processes are sufficiently strong, the abundances of primordial elements will significantly differ from the predictions of vanilla $\Lambda$CDM, which are known to match observations in pristine gas clouds and metal-poor stars fairly well. Any sizeable modification of their abundances due to efficient disintegration will therefore be in considerable tension with measurements.

The element abundances are typically measured with respect to the hydrogen abundance (with the exception of $^3\text{He}$), and we thus express the relevant element abundances as [D/H], [$^3\text{He}/\text{D}$], $Y_p = [^4\text{He}/\text{H}]$ and [$^7\text{Li}/\text{H}$]$^9$.

<footnote>$^9$We do not consider the abundances of other primordial elements such as $n, ^3\text{H}$ and $^7\text{Be}$ because they undergo further decay at late cosmological times and are totally depleted at the present time. We also do not consider $^6\text{Li}$, for which no reliable measurement of the primordial abundance exists.</footnote>
<table><thead><tr><th>Reaction</th><th><i>E</i><sub>th</sub> [MeV]</th></tr></thead><tbody><tr><td>D + <i>γ</i> → <i>n</i> + <i>p</i></td><td>2.22</td></tr><tr><td><sup>3</sup>H + <i>γ</i> → D + <i>n</i></td><td>6.26</td></tr><tr><td><sup>3</sup>H + <i>γ</i> → 2<i>n</i> + <i>p</i></td><td>8.48</td></tr><tr><td><sup>3</sup>He + <i>γ</i> → D + <i>p</i></td><td>5.49</td></tr><tr><td><sup>3</sup>He + <i>γ</i> → <i>n</i> + 2<i>p</i></td><td>7.12</td></tr><tr><td><sup>4</sup>He + <i>γ</i> → <sup>3</sup>H + <i>p</i></td><td>19.81</td></tr><tr><td><sup>4</sup>He + <i>γ</i> → <sup>3</sup>He + <i>n</i></td><td>20.58</td></tr><tr><td><sup>4</sup>He + <i>γ</i> → 2D</td><td>23.84</td></tr><tr><td><sup>4</sup>He + <i>γ</i> → D + <i>n</i> + <i>p</i></td><td>26.07</td></tr></tbody></table>

<table><caption>**Table 1.** Threshold energies of dissociation processes considered in this work.</caption><thead><tr><th>Reaction</th><th><i>E</i><sub>th</sub> [MeV]</th></tr></thead><tbody><tr><td><sup>6</sup>Li + <i>γ</i> → <sup>4</sup>He + <i>n</i> + <i>p</i></td><td>3.70</td></tr><tr><td><sup>6</sup>Li + <i>γ</i> → <sup>3</sup>H + <sup>3</sup>He</td><td>15.79</td></tr><tr><td><sup>7</sup>Li + <i>γ</i> → <sup>3</sup>H + <sup>4</sup>He</td><td>2.47</td></tr><tr><td><sup>7</sup>Li + <i>γ</i> → <sup>6</sup>Li + <i>n</i></td><td>7.25</td></tr><tr><td><sup>7</sup>Li + <i>γ</i> → <sup>4</sup>He + 2<i>n</i> + <i>p</i></td><td>10.95</td></tr><tr><td><sup>7</sup>Be + <i>γ</i> → <sup>3</sup>He + <sup>4</sup>He</td><td>1.59</td></tr><tr><td><sup>7</sup>Be + <i>γ</i> → <sup>6</sup>Li + <i>p</i></td><td>5.61</td></tr><tr><td><sup>7</sup>Be + <i>γ</i> → <sup>4</sup>He + <i>n</i> + 2<i>p</i></td><td>9.30</td></tr></tbody></table>

For most element abundances, the most up-to-date measurements can be found in [61], with the exception of $[^{3}\text{He}/\text{D}]$, for which we have computed our own improved estimation. To estimate $[^{3}\text{He}/\text{D}] = [[^{3}\text{He}/\text{H}]/[\text{D}/\text{H}]]$, we use the analysis of ref. [62], who estimate the protosolar abundance ratios based on their measures values today. To do so, we can use their estimate for protosolar $[^{3}\text{He}/\text{H}]$, and either an estimate for protosolar $[([^{3}\text{He} + \text{D})/\text{H}]$ or a corrected direct measurement of Jovian $[\text{D}/\text{H}]$. Unlike ref. [62], we use the weighted average of both methods, taking care when propagating the uncertainties. We find that $[^{3}\text{He}/\text{D}] = 0.82 \pm 0.11$. As the use of $[^{3}\text{He}/\text{D}]$ as a measurement is debated in the literature (see e.g. footnote 4 in ref. [63]), we only use it as an upper limit by assigning an optimal likelihood to values smaller than the estimate. With this value and the recommendations from ref. [61], the bounds that we use in this study are

$$ 
\begin{align}
[\text{D}/\text{H}] &= (2.547 \pm 0.025) \times 10^{-5} \nonumber \\
[^{3}\text{He}/\text{D}] &< 1.03 \text{ at } 95\% \text{ confidence level (CL)} \nonumber \\
Y_{\text{p}} &= 0.245 \pm 0.003 \nonumber \\
[^{7}\text{Li}/\text{H}] &= (1.6 \pm 0.3) \times 10^{-10}. \tag{3.10}
\end{align}
$$

### 3.2.1 Likelihood implementation and uncertainties

CosmoBit uses AlterBBN [64, 65] to compute an initial prediction of the primordial abundances of the light elements based on a given set of cosmological parameters. It then passes these abundances to ACROPOLIS [59]. The latter includes the effects of photodisintegration by computing the photon spectrum in the post-BBN thermal bath and solving the Boltzmann equation for each of the primordial element abundances. Details of the GAMBIT interface to ACROPOLIS can be found in appendix A.3.

Both AlterBBN and our interface to ACROPOLIS compute ratios of abundances to the hydrogen abundance $[\text{Y}/\text{H}]$, from which we calculate $[^{3}\text{He}/\text{D}] = [[^{3}\text{He}/\text{H}]/[\text{D}/\text{H}]]$ in order to match the measured value of $[^{3}\text{He}/\text{D}]$ in eq. (3.10). CosmoBit then computes a multi-dimensional Gaussian likelihood for all abundances, taking into account correlated uncertainties (see ref. [37] for details):

$$ -2 \ln \mathcal{L} = f(\mathcal{P}, \mathcal{O}) \mathcal{C}^{-1} f(\mathcal{P}, \mathcal{O})^T + (2\pi)^n \det(\mathcal{C}), \quad (3.11) $$

where $\mathcal{P}$ is the vector of predicted abundances, $\mathcal{O}$ the vector of observed abundances, and $\mathcal{C}$ the covariance matrix. The function $f(\mathcal{P}, \mathcal{O})$ is introduced to differentiate elements for
which the measurement is only an upper limit, such as $ ^{3}\text{He} $. Therefore, $ f(\mathcal{P}, \mathcal{O})_i = \mathcal{P}_i - \mathcal{O}_i $ for all elements except for $ ^{3}\text{He} $, for which $ f(\mathcal{P}, \mathcal{O})_{^{3}\text{He}} = \max(\mathcal{P}_{^{3}\text{He}} - \mathcal{O}_{^{3}\text{He}}, 0) $. This way, the likelihood becomes flat when the predicted $ [^{3}\text{He}/D] $ is smaller than the observed value, in line with our desire to implement the $ [^{3}\text{He}/D] $ measurement as an upper bound only.

We construct the covariance matrix $ \mathcal{C} $ by computing the uncertainties for each of the abundances independently, and then combining them with correlation information. We estimate the uncertainties on each abundance for each parameter point by taking the larger of the distances between the mean and the high/low operational modes of **AlterBBN**, i.e. $ \sigma_{Y_N} = \max(Y_N^{\text{high}} - Y_N^{\text{mean}}, Y_N^{\text{mean}} - Y_N^{\text{low}}) $. A detailed computation of the correlations is possible in **AlterBBN**, but is computationally too expensive for our purposes. Fortunately, the production of primordial elements in our model proceeds as in $ \Lambda $CDM. The correlations thus do not change significantly when varying the $ \Lambda $CDM parameters within their typical ranges. We therefore employ a constant correlation matrix through our analysis that matches the one computed by **AlterBBN** with default settings. We then propagate the covariance matrix when computing the corrections after photodisintegration with **ACROPOLIS**, using the transfer matrix provided (see appendix A.3 for details). Finally, in eq. (3.11) we use a corrected covariance matrix, which matches the change to the observed ratio of abundances for $ [^{3}\text{He}/D] $.

Besides their obvious use as predictions to compare with observations, the primordial abundances are also relevant for cosmological probes after BBN. In particular the primordial helium abundance $ Y_{\text{p}} $ affects the physics of recombination and thus is an input to the computation of the CMB power spectrum with **CLASS** (see section 3.1). As recombination happens at temperatures much lower than the energy threshold for the photodisintegration reactions for $ Y_{\text{p}} $, we make sure to include these effects on the predicted value of $ Y_{\text{p}} $ that we use for computing the spectrum of the CMB.

In addition to the BBN likelihood, we also include a nuisance likelihood for the neutron lifetime $ \tau_n $, as the reaction rates that produce the primordial elements are strongly affected by uncertainties in $ \tau_n $. In this analysis we use the laboratory measurement of $ \tau_n $ from ultra-cold neutron ("bottle") experiments, $ \tau_{n,\text{bottle}} = (879.4 \pm 0.6) \text{ s} $, due to their smaller uncertainties. Effects of the uncertainty on $ \tau_n $ on constraints from BBN, the CMB and other measurements were studied in ref. [37].

### 3.2.2 The $ ^{7}\text{Li} $ problem

The BBN predictions for the abundance of $ ^{7}\text{Li} $ are known to be in disagreement with observations [66]. This is usually attributed to either systematic uncertainties in the determination of the primordial lithium abundance from astrophysical systems or uncertainties in the cross sections of nuclear reactions used in the theoretical calculations. Measurements of $ ^{7}\text{Li} $ are therefore generally not considered robust enough to be used for fits of cosmological data, and we do not include them in our scans unless explicitly stated otherwise.

Nevertheless, it is still interesting to consider physics beyond the SM as a possible solution of the $ ^{7}\text{Li} $ problem. For example, ref. [9] suggests that the tension between predictions and observations can be resolved via ALP decays. The idea in this work is to consider ALPs in the mass range $ 3.17 \text{ MeV} < m_a < 4.44 \text{ MeV} $ such that the photons produced in their decay have enough energy for the photodisintegration of $ ^{7}\text{Be} $, but not enough for deuterium. We will return to this discussion in section 4.
## 3.3 Spectral distortions

While the CMB radiation closely follows a black-body spectrum with characteristic temperature $T_{\text{CMB}} \approx 2.725$ K [67, 68], spectral distortions (SDs) exist due to interactions taking place during the entire thermal history of the universe [4, 69–75]. As the relevance of the underlying physical processes that induce SDs changes over time, their redshift dependence has to be taken into account. This means that SDs are — at least in theory — sensitive to any known or new physics that alters the energy budget of the expanding universe, thus making SDs a useful, complementary probe of e.g. inflation, dark matter, or new particles (see refs. [5, 76] for recent reviews).

Spectral distortions affect the photon spectrum $f(x, z)$, where $x = \omega/T_{\text{CMB}}$ is the dimensionless frequency as observed today and $z$ is the redshift. This can be expressed as a deviation from the perfect black-body shape, i.e.

$$f(x, z) = \frac{1}{e^x + 1} + \Delta f(x, z). \qquad (3.12)$$

Historically, SDs have been classified into different types, such as y-type distortions (from the Kompaneets or Sunyaev-Zeldovich parameter $y$) or chemical potential-type $\mu$ distortions. However, it is sufficient to specify the heating rates of the primordial plasma as a function of redshift, and integrate their spectral intensity $\Delta I \equiv (xT_{\text{CMB}})^3\Delta f(x, z)/2\pi^2$ to compare to observations today. Given the frequency-integrated heating rate $\dot{Q}$ from the Boltzmann equation, we have (see ref. [5])

$$\Delta I \propto \frac{\Delta \rho_{\gamma}}{\rho_{\gamma}} = \int_{0}^{\infty} dz \frac{\dot{Q}}{(1+z)H\rho_{\gamma}}. \qquad (3.13)$$

### 3.3.1 ALP-induced spectral distortions

Axion-like particles can induce SDs through ALP-photon conversion in the (known) Galactic magnetic field or a putative large-scale magnetic field of primordial origin [77–84], or through ALP decays [19, 20, 85]. The former process is only relevant for ALPs with $m_a \ll$ eV. Since those ALPs are much lighter than the ones considered in this work, we only focus on the effect of decaying ALPs.

The effective energy injection rate of decaying ALPs, which contributes to the heating rate in eq. (3.13), is given by [5]

$$\frac{dQ_{\text{inj},a}}{dt} = f_{\text{eff}} \xi \rho_{\text{DM}} \Gamma_a e^{-\Gamma_a t}, \qquad (3.14)$$

where we set the efficiency factor $f_{\text{eff}} = 1$. This is a good assumption for decays before recombination [5]. In particular, $f_{\text{eff}} = 1$ is justified because the energy injection happens due to EM interactions, and we use the “on-the-spot approximation” of instantaneous energy injection at the time of the decay.

The ALP-induced SDs can in theory be calculated entirely via the unified approach of a full integration over the cosmological history outlined in eq. (3.13). However, this is in practise a very expensive computation. The authors of ref. [5] hence implemented a number of approximation schemes to simplify the calculation, which we briefly summarise in the next section. Examples for the ALP-induced spectral distortion signal in terms of $\Delta I$ are shown in the left panel of figure 3. There we choose ALP parameters that illustrate the power of using the general approach adopted in our analyses.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Left: predicted SDs (coloured lines) vs uncertainties of FIRAS observations (grey-shaded region). Right: predicted $\mu$ and $y$ distortions for ALPs with $m_a = 100$ MeV, $\xi = 5 \times 10^{-4}$, and $T_R = 5$ MeV as a function of $\tau_a$. Horizontal dotted lines indicate the COBE/FIRAS limits and PIXIE forecasted limits; SDs above these lines are excluded at 95% CL. The vertical dashed lines indicate the limits from $y$-type SDs (red) and combined amount of SDs (gold), respectively.</figcaption></figure>

However, it is still useful to understand how ALPs can induce SDs by altering the single [86] and double [87] Compton rates. Compton scattering is efficient at redshifts $z \gtrsim 5 \times 10^4$, where $\mu$-type distortions are sourced. The name is due to the fact that these SDs introduce a small effective chemical potential term, often denoted by $\mu$, in the blackbody spectrum of eq. (3.12). Compton scattering then becomes inefficient at redshifts $z \lesssim 5 \times 10^4$, leading to $y$-type distortions [88]. These arise as Compton scattering becomes too slow to maintain a thermal equilibrium, leading to photon energy distributions different from the $\mu$-type SDs discussed before.

Redshifts of $z \sim 5 \times 10^4$ correspond to ALP lifetimes of $\tau_a \sim 10^{10}$ s, indicating that the type of ALP-induced SDs can changes in different parts of the parameter space. Spectral distortions of a given type can thus only constrain the regions of parameter space where the underlying physical interactions are relevant.

If they apply, the magnitudes of ALP-induced $\mu$ [4] and $y$ distortions [19] (see also refs. [5, 20, 89]) can be approximated as

$$ \mu \approx \frac{1}{\kappa_{\mu}} \left( \frac{3\rho_a}{\rho_{\gamma}} - \frac{8n_a}{n_{\gamma}} \right) \Big|_{T=T_{\text{dec}}} , \quad y \approx \frac{1}{4} \ln \left( \frac{\rho_a}{\rho_{\gamma}} + 1 \right) \Big|_{T=T_{\text{dec}}} , \qquad (3.15) $$

with $\kappa_{\mu} \equiv 4\zeta(2)/\zeta(3) - 3\zeta(3)/\zeta(4) \approx 2.14$ and assuming that all ALPs decay at a temperature $H(T_{\text{dec}}) \sim \Gamma_a$.

If ALPs decay too fast or are too massive, they become inefficient in sourcing SDs, reducing the SDs in eq. (3.15). The more accurately predicted ALP-induced $\mu$ and $y$ distortions are shown in the right panel of figure 3. They should be compared to the constraints from COBE/FIRAS, which are $|\mu| < 0.9 \times 10^{-4}$ and $|y| < 1.5 \times 10^{-5}$ and at 95% CL [67, 90, 91]. We also included the expected limits for a null detection from a PIXIE-like experiment [5, table 1].
While µ-type SDs typically set stronger constraints, this need not always be the case, as shown in figure 3. Here we chose ALP parameters where only y-type SDs result in a limit of $\tau_a < 3.1 \times 10^{10}$ s (vertical dashed red line). However, by deriving constraints using the total distortion $\Delta I$ (shown in the left panel of figure 3), we can even deduce a stronger limit of $\tau_a < 2.1 \times 10^9$ s (vertical dashed gold line). This illustrates the power of the spectral likelihood approach presented in the next section.

### 3.3.2 Likelihood implementation

The authors of ref. [5] implemented the theoretical prediction of SDs as a part of CLASS v3.0 [50], incorporating adjusted versions of ExoCLASS [53], CosmoTherm [75, 92], and SZpack [93, 94]. They also provide a template likelihood for spectral distortions as a part of MontePython [95, 96], with specific realisations for COBE/FIRAS [67] data and PIXIE [97] forecasts.<fnref n="10" />

We pass the default options to MontePython for the SD likelihoods; in particular, this means that we use the 'exact' scheme for the SD branching approximation method. The 'exact' approximation scheme is based on the Green's function method and a principle component analysis, which requires frequency binning [92, 98].

After specifying a frequency binning for each experiment, the likelihood for the observed distortions of photon spectral intensities $\Delta I_{\text{obs}}(\omega_i)$ at frequencies $\omega_i$ in the $i$th bin is

$$ \ln \mathcal{L} = -\frac{1}{2} \sum_{i} \left( \frac{\Delta I_{\text{obs}}(\omega_i) - \Delta I_{\text{pred}}(\omega_i)}{\delta I(\omega_i)} \right)^2, \qquad (3.16) $$

where the $\Delta I_{\text{pred}}(\omega_i)$ are the theory predictions by CLASS and the $\delta I(\omega_i)$ are the expected sensitivities, which are the sum in quadrature of the experimental noise level and theory uncertainties. This includes foreground maps, which were accounted for in ref. [99] and largely based on ref. [100].

One of the strengths of implementing the likelihood with the computational framework described above is that we avoid the complications of using eq. (3.15) for $\mu$ and $y$ distortions separately. Instead, we constrain all types of SDs simultaneously, using essentially all spectral information contained in the data. Our results should thus supersede previous works such as refs. [20, 85].

## 3.4 ALP decays after SN1987A

The authors of ref. [101] set limits on $m_a$ and $g_{a\gamma}$ from ALPs emitted from core-collapse supernova SN1987A, using data [102] of the gamma-ray spectrometer onboard the Solar Maximum Mission (SMM) [103]. The (instantaneous) emission spectrum of ALPs with energy $E_a$ is

$$ \frac{dN_a}{dE_a} = C \frac{E_a^2}{e^{E_a/T_{\text{SN}}} - 1} \sigma_P(m_a, g_{a\gamma}, E_a), \qquad (3.17) $$

with $T_{\text{SN}} = 30.6$ MeV [101]. Here $C = 2.54 \times 10^{77}$ MeV$^{-1}$ is a fit to the results of ref. [104], and $\sigma_P$ is the Primakoff cross section [20], for which we use a Debye screening scale of $k_S = 16.8$ MeV [101]

<footnote n="10">¹⁰The PIXIE likelihood is based on the instrument's expected sensitivity and sets the data equal to the expected signal without any new, exotic physics (a null detection). For FIRAS, the authors of ref. [5] use the observed sensitivity but effectively also assumes a null detection. We expect this approximation to be reasonable as the observed residuals are small compared to the uncertainties shown in the left panel of figure 3.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Left: geometry for SN1987A ALPs decaying into two photons, emitted at angles $\alpha_{1/2}$ at a distance $L_a$ in the reference frame of the SMM satellite. Image credit for the SSM satellite to G. Nelson/NASA (JSC image library; public domain). Right: our profile likelihood and exclusion region of $m_a$ and $g_{a\gamma}$ (dashed black line) compared to that of ref. [101] (blue line) at $3\sigma$ CL, assuming 1 degree of freedom.</figcaption></figure>

The limits arise because the emitted ALPs can subsequently decay into two photons at some distance $L_a$, following an exponential distribution with mean distance

$$ \ell = \frac{\beta\gamma}{\Gamma_a} = \frac{64\pi}{g_{a\gamma}^2 m_a^3} \sqrt{\frac{E_a^2}{m_a^2} - 1}, \qquad (3.18) $$

where $\beta$ and $\gamma$ are the dimensionless ALP velocity and Lorentz boost factor, respectively.

To calculate the expected photon fluence from ALP decays, we estimate the “success rate”, i.e. the fraction of photons arriving at the SMM satellite, via Monte Carlo simulations. The expected fluence is then calculated by integrating eq. (3.17), dividing by the surface of the sphere centred on the supernova with the distance to the detector as the radius ($d_{SN} = 51.4$ kpc), and multiplying by the success rate. An illustration of the setup is provided in the left panel of figure 4.

Following ref. [101], we use a Gaussian likelihood function assuming a null result where the standard error is derived from the observed photon fluence in the energy range of 25–100 MeV. The resulting exclusion limit is shown in the right panel of figure 4. We provide details of our implementation in appendix C, as well as a table of the signal prediction that ships as part of the latest version of GAMBIT, and publicly available Python scripts for our Monte Carlo procedure.<fnref n="11" />

## 3.5 Stellar evolution and cluster counts

Due to their weak interactions with matter, ALPs can carry away energy from stellar interiors, affecting the evolution of e.g. Horizontal Branch (HB) stars by reducing their lifetimes [105–107]. This allows one to place a limit on $g_{a\gamma}$ by comparing the number of HB stars $N_{HB}$ to the number of stars on the Red Giant Branch (RGB) $N_{RGB}$ in globular clusters (GCs). The

<footnote n="11">¹¹Available on Github at https://github.com/marie-lecroq/ALP-fluence-calculation, tagged v1.0.</footnote>
associated “R parameter” observable $R_{GC} \equiv N_{HB}/N_{RGB}$ [108], is directly connected to the lifetimes and energy loss rates of these systems.

The theoretically expected R parameter can be estimated from stellar evolution simulations that include the energy loss from ALPs [24, 40, 108, 109],

$$R_{GC} \approx 0.022 - 0.443 (1 + 0.965 |g_{a\gamma}|)^{1/2} + 7.331 Y_{GC}, \quad (3.19)$$

where $Y_{GC} = 0.2515(17)$ [110] is an estimate for the Helium abundance in GCs. The observed value of the R parameter, $R_{GC,obs} = 1.39 \pm 0.03$ [40], is consistent with the axion-free value from (3.19), which is $R_{GC} = 1.42 \pm 0.01$.

Eq. (3.19) is only valid for ALP masses much smaller than the typical core temperature of an HB star, $m_a \ll T_{HB} \sim 10$ keV. For greater ALP masses, the Primakoff interaction cross section, which leads to the energy loss in HB stars, should be multiplied by a correction factor $F(m_a)$. Using eq. (A.5) of ref. [20], we obtain $F(m_a)$ normalised such that $F(m_a \ll T_{HB}) \simeq 1$. As the Primakoff cross section and the energy loss are proportional to $g_{a\gamma}^2$, it is sufficient to replace $g_{a\gamma} \mapsto g_{a\gamma} \times \sqrt{F(m_a)}$ in eq. (3.19) in order to apply the correction.

The R parameter likelihood is part of the DarkBit [111] module of GAMBIT. It was first introduced in ref. [112], and the correction for larger masses was added in ref. [113]. Also note that slightly stronger limits in the region around $m_a \sim 100$ keV may be obtained by considering photon coalescence [114]. As the R parameter constraints only marginally affect our results, we do not include these slightly improved constraints.

## 3.6 Other likelihoods

Finally, we make use of MontePython v3.5 to include the BOSS DR12 [115] likelihood for Baryon Accoustic Oscillations and the Pantheon [116] likelihood for type Ia supernovae. The latter likelihood requires the introduction of an additional nuisance parameter M corresponding to the absolute magnitude of a type Ia supernova, see ref. [37] for details.

## 4 Results

By combining all of the likelihoods introduced in section 3, we are now in the position to perform a global scan of ALP decays in the early universe. We first present the results from a frequentist analysis of the parameter space, which enables us to identify particularly interesting parameter regions, and then perform a Bayesian analysis, which provides complementary information on the preference between the ALP model and ΛCDM.

### 4.1 Frequentist scans

The parameters included in our scans are summarised in table 2. We emphasise that the parameter ranges for the ALP model are somewhat arbitrary and are chosen in such a way as to capture all relevant features.<fnref n="12" /> In other words, while viable parameter regions exist also beyond the scan ranges, they do not lead to any qualitative changes in our conclusions. The scan ranges for the ΛCDM and other nuisance parameters are chosen so that they fully

<footnote n="12">¹²We note that for the smallest ALP masses and the shortest lifetimes that we consider, the assumption of non-relativistic ALP decays is not a good approximation for ALPs produced via freeze-in, because it underestimates the energy density of ALPs. Nevertheless, we find that this parameter region is robustly excluded in our scans. A more accurate treatment would only strengthen this exclusion and is therefore unnecessary.</footnote>
<table><caption>**Table 2.** List of model and nuisance parameters included in our scans and their corresponding scan ranges.</caption><thead><tr><th>Model parameter</th><th></th><th>Scan range</th><th></th></tr></thead><tbody><tr><td>ALP mass</td><td>m<sub>a</sub></td><td>[0.001, 200]</td><td>MeV</td></tr><tr><td>ALP lifetime</td><td>τ<sub>a</sub></td><td>[10<sup>4</sup>, 10<sup>13</sup>]</td><td>s</td></tr><tr><td>ALP abundance</td><td>ξ</td><td>[10<sup>-12</sup>, 10<sup>2</sup>]</td><td></td></tr><tr><td>Baryon abundance</td><td>ω<sub>b</sub></td><td>[0.020, 0.024]</td><td></td></tr><tr><td>Dark matter abundance</td><td>ω<sub>DM</sub></td><td>[0.10, 0.13]</td><td></td></tr><tr><td>Hubble constant</td><td>H<sub>0</sub></td><td>[62, 74]</td><td>km s<sup>-1</sup> Mpc<sup>-1</sup></td></tr><tr><td>Redshift of reionisation</td><td>z<sub>reio</sub></td><td>[4.5, 9.5]</td><td></td></tr><tr><td>Primordial curvature</td><td>ln(10<sup>10</sup>A<sub>s</sub>)</td><td>[2.9, 3.2]</td><td></td></tr><tr><td>Scalar spectral index</td><td>n<sub>s</sub></td><td>[0.9, 1.1]</td><td></td></tr><tr><td>Neutron lifetime</td><td>τ<sub>n</sub></td><td>[875, 895]</td><td>s</td></tr><tr><td>Planck nuisance parameter</td><td>A<sub>Planck</sub></td><td>[0.9, 1.1]</td><td></td></tr><tr><td>Pantheon nuisance parameter</td><td>M</td><td>[-20, -18]</td><td></td></tr></tbody></table>

encompass the allowed parameter regions. Note that this statement is not in contradiction with our treatment of the reheating temperature, which we fix to $T_R = 5$ MeV. This choice is conservative since limits that depend on $T_R$ can typically be evaded in the region of low $T_R$. Thanks to the various parameter degeneracies of the ALP abundance — schematically captured in eq. (2.8) — and our choice of leaving $\xi$ as a free parameter, we can still access the parameter space studied in previous works such as ref. [35]. In terms of our numerical results, our choice of $T_R$ only has a mild effect: for example, going from $T_R = 5$ MeV to $T_R = 100$ MeV will change the lowest allowed ALP mass (at 95% CL) from $m_a \lesssim 0.25$ MeV to $m_a \lesssim 0.71$ MeV.

To perform comprehensive scans of the parameter space we employ the differential evolution sampler *Diver* v1.0.4 [117], using a population size of 20,000 and a convergence threshold of $10^{-8}$. We carried out additional scans targeted at specific regions of parameter space in order to ensure that the likelihood function is fully explored. As we perform a frequentist analysis based on the profile likelihood, we simply combine samples from different scans. The results shown below are based on a total of 33.5 million samples and took over 700,000 CPU hours on the Joliot-Curie HPC cluster at CEA (France).

In what follows, we use Wilks' theorem to determine the exclusion regions of the profile likelihood. To do so we assume that the number of degrees of freedom corresponds to the number of free parameters — i.e. two degrees of freedom for the two-dimensional plots shown below. It is well known that the assumptions of Wilks' theorem often do not hold exactly, and that care should be taken to ensure its validity when claiming a discovery [e.g. 118]. For the purpose of setting limits w.r.t the best-fitting point, however, we argue that Wilks' theorem typically presents a good approximation and avoids costly Monte Carlo simulations that do not alter significantly the findings of our analysis (see e.g. ref. [113, appendix D] for a demonstration in the context of global fits).
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Profile likelihood for various combinations of the parameters. Black lines correspond to the $1\sigma$, $2\sigma$, and $3\sigma$ contours around the best-fitting point (white star), whereas coloured lines are 99.7% CL limits from the individual likelihoods. Shaded grey regions are excluded from our study as they are beyond the scan ranges.</figcaption></figure>

Our main results are shown in figure 5, where we plot the profile likelihood as a function of different combinations of ALP model parameters. We remind the reader that there are three independent ALP model parameters, one of which is profiled out in each panel. We use grey shading to indicate parameter regions that are not explored in our scans due to the chosen scan ranges for $m_a$ and $\tau_a$. In particular, the upper bound on $\tau_a$ implies a lower bound on $g_{a\gamma}$ as a function of $m_a$; this, in turn, implies a lower bound on $\xi$ because of the irreducible freeze-in contribution, cf. eq. (2.8).

Each of the applied constraints is sensitive to a specific combination of model parameters and can therefore be understood most naturally in one of the panels in figure 5. For example, the astrophysical constraints are independent of $\xi$ and therefore directly constrain the combination ($m_a, g_{a\gamma}$) or, equivalently, ($\tau_a, g_{a\gamma}$). The cosmological constraints, on the other hand, are to first approximation independent of $m_a$. They are therefore most easily
expressed in terms of $(\xi, \tau_a)$. Coloured lines and labels indicate the locations of individual
constraints (at 99.7% CL) in the various panels.

We find that combining the various likelihoods excludes large parts of the ALP parameter
space. The cosmological constraints become stronger with increasing ALP lifetime and require
$\xi \ll 1$, unless $m_a \gtrsim 100$ MeV and $\tau_a \lesssim 10^5$ s. Moreover, we find that there are no viable
parameter regions for $m_a \lesssim 300$ keV as a result of the SN1987A constraint. This is because
it places an upper bound on $g_{a\gamma}$, and therefore a lower bound on $\tau_a$, which in turn leads to
strong cosmological constraints. We emphasise that this lower bound on $m_a$ is a direct result
of observational constraints and does not depend on the scan ranges that we have chosen.
However, when significantly expanding the scan ranges, we find that there is a second disjoint
region of allowed parameter space at smaller ALP masses, ALP lifetimes comparable to or
greater than the age of the universe, and a tiny abundance. The most promising way to probe
this parameter region would be through searches for X-ray lines [e.g. 20], which we do not
consider in the present work.

We observe that it is essential to combine different likelihoods in order to constrain
the parameter space from various directions. However, we find that no relevant constraints
are obtained from the modification of $N_{\text{eff}}$ and $\eta_b$ discussed in section 3.1.1. The reason is
that in the parameter regions allowed by all other constraints, the energy density of ALPs
relative to radiation is simply too small to change the photon temperature in an observable
way. Even for the shortest ALP lifetimes and largest ALP masses that we consider, eq. (2.9)
implies $\rho_a/\rho_{\text{tot}} < 10^{-3}$, which restricts the relative modification of $N_{\text{eff}}$ and $\eta_b$ to a similar
level. Furthermore, we find that once all other constraints are imposed, the likelihood for
$R_{\text{GC}}$ has no effect on the remaining parameter space, since in the allowed parameter region
the ALPs are too heavy to be produced in horizontal branch stars.

We find that large regions of parameter space feature a relatively constant likelihood,
which corresponds to the case that the ALP has negligible impact on cosmological observables.
However, we also find a small region that is slightly preferred by data. The best-fitting point
is found at $m_a = 126.1$ MeV, $\tau_a = 6.04 \times 10^9$ s (corresponding to $g_{a\gamma} = 3.31 \times 10^{-15}$ GeV$^{-1}$),
and $\xi = 1.48 \times 10^{-5}$ (compatible with the irreducible freeze-in contribution for $T_R \lesssim 4$ GeV)
and gives $-2\Delta \ln \mathcal{L} = 1.51$ compared to the case of negligible ALP abundance.

To understand the origin of this improvement, we show in figure 6 the allowed parameter
region in terms of the deuterium abundance D/H and the baryon abundance $\Omega_b h^2$. For
comparison, we show the measured $2\sigma$ CL interval for both quantities as horizontal and
vertical bands (see section 3.2 for the experimental measurement of D/H and ref. [51] for
that of $\Omega_b h^2$). On the one hand, we find that in ΛCDM these two quantities are tightly
correlated, and it is not possible to perfectly fit both measurements simultaneously. On the
other hand, in our ALP model it is possible to increase the deuterium abundance through
photodisintegration of $^4$He, which clearly improves the agreement with observations.<fnref n="13" />

We remind the reader that our scans include a theoretical uncertainty for the prediction
of D/H, which is calculated on the fly using AlterBBN. Indeed, for most points the relative
theoretical uncertainty is large compared to the observational uncertainty (see also ref. [37]).
This is why the allowed parameter region of our model extends well beyond the observed
confidence interval for D/H. The best-fitting point of ΛCDM is hence found to lie within
the $1\sigma$ allowed parameter region of the ALP model. This finding emphasises the need for an

<footnote n="13">¹³Note that the production of deuterium through photodisintegration of $^4$He dominates over the photodisintegration of deuterium itself, because $^4$He is much more abundant than deuterium and hence hard photons are more likely to be absorbed by the former than the latter.</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Left: profile likelihood in the plane of baryon density $\Omega_b h^2$ and deuterium abundance D/H. Black lines correspond to the $1\sigma$ and $2\sigma$ contours around the best-fitting point for the ALP model (white star). The brown line indicates the $2\sigma$ contour around the best-fitting point in $\Lambda$CDM (brown star). Horizontal and vertical lines mark the central values and the $2\sigma$ intervals of the measured $\Omega_b h^2$ [51] and D/H (see section 3.2) values, respectively. Right: time evolution of the ALP energy density (relative to the energy density of photons) as well as of the abundance ratios D/H and $^3$He/D for the best-fit point. Dashed (solid) red lines indicate the observed values ($2\sigma$ bounds) of the abundance ratios, the shaded bands represent theoretical uncertainties, which for the case of D/H are comparable to the measurement uncertainties.</figcaption></figure>

accurate treatment of uncertainties in the calculation of BBN abundances. Indeed, neglecting the theoretical uncertainty on D/H would significantly overestimate the preference of the ALP model over $\Lambda$CDM.

While the preference for our model over $\Lambda$CDM is certainly not significant, it is still interesting to understand which future observations may shed further light on this issue. We therefore performed another set of scans where we replace the COBE/FIRAS likelihood with a mock likelihood corresponding to the sensitivity of a PIXIE-like mission [97]. Such a PIXIE mock likelihood, where a null result for new physics was assumed to generate typical observed data, was first studied in ref. [5] and made publicly available through MontePython.<fnref n="14" /> The results are shown in figure 7 and are indicative of the improvements that can be expected from future measurements of CMB SDs. We find that a null result from PIXIE would indeed significantly tighten the constraints on our model and robustly exclude the best-fit point found in the previous scan. In other words, the best-fitting point of our main analysis including the COBE/FIRAS likelihood predicts CMB SDs that should be easily detectable

<footnote n="14">¹⁴Note that it is not currently possible in MontePython to combine the COBE/FIRAS likelihood and the PIXIE mock likelihood. However, since the latter is assumed to be substantially more constraining than the former, we expect results from a combination of both likelihoods to be very similar to the case with only the PIXIE mock likelihood. In any case, the PIXIE likelihood is subject to considerable unknowns at this point and should therefore only be considered an approximation.</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Profile likelihoods as in figure 5 but for the projected PIXIE limits instead of the FIRAS likelihood. Black lines again correspond to the $1\sigma$, $2\sigma$, and $3\sigma$ contours around the best-fit point. The likelihood in the allowed region is very flat such that the best-fit point is not unique; it is thus not shown here. The grey shaded region corresponds to the $2\sigma$ region around the best-fit point (grey star) of the scans with the FIRAS likelihood from figure 5.</figcaption></figure>

with future missions, such that PIXIE would be able to determine to high significance whether photodisintegration of $^4$He caused by ALP decays is responsible for the large value of D/H compared to the naive $\Lambda$CDM prediction. We emphasise that future measurements of CMB anisotropies (using for example CMB-S4 [119, 120]) would not lead to similar improvements in the constraints, since the predicted modifications of $\Delta N_{\text{eff}}$ and $\eta_b$ are simply too small to be observable.

Let us finally comment on the $^7$Li problem, which can potentially be addressed through the decays of ALPs with mass in the range $3.17 \text{ MeV} < m_a < 4.44 \text{ MeV}$ (see section 3.2.2). To understand if the ALP model considered in this work can solve this problem, we ran a separate scan in which the $^7$Li measurement is included as an additional likelihood. In agreement with the results from ref. [9], we find that the ALP abundances required to significantly modify
the ⁷Li abundance must be of order $\xi \sim 1$. Such large ALP abundances are however robustly excluded by CMB SDs, which were not considered in ref. [9]. We therefore conclude that, once all relevant constraints are taken into consideration, ALPs decaying into photons cannot resolve the ⁷Li problem.

To conclude, we note that we have checked explicitly that all results presented in this section are robust against additional contributions to $\Delta N_{\text{eff}}$ from further ultrarelativistic species. This is a direct consequence of the fact that the constraint on $N_{\text{eff}}$ does not dominate our exclusion limits.

## 4.2 Bayesian scans

The priors that we employ for Bayesian scans extend across the ranges quoted in table 2. For the $\Lambda$CDM parameters and the nuisance parameters, all of which are relatively tightly constrained by data, we employ flat priors on each interval.<fnref n="15" /> The ALP model parameters, on the other hand, can vary over many orders of magnitude. In order not to impose a specific scale, we employ logarithmic priors for these parameters.

These scans took place on the Cambridge CSD3 DiRAC cluster, with sampling by PolyChord [121] v1.20.1 with nlive= 1000, nprior= 40000, tol= $10^{-10}$ and logzero= $-10^{101}$. These settings ensure good sampling of the prior, the likelihood peak and the broader posterior, and give control over the points thrown away at the prior level.

We performed a $\Lambda$CDM run and an ALP run with the same priors on the cosmological parameters. The resulting posteriors from the two runs for the baseline cosmological parameters are indistinguishable and no correlations are found between the $\Lambda$CDM parameters and the ALP parameters. For the ALP run, we consider two effective priors on the ALP parameters: the default “box prior” given by the ranges in table 2, and a “Primakoff prior”, which results from throwing away points within the box prior where the irreducible ALP abundance from the Primakoff process (for $T_R$ = 5 MeV) is incompatible with the chosen value of $\xi$. We emphasise that in the case of the “box prior” we still impose the lower bound on $\xi$ at the level of the likelihood, so that the two approaches lead to indistinguishable posteriors.

Figure 8 shows the Bayesian equivalent of figure 5. The key difference here compared to figure 5 is that both the priors and the posterior are plotted in figure 8. Note that whilst these panels do show individual points, they do not show all points sampled. The points shown are an equal number of equal-weight draws from each distribution. This is in contrast to figure 5, which is based on profiling the likelihoods of all sample points within fixed-width bins, and interpolating the resulting maximum likelihood between them.

For a likelihood $\mathcal{L}$ and prior $\pi$, one may compute a Bayesian evidence (or marginal likelihood) $\mathcal{Z}$, posterior distribution $\mathcal{P}$ from Bayes theorem, and Kullback-Leibler divergence $\mathcal{D}_{\text{KL}}$ between the prior and posterior distribution [122]:

$$ \mathcal{P} = \frac{\mathcal{L}\pi}{\mathcal{Z}} \qquad \mathcal{Z} = \int \mathcal{L}\pi d\theta = \langle \mathcal{L} \rangle_{\pi} \qquad \mathcal{D}_{\text{KL}} = \int \mathcal{P} \ln \frac{\mathcal{P}}{\pi} d\theta = \left\langle \ln \frac{\mathcal{P}}{\pi} \right\rangle_{\mathcal{P}} \quad (4.1) $$

In order to quantitatively compare the two models, $\Lambda$CDM and $\Lambda$CDM +ALPs, we compute the Bayesian evidence, $\mathcal{Z}$ for the different prior definitions:

$$ \begin{aligned} \ln \mathcal{Z}_{\text{ALPs}}^{\text{box}} &= -1039.09 \pm 0.14 \\ \ln \mathcal{Z}_{\text{ALPs}}^{\text{Primakoff}} &= -1038.12 \pm 0.14 \\ \ln \mathcal{Z}_{\Lambda\text{CDM}} &= -1037.11 \pm 0.14 \end{aligned} \quad (4.2) $$

<footnote n="15">¹⁵In particular, we use a flat prior for $A_{\text{Planck}}$ and include a Gaussian likelihood in our scans to constrain it. This procedure leads to the same posteriors as when using a Gaussian prior with no additional likelihood.</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Bayesian scan analogue of figure 5. Light green points represent equally-weighted samples drawn from the default box prior, dark green points are samples drawn from the more restrictive Primakoff prior (see text for details), and light and dark blue points are samples drawn from the posterior using the COBE/FIRAS and PIXIE-type likelihoods, respectively.</figcaption></figure>

$\Lambda$CDM is thus marginally preferred to the ALP model in a Bayesian sense by $\Delta \ln \mathcal{Z}^{\text{box}} = -1.98 \pm 0.20$ or $\Delta \ln \mathcal{Z}^{\text{Primakoff}} = -1.02 \pm 0.20$ log units. The error bars in this section all arise from the nested sampling procedure [123, 124] as a result of the probabilistic integration used for estimating the evidence.

We can see where the difference in Bayesian evidence arises by using the Occam's razor equation to decompose the evidence into a goodness of fit (average log-likelihood) and Occam penalty (Kullback-Leibler divergence) [125]. By taking logarithms and rearranging Bayes theorem before posterior averaging, we recover the Occam's razor equation, which in our case
decomposes the evidences as

$$
\begin{align}
\ln \mathcal{Z} &= \langle \ln \mathcal{L} \rangle_p - \mathcal{D}_{\text{KL}} \nonumber \\
\ln \mathcal{Z}_{\text{ALPs}}^{\text{box}} &= -1012.27 - 26.82 \nonumber \\
\ln \mathcal{Z}_{\text{ALPs}}^{\text{Primakoff}} &= -1012.27 - 25.85 \nonumber \\
\ln \mathcal{Z}_{\text{ACDM}} &= -1012.38 - 24.72 \tag{4.3}
\end{align}
$$

These results demonstrate that all models in a Bayesian sense fit the data equally well, with $\Delta \langle \ln \mathcal{L} \rangle_p = 0.11 \pm 0.07$ relative to $\Lambda$CDM, consistent with zero. The ALPs are penalised by having additional constrained parameters, with an additional compression of $\Delta \mathcal{D}_{\text{KL}}^{\text{box}} = 2.09 \pm 0.19$ or $\Delta \mathcal{D}_{\text{KL}}^{\text{Primakoff}} = 1.12 \pm 0.19$, i.e. the parameter space contained in the posterior is a factor of $\exp(\Delta \mathcal{D}_{\text{KL}}^{\text{box}}) = 8.3 \pm 1.5$ or $\exp(\Delta \mathcal{D}_{\text{KL}}^{\text{Primakoff}}) = 3.1 \pm 0.6$ smaller than the prior. This compression factor measures how well the data constrain the model. Note that employing the more stringent Primakoff prior rules out a little more than half the parameter space in advance, and therefore reduces the Occam penalty. If the COBE/FIRAS data set is replaced with the PIXIE forecast, the only change is an increase in parameter space compression by approximately a factor of e: $\Delta \mathcal{D}_{\text{KL}}^{\text{box}} = 3.08 \pm 0.19$ and $\Delta \mathcal{D}_{\text{KL}}^{\text{Primakoff}} = 2.11 \pm 0.19$, i.e. the posterior is a factor of $\exp(\Delta \mathcal{D}_{\text{KL}}^{\text{box}}) = 22.2 \pm 4.2$ or $\exp(\Delta \mathcal{D}_{\text{KL}}^{\text{Primakoff}}) = 8.4 \pm 1.6$ smaller than the prior. The Likelihood contribution to the evidence remains the same.

Figure 9 shows that using the full scope of the nested sampling run, we can gain insight into the high likelihood region shown in figure 5. This plot showcases the Bayesian interplay between likelihood and parameter volume.

5 Conclusions

Axion-like particles (ALPs) are well-motivated candidates for light new physics and arise generically in many extensions of the Standard Model. In the present work we studied ALPs with masses in the keV–MeV range. We focused on ALPs with couplings to photons small enough to be inaccessible to accelerator experiments, but large enough to ensure that the ALPs decay before recombination and therefore do not contribute to the observed dark matter abundance. In this coupling range, ALPs can affect various cosmological and astrophysical observables, which makes it possible to search for them and to constrain their properties.

A central aspect of our analysis is that we did not assume that ALPs enter into thermal equilibrium in the early universe, such that their abundance is essentially an independent parameter. Nevertheless, one can place a lower bound on the abundance from the irreducible freeze-in production, and we provided a detailed calculation of this abundance as a function of the ALP mass and reheating temperature. At the same time, BBN and CMB constraints (both from anisotropies and spectral distortions) place tight upper bounds on the ALP abundance, typically requiring that these particles constitute only a tiny fraction of the dark matter abundance even before they decay.

To perform a global analysis of the ALP parameter space, we interfaced the GAMBIT framework with a number of external codes in order to combine state-of-the-art likelihood calculations for all relevant data sets. Most notably, for the first time we combine detailed calculations of the effects of photodisintegration of light elements with constraints on CMB spectral distortions, both of which give highly relevant constraints on ALP models. Our setup makes it possible to scan simultaneously over the ALP model parameters, the $\Lambda$CDM
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Top: live point locations in the $(\tau_a, \xi)$ plane during different stages of nested sampling (NS). The $i$th NS iteration may be parametrised by the log-prior volume $\ln X \approx -i/n_{\text{live}}$. Three critical stages are shown: the initial (non-Primakoff) prior points when $\ln X = 0$, at the “typical set” where $\ln X = -D_{\text{KL}}$, and at a much lower value around the point where NS terminates according to our selected tolerance (which in this case is around the best-fit point). The prior points have large volume and low likelihood, whereas the peak has high likelihood and low volume. The Bayesian “sweet spot” (i.e. majority contribution to the evidence) occurs around the typical set with approximate parameter volume given by $e^{D_{\text{KL}}}$. The higher-likelihood phase therefore occupies a parameter region $\approx 10^7$ times smaller than the usual posterior, so is not preferred in a Bayesian sense and does not appear in figure 8. Bottom: phase plot of log-likelihood against log-volume, showing the average log-likelihood $\langle \ln \mathcal{L} \rangle_{\mathcal{P}}$ over the posterior $\mathcal{P}$, the Kullback-Leibler divergence $D_{\text{KL}}$, and the maximum log-likelihood value $\ln \mathcal{L}_{\text{max}}$. The high likelihood peak can be identified as beyond the region where there is a small inflection point at $\ln X \approx -42$. All of this information is extractable in post-processing with the NS Python package anesthetic [126].</figcaption></figure>

parameters, and various nuisance parameters, in order to identify the most interesting regions of parameter space (see figure 5).

In particular, the decays of ALPs into photons can modify both the abundances of light elements and the baryon-to-photon ratio $\eta_b$. As a result, we find that [D/H] and $\Omega_b h^2$ are less tightly correlated in our model than in ΛCDM (see figure 6). This makes it possible to achieve a slightly better fit to observations in our model than in ΛCDM, although the Bayesian evidence does not prefer the ALP model over standard ΛCDM. However, we did not find any allowed parameter regions in our model that would address the $^7$Li problem.

Compared to previous results, we note that leaving the nonthermal abundance as a free parameter opens up parameter space for higher ALP masses, which are otherwise excluded by cosmological constraints (see figure 1). Nonetheless, bounds on the decays of ALPs
produced in SN1987A still provide relevant constraints in that region. This emphasises the importance of an analysis framework that allows for the inclusion of all relevant parameters and complementary constraints from astrophysics — both of which can be achieved with global-fitting tools such as GAMBIT.

Our findings also highlight the fact that measurements of the spectrum of the CMB from over 20 years ago can provide some of the leading cosmological constraints for ALPs with masses around the MeV scale. In particular, these constraints are stronger than those from the photodisintegration of light elements after BBN, which only constrain ALP masses $\gtrsim 5$ MeV. More accurate measurements from a future satellite mission such as PIXIE therefore promise substantially stronger constraints and the potential to discover weakly coupled MeV-scale ALPs via spectral distortions.

To conclude, we emphasise that our analysis was focused on a fairly restrictive set of ALP models, in which the ALP couplings to electrons and hadrons are assumed to be negligible. While this is a common assumption in the literature, there has been substantial, recent progress in understanding the effects of ALP-matter couplings on cosmology [127-130]. This research programme goes hand in hand with a detailed exploration of the importance of ALP-matter coupling for supernovae such as SN1987A [41, 131, 132]. Moreover, ALPs with hadronic couplings have been proposed as a possible solution of the $^7$Li problem [133]. It will therefore be very interesting to extend our setup and explore the interplay between such ALP models and observations, and the resulting constraints on more general ALP parameter spaces.

## Acknowledgments

We thank Deanna Hooper, Lloyd Knox, Julien Lesgourgues, Kai Schmidt-Hoberg, Edoardo Vitagliano, and all members of the GAMBIT community for discussions. We are thankful to Sebastian Wild for collaboration in the early stages of this work. We are furthermore very grateful to Marco Hufnagel for help with ACROPOLIS, to Matteo Lucca and Nils Schöneberg for help with the spectral distortions code, and to Joerg Jaeckel and Pedro C. Malta for helpful discussions on ref. [101] and letting us inspect the associated software code. CB received financial support from the Australian Research Council Discovery Projects DP180102209 and DP210101636. TEG, FK and PSt were funded by the Deutsche Forschungsgemeinschaft (DFG) through the Emmy Noether Grant No. KA 4662/1-1. SH was supported by the Alexander von Humboldt Foundation and the German Federal Ministry of Education and Research. DJEM is supported by an Ernest Rutherford Fellowship from the Science and Technologies Facilities Council (U.K.). PSc acknowledges funding support from the Australian Research Council under Future Fellowship FT190100814. PSt and SH thank FK for hospitality at RWTH Aachen, where part of this work was done. We thank PRACE for access to the Joliot-Curie HPC cluster at CEA. This work was performed using the Cambridge Service for Data Driven Discovery (CSD3), part of which is operated by the University of Cambridge Research Computing on behalf of the STFC DiRAC HPC Facility (https://dirac.ac.uk/). The DiRAC component of CSD3 was funded by BEIS capital funding via STFC capital grants ST/P002307/1 and ST/R002452/1 and STFC operations grant ST/R00689X/1. DiRAC is part of the National e-Infrastructure. This work used the Scientific Compute Cluster at GWDG, the joint data centre of Max Planck Society for the Advancement of Science (MPG) and University of Göttingen. This article made use of pippi v2.1 [134].
# A GAMBIT implementation

The computations and simulations perfomed in this work were done using the global fitting framework GAMBIT. In this appendix we document the modifications and additions we implemented for the purpose of this study. For full details on the structure of GAMBIT, we refer the reader to the GAMBIT core paper [38]. Further details on the CosmoBit extension can be found in refs. [37, 135].

## A.1 New models

The new models added to GAMBIT for this study are a generic parametrisation of a cosmological ALP, **GeneralCosmoALP**, two simplified ALP models **CosmoALP** and **CosmoALP_gg_tau**, a re-parametrisation of ΛCDM, **LCDM_zreio**, two new models for non-standard radiation, **etaBBN_rBBN_dNurBBN** and **rBBN_dNurBBN**, and two models of nuisance parameters, **cosmo_nuisance_BK15** and **cosmo_nuisance_SpectralDistortions**.

Various models of ALPs were already implemented in GAMBIT in ref. [112], including the **GeneralALP** model, which sits at the top of the model hierarchy. We generalised this base model even more to take into account cosmological effects, such as the thermal abundance of ALPs and the reheating temperature of the Universe, introducing a new “parent model” of **GeneralALP** named **GeneralCosmoALP**. We also added two other “child models” of **GeneralCosmoALP**: **CosmoALP** and **CosmoALP_gg_tau**. As ALPs may decay anytime before neutrino decoupling and the present day, translation functions are provided to the decaying DM model **DecayingDM_photon**, assuming that the primary decay mode is via photons. A “friendship relation” between these models ensures that energy injection to the CMB is taken into account in models with short ALP lifetimes.

**GeneralCosmoALP:** gagg, gaee, gaN, fa, ma0, Tchi, beta, thetau, f0_thermal, T_R

Generalization of the **GeneralALP** model from ref. [112], with the addition of the thermal ALP abundance $\xi$ (the f0_thermal parameter) and reheating temperature $T_R$.

**CosmoALP:** Cagg, fa, ma0, thetau, f0_thermal, T_R

Child model of **GeneralCosmoALP**, where the ALP-electron coupling and ALP-nucleon coupling are neglected, $g_{ae} = g_{aN} = 0$, and with temperature-independent mass, i.e. $\beta = 0$ and $T_\chi$ a large and irrelevant value. The ALP-photon coupling $g_{a\gamma}$ is replaced by its dimensionless equivalent $C_{a\gamma}$.

**CosmoALP_gg_tau:** tau, fa, ma0, Tchi, beta, thetau, f0_thermal, T_R

Child model of **GeneralCosmoALP** where the ALP-electron coupling and ALP-nucleon coupling are neglected, $g_{ae} = g_{aN} = 0$, and the ALP-photon coupling is replaced by the ALP lifetime $\tau_a$.

The standard parametrisation of ΛCDM in GAMBIT includes the baryon density $\omega_b$, the DM density $\omega_{DM}$, the Hubble constant $H_0$, the CMB temperature $T_{CMB}$, and the optical depth at reionisation $\tau_{reio}$. However, some of the external libraries we use, e.g. **CLASS**, internally use the redshift at reionisation $z_{reio}$ instead of $\tau_{reio}$. Typically obtaining $z_{reio}$ from $\tau_{reio}$ is not an issue, but in the presence of energy injection at the time of reionisation this is not trivial. For this reason we introduce a new model, **LCDM_zreio**.
**LCDM\_zreio:** omega_b, omega_cdm, H0, T_cmb, z_reio

Re-parametrisation of the ΛCDM model with the redshift at reionisation $z_{\text{reio}}$ instead of the optical depth at reionisation $\tau_{\text{reio}}$.

In GAMBIT, all non-standard radiation content that modifies the Hubble rate is parametrised by the five parameter model etaBBN_rBBN_rCMB_dNurBBN_dNurCMB and children thereof, cf. ref. [37, appendix B.4]. For this analysis we have implemented two additional models, etaBBN_rBBN_dNurBBN and rBBN_dNurBBN, which allow for scenarios in which the modifications to $N_{\text{eff}}$ at BBN are free model parameters that fully determine the respective values at recombination.

**etaBBN_rBBN_dNurBBN:** dNur_BBN, r_BBN, eta_BBN.

Parametrisation of non-standard radiation content with $\eta_{\text{BBN}}$, the baryon-to-photon ratio at the end of BBN; $\Delta N_{\text{ur, BBN}}$, the ultra-relativistic degrees of freedom at the end of BBN; and $r_{\nu, \text{BBN}}$, the neutrino temperature ratio at the end of BBN. A translation function is provided to map these parameters to those of etaBBN_rBBN_rCMB_dNurBBN_dNurCMB, requiring the new capability Neff_evolution, defined below.

**rBBN_dNurBBN:** dNur_BBN, r_BBN

Child model of etaBBN_rBBN_dNurBBN, where $\eta_{\text{BBN}}$ is inferred from the baryon-to-photon ratio today $\eta_0$, and thus the translation depends on the new capability eta_ratio, defined below.

Lastly, we introduced nuisance parameter models to include the systematic uncertainties of various experiments whose likelihoods are made available by MontePython. The nuisance model for the BICEP/Keck Array was updated to include the likelihood for the BK15 results [136] in the model cosmo_nuisance_BK15. Additionally, as this work is the first introduction of CMB spectral distortions in GAMBIT (cf. section 3.3), we implemented a model to include the nuisance parameters for the COBE/FIRAS and PIXIE experiments: cosmo_nuisance_SpectralDistortions.

**cosmo_nuisance_BK15:** BBbetadust, BBbetasync BBdust, BBsync, BBalphadust, BBTdust, BBalphasync, BBdustsynccorr, EEtoBB_dust, EEtoBB_sync, Delta_dust, Delta_sync, gamma_corr, gamma_95, gamma_150, gamma_220

Extension of the existing cosmo_nuisance_BK14 model including additional systematics.

**cosmo_nuisance_SpectralDistortions:** sd_delta_T, sd_T_D, sd_beta_D, sd_A_D, sd_T_CIB, sd Beta_CIB, sd_A_CIB, sd Alpha_sync, sd Omega_sync, sd A_sync, sd T_e, sd_EM, sd Nu_p_spin, sd_A_spin, sd_A_C0, sd_y_reio_nuisance

The spectral distortion nuisance model contains parameters that take into account various uncertainties on the COBE/FIRAS and PIXIE experiments.

## A.2 New capabilities provided in CosmoBit and DarkBit

In addition to the new models introduced for this study, we also added new capabilities, and modified existing ones, to work with the ALP models. Most of the new module functions that provide these capabilities were added to CosmoBit, as they had to do with the cosmological implications of ALPs, but one was added to DarkBit, as it deals with late-time astrophysical constraints. The full list of new capabilities, along with the module functions that provide
<table><caption>**Table 3.** Capabilities added to CosmoBit and DarkBit associated with this study. Here, `map_str_db1` is a `typedef` for `std::map<std::string, double>`, `map_str_int` for `std::map<std::string, int>` and `[pd]` is shorthand for `photodisintegration`. All functions are contained in CosmoBit except for `calc_PhotonFluence_SN1987A_Decay`, which is part of DarkBit.</caption><thead><tr><th>Capability</th><th>Function (type)</th><th>Dependencies [type] /<br>Backend reqs [type (args)]</th></tr></thead><tbody><tr><td>lifetime</td><td>lifetime_ALPagg (double)</td><td></td></tr><tr><td>minimum_abundance</td><td>minimum_abundance_ALP_analytical (double)<br>minimum_abundance_ALP_numeric (double)</td><td></td></tr><tr><td>minimum_fraction</td><td>minimum_fraction_ALP (double)</td><td>minimum_abundance [double]</td></tr><tr><td>DM_fraction</td><td>DM_fraction_ALP (double)</td><td>minimum_fraction [double]<br>lifetime [double]<br>RD_oh2 [double]</td></tr><tr><td>total_DM_abundance</td><td>total_DM_abundance_ALP (double)</td><td>DM_fraction [double]</td></tr><tr><td>external_dNeff_etaBBN</td><td>compute_dNeff_etaBBN_ALP (map_str_db1)</td><td>total_DM_abundance [double]<br>lifetime [double]</td></tr><tr><td>eta_ratio</td><td>eta_ratio_ALP (double)</td><td>external_dNeff_etaBBN [double]</td></tr><tr><td>Neff_evolution</td><td>Neff_evolution_ALP (map_str_db1)</td><td>external_dNeff_etaBBN [double]</td></tr><tr><td>tau_reio</td><td>get_tau_reio_classy (double)</td><td>class_get_tau_reio [double ()]</td></tr><tr><td>z_reio</td><td>get_z_reio_classy (double)</td><td>class_get_z_reio [double ()]</td></tr><tr><td>H_at_z</td><td>get_H_at_z_classy (daFunk::Funk)</td><td>class_get_Hz [double (double)]</td></tr><tr><td>time_at_z</td><td>get_time_at_z_classy (daFunk::Funk)</td><td>class_get_tz [double (double)]</td></tr><tr><td>age_universe</td><td>get_age_universe_from_time_at_z (double)</td><td>time_at_z [daFunk::Funk]</td></tr><tr><td>primordial_abundances_BBN</td><td>compute_primordial_abundances_BBN<br>(BBNcontainer)</td><td>AterBBNInput [map_str_db1]<br>call_nucl_err [int<br>(map_str_db1, double, double)]<br>call_nucl [int<br>(map_str_db1, double)]<br>get_NNUC [size_t ()]<br>get_abund_map_AlterBBN<br>[map_str_int ()]</td></tr><tr><td>primordial_abundances</td><td>primordial_abundances (BBNcontainer)<br>primordial_abundances_decayingDM<br>(BBNcontainer)</td><td>primordial_abundances_BBN<br>[BBNcontainer]<br>set_input_params [void<br>(bool, int, int, double)]<br>abundance_[pd]_decay [void<br>(double*, ...)]</td></tr><tr><td>PhotonFluence_SN1987A_Decay</td><td>calc_PhotonFluence_SN1987A_Decay (double)</td><td></td></tr></tbody></table>

them and their dependencies and backend requirements, can be seen in table 3. We give a quick summary of the most important ones below.

A few of the newly added capabilities involve computing properties of the ALPs that can be derived from the model parameters. The capability `lifetime` computes the ALP lifetime as a function of the mass and coupling. The `minimum_abundance` capability computes the minimal freeze-in abundance of ALPs, produced via Primakoff processes, either analytically or numerically (by reading in tables pre-computed with micrOMEGAs), and `minimum_fraction` computes the fraction of the total abundance of dark matter constituted by the minimal freeze-in abundance. The `total_DM_abundance` is the total abundance of ALPs, produced by realignment and freeze-in, and `DM_fraction` provides the fraction of the total DM abundance constituted by ALPs (assuming no ALP decays).
In decaying DM models such as ours, the baryon-to-photon ratio at BBN $\eta_{\text{BBN}}$, the neutrino-to-photon temperature ratio $r_\nu$, and the number of ultra-relativistic degrees of freedom at recombination $\Delta N_{\text{eff}}$, are provided by the capability `external_dNeff_etaBBN`. These values are returned as a `std::map`, and are can be accessed individually using the capabilities `eta_ratio` for $\eta_{\text{BBN}}$ and `Neff_evolution` for both $r_\nu$ and $\Delta N_{\text{eff}}$. Details of these computations can be found in appendix B.

As mentioned above in A.1, we have defined a re-parametrisation of the $\Lambda$CDM model in **GAMBIT** to more closely match the parametrisation in **CLASS**. Along with the new model, a few capabilities are defined to provide useful quantities. `tau_reio` and `z_reio` extract the values of the optical depth and redshift of reionisation from **CLASS**. The capabilities `H_at_z` and `time_at_z` return the Hubble parameter and time since the Big Bang from **CLASS** as `daFunk` functions of redshift (see [111] for more information on the `daFunk` library). Lastly `age_universe` computes the age of the Universe at a given redshift with **CLASS**.

To take into account the effects of the photodisintegration of light elements due to the EM spectrum injected by the ALP decays, we have restructured the capabilities relating to BBN. The capability `primordial_abundances_BBN` replaces the old `BBN_abundances` capability, which computes the abundances after BBN using **AlterBBN**. The new capability `primordial_abundances` computes the abundances at the present time. Two module functions provide this capability, `primordial_abundances`, which simply returns the abundances unchanged from BBN, and `primordial_abundances_decayingDM`, which uses **ACROPOLIS** to modify the computed abundances at BBN in order to include the effects of photodisintegration.

Finally, the only new capability added to **DarkBit** is `PhotonFluence_SN1987A_Decay`, which computes the $\gamma$ fluence from the supernova SN1987A as a result of the decay of ALPs into photons.

### A.3 New and updated backend interfaces

`GAMBIT` provides interfaces to a large collection of external tools, known as backends. The backend interfaces are documented along with the modules that use them [37, 111, 137-140]. Backend interfaces can also be auto-generated [141]. For this analysis we have updated some of the backends reported in [37] to their latest version, such as **CLASS 3.1** [50] and **MontePython 3.5** [96]. The updated version of **CLASS** brings along the possibility of including spectral distortions, as described in section 3.3, which also benefits from the new version of **MontePython**, as it comes with more likelihoods, including COBE/FIRAS and a mock PIXIE-like likelihood. Furthermore, we have modified the interface to **AlterBBN** [64, 65] in order to allow the computation of abundances without uncertainties, as well as to return three values of the abundances, corresponding to the central, high and low values computed by **AlterBBN**.

The only absolutely new backend interface implemented for this study is to the external tool **ACROPOLIS 1.2.1** [59]. Written in Python, **ACROPOLIS** computes the effects of photodisintegration of light elements on pre-computed abundances. In addition to the light element abundances, **ACROPOLIS** takes as input the mass of the decaying particle, $m_a$, its lifetime $\tau_a$, the ratio of the number density of the decaying particle today to that of photons, $N_{0,a}$, and the branching fractions to electrons and photons, respectively. From this set of input parameters, the backend convenience function `abundance_photodisintegration decay` performs various transformations on the input and output abundances and covariance matrix, and returns the final abundances after photodisintegration, and their associated covariance matrix.
The light element abundances in ACROPOLIS are taken as absolute abundances, as opposed to the treatment in GAMBIT and AlterBBN, where they are treated as ratios over the hydrogen abundance, $r_i = Y_i/Y_H$. Hence before feeding the abundances to ACROPOLIS, they are modified accordingly. Furthermore, the abundances computed by AlterBBN include the late-time decays of $^3$H and $^7$Be. As these happen after photodisintegration, and are in fact automatically added by ACROPOLIS, we revert the $^3$H and $^7$Be decays when converting to absolute abundances for sending to ACROPOLIS. Conversely, before returning the final abundances, they are converted back to ratios over $Y_H$.

ACROPOLIS takes the abundances as a fixed-size array ordered as {$n, p, D, ^3$H, $^3$He, $^4$He, $^6$Li, $^7$Li, $^7$Be} and, along with the input parameters, computes a transfer matrix $M$, which provides the abundances after photodisintegration as $Y_f = MY_i$. Besides the computation of the final abundances, this transfer matrix is used to propagate the covariance matrix as

$$ \Sigma_f = (\mathbf{J}_2\mathbf{M}\mathbf{J}_1)\Sigma_i(\mathbf{J}_2\mathbf{M}\mathbf{J}_1)^T, \qquad (A.1) $$

where $\mathbf{J}_1$ is the Jacobian of the transformation from ratios over $Y_H$ to absolute abundances, and $\mathbf{J}_2$ the converse. These Jacobians are different because the former includes the subtraction of late-time decays, whereas these are not included in the latter, but rather in the transfer matrix $M$.

Lastly, there are few internal options in ACROPOLIS that determine the precision and speed of the computations, among others. These can be provided as options in the YAML file, which can be seen below, including three configurations for default, fast and aggressive settings. The backend convenience function `set_input_params` ensures that the given options are passed to the relevant ACROPOLIS modules.

```yaml
# Options for ACROPOLIS
- capability: primordial_abundances
  function: primordial_abundances_decayingDM
  options:
    verbose: false
    NE_pd: 30 # default: 150, fast: 75, aggressive: 30
    NT_pd: 10 # default: 50, fast: 25, aggressive: 10
    eps: 1e-1 # default: 1e-3, fast: 1e-2, aggressive: 1e-1
```

## B Details of implementation of entropy injection between BBN and recombination

The evolution of the photon temperature including photon production through ALP decays $a \to \gamma\gamma$ is given by eq. (3.4). It is important to note that the functions $n_a(t)$ and $H(t)$ in eq. (3.4) both have an implicit dependence on the photon temperature $T_\gamma$ and neutrino temperature $T_\nu$. Both temperatures affect the energy density in radiation $\rho_{\text{rad}}$, see eq. (3.1), which in turn affects the expansion rate $H(t)$ and the scale factor $R(t)$ entering $n_a(t)$. Self-consistently integrating eq. (3.4) therefore requires taking the implicit temperature dependence into account. Doing so can, however, lead to numerical instabilities. In order to circumvent these complications, we take an alternative approach: instead of solving eq. (3.4) only once, while respecting the implicit dependence of $H(t)$ and $n_a(t)$ on $T_\gamma$, we solve this equation several times while neglecting the changes in $H(t)$ and $R(t)$ due to the changes in $T_\gamma$. After each iteration $H(t)$ and $R(t)$ are updated and the differential equation is solved again. This
means that in the ith iteration, the differential equation is given by

$$ \frac{dT_{\gamma}^{(i)}}{dt} = \frac{15}{4\pi^2} \frac{m_a n_a^{(i-1)}(t)}{\tau_a} \frac{1}{(T_{\gamma}^{(i)})^3} - H^{(i-1)}(t)T_{\gamma}^{(i)}, \quad (B.1) $$

where $H^{(i-1)}(t)$ and $n_a^{(i-1)}(t)$ are derived assuming the photon temperature of the previous iteration $T_{\gamma}^{(i-1)}$. This procedure is repeated until $T_{\gamma}$ has converged, i.e. $T_{\gamma}^{(i)}(t) = T_{\gamma}^{(i-1)}(t)$.

## B.1 Numerical implementation

Here we describe the numerical implementation of this procedure. We create a look-up table of the photon temperature $T_{\gamma}(t)$, the neutrino temperature $T_{\nu}(t)$, the Hubble rate $H(t)$, as well as the natural logarithm of the scale factor $\ln R(t)$, which enter $n_a(t)$ as given by eq. (3.5).<fnref n="16" />

When this table is initialised ($i=0$), it has logarithmic spacing in time between $t_{ini}$ and $t_{fin}$. Furthermore, we assume that there are no modifications to the photon temperature and neutrino temperature, meaning that initially the effective number of relativistic degrees of freedom $g_*(T)$ is assumed to be constant in the temperature range of interest. In this case $H = \dot{R}/R \sim T_{\gamma}^2$ and $T_{\gamma} \sim 1/R$ such that the initial Hubble rate $H^{(0)}(t)$ is given by

$$ H^{(0)}(t) = \frac{1}{2t}. \quad (B.2) $$

Using the right-hand side of eq. (2.4), the initial photon temperature $T_{\gamma}^{(0)}(t)$ is calculated from $H^{(0)}(t)$ and the initial neutrino temperature is given by assuming instant decoupling $T_{\nu}^{(0)} = (4/11)^{1/3} T_{\gamma}^{(0)}$. The last column of the table $\ln a^{(0)}(t)$ is initialised via

$$ \ln R^{(0)}(t) = \int_{t_{ini}}^{t} dt' H^{(0)}(t'). \quad (B.3) $$

After eq. (B.1) is solved the table has to be updated in order to be self-consistent. In particular, as $T_{\gamma}(t)$ is modified, $H(t)$ is modified as well. Moreover, it is important to note that the constant column in this table is not the time $t$ but the scale factor or its natural logarithm $\ln R(t)$. Hence, once the column containing $H$ is updated by using the respective values of $T_{\gamma}$ (and $T_{\nu}$) after solving the differential equation, the column containing the time coordinate $t$ has to be updated as well by using $H = d \ln R/dt$, such that $t$ becomes a function of $\ln R$ and is given by

$$ t(\ln R) = \int_{\ln R_0}^{\ln R} d \ln R' \frac{1}{H(\ln R')}. \quad (B.4) $$

Once $H$ and $t$ have been self-consistently updated, the differential equation can be solved again.

## C Details of the SN1987A limits from ALP decay

Here we give the full details of our SN1987A likelihood. While a description of the procedure and some of the relevant equations was already provided in the literature [101, 142, 143], we extend the discussion by adding more details on the reference frame transformations and Monte Carlo simulations.

<footnote n="16">¹⁶In the scenario studied in this paper, the neutrino temperature can be derived through eq. (3.6) such that there would be no need to include $T_{\nu}(t)$ in this table, as it can be directly derived through $\ln R(t)$. We include this, however, to allow for easy generalisation to scenarios in which the neutrino temperature is also modified.</footnote>
## C.1 Reference frame transformations

In figure 4, we show the ALP a emerging from SN1987A and decaying into two photons $\gamma_1$ and $\gamma_2$ with angles $\alpha_1$ and $\alpha_2$. For convenience, let the x-axis be parallel to the direction of the emergent ALP. We denote the ALP rest frame by “0” and the lab frame, i.e. the rest frame of the SMM satellite, by “L.” The Lorentz boost into the lab frame, for dimensionless velocity $\beta$ and $\gamma \equiv 1/\sqrt{1-\beta^2}$, is then given by

$$ \Lambda = \begin{pmatrix} \gamma & \beta\gamma & 0 & 0 \\ \beta\gamma & \gamma & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{pmatrix}. \qquad (C.1) $$

Boosting the 4-momentum of an ALP from the rest frame to the lab frame gives

$$ p_{a,0}^{\mu} = \begin{pmatrix} m_a \\ 0 \\ 0 \\ 0 \end{pmatrix} \quad \mapsto \quad p_{a,L}^{\mu} = \Lambda^{\mu\nu} (p_{a,0})_{\nu} = \begin{pmatrix} \gamma m_a \\ \beta\gamma m_a \\ 0 \\ 0 \end{pmatrix}. \qquad (C.2) $$

The two decay photons are emitted back-to-back in the ALP rest frame, which can be used to boost them into the lab frame:

$$ p_{1/2,0}^{\mu} = \frac{m_a}{2} \begin{pmatrix} 1 \\ \pm \cos \phi \sin \theta \\ \pm \sin \phi \sin \theta \\ \pm \cos \theta \end{pmatrix} \quad \mapsto \quad p_{1/2,L}^{\mu} = \frac{\gamma m_a}{2} \begin{pmatrix} 1 \pm \beta \cos \phi \sin \theta \\ \beta \pm \cos \phi \sin \theta \\ \pm \frac{1}{\gamma} \sin \phi \sin \theta \\ \pm \frac{1}{\gamma} \cos \theta \end{pmatrix}. \qquad (C.3) $$

Consequently, the two photon energies $E_1$ and $E_2$ in the lab frame are

$$ E_{1/2,L} = \frac{\gamma m_a}{2} (1 \pm \beta \cos \phi \sin \theta) = \frac{E_a}{2} (1 \pm \beta \cos \phi \sin \theta), \qquad (C.4) $$

while the photon angles $\alpha_1$ and $\alpha_2$ in the lab frame are defined geometrically as the angle between the 3-momenta and the x-direction $e_x$ such that

$$ \cos \alpha_{1/2,L} = \frac{e_x \cdot p_{1/2}}{|e_x| |p_{1/2}|} = \frac{E_a \beta \pm \cos \phi \sin \theta}{2 |p_{1/2}|}. \qquad (C.5) $$

Calculating the angles in eq. (C.5) requires some algebra, in particular we need that

$$ \frac{2}{E_a} |\mathbf{p}_{1/2}| = |1 \pm \beta \cos \phi \sin \theta| = 1 \pm \beta \cos \phi \sin \theta, \qquad (C.6) $$

where the absolute value in the last line can be dropped because the absolute values of $\beta$, $\cos \phi$, and $\sin \theta$ are all smaller than unity. The angles are therefore given by

$$ \alpha_{1/2,L} = \arccos \left( \frac{\beta \pm \cos \phi \sin \theta}{1 \pm \beta \cos \phi \sin \theta} \right). \qquad (C.7) $$

Equation (C.7) reveals additional rotational symmetry in the problem because only the product $\cos \phi \sin \theta$ is important; one could thus e.g. set $\theta = \pi$ for simplicity.
## C.2 Monte Carlo procedure

For each $(m_a, g_{a\gamma})$ pair we simulate $10^7$ events. This numbers reasonable since the naïve photon fluence in the parameter region of interest is of order $10^6 \text{ cm}^{-2}$ while the $3\sigma$ fluence limit is about $2 \text{ cm}^{-2}$. For example, ALPs with $m_a = 1 \text{ MeV}$ and $g_{a\gamma} = 10^{-11} \text{ GeV}^{-1}$ are just outside of the excluded region; the corresponding acceptance fraction of the MC procedure is about $2 \times 10^{-5}$, or $203 \pm 14$ events out of $10^7$.

We first draw ALP energies from the spectrum in eq. (3.17) by numerically calculating the inverse of the (normalised) spectrum's cumulative distribution. Similarly, we draw the decay lengths $L_a$ from an exponential distribution, with mean value parameter given by eq. (3.18). For ALPs to yield detectable decays, they have to:

a) decay at least $10^{10} \text{ m}$ away from SN1987A (the orange shell around SN1987A in figure 4), as photons from decays closer to the SN will be absorbed or deflected by its turbulent environment, and

b) decay before reaching the Earth, as almost no photons will reach the detector beyond this point (cf. discussion in ref. [101]).

For photons that satisfy these conditions, we draw a random emission angle for one of the two photons in the ALP rest frame from uniform distributions: $\phi \sim U(0, 2\pi)$ and $\cos\theta \sim U(-1, 1)$. The angles $\phi$ and $\theta$ of the other photon are the negative of the angles of the first photon, as the two photons are emitted back-to-back. Then, for each of the two photons, we calculate the energies $E_{1/2,L}$ according to eq. (C.6) in the lab frame and check that these energies are within the detector energy window of 25-100 MeV.

Finally, for the photons also making this cut, we calculate the photon angles $\alpha_{1/2}$ in the lab frame according to eq. (C.7). This allows us to compute the total length of their path, $L_{1/2}$, to the detector and hence their arrival times $t_{1/2}$, which are [101, section II.C]

$$L_{\gamma,1/2} \equiv -L_a \cos(\alpha_{1/2}) \sqrt{d_{\text{SN}}^2 - L_a^2 \sin^2(\alpha_{1/2})}, \qquad (C.8)$$

$$t_{1/2} = L_{1/2} = L_a/\beta + L_{\gamma,1/2} - d_{\text{SN}}. \qquad (C.9)$$

Demanding that photons arrive earlier than 223 s after the arrival of the first neutrino is the last check in the Monte Carlo procedure.

## C.3 Results

We tabulate the calculated fluence on a logarithmic grid of masses from 1 keV to 1 GeV and ALP-photon couplings from $10^{-13} \text{ GeV}^{-1}$ to $10^{-3} \text{ GeV}^{-1}$. We use a grid spacing of 0.05 dex in both dimensions, leading to a grid of 24,000 points. We chose this grid in order to close the gap between SN1987A and Horizontal Branch star limits [e.g. 101]. The resulting table is provided in GAMBIT as DarkBit/data/SN1987A_DecayFluence.dat.

The results are shown in figure 4, where we compare the excluded region from our implementation with that of ref. [101]. Our contours are smoother as we perform full Monte Carlo simulations without the high-mass cut used in ref. [101]. Apart from this difference, we get good qualitative agreement with the previous bounds. However, note that we typically find that our calculated flux is a factor $\lesssim 2$ larger than ref. [101] — in line with what ref. [41] observe using the non-probabilistic integration methods from ref. [144]. After inspecting the code of ref. [101] and its output, we suspect that apart from the wide binning used in ref. [101] there might also be some “double counting” in the implementation of the factors considered in their eq. (14).
### C.4 Technical details

Our scripts, written in Python v3.7.5 and using the ` numpy`, ` scipy`, ` time`, and ` mpi4py` packages, are available at https://github.com/marie-lecroq/ALP-fluence-calculation. The main script ` functions.py ` contains the function ` expected_photon_fluence`, which returns the photon fluence (in units of $cm^{-2}$) given an ALP mass $m$ (in units of eV) and ALP-photon coupling $g$ (in units of GeV$^{-1}$). We also included the script ` run_analysis_mpi.py ` to run our algorithm in parallel using MPI.

### References

[1] R. Diamanti, L. Lopez-Honorez, O. Mena, S. Palomares-Ruiz and A.C. Vincent, *Constraining Dark Matter Late-Time Energy Injection: Decays and P-Wave Annihilations*, JCAP 02 (2014) 017 [arXiv:1308.2578] [INSPIRE].

[2] T.R. Slatyer and C.-L. Wu, *General Constraints on Dark Matter Decay from the Cosmic Microwave Background*, Phys. Rev. D 95 (2017) 023010 [arXiv:1610.06933] [INSPIRE].

[3] V. Poulin, J. Lesgourgues and P.D. Serpico, *Cosmological constraints on exotic injection of electromagnetic energy*, JCAP 03 (2017) 043 [arXiv:1610.10051] [INSPIRE].

[4] W. Hu and J. Silk, *Thermalization constraints and spectral distortions for massive unstable relic particles*, Phys. Rev. Lett. 70 (1993) 2661 [INSPIRE].

[5] M. Lucca, N. Schöneberg, D.C. Hooper, J. Lesgourgues and J. Chluba, *The synergy between CMB spectral distortions and anisotropies*, JCAP 02 (2020) 026 [arXiv:1910.04619] [INSPIRE].

[6] V. Poulin and P.D. Serpico, *Nonuniversal BBN bounds on electromagnetically decaying particles*, Phys. Rev. D 91 (2015) 103007 [arXiv:1503.04852] [INSPIRE].

[7] M. Hufnagel, K. Schmidt-Hoberg and S. Wild, *BBN constraints on MeV-scale dark sectors. Part II. Electromagnetic decays*, JCAP 11 (2018) 032 [arXiv:1808.09324] [INSPIRE].

[8] M. Kawasaki, K. Kohri, T. Moroi, K. Murai and H. Murayama, *Big-bang nucleosynthesis with sub-GeV massive decaying particles*, JCAP 12 (2020) 048 [arXiv:2006.14803] [INSPIRE].

[9] P.F. Depta, M. Hufnagel and K. Schmidt-Hoberg, *Updated BBN constraints on electromagnetic decays of MeV-scale particles*, JCAP 04 (2021) 011 [arXiv:2011.06519] [INSPIRE].

[10] X.-L. Chen and M. Kamionkowski, *Particle decays during the cosmic dark ages*, Phys. Rev. D 70 (2004) 043502 [astro-ph/0310473] [INSPIRE].

[11] J.E. Kim, *Light Pseudoscalars, Particle Physics and Cosmology*, Phys. Rept. 150 (1987) 1 [INSPIRE].

[12] J. Jaeckel and A. Ringwald, *The Low-Energy Frontier of Particle Physics*, Ann. Rev. Nucl. Part. Sci. 60 (2010) 405 [arXiv:1002.0329] [INSPIRE].

[13] M. Cicoli, M. Goodsell and A. Ringwald, *The type IIB string axiverse and its low-energy phenomenology*, JHEP 10 (2012) 146 [arXiv:1206.0819] [INSPIRE].

[14] S. Weinberg, *A New Light Boson?*, Phys. Rev. Lett. 40 (1978) 223 [INSPIRE].

[15] F. Wilczek, *Problem of Strong P and T Invariance in the Presence of Instantons*, Phys. Rev. Lett. 40 (1978) 279 [INSPIRE].

[16] I. Brivio et al., *ALPs Effective Field Theory and Collider Signatures*, Eur. Phys. J. C 77 (2017) 572 [arXiv:1701.05379] [INSPIRE].

[17] M. Bauer, M. Neubert, S. Renner, M. Schnubel and A. Thamm, *The Low-Energy Effective Theory of Axions and ALPs*, JHEP 04 (2021) 063 [arXiv:2012.12272] [INSPIRE].
[18] E. Masso and R. Toldra, *On a light spinless particle coupled to photons*, Phys. Rev. D **52** (1995) 1755 [hep-ph/9503293] [INSPIRE].

[19] E. Masso and R. Toldra, *New constraints on a light spinless particle coupled to photons*, Phys. Rev. D **55** (1997) 7967 [hep-ph/9702275] [INSPIRE].

[20] D. Cadamuro and J. Redondo, *Cosmological bounds on pseudo Nambu-Goldstone bosons*, JCAP **02** (2012) 032 [arXiv:1110.2895] [INSPIRE].

[21] M. Millea, L. Knox and B. Fields, *New Bounds for Axions and Axion-Like Particles with keV-GeV Masses*, Phys. Rev. D **92** (2015) 023010 [arXiv:1501.04097] [INSPIRE].

[22] P. Arias, D. Cadamuro, M. Goodsell, J. Jaeckel, J. Redondo and A. Ringwald, *WISPy Cold Dark Matter*, JCAP **06** (2012) 013 [arXiv:1201.5902] [INSPIRE].

[23] M. Meyer, D. Horns and M. Raue, *First lower limits on the photon-axion-like particle coupling from very high energy gamma-ray observations*, Phys. Rev. D **87** (2013) 035027 [arXiv:1302.1208] [INSPIRE].

[24] M. Giannotti, I. Irastorza, J. Redondo and A. Ringwald, *Cool WISPs for stellar cooling excesses*, JCAP **05** (2016) 057 [arXiv:1512.08108] [INSPIRE].

[25] M. Bauer, M. Neubert and A. Thamm, *Collider Probes of Axion-Like Particles*, JHEP **12** (2017) 044 [arXiv:1708.00443] [INSPIRE].

[26] M.J. Dolan, T. Ferber, C. Hearty, F. Kahlhoefer and K. Schmidt-Hoberg, *Revised constraints and Belle II sensitivity for visible and invisible axion-like particles*, JHEP **12** (2017) 094 [Erratum ibid. **03** (2021) 190] [arXiv:1709.00009] [INSPIRE].

[27] D.J.E. Marsh, *Axion Cosmology*, Phys. Rept. **643** (2016) 1 [arXiv:1510.07633] [INSPIRE].

[28] I.G. Irastorza and J. Redondo, *New experimental approaches in the search for axion-like particles*, Prog. Part. Nucl. Phys. **102** (2018) 89 [arXiv:1801.08127] [INSPIRE].

[29] C. O'Hare, *AxionLimits: Data, plots and code for constraints on axions, axion-like particles and dark photons*, Zenodo (2020), current version of the repository available on Github at https://github.com/cajohare/AxionLimits [doi:10.5281/zenodo.3932430].

[30] J. Preskill, M.B. Wise and F. Wilczek, *Cosmology of the Invisible Axion*, Phys. Lett. B **120** (1983) 127 [INSPIRE].

[31] L.F. Abbott and P. Sikivie, *A Cosmological Bound on the Invisible Axion*, Phys. Lett. B **120** (1983) 133 [INSPIRE].

[32] M. Dine and W. Fischler, *The Not So Harmless Axion*, Phys. Lett. B **120** (1983) 137 [INSPIRE].

[33] M.S. Turner, *Coherent Scalar Field Oscillations in an Expanding Universe*, Phys. Rev. D **28** (1983) 1243 [INSPIRE].

[34] M.S. Turner, *Cosmic and Local Mass Density of Invisible Axions*, Phys. Rev. D **33** (1986) 889 [INSPIRE].

[35] P.F. Depta, M. Hufnagel and K. Schmidt-Hoberg, *Robust cosmological constraints on axion-like particles*, JCAP **05** (2020) 009 [arXiv:2002.08370] [INSPIRE].

[36] M. Bolz, A. Brandenburg and W. Buchmüller, *Thermal production of gravitinos*, Nucl. Phys. B **606** (2001) 518 [Erratum ibid. **790** (2008) 336] [hep-ph/0012052] [INSPIRE].

[37] GAMBIT COSMOLOGY WORKGROUP collaboration, *CosmoBit: A GAMBIT module for computing cosmological observables and likelihoods*, JCAP **02** (2021) 022 [arXiv:2009.03286] [INSPIRE].

[38] GAMBIT collaboration, *GAMBIT: The Global and Modular Beyond-the-Standard-Model Inference Tool*, Eur. Phys. J. C **77** (2017) 784 [Addendum ibid. **78** (2018) 98] [arXiv:1705.07908] [INSPIRE].
[39] A. Kvellestad, P. Scott and M. White, GAMBIT and its Application in the Search for Physics Beyond the Standard Model, arXiv:1912.04079 [INSPIRE].

[40] A. Ayala, I. Domínguez, M. Giannotti, A. Mirizzi and O. Straniero, Revisiting the bound on axion-photon coupling from Globular Clusters, Phys. Rev. Lett. 113 (2014) 191302 [arXiv:1406.6053] [INSPIRE].

[41] A. Caputo, G. Raffelt and E. Vitagliano, Muonic boson limits: Supernova redux, Phys. Rev. D 105 (2022) 035022 [arXiv:2109.03244] [INSPIRE].

[42] A. Caputo, H.-T. Janka, G. Raffelt and E. Vitagliano, Low-Energy Supernovae Severely Constrain Radiative Particle Decays, Phys. Rev. Lett. 128 (2022) 221103 [arXiv:2201.09890] [INSPIRE].

[43] A. Boyarsky, A. Neronov, O. Ruchayskiy, M. Shaposhnikov and I. Tkachev, Where to find a dark matter sterile neutrino?, Phys. Rev. Lett. 97 (2006) 261302 [astro-ph/0603660] [INSPIRE].

[44] A. Boyarsky, O. Ruchayskiy and M. Shaposhnikov, The Role of sterile neutrinos in cosmology and astrophysics, Ann. Rev. Nucl. Part. Sci. 59 (2009) 191 [arXiv:0901.0011] [INSPIRE].

[45] J.W. Foster et al., Deep Search for Decaying Dark Matter with XMM-Newton Blank-Sky Observations, Phys. Rev. Lett. 127 (2021) 051101 [arXiv:2102.02207] [INSPIRE].

[46] GAMBIT Community, Supplementary Data: Cosmological constraints on decaying axion-like particles: a global analysis, doi:10.5281/zenodo.6573347 (2022).

[47] J. Jaeckel, J. Redondo and A. Ringwald, 3.55 keV hint for decaying axionlike particle dark matter, Phys. Rev. D 89 (2014) 103511 [arXiv:1402.7335] [INSPIRE].

[48] P.F. de Salas, M. Lattanzi, G. Mangano, G. Miele, S. Pastor and O. Pisanti, Bounds on very low reheating scenarios after Planck, Phys. Rev. D 92 (2015) 123534 [arXiv:1511.00672] [INSPIRE].

[49] T. Hasegawa, N. Hiroshima, K. Kohri, R.S.L. Hansen, T. Tram and S. Hannestad, MeV-scale reheating temperature and thermalization of oscillating neutrinos by radiative and hadronic decays of massive particles, JCAP 12 (2019) 012 [arXiv:1908.10189] [INSPIRE].

[50] D. Blas, J. Lesgourgues and T. Tram, The Cosmic Linear Anisotropy Solving System (CLASS) II: Approximation schemes, JCAP 07 (2011) 034 [arXiv:1104.2933] [INSPIRE].

[51] PLANCK collaboration, Planck 2018 results. VI. Cosmological parameters, Astron. Astrophys. 641 (2020) A6 [Erratum ibid. 652 (2021) C4] [arXiv:1807.06209] [INSPIRE].

[52] PLANCK collaboration, Planck 2018 results. V. CMB power spectra and likelihoods, Astron. Astrophys. 641 (2020) A5 [arXiv:1907.12875] [INSPIRE].

[53] P. Stöcker, M. Krämer, J. Lesgourgues and V. Poulin, Exotic energy injection with ExoCLASS: Application to the Higgs portal model and evaporating black holes, JCAP 03 (2018) 018 [arXiv:1801.01871] [INSPIRE].

[54] GAMBIT collaboration, Thermal WIMPs and the scale of new physics: global fits of Dirac dark matter effective field theories, Eur. Phys. J. C 81 (2021) 992 [arXiv:2106.02056] [INSPIRE].

[55] J.J. Bennett, G. Buldgen, M. Drewes and Y.Y.Y. Wong, Towards a precision calculation of the effective number of neutrinos $N_{\text{eff}}$ in the Standard Model I: the QED equation of state, JCAP 03 (2020) 003 [Addendum ibid. 03 (2021) A01] [arXiv:1911.04504] [INSPIRE].

[56] K. Akita and M. Yamaguchi, A precision calculation of relic neutrino decoupling, JCAP 08 (2020) 012 [arXiv:2005.07047] [INSPIRE].

[57] J. Froustey, C. Pitrou and M.C. Volpe, Neutrino decoupling including flavour oscillations and primordial nucleosynthesis, JCAP 12 (2020) 015 [arXiv:2008.01074] [INSPIRE].
[58] J.J. Bennett et al., Towards a precision calculation of $N_{\text{eff}}$ in the Standard Model II: Neutrino decoupling in the presence of flavour oscillations and finite-temperature QED, JCAP 04 (2021) 073 [arXiv:2012.02726] [INSPIRE].

[59] P.F. Depta, M. Hufnagel and K. Schmidt-Hoberg, ACROPOLIS: A generic framework for Photodisintegration of Light elementS, JCAP 03 (2021) 061 [arXiv:2011.06518] [INSPIRE].

[60] M. Kawasaki and T. Moroi, Electromagnetic cascade in the early universe and its application to the big bang nucleosynthesis, Astrophys. J. 452 (1995) 506 [astro-ph/9412055] [INSPIRE].

[61] PARTICLE DATA GROUP collaboration, Review of Particle Physics, PTEP 2020 (2020) 083C01 [INSPIRE].

[62] J. Geiss and G. Gloeckler, Isotopic Composition of H, HE and NE in the Protosolar Cloud, Space Sci. Rev. 106 (2003) 3.

[63] P.F. Depta, M. Hufnagel, K. Schmidt-Hoberg and S. Wild, BBN constraints on the annihilation of MeV-scale dark matter, JCAP 04 (2019) 029 [arXiv:1901.06944] [INSPIRE].

[64] A. Arbey, AlterBBN: A program for calculating the BBN abundances of the elements in alternative cosmologies, Comput. Phys. Commun. 183 (2012) 1822 [arXiv:1106.1363] [INSPIRE].

[65] A. Arbey, J. Auffinger, K.P. Hickerson and E.S. Jenssen, AlterBBN v2: A public code for calculating Big-Bang nucleosynthesis constraints in alternative cosmologies, Comput. Phys. Commun. 248 (2020) 106982 [arXiv:1806.11095] [INSPIRE].

[66] B.D. Fields, The primordial lithium problem, Ann. Rev. Nucl. Part. Sci. 61 (2011) 47 [arXiv:1203.3551] [INSPIRE].

[67] D.J. Fixsen, E.S. Cheng, J.M. Gales, J.C. Mather, R.A. Shafer and E.L. Wright, The Cosmic Microwave Background spectrum from the full COBE FIRAS data set, Astrophys. J. 473 (1996) 576 [astro-ph/9605054] [INSPIRE].

[68] D.J. Fixsen, The Temperature of the Cosmic Microwave Background, Astrophys. J. 707 (2009) 916 [arXiv:0911.1955] [INSPIRE].

[69] Y.B. Zeldovich and R.A. Sunyaev, The Interaction of Matter and Radiation in a Hot-Model Universe, Astrophys. Space Sci. 4 (1969) 301 [INSPIRE].

[70] R.A. Sunyaev and Y.B. Zeldovich, The Interaction of matter and radiation in the hot model of the universe, Astrophys. Space Sci. 7 (1970) 20 [INSPIRE].

[71] Y.B. Zeldovich, A.F. Illarionov and R.A. Sunyaev, The Effect of Energy Release on the Emission Spectrum in a Hot Universe, Zh. Eksp. Teor. Fiz. 68 (1972) 1217.

[72] A.F. Illarionov and R.A. Siuniaev, Comptonization, characteristic radiation spectra, and thermal balance of low-density plasma, Sov. Astron. 18 (1975) 413.

[73] L. Danese and G. de Zotti, Double Compton process and the spectrum of the microwave background, Astron. Astrophys. 107 (1982) 39.

[74] C. Burigana, L. Danese and G. de Zotti, Formation and evolution of early distortions of the microwave background spectrum — A numerical study, Astron. Astrophys. 246 (1991) 49 [INSPIRE].

[75] J. Chluba and R.A. Sunyaev, The evolution of CMB spectral distortions in the early Universe, Mon. Not. Roy. Astron. Soc. 419 (2012) 1294 [arXiv:1109.6552] [INSPIRE].

[76] J. Chluba et al., Spectral Distortions of the CMB as a Probe of Inflation, Recombination, Structure Formation and Particle Physics: Astro2020 Science White Paper, Bull. Am. Astron. Soc. 51 (2019) 184 [arXiv:1903.04218] [INSPIRE].

[77] T. Yanagida and M. Yoshimura, Resonant Axion-Photon Conversion in the Early Universe, Phys. Lett. B 202 (1988) 301 [INSPIRE].
[78] A. Mirizzi, J. Redondo and G. Sigl, Constraining resonant photon-axion conversions in the Early Universe, JCAP 08 (2009) 001 [arXiv:0905.4865] [INSPIRE].

[79] T. Higaki, K. Nakayama and F. Takahashi, Cosmological constraints on axionic dark radiation from axion-photon conversion in the early Universe, JCAP 09 (2013) 030 [arXiv:1306.6518] [INSPIRE].

[80] H. Tashiro, J. Silk and D.J.E. Marsh, Constraints on primordial magnetic fields from CMB distortions in the axiverse, Phys. Rev. D 88 (2013) 125024 [arXiv:1308.0314] [INSPIRE].

[81] D. Ejlli and A.D. Dolgov, CMB constraints on mass and coupling constant of light pseudoscalar particles, Phys. Rev. D 90 (2014) 063514 [arXiv:1312.3558] [INSPIRE].

[82] S. Mukherjee, R. Khatri and B.D. Wandelt, Polarized anisotropic spectral distortions of the CMB: Galactic and extragalactic constraints on photon-axion conversion, JCAP 04 (2018) 045 [arXiv:1801.09701] [INSPIRE].

[83] S. Mukherjee, R. Khatri and B.D. Wandelt, Constraints on non-resonant photon-axion conversion from the Planck satellite data, JCAP 06 (2019) 031 [arXiv:1811.11177] [INSPIRE].

[84] S. Mukherjee, D.N. Spergel, R. Khatri and B.D. Wandelt, A new probe of Axion-Like Particles: CMB polarization distortions due to cluster magnetic fields, JCAP 02 (2020) 032 [arXiv:1908.07534] [INSPIRE].

[85] D. Cadamuro, S. Hannestad, G. Raffelt and J. Redondo, Cosmological bounds on sub-MeV mass axions, JCAP 02 (2011) 003 [arXiv:1011.3694] [INSPIRE].

[86] A.S. Kompaneets, The Establishment of Thermal Equilibrium between Quanta and Electrons, JETP 4 (1957) 730.

[87] A.P. Lightman, Double Compton emission in radiation dominated thermal plasmas, Astrophys. J. 244 (1981) 392.

[88] J. Chluba et al., New horizons in cosmology with spectral distortions of the cosmic microwave background, Exper. Astron. 51 (2021) 1515 [arXiv:1909.01593] [INSPIRE].

[89] J. Chluba, Refined approximations for the distortion visibility function and μ-type spectral distortions, Mon. Not. Roy. Astron. Soc. 440 (2014) 2544 [arXiv:1312.6030] [INSPIRE].

[90] J.C. Mather et al., Measurement of the Cosmic Microwave Background spectrum by the COBE FIRAS instrument, Astrophys. J. 420 (1994) 439 [INSPIRE].

[91] D.J. Fixsen, E. Dwek, J.C. Mather, C.L. Bennett and R.A. Shafer, The Spectrum of the extragalactic far infrared background from the COBE FIRAS observations, Astrophys. J. 508 (1998) 123 [astro-ph/9803021] [INSPIRE].

[92] J. Chluba, Green's function of the cosmological thermalization problem, Mon. Not. Roy. Astron. Soc. 434 (2013) 352 [arXiv:1304.6120] [INSPIRE].

[93] J. Chluba, D. Nagai, S. Sazonov and K. Nelson, A fast and accurate method for computing the Sunyaev-Zeldovich signal of hot galaxy clusters, Mon. Not. Roy. Astron. Soc. 426 (2012) 510 [arXiv:1205.5778] [INSPIRE].

[94] J. Chluba, E.R. Switzer, D. Nagai and K. Nelson, Sunyaev-Zeldovich signal processing and temperature-velocity moment method for individual clusters, Mon. Not. Roy. Astron. Soc. 430 (2013) 3054 [arXiv:1211.3206] [INSPIRE].

[95] B. Audren, J. Lesgourgues, K. Benabed and S. Prunet, Conservative Constraints on Early Cosmology: an illustration of the Monte Python cosmological parameter inference code, JCAP 02 (2013) 001 [arXiv:1210.7183] [INSPIRE].

[96] T. Brinckmann and J. Lesgourgues, MontePython 3: boosted MCMC sampler and other features, Phys. Dark Univ. 24 (2019) 100260 [arXiv:1804.07261] [INSPIRE].
[97] A. Kogut et al., *The Primordial Inflation Explorer (PIXIE): A Nulling Polarimeter for Cosmic Microwave Background Observations*, JCAP 07 (2011) 025 [arXiv:1105.2044] [INSPIRE].

[98] J. Chluba and D. Jeong, *Teasing bits of information out of the CMB energy spectrum*, Mon. Not. Roy. Astron. Soc. 438 (2014) 2065 [arXiv:1306.5751] [INSPIRE].

[99] N. Schöneberg, M. Lucca and D.C. Hooper, *Constraining the inflationary potential with spectral distortions*, JCAP 03 (2021) 036 [arXiv:2010.07814] [INSPIRE].

[100] M.H. Abitbol, J. Chluba, J.C. Hill and B.R. Johnson, *Prospects for Measuring Cosmic Microwave Background Spectral Distortions in the Presence of Foregrounds*, Mon. Not. Roy. Astron. Soc. 471 (2017) 1126 [arXiv:1705.01534] [INSPIRE].

[101] J. Jaeckel, P.C. Malta and J. Redondo, *Decay photons from the axionlike particles burst of type-II supernovae*, Phys. Rev. D 98 (2018) 055032 [arXiv:1702.02964] [INSPIRE].

[102] E.L. Chupp, W.T. Vestrand and C. Reppin, *Experimental Limits on the Radiative Decay of SN1987A Neutrinos*, Phys. Rev. Lett. 62 (1989) 505 [INSPIRE].

[103] D.J. Forrest et al., *The gamma ray spectrometer for the Solar Maximum Mission*, Sol. Phys. 65 (1980) 15.

[104] A. Payez, C. Evoli, T. Fischer, M. Giannotti, A. Mirizzi and A. Ringwald, *Revisiting the SN1987A gamma-ray limit on ultralight axion-like particles*, JCAP 02 (2015) 006 [arXiv:1410.3747] [INSPIRE].

[105] K. Sato and H. Sato, *Higgs Meson Emission from a Star and a Constraint on Its Mass*, Prog. Theor. Phys. 54 (1975) 1564 [INSPIRE].

[106] G.G. Raffelt, *Astrophysical methods to constrain axions and other novel particle phenomena*, Phys. Rept. 198 (1990) 1 [INSPIRE].

[107] G.G. Raffelt, *Stars as Laboratories for Fundamental Physics: The Astrophysics of Neutrinos, Axions, and Other Weakly Interacting Particles*, University Of Chicago Press, Chicago & London (1996).

[108] A. Buzzoni, F.F. Pecci, R. Buonanno and C.E. Corsi, *Helium abundance in globular clusters — The R-method*, Astron. Astrophys. 128 (1983) 94.

[109] G.G. Raffelt, *Core Mass at the Helium Flash From Observations and a New Bound on Neutrino Electromagnetic Properties*, Astrophys. J. 365 (1990) 559 [INSPIRE].

[110] E. Aver, K.A. Olive and E.D. Skillman, *The effects of He I λ10830 on helium abundance determinations*, JCAP 07 (2015) 011 [arXiv:1503.08146] [INSPIRE].

[111] GAMBIT DARK MATTER WORKGROUP collaboration, *DarkBit: A GAMBIT module for computing dark matter observables and likelihoods*, Eur. Phys. J. C 77 (2017) 831 [arXiv:1705.07920] [INSPIRE].

[112] S. Hoof, F. Kahlhoefer, P. Scott, C. Weniger and M. White, *Axion global fits with Peccei-Quinn symmetry breaking before inflation using GAMBIT*, JHEP 03 (2019) 191 [Erratum ibid. 11 (2019) 099] [arXiv:1810.07192] [INSPIRE].

[113] P. Athron et al., *Global fits of axion-like particles to XENON1T and astrophysical data*, JHEP 05 (2021) 159 [arXiv:2007.05517] [INSPIRE].

[114] P. Carenza, O. Straniero, B. Döbrich, M. Giannotti, G. Lucente and A. Mirizzi, *Constraints on the coupling with photons of heavy axion-like-particles from Globular Clusters*, Phys. Lett. B 809 (2020) 135709 [arXiv:2004.08399] [INSPIRE].

[115] BOSS collaboration, *The clustering of galaxies in the completed SDSS-III Baryon Oscillation Spectroscopic Survey: cosmological analysis of the DR12 galaxy sample*, Mon. Not. Roy. Astron. Soc. 470 (2017) 2617 [arXiv:1607.03155] [INSPIRE].
[116] PAN-STARRS1 collaboration, *The Complete Light-curve Sample of Spectroscopically Confirmed SNe Ia from Pan-STARRS1 and Cosmological Constraints from the Combined Pantheon Sample*, Astrophys. J. 859 (2018) 101 [arXiv:1710.00845] [INSPIRE].

[117] GAMBIT collaboration, *Comparison of statistical sampling methods with ScannerBit, the GAMBIT scanning module*, Eur. Phys. J. C 77 (2017) 761 [arXiv:1705.07959] [INSPIRE].

[118] S. Algeri, J. Aalbers, K. Dundas Morå and J. Conrad, *Searching for new phenomena with profile likelihood ratio tests*, Nature Rev. Phys. 2 (2020) 245 [arXiv:1911.10237] [INSPIRE].

[119] CMB-S4 collaboration, *CMB-S4 Science Book, First Edition*, arXiv:1610.02743 [INSPIRE].

[120] K. Abazajian et al., CMB-S4 Science Case, Reference Design, and Project Plan, arXiv:1907.04473 [INSPIRE].

[121] W.J. Handley, M.P. Hobson and A.N. Lasenby, *polychord: next-generation nested sampling*, Mon. Not. Roy. Astron. Soc. 453 (2015) 4385 [arXiv:1506.00171] [INSPIRE].

[122] R. Trotta, *Bayes in the sky: Bayesian inference and model selection in cosmology*, Contemp. Phys. 49 (2008) 71 [arXiv:0803.4089] [INSPIRE].

[123] J. Skilling, *Nested sampling for general bayesian computation*, Bayesian Anal. 1 (2006) 833.

[124] G. Ashton et al., *Nested sampling for physical scientists*, Nat. Rev. Meth. Primers 2 (2022) 39 [arXiv:2205.15570] [INSPIRE].

[125] L.T. Hergt, W.J. Handley, M.P. Hobson and A.N. Lasenby, *Bayesian evidence for the tensor-to-scalar ratio r and neutrino masses m$_{ν}$: Effects of uniform vs logarithmic priors*, Phys. Rev. D 103 (2021) 123511 [arXiv:2102.11511] [INSPIRE].

[126] W. Handley, *anesthetic: nested sampling visualisation*, J. Open Source Softw. 4 (2019) 1414 [arXiv:1905.04768] [INSPIRE].

[127] D. Baumann, D. Green and B. Wallisch, *New Target for Cosmic Axion Searches*, Phys. Rev. Lett. 117 (2016) 171301 [arXiv:1604.08614] [INSPIRE].

[128] F. D'Eramo, R.Z. Ferreira, A. Notari and J.L. Bernal, *Hot Axions and the $H_0$ tension*, JCAP 11 (2018) 014 [arXiv:1808.07430] [INSPIRE].

[129] D. Ghosh and D. Sachdeva, *Constraints on Axion-Lepton coupling from Big Bang Nucleosynthesis*, JCAP 10 (2020) 060 [arXiv:2007.01873] [INSPIRE].

[130] D. Green, Y. Guo and B. Wallisch, *Cosmological implications of axion-matter couplings*, JCAP 02 (2022) 019 [arXiv:2109.12088] [INSPIRE].

[131] G. Lucente, P. Carenza, T. Fischer, M. Giannotti and A. Mirizzi, *Heavy axion-like particles and core-collapse supernovae: constraints and impact on the explosion mechanism*, JCAP 12 (2020) 008 [arXiv:2008.04918] [INSPIRE].

[132] F. Calore, P. Carenza, M. Giannotti, J. Jaeckel and A. Mirizzi, *Bounds on axionlike particles from the diffuse supernova flux*, Phys. Rev. D 102 (2020) 123005 [arXiv:2008.11741] [INSPIRE].

[133] A. Goudelis, M. Pospelov and J. Pradler, *Light Particle Solution to the Cosmic Lithium Problem*, Phys. Rev. Lett. 116 (2016) 211303 [arXiv:1510.08858] [INSPIRE].

[134] P. Scott, *Pippi — painless parsing, post-processing and plotting of posterior and likelihood samples*, Eur. Phys. J. Plus 127 (2012) 138 [arXiv:1206.2245] [INSPIRE].

[135] GAMBIT COSMOLOGY WORKGROUP collaboration, *Strengthening the bound on the mass of the lightest neutrino with terrestrial and cosmological experiments*, Phys. Rev. D 103 (2021) 123508 [arXiv:2009.03287] [INSPIRE].
[136] BICEP2 and KECK ARRAY collaborations, BICEP2/Keck Array x: Constraints on Primordial Gravitational Waves using Planck, WMAP, and New BICEP2/Keck Observations through the 2015 Season, Phys. Rev. Lett. 121 (2018) 221301 [arXiv:1810.05216] [INSPIRE].

[137] GAMBIT collaboration, ColliderBit: a GAMBIT module for the calculation of high-energy collider observables and likelihoods, Eur. Phys. J. C 77 (2017) 795 [arXiv:1705.07919] [INSPIRE].

[138] GAMBIT FLAVOUR WORKGROUP collaboration, FlavBit: A GAMBIT module for computing flavour observables and likelihoods, Eur. Phys. J. C 77 (2017) 786 [arXiv:1705.07933] [INSPIRE].

[139] M. Chrzaszcz, M. Drewes, T.E. Gonzalo, J. Harz, S. Krishnamurthy and C. Weniger, A frequentist analysis of three right-handed neutrinos with GAMBIT, Eur. Phys. J. C 80 (2020) 569 [arXiv:1908.02302] [INSPIRE].

[140] GAMBIT MODELS WORKGROUP collaboration, SpecBit, DecayBit and PrecisionBit: GAMBIT modules for computing mass spectra, particle decay rates and precision observables, Eur. Phys. J. C 78 (2018) 22 [arXiv:1705.07936] [INSPIRE].

[141] S. Bloor et al., The GAMBIT Universal Model Machine: from Lagrangians to likelihoods, Eur. Phys. J. C 81 (2021) 1103 [arXiv:2107.00030] [INSPIRE].

[142] M. Giannotti, L.D. Duffy and R. Nita, New constraints for heavy axion-like particles from supernovae, JCAP 01 (2011) 015 [arXiv:1009.5714] [INSPIRE].

[143] A. Caputo et al., Axionlike Particles from Hypernovae, Phys. Rev. Lett. 127 (2021) 181102 [arXiv:2104.05727] [INSPIRE].

[144] L. Oberauer, C. Hagner, G. Raffelt and E. Rieger, Supernova bounds on neutrino radiative decays, Astropart. Phys. 1 (1993) 377 [INSPIRE].

[145] GAMBIT collaboration, GAMBIT: The Global and Modular Beyond-the-Standard-Model Inference Tool. Addendum for GAMBIT 1.1: Mathematica backends, SUSYHD interface and updated likelihoods, Eur. Phys. J. C 78 (2018) 98.