#  Probing New Physics with High-Redshift Quasars:
Axions and Non-standard Cosmology

Chen Sun<EMAIL_ADDRESS>Theoretical Division, Los Alamos National
Laboratory, Los Alamos, NM 87545, USA  Manuel A. Buen-Abad<EMAIL_ADDRESS>Maryland Center for Fundamental Physics, Department of Physics, University of
Maryland, College Park, MD 20742, U.S.A.  Dual CP Institute of High Energy
Physics, C.P. 28045, Colima, México  JiJi Fan<EMAIL_ADDRESS>Department
of Physics, Brown University, Providence, RI, 02912, USA  Brown Theoretical
Physics Center, Brown University, Providence, RI, 02912, U.S.A.

###### Abstract

The Hubble diagram of quasars, as candidates to “standardizable” candles, has
been used to measure the expansion history of the Universe at late times, up
to very high redshifts ($z\sim 7$). It has been shown that this history, as
inferred from the quasar dataset, deviates at $\gtrsim 3\sigma$ level from the
concordance ($\Lambda$CDM) cosmology model preferred by the cosmic microwave
background (CMB) and other datasets. In this article, we investigate whether
new physics beyond $\Lambda$CDM (${\mathrm{B}\Lambda\mathrm{CDM}}$) or beyond
the Standard Model (BSM) could make the quasar data consistent with the
concordance model. We first show that an effective redshift-dependent relation
between the quasar UV and X-ray luminosities, complementing previous
phenomenological work in the literature, can potentially remedy the
discrepancy. Such a redshift dependence can be realized in a BSM model with
axion-photon conversion in the intergalactic medium (IGM), although the
preferred parameter space could be in mild tension with various other
astrophysical constraints on axions, depending on the specific assumptions
made regarding the IGM magnetic field. We briefly discuss a variation of the
axion model that could evade these astrophysical constraints. On the other
hand, we show that models beyond $\Lambda$CDM such as one with a varying dark
energy equation of state (${w\mathrm{CDM}}$) or the phenomenological
cosmographic model with a polynomial expansion of the luminosity distance,
cannot alleviate the tension. The code for our analysis, based on emcee
Foreman_Mackey_2013  and corner.py corner , is publicly available at
github.com/ChenSun-Phys/high_z_candles.

††preprint: LA-UR-23-29579

###### Contents

  1. I Introduction
  2. II Cosmic Distance Inference with Quasar Data Sets
    1. II.1 A Lightning Review
    2. II.2 How New Physics Biases the Distance Inference
  3. III An Effective Evolution of the UV–X-ray Relation
  4. IV Implications for B$\Lambda$CDM and BSM Models
    1. IV.1 Beyond $\Lambda$CDM
    2. IV.2 Beyond SM
    3. IV.3 Evaluation of the Fits
    4. IV.4 Comment on the Axion Model
    5. IV.5 Effective Evolution of $\beta$
  5. V Conclusion
  6. A Quasar Dataset
  7. B Methodology
  8. C Posterior of the Fits

##  I Introduction

Quasars, or quasi-stellar objects (QSOs), serve as probes in the ultraviolet
(UV) and infrared (IR) frequencies up to high redshifts ($z\sim 7$). With some
theoretical modeling of their intrinsic luminosities, they can be used to
measure luminosity distances as a function of redshift. Recently there has
been a renewed interest in using them as “standardizable” candles to measure
the expansion history of the Universe at late times Risaliti:2018reu ;
Lusso:2019akb ; Yang:2019vgk ; Velten:2019vwo ; Dainotti:2023ebr ;
Dainotti:2022rfz ; Lenart:2022nip ; Bargiacchi:2023jse . This is of particular
importance in light of the recent Hubble tension Riess:2016jrr ; Riess:2019cxk
; Aylor:2018drw ; Knox:2019rjx ; Abdalla:2022yfra . However, this hope of
having a new standard candle has been confronted by a series of challenges in
the form of various of consistency checks. In early attempts at its
application as a tool to constrain cosmological models, various groups have
confirmed that the quasar data prefers an expansion history that stands in
stark tension (at $\gtrsim 3\sigma$ level) with the one inferred from the
cosmic microwave background (CMB) Aghanim:2018eyx , which is in turn
consistent with type Ia supernovae (SNIa) Pan-STARRS1:2017jku , and Baryon
Acoustic Oscillation (BAO) BOSS:2016apd ; BOSS:2016hvq  measurements. This
discrepancy persists even after the application of stringent cuts to the data
that correct for various biases, such as those stemming from dust reddening,
X-ray absorption, and Eddington bias Lusso:2017hgz . On the other hand, purely
data-driven analyses have been implemented in Refs. Dainotti:2022rfz ;
Lenart:2022nip ; Bargiacchi:2023jse , which correct the luminosity-redshift
correlations in the UV and X-ray bands separately. These analyses show that
the UV–X-ray relation is indeed robust, and that it does not arise from any
luminosity-redshift correlation potentially caused by selection bias.
Furthermore they find that, for a given cosmology, the correlation-corrected
UV and X-ray luminosities deduced from observations each presents a different
redshift evolution, which has been previously unaccounted for. These analyses
remain agnostic as to the origin of said evolutions, and limit themselves to
characterize their size and impact on the quasar data.

In this paper we take a strategy of reverse-engineering: rather than claiming
that quasar data favors a cosmology in tension with that from the combined
CMB+SNIa+BAO datasets, we take it as input in our search for a plausible
explanation of the apparent redshift evolution of observed quasar fluxes, in
order to restore cosmic concordance. We define the concordance cosmology as
the ${\Lambda\mathrm{CDM}}$ model, whose late-time expansion history is fixed
by the dark energy density parameter $\Omega_{\Lambda}$ and the scaling factor
$h$ of the Hubble expansion rate $H_{0}=100h\,\mathrm{km/s/Mpc}$
Aghanim:2018eyx ,

$\displaystyle\Omega_{\Lambda}=0.6847\pm 0.0073~{},~{}~{}~{}h=0.6736\pm
0.0054~{}.$ (1)

While it is possible that the aforementioned luminosity evolution is of a
purely astrophysical nature (e.g., coming from a poor understanding of the
quasar luminosity itself, a possible temporal bias in quasar formation
history, or unaccounted-for propagation effects Dainotti:2022rfz ;
Lenart:2022nip ; Bargiacchi:2023jse ), in this paper we take a step further
and consider whether its origin could be in part coming from new physics. We
examine the implications that this evolution has for alternative models beyond
$\Lambda$CDM (${\mathrm{B}\Lambda\mathrm{CDM}}$) or beyond the Standard Model
(SM) of particle physics (BSM). Since the luminosity distance can be affected
by both the expansion history of the Universe and photon attenuation, we test
three benchmark models: $w$CDM and the cosmographic model Risaliti:2018reu ;
Lusso:2019akb  as examples of ${\mathrm{B}\Lambda\mathrm{CDM}}$, and axion-
like particles (ALPs, or axions for short) coupled to photons as an example of
BSM Csaki:2001yk ; Csaki:2001jk ; Buen-Abad:2020zbd . We show that while all
these alternatives provide good fits to the QSO dataset, both $w$CDM and the
cosmographic model do not resolve the tension between it and the CMB+SNIa+BAO
datasets. On the other hand, axion-photon couplings allow for frequency-
dependent photon disappearance, which greatly reduces the tension between all
the datasets; although it is still in tension with some other astrophysical
constraints.

This paper is organized as follows. We briefly review the UV–X-ray relation in
quasars and how that is affected by new physics in Sec. II. We next fit
${\Lambda\mathrm{CDM}}$ to the QSO data with a flexible UV–X-ray relation
parametrization in Sec. III, where we reveal the tension between it and other
datasets. We demonstrate that, when allowed to change, the quasar data prefers
an apparent redshift evolution in the UV–X-ray flux relation, which naturally
resolves the aforementioned tension. In Sec. IV we study the performances of
two ${\mathrm{B}\Lambda\mathrm{CDM}}$ models (${w\mathrm{CDM}}$ and
cosmographic model) and one BSM model (axion) in alleviating the tension
between the QSO and other datasets (SNIa, BAO, CMB). We show that the axion
model has the best performance in resolving the tension. We conclude in Sec.
V. The code we used in our numerical analysis, based on emcee
Foreman_Mackey_2013  and corner.py corner , is publicly available at
github.com/ChenSun-Phys/high_z_candles.

##  II Cosmic Distance Inference with Quasar Data Sets

In this section we describe the procedure through which QSO data can be used
to determine cosmic luminosity distances as a function of redshift, as well as
how new physics can change said procedure.

###  II.1 A Lightning Review

The observed power-law behavior of the X-ray spectrum of the quasars has
motivated the so-called “two-phase model” Haardt:1991tp ; Haardt:1994tx  of
the environment surrounding the supermassive black holes (SMBH), which are
believed to power the accretion of the active galactic nucleus (AGN). In the
two-phase model, ultraviolet (UV) photons are emitted by the relatively cold,
optically thick accretion disk around the SMBH; while X-ray photons are
produced through up-scattering of the UV photons by the hot, optically thin
corona of the SMBH.111There is a third component from the reflection of the
corona photon into the accretion disk peaks around 30 keVMartocchia:2017bou ;
Salvestrini:2019thn .

Since both the UV and X-ray photons are related to the SMBH mass and the
accretion rate, their respective luminosities can be formally written as
functions of these quantities, namely $L_{\rm X}=f_{1}(M_{\rm BH},\dot{M}_{\rm
BH})$, $L_{\rm UV}=f_{2}(M_{\rm BH},\dot{M}_{\rm BH})$ Lusso:2017hgz . This
implies a relation between $L_{\rm X}$ and $L_{\rm UV}$, as discussed widely
in the quasar/AGN literature; see for example Refs. zotero-21229 ;
Vignali:2003it ; Lusso:2009nq ; Steffen:2006px ; zotero-21239 . A common
parametrization of this relation is the so-called “Risaliti-Lusso (RL)
relation”, $L_{\rm X}=10^{\beta}\,L_{\rm UV}^{\gamma}$, with $L_{\rm X,UV}$
normalized to $\mathrm{erg/s/Hz}$. The parameter $\gamma$ has been confirmed
to have negligible redshift evolution by binning the flux measurement in
redshift Risaliti:2018reu ; Lusso:2017hgz . In addition, a particular model
Lusso:2017hgz , based on the two-phase model Haardt:1994tx , predicts a
redshift-independent $\beta$ as well.

The RL quasar luminosity relation can also be understood in terms of the flux
in the UV and X-ray band. Under the standard assumption that photon number is
conserved, the photon flux $F$ measured by an observer at a luminosity
distance $D_{L}$ can be related to the luminosity $L$ of the photon source as
follows:

$\displaystyle F(z;\omega)$ $\displaystyle=\frac{L(z;\omega)}{4\pi
D_{L}(z)^{2}}\ ,$ (2)

where $z$ is the redshift and $\omega$ the photon energy in question.222We
denote the photon flux $F(z;\omega_{x})$ at a specific energy $\omega_{x}$ by
$F_{x}(z)$.333Unless otherwise stated, we work in natural units. The RL
relation between the quasar UV and X-ray luminosities can then be written in
terms of the measured flux as follows:444We use $\log()$ and $\ln()$ to denote
the base-10 and natural logarithms, respectively.

$\displaystyle\log\left(\frac{F_{\rm
X}(z)}{\mathrm{erg/s/Hz/cm^{2}}}\right)-\gamma\log\left(\frac{F_{\rm
UV}(z)}{\mathrm{erg/s/Hz/cm^{2}}}\right)$ (3)
$\displaystyle=2(\gamma-1)\log\left({D_{L}(z)}/{\mathrm{cm}}\right)+\beta+(\gamma-1)\log(4\pi)\
,$

where $\gamma$ and $\beta$, which depend on quasar properties and dynamics,
are treated as nuisance parameters. The left hand side depends on the directly
observed quasar UV and X-ray fluxes, which means that it can be taken as
observational data, modulo the nuisance parameter $\gamma$. The right hand
side is a function of $D_{L}(z)$, which is a prediction of the underlying
cosmological model, and of the QSO nuisance parameters $\beta$ and $\gamma$.
Thus it can be treated as the theory input.

###  II.2 How New Physics Biases the Distance Inference

Both astrophysical processes and new physics could alter the flux relation in
Eq. (3) in a redshift-dependent way, which can be understood in terms of an
effective $\beta_{\mathrm{eff}}(z)$.

If the effective evolution in $\beta_{\rm eff}(z)$ has an astrophysical
origin, it must be corrected before the UV–X-ray relation can be reliably
utilized for cosmological inferences. This is indeed the subject of Refs.
Dainotti:2022rfz ; Lenart:2022nip ; Bargiacchi:2023jse , where the authors
eliminate the luminosity-redshift correlation assuming it takes the form of a
given function, and making use of the Efron-Petrosian method Efron:1992sch ,
before applying the UV–X-ray RL relation. As mentioned in the introduction,
these studies showed that the RL relation is not an artifact of any possible
selection bias effects, and that there is reason to believe that the UV and
X-ray quasar luminosities, as derived from observations and within the context
of a specific cosmic history, present each a different, non-negligible
redshift evolution. In these investigations, the luminosities in UV and X-ray
are assumed to be each corrected by a power-law form,555More sophisticated
functional forms were also tested, which showed no significant difference.
$L_{\rm X,UV}\rightarrow L_{\rm X,UV}(1+z)^{k_{\rm X,UV}}$, to account for the
apparent luminosity-redshift evolution. This effectively results in the
substitution of the constant $\beta$ with the redshift-dependent
$\beta_{\mathrm{eff}}(z)$:

$\displaystyle\beta\rightarrow\beta_{\rm eff}(z)=\beta+(k_{\rm X}-\gamma
k_{\rm UV})\log(1+z)~{}.$ (4)

However, simply correcting for the luminosity-redshift correlation using data-
driven statistical methods, while at the same time remaining agnostic as to
its origin, is not suited to our purposes. Indeed, any new physics that may be
the source of even just part of this correlation would not be picked up by
these methods, but would instead be discarded away during the correction
process along with any other less exotic sources. Therefore, we follow a
different approach by instead fitting the quasar dataset with different
$\beta_{\rm eff}(z)$ functional forms directly predicted by new physics
models.

We consider two classes of new physics models in which this
$\beta_{\mathrm{eff}}(z)$ modification of $\beta$ can arise. The first class
consists of deviations from the standard $\Lambda$CDM cosmological expansion
history. The measured flux (left-hand side of Eq. (3)) constraints the
expansion history of the Universe through $D_{L}(z)$. Anchoring the cosmology
to a given concordance history $D_{L,c}(z)$, $\beta_{\mathrm{eff}}$ can be
seen as a deviation from this concordance. Indeed, $\beta$ in Eq. (3) is
modified to

$\displaystyle\beta\rightarrow\beta_{\mathrm{eff}}(z)=\beta+2(\gamma-1)\log[D_{L}(z)/D_{L,c}(z)]\
.$ (5)

We would like to point out that such a modification of the expansion history
of the Universe would manifest itself in many other independent measurements
beyond those of quasars, including CMB anisotropy Aghanim:2018eyx , SNeIa Pan-
STARRS1:2017jku , and BAO BOSS:2016apd ; BOSS:2016hvq ; Beutler:2011hx . As we
show below, the quasar data prefers a non-standard cosmology at odds with that
favored by these other observations, which effectively limits this
interpretation of the tension.

The second class of models involves frequency-dependent propagation effects
that result in photon disappearance, which can also break the flux relation in
Eq. (3) without violating the relation in the luminosities $L_{\rm X}$ and
$L_{\rm UV}$:

$\displaystyle F(z;\omega)$
$\displaystyle=P_{\gamma\gamma}(z;\omega)\frac{L(z;\omega)}{4\pi
D_{L}(z)^{2}}\ ,$ (6)

where $P_{\gamma\gamma}(z;\omega)$ is the photon survival probability. This
effectively introduces a redshift dependence in $\beta$ as follows:

$\displaystyle\beta\rightarrow\beta_{\mathrm{eff}}(z)=\beta+\log\,P_{\rm\gamma\gamma}(z;\omega_{\rm
X})-\gamma\log\,P_{\rm\gamma\gamma}(z,\omega_{\rm UV})\ .$ (7)

In other words, if there is unaccounted-for extra attenuation in the measured
photon flux, a preference for non-zero en-route photon disappearance will
arise within the context of a given concordance cosmological expansion
history, such as $\Lambda$CDM; and it will manifest itself as a redshift
dependence of $\beta$ in the quasar data. In the scenario where this photon
disappearance is caused by new particles, the observed flux’s dependence on
$P_{\gamma\gamma}(z;\omega)$ effectively turns quasars into a probe of BSM
physics.

In this paper we take axion-photon interactions as a benchmark scenario of BSM
physics responsible for photon disappearance, and test for axion-photon
conversion in inter-galactic medium (IGM) using the QSO dataset. We will see
how this model, in conjunction with the standard $\Lambda$CDM, provides a very
good fit to the data. While this fit is in tension with some astrophysical
constraints on the axion-photon coupling, axion-photon interactions provide a
concrete physics model that both greatly improves the fit to the QSO data and
restores the consistency between the expansion history preferred by this data
and that favored by CMB, SNIa, and BAO. We also compare its goodness of fit
with a few ${\mathrm{B}\Lambda\mathrm{CDM}}$ models.

To summarize, our approach differs from and complements the analyses in Refs.
Dainotti:2022rfz ; Lenart:2022nip ; Bargiacchi:2023jse  in a few ways.
Firstly, Refs. Dainotti:2022rfz ; Lenart:2022nip ; Bargiacchi:2023jse
designed a data-driven method to empirically correct any luminosity-redshift
corrections. This results in the elimination of both any selection bias and
any new physics possibly contributing to this redshift evolution. Thus, this
method could potentially hide any new physics–induced flux evolution, should
there be any. Secondly, it is not our purpose to understand all luminosity-
redshift correlation present in the quasar data. By fitting the data with new
physics-induced $\beta_{\rm eff}$, we perform a critical examination of whethe
part of the flux-redshift evolution may be caused by new physics, if the
concordance cosmological history encoded in $\Lambda$CDM and Eq. (1) indeed
describe nature. Lastly, while Refs. Dainotti:2022rfz ; Lenart:2022nip ;
Bargiacchi:2023jse  are of critical importance in establishing the UV–X-ray RL
relation on firmer ground, the templates of the luminosity-redshift evolution
adopted therein only allow for smooth apparent luminosity changes across a
relatively large range of redshifts. As will become clear later in this paper,
the new physics model we test cannot be captured by the simpler functional
forms tested in Ref. Dainotti:2022rfz . In particular, the axion-induced flux
attenuation generates a somewhat abrupt change in both the UV and the X-ray
bands, more akin to a step function than to a slowly-varying transition.

##  III An Effective Evolution of the UV–X-ray Relation

In this section, we first discuss the tension between QSO dataset and
CMB+SNIa+BAO within the context of $\Lambda$CDM. We then show how an effective
$\beta_{\mathrm{eff}}(z)$ with redshift dependence is preferred by quasars,
once the cosmological expansion history is itself anchored by CMB+SNIa+BAO.
The datasets we use include:

  * •

$\mathcal{B}$: Our baseline datasets. These include SNIa: Pantheon Pan-
STARRS1:2017jku ; BAO: 6dFGS Beutler:2011hx , SDSS using the MGS galaxy sample
Ross:2014qpa , CMASS and LOWZ galaxy samples of SDSS-III DR12 Alam:2016hwk .

  * •

$\mathcal{G}$: Gaussian priors including $H_{0}=67.36\pm
0.54\,\mathrm{km/s/Mpc}$ and $r_{s}=147.09\pm 0.26\,\mathrm{Mpc}$ from Planck
2018 Aghanim:2018eyx .666We also tested the Gaussian prior of $H_{0}$ from
SH0ES, which is in significant tension with the Planck results. We find the
same results. This is expected as only the shape of the Hubble diagram matters
in our analysis.

  * •

$\mathcal{Q}$: QSO dataset from Ref. Lusso:2020pdb . We further parametrize
the nuisance $\beta$ parameter in three possible ways:

    * –

$\mathcal{Q}_{\beta_{0}}$: $\beta_{\mathrm{eff}}=\beta_{0}$ is a constant
across all redshifts. This is the most common parametrization in the quasar
literature.

    * –

${\mathcal{Q}_{\beta_{\sharp}}}$: $\beta_{\mathrm{eff}}$ is a step function.
At redshift $z_{0}$, it sharply transitions from a value $\beta_{0}$ to
$\beta_{1}$.

    * –

${\mathcal{Q}_{\beta_{\rm t}}}$: $\beta_{\mathrm{eff}}$ goes through a smooth
transition, parametrized with a $\tanh$ function as follows

$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\beta_{\mathrm{eff}}(z)=\beta_{0}+\frac{\beta_{0}-\beta_{1}}{2}\left[\tanh\left(\frac{z-z_{0}}{\delta
z}\right)+1\right].$ (8)

According to this parametrization there are five nuisance parameters in the
QSO dataset in total: $\gamma$, $\beta_{0}$, $\beta_{1}$, $z_{0}$, and $\delta
z$.

The details of the fits can be found in Appendices A and B. The 1D posteriors
can be found in Appendix C. We list the log-likelihood of the best fit points
and the posteriors of dark energy density parameter, $\Omega_{\Lambda}$, in
seven different runs with $\Lambda$CDM in Tab. 1.

|  QSO |  SNIa |  BAO |  $\Omega_{\Lambda}$
---|---|---|---|---
$\mathcal{B}+\mathcal{G}$ |  |  -1177.4 |  -9.2 |  $0.68\pm 0.02$
$\mathcal{Q}_{\beta_{0}}+\mathcal{G}$ |  -106.6 |  |  |  $0.05^{+0.07}_{-0.03}$
${\mathcal{Q}_{\beta_{\sharp}}}+\mathcal{G}$ |  -169.2 |  |  |  $0.10^{+0.12}_{-0.07}$
${\mathcal{Q}_{\beta_{\rm t}}}+\mathcal{G}$ |  -179.1 |  |  |  $0.46^{+0.26}_{-0.27}$
$\mathcal{Q}_{\beta_{0}}+\mathcal{B}+\mathcal{G}$ |  -66.62 |  -1174.9 |  -9.4 |  $0.66\pm 0.02$
${\mathcal{Q}_{\beta_{\sharp}}}+\mathcal{B}+\mathcal{G}$ |  -149.2 |  -1177.0 |  -9.4 |  $0.67\pm 0.02$
${\mathcal{Q}_{\beta_{\rm t}}}+\mathcal{B}+\mathcal{G}$ |  -179.1 |  -1177.3 |  -8.9 |  $0.68\pm 0.02$
Table 1:  The log-likelihood of the best fit points in seven different runs
with $\Lambda$CDM.

We make a few comments on the results below. First, the aforementioned tension
between the QSO and SNIa+BAO datasets can be easily seen in the posterior of
$\Omega_{\Lambda}$ from the $\Lambda$CDM fit to $\mathcal{B}+\mathcal{G}$
($\Omega_{\Lambda}=0.68\pm 0.02$) and to $\mathcal{Q}_{\beta_{0}}+\mathcal{G}$
($\Omega_{\Lambda}=0.05^{+0.07}_{-0.03}$). Although
$\mathcal{Q}_{\beta_{0}}+\mathcal{B}+\mathcal{G}$ leads to a posterior in
$\Omega_{\Lambda}$ that is consistent with the concordance value, this is
mostly due to the addition of SNIa+BAO. This can be seen by comparing the log-
likelihood between $\mathcal{Q}_{\beta_{0}}+\mathcal{B}+\mathcal{G}$ and
$\mathcal{Q}_{\beta_{0}}+\mathcal{G}$. The fit to quasars is degraded
significantly once the baseline data combination SNIa and BAO is added, with a
change of $\Delta\chi^{2}=+40.0$ for the QSO dataset alone. What is more, the
fit to SNIa within the $\mathcal{Q}_{\beta_{0}}+\mathcal{B}+\mathcal{G}$ run
is also slightly worse than that within $\mathcal{B}+\mathcal{G}$,
$\Delta\chi^{2}=+2.5$. This is due to quasars pulling the parameters away from
the minimum of SNIa. Put together, these facts are indications of the
incompatibility between the QSO dataset assuming a constant $\beta$ and
SNIa+BAO.

Second, the tension is greatly relieved when considering a redshift-dependent
$\beta_{\mathrm{eff}}(z)$:

  * •

In the ${\mathcal{Q}_{\beta_{\sharp}}}+\mathcal{G}$ run, there is a clear
preferred transition point at $z_{0}=1.65^{+0.01}_{-0.01}$. While the
posterior in $\Omega_{\Lambda}$ is still largely off compared with that from
$\mathcal{B}+\mathcal{G}$, there is huge improvement in the fit to the QSO
dataset with $\Delta\chi^{2}=-62.6$, at the cost of only two more parameters.
Furthermore, by comparing
${\mathcal{Q}_{\beta_{\sharp}}}+\mathcal{B}+\mathcal{G}$ and
${\mathcal{Q}_{\beta_{\sharp}}}+\mathcal{G}$, we see the former has only a
milder deterioration in the log-likelihood of quasars, $\Delta\chi^{2}=20.0$,
while the value of $\Omega_{\Lambda}$ is consistent with that from
$\mathcal{B}+\mathcal{G}$.

  * •

With ${\mathcal{Q}_{\beta_{\rm t}}}$, the tension is completely resolved.
Comparing ${\mathcal{Q}_{\beta_{\rm t}}}+\mathcal{G}$ with
$\mathcal{Q}_{\beta_{0}}+\mathcal{G}$, we see an improvement in the fit to the
quasar data with $\Delta\chi^{2}=-72.5$ between the two fits. In addition, the
fit to quasars is as good in ${\mathcal{Q}_{\beta_{\rm
t}}}+\mathcal{B}+\mathcal{G}$ as in ${\mathcal{Q}_{\beta_{\rm
t}}}+\mathcal{G}$. More importantly, ${\mathcal{Q}_{\beta_{\rm
t}}}+\mathcal{G}$ leads to a value of $\Omega_{\Lambda}$ that is compatible
with the concordance cosmology. We show the corresponding shape of
$\beta_{\mathrm{eff}}(z)$ in Sec. IV.5.

##  IV Implications for B$\Lambda$CDM and BSM Models

In this section, we discuss several specific ${\mathrm{B}\Lambda\mathrm{CDM}}$
and BSM models that realize an effective redshift-dependent
$\beta_{\mathrm{eff}}(z)$, and whether they could alleviate the tension
between the QSO dataset and other cosmological datasets.

###  IV.1 Beyond $\Lambda$CDM

The luminosity distance can be computed once a cosmological expansion history
is given:

$\displaystyle D_{L}(z)$
$\displaystyle=(1+z)\,\int_{0}^{z}\,dz^{\prime}\frac{c}{H(z^{\prime})}~{},$
(9)

where $c$ is the speed of light. This affects the UV–X-ray relation in the
measured fluxes as shown in Eq. (3), which causes an effective modification of
the QSO nuisance parameter $\beta$ given in Eq. (5).

First we consider the case of $w$CDM, where the equation of state of the dark
energy, $w$, deviates from $-1$. A phenomenological parametrization commonly
used in the literature is

$\displaystyle w(a)$ $\displaystyle=w_{0}+w_{a}(1-a)~{},$ (10)

where $a$ is the scale factor. This function smoothly interpolates the current
equation of state (EOS) $w_{0}$ to its value in the early Universe
$(w_{0}+w_{a})$. The parameters $w_{0},w_{a}$ are allowed to vary between -1
and +1. The $w$CDM model has therefore four theory parameters to be fitted,
$\Omega_{\Lambda},H_{0},w_{0},w_{a}$.

We also follow Refs. Risaliti:2018reu ; Lusso:2019akb  and perform a free-form
polynomial expansion of the luminosity distance in order to study potential
alternative cosmologies in a more general way. If one fixes the low redshift
luminosity distance to $c/H_{0}$, the so-called cosmographic approach is an
effective expansion of the form:

$\displaystyle D_{L}(z)$
$\displaystyle=\frac{c}{H_{0}}\bigg{(}\ln(1+z)+\sum_{i=2}a_{i}\ln^{i}(1+z)\bigg{)}~{}.$
(11)

We truncate the expansion up to the forth order. Therefore, the model
parameters are $H_{0},a_{2},a_{3},a_{4}$. Note that Ref. Yang:2019vgk  showed
that the convergence of the series is poor when they are mapped to physical
models ($\Lambda$CDM, $w$CDM, etc.). Therefore, we use this only as a
benchmark phenomenological parametrization to compare with the literature, and
refrain from making any statements about the underlying cosmological model.

###  IV.2 Beyond SM

We now devote our attention to a concrete particle physics model that can
leave imprints in the UV–X-ray relation of Eq. (3). If axions exist and couple
to photons, the number of photons may not be conserved when they propagate
through a static magnetic field background, such as the IGM magnetic field.
The photon disappearance probability $P_{0}$ (within a single magnetic domain)
due to an axion $a$ with a mass $m_{a}$ and a coupling to photons
$g_{a\gamma}aF\tilde{F}/4$, where $g_{a\gamma}$ is the constant coupling with
energy dimension $-1$, and $F$ ($\tilde{F}$) is the (dual) electromagnetic
field strength, is given by the well-known formula Georgi:1983sy ;
Sikivie:1983ip ; Raffelt:1987im ; Csaki:2001yk :

$\displaystyle
P_{0}=\frac{(2\Delta_{a\gamma})^{2}}{k^{2}}\sin^{2}\left(\frac{kx}{2}\right)\
,$ (12)

where $x$ is the distance traveled by the photon, and

$\displaystyle k$
$\displaystyle\equiv\sqrt{(2\Delta_{a\gamma})^{2}+\left(\Delta_{a}-\Delta_{\gamma}\right)^{2}}\
,$ (13) $\displaystyle\Delta_{a\gamma}$
$\displaystyle\equiv\frac{g_{a\gamma}B}{2}\
,\quad\Delta_{a}\equiv\frac{m_{a}^{2}}{2\omega}\
,\quad\Delta_{\gamma}\equiv\frac{m_{\gamma}^{2}}{2\omega}\ ,$ (14)

with $B$ the IGM magnetic field transverse to the photon trajectory, and
$\omega$ the photon energy. Here $m_{\gamma}^{2}\equiv\frac{4\pi\alpha
n_{e}}{m_{e}}$, where $m_{e}$ and $\alpha$ are the electron mass and fine-
structure constant respectively, is the effective photon mass squared in the
presence of an ionized plasma with an electron number density $n_{e}$. In our
analysis, we take $n_{e}\simeq 1.6\times 10^{-8}\,\mathrm{cm}^{-3}$
Nicastro:2018eam ; Martizzi:2018iik . There are therefore a total of four
theory parameters to be fitted $\Omega_{\Lambda},H_{0},m_{a},g_{a\gamma}$.

We adopt the cell model for the IGM magnetic field, described in Refs.
Csaki:2001yk ; Csaki:2001jk ; Grossman:2002by ; Buen-Abad:2020zbd . 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}$
Grossman:2002by , whose component perpendicular to the photon’s path is
assumed to be the same in each domain. With these simplifications, the
resulting net probability of photon-axion conversion over many domains is then
given by Avgoustidis:2010ju

$\displaystyle
P_{a\gamma}(y)=(1-A)\left(1-\prod\limits_{i=1}^{N}\left(1-\frac{3}{2}P_{0,i}\right)\right)\
,$ (15)

where $A\equiv\frac{2}{3}\big{(}1+\frac{I_{a}^{0}}{I_{\gamma}^{0}}\big{)}$
depends on the ratio of the initial intensities of axions and photons coming
from the source, denoted by $I_{a}^{0}$ and $I_{\gamma}^{0}$ respectively; and
$P_{0,i}$ is the conversion probability in the i-th magnetic domain, which can
be obtained from Eq. (12) for $x=L_{i}$. Since $N$ is very large, Eq. (15) 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 $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=Ns$. Under these assumptions, we have

$\displaystyle
P_{a\gamma}(y)=(1-A)\left(1-\exp\left[\frac{1}{s}\int\limits_{0}^{y}\mathrm{d}y^{\prime}~{}\ln\left(1-\frac{3}{2}P_{0}(y^{\prime})\right)\right]\right)\
.$ (16)

We assume $A=2/3$ or equivalently $I_{a}^{0}=0$ throughout the paper. The
ratio of the observed photon flux and the emitted photon flux from the source
is then given by $P_{\gamma\gamma}=1-P_{a\gamma}$.

Let us now briefly comment on the coherent length of the magnetic field
domain. While the existence of IGM magnetic fields has been indirectly
attested (from the absence of $\gamma$-ray cascade emission) and there are
upper bounds as well (from CMB anisotropy as well as turbulence decay),
observational constraints on the coherent length are lacking Durrer:2013pga .
From the theoretical perspective, the coherent length of a magnetic field is
related to the structures that support the magnetic field. Therefore, it is
reasonable to believe that the magnetic fields at the intersections of
filaments has a domain size of $\sim$Mpc scale. Nevertheless, the IGM magnetic
field domains can be much larger due to characteristic size of other
structures in the IGM. Because voids and sheets make up the majority of the
IGM volume, it is natural to suppose that the relevant length scale associated
with their corresponding magnetic fields should be related to the size of
these structures. From the IllustrisTNG simulationMartizzi:2018iik , one can
see that the void size at low redshift ($z\sim 1$) can reach $\mathcal{O}(10)$
Mpc. In Sec. IV.4, we fit the axion model to our datasets with different
choices of the domain size. Elsewhere, and unless stated otherwise, we use as
a benchmark a magnetic domain size of $1~{}\mathrm{Mpc}$ and a magnitude field
of $B=1\,\mathrm{nG}$ throughout this paper.

At last, we stress that the impact of axion-photon interactions only manifests
itself in distance measurements that involve photon flux, such as the
luminosity distances of quasars. It is otherwise largely invisible in
experiments based on direct measurements of the angular diameter sizes (e.g.
CMB observations.) In addition, for the axion-photon couplings in which we are
interested, the impact on SNIa luminosity distances (sensitive to only
$\mathrm{eV}$-energy photons) is negligible, most effects being limited to the
UV photons at $z>1.5$.

###  IV.3 Evaluation of the Fits

We fit each model ($w$CDM, cosmographic, and axion models) to either quasars
alone $(\mathcal{Q}_{\beta_{0}}+\mathcal{G})$ or to quasars and our baseline
datasets ($\mathcal{Q}_{\beta_{0}}+\mathcal{B}+\mathcal{G}$), including in
both cases the Gaussian prior ($\mathcal{G}$) on $H_{0}$ and $r_{s}$. See
Appendices A and B for details on our methodology, including the QSO dataset,
the priors, and the fitting procedure. We summarize the resulting log-
likelihoods of the best fit points in Tab. 2. The 1D and 2D posteriors are
shown in Appendix C.

|  QSO |  SNIa |  BAO |  $\chi^{2}_{c}$
---|---|---|---|---
${w\mathrm{CDM}}(\mathcal{Q}_{\beta_{0}}+\mathcal{G})$ |  -171.1 |  |  |  1247.5
${w\mathrm{CDM}}(\mathcal{Q}_{\beta_{0}}+\mathcal{B}+\mathcal{G})$ |  -91.7 |  -1172.0 |  -9.1 |  743.0
${a_{2}a_{3}a_{4}}(\mathcal{Q}_{\beta_{0}}+\mathcal{G})$ |  -179.4 |  |  |  996.9
${a_{2}a_{3}a_{4}}(\mathcal{Q}_{\beta_{0}}+\mathcal{B}+\mathcal{G})$ |  -140.5 |  -1167.7 |  -8.9 |  822.5
$\mathrm{axion}(\mathcal{Q}_{\beta_{0}}+\mathcal{G})$ |  -169.7 |  |  |  609.0
$\mathrm{axion}(\mathcal{Q}_{\beta_{0}}+\mathcal{B}+\mathcal{G})$ |  -147.8 |  -1176.5 |  -9.3 |  514.2
Table 2:  The log-likelihood of the best fit points in six tests with beyond
$\Lambda$CDM and beyond SM models. The cosmographic model is labeled as
${a_{2}a_{3}a_{4}}$. In the last column, we show the “theory distance” between
the best fit points and the concordance cosmological model by comparing their
$\chi^{2}$ difference from the quasars likelihood; see Eq. (24). The smaller
the value is, the closer the preferred cosmology is to the concordance model.

Tension between QSO and SNIa+BAO. We start by noting that the axion model
significantly reduces the tension between QSOs ($\mathcal{Q}_{\beta_{0}}$) and
SNIa+BAO ($\mathcal{B}$). In comparison, the stronger tension between
$\mathcal{Q}_{\beta_{0}}$ and $\mathcal{B}$ remains in $w$CDM and the
cosmographic model.

This is reflected in the goodness of fit to $\mathcal{Q}_{\beta_{0}}$ of each
model. When SNIa and BAO are included, the goodness of fit to QSOs is degraded
by $\Delta\chi^{2}=79.4$ for $w$CDM, by $\Delta\chi^{2}=38.9$ for the
cosmographic model, and by $\Delta\chi^{2}=21.9$ for the axion model. This
shows that both $w$CDM and the cosmographic model take a significant deviation
from the concordance cosmology over a large range of redshifts in order to
accommodate the QSO dataset, with such changes disfavored by SNIa+BAO. As a
result, including SNIa+BAO limits how close these models can reach the minimum
of the quasar likelihood.

On the other hand, the axion model allows a modification of the UV and optical
photon flux starting at $z\gtrsim 1.5$ while leaving the low-$z$ part mostly
unchanged. This minimizes the impact on SNIa, and BAO is not affected by it at
all. As a result, including SNIa+BAO only has a mild effect on the goodness of
fit to QSO data. Such “flexibility” of the axion model compared to
${w\mathrm{CDM}}$ and the cosmographic model is the key to reduce the tension
between the $\mathcal{Q}_{\beta_{0}}$ and $\mathcal{B}$ datasets. Meanwhile,
the flexibility of the axion model comes at no extra cost in terms of the
introduction of more parameters than the ${w\mathrm{CDM}}$ and cosmographic
models.

It is interesting to note that the fit to quasar in
${a_{2}a_{3}a_{4}}(\mathcal{Q}_{\beta_{0}}+\mathcal{G})$ and
$\mathrm{axion}(\mathcal{Q}_{\beta_{0}}+G)$ only differ by less than
$\Delta\chi^{2}=10$, with the ${a_{2}a_{3}a_{4}}$ model performing slightly
better. However, when $\mathcal{B}$ is included, the axion model is preferred
with $\Delta\chi^{2}\approx-7.3$, due to precisely the aforementioned
flexibility in the axion model that effectively decouples the impact on quasar
flux, SNIa flux, and BAO.

Restoration of the Concordance Cosmology. Resolving the tension means the
minimum of $\mathcal{Q}_{\beta_{0}}$ should get closer to that of
$\mathcal{B}$ in the theory space. However, this does not guarantee a
restoration of concordance cosmology. For example, the minima of the two
datasets can both drift toward a theory point that is far away from the
concordance model. Therefore, aside from the relative distance between the
minima of $\mathcal{Q}_{\beta_{0}}$ and $\mathcal{B}$, we would also like to
evaluate the distance between the minima and a fixed point anchored by the
concordance model, which is essentially the minimum of the Planck likelihood.

To quantify their agreement to the well-established concordance cosmology
given by Eq. (1), we compare the distance modulus computed from the best fit
points of the six runs
$\\{{a_{2}a_{3}a_{4}},{w\mathrm{CDM}},\mathrm{axion}\\}\otimes\\{\mathcal{Q}_{\beta_{0}}+\mathcal{G},\mathcal{Q}_{\beta_{0}}+\mathcal{B}+\mathcal{G}\\}$
with that computed from Eq. (1). More concretely, we define the concordance
distance modulus as

$\displaystyle\mu_{c}$
$\displaystyle=5\;\log_{10}\left(\frac{D_{L,c}}{10\,\mathrm{pc}}\right),$ (17)

where $D_{L,c}$ is the luminosity distance computed from Eq. (9) with the
Hubble constant taken from Eq. (1). The error associated with $\mu_{c}$ is
given by

$\displaystyle\Delta\mu_{c}=\left[\left(\frac{\partial\mu_{c}}{\partial\Omega_{\Lambda}}\right)^{2}\left(\Delta\Omega_{\Lambda}\right)^{2}+\left(\frac{\partial\mu_{c}}{\partial
h}\right)^{2}\left(\Delta h\right)^{2}\right]^{1/2}.$ (18)

The distance modulus at each best fit point of the six runs can be computed
with the $D_{L}$ inferred from the cosmological parameters
[$(a_{2},a_{3},a_{4})$ in the cosmographic models,
$(\Omega_{\Lambda},h,w_{0},w_{a})$ in ${w\mathrm{CDM}}$, and
$(\Omega_{\Lambda},h)$ in $\mathrm{axion}$]. Since the best fit point fits the
data quite well in each of the six tests, one can also use the measurement to
represent the best fit theory point for simplicity,

$\displaystyle\mu_{\rm bf}$
$\displaystyle=\frac{5}{2(\gamma-1)}\left[\log\left(\frac{f_{\rm X}}{f_{\rm
UV}^{\gamma}}\right)-\beta\right]$ (19)
$\displaystyle~{}~{}~{}~{}-5\,\log[(4\pi)^{1/2}]-5\,\log\left[{10\,\mathrm{pc}}/{\mathrm{cm}}\right],$
(20)

where $f_{\rm X,UV}$ are the X-ray and UV band fluxes normalized in units of
$\mathrm{erg/s/Hz/cm^{2}}$.777In the case of the axion model, they are
corrected for extra attenuation due to the axion-photon conversion: $f_{\rm
X,UV}=\frac{F_{\rm
X,UV}(z)/(\mathrm{erg/s/Hz/cm^{2}})}{P_{\gamma\gamma\rm~{}X,UV}(z)}$. The
error of $\mu_{\rm bf}$ is computed as follows.

$\displaystyle\Delta$ $\displaystyle\mu_{\rm
bf}=\bigg{[}\left(\frac{\partial\mu}{\partial\log F_{\rm
X}}\right)^{2}\Delta\log F_{\rm X}^{2}+\left(\frac{\partial\mu}{\partial\log
F_{\rm UV}}\right)^{2}\Delta\log F_{\rm UV}^{2}$ (21)
$\displaystyle+\left(\frac{\partial\mu}{\partial\log
m_{a}}\right)^{2}\Delta\log(m_{a})^{2}+\left(\frac{\partial\mu}{\partial\log
g_{a\gamma}}\right)^{2}\Delta\log(g_{a\gamma})^{2}$ (22)
$\displaystyle+\left(\frac{\partial\mu}{\partial\beta_{0}}\right)^{2}\Delta\beta_{0}^{2}+\left(\frac{\partial\mu}{\partial\gamma}\right)^{2}\Delta\gamma^{2}+\left(\frac{\partial\mu}{\partial\delta}\right)^{2}\Delta\delta^{2}\bigg{]}^{1/2},$
(23)

where $\delta$ is fitted as a nuisance parameter to account for the intrinsic
scattering in the quasar dataset. See Appendix A for the explicit likelihood
function. We define a distance between the best fit ($\mu_{\rm bf}$) and the
concordance cosmology ($\mu_{\rm c}$) as

$\displaystyle\chi^{2}_{\rm c}$ $\displaystyle=\sum_{i}\left(\frac{\mu_{\rm
bf,i}-\mu_{\rm c,i}}{\Delta_{i}}\right)^{2},$ (24)

where $i$ is the index of each QSO data point, and
$\Delta=\sqrt{\Delta\mu_{\rm bf}^{2}+\Delta\mu_{\rm
c}^{2}}\approx\Delta\mu_{\rm bf}$, since $\Delta\mu_{\rm bf}\gg\Delta\mu_{\rm
c}$. In practice, we neglect the contribution to $\Delta\mu_{\rm bf}$ from
$\Delta\log(m_{a})$ and $\Delta\log(g_{a\gamma})$ for simplicity. Since the
purpose of $\chi_{c}^{2}$ is to quantify the compatibility between the best
fit point and the concordance cosmology, neglecting both makes the estimate
conservative, i.e. it is easier to spot any potential incompatibility.

We show the theory distance as defined above in the last column of Tab. 2.
Larger value of $\chi^{2}_{c}$ show a larger incompatibility between the best
fit point and the concordance cosmology expansion histories. In $w$CDM and the
cosmographic model, the incompatibility is much larger than the axion model
even after adding SNIa and BAO. This is because the axion model mitigates some
of the incompatibility by correcting the quasar flux evolution with the axion-
induced attenuation $P_{\gamma\gamma,\rm UV,X}(z)$. After this correction the
inferred expansion history becomes similar to the concordance cosmology. By
contrast, in the other two models the quasar data bends the luminosity
distance $D_{L}(z)$ away from the concordance ${\Lambda\mathrm{CDM}}$
directly, leading to a larger incompatibility between the best fit points and
the concordance cosmology.

###  IV.4 Comment on the Axion Model

Figure 1: 2D posteriors of the axion mass ($m_{a}$) and coupling to photons
($g_{a\gamma}$) in the $\mathrm{axion}(\mathcal{Q}_{\beta_{0}}+\mathcal{G})$
model, with a magnetic domain size fixed at a comoving length of
$s=1~{}\mathrm{Mpc}$ (blue) and $s=100~{}\mathrm{Mpc}$ (orange). The gray
curves correspond to other astrophysical constraints, the super star clusters
Dessert:2020lil , M87 Marsh:2017yvc , NGC 1275Reynolds:2019uqt , and H1821+643
Reynes:2021bpe .

As demonstrated in the last section, the axion model can alleviate the tension
between the quasar and other cosmological datasets, and its best fit point is
more consistent with the concordance cosmology. This motivates us to examine
further the posterior of the axion model fit.

In Fig. 1, we compare the 2D posterior of the axion parameters with other
independent astrophysical constraints Reynes:2021bpe ; Reynolds:2019uqt ;
Marsh:2017yvc ; Dessert:2020lil . We observe that such large coupling is
disfavored by those constraints. See Fig. 4 for more detailed 2D posterior
information, and Tab. 4 for the full 1D posterior.

On the other hand, since these constraints rely on galactic magnetic field or
the magnetic field in the intracluster medium, the uncertainty in the domain
size of the IGM magnetic field does not affect these bounds. As we discussed
in Sec. IV.2, the allowed domain size of the IGM magnetic field can easily
vary from $1~{}\mathrm{Mpc}$ to $100~{}\mathrm{Mpc}$. Because of this, we also
fit the axion model with a comoving IGM magnetic field domain length set to
$s=100~{}\mathrm{Mpc}$. This results in a more efficient photon-to-axion
conversion, as we showed in Ref. Buen-Abad:2020zbd . Therefore, the model
requires a smaller $g_{a\gamma}$ to generate the redshift-dependent feature
preferred by the QSO dataset. The results from taking $s=100~{}\mathrm{Mpc}$
are still in $\sim 2\sigma$ tension with H1821+643 Reynes:2021bpe  and NGC
1275 Reynolds:2019uqt ,888The intracluster medium magnetic field modeling for
NGC 1275 bound is questioned in Libanov:2019fzq . although they are still
allowed by M87 Marsh:2017yvc  and Super Star Cluster Dessert:2020lil . Another
strong constraint on the axion-photon coupling comes from the CMB spectral
distortion Mirizzi:2009nq ; Tashiro:2013yea . While a large swath of axion
mass values is thereby ruled out, the region $10^{-14}\,\mathrm{eV}\lesssim
m_{a}\lesssim 10^{-12}\,\mathrm{eV}$ can be allowed if (i) multiple resonant
conversions, and (ii) tuning of the single resonant conversion probabilities
both take place Tashiro:2013yea .

As an extra test to understand the axion-induced photon attenuation, we plot
out the photon-to-axion conversion probability at the best fit point of
axion($\mathcal{Q}_{\beta_{0}}+\mathcal{G}$) in Fig. 2, for different
frequencies.

Figure 2: Photon-to-axion conversion probability, at redshifts corresponding
to all quasars, in the X-ray band ($\omega=2~{}\mathrm{keV}$, blue), UV band
($\omega=4.96~{}\mathrm{eV}$, orange), and optical band
($\omega=1\,\mathrm{eV}$, black), for the
axion($\mathcal{Q}_{\beta_{0}}+\mathcal{G}$) best fit point. We also recompute
the UV photon survival probability by increasing $g_{a\gamma}$ by 16 times and
$m_{a}$ by 4 times in the orange curve. This shows that $P_{\rm UV}(z)$ is
almost invariant along the curve $(g_{a\gamma}/g_{a\gamma,\rm
bf})=(m_{a}/m_{a,\rm bf})^{2}$, where $(\cdot)_{\rm bf}$ is the best fit point
of the parameter. The domain size is chosen to be $s=1~{}\mathrm{Mpc}$.

The conversion is very effective in the X-ray frequency range, which
transitions from zero to 1/3, the saturated value, mostly below $z\sim 0.3$.
The conversion in the UV, on the other hand, is almost negligible below $z\sim
1$, merely at the (sub)percent level, while it has a large redshift dependence
at $1\lesssim z\lesssim 2.5$. Above $z\sim 2.5$ the UV conversion also
saturates at 1/3. Put together, the change in $F_{\rm X}$ and $F_{\rm UV}$
leads to a modification of Eq. (3), i.e. the UV–X-ray flux relation, with Eq.
(7). Such a modification is favored by the QSO data due to its tendency to
prefer what amounts to an effective change in $\beta$ around redshift $z\sim
2$, as shown by the
${\Lambda\mathrm{CDM}}({\mathcal{Q}_{\beta_{\sharp}}}+\mathcal{G})$ and
${\Lambda\mathrm{CDM}}({\mathcal{Q}_{\beta_{\sharp}}}+\mathcal{B}+\mathcal{G})$
tests in Sec. III; we discuss this point in more detail in the next section.
We further test that the axion-to-photon conversion in the UV band is in the
non-linear minimal mixing regime (i.e. $k\approx\Delta_{a}\gg 1/s$ in Eq. (12)
with $s$ being the magnetic domain size). In this regime, the conversion
probability has a power-law dependence on both the axion mass and the
coupling, $P_{a\gamma}\propto g_{a\gamma}^{2}/m_{a}^{4}$. To verify this, we
shift the axion parameters in the direction of $g_{a\gamma}\propto m_{a}^{2}$,
away from the best fit point. We show that the axion-to-photon conversion is
almost invariant along this curve in the $m_{a}-g_{a\gamma}$ plane, as shown
by the green and orange curves in Fig. 2. Note that this behavior corresponds
to the degeneracy exhibited by the 2D posterior contours in the
$m_{a}-g_{a\gamma}$ parameter space of Fig. 1.

We can also understand why the axion island is allowed even after SNIa is
added. In Fig. 2 we show the impact of the photon-to-axion conversion in the
optical band (black points). Below $z\sim 2$, which is the maximum redshift to
which the Pantheon dataset extends, the modification of the photon flux in
optical band is minimal, at most at the $\sim 2\%$ level. Therefore, this
modification hardly generates any perceivable changes in the SNIa dataset even
though the coupling is relatively large.

Lastly, we want to comment on a variation of the axion model, which might
evade the astrophysical constraints listed in Fig. 1. It is also possible to
have axion-photon conversion in a cosmic dark magnetic field $B^{\prime}$, in
the scenario with a dark photon, $\gamma^{\prime}$, and an axion-photon-dark
photon coupling $g_{a\gamma^{\prime}}aF\tilde{F^{\prime}}$, where
$\tilde{F^{\prime}}$ is the dual field strength of the gauged dark
$U(1)^{\prime}$. In this case, the conversion probability in a single magnetic
domain is similar to Eq. (12) and (13) with $g_{a\gamma}B$ replaced by
$g_{a\gamma^{\prime}}B^{\prime}$. Assuming that the coherent length of
$B^{\prime}$ is similar to that of the IGM magnetic field ($\sim{\cal
O}(1-100)$ Mpc), from Fig. 1, we need $g_{a\gamma^{\prime}}B^{\prime}\gtrsim
10^{-12}$ GeV${}^{-1}\times 1$ nG to restore the cosmic concordance. The
astrophysical constraints on $g_{a\gamma^{\prime}}$ and $B^{\prime}$ are
weaker compared to those on $g_{a\gamma}$ and the IGM $B$ field. More
specifically, $B^{\prime}$ could be as large as micro-Gauss, with the
constraints mostly from $N_{\rm eff}$ while $g_{a\gamma^{\prime}}\lesssim
5\times 10^{-10}$ GeV-1 due to star cooling Choi:2018mvk ; Choi:2019jwx . This
implies a potential new physics explanation for the cosmic concordance of the
quasar data, without conflicts with the other constraints. We leave this for
further investigation in future work.

###  IV.5 Effective Evolution of $\beta$

As we discuss in Sec. II.2, new physics modifications can be parametrized as
an effective redshift evolution in $\beta$. We summarize the results by
showing this effective redshift evolution, $\beta_{\rm eff}(z)$, in Fig. 3.

Figure 3: Effective redshift dependence of the nuisance parameter
$\beta_{\mathrm{eff}}(z)$. We show the amount of change in
$\beta_{\mathrm{eff}}(z)$ by measuring it with respect to its value at large
$z$ (in this plot, $z=8$). In the upper plot we show this difference for the
best fit points of the runs for
${\Lambda\mathrm{CDM}}({\mathcal{Q}_{\beta_{\sharp}}}+\mathcal{G})$,
${\Lambda\mathrm{CDM}}({\mathcal{Q}_{\beta_{\rm t}}}+\mathcal{G})$,
${a_{2}a_{3}a_{4}}(\mathcal{Q}_{\beta_{0}}+\mathcal{G})$,
axion$(\mathcal{Q}_{\beta_{0}}+\mathcal{G})$, and
${w\mathrm{CDM}}(\mathcal{Q}_{\beta_{0}}+\mathcal{G})$. In the lower panel we
do the same but including SNIa and BAO in the fits ($+\mathcal{B}$ dataset).
We find a good agreement between our results that make use of
${\mathcal{Q}_{\beta_{\rm t}}}$, and the effective evolution $\beta_{\rm eff}$
of Eq. (4) with parameter values taken from Ref. Dainotti:2022rfz .

In Sec. III, when fitted to the QSO data alone, we see roughly two classes of
solutions: a shallow but abrupt change in $\beta$ around $z\sim 1.5$ and a
deep and slow change of $\beta$ spanning the whole redshift range. The axion
model leads to a relatively quick change in redshift through
$P_{\gamma\gamma,\rm UV,X}(z)$. This serves as a good physical realization of
the sharp change given by the effective parametrization in
${\mathcal{Q}_{\beta_{\sharp}}}$. As we see in Sec. IV.4, the effective change
in $\beta$ around $z=1.5$ is due to the UV photon conversion to axions.
Interestingly, the maximum amount of change allowed in the UV is given by

$\displaystyle\Delta\beta=\gamma\log(2/3)\approx-0.12.$ (25)

This size of the modification in the UV-X-ray relation happens to be what the
data prefers, as evidenced by the run with the effective parametrization with
a sharp transition in $\beta$ ($\beta_{\sharp}$ in the plot). We would like to
highlight the fact that when the QSO dataset is allowed to have a sharp
transition in the UV-X-ray relation (i.e. ${\mathcal{Q}_{\beta_{\sharp}}}$),
the fits prefer this transition to take place at around the same redshift as
our axion model does. The striking agreement between both can be seen by
comparing ${\Lambda\mathrm{CDM}}({\mathcal{Q}_{\beta_{\sharp}}}+\mathcal{G})$
with $\mathrm{axion}(\mathcal{Q}_{\beta_{0}}+\mathcal{G})$, or by comparing
${\Lambda\mathrm{CDM}}({\mathcal{Q}_{\beta_{\sharp}}}+\mathcal{B}+\mathcal{G})$
and $\mathrm{axion}(\mathcal{Q}_{\beta_{0}}+\mathcal{B}+\mathcal{G})$. Adding
to its appeal, the axion model outperforms non-standard cosmological models,
even those as flexible as the cosmographic model, as we show in Sec. IV.3.

By contrast, the $w$CDM and cosmographic models have large changes when
$\mathcal{B}$ set is added. One should keep in mind that in the
${\mathrm{B}\Lambda\mathrm{CDM}}$ models, the effective evolution in $\beta$
is achieved at the price of directly modifying the Universe expansion history.
Therefore, SNIa+BAO will constrain the low-$z$ part ($z\lesssim 1.5$) of
$\beta_{\rm eff}(z)$ to be as flat as possible. This has consequences in both
the $w$CDM and the cosmographic model. In $w$CDM, having a flat low-$z$ part
of $\beta_{\rm eff}(z)$ restricts the overall amount of change one can achieve
from $z=0$ to $z\sim 8$. This is reflected by the change in the curve for
${w\mathrm{CDM}}(\mathcal{Q}_{\beta_{0}}+\mathcal{G})$ and
${w\mathrm{CDM}}(\mathcal{Q}_{\beta_{0}}+\mathcal{G}+\mathcal{B})$: adding
$\mathcal{B}$ significantly reduces the amount of change in
$\beta_{\mathrm{eff}}(z)$. On the other hand, the cosmographic model has
enough degrees of freedom to guarantee a relatively flat low-$z$ $\beta_{\rm
eff}(z)$. Therefore, while restricting the flatness at low-$z$ changes the
shape of the effective $\beta$ evolution, the large modification in $\beta$
survives. However, this large $\beta$ change itself comes at the price of a
large deviation from the concordance ${\Lambda\mathrm{CDM}}$ at $1\lesssim
z\lesssim 8$. It also remains unclear what physical models can achieve such a
modification in $D_{L}(z)$ as indicated by the cosmographic model.

##  V Conclusion

Although quasar data can be used to extract information about the late-time
cosmological history of the Universe, it stands in need of standardization. In
light of the recent efforts along these lines, we perform a critical
examination on the imprints that new physics may leave on the quasar data, in
the form of unaccounted-for redshift evolution of their fluxes. This is a
timely endeavor, given the ongoing Hubble tension, as well as the discrepancy
between the cosmological expansion inferred from quasars and other
observations.

We use a flexible fitting template to identify what is required to modify the
quasar flux in such a way that the concordance cosmological history given by
in Eq. (1) can be restored. We test ${w\mathrm{CDM}}$ and the cosmographic
model, two cosmological models beyond ${\Lambda\mathrm{CDM}}$, as well as the
axion model, a mechanism beyond the SM that can present extra photon
attenuation. We find that the axion model has an advantage over the other two,
restoring the consistency between the best-fit parameters emerging from the
quasar dataset, and the concordance cosmology. It also outperforms the
${\mathrm{B}\Lambda\mathrm{CDM}}$ models in improving the goodness of fit to
the quasar dataset alone.

Lastly we note that although the best-fit axion model parameters are in mild
tension with several astrophysical constraints, this tension can be relieved
by allowing for larger magnetic domain sizes in the IGM, or for the axion-
photon conversion to take place in a dark magnetic field. These variants of
our baseline axion model may evade other constraints while at the same time
providing the redshift evolution of the quasar fluxes required by observation,
maintaining the validity of the RL relation, and restoring the consistency of
the quasar data with the concordance cosmological model.

Acknowledgment We thank Daniele Alves, Michael Graesser, Hui Li, Nicole Lloyd-
Ronning, and Elisabeta Lusso, for useful comments. We thank Tomer Volansky for
the access to Tel Aviv Univeristy High-Performance Computing, where part of
the work was performed. This work was supported by the U.S. Department of
Energy through the Los Alamos National Laboratory. Los Alamos National
Laboratory is operated by Triad National Security, LLC, for the National
Nuclear Security Administration of U.S. Department of Energy (Contract No.
89233218CNA000001). Research presented in this article was supported by the
Laboratory Directed Research and Development program of Los Alamos National
Laboratory under project number 3W200A-XXG900, task SUN00000. MBA is supported
in part by the National Science Foundation under Grant Number PHY-2210361, and
the Maryland Center for Fundamental Physics. JF is supported by the NASA grant
80NSSC22K081 and the DOE grant DE-SC-0010010.

Figure 4: 2D posteriors of the parameters of the cosmographic (top panel),
${w\mathrm{CDM}}$ (middle panel), and axion (bottom panel) models. The
magnetic domain has a comoving size $s=1~{}\mathrm{Mpc}$.

##  Appendix A Quasar Dataset

We use the quasar data sample from Ref. Lusso:2020pdb . It contains seven
different groups: the XMM-Newton $z\sim 3$ sample, the new XMM-Newton $z\sim
4$ quasars, the High-$z$ sample, XXL, SDSS - 4XMM, SDSS - Chandra, and local
AGN.

The UV–X-ray relation is expected to hold up to a certain amount of scattering
that cannot be eliminated by the improvement of the measurement. Since this
intrinsic scattering of the quasars is unknown, we add a third nuisance
parameter, $\delta$, describing the standard deviation of the intrinsic
scattering, along with $\gamma$ and $\beta$. As a result, the log-likelihood
function is given by

$\displaystyle-2\,\ln(\mathrm{likelihood})$
$\displaystyle=\sum_{i}\left[\left(x^{\rm th}-x^{\rm
obs}\right)^{2}/\sigma_{i}^{2}+\ln(2\pi\sigma^{2})\right],$ (26)

where $\delta$ is combined into the error of each quasar flux measurement as
$\sigma_{i}^{2}=\sigma_{i,\rm\log F_{\rm UV,X}}^{2}+\delta^{2}$. By adopting
Eq. (26), $\delta$ can be treated as a nuisance parameter that gets fitted
together in the MCMC. Therefore, specific to the QSO dataset, we have the
following nuisance parameters: $\beta,\gamma,\delta$, with $\beta$ itself
parametrize in the various ways described in Sec. III.

##  Appendix B Methodology

We use emcee Foreman_Mackey_2013  to fit the three models (${w\mathrm{CDM}}$,
cosmographic, and $\mathrm{axion}$) to two data combinations
($\mathcal{Q}_{\beta_{0}}+\mathcal{G}$ and
$\mathcal{Q}_{\beta_{0}}+\mathcal{B}+\mathcal{G}$). Our code is publicly
available at github.com/ChenSun-Phys/high_z_candles. We use 100 walkers and a
chain length of 40000. All runs pass our MCMC convergence condition of the
chain length to be 50 times longer than the auto-correlation length, except
for ${w\mathrm{CDM}}(\mathcal{Q}_{\beta_{0}}+\mathcal{B}+\mathcal{G})$, where
we loosen the condition to 12 times the auto-correlation length due to its
strong bimodal posterior.

The results are analyzed with corner.py corner . Since different parameters
have different auto-correlation times, we choose the burn-in to be equal to
twice the largest one, and choose the thinning-length to be equal to half of
the smallest one.

We use the following flat priors for the theory parameters.

$\displaystyle\Omega_{\Lambda}\in(0,1)~{}~{}~{}h\in(0.6,0.8)$ (27)
$\displaystyle w_{0}\in(0,1)~{}~{}~{}w_{a}\in(0,1)$ (28) $\displaystyle
a_{2}\in(1,5)~{}~{}a_{3}\in(-6,6)~{}~{}a_{4}\in(-10,6)$ (29)
$\displaystyle\log\left[\frac{g_{a\gamma}}{\mathrm{GeV}^{-1}}\right]\in(-18,-8)~{}~{}\log\left[\frac{m_{a}}{\mathrm{eV}}\right]\in(-17,-11).$

We set the following priors for the nuisance parameters when the corresponding
datasets are used.

$\displaystyle\mathrm{SNIa}:M_{0}\in(-21,-18)$ (30)
$\displaystyle\mathrm{BAO}:r_{s}/\mathrm{Mpc}\in(120,160)$ (31)
$\displaystyle\mathrm{QSO}:\gamma\in(0.1,1)~{}~{}\delta\in(0.05,0.6)$ (32)
$\displaystyle~{}\mathcal{Q}_{\beta_{0}}:\beta_{0}\in(0,10)$ (33)
$\displaystyle~{}{\mathcal{Q}_{\beta_{\sharp}}}:\beta_{0}\in(0,10)~{}~{}\beta_{1}\in(0,10)~{}~{}z_{0}\in(0,9)$
(34) $\displaystyle~{}{\mathcal{Q}_{\beta_{\rm
t}}}:\beta_{0}\in(0,10)~{}~{}\beta_{1}\in(0,10)~{}~{}z_{0}\in(0,9)~{}~{}\delta
z\in(0.01,10)~{}.$

| $\Omega_{\Lambda}$ | $h$ | $\beta_{0}$ | $\beta_{1}$ | $z_{0}$ | $\delta z$ | $\gamma$ | $\delta$ | $M_{0}$ | $r_{s}\mathrm{[Mpc]}$
---|---|---|---|---|---|---|---|---|---|---
$\mathcal{B}+\mathcal{G}$ | $0.68^{+0.02}_{-0.02}$ | $0.676^{+0.004}_{-0.004}$ |  |  |  |  |  |  | $-19.41^{+0.01}_{-0.01}$ | $147.12^{+0.26}_{-0.26}$
$\mathcal{Q}_{\beta_{0}}+\mathcal{G}$ | $0.05^{+0.07}_{-0.03}$ | $0.674^{+0.005}_{-0.005}$ | $7.04^{+0.24}_{-0.24}$ |  |  |  | $0.639^{+0.008}_{-0.008}$ | $0.228^{+0.004}_{-0.003}$ |  |
${\mathcal{Q}_{\beta_{\sharp}}}+\mathcal{G}$ | $0.10^{+0.12}_{-0.07}$ | $0.674^{+0.005}_{-0.005}$ | $8.20^{+0.28}_{-0.28}$ | $8.30^{+0.29}_{-0.29}$ | $1.65^{+0.01}_{-0.01}$ |  | $0.599^{+0.009}_{-0.009}$ | $0.225^{+0.003}_{-0.003}$ |  |
${\mathcal{Q}_{\beta_{\rm t}}}+\mathcal{G}$ | $0.46^{+0.26}_{-0.27}$ | $0.674^{+0.005}_{-0.005}$ | $8.38^{+0.37}_{-0.36}$ | $9.37^{+0.35}_{-0.36}$ | $4.91^{+3.27}_{-2.05}$ | $1.24^{+1.13}_{-0.86}$ | $0.580^{+0.011}_{-0.010}$ | $0.224^{+0.003}_{-0.003}$ |  |
$\mathcal{Q}_{\beta_{0}}+\mathcal{B}+\mathcal{G}$ | $0.66^{+0.02}_{-0.02}$ | $0.674^{+0.004}_{-0.004}$ | $6.40^{+0.23}_{-0.23}$ |  |  |  | $0.662^{+0.008}_{-0.008}$ | $0.230^{+0.004}_{-0.004}$ | $-19.41^{+0.01}_{-0.01}$ | $147.09^{+0.26}_{-0.26}$
${\mathcal{Q}_{\beta_{\sharp}}}+\mathcal{B}+\mathcal{G}$ | $0.67^{+0.02}_{-0.02}$ | $0.675^{+0.004}_{-0.004}$ | $7.79^{+0.27}_{-0.27}$ | $7.90^{+0.28}_{-0.28}$ | $1.65^{+0.01}_{-0.01}$ |  | $0.616^{+0.009}_{-0.009}$ | $0.226^{+0.003}_{-0.003}$ | $-19.41^{+0.01}_{-0.01}$ | $147.10^{+0.26}_{-0.26}$
${\mathcal{Q}_{\beta_{\rm t}}}+\mathcal{B}+\mathcal{G}$ | $0.68^{+0.02}_{-0.02}$ | $0.676^{+0.004}_{-0.004}$ | $8.40^{+0.39}_{-0.40}$ | $9.40^{+0.38}_{-0.37}$ | $4.30^{+2.75}_{-1.52}$ | $0.92^{+0.76}_{-0.64}$ | $0.579^{+0.011}_{-0.011}$ | $0.224^{+0.003}_{-0.003}$ | $-19.41^{+0.01}_{-0.01}$ | $147.12^{+0.26}_{-0.26}$
Table 3:  The mean and $1\sigma$ values of the ${\Lambda\mathrm{CDM}}$
parameters, fitted to various data combinations.

cosmographic | $h$ | $a_{2}$ | $a_{3}$ | $a_{4}$ | $\beta_{0}$ | $\gamma$ | $\delta$ | $M_{0}$ | $r_{s}\mathrm{[Mpc]}$
---|---|---|---|---|---|---|---|---|---
$\mathcal{Q}_{\beta_{0}}+\mathcal{G}$ | $0.674^{+0.005}_{-0.005}$ | $3.66^{+0.71}_{-0.86}$ | $-4.35^{+1.97}_{-1.19}$ | $1.61^{+1.25}_{-1.69}$ | $8.81^{+0.32}_{-0.32}$ | $0.579^{+0.011}_{-0.011}$ | $0.224^{+0.004}_{-0.003}$ |  |
$\mathcal{Q}_{\beta_{0}}+\mathcal{B}+\mathcal{G}$ | $0.677^{+0.005}_{-0.005}$ | $2.95^{+0.15}_{-0.12}$ | $5.19^{+0.55}_{-0.78}$ | $-8.01^{+1.05}_{-0.82}$ | $7.67^{+0.27}_{-0.27}$ | $0.619^{+0.009}_{-0.009}$ | $0.226^{+0.003}_{-0.003}$ | $-19.42^{+0.01}_{-0.01}$ | $147.13^{+0.26}_{-0.26}$

${w\mathrm{CDM}}$ | $\Omega_{\Lambda}$ | $h$ | $w_{0}$ | $w_{a}$ | $\beta_{0}$ | $\gamma$ | $\delta$ | $M_{0}$ | $r_{s}\mathrm{[Mpc]}$
---|---|---|---|---|---|---|---|---|---
$\mathcal{Q}_{\beta_{0}}+\mathcal{G}$ | $0.28^{+0.28}_{-0.13}$ | $0.674^{+0.005}_{-0.005}$ | $0.77^{+0.17}_{-0.29}$ | $0.83^{+0.12}_{-0.24}$ | $8.61^{+0.30}_{-0.31}$ | $0.58^{+0.01}_{-0.01}$ | $0.225^{+0.003}_{-0.003}$ |  |
$\mathcal{Q}_{\beta_{0}}+\mathcal{B}+\mathcal{G}$ | $0.96^{+0.02}_{-0.08}$ | $0.675^{+0.005}_{-0.004}$ | $-0.86^{+0.03}_{-0.06}$ | $0.96^{+0.03}_{-0.10}$ | $6.68^{+0.25}_{-0.25}$ | $0.653^{+0.008}_{-0.008}$ | $0.229^{+0.004}_{-0.003}$ | $-19.43^{+0.015}_{-0.014}$ | $147.10^{+0.26}_{-0.25}$

axion | $\Omega_{\Lambda}$ | $h$ | $\log(m_{a})$ | $\log(g_{a\gamma})$ | $\beta_{0}$ | $\gamma$ | $\delta$ | $M_{0}$ | $r_{s}\mathrm{[Mpc]}$
---|---|---|---|---|---|---|---|---|---
$\mathcal{Q}_{\beta_{0}}+\mathcal{G}$ | $0.09^{+0.12}_{-0.06}$ | $0.674^{+0.005}_{-0.005}$ | $-12.85^{+0.58}_{-0.46}$ | $-9.77^{+1.21}_{-0.80}$ | $8.18^{+0.25}_{-0.28}$ | $0.606^{+0.009}_{-0.008}$ | $0.225^{+0.004}_{-0.003}$ |  |
$\mathcal{Q}_{\beta_{0}}+\mathcal{B}+\mathcal{G}$ | $0.67^{+0.02}_{-0.02}$ | $0.675^{+0.005}_{-0.004}$ | $-12.71^{+0.47}_{-0.37}$ | $-9.39^{+0.95}_{-0.77}$ | $7.68^{+0.23}_{-0.22}$ | $0.625^{+0.007}_{-0.007}$ | $0.226^{+0.004}_{-0.003}$ | $-19.41^{+0.01}_{-0.01}$ | $147.11^{+0.26}_{-0.26}$

Table 4:  The mean and $1\sigma$ values of the parameters of the cosmographic
(top), ${w\mathrm{CDM}}$ (middle), and axion (bottom) models, fitted to
various data combinations.

##  Appendix C Posterior of the Fits

We show the 1D, $1\sigma$ posterior range of the ${\Lambda\mathrm{CDM}}$
parameters in Tab. 3. The corresponding results for the
${\mathrm{B}\Lambda\mathrm{CDM}}$ and BSM fits are shown in Tab. 4. The 2D
posterior of the ${\mathrm{B}\Lambda\mathrm{CDM}}$ and BSM parameters are
shown in Fig. 4.

## References

  * (1) D. Foreman-Mackey, D. W. Hogg, D. Lang, and J. Goodman, “emcee: The mcmc hammer,” Publications of the Astronomical Society of the Pacific 125 no. 925, (Mar, 2013) 306–312. http://dx.doi.org/10.1086/670067.
  * (2) D. Foreman-Mackey, “corner.py: Scatterplot matrices in python,” The Journal of Open Source Software 1 no. 2, (Jun, 2016) 24. https://doi.org/10.21105/joss.00024.
  * (3) G. Risaliti and E. Lusso, “Cosmological constraints from the Hubble diagram of quasars at high redshifts,” Nature Astron. 3 no. 3, (2019) 272–277, arXiv:1811.02590 [astro-ph.CO].
  * (4) E. Lusso, E. Piedipalumbo, G. Risaliti, et al., “Tension with the flat {\}Lambda{}CDM model from a high redshift Hubble Diagram of supernovae, quasars and gamma-ray bursts,” Astronomy & Astrophysics 628 (Aug., 2019) L4, arXiv:1907.07692 [astro-ph]. http://arxiv.org/abs/1907.07692.
  * (5) T. Yang, A. Banerjee, and E. Ó. Colgáin, “Cosmography and flat $\Lambda$CDM tensions at high redshift,” Physical Review D 102 no. 12, (Dec., 2020) 123532, arXiv:1911.01681 [astro-ph, physics:gr-qc]. http://arxiv.org/abs/1911.01681.
  * (6) H. Velten and S. Gomes, “Is the Hubble diagram of quasars in tension with concordance cosmology?,” Physical Review D 101 no. 4, (Feb., 2020) 043502, arXiv:1911.11848 [astro-ph]. http://arxiv.org/abs/1911.11848.
  * (7) M. G. Dainotti, G. Bargiacchi, M. Bogdan, et al., “Reduced uncertainties up to 43% on the Hubble constant and the matter density with the SNe Ia with a new statistical analysis,” arXiv:2303.06974 [astro-ph.CO].
  * (8) M. G. Dainotti, G. Bardiacchi, A. L. Lenart, et al., “Quasar Standardization: Overcoming Selection Biases and Redshift Evolution,” Astrophys. J. 931 no. 2, (2022) 106, arXiv:2203.12914 [astro-ph.HE].
  * (9) A. L. Lenart, G. Bargiacchi, M. G. Dainotti, et al., “A Bias-free Cosmological Analysis with Quasars Alleviating H 0 Tension,” Astrophys. J. Suppl. 264 no. 2, (2023) 46, arXiv:2211.10785 [astro-ph.CO].
  * (10) G. Bargiacchi, M. G. Dainotti, S. Nagataki, and S. Capozziello, “Gamma-Ray Bursts, Quasars, Baryonic Acoustic Oscillations, and Supernovae Ia: new statistical insights and cosmological constraints,” arXiv:2303.07076 [astro-ph.CO].
  * (11) A. G. Riess et al., “A 2.4% Determination of the Local Value of the Hubble Constant,” Astrophys. J. 826 no. 1, (2016) 56, arXiv:1604.01424 [astro-ph.CO].
  * (12) A. G. Riess, S. Casertano, W. Yuan, et al., “Large Magellanic Cloud Cepheid Standards Provide a 1% Foundation for the Determination of the Hubble Constant and Stronger Evidence for Physics beyond $\Lambda$CDM,” Astrophys. J. 876 no. 1, (2019) 85, arXiv:1903.07603 [astro-ph.CO].
  * (13) K. Aylor, M. Joy, L. Knox, et al., “Sounds Discordant: Classical Distance Ladder $\&$ $\Lambda$CDM -based Determinations of the Cosmological Sound Horizon,” Astrophys. J. 874 no. 1, (2019) 4, arXiv:1811.00537 [astro-ph.CO].
  * (14) L. Knox and M. Millea, “Hubble constant hunter’s guide,” Phys. Rev. D 101 no. 4, (2020) 043533, arXiv:1908.03663 [astro-ph.CO].
  * (15) E. Abdalla, G. F. Abellán, A. Aboubrahim, et al., “Cosmology Intertwined: A Review of the Particle Physics, Astrophysics, and Cosmology Associated with the Cosmological Tensions and Anomalies,” Journal of High Energy Astrophysics 34 (June, 2022) 49–211, arxiv:2203.06142 [astro-ph, physics:hep-ph]. http://arxiv.org/abs/2203.06142.
  * (16) Planck Collaboration, N. Aghanim et al., “Planck 2018 results. VI. Cosmological parameters,” Astron. Astrophys. 641 (2020) A6, arXiv:1807.06209 [astro-ph.CO].
  * (17) D. M. Scolnic, D. O. Jones, A. Rest, et al., “The Complete Light-curve Sample of Spectroscopically Confirmed Type Ia Supernovae from Pan-STARRS1 and Cosmological Constraints from The Combined Pantheon Sample,” The Astrophysical Journal 859 no. 2, (May, 2018) 101, arxiv:1710.00845. http://arxiv.org/abs/1710.00845.
  * (18) A. J. Ross, F. Beutler, C.-H. Chuang, et al., “The clustering of galaxies in the completed SDSS-III Baryon Oscillation Spectroscopic Survey: Observational systematics and baryon acoustic oscillations in the correlation function,” Monthly Notices of the Royal Astronomical Society 464 no. 1, (Jan., 2017) 1168–1191, arxiv:1607.03145. http://arxiv.org/abs/1607.03145.
  * (19) F. Beutler, H.-J. Seo, A. J. Ross, et al., “The clustering of galaxies in the completed SDSS-III Baryon Oscillation Spectroscopic Survey: Baryon Acoustic Oscillations in Fourier-space,” Monthly Notices of the Royal Astronomical Society 464 no. 3, (Jan., 2017) 3409–3430, arxiv:1607.03149. http://arxiv.org/abs/1607.03149.
  * (20) E. Lusso and G. Risaliti, “Quasars as standard candles I: The physical relation between disc and coronal emission,” Astronomy & Astrophysics 602 (June, 2017) A79, arXiv:1703.05299 [astro-ph]. http://arxiv.org/abs/1703.05299.
  * (21) C. Csaki, N. Kaloper, and J. Terning, “Dimming supernovae without cosmic acceleration,” Phys. Rev. Lett. 88 (2002) 161302, arXiv:hep-ph/0111311 [hep-ph].
  * (22) C. Csaki, N. Kaloper, and J. Terning, “Effects of the intergalactic plasma on supernova dimming via photon axion oscillations,” Phys. Lett. B535 (2002) 33–36, arXiv:hep-ph/0112212 [hep-ph].
  * (23) M. A. Buen-Abad, J. Fan, and C. Sun, “Constraints on axions from cosmic distance measurements,” JHEP 02 (2022) 103, arXiv:2011.05993 [hep-ph].
  * (24) F. Haardt and L. Maraschi, “A two-phase model for the X-ray emission from Seyfert galaxies,” The Astrophysical Journal 380 (Oct., 1991) L51. http://adsabs.harvard.edu/doi/10.1086/186171.
  * (25) F. Haardt, L. Maraschi, and G. Ghisellini, “A Model for the X-Ray and UV Emission from Seyfert Galaxies and Galactic Black Holes.” http://arxiv.org/abs/astro-ph/9405059, May, 1994.
  * (26) S. Martocchia, E. Piconcelli, L. Zappacosta, et al., “The WISSH quasars project: III. X-ray properties of hyper-luminous quasars,” Astronomy & Astrophysics 608 (Dec., 2017) A51. http://www.aanda.org/10.1051/0004-6361/201731314.
  * (27) F. Salvestrini, G. Risaliti, S. Bisogni, et al., “Quasars as standard candles II: The non linear relation between UV and X-ray emission at high redshifts,” Astronomy & Astrophysics 631 (Nov., 2019) A120, arXiv:1909.12309 [astro-ph]. http://arxiv.org/abs/1909.12309.
  * (28) Y. Avni and H. Tananbaum, “On the cosmological evolution of the X-ray emission from quasars,” The Astrophysical Journal 262 (Nov., 1982) L17–L21. https://ui.adsabs.harvard.edu/abs/1982ApJ...262L..17A.
  * (29) C. Vignali, W. N. Brandt, D. P. Schneider, et al., “Chandra and XMM-Newton Observations of the First Quasars: X-Rays from the Age of Cosmic Enlightenment,” The Astronomical Journal 125 no. 6, (June, 2003) 2876. https://iopscience.iop.org/article/10.1086/374998/meta.
  * (30) E. Lusso, A. Comastri, C. Vignali, et al., “The X-ray to optical-UV luminosity ratio of X-ray selected type 1 AGN in XMM-COSMOS,” Astronomy and Astrophysics 512 (Mar., 2010) A34. https://ui.adsabs.harvard.edu/abs/2010A&A...512A..34L.
  * (31) A. T. Steffen, I. Strateva, W. N. Brandt, et al., “The X-Ray-to-Optical Properties of Optically Selected Active Galaxies over Wide Luminosity and Redshift Ranges,” The Astronomical Journal 131 no. 6, (June, 2006) 2826. https://dx.doi.org/10.1086/503627.
  * (32) E. Lusso and G. Risaliti, “THE TIGHT RELATION BETWEEN X-RAY AND ULTRAVIOLET LUMINOSITY OF QUASARS,” The Astrophysical Journal 819 no. 2, (Mar., 2016) 154. https://dx.doi.org/10.3847/0004-637X/819/2/154.
  * (33) B. Efron and V. Petrosian, “A Simple Test of Independence for Truncated Data with Applications to Redshift Surveys,” The Astrophysical Journal 399 (Nov., 1992) 345. https://ui.adsabs.harvard.edu/abs/1992ApJ...399..345E.
  * (34) F. Beutler, C. Blake, M. Colless, et al., “The 6dF Galaxy Survey: Baryon Acoustic Oscillations and the Local Hubble Constant,” Monthly Notices of the Royal Astronomical Society 416 no. 4, (Oct., 2011) 3017–3032, arxiv:1106.3366. http://arxiv.org/abs/1106.3366.
  * (35) A. J. Ross, L. Samushia, C. Howlett, et al., “The clustering of the SDSS DR7 main Galaxy sample –I. A 4 per cent distance measure at $z~{}=~{}0.15$,” Mon. Not. Roy. Astron. Soc. 449 no. 1, (2015) 835–847, arXiv:1409.3242 [astro-ph.CO].
  * (36) BOSS Collaboration, S. Alam et al., “The clustering of galaxies in the completed SDSS-III Baryon Oscillation Spectroscopic Survey: cosmological analysis of the DR12 galaxy sample,” Mon. Not. Roy. Astron. Soc. 470 no. 3, (2017) 2617–2652, arXiv:1607.03155 [astro-ph.CO].
  * (37) E. Lusso, G. Risaliti, E. Nardini, et al., “Quasars as standard candles III. Validation of a new sample for cosmological studies,” Astronomy & Astrophysics 642 (Oct., 2020) A150, arXiv:2008.08586 [astro-ph]. http://arxiv.org/abs/2008.08586.
  * (38) H. Georgi, P. H. Ginsparg, and S. L. Glashow, “Photon Oscillations and the Cosmic Background Radiation,” Nature 306 (1983) 765–766.
  * (39) P. Sikivie, “Experimental Tests of the Invisible Axion,” Phys. Rev. Lett. 51 (1983) 1415–1417. [Erratum: Phys.Rev.Lett. 52, 695 (1984)].
  * (40) G. Raffelt and L. Stodolsky, “Mixing of the Photon with Low Mass Particles,” Phys. Rev. D 37 (1988) 1237.
  * (41) F. Nicastro et al., “Observations of the MIssing Baryons in the warm-hot intergalactic medium,” Nature 558 (2018) 406, arXiv:1806.08395 [astro-ph.GA].
  * (42) D. Martizzi et al., “Baryons in the Cosmic Web of IllustrisTNG – I: gas in knots, filaments, sheets, and voids,” Mon. Not. Roy. Astron. Soc. 486 no. 3, (2019) 3766–3787, arXiv:1810.01883 [astro-ph.CO].
  * (43) Y. Grossman, S. Roy, and J. Zupan, “Effects of initial axion production and photon axion oscillation on type Ia supernova dimming,” Phys. Lett. B 543 (2002) 23–28, arXiv:hep-ph/0204216.
  * (44) A. Avgoustidis, C. Burrage, J. Redondo, et al., “Constraints on cosmic opacity and beyond the standard model physics from cosmological distance measurements,” JCAP 10 (2010) 024, arXiv:1004.2053 [astro-ph.CO].
  * (45) R. Durrer and A. Neronov, “Cosmological Magnetic Fields: Their Generation, Evolution and Observation,” Astron. Astrophys. Rev. 21 (2013) 62, arXiv:1303.7121 [astro-ph.CO].
  * (46) C. Dessert, J. W. Foster, and B. R. Safdi, “X-ray Searches for Axions from Super Star Clusters,” arXiv:2008.03305 [hep-ph].
  * (47) M. D. Marsh, H. R. Russell, A. C. Fabian, et al., “A New Bound on Axion-Like Particles,” JCAP 12 (2017) 036, arXiv:1703.07354 [hep-ph].
  * (48) C. S. Reynolds, M. D. Marsh, H. R. Russell, et al., “Astrophysical limits on very light axion-like particles from Chandra grating spectroscopy of NGC 1275,” arXiv:1907.05475 [hep-ph].
  * (49) J. S. Reynés, J. H. Matthews, C. S. Reynolds, et al., “New constraints on light Axion-Like Particles using Chandra Transmission Grating Spectroscopy of the powerful cluster-hosted quasar H1821+643,” Monthly Notices of the Royal Astronomical Society 510 no. 1, (Dec., 2021) 1264–1277, arxiv:2109.03261. http://arxiv.org/abs/2109.03261.
  * (50) M. Libanov and S. Troitsky, “On the impact of magnetic-field models in galaxy clusters on constraints on axion-like particles from the lack of irregularities in high-energy spectra of astrophysical sources,” Phys. Lett. B 802 (2020) 135252, arXiv:1908.03084 [astro-ph.HE].
  * (51) A. Mirizzi, J. Redondo, and G. Sigl, “Constraining resonant photon-axion conversions in the Early Universe,” JCAP 0908 (2009) 001, arXiv:0905.4865 [hep-ph].
  * (52) H. Tashiro, J. Silk, and D. J. E. Marsh, “Constraints on primordial magnetic fields from CMB distortions in the axiverse,” Physical Review D 88 no. 12, (Dec., 2013) 125024, arxiv:1308.0314. http://arxiv.org/abs/1308.0314.
  * (53) K. Choi, S. Lee, H. Seong, and S. Yun, “Gamma-ray spectral modulations induced by photon-ALP-dark photon oscillations,” Phys. Rev. D 101 no. 4, (2020) 043007, arXiv:1806.09508 [hep-ph].
  * (54) K. Choi, H. Seong, and S. Yun, “Axion-photon-dark photon oscillation and its implication for 21 cm observation,” arXiv:1911.00532 [astro-ph, physics:hep-ph] (Nov., 2019) , arxiv:1911.00532 [astro-ph, physics:hep-ph]. http://arxiv.org/abs/1911.00532.