Axion global fits with Peccei-Quinn symmetry breaking before inflation using GAMBIT

We present global fits of cosmologically stable axion-like particle and QCD axion models in the mass range 0.1 neV to 10 eV. We focus on the case where the Peccei-Quinn symmetry is broken before the end of inflation, such that the initial value of the axion field can be considered to be homogeneous throughout the visible Universe. We include detailed likelihood functions from light-shining-through-wall experiments, haloscopes, helioscopes, the axion relic density, horizontal branch stars, supernova 1987A, white dwarf cooling, and gamma-ray observations. We carry out both frequentist and Bayesian analyses, with and without the inclusion of white dwarf cooling. We explore the degree of fine-tuning present in different models and identify parameter regions where it is possible for QCD axion models to account for both the dark matter in the Universe and the cooling hints, comparing them to specific DFSZ- and KSVZ-type models. We find the most credible parameter regions, allowing us to set (prior-dependent) upper and lower bounds on the axion mass. Our analysis also suggests that QCD axions in this scenario most probably make up a non-negligible but sub-dominant component of the dark matter in the Universe.


Introduction
QCD axions [1][2][3][4] and axion-like particles (ALPs) are perhaps among the most intriguing classes of hypothetical particles. They have been extensively studied in the literature and experiments are expected to probe the relevant parameter space for many axion models in the near future (see Ref. [5] for an up-to-date review on axion searches). Because QCD axions can behave as cold dark matter (DM) [6][7][8][9], the particularly interesting regions of the parameter space are where they contribute significantly to the DM density of the Universe. If so, they would solve at least two outstanding problems in physics at the same time (the other being the Strong CP problem [1,2]).
There are now many complementary searches for axions 1 underway, and many new search strategies have emerged in recent years (see Ref. [5]). Axions could also explain some apparent anomalies in the cooling of white dwarf stars [10][11][12][13][14][15][16], or the transparency of the Universe to gamma rays [17][18][19][20][21][22]. It is therefore crucial to combine all available results in order to extract the maximum information from the data. In doing so, we can learn more about the parameter space of different models, help guide the planning of future searches towards the most promising search areas and -if axions do indeed exist -find them in the correlated signals of several experiments and determine their properties.
Some of these goals can be satisfactorily achieved by over-plotting exclusion limits from several experiments. However, such simple exclusion plots generally have some shortcomings. They make assumptions about the relative importance of interactions with other matter, and about important inputs such as solar physics, the local density of DM and theoretical uncertainties. Incompatibilities between the assumptions can undermine the validity of a simple combination. They also do not allow for a quantitative model comparison between different theories that include axions.
A solution is to perform a consistent global statistical analysis of all available constraints, accounting for the leading theoretical and experimental uncertainties. Implementing such an analysis requires careful design, along with compromises relating to the availability of data and the total computational runtime required. An example of such a computational framework for BSM physics is GAMBIT [23]. Full details of GAMBIT's features can be found in the relevant publications [23][24][25][26][27][28] and physics analyses [29][30][31][32][33][34]. The most relevant features for this work are GAMBIT's modularity, and the options that it offers for carrying out both Bayesian and frequentist analyses. GAMBIT's structure allows easy integration of new components such as models, theory calculations, likelihoods and sampling algorithms. GAMBIT contains a variety of advanced samplers for both Bayesian and frequentist analyses, which are particularly useful for including nuisance parameters and assessing fine-tuning.
In the following section, we review some aspects of axion physics, including specific issues that should be accounted for in global fits. In Sec. 3, we turn to axion models, their corresponding effective field theories and the family tree of axion models available in GAMBIT. Sec. 4 describes our observables and likelihood functions, including their implementation, incorporated experimental data, potential caveats and restrictions. Results and discussion of the first global scans of axion models can be found in Sec. 5, and Sec. 6 summarises our results.
The axion routines developed for this paper are available in the DarkBit [25] module of GAMBIT 1.3.0, available at https://gambit.hepforge.org under the 3-clause BSD license. 2 Likelihood and posterior samples from this study can be downloaded from Zenodo [35].

Axion physics
Here we present a brief overview of axion physics, highlighting caveats of our current implementation and pointing out opportunities for future extensions. More details can be found in the literature [5,[36][37][38][39][40].

The QCD axion
The Peccei-Quinn (PQ) mechanism [1,2] is a proposed solution to the Strong CP problem that gives rise to a pseudo-scalar particle: the QCD axion [3,4]. The QCD axion is an excellent DM candidate [6][7][8][9], as it can account for the entire cosmological abundance of DM via the vacuum misalignment or realignment mechanism (Sec. 2.3.1). The original QCD axion model also inspired further archetypical models, such as the KSVZ [41,42] and DFSZ [43,44] axion models, which we will introduce in Secs 3.2.1 and 3.2.2.
The symmetries of the Standard Model (SM) Lagrangian permit a term of the form where G a µν is the gluon field strength, a is the SU (3) gauge index and α S is the strong coupling constant. The angle θ QCD ∈ [−π, π] is an unknown parameter and G µν,a = µνκλ G a κλ /2. A contribution is also generated by chiral transformations due to the chiral anomaly, which replaces θ QCD by an effective angle, (2.2) where Y d and Y u are the down-and up-type Yukawa matrices, respectively [45,Sec. 29.5].
The G G-term is anti-symmetric under the discrete parity (P ) and charge-parity (CP ) transformations. Due to the presence of this term, one would naïvely expect the strong interaction to show some CP -violating effects -especially because weak interactions are known to violate P maximally and CP mildly. There is no CP -violation if and only if θ eff vanishes (modulo the periodicity). Experiments trying to measure the CP-violating electric dipole moment of the neutron (nEDM) can place limits on the value of θ eff . Within 40-50% uncertainty, the dipole moment induced by the G G-term is given by [46][47][48] |d n | ≈ 2.4 × 10 −16 e cm |θ eff | . (2. 3) The current limit on the nEDM is |d n | < 3.6 × 10 −26 e cm at 95% confidence level [49], resulting in |θ eff | ∼ < 10 −10 . This observation poses a fine-tuning issue in the SM, commonly referred to as the Strong CP problem. The mechanism suggested by Peccei and Quinn [1,2] solves this problem by introducing a new global, axial U (1) symmetry. This so-called PQ symmetry is spontaneously broken by the vacuum expectation value v of a complex scalar field. The resulting Nabu-Goldstone boson is the QCD axion, a pseudo-scalar field denoted by a(x). It adds another contribution to θ eff of the form N a(x)/v, where the non-zero integer N is the colour anomaly of the PQ symmetry. The associated shift symmetry can then be used to cancel the θ eff term by driving θ eff + N a(x)/v to zero. In fact, it can be shown by the Vafa-Witten theorem [50,51] that the axion dynamically and asymptotically relaxes θ eff to the CP -conserving minimum. The Strong CP problem is therefore effectively solved by promoting θ eff to a dynamical degree of freedom.
The continuous shift symmetry of the PQ U (1) phase forbids any mass terms for the axion at the Lagrangian level. In the presence of an a G G term in the Lagrangian, however, this symmetry is broken after the QCD phase transition due to topologically non-trivial fluctuations of the gluon fields, leaving only a discrete shift symmetry of size 2πN v.
The resulting effective potential can be written as 3 where m a is the temperature-dependent axion mass and f a ≡ v/N . The axion is therefore practically massless until the time of the QCD phase transition, where it picks up a small mass that can be calculated using the chiral Lagrangian formalism.
The zero-temperature mass of the QCD axion was initially calculated by relating v to the scale of weak interactions [3]. However, this model was ruled out and the parameter space for f a was opened up in order to make the axion "invisible" -in the sense that it evaded experimental constraints at the time. Neglecting suppressed quark mass ratios, the QCD axion mass at zero temperature is given by [3] m a,0 where z d = m u /m d is the up-and down-quark mass ratio, and f π 0 and m π 0 are the decay constant and mass of the pion respectively. Using next-to-leading order chiral perturbation theory, a recent study found the QCD axion mass at zero temperature to be [53] m a,0 = 5.70 (7) µeV 10 12 GeV f a . (2.6) At temperatures exceeding the QCD scale, the axion becomes increasingly light as the shift symmetry is restored. Numerical estimates of the temperature dependence can be obtained directly from recent lattice QCD simulations [55,56]. These results are well-approximated at higher temperatures by a power law m a (T ) ∝ T −β/2 for some β > 0, and by a constant axion mass below some transition temperature T χ . This also agrees with the behaviour predicted from the finite-temperature dilute instanton gas approximation [57].
As we solve all relevant equations numerically in this paper, it would be straightforward to include the full lattice QCD results and calculate the axion mass at every given temperature. However, including the statistical uncertainties 4 is most easily achieved by having a parametrised form of the temperature dependence of the QCD axion mass. This is commonly done by introducing two parameters, β and T χ , which we fit to lattice QCD results. More details about the implementation can be found in Sec. 3.2.

Axion-like particles
Apart from the original QCD axion, pseudo-Nambu-Goldstone bosons with fundamental shift symmetries appear in many other contexts. These are usually referred to as ALPs (axion-like particles), and are appealing for theorists and model builders because of the shift symmetry, and the ability of their field values to undergo relaxation via the realignment mechanism in the same manner as axions. They do not necessarily solve the Strong CP problem. However, under the right circumstances, they can be cold DM candidates, or play an important role in solving other physics puzzles. ALPs typically arise from the breaking of a U (1) symmetry at some scale f a and generate a mass from explicit breaking of the residual symmetry at scale Λ. As a result, they can be fairly light, with masses of the order m a ∼ Λ 2 /f a , and have suppressed couplings to the SM. Due to the lack of a direct relation between the two scales involved (f a and m a ), they occupy a larger parameter space than QCD axions. More details on the theory and phenomenology of these particles can be found in the literature [5,36,58].

Axion creation mechanisms
QCD axions (with roughly µeV to meV masses) and many ALP models are expected to be lighter than other hypothetical particles, such as WIMPs, which typically have GeV to TeV-scale masses [59]. For thermal DM, this would be problematic, as DM would not be sufficiently cold today to reproduce the observed large-scale structure of the Universe. However, although a small population of axions is produced thermally, the relic abundance is typically dominated by non-thermal mechanisms, namely the realignment mechanism and topological defects. In this work, we focus on the realignment mechanism, which allows axions to be both ultra-cold and very light at the same time.

Realignment mechanism
The equation of motion for a homogeneous QCD axion or ALP field θ(t) = a(t)/f a with potential V (θ) in a Friedmann-Robertson-Walker-Lemaître universe reads [60, app. B.12] θ + 3H(t)θ + 1 where H ≡ȧ/a is the Hubble parameter. For the canonical axion potential, The general form of this equation does not possess analytic solutions. In the early Universe m a H and the system is an overdamped oscillator. We can integrate the differential equation with boundary conditions θ(t i ) = θ i andθ(t i ) = 0, where the value θ i of the axion field is called the initial misalignment angle.
At later times, around the time when m a ∼ H, the system becomes critically damped and the field starts to oscillate. In the regime of m a H, the axion field oscillations are adiabatic while their amplitude continuously decreases due to Hubble damping. In other words, at some time t , the axion field evolution is guaranteed to eventually enter an epoch where |θ(t)| 1 for all t > t . In this harmonic limit, (2.9) is a damped harmonic oscillator, and the field evolution is very well described by the WKB approximation [38] θ(t) ∝ m a (t) a 3 (t) where we can match the solution at time t with the phase δ . This limiting form of the axion scalar field oscillations can be used to demonstrate that they behave like cold DM at sufficiently late times by averaging (2.10) over an oscillation period to find the effective behaviour of the energy density, ε, such that [38] ε ∝ m a a −3 . (2.11) Once the mass becomes constant, this is the scaling behaviour of cold matter. We can rephrase this statement for non-relativistic axions to see that the (averaged) comoving axion number density is conserved: (2. 12) In fact, there is an adiabatic invariant in the harmonic limit of sin(θ) θ, which has been studied in the literature for QCD axions and ALPs [7,9,61]. However, if the initial misalignment angle θ i ∼ O(1), then the adiabatic invariant is not sufficient to describe the system, and anharmonic effects arising from the full potential must be taken into account. These have been estimated and calculated by several authors [9,[62][63][64][65].
The energy density in axions today can hence be parametrised by an overall transfer function F , with the properties that F is bounded from below, and becomes constant in the harmonic limit |θ i | 1 [9,63,64]. This allows us to conveniently write the axion energy density in the Universe today as where we explicitly denote the dependence of the transfer function F on the values of the axion parameters. For a given set of these parameters, the transfer function -and hence the axion relic density -could be determined once and used in subsequent calculations. For sampling over the parameter space, it is necessary to repeat this process for each set of parameters (or to tabulate the results on a sufficiently dense grid ahead of the scan).
We have thus far only considered the homogeneous field equation (2.9). It is not obvious that this is sufficient for determining the energy density in axions today. Let us consider the two possible scenarios in a cosmology with inflation: either PQ symmetry breaks after inflation, or it breaks before inflation ends.
In the first case, the Universe consists of a large number of causally disconnected "bubbles" where the initial misalignment angles take random values from a uniform distribution on the interval −π to π [9]. As a consequence, the initial misalignment angle θ i effectively becomes a function of space. Nonetheless, the resulting overall energy density in axions is fixed, because it can be calculated as the average over all the misalignment angles in all the patches. As θ i is constant within the patches, we may use (2.13) to calculate the result: In the second case, where the PQ symmetry breaks before the end of inflation, one causally connected patch gets blown up to at least the size of the observable Universe. The axion field is therefore homogeneous, with a random initial value in the interval (−π, π]. The stochastic nature of the initial misalignment angle gives rise to a physically motivated prior probability for θ i in a Bayesian analysis, which will be discussed in Sec. 5.
A major issue with this scenario is perturbations from inflation, which will affect the initial misalignment angle by some amount that depends on the energy scale of inflation [see e.g. 62]. It has been argued that this scenario is therefore finely-tuned to a degree that can be considered "worse" than the Strong CP problem itself [66,67]. Since the necessary additional implementation of inflationary models is beyond the scope of this work, we neglect the issue of field fluctuations during inflation.
Some authors have parametrised the topological defect contribution as Ω a = Ω realign where Ω realign a is the contribution from realignment. Including α as a nuisance parameter in our work would remove much of the predictability of the model, and would assume similar scaling of the energy density contributions. As we consider PQ symmetry to break before inflation, we assume that contributions from topological defects are diluted away during inflation [87].

Thermal creation
Axions are also thermally produced in the early Universe [69,88]. The resulting abundance is however rather dependent on the additional new field content associated with the axion.
For the QCD axion, the Boltzmann equations for the least model-dependent processes give useful estimates. The contributions of low-temperature processes such as π + π π + a have been computed [69,[88][89][90], including all combinations of pionic states π. Thermal axion production during reheating has also been studied [91][92][93]. As hadronisation would not yet have taken place by the time of reheating, these calculations consider processes like axion-gluon interactions, i.e. g + g g + a. Recent calculations also considered modeldependent processes with quarks q, i.e. q +q g + a and q + g q + a, which can be dominant in the early Universe [94].
A significant abundance of thermal axions modifies the effective number of relativistic species, which can be used to set limits on the axion mass [95][96][97][98][99][100][101][102][103]. Generally speaking, thermally-produced QCD axions with a mass of more than O(eV) are hot DM with a relic abundance comparable to that of neutrinos and photons. Those axion models are therefore excluded, as they would exceed the observational bounds on the fraction of DM that can be hot [e.g. 40,69,87]. However, such bounds have some dependence on the choice of cosmological datasets [102], and the limits can be relaxed if the Universe had a non-standard thermal history [104].

Axion models and theoretical uncertainties
Let us now discuss the various axion models that we consider in this work with a focus on the interactions between axions and other particles. We already saw that axions acquire a mass, so that they interact gravitationally. QCD axions also interact via the strong interaction. For generic ALPs, however, it is not clear a priori what the PQ charges of SM particles are, nor whether the axion interacts with a given particle at tree level or only at higher order in perturbation theory.  In order to study axion phenomenology and to identify useful observables, axions have been studied in an effective field theory (EFT) framework [105][106][107], which can be adjusted to a given energy scale and scenario. One can then apply exclusion limits on the effective interactions to specific axion models that establish a relation between the effective couplings and the fundamental axion parameters; examples of such models include the so-called KSVZ model [41,42] and the DFSZ model [43,44].
One important consequence of the fundamental shift symmetry is that axions can only directly interact with matter via derivative couplings. Below the electroweak scale, the effective Lagrangian for axion-SM interactions is where a is the axion field, F and G are the duals of electromagnetic and strong field strengths and g af f and g aγγ are effective coupling constants of mass dimension −1.
In principle, f runs over all SM fermions with mass m f , but the couplings most relevant for axion searches are those to electrons and nucleons, the latter arising from matching the EFT to chiral perturbation theory [107]. 5 Also note that the Lagrangian (3.1) is flavour-diagonal, which is not necessarily the case in all models. While we do not consider the possibility of flavour non-diagonal interactions in this paper, the resulting models can be phenomenologically interesting, rather predictive, and within reach of existing and future experimental searches [e.g. [109][110][111]. 5 The same matching also gives rise to interactions between axions and mesons, which can in particular induce flavour-changing rare decays. While of phenomenological relevance for heavier ALPs [108], these processes do not lead to relevant constraints on the parameter space considered in the present work.
We now turn to the specific axion models implemented in GAMBIT (see Fig. 1). We describe these models in the following subsections, while we leave the discussion of the parameter ranges and prior distributions for Sec. 5. Note that the models presented here are only a subset of the many general and more specific ALP models studied in the literature; detailed overviews can be found in review articles [5,40].
Models in GAMBIT are defined as collections of named parameters. All relevant observables for a model must be computable from these parameters. Models can have relationships to other models, allowing parameter combinations in one model to be translated to equivalent combinations in another model or to an alternative parameterisation of the same model. This is achieved by adding new "children" to the family tree of a more general model and defining a translation between the two models. This translation may, e.g., fix the values of some of the "parent" parameters or compute them as functions of (possibly new) parameters in the child model. This ensures flexibility in the choice of the independent model parameters to work with in any given theory calculation, allowing the most convenient definitions to always be used for any calculation. More details on the general implementation of models in GAMBIT can be found in Ref. [23].

General ALP model
We define a new family of GAMBIT models for axions and ALPs. On top of the family tree is the GeneralALP model, with seven parameters. This model describes an effective Lagrangian that is a subset of (3.1), In this study we do not include couplings of axions to nucleons. However, future versions of GAMBIT may include additional interactions such as these, subject to the availability of interesting observables and constraints. We employ the rescaled field value θ = a/f a , assume the canonical cosine potential (2.8) for all types of axions, 6 and take A summary of all the model parameters can be found in Table 1.

QCD axion models
The QCDAxion model is a child of the GeneralALP model. It is inspired by the original QCD axion models, which solve the Strong CP problem. The scale of the explicit breaking of the shift symmetry by instanton-like effects is therefore linked to the QCD scale. This connection can be exploited to uniquely determine the parameters m a,0 , β and T χ . However, there are uncertainties from theory, experiment and lattice QCD simulations that should be taken into account. We treat these as nuisance parameters; Table 2 provides an overview.  Exploiting the connection to QCD, we can replace the parameter m a,0 of the GeneralALP model with the energy scale Λ χ , defined via the topological susceptibility at zero temperature, To determine Λ χ , we use first principle calculations of the zero-temperature axion mass [53]. We include these results via a Gaussian likelihood whereΛ χ = 75.5 MeV and σ Λχ = 0.5 MeV are the most likely values for Λ χ and its uncertainty from Ref. [53]. 7 This is shown in the left panel of Fig. 2

Figure 2:
Nuisance likelihoods for the scale Λ χ (left, from direct theory calculations [53]) and β and T χ (right, from lattice QCD [55]). Note that β is correlated with T χ when the full lattice QCD results are taken into account.
The two parameters β and T χ can then be constrained by using, e.g., the full lattice QCD results from Table S.7 in the supplementary material of Ref. [55]. We construct a likelihood function by performing a chi-squared fit to the N = 20 data points, where X = log 10 [χ(T | β, T χ )/χ(T = 0)] is the logarithm of the normalised topological susceptibility, andX i and σ 2 X i are its value and uncertainty for the ith data point. 8 In Fig. 2, we show the two-dimensional profile likelihood for β and T χ . These parameters show a clear correlation. This is expected, as a higher transition temperature implies a steeper slope in the temperature-dependent branch of the axion mass (corresponding to larger β) in order to maintain a good fit to the shape of m a (T ), as determined by lattice QCD. We provide the best-fit values in Table 2.
We note that the fit to our functional form for the QCD axion mass (3.3) captures the temperature dependence established by lattice QCD well everywhere except in the region around T = T χ . This is because (3.3) is not smooth there and the overall fit is poor (χ 2 = 55.7 for 18 d.o.f., which corresponds to a p-value of about 10 −5 ). However, the disagreement stems only from a narrow temperature range, and has no impact on any of our results. Excluding the only data point in that region improves the fit to an acceptable level (χ 2 = 21.6 for 17 degrees of freedom, which corresponds to a p-value of about 0.2). tice QCD methods. 8 As the uncertainty quoted in Ref. [55] also includes an uncertainty on the value of Λχ, we divide the topological susceptibility χ = f 2 a m 2 a by their best-fit value for χ(T = 0), and remove this uncertainty by assuming simple error propagation.
We also replace the axion-electron coupling with the model-dependent factor C aee g aee = m e f a C aee , (3.7) and the axion-photon coupling with the model-dependent ratio of the electromagnetic and colour anomalies E/N where C aγγ is the model-independent contribution from axion-pion mixing. We use the value obtained in Ref. [53] for C aγγ and include it as a simple nuisance likelihood, with X = C aγγ andX and σ X again the most likely value and its uncertainty (see Table 2). Finally, we want to emphasise some subtle considerations about the coupling strengths of the QCD axion model. The possible ratios E/N are rational numbers arising from group theoretical considerations. In this study, we sample over E/N as if it were a continuous parameter, as the possible rational numbers that it can take on are quite densely spaced along the real line, at least over the range of values that we consider. The sometimes so-called "classical axion window" considers a canonical, small and somewhat arbitrarilydefined range of couplings for the prototypical axion models [105,[113][114][115], arising from only quite a small range in E/N . It has recently been pointed out that the range of choices can, indeed, be extended to include more possibilities [116][117][118][119]. To assess the whole range of various axion models, we use the minimum and maximum values for E/N from Table IV in Ref. [117], so that E/N ∈ [−4/3, 524/3]. These values arise from a systematic study of DFSZ-and KSVZ-type axion models, where the additional heavy quarks in KSVZ-type models have cosmologically safe lifetimes and do not introduce Landau poles below the Planck scale. Furthermore, in DFSZ-type models, the number of Higgs doublets may go up to the maximum of nine.
Note that (3.8) implies the possibility of having g aγγ < 0 within the valid range for E/N . Note however that all the likelihood functions that we use in this paper only depend on the absolute value of g aγγ . We therefore plot only |g aγγ | in our results, even though we do scan over parameter values that lead to negative couplings in the range −3.25 ∼ < E/N − C aγγ ≤ 0.

KSVZ models
The archetypical axion model is the KSVZ model [41,42], where the SM is supplemented by one or more electrically neutral, heavy quarks. In our implementation, the KSVZAxion is a child model of the QCDAxion, where the anomaly ratio, E/N , is a free parameter. In these models, axions have no tree-level interactions with fermions. However, there is a loop-induced coupling to electrons due to the axion-photon interaction, which -in the absence of a leading order contribution -must be taken into account [106]: whereΛ ∼ 1 GeV is the QCD confinement scale. Several previous works have used this expression withΛ = 1 GeV [16,120], even though this quantity is not uniquely defined. We too assumeΛ = 1 GeV, relying on the fact that any deviations enter only logarithmically. Although the anomaly ratio in the original KSVZ model was E/N = 0 [41,42], other assignments are possible. As in Ref. [16], we will consider four different KSVZAxion models: E/N = 0, 2/3, 5/3, and 8/3.

DFSZ models
In contrast to the KSVZ model, DFSZ models are obtained by adding an additional Higgs doublet to the SM [43,44]. This results in direct axion-electron interactions. One can define two manifestations of this model, often called DFSZAxion-I and DFSZAxion-II. The couplings in these two models are given by where tan(β ) is the ratio of the two Higgs vacuum expectation values [106]. Perturbative unitary requires 0.28 < tan(β ) < 140. It is therefore convenient to replace the parameter C aee in the QCDAxion model by tan(β ).

ALP models with constant mass
As a template model for ALPs that arise as pseudo-Nabu-Goldstone bosons but do not have a temperature-dependent mass, we define the ConstantMassALP model. This is mainly for convenience in studies where we want to parametrically explore the coupling space whilst keeping the inverse dependence on f a , but are not interested in a temperature-dependent mass.
For ConstantMassALP models, T χ is irrelevant because β = 0, reducing the total number of parameters to five. Similar to QCDAxion models, we replace the ALP mass with a pseudo-Nambu-Goldstone scale Λ and introduce dimensionless coupling constants C aγγ and C aee , consistent with the other models: Table 3 gives an overview of the observables and likelihood functions that we use in this paper. In what follows, we give details of the experimental data, computational methods and likelihood implementations that they employ.

Light-shining-through-wall experiments
Light-shining-through-wall (LSW) experiments shine laser light through a magnetic field onto an opaque material, and attempt to detect it on the other side. A photon may convert into an axion within the magnetic field before the material, pass through it as an axion, Table 3: Overview of the likelihood functions that we employ in this paper (in the order they are discussed).

Likelihood/Observable Comments References
QCD nuisance parameters

Stellar cooling hints (optional likelihood)
White dwarf cooling hints only considered in Sec. 5.5 [11,12,14,15]; Sec. 4.5.4 and convert back into a photon in the magnetic field on the other side [137][138][139]. Examples of experiments based on this technique are ALPS [121,140] and OSQAR [141,142]. The predicted number of photons on the opposite side of the wall to the laser source is where tot is the detector efficiency, P is the laser power, ω γ the laser energy, and t obs the observation period. P (γ → a → γ) is the probability of a photon converting into an axion and back, in an appropriately aligned magnetic field of length L and strength B. It is given by P (γ → a → γ) = P 2 (γ → a) with [121] where M 2 ≡ m 2 a /2 + ω 2 γ (n − 1), sinc(x) ≡ sin(x)/x, m a is the axion mass, ω γ is the photon energy, and n is the refractive index of the medium in the experimental setup (n = 1 for vacuum).
Our LSW likelihood is based on the final results of the ALPS-I experiment, using data for both evacuated and gas-filled magnets [121]. The ALPS Collaboration took data in "frames" of t f = 1 h each, binning physical pixels of the detector into 3 × 3-pixel blocks. ALPs refer to these simply as "pixels" of area 42 µm × 42 µm. The collaboration searched their data frames for cosmic rays and signatures of other systematics, over a wide region around the single pixel where most of the laser light would fall in the absence of a wall, referred to as the "signal pixel". We show the 90%, 95% and 99% C.L.s compared to the envelope of their strongest vacuum or gas results. The difference between our likelihood and the limit published by ALPS arises mostly due to the fact that we combine both likelihoods, rather than taking the envelope.
ALPS adjusted the raw ADU values (electron counts) of the signal pixel in an attempt to account for the average background in surrounding pixels. These reduced ADU (ADUred) values are obtained by subtracting the average ADU values across all pixels in the region surrounding the signal pixel, from the ADU value of the signal pixel. Doing this for every frame where the laser was on ("signal frames") and off ("background frames"), ALPS constructed histograms of ADUred values for both signal and background. By fitting Gaussian functions to these histograms, they were able to estimateb and σ b , the expected value and standard deviation of ADUred for the background, as well asô and σ o , the equivalent quantities for signal frames. An example and more details can be found in Fig. 2 of Ref. [121] (see also Ref. [140]). From (4.1), the expected signal s from photon-axion-photon production per frame of ALPS-I data (with B = 4.98 T, L = 4.2 m, and ω γ = 2.33 eV) is Apart from fluctuations in the experimental performance tot , which amounts to an uncertainty of ∆ tot / tot ≈ 6%, the experimental parameters are known to sufficient precision to be fixed to their reference values. We incorporate the uncertainty on tot into the estimate of the signal prediction, σ s = s∆ tot / tot . As a result of the ADU reduction procedure, the measured signal and background estimators are non-integer numbers. We test the signal-plus-background hypothesis with a Gaussian likelihood: We add together the log-likelihoods for data sets i = 1 to 3, where two data sets (five and six frames compared to 122 and 47 background frames, respectively) consist of vacuum data and one data set (eight frames compared to 155 background frames) consists of argon gas data. 9 We show the resulting exclusion limits in Fig. 3 and compare with the envelope of the strongest vacuum or argon gas limits in Ref. [121]. The differences between the published results and our implementation are due to fact that we combine the likelihoods instead of just adopting the more constraining of the two, and also because the authors of Ref. [121] used the Feldman-Cousins method [143], assuming Gaussianity and physical signals s > 0. Considering the vacuum data alone, the method of Ref. [143] gives a slightly stronger limit than our log-likelihood ratio method (6.5 × 10 −8 GeV −1 at 90% confidence limit (C.L.) in the low mass limit, as compared to 6.9 × 10 −8 GeV −1 in our implementation). By combining the vacuum and argon likelihoods however, our final limit is somewhat stronger: g aγγ < 5.8 × 10 −8 GeV −1 at 90% C.L. at low masses.

Helioscopes
Axion helioscopes attempt to detect axions produced by interactions in the Sun by observing the solar disc with a "telescope" consisting of a long magnet contained in an opaque casing [144][145][146]. Solar-produced axions would pass through the casing, convert into photons in the field of the magnet, and be observed in a detector behind the magnet.
Multiple processes in the Sun can produce axions: resonant production in the oscillating electric field of the solar plasma, non-resonant production in solar magnetic fields, and emission from the interaction of electrons with photons, nuclei or one another. The predicted axion-induced photon flux at Earth therefore depends on the assumed solar model, and on the couplings of the axion to both photons and electrons [see e.g. 147].
The dominant process for axion-photon interactions is Primakoff production [148,149], where photons are resonantly converted into axions in the presence of an atomic nucleus. The rate at which axions of energy E can be produced in a plasma from photons of the same energy is given by [150,151] Here, the inverse screening length is given in the Debye-Hückel approximation by with n e denoting the electron number density, and n j and Z j representing the number density and charge, respectively, of the jth nucleus. Note that the number densities and temperature, and hence the conversion rate, vary with the radial position r. Using the expression for the axion-photon conversion probability (4.2), the differential photon flux at the detector is [122] dΦ with ρ a dimensionless radial co-ordinate in the Sun, and r a dimensionless radial co-ordinate on the solar disc on the sky. The quantity D is the (average) Sun-Earth distance, which we take to be one astronomical unit, and ω 2 pl (ρ) = 4πα EM n e (ρ)/m e is the plasma frequency calculated from the electron number density n e (ρ) and electron mass m e . The upper limit of outer integral, r s , controls how much of the inner part of the image of the solar disc on the detector is included in the analysis. This need not always be 1, as the signal-to-noise ratio can be maximised by introducing a cut-off r s < 1 [see e.g . 122].
The contribution to the solar axion flux from axion-electron interactions can be taken into account by including the additional interaction rate Γ e→a (E, ρ) in (4.8). However, these contributions are not so straightforward to calculate from first principles. This is due to narrow free-bound transition lines [152][153][154][155], axion bremsstrahlung [150,156,157] and Compton scattering [158][159][160]. To include these contributions in our signal prediction, we use tabulated data for the axion-electron spectrum provided in Ref. [147].
The spectrum in Ref. [147] was computed with the 2009 iteration of the AGSS09met solar model [161], which is based on photospheric abundances for non-refractory species and meteoritic abundances for refractory elements [162]. The AGSS09met model is thus the default solar model in GAMBIT, and the one that we use throughout the rest of this paper. We however make use of its latest iteration [163] in preference to the earlier version wherever possible, such as when computing the axion flux from axion-photon interactions. In the limit of small axion mass, the predicted flux from axion-photon interations deviates by no more than 4.4% between solar models, with the greatest difference ocurring between the GS98 [163,164] and most recent AGSS09met models [163].
Full details of our integration routines for the axion-photon and axion-electron contributions, as well as the options available for the inclusion of solar models for axion physics in GAMBIT, can be found in Appendix A.

2007 CAST results
The first of our two CAST likelihoods is based on the CCD results published in 2007 [122] (CCD detector data from the 2004 vacuum run in Table 1 of Ref. [122]). The other data in Ref. [122] are less constraining; including them only improves the upper limit by 1% [122]. We therefore do not provide separate likelihoods for the other runs.
Although the 2007 CAST analysis was based on the solar model of Ref. [165], and a follow-up analysis [166] on axion-electron interactions was based on a different model [167],  [122], including the 2013 re-interpretation in terms of axion-electron couplings [166]. Left: Exclusion limits for the axion-photon coupling [122]. Right: Exclusion limits for the product of axion-photon and axion-electron coupling, assuming that axion-electron interactions dominate the axion production inside the Sun. We only make this assumption here to compare our implementation to Ref. [166].
here we use the AGSS09met model of Ref. [163]. For both the axion-photon and axionelectron interactions, we integrate the total flux over all 20 energy bins (from 0.8 to 6.8 keV), taking into account the observation time, effective area, and detector efficiency (see Appendix A). In this case r s ≈ 0.23 for the axion-photon contribution. As we use the interpolated spectrum for the axion-electron contribution (calculated for r s = 1), we can only rescale the resulting flux by an overall factor in order to estimate the flux inside r s = 0.23. We take this number to be 0.826, assuming the findings from the axion-photon interaction also apply in the axion-electron case [122]. Our implementation follows the CAST analysis [122,166], using a Poisson likelihood in each of the 20 energy bins where o i , s i and b i are respectively the observed number of photons, the number of expected signal photons, and the expected number of background photons based on observations away from the Sun, in the ith energy bin. In total, 26 photons were observed in the detector during data-taking compared to 30.9 expected background events. The resulting exclusion limits can be found in Fig. 4.

2017 CAST results
Our implementation of the latest CAST results [123] is analogous to that of the 2007 results, using the signal and expected background counts, exposure and detector efficiencies for the 2017 data. 10  [123]. Differences are mainly due to our simplified implementation of the likelihood function, which does not employ event-level information.
To calculate the signal predictions, we integrate the axion-photon and axion-electron fluxes over each of the 10 energy bins from 2 to 7 keV, and then scale the predictions by the effective exposure for each of the ten datasets in Ref. [123].
The exclusion limits presented in Ref. [123] are based on an unbinned likelihood. Here, we treat each energy bin in each of the ten datasets as a separate counting experiment, combining them into a binned Poisson likelihood (4.10) Here o are respectively the observed number of photons, the expected number of signal events and the expected number of background events in the ith energy bin of the jth experiment. In total, 226 photons were observed in the detector during data taking compared to 246.6 expected background events. Our slightly different choice of likelihood function to the original CAST analysis is significantly simpler, because it does not require event-level information -but it still reproduces the published exclusion limits rather well (see Fig. 5).

Haloscopes (cavity experiments)
Axion haloscopes are designed to detect DM axions by resonant axion-photon conversion inside a tunable cavity [144,145]. Microwave cavities are the most sensitive axion experiments in existence, but only cover a small mass range compared to other techniques. The ability of haloscopes to detect axions therefore directly depends on their cosmological abundance, and to a lesser extent, their velocities in the Galactic halo [168]. Here we consider only the case where axions are fully virialised within the halo.  [124,125], UF [126] and ADMX [127-132, 170, 171] experiments. Right: Magnified details of the latest ADMX results [132].
The power expected to be converted from axions to photons in a cavity is [144,145,169] where C is the form factor (a dimensionless integral over the E-and B-field configuration of the cavity), B 0 is magnetic field strength in the cavity, V is its volume, and Q and Q a are the quality factors of the cavity and axions, respectively. In this context, Q describes the ratio of stored vs dissipated energy of the cavity, while Q a is proportional to the axion velocity dispersion (just as Q effectively characterises the bandwidth of the cavity).
The signal prediction also depends on the local DM density in axions, which we obtain by rescaling the local DM density as Obtaining exclusion limits from cavity experiments is often quite involved, generally requiring simulation of the selection procedure of the detector [e.g . 128]. Without access to this information, we must approximate the underlying likelihood functions based on the publicly available limits and publications.
In the following, we describe our likelihoods for three different haloscope experiments. An overview of the resulting exclusion limits can be found in Fig. 6.

RBF results
The Rochester-Brookhaven-Fermi (RBF) collaboration performed a search for axions using several cavities [124,125]. Table I, Eq. (26) and Fig. 14 in Ref. [125] provide useful information for approximating their results.
By (4.11), the axion-induced power in a haloscope is proportional to ρ a, local g 2 aγγ , which we define as the "reduced signal" s ≡ ρ a, local g 2 aγγ . (4.13) The remaining factors from (4.11) are effectively constant across all frequency/mass bins and detectors. The signal is expected to occur in a single frequency bin i, which satisfies m a,0 ∈ [ω i , ω i+1 ). Using this definition, we adopt an ansatz for the likelihood function inspired by Eq. (26) of Ref. [125]: Here a i is a threshold parameter, b i effectively corresponds to an expected standard deviation of the reduced signal, and Θ is the Heaviside step function. The threshold values a i arise because RBF manually inspected all candidate frequencies in their data over a certain significance level. The two parameters are related as a i = N b i , where N is the number of standard deviations required for manual inspection of a candidate signal. Although Table I in Ref. [125] would allow us to determine b i , using the central frequency of the bin as well as the 95% C.L. on the coupling strength, the resulting bins are not quite identical to what is shown in Fig. 14 of the same paper. We therefore determine a i and b i for each frequency bin from the limits in Fig. 15 of that paper, assuming N = 4 in all cases (cf. Table I in the same paper). This leads to limits in 14 bins ranging from m a,0 = 4.4 to 10.1 µeV and from m a,0 = 11.2 to 16.2 µeV.

UF results
While the results from the University of Florida (UF) collaboration [126] could be implemented in the same way as the RBF experiment, the published data do not allow us to infer the threshold parameter a i for the one mass bin. However, Eq. (6) in Ref. [126] quotes the "noise power fluctuation", which we use as a standard deviation σ P ≈ 2.86 × 10 −22 W for the expected power P . We obtain the expected power for each axion model using the information provided in Ref. [126], and check that the quoted limit is comparable to the expected sensitivity (which we obtain by assuming that the observed data is equal to the background expectation). The corresponding likelihood function for the single bin from m a,0 = 5.4 to 5.9 µeV and signal s from (4.13) is given by (4.15)

ADMX results 1998 -2009
The procedure used by the ADMX Collaboration for setting limits in the absence of a detection is highly customised for the experiment [128]. Unfortunately, the necessary information for fully implementing their numerous results [127][128][129][130][131]170] is not available. Similar to the RBF likelihood, we therefore use the reduced signal (4.13) and the following ansatz for the likelihood: Here s i is the known limit in the ith frequency/mass bin, and a and b are free parameters that have to be determined by a fit to published exclusion curves. We do this using the exclusion limits at three confidence levels published in Ref. [127], and by assuming that the shape of the likelihood function is representative of the shape in all other ADMX bins over the range from m a,0 = 1.90 to 3.65 µeV. Doing so results in a = 0.0131 and b = 0.455.

ADMX results 2018
In their recent publication, the ADMX collaboration increased the sensitivity of their setup, which is now able to rule out some DFSZ-type models [132] in the range 2.66 ≤ m a,0 ≤ 2.81 µeV. We approximate the likelihood of this result using the 90% C.L. limits in Fig. 4 of Ref. [132], for the Maxwell-Boltzmann velocity distribution (consistent with the model of the halo velocity distribution that we assume for analysing the results of other searches).
Because the experimental setup changed compared to the previous runs, we do not employ the shape parameters from Sec. 4.3.3 for the 2018 dataset. Instead, as with the UF experiment we assume that ADMX saw no signal events, approximating the 2018 likelihood as Here, the effective standard deviation is given by σ 2 eff = σ 2 stat (m a,0 ) + σ 2 sys (s). We infer the statistical contribution σ stat by setting the log-likelihood at the observed values of the limits in Fig. 4 of Ref. [132] to that corresponding to a 90% C.L. for one degree of freedom. We read off σ stat at 194 different masses, and interpolate between them linearly for intermediate mass values (ignoring the narrow region from 2.7302 to 2.7307 µeV where the ADMX limits do not apply, due to radio interference [132]). We add the systematic uncertainty of 13% of the signal prediction (quoted in Ref. [132]) in quadrature with the statistical uncertainty.

Dark matter relic density
The realignment mechanism gives an axion contribution to the observed cold DM relic density. We solve (2.9) numerically as a function of temperature to obtain Ω a , the fraction of the critical energy density in axions today; details can be found in Appendix B. The resulting energy densities for the QCDAxion are shown in Fig. 7, similar to the presentation in other works [172, e.g.]. We reiterate that in this paper, we do not consider other contributions to the relic density than vacuum realignment.
The relic density of DM is very well constrained by the most recent Planck analysis [133]. We employ a Gaussian likelihood with the central value and standard deviation from [133] (Ω DM h 2 = 0.1188, σ exp = 0.0010), combining the experimental uncertainty in quadrature with a further 5% theory uncertainty, 11 σ theo = 0.05 Ω a h 2 , to give Band of m a,0 -θ i combinations (from Diver) to get the correct DM density (including β, T χ , and Λ χ as nuisance parameters). We show the results of Ref. [55] for comparison and also include the hypothetical case of a "temperature-independent QCD axion" (with Λ χ as a nuisance parameter, but β = 0 and T χ irrelevant) as an example.
GAMBIT offers two options for this likelihood: a detection or an upper limit. These allow us to demand either that axions are responsible for all DM, or only a fraction. For the upper limit, we simply set ln (L) = 0 for Ω a < Ω DM in (4.18). Except where we state otherwise, we show results based on the upper limit option.

Astrophysical probes
Astrophysical systems can provide significant additional constraints on axions, especially the axion-electron and axion-neutron couplings, which are not well constrained by other probes. Due to their weak interactions with matter, axions can efficiently transport energy across large distances in free space, or through stellar matter, thereby influencing stellar structure and evolution. Intriguingly, a number of astrophysical systems appear to exhibit an unexplained mechanism of energy transport, which might be due to ALPs: white dwarfs display apparently anomalous cooling rates (Sec. 4.5.4) or deviations in the shape of their luminosity function [182,183] and highly-energetic gamma rays seem to experience significantly less attenuation through intergalactic space than might be expected [17][18][19][20][21][22]. Unfortunately, the systematic uncertainties associated with these potential hints of new physics are quite difficult to quantify. Nonetheless, if the observations and associated theoretical uncertainties turn out to be robust, ALPs can indeed explain the observed deviations from expectation.
The probability of photon-axion conversion in a domain of size filled with a suitably aligned magnetic field B eff and plasma with electron number density n e is given by [184,191] where ω pl and E crit are the plasma frequency and critical energy, respectively: The quantity E crit describes the energy scale at which photons will efficiently convert into axions in the extragalactic magnetic field domains. In the absence of dust or any other photon absorber, after traversing N such domains of size , the remaining fraction of photons is [185,186,191,198] Equations (4.19) and (4.21) reveal that we do not expect to see any effect due to axions for E E crit , because p 1 0. For a given photon energy E, this happens for large axion masses m a and small axion-photon couplings g aγγ .
For E E crit , on the other hand, conversion is very efficient, but the observed spectrum would simply decrease by a constant factor over the entire energy range. In this case it is also not possible to test the axion hypothesis: the expected spectral normalisation of any source is not well constrained and it therefore has to be a free fitting parameter in the analysis.
It is only possible to constrain models where the critical energy lies within the spectral window of the instrument, such that one end of the spectrum is suppressed, but the other is not [192,193]. Limits from the distortion of gamma-ray spectra are therefore strongest at axion masses that lead to threshold energies similar to the photon energies observed by the experiment. This explains the characteristic shape of the limits (in particular why this method is not sensitive to axion masses below a certain value; see Fig. 8).
The H.E.S.S. Collaboration applied this technique to the spectrum of the active galactic nucleus PKS 2155-304, using data obtained with their Cherenkov telescope array [134]. Unfortunately, their signal prediction requires Monte Carlo simulation of magnetic field realisations, which are no longer available. 12 We therefore approximate their likelihood function for the galactic cluster magnetic field using Figs 6 and 7 of Ref. [134], based on a scheme that we describe in detail in Appendix C. The main idea is to use common interpolation methods inside the published exclusion contours, and to extrapolate to likelihood values outside the known contours using a method that mimics the shape of the known iso-likelihood contours and preserves the mathematical   properties of the likelihood. Our approximation procedure is of course somewhat arbitrary. The advantages of our scheme are that it exactly reproduces the known exclusion curves by construction, and that the general likelihood function is well-behaved. The obvious downsides are that we can neither guarantee that the likelihood in the outermost and innermost regions are completely accurate, nor can we extend it to larger couplings than the values shown in Fig. 7 of Ref. [134].

Supernova 1987A
Supernovae are excellent particle laboratories. Supernova 1987A (SN1987A) provides significant constraints on axion parameters, based on the neutrino burst duration [199][200][201][202] and axion-photon interaction in magnetic fields external to the supernova [203,204]. Our likelihood for SN1987A is based on the results from Ref. [135]. The authors of that study derived limits based on the absence of a coincident gamma-ray burst from SN1987A, which should have been observed by instruments on board of the Solar Maximum Mission [205] if axions were produced in the explosion and converted to gamma rays in the Galactic magnetic field.
The gamma-ray spectrum of photons with energy E expected at Earth per unit time from axions produced in SN1987A is where d is the distance to the supernova, P (a → γ) is the conversion probability (4.19), N is the number of axions created in the supernova and dṄ /dE is the axion spectrum at the source, as predicted from a supernova model. To obtain the measured photon fluence at Earth (gamma rays per unit area during the observation), F a→γ , equation (4.22) has to be integrated over the energy range and time duration of the observation, while modelling the transport of the axions through the Galactic magnetic field. 13 From Fig. 6 in Ref. [135], we can see that the photon fluence at Earth for a given axion-photon coupling becomes constant below a certain mass scale m * and rapidly decreases for bigger masses. This is not surprising, given that most axion experiments lose sensitivity at large masses due to the loss of coherence in axion-photon conversion. We therefore make the following ansatz for the fluence, where F ≈ 0.57 × 10 −12 cm −2 is the fluence for small axion masses at the reference value of g aγγ ≈ 5.3 × 10 −10 GeV −1 ; we obtain this value by integrating Eq. (2.11) in Ref. [135]. We determine the best-fit values for m * and the exponent b from the higher-mass region (m a,0 > 6.0 × 10 −10 eV) via a least-squares fit to the fluence contour in Fig. 6 of Ref. [135], givingm * ≈ 5.43 × 10 −10 eV andb ≈ 4.02. The value ofb that we obtain can be qualitatively understood by examining the axion-photon conversion probability given in (4.19), which can be written as and, since the oscillatory part is washed out by the turbulent magnetic fields, the conversion probability is effectively suppressed by a factor of m 4 a .
The likelihood for s = F a→γ in the absence of a photon signal with background fluctuations σ 2 F is then given by (4.24) Fig. 9 shows how the limits obtained from our approximation compare to the original reference.

Horizontal Branch stars and R parameter
Weakly-interacting particles can influence stellar evolution by providing an additional energy loss mechanism, cooling stars over the course of their evolution [208][209][210]. The so-called R parameter, R = N HB /N RGB , is the ratio between the number of Horizontal Branch (HB) stars, N HB , and upper Red Giant Branch (RGB) stars, N RGB , in Galactic globular clusters. Its value depends on the relative time that stars spend on each branch, which is sensitive to the details of stellar evolution and cooling. Axions are expected to be produced in the cores of both types of stars, but would remove heat more efficiently from the cores of HB stars, reducing the time that they spend on the HB and leading to a reduction in R.
Based on a weighted average of a selection of cluster count observations [211], the observed value is R obs = 1.39 ± 0.03 [136]. The dependence of the predicted value of the R parameter on the properties of axions can be approximated as [13,136,[212][213][214] where Y is the Helium abundance, x aγγ ≡ g aγγ /10 −10 GeV −1 and x aee ≡ g aee /10 −13 . The equation above is valid only if axions are sufficiently light compared to the typical temperatures of the stellar interior, which are T ∼ 10 8 K ≈ 10 keV [210], i.e. much higher than the axion masses we consider. Our R parameter likelihood is then simply where σ pred and σ obs are the uncertainties of the predicted and observed values. The authors of Ref.
[13] adopted a Helium abundance of Y = 0.255 ± 0.002 [215], leading to a predicted value from standard (axion-free) stellar evolution calculations of R theo = 1.45 ± 0.01, almost 2σ higher than the observed value. We adopt the updated value for such low-metallicity environments of Y = 0.2515 ± 0.0017 [216], leading to R theo = 1.42 ± 0.01, entirely consistent with the observed R parameter. The effect of the uncertainty of Y on R pred can be estimated according to (4.25), as σ pred = 7.33 × 0.0017 ≈ 0.012. A comparison to the exclusion curves in Ref. [13] can be found in Fig. 10.

White Dwarf cooling hints
White dwarfs (WDs) are a particularly interesting environment in which to study axionelectron interactions, due to their electron-degenerate cores [10,[217][218][219]. Current observations can be interpreted as indicating a need for additional cooling in WDs compared   [13]. Note that our adopted likelihood leads to limits only, whereas the results of Ref. [13] indicate an almost 2σ preference for signal (hence the presence of 1σ upper and lower limits in the results that we plot from Ref. [13]). The difference is due to the updated He abundance that we employ here.
to standard models. The coupling necessary to explain the cooling with axions has been estimated to be g aee ∼ O 10 −13 [13]. A more recent analysis also considered these cooling hints in a global fitting framework [16]. Whilst the systematics of such analyses are still a matter of debate, and alternative explanations (whether involving BSM physics or not) are certainly still possible, it is intriguing to investigate the impact on axion global fits of including the WD cooling hints. Due to the speculative nature of theses hints, we will do so in a separate (alternative) analysis in Sec. 5.5, presented alongside our main global fits. WDs typically pulsate, allowing the oscillation of their radii and luminosity to be used to probe their internal structure via astroseismology. The periods, Π, of their pulsations decrease with time, with a rate X = dΠ/dt, which can be related to the energy loss in the system.
Refs [11,12,14,15] simulated the evolution of WDs with and without axions, predicting the period decrease dΠ/dt in each case. For our predictions of WD cooling rates, we interpolate these results and their stated uncertainties, using natural splines. The specific figures and objects from those papers that our implementation is based on are listed in Table 4. 14 Note that the plots in Refs [11,12,14,15] show the quantity m a,0 cos 2 (β ) which, for the DFSZAxion-I axion model they consider, is proportional to the more fundamental parameters C aee or g aee (cf. Sec. 3.2.2). For values of the axion-electron coupling larger than considered in the simulations (i.e. g aee /10 −13 > 5.6 or 8.4), we assign the likelihood corresponding to the largest simulated coupling. This is a conservative assumption, as Table 4: Overview of available cooling hints for WD variables of spectral types DA and DB. We list the couplings allowed at 2σ confidence (whether limits or intervals). The numerical values necessary for constructing the corresponding likelihoods were kindly provided by A. H. Córsico and T. Battich. In the case of R548, we have two more data points compared to what is shown in Fig. 1 of Ref. [12].

Object
Type   the disagreement between prediction and observation will in reality only worsen as the coupling increases further (until the WDs become opaque to axions -but this would occur at couplings well beyond what we consider).
For each WD listed in Table 4, we use a simple Gaussian likelihood function for the observed (X obs,i ) and theoretically-expected (X pred,i ) period decrease, such that our total WD cooling likelihood is where, again, the predictions and corresponding uncertainties σ pred,i are taken from the respective figures in the references listed in Table 4 and interpolated via natural splines.
The resulting individual and combined likelihoods can be found in Fig. 11.
We emphasise that the interpretation of WD cooling is subject to a number of assumptions and caveats. Statistical and systematic uncertainties associated with the inputs and algorithms of stellar models were considered by the authors of Refs [11,12,14,15], but a number of other potential issues remain. These include theoretical modelling of the transition from the main sequence to the WD phase, and the accuracy of the observed period decrease of PG 1351+489. Despite these problems, and in contrast to Ref. [16], we include the PG 1351+489 system in our discussion. The authors of Ref. [16] exclude this object due to its similarity to R548, and the uncertainties associated with R548 being "more conservative". However, the different estimated uncertainties in these two systems are in fact due to a real physical effect, namely the difference in the influence of trapped vs non-trapped oscillation modes in the two systems. While the latter might give rise to concerns regarding the understanding of different modes (cf. Table 4), we do not conclude that the arguments in favour of excluding PG 1351+489 are strong enough to do so.

Results and discussion
In this section, we present the central findings from our global fits of various axion models, identifying the most promising regions in parameter space and comparing the various models. We present frequentist and Bayesian results side-by-side, after discussing our choice of priors for the model parameters.

Sampling algorithms and settings
We use the differential evolution sampler Diver [28] to sample the composite likelihood function and T-Walk [28] to sample the posterior distributions. We employ MultiNest [224][225][226] primarily to compute Bayesian evidences.
We use the sampler settings established in an earlier study by the GAMBIT Collaboration as starting points [28]. For Diver, we generally use a population size (NP) of 2 × 10 4 and a tolerance (convthresh) of 10 −4 , and turn off the lambdajDE optimisation, preferring to use regular jDE for its slightly less aggressive optimisation. In addition to combining samples from various runs, where necessary to resolve fine-tuned regions we increase NP to 3 × 10 4 or 5 × 10 4 and/or reduce convthresh to 10 −5 . For T-Walk, we use the default settings for 340 or 544 MPI processes, until reaching a tolerance (sqrtR − 1) of 0.01 or 0.005. All initialisation YAML files that we use in this study are available on Zenodo [35]; the exact scans for which we use the different settings can be ascertained by examining the input files. Because we use MultiNest primarily for estimating Bayesian evidences, we set the sampling efficiency (efr) to 0.3, as recommend for this task [225], and use 2 × 10 4 live points (nlive) with a tolerance (tol) of 10 −4 .

General ALP models
Starting at the top of the model hierarchy (see Fig. 1), we first consider the GeneralALP model. This is a phenomenological model; parameter combinations in this model need not correspond to physical models, as their couplings do not depend on the inverse of f a . The main purpose of the GeneralALP is to provide a straightforward, universal connection to observables and to compare to results in the literature.
Parameter ranges. The parameter ranges and scales that we use for GeneralALP models are given in Table 5. Because the axion potential is periodic, all normalised field values are equivalent to a value in the interval (−π, π]. For g aγγ , g aee , and m a,0 , there is no obvious range to choose; we adopt parameter ranges encompassing values informed by previous studies and phenomenology. Recall that g aγγ could be negative but the likelihood functions in the present work only depend its absolute value. We therefore scan only over positive values of g aγγ , but label plots with |g aγγ |, to make it explicit that the results are equally valid for the corresponding negative values. For the local DM density ρ 0 , we adopt the same range as in earlier GAMBIT studies [29][30][31][32] (and implement the same likelihood function).
The appropriateness of different ranges on f a , β, and T χ depends on the fundamental properties of the symmetries and scales of the underlying ALP model. Although there are theoretical arguments for the existence of ALPs from e.g. string theory [38,227], they do not provide any quantitative guidance. 15 Apart from the likely case that our calculations become meaningless for f a ∼ > m Pl , we can only impose β > 0, as the axion mass should become smaller as the underlying symmetry is restored at higher temperatures. The ranges in Table 5 are therefore an attempt to include a variety of cases around the values known or preferred for the QCD axion.
We do not produce Bayesian results for the GeneralALP model, as no strong physical arguments exist for any particular choice of prior on most of its parameters. The only exception is the initial misalignment angle θ i , due to the causal structure of the early-Universe cosmology mentioned in Sec. 2.3.1. Frequentist results. We first scan the GeneralALP model assuming ALPs to be all of DM. The resulting limits on the axion-photon coupling (Fig. 12) are comparable to summary plots elsewhere in the literature [e.g. 115]. However, we would like to stress that unlike overplotted exclusion limits, the exclusion curve in Fig. 12 arises from a composite likelihood, and profiling takes into account uncertainties in the local DM density (the only relevant nuisance parameter here). The left panel of Fig. 13 shows that the joint constraints on the two coupling parameters are essentially dictated by the constraint on the R parameter (Fig. 10). As a consequence, the axion-electron equivalent of Fig. 12 would show that values of g aee ∼ > few × 10 −13 are excluded across the entire mass range. The right panel of Fig. 13 shows the possible combinations of f a and |θ i | that allow the GeneralALP to be all of DM. The extent and shape of this region is mostly due to the limited ranges of m a,0 , β, and T χ . The axion potential, and therefore the initial energy density in axions, is proportional to f a m a,0 . The observed DM abundance can only be achieved if f a m a,0 is large enough, because |θ i | < π. On the other hand, the axion starts to oscillate when H ∼ m a (Sec. 2.3.1). The associated temperature scale depends on m a,0 , β, and T χ and sets the amount by which the axion energy density is red-shifted up to the present day. To obtain the correct abundance in axions today while e.g. going to lower values of f a , the values of m a,0 , β, or θ i must be increased or the value of T χ decreased.
Within our selected parameter ranges, the profile likelihood does not identify preferred   Fig. 14 is a consequence of how the energy density scales right after the axion field begins to oscillate. For the allowed range of values for m a,0 , the corresponding band of possible f a -given a value for the initial misalignment angle -is also different.

QCD axions
QCD axions are the most well-studied type of axions to date. Unlike in previous studies, here we take into account the uncertainties due to nuisance parameters (see Sec. 3.2) in every part of the analysis. We also consistently scale the local DM density in axions according to their cosmological abundance. This affects the limits on the axion-photon interaction from haloscope experiments such as ADMX, as the detector signal in (4.13) is proportional to ρ a,local g 2 aγγ . The limits on g aγγ therefore scale with 1/ √ ρ a, local . We again cap the local axion abundance at 100% of the local DM abundance, and penalise models that predict too much DM via the Planck likelihood for Ω DM h 2 (Sec. 4.4).

Prior choices.
The priors that we apply to the QCDAxion parameters can be found in Table 6. This model imposes a number of relations between the phenomenological parameters of its parent GeneralALP model, which depend on nuisance parameters, i.e. quantities determined by simulations, theory or experiments, which are only known within an appreciable uncertainty. While we are generally not interested in inference on such parameters, their uncertainties can affect results for the actual parameters of interest. The additional nuisance parameters for QCDAxion models are C aγγ , Λ χ , β, and T χ . For C aγγ and Λ χ , the nuisance likelihood is given by a 1D Gaussian for each parameter, whereas the likelihood for β and T χ takes into account correlations between the two parameters (Sec. 3.2). We choose flat priors from about −5σ to +5σ around the respective central values for all four nuisance parameters.
The range of values that we choose for the anomaly ratio, E/N , is inspired by the selection criteria and range established in phenomenological studies of axion models (cf. Sec. 3.2). While the different preferred models presented in Ref. [117] form a discrete set, we assume that there is a continuous band of possible axion models, spanning a range from the lowest (E/N = −4/3) to the highest (E/N = 524/3) possible value of the anomaly ratio. Given that the number of possibilities grows very quickly if we allow for an arbitrary number of new heavy quarks in KSVZ-type models (where, however, E/N ≤ 170/3), it is not inconceivable that such a band exists. We treat each value inside the band as equally probable before contact with data, employing a flat prior for E/N . 16 Note that this necessarily encompasses negative values for g aγγ , as discussed in Sec. 3.1; whilst this does not impact our likelihoods (which depend only on |g aγγ |), it does imply an asymmetric effective prior on the two signs of g aγγ .
Assigning priors to f a and C aee is more difficult. For f a , we choose a range that corresponds to our region of interest in mass: from the largest masses allowed by bounds on hot DM, to highly fine-tuned regions with very small masses and f a somewhat below the Planck scale. The logarithmic prior reflects our ignorance about the scale of new physics, given that the ability of the original QCD axion to solve the Strong CP problem does not depend on the value of f a . We choose a generous range for C aee , taking a logarithmic prior around values of order unity, which may be considered the most natural value for C aee . Note that the lower end roughly corresponds to the minimum value that can be constrained by the R parameter for the highest QCDAxion masses we consider. Values any lower will be effectively indistinguishable.
Our choice of priors on ρ 0 and θ i follow the logic from the previous section on GeneralALP models. Because ρ 0 is rather well constrained by data, the choice of log or flat prior has little impact on the final results. For θ i , the causal structure of the early-Universe cosmology mentioned in Sec. 2.3.1 means that all values of the initial misalignment angle are equally likely, so a flat prior is most appropriate.
Frequentist results. First, let us consider statistical inference on the axion coupling strengths. There are essentially only upper limits on the axion couplings or the associated model parameters. We begin by focusing on the axion-photon interactions, as determined by the anomaly ratio E/N , shown in the upper row of Fig. 15. In the left panel we impose the relic density constraint as an upper limit, while in the right panel we demand that axions be all of DM. A notable difference between these two assumptions is that haloscopes (UF, RBF and ADMX) only provide strong limits in the latter case.
If axions are not required to be all of DM, the high-mass (low-f a ) region is excluded by the R parameter and CAST likelihoods (cf. Fig. 12) except at very low values of E/N . If axions constitute all of the DM in the Universe, these constraints are not relevant because the realignment mechanism cannot produce enough DM when |θ i | ≤ 3.14159 and m a,0 ∼ > 1 meV (cf. right panel of Fig. 7), so the high-mass region is excluded.
We also see slightly lower profile likelihood values for masses below about 0.1 µeV. This is due to the role of the axion-electron coupling in the R parameter likelihood: while 16 The assignment of weights to the different discrete values or to the different parts of the band is not trivial; it becomes complicated quickly if we consider the general QCDAxion family instead of specific DFSZ-type and KSVZ-type models, because the number of additional components (Higgs doublets or heavy quarks) is not fixed. Although it could be argued that all possible values of E/N are equally likely within each class of model, models with more additional particles might be considered more "contrived", and hence less probable a priori. This is relevant because it particularly affects the higher values of E/N , which cannot be achieved in the simpler models with only one new quark or two Higgs doublets. Creating a probability density function based on how often the values of E/N occur for all the different cases might hence not reflect the a priori probability for each version of the QCDAxion model. There are also significant practical challenges to computing all possible values of E/N when the number of quarks becomes very large, as well as in the most general versions of DFSZ-type models. our updated value for the helium abundance reduces the tension between theory and observations, there is still a slight preference for g aee = 0. For small masses, however, the maximum allowed value for the axion-electron coupling, C aee ≤ 10 4 , is still not large enough to satisfy this small preference.
In the bottom row of Fig. 15, we show the allowed values for the magnitude of the initial misalignment angle, with and without the assumption that axions constitute all of DM. Due to the influence of the various nuisance parameters and the relic density likelihood, the allowed region in the right panel is not simply a line, but a band of parameter combinations that reproduce the observed DM density within the allowed uncertainties. This panel also illustrates the well-known result that the initial misalignment angle needs to be fine-tuned, i.e. |θ i | 1, for QCD axion masses of m a,0 ∼ < 0.1 µeV. 17 We will investigate this issue in more detail below using a Bayesian analysis. Bayesian results. Breaking the PQ symmetry before inflation effectively results in a single, homogeneous value for the misalignment angle in the entire observable Universe. This gives a physical motivation for choosing a flat prior on θ i . Parameter regions in which θ i must be very small to avoid axion overproduction are hence theoretically less appealing. In a Bayesian analysis, we can see and quantify these fine-tuning issues in the (marginalised) posterior distributions, which quantify the degree of belief in certain values of the parameters given data and prior information.
We show marginalised posteriors for the QCDAxion model in Fig. 16, once again without (left) and with (right) the requirement that QCD axions are all of DM. As a consequence of fine-tuning in θ i , the low-mass (high-f a ) region of the parameter space in Fig. 16 is disfavoured, even when taking the DM relic density as an upper limit only. This is because in the low-mass region, large absolute values of the initial misalignment angles have a small likelihood. An O(1) value for the magnitude of the initial misalignment angle is a priori more probable than finding a value close to zero, due to the flat prior. This conflict leads to fine-tuning becoming increasingly necessary as the axion mass decreases, which is penalised in the Bayesian analysis. Although such parameter combinations might still give valid solutions that evade all constraints, they are not as probable as others.
A similar logic applies to large axion-photon coupling, i.e. large E/N . Due to the fine-tuning in E/N necessary to evade the helioscope and R parameter constraints at large axion mass (cf. the corresponding profile likelihood in the top left panel of Fig. 15), the large-m a,0 (low-f a ) region in the top left panel of Fig. 16 is disfavoured in the Bayesian posterior.
If we demand that axions explain all of DM, the consequences are even more dramatic. The most probable axion models are confined to the narrow band in m a,0 , visible in the upper right panel of Fig. 16. This mass range presents a feasible target for haloscope searches, and ADMX in particular is already beginning to cut into these models from the left (low-mass end). This also explains why the band of m a,0 -θ i values in the bottom right panel of Fig. 16 is not continuous, but disrupted around two points. These correspond to the ADMX and RBF/UF haloscope searches, respectively. While the RBF and UF haloscopes cannot reach as far down into the coupling space as ADMX, they do still constrain a significant fraction of the coupling range.
The fact that the Bayesian analysis singles out a well-defined range for the QCD axion mass becomes even more apparent in Fig. 17, where we compare one-dimensional profile likelihoods and marginalised posteriors for the axion mass. The frequentist approach does not yield a clear preference for any mass range, but the posterior distributions are strongly peaked around m a,0 ∼ 100 µeV. Clearly, such a result is not completely prior-independent, and we discuss the impact of adopting different priors in Appendix D. Nevertheless, it is appealing that a Bayesian analysis can identify a preferred region of m a,0 which, intriguingly, falls into the range that can be covered by experimental searches. Indeed, the impact of ADMX and other haloscopes already manifests itself as dips in the right panel of Fig. 17. for Ω a h 2 in QCDAxion models with upper limits (red shading) and matching condition (blue shading) for the DM relic density.
The marginalised posterior in Fig. 17 allows us to infer a preferred QCDAxion mass range. When demanding that axions explain all of DM, we find that the 95% equal-tailed credible interval for the axion mass is 0.12 µeV ≤ m a,0 ≤ 0.15 meV; allowing them to constitute a fraction of DM, this becomes 0.48 µeV ≤ m a,0 ≤ 3.8 meV. These numbers have minimal dependence on the adopted prior for E/N , but a stronger dependence on the choice of priors for C aee and f a (Appendix D).
We also note that if the PQ symmetry is broken after inflation, the preferred axion mass range will generally shift to larger values due to averaging of the energy density and inclusion of topological defects (cf. Sec. 2.3). The lower bounds on m a,0 that we quote can therefore be viewed as robust against changes of assumptions about inflation.
Finally, we also show the one-dimensional profile likelihoods and marginalised posteriors for the QCDAxion relic density in Fig. 18. Demanding that axions be all of the DM effectively dominates the outcome of this analysis. Using the DM relic density as an upper limit causes the profile likelihood to essentially follow the relic density likelihood function (left panel). In a Bayesian analysis, however, we immediately see that QCDAxions are not expected to generally provide all of the DM in the Universe, given our definition of the parameter space and priors. Imposing the DM relic density as an upper limit, the median axion relic density is 6.5 × 10 −3 , or about 5% of the observed DM abundance. The 95% credibility equal-tailed preferred range is 6.8 × 10 −6 ≤ Ω a h 2 ≤ 0.10, which corresponds to between about 0.006% and 90% of the cosmological density of DM. This demonstrates that in the pre-inflationary PQ symmetry-breaking scenario, although QCDAxions can provide a sizeable contribution to the DM density of the Universe, they probably do not contribute all of DM. Again, we stress that these statements are sensitive to the adopted prior on C aee and f a (Appendix D).

DFSZ-and KSVZ-type models
The DFSZ-type (DFSZAxion-I, DFSZAxion-II) and KSVZ-type (KSVZAxion) models differ from their parent model, the QCDAxion, in that they specify the axion-photon and axionelectron coupling strengths, or at least limit them to a well-defined range for a given axion mass. They are but a few of the many possible phenomenologically-inspired models, but they serve as interesting archetypes of their respective subclasses to compare with more general QCDAxion models.
Prior choices. Our prior choices for the DFSZ and KSVZ model can be found in Table 7.
For most of them, the rationale is the same as for QCDAxions presented in Sec. 5.3. The only differences are in the parameters related to couplings. We fix E/N to some typical values considered previously in the literature [e.g. 16]. The range that we choose for tan(β ) in DFSZ-type models reflects the values allowed by perturbativity bounds [229]. Our choice of a log prior for tan(β ) reflects the assumption that each possible Higgs vacuum expectation value is equally likely; indeed, any sensible prior choice for this parameter should reflect the fact that the two Higgs doublets may be interchanged, and the prior should be invariant under inversion of the ratio of vacuum expectation values.
Frequentist results. The DFSZ-and KSVZ-type models are essentially restrictions of the allowed QCD axion couplings. We therefore only expect to see qualitative differences in the results from the different models where the DFSZ and KSVZ interaction strengths cannot be tuned sufficiently to evade constraints from haloscopes (if axions make up all of DM) and the R parameter. Figure 19 shows the profile likelihood constraints on various axion models, imposing the DM density as an upper limit. We can see that the upper limit on the axion mass depends on the value of E/N in a given model (cf. Fig. 12), giving KSVZAxion models with E/N = 5/3 the largest allowed parameter space out of all the models compared here. The  Demanding that axions explain all of DM, Fig. 19 would change slightly. All models with m a,0 ∼ > 0.1 meV would be ruled out (not being able to provide all of the DM through the realignment mechanism), and ADMX would partially constrain all models except those with E/N = 5/3 (cf. Fig. 12).
The relation between the DFSZ-and KSVZ-type models and their parent QCDAxion model determines their allowed axion-electron couplings. Figure 20 shows how the allowed parameter space in the m a,0 -C aee parameter plane of the QCDAxion model is constrained further by imposing additional relations between the different model parameters in the KSVZ-and DFSZ-type models. This is most striking in the case of the KSVZ-type models, for which C aee is only induced at the loop level and depends directly on m a,0 (3.10). Note that the ordering of the KSVZAxion regions is also non-monotonic in E/N due to the difference term in (3.10). The finite sizes of the allowed parameter regions are simply a result of the nuisance parameters included in the relation between C aee and m a,0 . For DFSZ-type models, C aee depends on the additional parameter tan(β ), which makes it possible to accommodate a wide range of axion-electron couplings. However, the parameter space is also more constrained in this case, as very large values of C aee cannot be realised given other constraints on tan(β ). Also note that, due to the different coupling structure in DFSZAxion-I and DFSZAxion-II models (3.11), the same range for tan(β ) translates into a lower minimal value of C aee in DFSZAxion-II models than in DFSZAxion-I models. The resulting possible range for C aee in DFSZAxion-II models (and KSVZAxion models with E/N = 5/3) also extends to slightly lower values than the prior box that we chose for QCDAxions.
Bayesian results. We use the nested sampling package MultiNest to estimate the Bayesian evidences Z(M) for each model M. From these we construct the Bayes factor [230][231][232] where M 1 and M 2 are the two models under investigation, θ 1 and θ 2 are their parameters, P 1 and P 2 are their priors and L is the likelihood. The Bayes factor is connected to the odds, i.e. the ratio of posterior probabilities, of the models being correct: In this paper, we assign equal prior probabilities to both models being correct, i.e. choose P(M 1 )/P(M 2 ) = 1, so the Bayes factor is the same as the posterior odds ratio. Using MultiNest's nested sampling (as opposed to importance nested sampling) estimates for evidences, we calculated the odds in favour of KSVZ-and DFSZ-type axion models over the QCDAxion model. In terms of the commonly used scale for Bayes factors [230,231], we find that there is generally no noticeable evidence for or against any of these models, which would require an odds ratio of more than 3:1 (or less than 1:3).
Imposing the relic DM density as an upper limit, the odds in favour of any KSVZor DFSZ-type axion models, compared to QCDAxions, are 2:1. If we demand that axions constitute all of DM, the odds reduce to 1:1.
The outcome of the model comparison is not surprising, as we have not included any positive evidence for axions at this stage. We will discuss in the following section how these conclusions change when including WD cooling hints.

Cooling hints
Observables related to stellar cooling offer a unique opportunity to constrain the axion parameter space. If future observations confirm the need for additional cooling channels  to explain the observed decreases in WD pulsation periods, we may be able to use WDs to measure the axion mass and coupling strengths. In this section we add the likelihoods related to the WD cooling hints to our analysis, emphasising once again the caveats and difficulties associated with assigning uncertainties to the model predictions (cf. Sec. 4.5.4). Here our prior choices for each model are the same as in the preceding sections. A detailed numerical comparison of our results to previous works [13,16] is not meaningful due to differences in the choice of WD likelihood function, but the findings are qualitatively similar.

QCD axions
Previous studies have mostly considered the phenomenological couplings g aγγ and g aee or specific QCD axion models with fixed E/N . Here, we investigate which parts of the broader QCDAxion parameter space can explain the cooling hints.
Frequentist results. Figure 21 is the cooling-hint equivalent of Fig. 15, summarising the allowed anomaly ratio and magnitude of the initial misalignment angle. The only notable difference is at m a,0 ∼ < 0.1 µeV (f a ∼ > 3 × 10 13 GeV), where none of the possible values for C aee under consideration is large enough to fully account for the anomalous cooling. Recall that g aee ∝ C aee m a,0 (3.7) and that the cooling hints point towards a relatively narrow range of couplings g aee (Fig. 11). The overall effect of the cooling hints is therefore to disfavour lower masses. Had we chosen the range of possible values for C aee to be smaller, these constraints would extend to even larger values of m a,0 (and vice versa if we had permitted even larger values of C aee ).
The right panels of Fig. 21 show that QCDAxion models can satisfy the cooling hints and be all of the DM in the interval 0.1 ∼ < m a,0 /µeV ∼ < 300. The lower bound on this mass region depends on the largest allowed value for C aee .
It is interesting to consider how good the fit of the QCDAxion model is in an absolute sense. Most constraints are easily satisfied by the best-fit point, such that the corresponding partial likelihoods give p-values of order one, which we will not discuss further. 18 One exception is the fit to the temperature dependence of the QCD axion mass, which gives a p-value of order 10 −5 (see Sec. 3.2). Ignoring this likelihood (and the two model parameters constrained by it) we are left with 7 model d.o.f. and 48 data d.o.f. when including the WD cooling hints; without the cooling hints, the data d.o.f. is 43. The corresponding p-value is 0.30 with cooling hints included, and 0.60 without. The decrease in p-value when including the WD cooling hints results from the slight discrepancies between the cooling hints themselves (cf. Fig. 11) and their slight tension with the R parameter likelihood.
Bayesian results. Selected results from the Bayesian analysis of QCD axions in combination with cooling hints can be found in Fig. 22. Compared to the Bayesian results without cooling hints in Fig. 16, we can see that the preferred mass regions in the m a,0 -θ i plane get narrowed down slightly when we impose the DM relic density constraints as an upper limit (bottom left plot). However, for the anomaly ratio, E/N , this is not the case (top left plot). Generally speaking, these results identify the most credible regions for a compromise between QCD axions fitting the cooling hints and "naturally" not overproducing DM (which prefers masses of O(1 to 100 µeV), cf. Fig. 16). Despite the slight differences, which might also depend on the adopted priors, the overall results with and without cooling hints are remarkably similar. This is mainly due to the influence of the R parameter likelihood included in both cases and its slight preference for non-zero couplings.
The influence of the cooling hints is illustrated further in Fig. 23, which shows the regions of highest posterior probability in the m a,0 -C aee parameter plane with and without the inclusion of cooling hints. Since the cooling hints strongly require g aee ∼ 3 × 10 −13 (cf. right panel of Fig. 11), we find the highest posterior probabilities along a line of constant C aee m a,0 . The chosen range of C aee then implies that the cooling hints can only be explained for m a,0 ∼ > 0.3 µeV. For values of C aee ∼ 1, the cooling hints would point towards meV-scale axions, which is incompatible with the requirement Ω a ∼ Ω DM . This regime is therefore disfavoured in the right panel of Fig. 23. Not including the cooling hints  essentially only results in a upper limit in the most credible mass regions close to where the line of constant g aee was (grey contours in Fig. 23), due to the R parameter likelihood.
As mentioned before, QCD axions can account for both the cooling hints and all of the DM in the Universe (cf . Figs 21 and 22). However, because the posterior probability in Fig. 22 is normalised, one cannot infer from this plot if these solutions occur naturally or if considerable fine-tuning is required. Figure 24 gives an idea of the "naturalness" of QCDAxion DM by showing the marginalised posterior as a function of m a,0 and Ω a h 2 (with and without the WD cooling hints). In the colour density plots of both panels in Fig. 24, we can see that the scan finds credible parts of the parameter space where axions account for a sizeable fraction of the DM while being consistent with all experiments and observations. The differences between including and not including the cooling hints regarding the preferred regions of Ω a h 2 are rather small, consistent with the other plots.
Similar to the discussion at the end of Sec. 5.3, we can infer the most credible regions for the relic abundance of axions as well as for m a,0 . The preferred axion mass is very   for QCDAxion models with upper limits (density plots and black contour lines) and matching condition (grey contour lines) for the observed DM relic density. We show the constraints on the energy density in axions today, Ω a h 2 , without (left) and with (right) the inclusion of cooling hints. similar with or without the inclusion of the cooling hints. Imposing the DM relic density as an upper limit, we find 0.70 µeV ≤ m a,0 ≤ 2.8 meV at 95% credibility (equal-tailed interval); demanding that axions be all of the DM, this becomes 0.41 µeV ≤ m a,0 ≤ 0.14 meV. Including the cooling hints slightly modifies the preferred range for Ω a h 2 . At 95% credibility (equal tails), 1.0 × 10 −5 ≤ Ω a h 2 ≤ 0.10, corresponding to 0.009-90% of DM. The median value is Ω a h 2 = 9.3 × 10 −3 , or about 8% of the observed DM.
Finally, let us return to the m a,0 -g aγγ parameter plane, as discussed in the GeneralALP model without cooling hints (see Fig. 12). In Fig. 25, we contrast the naïve bounds on the parameter space (from the phenomenological constraints on GeneralALP models and the maximum possible value of E/N ) with the regions preferred by a Bayesian analysis. These regions are not only determined by the constraints from data (satisfying the cooling hints in both panels and matching the DM density in the right panel), but also by the fine-tuning in some parts of the parameter space. Fine tuning is necessary for avoiding overproduction of DM at small m a,0 , and for achieving low values of g aγγ through cancellations between E/N and C aγγ at large m a,0 (cf. 3.8). For our adopted priors, the most credible parameter regions correspond to a few orders of magnitude around m a,0 ∼ 10 µeV and g aγγ ∼ 10 −12 GeV −1 . In Appendix D we discuss how choosing different priors may affect these conclusions.

DFSZ-and KSVZ-type axions
DFSZ-type models have intrinsically larger coupling to electrons than KSVZ-type models, which only obtain their interactions with electrons at loop level. DFSZ models are therefore the natural choice to account for the potential WD cooling anomalies [16], by way of an axion-electron coupling of g aee ∼ 3 × 10 −13 .
Frequentist results. Fig. 26 shows the profile likelihood for the DFSZAxion-I and DFSZAxion-II models, compared to the band that maximises the profile likelihood for the QCDAxion model. Clearly both DFSZ-type models can accommodate the cooling hints with large axion-electron couplings. We already saw in Fig. 20 that even without the cooling hint likelihood, the KSVZAxion fails to achieve large C aee values; with cooling hints included, the highest likelihood regions for KSVZAxion models are therefore essentially the same as in Fig. 20. Table 8 gives the best-fit values for the six classic axion models that we consider in this section, under the requirement that they do not exceed the observed abundance of DM. We do not report best-fit values for |θ i |, as even with the maximum value included in  Table 8: Best-fit values for DFSZAxion-I, DFSZAxion-II, and KSVZAxion models when imposing an upper limit on the DM relic density, Ω a ∼ < Ω DM . In the final column, we compare the likelihood of the respective best-fit points to QCDAxion models. Note that the QCDAxion model has two more degrees of freedom than the KSVZ models, and one more than the DFSZ-type models.  our scans, axions only account for a few percent of the observed DM abundance. For each model, we calculate ∆ ln(L), the logarithm of the ratio of the best-fit likelihood relative to the QCDAxion model. As anticipated, DFSZ-type models perform much better than KSVZAxions (at the expense of having one additional degree of freedom). This is because DFSZ-type models can easily reach the required axion-electron coupling to fit the cooling hints with masses of order m a,0 ∼ 10 meV (as noted previously [16]). For KSVZ-type axions, the masses required to naturally fit the cooling hints are about three to four orders of magnitude larger (see e.g. 3.10), and the associated axion-photon coupling is therefore in conflict with the R parameter likelihood (as well as hot DM bounds, which are not included). Nevertheless, even for KSVZAxion models there is still a slight preference compared to having no axion at all, which corresponds to −2∆ ln(L) = 10.54.

E/Nm
For DFSZ-type models, the DFSZAxion-I model gives a better fit than the DFSZAxion-II model. This is due to the influence of the R parameter likelihood, which combines with the cooling hints to force the axion-photon coupling to g aγγ ∼ < 2 × 10 −11 GeV −1 (at 95% C.L.). The best-fit point is therefore a balance between reaching high enough g aee and minimising g aγγ . The maximum value of C aee for DFSZAxion-I models is about a factor of 1.08 larger than for DFSZAxion-II models, due to perturbativity constraints on tan(β ), whereas the axion-photon couplings are about a factor of 0.6 lower, yielding a better fit to the R parameter likelihood at any given mass.
Next, let us briefly consider the case where we demand that the classic axion models provide all of the DM. This results in much poorer maximum likelihood values of around −2∆ ln(L) ≈ 10.5 compared to QCDAxion models for all KSVZ-and DFSZ-type models that we consider in this paper. This is not surprising because, unlike some other QCDAxion models, they cannot account for both the cooling anomalies and DM. 19 The maximumlikelihood regions are also highly degenerate in this case, as none of the fits is actually "good" (with a maximum likelihood comparable to the case without any axion).
Finally, thanks to the model hierarchy shown in Fig. 1, we can perform nested hypothesis tests to determine whether or not the more constrained, specific DFSZAxion and KSVZAxion models are disfavoured compared to the broader class of QCDAxion models. Without the cooling hints, QCDAxion, DFSZAxion and KSVZAxion models cannot be discriminated. This situation changes if the WD cooling hints are taken into account. Imposing the measured DM density as an upper limit only, KSVZAxion models (null hypothesis) can be rejected with respect to QCDAxion models (alternative hypothesis) with a p-value of 8.3 × 10 −3 . DFSZAxion models as a null hypothesis, on the other hand, cannot be rejected (p ≥ 0.26). If we instead demand that axions be all of DM, both DFSZAxion and KSVZAxion models can be rejected with respect to the QCDAxion model with p-values of 1.2 × 10 −3 and 5.3 × 10 −3 , respectively.
Bayesian results. Figure 27 shows posteriors for the DFSZ-type models compared to the QCDAxion, in the m a,0 -C aee parameter plane (see also Fig. 23). These regions resemble the ones found in Fig. 26. The log prior on tan β enables both models to achieve sufficiently large values of C aee quite naturally, despite their structural differences. At first sight, DFSZAxion-II models appear to be able to occupy a larger region of parameter space in the 68.27% C.R. However, this is mainly a reflection of the effective prior on C aee , which for the DFSZAxion-I model is strongly peaked towards the upper boundary of the accessible parameter space in Fig. 27. In contrast, the effective prior in C aee is almost flat in the DFSZAxion-II model. In terms of the fundamental parameters (i.e. tan β ), the credible parameter region in the DFSZAxion-I model is in fact larger, because the large values of C aee needed to explain the cooling hints can be achieved more easily.   Table 9: Odds ratios in favour of DFSZ-and KSVZ-type models, compared to the parent QCDAxion, as calculated from the nested sampling evidence estimates in MultiNest, and including the WD cooling hints. We either impose the DM relic density as an upper limit (Ω a ∼ < Ω DM ) or demand that axions be all of DM (Ω a ∼ Ω DM ). Note that the estimated uncertainties on the evidence values are small enough that the corresponding uncertainties on the odds ratios are negligible.

Model
DFSZAxion-I DFSZAxion-II KSVZAxion To investigate the consequences of fine-tuning in more detail, we consider all models in a Bayesian model comparison. The resulting odds ratios can be found in Table 9. If we impose the DM relic density as an upper limit, the odds ratios are still mostly inconclusive. However, for the DFSZAxion-II and KSVZAxion models, the trend swings in favour of the broader class of QCDAxion models when cooling hints are added to the analysis. In contrast, the DFSZAxion-I model fares better than the QCDAxion, with an odds ratio of 3:1. If we combine this 3:1 odds ratio with the 2:1 preference for the QCDAxion model over KSVZAxion models, there is a 6:1 positive preference for the DFSZAxion-I model over all KSVZAxion models. This preference is not surprising, given the differences in the natural axion-electron coupling strength in these models.
If we demand that axions solve the cooling hints and constitute all of DM, model comparison confirms quantitatively (with ratios of ≤ 1:5) that there is a positive preference for the QCDAxion over the DFSZ-and KSVZ-type models. This is because we allowed for much larger values of C aee with QCDAxion models than with DFSZ-and KSVZ-type models.
In fact, the one-dimensional marginal posterior for the QCDAxion electron coupling peaks at C aee ∼ 100 (C aee ∼ 50 if we impose the DM relic density as an upper limit), whereas the DFSZ and KSVZ models are limited to couplings of less than one. However, it should be noted that QCDAxion models with such a large coupling are not expected from traditional axion models, and may therefore pose a challenge for model building.

Conclusions
In this study we presented the first global fits of axion models in the pre-inflationary PQ symmetry-breaking scenario, using frequentist and Bayesian methods. We identified the most viable regions of parameter space for these models and discussed the effect of adding cooling hints seen in white dwarfs. We not only considered the phenomenological parameter space, but also the underlying parameters in a generic QCD axion model and six specific DFSZ-and KSVZ-type models. We extended previous results in the literature by including various nuisance parameters and by quantitatively considering the fine-tuning of the initial misalignment angle.
We found, in agreement with previous work, that QCD axion models are viable as explanations for both the observed cold dark matter and the white dwarf cooling hints. We showed that, for a broader class of QCDAxion models than what is traditionally considered, these can be achieved simultaneously. We also quantitatively confirmed that this is not possible for six specific DFSZ-and KSVZ-type models. However, if the condition of being all the dark matter in the Universe is relaxed, we found that DFSZAxion-I models are positively preferred over DFSZAxion-II models, and over all four KSVZAxion models that we investigated.
We determined the most credible predicted ranges for the QCD axion mass and its cosmological abundance in the Bayesian statistical framework. Although these results are somewhat prior-dependent, QCDAxions appear likely to be a cosmologically relevant (but probably not dominant) component of dark matter. Moreover, the most credible axion mass range is within reach of current and planned haloscope experiments.
Global fits of QCD axions and axion-like particles have the potential to confirm and refine previously known phenomenological statements about the relevant parameter spaces (e.g. exclusion limits, existence of fine-tuning, compatibility with cooling hints and dark matter). In other cases, they offer new and more rigorous insights (e.g. model comparison, most credible parameter regions). This is true in particular for the Bayesian analyses that we performed in this paper, because they inherently take fine-tuning into account. Due to the orthogonality of constraints and the insufficient sensitivity of most experiments, it is not (yet) possible to decisively probe the axion parameter space in the pre-inflationary PQ symmetry-breaking scenario, but it is realistically possible to target the most likely versions of QCD axion and axion-like particle models in the near future (see Ref. [5] for a recent review of upcoming axion searches).
The axion routines that we developed for this paper will be publicly available in DarkBit within GAMBIT 1.3.0, to be released shortly. Statistical samples and input files from this study are freely available from Zenodo [35].

Acknowledgments
We are grateful to the UK Materials and Molecular Modelling Hub for computational resources, which is partially funded by EPSRC (EP/P020194/1), PRACE for awarding us access to Marconi at CINECA, Italy, and the Imperial College Research Computing Service as well as to all those individuals and experimental collaborations who kindly provided their data and results. We thank A. Serenelli

A Integrating solar models for the signal prediction in CAST
To obtain the axion flux at Earth from axion-photon interactions, we need four inputs from a solar model: the solar radius R , the temperature T (r), the plasma frequency ω pl (r) and the screening scale κ s (r) (see Eqs. 4.5-4.8). While T (r) can be obtained directly from solar model files, ω pl (r) and κ s (r) have to be calculated using the mass density ρ(r) and the mass fractions X i (r) for each ion/atom with label i using (4.6). Assuming that the plasma is fully ionised, we can recast the equations into the following form [210]: and where m u is the atomic mass unit, r is the distance from the centre of the Sun in units of R , and Z i and A i are the charge and atomic weight of the ith element. Note that the approximation of full ionisation is justified almost everywhere inside the Sun, i.e. r 0.95. 20 As the largest contribution to the axion flux comes from the innermost region (r ∼ < 0.2 [122]), we can safely employ this assumption. For elements tracked by the solar model in bulk (i.e. without isotopic information), we calculate the mean atomic weights A i using the isotopic composition of Ref. [162], with values based on the terrestial composition updated to use more recent data [233]. The expected number of photons in the energy range [E j , E j+1 ] is where E is the effective exposure and dΦ/dE is the combined photon spectrum from axion-electron and axion-photon contributions. We neglect the energy dispersion of the CAST detector, as it is always less than about 0.2 keV [122], and therefore smaller than or comparable to the bin width of the CAST analyses (0.3 keV and 0.5 keV respectively for the 2007 and 2017 analyses). We provide tabulated data within DarkBit for the effective exposure E; 21 this can be found in DarkBit/data as dataset_EffectiveExposure.dat, where dataset is either CAST2007 or CAST2017_X, where X corresponds to the data sets A to L in Ref. [123]. All other files that we mention in the following can be found in the same folder. Note that performing the energy integration in (A.3) after the density integral in (4.8) is only possible because all contributions to the integral are for energies greater than ω pl (0) ≈ 0.3 keV. For energies lower than ω pl , axions cannot be produced from the plasma and the square root in (4.8) becomes ill-defined.
To obtain the contribution to dΦ/dE from axion-photon interactions, we integrate (4.7) and (4.8) over r and ρ, using the adaptive 51 point Gauss-Kronrod rule gsl_integration_ qag, from the gsl library. We obtain quantities necessary for these integrations (temperature, etc.) by interpolating the solar model linearly in radius. For the contribution to dΦ/dE from axion-electron interactions, we use the spectrum published in Ref. [147] and redistribute it as Axion_Spectrum_AGSS09met_old_gaee.dat.
The peaks in the spectrum from axion-electron interactions specifically require using an algorithm that takes them into account as singularities. We therefore compute the contribution from axion-electron interactions to the signals s j (A.3) using the gsl_ integration_qagp integrator of the gsl library. 22 For the axion-photon contribution to the signal, we again use the 51-point Gauss-Konrod rule. We perform all signal calculations to a relative accuracy of 10 −6 .
We compute signals at reference values of g aγγ = 10 −10 GeV −1 and g aee = 10 −13 , at 183 mass values ranging from 10 −3 to 10 2 eV. We use unequally spaced mass values, as the density of points needs to be higher in certain regions to obtain good interpolating functions. We provide these pre-calculated signal count files as dataset_ReferenceCounts_ solarmodel_coupling.dat, where dataset refers to CAST2007 or CAST2017_X results, solarmodel is the name of the solar model, and coupling is either the contribution from axion-photon interactions (gagg) or axion-electron interactions (gaee).
For the axion-photon contribution (coupling = gagg), we provide these files for the If these options are left unspecified, the defaults will be chosen. For the axion-photon coupling, this is AGSS09met [163]; for the axion-electron coupling, the default is AGSS09met_ old [161]. If the corresponding dataset_ReferenceCounts_solarmodel_coupling.dat file cannot be found, DarkBit will attempt to create such a file by solving (A.3). For the axionphoton coupling, it will attempt to do this on the basis of a provided solar model file, which should be called SolarModel_solarmodel.dat. We have included the solar model file SolarModel_AGSS09met.dat from Ref. [163] as an example of the expected formatting.
For the axion-electron coupling, DarkBit will attempt to generate the relevant dataset_ ReferenceCounts_solarmodel_gaee.dat file on the basis of a provided spectrum file Axion_ Spectrum_solarmodel_gaee.dat. If the relevant solar model or spectrum file cannot be located, or has the wrong format, GAMBIT will terminate. Note that this will currently happen for any choice of solar_model_gaee except for the default value. However, if the user can provide the spectrum file for any other model, the calculation will go ahead.

B Numerical implementation of the solution to the axion field equation
The most general equation of motion for a scalar field φ(t, x) in the background of a Friedmann-Robertson-Walker-Lemaître universe with Hubble parameter H ≡ȧ/a and potential V is given by [60, app. B.12] φ + 3H(t)φ + ∆φ a 2 (t) while the energy density associated with φ can be calculated using: To obtain the energy density in axions today, it is sufficient to solve the homogeneous field equation for θ ≡ a/f a , the normalised axion field: We assume the potential for QCD axions and ALPs alike. For QCD axions, this is a reasonable approximation to the full potential [52]. The axion energy density ε a ρ a from (B.2) is therefore given by The initial conditions for (B.3) are some initial misalignment angle θ(0) ≡ θ i and vanishing initial derivativesθ(0) ≡ 0. 23 For early times, the solution of (B.3) is a constant value, θ(t) = θ i . The behaviour only changes around the oscillation time T osc , defined as the point when 3H(T osc ) = m a (T osc ). 24 Assuming that the scale factor a is a monotonically increasing function of cosmic time t, and that the temperature T is a monotically decreasing function of t, we may re-write (B.3) as where X = t, α, τ with α ≡ a/a osc , τ ≡ T /T osc and where × denotes derivatives w.r.t. τ . Conservation of entropy can subsequently be used to relate α and temperature: Equation (B.9) has the advantage that an explicitly temperature-dependent axion mass or the effective relativistic degrees of freedom, g and g S , can be easily incorporated. The various terms can be further simplified as follows: Rescaling ϑ ≡ θ/θ i , we may use (B.9) to rewrite (B.6) as where the auxiliary functions F i for i = 1, 2 are given by 14) 23 The initial value for the derivative can be taken as zero assuming that |θ/H| 1 at early times [69, ch. 10.3.2]. The consequence of that assumption is essentially independent of the cosmological history and equation of state [180]. 24 We will solve the field equation numerically around this point and the exact definition of the onset of oscillation is therefore irrelevant. See e.g. Ref. [40] for a more detailed discussion of the general relevance of this definition.
where we assume a radiation-dominated universe, H ∝ g(T ) T 2 , for the last line. We use interpolated values for (B.14) and (B.16) based on the analytic forms for the effective degrees of freedom in Ref. [65]. The tables used for this procedure are included in GAMBIT.
Note that if the changes in the effective degrees of freedom are not important, F 1 (τ ) ≈ −1, F 2 (τ ) ≈ 0, and the field equation (B.13) reads in the harmonic limit of |θ i | 1: i.e. we obtain the approximate behaviour of a (damped) harmonic oscillator. In that limit, the comoving axion number density n a ≡ ρ a a 3 /m a is conserved on average. We can therefore stop integrating the differential equation once specified conditions for m a /3H and θ are met. Because the solution θ(t) oscillates as a function of time or temperature, we define a peak amplitudeθ in order to check that |θ(t )| 1 from some time t : .
In the harmonic limits, (B.18) indeed coincides with the amplitude of oscillation when all the energy of the field is stored in the potential part of the energy density. 25 Once the axion field starts to oscillate, fulfilment of the conditionθ 1 indicates that the harmonic limit has been achieved. This provides much clearer limiting behaviour than e.g. |θ|.
We use the gsl implementation of Brent's method to solve for T osc to a precision of 10 −6 . We then impose the initial conditions at some relatively high temperature τ start = 10 5 , and evolve the differential equation with gsl_odeiv2_step_bsimp until some relative temperature τ stop for which m a /3H > 10 3 andθ < 10 −2 . We then calculate the energy density from realignment axions today using the conservation of comoving number density where T CMB is the CMB temperature.

C H.E.S.S. likelihood implementation details
Our approximation of the H.E.S.S. likelihood is based on the information provided in Figs 6 and 7 of Ref. [134] on limits on axion-photon conversion in galactic cluster magnetic fields. We choose 13 lines of constant axion-photon coupling (corresponding to specific values of 25 Note that due to the inequality 1 − cos(x) ≤ x 2 /2, we can apply this procedure even if θ is not much smaller than unity. In fact, it should be valid even if the true axion potential is not cosine-shaped, as long as the true potential is also harmonic for small values of θ.  the axion decay rate Γ in their figure) that intersect with the exclusion curves at least five times (dashed black and blue horizontal lines and blue circles in Fig. 28). We subsequently define a family of one-dimensional natural cubic splines along each dashed line by including two additional points outside of the known data. We optimise each of those splines by iteratively modifying the locations of the additional points along the dashed lines, until the values and first derivatives of the splines at the additional points go to zero (blue diamonds in Fig. 28). This forcibly prevents ringing, ensuring that the interpolated log-likelihood is negative everywhere. To evaluate the log-likelihood at a given axion mass on one of the horizontal dashed lines, we use the resultant interpolating cubic spline from the fitting procedure. If a coupling value of interest sits between two horizontal dashed lines, we interpolate linearly between the values at a given mass on the adjacent dashed lines. We assign zero log-likelihood to all points outside of the area defined by the exterior points (blue diamonds).
This leaves two more regions in Fig. 28. The first is at higher couplings, where Fig. 6 of Ref. [134] does not show any limits, but Fig. 7 does (dotted horizontal lines and blue triangles in Fig. 28). To complete the likelihood curve families in this region, we infer additional points (by eye; light blue squares in Fig. 28) and follow the procedure described in the previous paragraph. We assign zero log-likelihood to all couplings above the uppermost horizontal line in Fig. 28.
At low couplings, less than five contours are available for interpolating (inverted triangles in Fig. 28). Here, we construct a single cubic spline in the vertical direction (blue stars and vertical dot-dash line in Fig. 28), in order to determine an exterior point at low coupling (the diamond at the bottom of the plot). We define three parabolae, based on nine points: the set of two diamonds and four outer circles on the dashed black line around g aγγ = 10 −10.6 GeV −1 , the two stars below this line, and the one diamond at the bottom of the plot. We designate the lower parabola, passing through the exterior point (bottommost diamond) at low coupling and drawn as a dashed black curve, as the zero-log-likelihood contour. We use the middle parabola to model the 99% C.L. exclusion curve, and the upper one to model the 95% C.L. curve. These three parabolae allow us to map any point (represented by the white star) inside the lower parabola to a likelihood, by interpolating between the three parabolae. We do this using the one-dimensional spline constructed along the central vertical axis, using as input to the interpolation the position of the point in question between the central star and the zero-log-likelihood parabola. We carry out all interpolations in log space for both the axion mass and the axion-photon coupling.

D Prior-dependence of results
For QCD axions some prior choices are relatively straightforward (such as the use of a flat prior for θ i ) or do not have a significant impact on the results (such as the prior assignments for the well-constrained nuisance parameters). In general, however, the prior beliefs can have a significant impact on the posteriors and Bayes factors, especially in the absence of constraining data. In this appendix, we investigate different prior choices for the parameters f a , C aee , and E/N in order to assess the robustness of our results and to understand how our conclusions may be affected by a variation of priors. All the new results that we present in this appendix are based on MultiNest runs only, as we need the evidence values for examining the impacts of priors on model comparison, but are not sufficiently interested in the finer details of posterior maps to warrant also running T-Walk.
First, let us revisit the determination of the most credible regions for the axion mass in QCDAxion models from Sec. 5.3. The left panel of Fig. 29 shows the marginalised posterior distribution for the axion mass, m a,0 , for a number of alternative priors. The grey line and shading correspond to the prior choice made in the main text (right panel of Fig. 17). We can see that choosing a flat prior on either C aee or f a shifts the most credible region for the axion mass to lower masses. The effect is more extreme for a flat prior on f a than for a flat prior on C aee . On the other hand, restricting the prior of the anomaly ratio to E/N ≤ 170/3 (the highest values achievable in KSVZ-type models) does not have a significant impact on the result.
The different prior choices also alter the most credible regions of the axion-photon couplings (cf. Fig. 25). The resulting posterior distributions are shown in the right panel of Fig. 29. As expected, priors that imply a preference for smaller axion masses also lead to a preference for smaller values of g aγγ , while changes in prior that do not significantly affect the most probable region for m a,0 also do not impact the findings for g aγγ very much. In particular, for the alternative priors on E/N and C aee , the regions of highest posterior density still overlap with the region that we obtained with our adopted prior, whereas the regions resulting from a flat prior on f a are almost completely disjoint from the others. 68  Here we impose the observed DM relic density as an upper limit. We show the one-dimensional posterior for the axion mass without the cooling hints (left; to be compared to the right panel of Fig. 17) and the two-dimensional posterior for m a,0 vs g aγγ , including the cooling hints (right; to be compared with Fig. 25). The priors adopted in the main part of this study are marked in grey in the left panel, and in the colour scale and black contours in the right panel.
The preference for very small axion masses for the case of a flat prior on f a is a direct consequence of the relation between f a and m a,0 , which implies that the prior probability for m a,0 is strongly peaked at the lowest possible axion masses. This prior then overwhelms the preference for large m a,0 from the combination of tuning in the initial misalignment angle and the relic density requirement. In other words, the preferred mass range in this case is dictated almost entirely by the assumed prior. A log prior on f a , on the other hand, essentially corresponds to a log prior on m a,0 , which has the important advantage that any preference in the axion mass range is the result of phenomenological requirements on the model. A logarithmic prior for f a is therefore less informative than a flat one, which favours a particular scale of new physics near the upper end of the prior range.
However, it is also important to determine how changing the range of the logarithmic prior impacts our results. Taking the example of the grey curve in the left panel of Fig. 29 (canonical priors, no cooling hints, Ω a ∼ < Ω DM ), we find that restricting the log prior to f a ∈ [10 7 , 10 15 ] GeV results in a 95% equal-tailed credible interval for the axion mass of 0.50 µeV ≤ m a,0 ≤ 3.8 meV, compared to 0.48 µeV ≤ m a,0 ≤ 3.8 meV for the prior range that we adopted in the main paper (f a ∈ [10 6 , 10 16 ] GeV). Narrowing this range further to f a ∈ [10 8 , 10 14 ] GeV only raises the lower edge of the interval slightly: 0.63 µeV ≤ m a,0 ≤ 3.6 meV. Clearly the range of the logarithmic prior on f a has little impact on our results. The median values and credible intervals of Ω a are even more stable, with the lower bound of the 95% equal-tailed credible interval slightly increasing from 6.8 × 10 −6 to 7.3 × 10 −6 for f a ∈ [10 8 , 10 14 ] GeV (corresponding to about 0.006% of the cosmological density of DM, regardless of the prior range). Table 10: Odds ratios in favour of DFSZ-and KSVZ-type models, compared to the parent QCDAxion model, as calculated from nested sampling evidence estimates by MultiNest. We impose the DM relic density as an upper limit (Ω a ∼ < Ω DM ). Note that the estimated uncertainties on the evidence values are small enough that the corresponding uncertainties on the odds ratios are negligible.

Model
DFSZAxion-I DFSZAxion-II KSVZAxion We also calculate the Bayes factors for KSVZ-and DFSZ-type models compared to QCDAxion models, adopting the alternative priors for the appropriate models. The resulting odds ratios are given in Table 10, in addition to the values for the priors that we adopt in Sec. 5. Fig. 29, adopting a reduced range for the anomaly ratio E/N has little impact, slightly increasing the odds in favour of KSVZ-type models. Choosing the alternative prior on the axion-electron coupling constant, C aee , the odds in favour of KSVZ-and DFSZ-type models increase by an order of magnitude (if cooling hints are not included). This is because the flat prior on C aee causes a preference for large f a (smaller m a,0 ). This, in turn, implies a considerable fine-tuning in θ i in order to avoid overproducing DM. After including the cooling hints, axion-electron interactions are much more tightly constrained by data, and the prior-dependence of the odds ratios is reduced.

As in
Finally, a flat prior on f a favours the lowest allowed axion masses in all models. This results in even more fine-tuning of the initial misalignment angle than using a flat prior on C aee for QCDAxion models. However, although the Bayesian evidence in QCDAxion models is drastically reduced, the odds ratios in favour of KSVZ-and DFSZ-type models turn out to be 1:1. This is because the alternative prior is applied to all models: the fine-tuning required to overcome the prior pushes all evidences equally low, washing out any preference one way or another from the data.
We see that prior choice can indeed significantly influence the results of our analysis. This is true in particular for a flat prior on f a , where the prior-dependence dominates the results. On the other hand, we also saw that data with a strong preference for certain model parameter values, such as the cooling hints, can reduce the prior-dependence. In the face of more such data, the impact of prior choices tends to become less pronounced, ideally leaving only the physically-motivated prior on the initial misalignment angle able to influence final results.

E Overview of new capabilities
For reference, in Table 11 we provide a complete list of the new capabilities, dependencies and options that we have added to DarkBit whilst preparing this paper.