Constraints on Axions from Cosmic Distance Measurements

Axion couplings to photons could induce photon-axion conversion in the presence of magnetic fields in the Universe. This conversion could impact various cosmic distance measurements, such as luminosity distances to type Ia supernovae and angular distances to galaxy clusters, in different ways. In this paper we consider different combinations of the most up-to-date distance measurements to constrain the axion-photon coupling. Employing the conservative cell magnetic field model for the magnetic fields in the intergalactic medium (IGM) and ignoring the conversion in the intracluster medium (ICM), we find the upper bounds on axion-photon couplings to be around $5 \times 10^{-12}$ (nG/$B$) $\sqrt{\mathrm{Mpc/s}}$ GeV$^{-1}$ for axion masses $m_a$ below $10^{-13}$ eV, where $B$ is the strength of the IGM magnetic field, and $s$ is the comoving size of the magnetic domains. When including the conversion in the ICM, the upper bound is lowered and could reach $5 \times 10^{-13}\, $GeV$^{-1}$ for $m_a<5 \times 10^{-12}$ eV. While this stronger bound depends on the ICM modeling, it is independent of the strength of the IGM magnetic field, for which there is no direct evidence yet. These constraints could be placed on firmer footing with an enhanced understanding and control of the astrophysical uncertainties associated with the IGM and ICM. All the bounds are determined by the shape of the Hubble rate as a function of redshift reconstructable from various distance measurements, and insensitive to today's Hubble rate, of which there is a tension between early and late cosmological measurements. As an appendix, we discuss the model building challenges of the use of photon-axion conversion to make type Ia supernovae brighter to alleviate the Hubble problem/crisis.


Introduction
Axions, as periodic scalar fields, arise ubiquitously in both low-energy phenomenological models [1][2][3][4][5][6][7][8] and quantum gravity theories [9]. They serve as an important benchmark of feebly-coupled light particles beyond the Standard Model (SM). In particular, one of the most active experimental and observational targets is the coupling of an axion, a, to photons, which takes the form where F µν is the electromagnetic field strength,F µν = 1 2 µνρσ F ρσ is its dual field strength, and E and B are the electric and magnetic fields. The axion-photon coupling coefficient g aγγ has mass dimension −1 and is inversely proportional to a high energy scale. Various introductions to axion physics basics can be found in [10][11][12][13].
On the other hand, the past 20 years have seen an increased interest and corresponding progress in the efforts to measure various cosmic distances to chart out the expansion history of our Universe. One outstanding example is the measurements of luminosity distances (LD), D L , to Type Ia supernovae (SNIa). SNIa are used as "standard candles" in the Universe given their very similar peak brightnesses. A large number of SN surveys have brought about the Pantheon dataset, which is the largest and most accurate SNIa compilation at present [14]. It consists of a total of 1048 SNIa in the redshift range of 0.01 < z < 2.3, which can be used to constrain D L as a function of redshift z. Since D L (z) is determined by the Hubble expansion rate H(z), the Pantheon sample could consequently determine the shape of H(z), at late times. Type Ia supernovae samples are also a crucial input for late-time measurement of today's expansion rate H 0 , SH0ES [15,16], which is seriously at odds with the early time determination using the CMB data collected by the Planck satellite [17]. This is dubbed the "Hubble problem" or "Hubble crisis" (see [18,19] and references therein). In addition to LD measurements, there have been several kinds of precise measurements of angular diameter distance (ADD), which is defined as D A = d/θ for an astrophysical object of physical size d and angular size θ. Two examples we will use in this paper are from Baryonic Acoustic Oscillations (BAO) [20][21][22] and galaxy clusters [23,24].
These two seemingly unrelated subjects (axions and their couplings to photons in particle physics on the one hand, and cosmic distances in cosmology on the other) have an intriguing connection. The coupling of axions to photons in Eq. (1.1) suggests that in the presence of an external magnetic field photons could convert into axions, and vice versa. Indeed, since there could exist non-negligible magnetic fields in the intergalactic medium (IGM) and/or in the intracluster medium (ICM), the propagation of photons from astrophysical sources could be affected by their conversion into light axions in certain regions of the axion parameter space. This in turn could affect the inference of various distance observables. Take SNIa for example. Assuming the standard model of cosmology, ΛCDM, photons in the optical band converting into axions could result in a significant dimming of SNIa at higher z's, while SNIa at lower z's are less affected or not affected at all. An effective luminosity distance D L (z) which takes into account photon-axion conversion could then be constrained by the Pantheon sample, which is consistent with the prediction of pure ΛCDM. Other distance observables may be modified by photon-axion conversion as well, albeit in different ways. One such example is the angular diameter distance D A (z) to galaxy clusters. In summary, conversion of photons into axions is effectively equivalent to a departure of the Hubble diagram, H(z), from that of ΛCDM at late times, which could be constrained by various combinations of cosmic distance measurements.
In this work, we consider different combinations of cosmic distance measurements and carry out statistical analyses to map out the allowed parameter space in the plane of the axion's mass and coupling to photons, m a and g aγγ respectively. Analyses using cosmic distance measurements have been performed before, e.g., in Refs. [25][26][27] with an older and smaller dataset of SNIa. In addition to including more and updated datasets, our analyses differ from previous ones in the choice of observables. Earlier works usually interpret the constraints as arising from violations of the "Etherington relation" [28], the distance duality relation between D L and D A : D L (z) = (1 + z) 2 D A (z). In other words, the chosen observable is the ratio D L /D A . It is then (implicitly) assumed that the violation due to photon-axion conversion could be parametrized by a single parameter, e.g., such that D L = D A (1 + z) 2+ . As we will discuss in detail, depending on the datasets involved, D L and D A could be affected by photon-axion conversion in very different ways and the photonaxion conversion may not be encoded in a single function of z or a single parameter. Instead, we simply choose the observables to be those quantities directly measured or inferred in each dataset, such as the apparent magnitude of SNIa, the ADDs to galaxy clusters, and the characteristic angular scale of the matter two-point correlation function for BAO; and build corresponding likelihood functions.
The paper is organized as follows: in Sec. 2, we discuss the basic formalism of axionphoton conversion in the IGM or the ICM. We also explain how D L and D A could be affected by the conversion. In Sec. 3, we describe the datasets included in our analyses and the statistical method we use. In Sec. 4, we present and discuss the results as constraints on the axion parameter space. We conclude in Sec. 5. Throughout the paper, we assume that there is negligible axion production at SNIa and that, consequently, the effect of photon-axion conversion in IGM is to dim the SNs. In Appendix A, we will entertain the readers with the possibility of resonant axion production at SNIa, which might open up the possibility of brightening SNIa through IGM conversion. We will explain the related model building challenges and why this could not work as a solution to the current Hubble problem/crisis.

Axion-photon conversion
In this section, we will first review the basic formulas that describe photon-axion conversion in a magnetic field. We will then discuss the models of the two media, IGM and ICM, in which the conversion takes place. Lastly, we will discuss how various cosmic distances, D L to SNIa and D A to galaxy clusters, could be affected by the conversion in different ways.
Throughout the rest of this paper we will make frequent reference to the parameters that describe axion-photon conversion in a flat ΛCDM cosmological setting. As a shorthand, we denote these parameters as θ = {Ω Λ , H 0 , m a , g aγγ } and Θ = θ∪{M, r drag s }; where Ω Λ is the fraction of today's energy density in the cosmological constant, H 0 is the Hubble parameter, m a is the axion mass, g aγγ is the axion-photon coupling, M is the absolute magnitude of the SNIa standard candles, and r drag s is the comoving sound horizon size at the time of baryon drag.

Basic formulas
In the presence of external magnetic fields, the operator in Eq. (1.1) implies that the propagation eigenstates of the photon-axion system are mixtures of axion and photon states. As a result, there is a non-zero probability P 0 that a photon oscillates and converts into an axion while traveling through the magnetic field, effectively resulting in photon number violation. When birefringence and Faraday rotation effects are small, as is the case with propagation in the IGM [29], the axion mixes only with the photon polarization parallel to the component of the magnetic field B T , which is transversal to the direction of motion. In the simple case of photons with energy ω propagating in a constant and homogeneous magnetic field with B = |B T |, the axion-photon conversion probability is given by the well-known formula [30][31][32][33]: where x is the distance traveled by the photon, and is the effective photon mass squared in the presence of an ionized plasma with an electron number density n e . 1 The photons associated with typical observables travel through various environments, such as the IGM or the ICM, traversing a large number of magnetic domains. In order to quantitatively describe this phenomenon, some simplifying assumptions are made about the configuration of the magnetic fields in these environments and about the path traveled by the photons. We adopt the simple cell magnetic field model, first introduced in [33] and further developed in [25,36]. In this model the magnetic field is assumed to be split into domains (cells) in which it can be taken to be homogeneous. The photon path, extending from a source at some distance y to the observer, is assumed to cross a large number N of these magnetic domains. Each i-th domain has a physical size L i and a randomly oriented magnetic field of strength B i [36], whose component perpendicular to the photon's path is the same in each domain. With these simplifications, the resulting net probability of photon-axion conversion over many domains is then given by [25] depends on the ratio of the initial intensities of axions and photons coming from the source, denoted by I 0 a and I 0 γ respectively; and P 0,i is the conversion probability in the i-th magnetic domain, which can be obtained from Eq. (2.1) for x = L i .
Since N is very large, Eq. (2.4) can be rewritten as an integral. In order to do this, we further assume that y is a distance that scales linearly with N , such that s = y/N remains constant as N goes to infinity. For example, for IGM propagation the domains are typically assumed to be evenly distributed in comoving space, which means that each domain has comoving size s and the distance to the source is a comoving distance y = N s. Under these assumptions, we have The ratio of the observed photon flux and the emitted photon flux from the source is then given by P γγ = 1 − P aγ . (2.6)

Intergalactic medium propagation
We will consider the propagation of photons in different media. In this section, we focus on the IGM first. The IGM, more precisely the space between large scale structures, could be home to primordial magnetic fields, which serve as "seeds" for the observed magnetic fields in astronomical sources of different sizes, from stars to galaxy clusters. They could be generated during the preheating/reheating epochs immediately after inflation or during cosmological phase transitions before the formation of CMB. Magnetic fields produced at late times (at redshifts z < 10) from outflows of already formed galaxies could also reside in IGM. For a review of the generation mechanisms of IGM magnetic fields, see [37]. At the moment, there is no direct evidence of the IGM magnetic field. Instead there are observational upper and lower bounds on the amplitude of the magnetic field in IGM. CMB anisotropies set upper limits about nG on the present value of primordial magnetic field [38][39][40][41]. Other methods, such as the non-observation of Faraday rotation of the polarization plane of radio emission from distant quasars, set a similar upper limit [37]. 2 On the other hand, the non-observation of very high energy γ-ray cascade emission sets a lower bound on the magnetic field B IGM 10 −16 G for a coherent length above Mpc and becomes more stringent at smaller coherent lengths. For recent reviews of the constraints, see [37,43,44]. In our paper, we will adopt the cell magnetic field model in comoving space for the IGM. More concretely, we take 1 nG as a convenient benchmark value for the component of the comoving magnetic field perpendicular to the line of sight, as well as comoving coherent length (s IGM ) benchmarks of 0.1 Mpc, 1 Mpc, and 10 Mpc for the domain sizes. The cell magnetic field model is a very simple approximation to the structure of astrophysical magnetic fields. However, it has been shown to give the same results for axion-photon conversion as other more refined methods (such as power spectrum models) at high photon frequencies, while underestimating the conversion probability at lower frequencies [25,45,46]. The cell magnetic field model thus leads to conservative bounds on the axion parameter space. Finally, we caution the reader that the bounds we derive on axion coupling from datasets that are sensitive to the photon propagation in IGM should be understood as an upper bound on g aγγ × B IGM 1 nG for a fixed coherent length. We leave the discussion on more general combinations of B IGM and s IGM to Sec. 4.
Another important quantity of IGM that matters in our analysis is the electron density n e , which determines the plasma photon mass. At low redshifts, most of the baryons are in photoionized diffuse intergalactic gas (Lyman-α forest) and warm-hot intergalactic matter [47]. Among these two structures, Lyman-α forest contributes 28 ± 11% of the total mass (at z < 0.5) [47] but occupies 90% of the total volume [48]. The other structures, including warm-hot intergalactic matter, are more condensed and take up a much smaller volume. Thus, what matters more for the photon propagation is the Lyman-α forest. The average electron density of Lyman-α forest is about 6.5 × 10 −8 cm −3 , assuming its mass fraction to be the central value 28%. This, however, is not the entire story. Recent simulations show that for diffuse gas, most of the volume is occupied by cosmic voids (large under-dense patches) and sheets (two-dimensional structures of matter), which constitute approximately ∼ 30% and ∼ 40% of the entire volume at z 2 respectively [48]. Their mass fractions are significantly smaller, however, 8% and 20% respectively [48]. Based on this, the electron density of the sheet component of the Lyman-α forest is about half of the average one over all components, 3 × 10 −8 cm −3 ; while the electron density of the void component is about 1/4 of the average, 1.6 × 10 −8 cm −3 at z = 0. We will take these two values as benchmarks of the plasma electron density in IGM in our analysis.
For photons traveling through the IGM, Eqs. (2.5) and (2.6) could be rewritten in terms of z as The benchmark values of B IGM , s IGM , and n e,IGM are discussed and explained above.

Intracluster medium propagation
The angular diameter distances to galaxy clusters, as we will see in Sec. 2.4, rely on measurements of cluster X-ray brightness. The X-ray photons are produced throughout the cluster via Bremsstrahlung and line-emission involving the ionized plasma composing the ICM. These photons travel first through the ICM and then the IGM to reach the detector. Faraday rotation measurements in long wavelengths have shown [49][50][51][52] that ICM has magnetic fields with a strength of order O(µG). Therefore a fraction of the X-ray photons could convert into axions. This possibility has been studied in the literature and yields some of the strongest limits on couplings of very low mass axions to photons [53][54][55][56][57]. We devote the rest of this section to the computation of the effect ICM propagation has on X-rays photons as they leave the cluster. In order to perform this calculation, we need prescriptions for the ICM's electron number density n e,ICM and magnetic field B ICM .
We model n e,ICM with the double-β profile [24,58] n e,ICM (r) = n e,0 where n e,0 is the central density, r c1 , r c2 are the two core radii, f is the fractional contribution from the inner core, and β is the slope. Eq. (2.8) allows us to compute the photon plasma mass m γ , necessary for the calculation of the axion-photon conversion probability in Eq. (2.1). The values of the parameters for the double-β profiles of the clusters we use in this work can be found in [24].
For the magnetic field we follow previous literature [51,52,57,59] and assume the magnetic field follows a power law on the number density: where r ref is some reference radius from the cluster's center, B ref is the magnetic field value at that point, and η some power. We will take the two models of the ICM magnetic field of the Perseus cluster found in [57] and the one for the magnetic field of the Coma cluster in [51] as benchmarks for our analysis of the ICM effect: At small radii in Model A, r < 10 kpc, the electron number density is underestimated and spherical symmetry is unjustified. Therefore, we exclude the photon-axion conversion in the region of r < 10 kpc, following the treatment in [57]. We take L ICM = 6.08 kpc to be the (uniform) size of the magnetic domains, which is the mean of the L −1.2 distribution between 3.5−10 kpc proposed in [57]. 3 We take the virial radius of the cluster to be R vir = 1.8 Mpc, that of the Perseus cluster. 4 In treating the orientation of the magnetic field, we assume B ref to be the magnetic field value in the transverse direction, perpendicular to the photon's propagation direction. If we take the direction of the magnetic field to distribute uniformly between 0 and π, it is equivalent to substitute , witĥ B the direction of the magnetic field andk that of the photon propogation. This will make the bound on g aγγ derived from the ICM effect about a factor of √ 2 weaker. This is verified numerically by assigning a random orientation in each domain. 5 For X-rays originating at a radius r in the cluster, we can then approximate the ratio of outgoing to initial X-ray photon flux after axion-photon conversion, following Eq. (2.4), as: 13) where N (r) = (R vir − r)/L ICM is the number of domains with size L ICM from origin point r to the virial radius of the cluster R vir ; P 0 (r i ) is the axion-photon probability conversion at the center r i of the i-th domain, given by Eq. Finally we want to comment on the uncertainties associated with the ICM magnetic fields. Recently, for example, just how coherent or turbulent these fields are has been the subject of some concern, and the stringent bounds on the axion-photon coupling relying on the ICM propagation in [57] has been questioned [60]. As we will discuss below in more detail, we sidestep this issue by computing and comparing bounds based on various settings that, either include ICM photon-axion conversion with one of the benchmarks in Eqs. (2.10)-(2.12), or ignore the ICM conversion altogether.

Effects on distance observables
In the absence of a significant initial axion flux from the sources, non-negligible axionphoton mixing results in dimming: the brightness of a distant source will be decreased by a factor of P γγ . Historically, this effect was used in an early attempt at explaining away the cosmological constant [33]. While the cosmological constant has since been vindicated, the observation (or lack thereof) of dimming of distant sources can be used to constrain the axion parameter space.
In this section, we discuss the impact of dimming via axion-photon mixing on luminosity distances to type Ia SN and on angular diameter distances to galaxy clusters.

Luminosity distances and distance moduli
The flux F from a source of luminosity L located at redshift z is given by: where P γγ accounts for the possible non-conservation of photon flux between the observer and the source, and the luminosity distance D L is: .
Flux measurements of distant sources such as SNIa, our primary concern in this section, are usually expressed in terms of the source's apparent magnitude m, which is conventionally written as: where M (F 10 ) is the absolute magnitude (flux) of the source, defined at a distance of 10 pc from it; and µ is called the distance modulus which, from Eq. (2.14), can be rewritten as: is the effective LD in the presence of axion-photon conversion, and we have taken P γγ to be 1 at a distance of 10 pc from the source.
Putting everything together and making explicit the dependence on the parameters θ in our analysis, the effective apparent magnitude of the SNIa located at redshift z is with D L (z; Ω Λ , H 0 ) given by Eq. (2.15) and P γγ (z; θ) by Eq. (2.7). We will take A in Eq. (2.7) to be 2/3 since the initial axion flux from SNIa is negligible [36]. Note that Ref. [36] didn't consider the possibility of resonant production of axions at SNIa, which we will entertain in App. A and show that it doesn't modify the conclusion. Finally, we want to comment on the energy dependence of the photons. The photons from the SNIa are in the optical band with ω ≈ 1 eV. If the source is observed in various frequencies, its magnitude or flux will in general undergo spectral distortion (also called chromaticity) as a result of the photon energy dependence of P γγ , which can in principle be used to further constrain the axion parameter space. In the parameter space we are interested in, however, this distortion is negligible. This can be estimated by comparing the oscillation probabilities for the B band (4.3 eV) and V band (3.4 eV) photons respectively. The achromaticity requirement from data [61] could be translated into a constraint of P B γγ − P V γγ 0.03 [62]. For illustrative purposes let us assume that the photons transverse N = 3000 magnetic domains and g aγγ = 10 −11 GeV −1 . The probability difference computed using Eq. (2.5) and Eq. (2.6) is presented in Fig. 1. One can see from the figure that even with this relatively large g aγγ (which is already excluded by SN1987a [63]), the monochromaticity requirement only constrains m a in a tiny range around 10 −14 eV. Therefore, we do not consider achromaticity in the SNIa observations from photon-axion conversion further in our analysis.

Angular diameter distances of galaxy clusters
The angular diameter distance D A (ADD) is defined as the ratio of an astrophysical object's physical size d to the arc θ that it subtends in the sky. It can be shown to be equal to: . (2.20) In general ADDs are unaffected by the axion-photon conversion since they do not rely on brightness measurements of any kind. This is the case, for example, for observations that measure the imprint of the comoving sound horizon on the galaxy two-point correlation function, such as the BAO measurements. However, microwave and X-ray surveys can be used to determine ADDs to galaxy clusters, which would then be impacted by axion-photon conversion. Indeed, it has been shown in [64] that ADDs to galaxy clusters can be obtained from measurements of the clusters' X-ray surface brightness S X , due to Bremsstrahlung and line-emission resulting from ion-electron collisions in the ICM; combined with observations of the brightness temperature decrement ∆T SZ from CMB photons undergoing inverse Compton scattering with the same ICM, the so-called Sunyaev-Zeldovich effect (SZE) 6 : As we have seen in Sec. 2.3, X-rays originate in the cluster's ICM and are affected by axion-photon conversion as they travel through the ICM magnetic fields. Eq. (2.13) 6 Note that any prior photon-axion conversion effects in IGM, as the photons travel from the surface of last scattering to the clusters, would leave the brightness temperature decrement ∆TSZ unaffected. This is because this quantity is concerned with the relative brightness of the CMB photons that pass through clusters compared to those that do not. However, the absolute value of CMB distortion effects due to photon-axion conversion from the surface of last scattering can be used to constrain the axion parameter space, see [65,66].
describes the ratio of outgoing to produced X-rays flux, for photons traveling radially outwards from their oirigin at a distance r from the cluster's center. Since clusters are extended objects in the sky, the change in the X-ray brightness S X is not exactly described by Eq. (2.13), there being photons whose path from the interior to the exterior of the cluster and from there to the observer is not radial. Nevertheless, we expect that a good proxy for the exact ICM effect across the surface of the cluster is to weigh P γγ by the integrand sourcing the brightness S X and average over the line of sight. Indeed, S X ∝ dl n 2 e,ICM Λ ee [64], where l is the line-of-sight variable, and Λ ee the cluster cooling rate, which scales like the square root of the ICM temperature, ∼ T 1/2 e . From [24], we see the ICM temperature is approximately constant throughout most clusters, except for only a handful of them where the temperature fluctuates at most by a factor of four with small angular variation, such as RX J1347.5-1145, Abell 1835, Abell 2204, and Abell 1914. The resulting factor of two change in Λ ee leads to a percentage level change in the weighted average of P γγ . Therefore, we approximate the suppression on the X-ray brightness S X due to the ICM effect with: where the integral is taken from some initial radius r ini . For the ICM magnetic field Model A we follow [57] and take r ini = 10 kpc, whereas for models B and C we take r ini = 0 kpc. The suppression described in Eq. (2.22) immediately implies that there is a fraction 1 − P ICM γγ of the initial X-ray flux that has converted into axions. This means that the ratio of axions to photons outside the cluster is given by:  .7), for photons of both CMB and X-ray energies propagating in the IGM, in order to finally arrive at the effective ADDs to clusters:

Data and methodology
Having explained described the effects that axion-photon conversion has on various cosmological observables, we devote this section to describing the datasets and methodology we have used for our model fits.

Datasets
For our analysis we consider data from the following experiments, which we will combine in different ways: • Pantheon: the Pantheon dataset [14], consisting of apparent magnitude measurements of 1048 SNIa; • Clusters: a set of ADDs measurements for 38 galaxy clusters [24]; • SH0ES: measurements of the absolute magnitudes of 19 SNIa by the SH0ES collaboration [15]; • TDCOSMO: the Hubble parameter as measured by the TDCOSMO collaboration using strong lensing [67]; • BAO: the measurements of the imprint of baryon acoustic oscillations in galaxy distributions [20][21][22]; • Planck: The value of the comoving sound horizon at the epoch of baryon drag, given by the Planck collaboration's observation of CMB anisotropies [17].
In the rest of this section, we will describe in more detail these datasets and provide their corresponding likelihoods, some of which are inspired by MontePython [68,69]. We will then use these likelihoods to constrain the axion parameter space.

Pantheon
As we have seen in the previous section, axion-photon conversion impacts those cosmic distance measurements that rely on the brightness of astrophysical sources. The observation of the brightness of SNIa is one such kind of measurement. The Pantheon dataset is the most up-to-date collection of apparent magnitude measurements for 1048 SNIa in the redshift range of 0.01 < z < 2.3 [14]. The corresponding likelihood we use is given by where C Pan is the Pantheon inverse covariance matrix; m Pan i is the Pantheon measurements for the apparent magnitudes of the SNIa located at redshift z i while m eff (z i ; θ, M ) is the corresponding theory prediction given by Eq. (2.18) in the axion-photon conversion model. We take M as a free parameter in our MCMC runs and fit together with the model parameters.
For the SNIa in the Pantheon set, we will take the energy of their optical photons to be ω = 1 eV, the IGM magnetic field B IGM = 1 nG, the comoving size of the magnetic fields s IGM ∈ {0.1, 1, 10} Mpc, and the IGM electron number density n e,IGM either 1.6 × 10 −8 cm −3 or 3.0 × 10 −8 cm −3 . All benchmarks are in accordance with the discussion in Sec. 2.2.

Cluster angular diameter distances
Measurements of angular diameter distances (ADDs) to galaxy clusters, inferred from SZE and X-ray cluster data, are also sensitive to axion-photon conversion in a manner described in Sec. 2.4, and can therefore be used to constrain the axion parameter space. In our present work we use the sample of 38 clusters from [24] as listed in their Table 2, which assumes spherically symmetric clusters in hydrostatic equilibrium. 7 We then construct a likelihood for the ADD measurements from this dataset taking into account the statistical and systematic uncertainties enumerated in Table 3 of [24], which we add in quadrature. The likelihood is given by: where D eff A (z i ; θ) is given by Eq. (2.24). The data in [24] provides not only the redshifts and ADDs to these clusters but also the means and uncertainties of the n e,0 , f , r c1,c2 , and β parameters for the double-β profile of Eq. (2.8) describing the ICM electron number density n e,ICM . We have computed the impact of the uncertainties in these parameters on the line of sight-averaged ICM survival probability of Eq. (2.22). We found that these uncertainties only lead to a variation below 4% on top of the result using the mean values. We therefore ignore these subdominant effects, and restrict ourselves to the mean values provided by [24].
For the factors in Eq. (2.24) related to IGM propagation we assume the same benchmark quantities as for the SNIa Pantheon dataset. For the factors dealing with the ICM effect, we use the three benchmark magnetic field models described in Eqs. (2.10)-(2.12). We take the CMB photons to have energy ω CMB = 2.4 × 10 −4 eV. We average the X-ray photon energy in the band 0.7 − 7 keV using the measured temperature of each cluster [24], and the resulting photon effective energy is around ω X = 5 keV, which we use for our fits. In light of the uncertainties in the axion-photon conversion for X-rays in the ICM discussed in Sec. 2.3, we also perform fits to the ADD data ignoring the ICM effect.

BAO
Galaxy surveys can determine the imprint of baryon acoustic oscillations on matter distribution and then ADDs at low redshifts. More concretely, they measure ratios of the comoving sound horizon at the epoch of baryon drag r drag s to either the comoving angular We use the recent observations of r drag s /D V at z = 0.106 by 6dFGS [20], of D V /r drag s at z = 0.15 by SDSS using the MGS galaxy sample [21], and of both D M /r drag s and r drag s /D H 7 The sphericity requirement is relaxed in the sample of 25 clusters studied in [23], where an elliptical morphology is assumed instead. In it, however, the values of the ΛCDM cosmological parameters are fixed, since they are highly degenerate with the shape parameters, whose determination is the main goal of the paper. Since we are interested in fitting the cosmological parameters along with those of the axion-photon system, we use the dataset in [24] instead. Note that [24] quantifies an error of 15% arising from the sphericity assumption. at z = 0.38, 0.51, and 0.61 by BOSS, from the CMASS and LOWZ galaxy samples of SDSS-III DR12 [22]. We use the covariance matrix to take care of the correlation between the three redshift bins from BOSS as the middle one completely overlaps with the other two. There is no correlation between 6dFGS, MGS sample, and BOSS since BOSS only contains data with z > 0.2.
Note that since none of these surveys rely on the brightness of sources, these measurements are insensitive to axion-photon conversion effects, and therefore can be used to constrain the cosmological parameters {H 0 , Ω Λ } of ΛCDM. Since these measurements depend on r drag s , whenever we use these datasets we include r drag s as an extra parameter to our model. Schematically, then, the BAO likelihood is given by: where C BAO is the inverse covariance matrix of the BAO measurements. Q BAO i is the quantity being measured at redshift z i , and Q ΛCDM (z i ; Ω Λ , H 0 , r drag s ) is the model's prediction, which depends only on the cosmological parameters {H 0 , Ω Λ , r drag s } and is therefore identical to that of ΛCDM.

SH0ES
The SH0ES collaboration used parallax to deduce the distances to standard candles such as Cepheid variables in order to determine the absolute magnitude of 19 accompanying SNIa [15]. We then construct the corresponding likelihood: (3.6) Note that we are using the SH0ES collaboration's determination of the absolute magnitude M and not their value for H 0 , since this was determined under the assumption of photon flux conservation, which is not true in the axion-photon conversion framework.

TDCOSMO
The Hubble parameter can be determined through strong lensing. A sample of 7 such lenses was used by the TDCOSMO collaboration to determine a value of H 0 = 74.5 +5.6 −6.1 km sec −1 Mpc −1 [67]. We note that this measurement is independent of photon brightness and thus constrains H 0 only, not the axion parameter space. The likelihood we use is therefore: where for simplicity we take the symmetrized error σ TD = 5.85 km sec −1 Mpc −1 .

Planck
The use of the BAO likelihood defined in Eq. (3.4) requires the introduction of the comoving sound horizon at baryon drag r drag s as an extra parameter in our model. There is enough constraining power in the late-Universe data from Pantheon+SH0ES+TDCOSMO to determine the value of this parameter. However a different possibility is to use early-Universe data from Planck's observations of the CMB anisotropies [17], which yield r drag,obs s = 147.09 ± 0.26 for TT,TE,EE+low-E+lensing measurements. The likelihood we use is then simply given by: (3.8) In the next section we describe how we deal with the so-called Hubble crisis and the discrepancies between SH0ES, TDCOSMO, and Planck.

Methodology
We consider various combinations of the datasets described in Sec. 3.1, as well as different assumptions regarding the ICM and IGM, with the goal of deriving and comparing bounds on the axion parameter space (m a , g aγγ ) in different cases: • Early vs. Late: In light of the Hubble crisis and the disagreement regarding H 0 and r drag s between the SH0ES collaboration on the one hand and the Planck collaboration on the other [18,19], we split our datasets into two subsets with likelihoods given by L early ≡ L Pan · L cl · L BAO · L Pl , (3.9) L late ≡ L Pan · L cl · L BAO · L SH0ES · L TD , (3.10) and we fit to each likelihood separately.
• with vs. without ICM propagation: Given the debate surrounding the robustness of bounds on axion-photon interactions obtained from ICM propagation effects [57,60], we perform fits both with and without this effect. For the analyses that include the ICM effect, we assume three different models for the ICM magnetic field: A, B, and C, given by Eqs. (2.10)-(2.12).
• IGM magnetic domain sizes: We take three benchmarks for the coherent length of the IGM magnetic domains: s IGM ∈ {0.1, 1, 10} Mpc, described in Sec. 2.2.
We employ the likelihood-ratio test in order to find the 95% confidence level (C.L.) one-sided upper limits on the axion parameter space, putting constraints on the axion's coupling to photons at fixed masses. We will describe the procedure below, following Ref. [70]. We first take the total likelihood L tot (Θ) from one of the options described in the itemized cases above, based on either Eq. (3.9) or Eq. (3.10), for each of the ICM and IGM assumptions listed above. L tot (Θ) is a function of the combined set of cosmological and axion parameters Θ = {H 0 , Ω Λ , M, r drag s , m a , g aγγ }. We then scan over the parameter space of ξ ≡ {H 0 , Ω Λ , M, r drag s , g aγγ } to find the maximal likelihood at a fixed axion mass m a , max ξ L tot (m a ) ≡ L max (m a ). We allow g aγγ to take imaginary values as well, in order to cover the case with negative signal strength g 2 aγγ < 0. We then follow the prescription in Ref. [70] to construct the test statistic in Eq. (3.11). Instead of taking an evenly spaced grid, we make use of the public MCMC code emcee [71] to speed up the maximization of the likelihood. After that, at each point in the {m a , g aγγ } axion parameter space, we scan over the set of cosmological parameters only α ≡ {H 0 , Ω Λ , M, r drag s } to find max α L tot ≡ L max (m a , g aγγ ). Then we compute the log-likelihood ratio, We then set ∆χ 2 tot = 2.71 (the value for 95% C.L. one-sided upper limits with one degree of freedom [70]) to find the 95% C.L. upper bound. The range of parameters we scan is: 0.6 < Ω Λ < 0.75, 60 < H 0 /(km sec −1 Mpc −1 ) < 80, −21 < M < −18, 120 Mpc < r drag s < 160 Mpc, 10 −17 < (m a /eV) < 10 −11 , 10 −18 < (g aγγ GeV) < 10 −8 . This range is motivated by physical considerations (e.g., g −1 aγγ < 10 18 GeV, the Planck scale) 8 and sufficiently broad to include the points that maximize the likelihoods.

Results
In this section, we will present our results, based on the datasets and methodology described in the previous section. Our python code, which implements both the physics of axionphoton conversion and the MCMC analysis of its parameters using emcee, is publicly available at github.com/ManuelBuenAbad/cosmo axions.
First, we vary both m a and g aγγ to find the best fits for a given dataset. We find no evidence for axion-photon conversions, for any combinations of datasets, with the IGM and ICM assumptions considered and listed in the previous section. In order to investigate whether the axion-photon conversion and its impact on physical observables has been mismodelled, or whether there are unaccounted systematics in our datasets, we need to compare the best fit χ 2 values with non-negative signal strengths (i.e. real g aγγ ) to those where negative signal strengths (i.e. imaginary g aγγ ) are allowed. When only non-negative signal strengths are considered, we find that the best fits are consistent with the null hypothesis, or equivalently, ΛCDM with no axion-photon conversion. Now we allow negative signal strengths. For our analyses involving L early and no ICM conversion, g aγγ 's at best fits take imaginary values with the absolute values in the range between 10 −12 -10 −11 GeV −1 , and a |∆χ 2 | difference from that of the null hypothesis in the range 3.2-3.8. For cases with L early and ICM conversion, g aγγ 's at best fits are around imaginary 10 −12 GeV −1 , with |∆χ 2 | 2.2 with respect to the null hypothesis. On the other hand, fitting with L late yields |g aγγ | ∼ 10 −13 -10 −11.5 GeV −1 (some real, some imaginary) with |∆χ 2 | < 1, independent of whether ICM conversion is implemented or not. In most of these cases, the |∆χ 2 | difference between the best fits allowing only non-negative signal strengths and those allowing negative signal strengths is below 1σ. In all of the cases, this difference lies below the 2σ level (the two-sided χ 2 thresholds for two degrees of freedom are 2.28 and 5.99 for 68% and 95% respectively). This indicates both a consistent modeling of the axion-photon conversion and a correct accounting of the uncertainties in the dataset.
We then use the log-likelihood ratio in Eq. (3.11) to obtain 95% C.L. upper bound in the (m a , g aγγ ) plane assuming B IGM = 1 nG from both L early and L late , which is shown in Fig. 2. We also show constraints from varying the electron density in IGM and models to describe possible ICM effects on D A to galaxy clusters as described in Sec. 2.3 and Sec. 2.4.2. From all the numerical results, we learn that • The results are very similar for both L early and L late , for a given ICM model. The datasets used in L early and L late mainly differ in H 0 and r drag s anchors for BAO. Yet the constraints on g aγγ are mainly due to the shape of H(z) at late times constructed from various distance measurements, which are used for both L early and L late . In other words, our bounds do not depend on the resolution of Hubble crisis.
• The bounds are insensitive to the precise values of n e,IGM .
• We assume B IGM = 1 nG. For our results without ICM effects (red curves in all the figures of this section), the bounds should be understood as being constraints on g aγγ × B IGM 1 nG , for fixed values of the IGM coherent length s IGM . A better understanding of B IGM could help improve these bounds. For small axion masses, the conversion probability P 0 , within a single IGM magnetic domain, is in its linear regime, i.e. k s IGM in Eq. (2.1) and P 0 ∼ (g aγγ B IGM s IGM ) 2 . Eq. (2.5) in turn implies P aγ ∼ y s IGM (g aγγ B IGM ) 2 for sources at a comoving distance y. Consequently the bounds on g aγγ scale with ∝ 1/ √ s IGM for small axion masses.
On the other hand, for large axion masses P 0 oscillates rapidly within a single IGM domain and sin 2 (kx/2) averages to 1/2. Therefore P 0 ∼ (g aγγ B IGM ω) 2 /m 4 a and P aγ ∼ y s IGM (g aγγ B IGM ω) 2 /m 4 a . This means that the bounds on g aγγ scale with ∝ √ s IGM for large axion masses. The transition between both regimes occurs when • Including the ICM photon-axion conversion effects on the X-ray propagation used for inferring D A to galaxy clusters will make the upper limit stronger, for all the ICM magnetic field models we consider. In particular, for Model A in Eq. (2.10), the upper bound on g aγγ could be improved by one order of magnitude compared to the bound assuming no ICM effect. What is more, these ICM effects completely overshadow those from IGM propagation. Even choosing as a benchmark the smallest possible IGM magnetic field, B IGM = 10 −16 G, the constraints from ICM conversions remain the same, as strong as the ones presented in Fig. 2. In other words, constraints that include ICM effects are independent of the magnetic field strength of IGM.
• The galaxy cluster ADD measurements drive the likelihood-ratio, with subdominant contributions from the Pantheon dataset, both in the case where we ignore or include ICM X-ray conversion effects. The Pantheon SNIa dataset by itself does place constraints in the (m a , g aγγ ) parameter space, albeit somewhat weaker ones.
To further illustrate this last point, Fig. 4 shows the residuals for the Pantheon apparent magnitude (left panel) and cluster ADD (right panel) data, compared to a ΛCDM model with Ω Λ = 0.69, H 0 = 69 km sec −1 Mpc −1 , and M = −19.39. We also plot the effects of the axion-photon conversion on those observables, for n e,IGM = 1.6 × 10 −8 cm −3 , m a = 10 −16 eV, and both g aγγ = 6 × 10 −13 GeV −1 (orange) and g aγγ = 6 × 10 −12 GeV −1 (green). Note that in both panels the disagreement due to IGM conversion grows with redshift, as the effect gets stronger for the more distant sources. In the right panel, for the cluster ADD, we consider both cases where we only keep the IGM conversion and ignore ICM effects (lines), and where we include ICM effects with the model A from Eq. (2.10) for the ICM magnetic field (diamonds). Note that since each cluster has different parameters for its double-β profile, the ADD with ICM effects is different for each cluster. Also note that the presence of ICM conversion is the dominant contribution to the modification of the ADD distances to clusters, overshadowing the z-dependent IGM effect behind, and making the bounds independent of B IGM .  Lastly, we want to compare our results with existing studies in the literature, which is shown in Fig. 5. For readability, we only show the 95% C.L. upper limits from either assuming no ICM conversion effects on the galaxy cluster data or assuming model A in Eq. (2.10) for the effect. The upper limits for model B and C in Eqs. (2.11) and (2.12) are in between them. In the figure, we also show several other strong bounds on g aγγ in the same mass range from CAST [72], SN1987a [63] (note that [73] proposes a looser bound, due to an alternative modeling of the neutrino emission), X-ray searches from super star cluster [74] and X-ray spectroscopy from AGN NGC 1275 [57] (note that the ICM magnetic field modeling for NGC 1275 bound is questioned in [60]). We could see that, • the weakest limit we have, assuming that no X-ray photon-axion conversion in ICM, leads to a bound comparable to existing ones from SN1987a and super star cluster: g aγγ (4 − 5) × 10 −12 GeV −1 for m a 5 × 10 −13 eV, assuming B IGM = 1 nG and s IGM = 1 Mpc. For other IGM benchmarks, the bounds should be scaled by (nG/B IGM ) Mpc/s IGM accordingly for light axions.
• if the magnetic field in ICM is described by model A in Eq. (2.10), the strongest limit we have pushes g aγγ (5 − 6) × 10 −13 GeV −1 for m a 5 × 10 −12 eV. As mentioned above, these bounds are independent of B IGM and s IGM . Note that to avoid a busy plot, we do not show the bounds assuming model B and C in Eqs. (2.11) and (2.12). They are weaker than the one from model A but still stronger than the weakest limit assuming only IGM conversion.
Note that axions in the narrow mass range (6 × 10 −13 − 10 −11 ) eV are ruled out by superradiance of stellar black holes [75] and for even lighter axions with mass around or below 10 −20 eV, there exists interesting constraints on g aγγ from AGN [76], protoplanetary disk polarimetry [77] and CMB birefringence [78], which we do not show in the figure.
In addition, the distortion of CMB spectrum due to γ − a conversion only places strong bounds at m a > 10 −14 eV [65,66], which scales with B IGM .
It has been noted in [79] that for ultralight axions, cosmological considerations requiring axions to have a matter-power spectrum that matches that of cold dark matter constrains the magnitude of the axion couplings to the visible sector. As a result, at least part of the parameter space the cosmic distance measurements could probe is associated with non-trivial axion models, in which axions have an abnormally large coupling to photons, as constructed in [79][80][81][82].  When axion-photon conversion in ICM is neglected, which serves as a conservative benchmark to avoid the uncertainty of the magnetic field in ICM, we derive a bound comparable to existing bounds from SN1987a and super star clusters. These bounds are effectively constraints on (g aγγ × B IGM 1 nG × s IGM 1 Mpc ) for small axion masses, and therefore a direct measurement of B IGM and s IGM would fix exactly where this bound lies in the (m a , g aγγ ) parameter space. On the other hand, the inclusion of X-ray axion conversion in ICM makes the bound even stronger, no matter what ICM magnetic model we choose, and these bounds are entirely independent of the IGM parameters. In particular, model A of the magnetic field in ICM pushes the bound an order of magnitude stronger. A better understanding of the magnetic field in ICM could help reduce the uncertainties associated with its modeling. In addition, positive detection of axion-photon coupling from future experiments probing axion-photon coupling in this mass range [83][84][85][86] could help fix these bounds. Lastly, with future improvements in the precision of cosmic distance measurements, a better determination of late-time Hubble diagram H(z) is expected, which could further improve the sensitivity to possible departures from the ΛCDM prediction due to photon-axion conversion.

Acknowledgments
We thank David Pinner for early collaboration of this project. We thank Prateek Agrawal, Michael Geller, Dan Hooper, Matt Reece, Martin Schmaltz, Yu-Dai Tsai and Tomer Volansky for discussions at different stages of the project, as well as the anonymous referee, whose suggestions helped improve this paper. MBA

A Brightening Supernovae with axions and the Hubble crisis
In this appendix, we will discuss the interesting possibility of using axions to solve the Hubble crisis between early and late time measurements. This is not directly related to the main goal of our paper but has some similar ingredients, such as photon-axion conversion in IGM due to the magnetic fields. We first discuss some minimum requirements for this possibility and demonstrate why it does not work, at least for some minimal models. Ref. [19] also briefly discusses this possibility and comments on the potential observational challenges it faces, e.g., to explain other late-time datasets such as strong lensing [67]. We will provide a simple argument why this idea could not work even if we simply try to reconcile the SH0ES and Planck results.
The basic idea is that SNIa's further away on the cosmic distance ladder actually appear brighter than they would be in pure standard ΛCDM, because they also produce axions, which convert to photons en route and increase the net photon flux observed. Without correcting for the axion effects, the SNIa's further away will appear to be closer to us than they actually are. Thus the deduced D L 's of brightened SNIa's are shorter, resulting in a larger deduced H 0 , compared to its true value. More precisely, the effective luminosity distance D eff L from the observed flux of photons, F obs γ in SH0ES is given by where L SN is the luminosity of SNIa's. The Hubble value today measured by SH0ES is related to that inferred from Planck data as H SH0ES 0 = H Planck 0 (1 + ), where ∼ 10%. Therefore, assuming H Planck 0 is the true value of the Hubble rate today and to reconcile the late-time and early-time measurements, we need the observed photon flux to be enhanced by ∼20% compared to the flux without contribution from axions converting into photons. Using the formalism in Sec. 2, we have observed photon intensity from SNIa further away (e.g. at redshift z ∼ 0.1, or a distance of y ∼ 1 Gpc away), enhanced by a factor of about 1. Thus we need an initial axion flux I 0 a almost as large as the photon flux I 0 γ emitted by SNIa further away to solve the Hubble crisis in this scenario! This poses the first challenge to this potential solution. As shown in Ref. [36], the initial axion flux is negligible considering direct axion productions, namely, non-resonant conversions of photons in the SNIa's magnetic fields and in the magnetic fields of their host galaxies. One possibility that was ignored is the resonant conversion of photons to axions. In general, it is not easy to generate a large initial axion flux through resonant conversions, of which the general conditions required could be found in [87,88]. One necessary but not sufficient condition is to have a resonant shell in or near the SN, at which m a matches the plasma photon mass m γ . In a SNIa with about one solar mass and a radius of order 10 15 cm (the characteristic radius at 10 days when SNIa reaches its peak luminosity after the explosion of its progenitor white dwarf), the average electron density corresponds to a plasma photon mass ∼ 10 −5 eV. In the interstellar medium of the host galaxy outside SN, the plasma photon mass is of order 10 −11 eV. Thus, to have resonant conversions inside or near SN, the axion mass has to be m a 10 −11 eV.
On the other hand, for axion masses m a 10 −11 eV, the photon-axion conversion probability is negligible in IGM. For this axion mass range, the conversion probability in a single magnetic domain is approximately (A.4) The probability of axion-photon conversion remains tiny after photons/axions travel over 10 3 − 10 4 domains from a source Gpc away. This is consistent with our discussion in the main text. We only see a strong bound for m a 10 −13 eV, in which the photon-axion conversion in IGM becomes non-negligible.
In summary, to have axions brighten SNIa, we need resonant conversions inside or near SN in order to generate an initial axion flux as large as the initial photon flux. We also need more axions converting into photons in the IGM rather than the other way around. Yet as we show above by considering some simple necessary conditions for the scenario to work, the two requirements mentioned point towards very different axion mass ranges.
As bold model builders, we could consider more complicated scenarios, e.g, a photondark photon-axion system, similar to the setup in Ref. [89] for a different purpose. Then in the IGM, it is the dark magnetic field, which could be much larger than the ordinary magnetic field, that converts axions into photons or vice versa. Yet even considering a large dark magnetic field of order 10µG as considered in [89], we could see from Eq. (A.4) that the photon-axion conversion probability is still tiny for axion mass above 10 −11 eV. We will leave it for interested readers to explore further whether there are loopholes in our arguments.