Also at ]

# Equality scale-based and sound horizon-based analysis of the Hubble tension

Zhihuan Zhou<EMAIL_ADDRESS>[  Yuhao Mu  Gang Liu  Lixin Xu
<EMAIL_ADDRESS>Institute of Theoretical Physics
School of Physics
Dalian University of Technology
Dalian 116024, People’s Republic of China  Jianbo Lu  Department of Physics
Liaoning Normal University
Dalian 116029, People’s Republic of China

(August 28, 2024)

###### Abstract

The Hubble horizon at matter-radiation equality ($k^{-1}_{\rm{eq}}$) and the
sound horizon at the last scattering surface ($r_{s}(z_{*})$) provides
interesting consistency check for the $\Lambda$CDM model and its extensions.
It is well known that the reduction of $r_{s}$ can be compensated by the
increase of $H_{0}$, while the same is true for the standard rulers
$k_{\rm{eq}}$. Adding extra radiational component to the early universe can
reduce $k_{\rm{eq}}$. The addition of early dark energy (EDE), however, tends
to increase $k_{\rm{eq}}$. We perform $k_{\rm{eq}}$\- and $r_{s}$-based
analyses in both the EDE model and the Wess-Zumino Dark Radiation (WZDR)
model. In the latter case we find $\Delta H_{0}=0.4$ between the $r_{s}$\- and
$k_{\rm{eq}}$-based datasets, while in the former case we find $\Delta
H_{0}=1.2$. This result suggests that the dark radiation scenario is more
consistent in the fit of the two standard rulers ($k_{\rm{eq}}$ and $r_{s}$).
As a forecast analyses, we fit the two models with a mock $k_{\rm{eq}}$ prior
derived from _Planck_ best-fit $\Lambda$CDM model. Compared with the best-fit
$H_{0}$ in baseline $\Lambda$CDM model, we find $\Delta H_{0}=1.1$ for WZDR
model and $\Delta H_{0}=-2.4$ for EDE model.

##  I Introduction

The Cepheid-calibrated supernovae Ia (SnIa) ($z\lesssim 1.6$) Riess  _et al._
(2021, 2022), the acoustic scale extracted from CMB ($z\sim 1100$) and the BAO
feature measured in LSS surveys Abbott  _et al._ (2018); Alam  _et al._
(2017); Hildebrandt  _et al._ (2020); Hikage  _et al._ (2019) ($z\sim 0.4$)
are amongst the most precise cosmic distance anchors. Mismatches between these
anchors give rise to the well-known $H_{0}$ discrepancy Di Valentino  _et al._
(2021); Perivolaropoulos and Skara (2021a); Abdalla  _et al._ (2022).
Specifically, the measured $H_{0}$ from local distance ladder111The local
distance ladder depends on the choice of local calibrator, i.e., the Cepheid-
calibrated SnIa from SH0ES found $H_{0}=73.04\pm 1.04$Riess  _et al._ (2022),
whilst $H_{0}=69.8\pm 1.9$Freedman (2021) when the local distance ladder is
calibrated by the tip of the red giant branch (TRGB) methods. are
significantly higher than that inferred from the CMB Aghanim  _et al._
(2020a); Chen  _et al._ (2019) and BAO Alam  _et al._ (2017); Gil-Marín  _et
al._ (2018); Bautista  _et al._ (2017), and both of the latter two relies on
the determination of the sound horizon scale at either the last scattering
surface ($r_{s}(z_{*})$) or the drag epoch ($r_{s}(z_{d})$).

Direct solutions to the $H_{0}$ discrepancy from a theoretical point of
view222It is possible that the $H_{0}$ tension is caused by unknown
systematics. Many recent works are focusing on this issue Mortsell  _et al._
(2021); Perivolaropoulos (2022); Perivolaropoulos and Skara (2021b). can
roughly be categorized as as the late-time $(z\lesssim 2)$ Benevento  _et al._
(2020); Alestas  _et al._ (2022); Keeley  _et al._ (2019); Zhou  _et al._
(2022); Alestas  _et al._ (2021a); Bassett  _et al._ (2002); Li and Shafieloo
(2019); Hernández-Almada  _et al._ (2020); Efstathiou (2021); Perivolaropoulos
and Skara (2021b); Vagnozzi (2021); Cai  _et al._ (2022) and the early-time
($z\gtrsim 1100$) deformations Bernal  _et al._ (2016); Poulin  _et al._
(2019); Smith  _et al._ (2020); Moss  _et al._ (2021); Berghaus and Karwal
(2020); Hill and Baxter (2018); Murgia  _et al._ (2021), as well as other
approaches Carneiro  _et al._ (2019); Feng  _et al._ (2022); Alestas  _et al._
(2021b); Marra and Perivolaropoulos (2021); Pandey  _et al._ (2020); Becker
_et al._ (2021); Velten  _et al._ (2021); Harvey  _et al._ (2015); Krall  _et
al._ (2017); Buen-Abad  _et al._ (2018); Pan  _et al._ (2018); Escudero  _et
al._ (2018). The late solutions include a phantom transition of DE which
raises the CMB predicted $H_{0}$ value while preserving the overall shape of
the CMB spectra (including $r_{s}$) Li and Shafieloo (2019); Hernández-Almada
_et al._ (2020); Zhou  _et al._ (2022), however, such modifications would
violate the BAO predictions. To shift both the CMB and BAO anchors at the same
time, one may consider reducing the sound scale at last scattering,
$r_{s}(z_{*})$ by adding extra component, e.g., early dark energy Poulin  _et
al._ (2019); Hill and Baxter (2018); Murgia  _et al._ (2021); Hill  _et al._
(2020); Poulin  _et al._ (2018); Gogoi  _et al._ (2020); Smith  _et al._
(2022) and dark radiation Aloni  _et al._ (2021); Bernal  _et al._ (2016);
Seto and Toda (2021); Dvorkin  _et al._ (2014); Aboubrahim  _et al._ (2022);
Steigman (2013); Ghosh  _et al._ (2022); Vagnozzi (2020). In this case, the
$H_{0}$ tension can still not be fully resolved Jedamzik  _et al._ (2021). One
reason is that the additional component, e.g., early dark energy (EDE), would
bring in new tension with the density fluctuation amplitude, $\sigma_{8}$ Hill
_et al._ (2020). Another is that given a shift in $r_{s}$, the shift in
$H_{0}$ needed to match CMB will be different from the one needed to match the
LSS D’Amico  _et al._ (2021).

Historically, much of the focus has been put on $r_{s}$. To shed light on the
$H_{0}$ discrepancy, additional standard rulers, especially those in the early
universe, are urgently needed. In this work, we focus on the standard ruler
$k^{-1}_{\rm{eq}}$, i.e., the Hubble horizon at matter-radiation equality,
which can be derived from the turnover scale of the Full Shape (FS) power
spectrum. The isolation of $k_{\rm{eq}}$ based information in Full shape (FS)
galaxy power spectrum Philcox  _et al._ (2021); Farren  _et al._ (2022) (or
CMB lensing power spectrum Baxter and Sherwin (2021)) require
$r_{s}$-marginalization operations, which can be achieved either by excluding
the prior on the baryon density $\omega_{b}$ Philcox  _et al._ (2021) or
removing the BAO features in the matter power spectrum via rescaling of the
BAO wiggles Farren  _et al._ (2022). It is shown in Ref. Philcox  _et al._
(2022) that more information can be extracted from the features FS power
spectrum by adopting the latter approach. When combined with CMB lensing and
$\Omega_{m}$ prior from the recent Pantheon+ analysis Brout  _et al._ (2022),
the standard ruler $k_{\rm{eq}}$ is able to deliver precise cosmological
constraints. Especially, $k_{\rm{eq}}$ probes the features of new physics
around the equality $a_{\rm{eq}}$, which is otherwise difficult to detect.

The additional EDE component which peaked around the equality scale
($\log_{10}z_{c}\sim 3.5$) would increase the expansion rate
$H_{\rm{eq}}\equiv H(a_{\rm{eq}})$ as well as the equality scale
$k_{\rm{eq}}\equiv a_{\rm{eq}}H_{\rm{eq}}$, while reducing the sound horizon
scale $r_{s}(z_{*})$. In mock EDE cosmology333Euclid like mock data are
generated assuming EDE model. See Ref. Farren  _et al._ (2022) for detailed
discussions., one would notice a significant peak shift in $H_{0}$ measured
from $k_{\rm{eq}}$\- to $r_{s}$\- based analysis if $\Lambda$CDM model is
assumed Farren  _et al._ (2022). This method provides a novel consistency
check of the $\Lambda$CDM cosmology. On the other hand, the $H_{0}$ value
measured by the two rulers should be consistent as long as the assumed model
is correct. Thus, it is of particular interest that can the modifications to
the early universe shift both of the rulers ($k_{\rm{eq}}$, $r_{s}$) so that
the inferred $H_{0}$ from each of the rulers be consistent with late-time
measurements? Note that the redshift of equality $z_{\rm{eq}}$ is sensitive to
the change in the radiational (matter) component. An addition of radiation
component before the equality will reduce $z_{\rm{eq}}$ and thus shift
$k_{\rm{eq}}$ to smaller value, which is in contrast with the shift brought by
the extra DE component. Given $k_{\rm{eq}}$ (measured in unit of
$h\rm{Mpc}^{-1}$) and a probe of $\Omega_{m}$, one can solve for the Hubble
constant Philcox  _et al._ (2021). We show in Sec. III that the reduction of
$k_{\rm{eq}}$ will lead to a higher inferred value of $H_{0}$ compared with
the best-fit of the $\Lambda$CDM model.

One example of adding the radiational component is the recently investigated
Wess-Zumino Dark Radiation (WZDR) model. The cosmological features of the WZDR
model can be effectively described by a step transition of the effective
number of neutrinos, $N_{\rm{eff}}$, when the interacting dark radiation
density increases as the mediator deposits its entropy into the lighter
species. Compared with the reference self-interacting dark radiation (SIDR)
model, the transition of radiation component in the WZDR model will shift the
position of the high-$\ell$ modes, which enter the horizon before the step. It
is shown in Ref. Aloni  _et al._ (2021) that including such "stepped" dark
radiation does significantly improve the combined fit to CMB, BAO, and SH0ES
data while not degrading the fit to the CMB and BAO.

The outline of this paper is as follows: in Sec. II we take a brief review of
the basic equations associated with the WZDR model. in Sec. III, we show how
EDE and WZDR model affects the standard ruler $k_{\rm{eq}}$. in Sec. IV, we
describe our numerical implementation of both models and the datasets used in
our analysis. The numerical results of the Markov chain Monte Carlo (MCMC)
analysis are presented in Sec. V. The discussion and conclusions are presented
in Sec. VI.

Figure 1: The relative difference between the best-fit parameters
($k_{\rm{eq}},r_{s},z_{\rm{eq}}$) calculated in fiducial $\Lambda$CDM model
and its extensions (EDE, WZDR).

##  II Models

In this section, we present the basic features of the WZDR model Wess and
Zumino (1974) for cosmological implementations. We refer the reader to
Appendix A in ref. Aloni  _et al._ (2021) for detailed discussions. The
interacting dark radiation consists of a Weyl fermion $\psi$ and complex
scalar $\phi$, which interact through a Yukawa coupling,
$\lambda\psi^{2}\phi$. The Lagrangian reads

$\mathcal{L}_{\text{WZ}}=\lambda\phi\psi^{2}+\lambda^{2}(\phi^{*}\phi)^{2}.$
(1)

This model is referred to as the “Wess-Zumino Dark Radiation” (WZDR). When the
temperature of the coupled fluid drops below the scalar mass $m_{\phi}$, the
scalars deposit their entropy into the lighter (massless) $\psi$ species.
According to the entropy conservation, during this process, the lighter
species would be heated, while the energy density of interacting fluid
transitions as $N_{\rm{wzdr}}\rightarrow N_{\rm{wzdr,early}}$ at a redshift
$z_{t}$, where

$N_{\rm{wzdr}}=\frac{\rho_{\rm{wzdr}}}{\rho_{1\nu}}$ (2)

is the effective number of WZDR in the present time, with $\rho_{1\nu}(a)$ the
energy density of a single neutrino in standard $\Lambda$CDM cosmology. In
_Planck_ best-fit WZDR model, the transition of $N_{\rm{dr}}$ happens just
before the matter radiation equality, i.e, $\log_{10}z_{t}\approx 4.3\pm 0.2$,
which implies $m_{\phi}$ is of order of $1\rm{eV}$ or $10\rm{eV}$ The fluid
consists of both massive particle $\phi$ and massless particles $\psi$ during
the transition. After the transition, the massive particle $\phi$ have
annihilated away so that the energy density of fluid in the present time is
dominated by massless particle $\psi$, with

$\rho_{\psi}=\frac{7}{4}\frac{\pi^{2}}{30}T^{4}_{d}=\rho_{1\nu}f^{4}_{T},$ (3)

where $f_{T}\equiv T_{\nu}/T_{\rm{dr}}$ is the temperature ratio between
neutrinos and WZDR. Accordingly, The energy and pressure density of WZDR can
be written as

$\displaystyle\rho_{\rm{wzdr}}\equiv\rho_{\psi}+\rho_{\phi}=\rho_{1\nu}[1+r_{g}\hat{\rho}(x)]f^{4}_{T}$
(4a) $\displaystyle p_{\rm{wzdr}}\equiv
p_{\psi}+p_{\phi}=p_{1\nu}[1+r_{g}\hat{p}(x)]f^{4}_{T}$ (4b)

where we have assumed

$\rho_{\phi}=r_{g}\rho_{\psi}\hat{\rho}(x),\quad
p_{\phi}=r_{g}p_{\psi}\hat{p}(x),$ (5)

with $x\equiv m_{\phi}/T_{\rm{dr}}$, and $r_{g}=8/7$ the density ratio of
bosons to fermions for a single degree of freedom. Note that the dimensionless
integrals $\hat{\rho}$ and $\hat{p}$ can analytically be expressed in terms of
the Bessel functions of the second kind. According to the entropy conservation
in the dark fluid

$S\propto
a^{3}\frac{\rho_{\rm{dr}}+p_{\rm{dr}}}{T_{\rm{dr}}}=\frac{\rho_{\rm{dr,0}}+p_{\rm{dr,0}}}{T_{\rm{dr,0}}},$
(6)

with $\rho_{\rm{dr,0}}=3\>p_{\rm{dr,0}}=N_{\rm{dr,0}}\rho_{1\nu,0}$. Combing
Eq. (3), Eq. (4a, 4b) and Eq. (6) one can solve $x(a)$ numerically, while
$f_{T}$ can be expressed in terms of $x(a)$, i.e., can be expressed as

$\displaystyle f_{T}$
$\displaystyle=f_{T,0}\frac{a}{xa_{t}}=f_{T,0}\left(1+\frac{r_{g}}{4}[3\hat{\rho}(x)+\hat{p}(x)]\right)^{-1/3}$
(7)

where $f_{T,0}=N_{\rm{dr,0}}^{\frac{1}{4}}$, and $a_{t}\equiv
T_{dr,0}/m_{\phi}$ is the transition scale factor.

In this work, we assume one massive neutrino species with a mass of
$0.06\rm{eV}$ and $N_{\rm{ur}}$ to be the effective number of ultra-
relativistic (massless neutrino) species. We use the fractional contribution
of WZDR to total ultra-relativistic species,

$f_{\rm{wzdr}}\equiv\frac{N_{\rm{wzdr}}}{N_{\rm{wzdr}}+N_{\rm{ur}}},$ (8)

to parameterize the WZDR model (denote as (f)WZDR). This is analogous to the
EDE model, which is parameterized by maximal fractional contribution of the
EDE component to the energy density of the universe,
$f_{\rm{EDE}}\equiv\rho_{\rm{EDE}}(z_{c})/\rho_{\rm{tot}}(z_{c})$ and the
critical redshift $z_{c}$. The WZDR model can also be effectively described by
the effective number of additional neutrino species $N_{\rm{wzdr}}$
(equivalent to $N_{\rm{IR}}$ in Ref. Aloni  _et al._ (2021)). See Appendix A
for the cosmological constraints on both of the parameterizations described
above444$N_{\rm{ur}}$ can either a free parameter or be fixed at given value.
In (f)WZDR parameterization, we find $N_{\rm{ur}}=2.00\pm 0.27$, which closely
matches the value given in standard cosmology. Consequently, we have fixed
$N_{\rm{ur}}$ to 2.0328 in the following analysis..

##  III Impact on Standard Rulers

Figure 2: The derived $H_{0}$ in EDE and (f)WZDR scenarios with $k_{\rm{eq}}$
and $\Omega_{m}$ fixed at _Planck_ best-fit value assuming $\Lambda$CDM model.

In $\Lambda$CDM model, the equality scale (in unit of $\rm{Mpc}^{-1}$) can
simply be approximated as

$\displaystyle k_{\rm{eq}}$
$\displaystyle=(2\Omega_{m}H^{2}_{0}z_{\rm{eq}})^{1/2},$ (9a) $\displaystyle
z_{\rm{eq}}$ $\displaystyle=2.5\times 10^{4}\Omega_{m}h^{2}\Theta^{-4}_{2.7}$
(9b)

with $\Omega_{m}$ the matter density fraction and $\Theta_{2.7}\equiv
T_{\rm{CMB}}/(2.7K)$ the CMB photon temperature. The addition of DE component
in the early time does not change the ratio between matter and radiation. Thus
Eq. (9b) still holds in the EDE model (see the middle panel of Fig.1).
However, Eq. (9a) is no longer valid in this scenario, because the the extra
energy component increases the expansion rate $H(z_{\rm{eq}})$ while reduces
the comoving hobble horizon at equality scale ($k^{-1}_{\rm{eq}}$). On the
contrary, Eq. (9a) is valid in WZDR scenario whereas Eq. (9b) is not. The
additional radiation component $N_{\rm{wzdr}}$ shifts $z_{\rm{eq}}$ to smaller
value, and $k_{\rm{eq}}$ is reduced as a consequence. As is shown in Fig. 1,
the addition of EDE and WZDR shifts the standard ruler $k_{\rm{eq}}$
differently. While the sound horizon $r_{s}(z_{*})$ is reduced in both
scenarios, which increases the inferred value of $H_{0}$ in $r_{s}$-based
analysis.

$H_{0}$ is constrained by measuring the angular scale of the cosmological
horizon at matter-radiation-equality, i.e. $k_{\rm{eq}}/D_{A}(z)\propto
k_{\rm{eq}}/h$ Philcox  _et al._ (2022), which is equivalent to the value of
$k_{\rm{eq}}$ in unit of $h\>\rm{Mpc}^{-1}$ (If not stated otherwise,
$k_{\rm{eq}}$ discussed in this paper is in unit of $h\>\rm{Mpc}^{-1}$). As
can be seen in Eq. (9a, 9b), given $k_{\rm{eq}}$ and a probe of $\Omega_{m}$,
one can solve for the Hubble constantPhilcox  _et al._ (2021). We have
discussed how different models shift standard ruler $k_{\rm{eq}}$ compared
with the fiducial $\Lambda$CDM model. Next, we fix $k_{\rm{eq}}$ and
$\Omega_{m}$ to the value inferred in _Planck_ best-fit $\Lambda$CDM
model555The baseline $\Lambda$CDM cosmology
$\\{h=0.6821,\omega_{c}=0.1177,\omega_{b}=0.02253\\}$. We refer the readers to
Section III A in Ref. Farren  _et al._ (2022) for detailed discussions. ,i.e.,
$k_{\rm{eq}}=1.545\>h\rm{Mpc}^{-1},\Omega_{m}=0.311$, to see the best-fit of
$H_{0}$ in these scenarios by “shooting” the targeted $k_{\rm{eq}}$ value. We
show in Fig. 2 the best-fit $H_{0}$ value as a function of the fractional
density of EDE and WZDR. As can be seen, the increase of $k_{\rm{eq}}$ by EDE
component can be compensated by a lower inferred value of $H_{0}$, while the
decrease of $k_{\rm{eq}}$ is compensated by a higher inferred value of
$H_{0}$. We will further verify this point by adopting a mock prior on
$k_{\rm{eq}}$ in MCMC analysis (See Section V.2).

Dataset | FS without $r_{s}$ |  $k_{\rm{eq}}$ prior
---|---|---
Model |  (f)WZDR | EDE | (f)WZDR | EDE
$H_{0}$ | $69.7\pm 1.1$ | $70.6\pm 1.0\quad$ | $70.3\pm 2.4$ | $65.8^{+2.0}_{-2.3}$
$10^{2}k_{\rm{eq}}\quad$ | $1.517^{+0.055}_{-0.050}$ | $1.673^{+0.048}_{-0.055}\quad$ | $\quad 1.545\pm 0.020$ | $1.545\pm 0.021$
$S_{8}$ | $0.809\pm 0.027$ | $0.820^{+0.025}_{-0.029}\quad$ | $0.825\pm 0.035$ | $0.725^{+0.026}_{-0.030}$
$\Omega{}_{m}$ | $0.310^{+0.011}_{-0.0094}$ | $0.3129\pm 0.0088\quad$ | $0.310^{+0.011}_{-0.0094}$ | $0.3129\pm 0.0088$
$\chi^{2}_{\rm{min}}$ | $1072.06$ | $1073.42$ | $0.06$ | $0.08$
Table 1: The the mean $\pm 1\sigma$ constraints on the derived parameters in
EDE and $f$WZDR model, as inferred from base datasets $\mathcal{B}$ \+ FS
power spectrum with $r_{s}$ marginalized over (left panel), and
$\mathcal{B^{\prime}}$ \+ $k_{\rm{eq}}$ prior derived from _Planck_ best-fit.
Upper and lower bounds correspond to the $68\%$ C.L. interval.

##  IV datasets and numerical method

Following in Ref. Philcox  _et al._ (2022), we consider the following
combination as our base datasets (denote as $\mathcal{B}$), which includes:

  * •

Supernovae (SNe): The derived $k_{\rm{eq}}/h$ (in unit of $\rm{Mpc}^{-1}$) is
degenerate with the matter density, i.e.,
$k_{\rm{eq}}/h\propto\Omega_{m}^{1/2}h$ (see Eq. 9a, and Eq. 9b), thus, the
additional priors on $\Omega_{m}$ can improve the constraints on $H_{0}$ from
the equality scale. We consider the Gaussian prior $\Omega_{m}=0.338\pm 0.018$
from the most recent Pantheon analysis Aboubrahim  _et al._ (2022).

  * •

Big Bang Nucleosynthesis (BBN): Including the $\omega_{b}$ prior from BBN
breaks the $H_{0}-\omega_{b}$ degeneracy and will strengthen the sources of
$H_{0}$ information from the equality scale. Consequently, a prior on the
physical baryon density $\omega_{b}=0.02268\pm 0.00036$ is imposed to maximize
the extracted information, following [28].

  * •

Spectral Index and Primordial Amplitude: Following from previous works [41],
we add a flattened Gaussian prior on $A_{s}$ and $n_{s}$ centred at _Planck_
best-fit Aghanim  _et al._ (2020a) i.e., $n_{s}=0.965\pm 0.02$, and
$A_{s}=(2.11\pm 0.108)\times 10^{-9}$.

The base datasets $\mathcal{B}$ can optionally be combined with:

  * •

Galaxy Power Spectra: We adopt the method described in Ref. Farren  _et al._
(2022) to marginalized over $r_{s}$ information in galaxy power spectrum. This
method can be applied to the current version of the BOSS DR12 galaxy survey
dataset Beutler  _et al._ (2017), which is part of SDSS-III Alam  _et al._
(2015); Dawson  _et al._ (2013). This data split across two redshift bins (at
z = 0.38, 0.61) in each of the Northern and Southern galactic caps (NGC and
SGC). Following from previous work, we utilize the CLASS-PT implementation
Chudaykin  _et al._ (2020) with publicly available likelihoods666Available at
https://github.com/oliverphilcox/full_shape_likelihoods., in which the power
spectrum are modelled with the Effective Field Theory of Large Scale Structure
Ivanov  _et al._ (2020) including one-loop perturbation theory and etc.

  * •

Baryon Acoustic Oscillation: measurements of the BAO signal at $z=0.38,0.51$
and $0.61$ from BOSS DR12 galaxy power spectrum Alam  _et al._ (2017).

  * •

$S_{8}$ prior: We adopt the prior $S_{8}=0.800\pm 0.029$ from KiDS-450+GAMA
van Uitert  _et al._ (2018), similar results can be found in the the
measurement from the first-year data of HSC SSP20, for which $S_{8}=0.804\pm
0.032$ Hamana  _et al._ (2020) .

We are interested in how the extra radiation (DE) component affects the Hubble
constant inferred from the standard ruler $k_{\rm{eq}}$. Hence, throughout the
analysis we fix the fractional density of additional components to best-fit
value obtained in cosmological analysis (see Appendix A for a detailed
discussion of the best-fit (f)WZDR model), i.e., $f_{\rm{wzdr}}=0.21$ and
$f_{\rm{EDE}}=0.122$777We adopt the best-fit values of EDE model from Table I
of Ref. Smith  _et al._ (2020). We sample over the following set of
cosmological parameters:

$\\{H_{0},\omega_{\rm b},\omega_{\rm cdm},\log 10^{10}A_{s},n_{s},\sum
m_{\nu}\\}.$ (10)

The optical depth of reionization, $\tau_{\rm reio}$, is fixed to $0.055$. We
implement the EDE and the WZDR scenarios as modifications to the publicly
available Einstein-Boltzmann code CLASS Lesgourgues (2011); Blas  _et al._
(2011) package. The non-linear matter power spectrum required by redshift-
space distortion (RSD) likelihoods are computed using the “HMcode” Mead  _et
al._ (2015); Mead (2017); Mead  _et al._ (2020) implemented in CLASS. The MCMC
analyses are performed using the publicly available code MontePython
Brinckmann and Lesgourgues (2019).

##  V numerical results

(a) FS-$r_{s}$

(b) $\>k_{\rm{eq}}$

Figure 3: Cosmological parameters constraints for EDE model (blue contour) and
(f)WZDR model (red contour) from the dataset: $\mathcal{B}$ \+ FS power
spectrum with $r_{s}$ marginalized over (denoting as FS - $r_{s}$ in the left
panel), and fitting to $\mathcal{B^{\prime}}$ \+ $k_{\rm{eq}}$ prior derived
from _Planck_ best-fit $\Lambda$CDM model (right panel).

We show in Fig. 3(a) and Tab. 1 the $r_{s}$-independent constraints from FS
galaxy power spectra, supplemented by the base datasets $\mathcal{B}$ which
includes priors on $\Omega_{m}$ from Pantheon+ and $\omega_{b}$ from BBN.

For comparison between constraints with and without $r_{s}$ information, we
show in Table 2 the parameter constraints obtained from the base datasets
$\mathcal{B}$ combined with three choices of likelihoods: (1), the FS
likelihoods with $r_{s}$ marginalized over, (2), the FS likelihoods including
the BAO information, and (3), the BAO likelihoods combined with $S_{8}$ prior
derived from lensing experiments. The triangular plot with the 1D posterior
distributions and the 2D contour plots for the constrained parameters are
shown in Fig. 4 and Fig. 5.

###  V.1 Derived equality scale

In the discussion of Sec. III, we have shown that the addition of the WZDR
reduces $k_{\rm{eq}}$ while the opposite is true for the EDE component. As is
shown in Table 2, the $1\sigma$ constraint of $k_{\rm{eq}}$ from $r_{s}$-based
analysis (BAO+$S_{8}$) in (f)WZDR model, $k_{\rm{eq}}=1.512\pm 0.041$, is
$2.4\sigma$ lower than that in EDE model, $k_{\rm{eq}}=1.659\pm 0.045$. Still,
there is a $2.2\sigma$ tension between the EDE and (f)WZDR best-fit
$k_{\rm{eq}}$ when fitting with the $k_{\rm{eq}}$-based analysis
($r_{s}$-marginalized FS power spectra). While in each model, $k_{\rm{eq}}$
derived from the BAO+$S_{8}$ prior is close to that from the
$r_{s}$-marginalized FS power spectrum. These results suggest that
$k_{\rm{eq}}$ can not yet be tightly (model independently) constrained by the
Equality-based measurements. Currently, the information extracted from the
$r_{s}$-marginalized FS power spectrum is dominated by information concerning
structure growth ($f\sigma_{8}$, $S_{8}$) rather than the angular scale of the
cosmological horizon at matter-radiation-equality. Otherwise, the derived
$k_{\rm{eq}}$ in $k_{\rm{eq}}$-based analysis from different models should be
more consistent.

###  V.2 $k_{\rm{eq}}$\- and $r_{s}$-based inferrance of $H_{0}$

In WZDR model, we find a minor peak shift of $DeltaH_{0}=0.2$ between the fit
of FS power spectrum with and without $r_{s}$ information, and a peak shift of
$DeltaH_{0}=0.4$ from the $k_{\rm{eq}}$\- to $r_{s}$-based analyses (see Table
2). While in EDE model we have $DeltaH_{0}=1.0$ between FS and FS - $r_{s}$,
and $DeltaH_{0}=1.2$ between $k_{\rm{eq}}$\- and $r_{s}$-based analyses. Since
the addition of radiation component shifts the two standard rulers more
coherently (see the discussion in Sec. III), the best-fit $H_{0}$ from
$r_{s}$\- and $k_{\rm{eq}}-$based analyses in the WZDR model are more
consistent than that in the EDE model. However, the above results are
inconclusive since the Hubble constant measured by the two rulers is
consistent (within $2\sigma$) for both EDE and WZDR models. As is discussed
above, $k_{\rm{eq}}$-based datasets are dominated by the information
$f\sigma_{8}$ ($S_{8}$) rather than the angular scale of the cosmological
horizon at matter-radiation-equality, which makes it difficult for us to
detect the inconsistency.

To see how the standard ruler $k_{\rm{eq}}$ affects the infferance of $H_{0}$,
we replace the $r_{s}$-marginalized FS power spectra with a mock prior:
$k_{\rm{eq}}=1.545\pm 0.020\times 10^{-2}h\rm{Mpc}^{-1}$ derived from _Planck_
best-fit $\Lambda$CDM model (Accordingly, the $\Omega_{m}$ prior in the
supplement datasets $\mathcal{B}$ are substituted by the emphPlanck best-fit
value, $\Omega_{m}$, denoting as $\mathcal{B^{\prime}}$). The results can be
found in the right panel of Tab. 1. As is discussed in Sec. III, to fit with
the standard rulers fixed by the fiducial $\Lambda$CDM model, the increase of
$k_{\rm{eq}}$ should be compensated by a lower $H_{0}$. Accordingly, its
decrease should be compensated by a higher $H_{0}$. In this case, the inferred
value of $H_{0}$ in the WZDR model is significantly higher than in the EDE
model. Note that we have also adopted a tight prior on $\Omega_{m}$, so higher
$H_{0}$ will result in higher $\omega_{c}=\Omega_{m}h^{2}-\omega_{b}$ and
higher $S_{8}$. The inferred $S_{8}$ in the EDE model is much lesser than that
of the WZDR model.

Model | Dataset | $H_{0}$ | $10^{2}k_{\rm{eq}}$ | $S_{8}$ | $\Omega_{m}$ | $\ln 10^{10}A_{s}$ | $n_{s}$ | $\chi^{2}_{\rm{min}}$
---|---|---|---|---|---|---|---|---
(f)WZDR |  $\mathcal{B}$+FS - $r_{s}$ | $69.7\pm 1.1$ | $1.517^{+0.055}_{-0.050}$ | $0.809\pm 0.027$ | $0.310^{+0.011}_{-0.0094}$ | $3.024^{+0.039}_{-0.031}$ | $0.973\pm 0.019$ | 1072.06
|  $\mathcal{B}$+FS | $69.90^{+0.90}_{-0.79}$ | $1.515^{+0.047}_{-0.029}$ | $0.817^{+0.020}_{-0.016}$ | $0.308^{+0.011}_{-0.0042}$ | $3.036^{+0.023}_{-0.034}$ | $0.971^{+0.021}_{-0.026}$ | 1072.26
|  $\mathcal{B}$+BAO+$S_{8}$ | $70.13\pm 0.96$ | $1.512\pm 0.041$ | $0.813^{+0.019}_{-0.021}$ | $0.3042\pm 0.0082$ | $3.037\pm 0.033$ | $0.963^{+0.016}_{-0.022}$ | 7.78
EDE |  $\mathcal{B}$+FS - $r_{s}$ | $70.6\pm 1.0$ | $1.673^{+0.048}_{-0.055}$ | $0.820^{+0.025}_{-0.029}$ | $0.3129\pm 0.0088$ | $3.030\pm 0.036$ | $0.975\pm 0.020$ | 1073.42
|  $\mathcal{B}$+FS | $71.6\pm 0.9$ | $1.692^{+0.053}_{-0.044}$ | $0.820\pm 0.023$ | $0.3134\pm 0.0085$ | $3.000^{+0.035}_{-0.042}$ | $0.969^{+0.018}_{-0.020}$ | 1073.58
|  $\mathcal{B}$+BAO+$S_{8}$ | $71.8\pm 1.1$ | $1.659\pm 0.045$ | $0.819\pm 0.022$ | $0.3043^{+0.0076}_{-0.0091}$ | $3.032\pm 0.040$ | $0.964^{+0.015}_{-0.020}$ | 8.18
Table 2: The mean and $1$-$\sigma$ constraints on cosmological parameters in
the (f)WZDR model (top panel) from the base datasets $\mathcal{B}$ combined
with three choices of likelihoods: (1), the FS likelihoods with $r_{s}$
marginalized over (denoting as FS - $r_{s}$), (2), the FS likelihoods
including the BAO information, and (3), the BAO likelihoods combined with
$S_{8}$ prior derived from lensing experiments. Upper and lower bounds
correspond to the $68\%$ C.L. interval. Figure 4: Constraints on cosmological
parameters in the WZDR scenario from the base datasets $\mathcal{B}$ combined
with three choices of likelihoods: (1), the FS likelihoods with $r_{s}$
marginalized over, (2), the FS likelihoods including the BAO information, and
(3), the BAO likelihoods combined with $S_{8}$ prior derived from lensing
experiments. Figure 5: Constraints on cosmological parameters in the EDE
scenario from three choices of likelihoods discussed above.

##  VI CONCLUSIONS

The Hubble horizon at matter-radiation equality, $k^{-1}_{\rm{eq}}$, provides
valuable test of new physics around (prior to) $z_{\rm{eq}}$, while the
commonly used standard ruler, $r_{s}$, is sensitive to physics prior to the
last scattering. These two rulers provide an interesting consistency check for
the $\Lambda$CDM model and its extensions, i.e., the inferred value of $H_{0}$
from both the rulers should be consistent. We find that the addition of EDE
component which peaked around the equality scale ($\log_{10}z_{c}\sim 3.5$)
increases equality scale, $k_{\rm{eq}}$. On the contray, the addition of
radiational component before the equality reduces $z_{\rm{eq}}$ and
$k_{\rm{eq}}$. We show in Fig. 2 that the increase of $k_{\rm{eq}}$ by EDE
component can be compensated by an upward shift of $H_{0}$, while the decrease
of $k_{\rm{eq}}$ is compensated by a higher inferred value of $H_{0}$. The two
standard rulers should be shifted coherently to relieve the $H_{0}$
discrepancy between the early and late universe.

Our numerical results shows that, the best-fit $H_{0}$ values obtained from
$r_{s}$\- and $k_{\rm{eq}}-$based analyses are consistent within $1\sigma$ in
(f)WZDR model ($\Delta H_{0}=0.4$). In the same analysis assuming the EDE
model, we observed a peck-shift of $\Delta H_{0}=1.2$. There is a significant
tension between the derived $k_{\rm{eq}}$ in EDE and the (f)WZDR model.
Meanwhile, in each model, $k_{\rm{eq}}$ derived from the BAO+$S_{8}$ prior is
close to that from the $r_{s}$-marginalized FS power spectrum. These results
suggest that, currently, the $r_{s}$ marginalized FS power can not tightly
contain the standard ruler $k_{\rm{eq}}$. The extracted signal consists mainly
of the information concerning the structure growth ($f\sigma_{8}$, $S_{8}$)
amplitude rather than the angular scale of the cosmological horizon at matter-
radiation-equality.

We expect more information on $k_{\rm{eq}}$ to be extracted from next
generation of galaxy survey, e.g., Euclid, SKA. As a forecast analysis, we
adopt a mock prior derived from _Planck_ best-fit $\Lambda$CDM model to see
how the standard ruler $k_{\rm{eq}}$ affects the inference of $H_{0}$. In WZDR
model, we find $\Delta H_{0}=1.1$ compared with the best-fit $H_{0}$ in
baseline $\Lambda$CDM model, while in the EDE model the pair difference is
$\Delta H_{0}=-2.4$. This result further suggests that the two standard rulers
should be reduced simultaneously to relieve the tension with the local
distance ladder.

###### Acknowledgements.

This work is supported in part by National Natural Science Foundation of China
under Grant No.12075042, Grant No.11675032 (People’s Republic of China), Grant
No.12175095, and supported by LiaoNing Revitalization Talents Program
(XLYC2007047).

##  Appendix A Best-fit WZDR

We repeat the cosmological constraints on the WZDR model with the same
datasets discussed in Ref. Schöneberg  _et al._ (2021). Meanwhile, we adopt a
new parameterization with $N_{\rm{ur}}$ being a free parameter (denote as
(f)WZDR). This datasets includes: the Planck _Planck_ 2018 low-$\ell$ TT+EE
and high-$\ell$ TT+TE+EE temperature and polarization power spectrum Aghanim
_et al._ (2020b, a), as well as CMB lensing Aghanim  _et al._ (2020a), the BAO
measurements Alam  _et al._ (2017); Beutler  _et al._ (2011); Ross  _et al._
(2015), the PANTHEON supernovae data Scolnic  _et al._ (2018), as well as
measurements from SH0ES Riess  _et al._ (2022) via a prior on the intrinsic
magnitude of supernovae $M_{b}$ Camarena and Marra (2021). Constraints at
$68\%$ C.L. on the cosmological parameters and $\chi^{2}$ statistics can be
found in Table 3. The posterior distributions and for WZDR and (f)WZDR
parameterization are shown in Fig. 6 and 7, respectively.

Parameter | WZDR | (f)WZDR
---|---|---
$10^{-2}\omega{}_{b}$ | $2.286\pm 0.018$ | $2.279\pm 0.016$
$\omega{}_{c}$ | $0.1263^{+0.0041}_{-0.0033}$ | $0.1256^{+0.0026}_{-0.0029}$
$\ln 10^{10}A_{s}$ | $3.046\pm 0.014$ | $1.04241\pm 0.00086$
$n_{s}$ | $0.9725^{+0.0036}_{-0.0040}$ | $3.047^{+0.021}_{-0.018}$
$\tau_{\rm{reio}}$ | $0.0585^{+0.0059}_{-0.0069}$ | $0.970^{+0.010}_{-0.0092}$
$N_{\rm{idr}}$ | $0.55^{+0.26}_{-0.19}$ | $-$
$f_{\mathrm{wzdr}}$ | $-$ | $0.21^{+0.11}_{-0.13}$
$N_{\mathrm{ur}}$ | $(2.0328)$ | $2.00\pm 0.27$
$log10z_{t}$ | $4.26^{+0.18}_{-0.20}$ | $4.29^{+0.16}_{-0.19}$
$H_{0}$ | $71.0^{+1.3}_{-0.87}$ | $70.68\pm 0.89$
$S_{8}$ | $0.8185\pm 0.0059$ | $0.818^{+0.011}_{-0.010}$
$N_{\rm{eff}}$ | $3.47^{+0.20}_{-0.15}$ | $3.43^{+0.14}_{-0.16}$
Table 3: The the mean $\pm 1\sigma$ constraints on the cosmological parameters
and derived parameters in WZDR and (f)WZDR model from the datasets discussed
in Appendix A. Upper and lower bounds correspond to the $68\%$ C.L. interval.
Figure 6: Cosmological parameters constraints in WZDR model (with
$N_{\rm{wzdr}}$ being the only extended parameter) from the datasets discussed
in Appendix A. Figure 7: Cosmological parameters constraints in (f)WZDR model
(with $f_{\rm{wzdr}}$ and $N_{\rm{ur}}$ being free parameters) from the
datasets discussed in Appendix A.

## References

  * Riess  _et al._ (2021) A. G. Riess, S. Casertano, W. Yuan, J. B. Bowers, L. Macri, J. C. Zinn, and D. Scolnic, Astrophys. J. Lett. 908, L6 (2021), arXiv:2012.08534 [astro-ph.CO] .
  * Riess  _et al._ (2022) A. G. Riess _et al._ , Astrophys. J. Lett. 934, L7 (2022), arXiv:2112.04510 [astro-ph.CO] .
  * Abbott  _et al._ (2018) T. Abbott _et al._ (DES), Phys. Rev. D 98, 043526 (2018), arXiv:1708.01530 [astro-ph.CO] .
  * Alam  _et al._ (2017) S. Alam _et al._ (BOSS), Mon. Not. Roy. Astron. Soc. 470, 2617 (2017), arXiv:1607.03155 [astro-ph.CO] .
  * Hildebrandt  _et al._ (2020) H. Hildebrandt _et al._ , Astron. Astrophys. 633, A69 (2020), arXiv:1812.06076 [astro-ph.CO] .
  * Hikage  _et al._ (2019) C. Hikage _et al._ (HSC), Publ. Astron. Soc. Jap. 71, Publications of the Astronomical Society of Japan, Volume 71, Issue 2, April 2019, 43, https://doi.org/10.1093/pasj/psz010 (2019), arXiv:1809.09148 [astro-ph.CO] .
  * Di Valentino  _et al._ (2021) E. Di Valentino, O. Mena, S. Pan, L. Visinelli, W. Yang, A. Melchiorri, D. F. Mota, A. G. Riess, and J. Silk, (2021), arXiv:2103.01183 [astro-ph.CO] .
  * Perivolaropoulos and Skara (2021a) L. Perivolaropoulos and F. Skara, (2021a), arXiv:2105.05208 [astro-ph.CO] .
  * Abdalla  _et al._ (2022) E. Abdalla _et al._ , JHEAp 34, 49 (2022), arXiv:2203.06142 [astro-ph.CO] .
  * Freedman (2021) W. L. Freedman, (2021), arXiv:2106.15656 [astro-ph.CO] .
  * Aghanim  _et al._ (2020a) N. Aghanim _et al._ (Planck), Astron. Astrophys. 641, A6 (2020a), [Erratum: Astron.Astrophys. 652, C4 (2021)], arXiv:1807.06209 [astro-ph.CO] .
  * Chen  _et al._ (2019) L. Chen, Q.-G. Huang, and K. Wang, JCAP 02, 028 (2019), arXiv:1808.05724 [astro-ph.CO] .
  * Gil-Marín  _et al._ (2018) H. Gil-Marín _et al._ , Mon. Not. Roy. Astron. Soc. 477, 1604 (2018), arXiv:1801.02689 [astro-ph.CO] .
  * Bautista  _et al._ (2017) J. E. Bautista _et al._ , Astron. Astrophys. 603, A12 (2017), arXiv:1702.00176 [astro-ph.CO] .
  * Mortsell  _et al._ (2021) E. Mortsell, A. Goobar, J. Johansson, and S. Dhawan, (2021), arXiv:2105.11461 [astro-ph.CO] .
  * Perivolaropoulos (2022) L. Perivolaropoulos, Universe 8, 263 (2022), arXiv:2201.08997 [astro-ph.EP] .
  * Perivolaropoulos and Skara (2021b) L. Perivolaropoulos and F. Skara, Phys. Rev. D 104, 123511 (2021b), arXiv:2109.04406 [astro-ph.CO] .
  * Benevento  _et al._ (2020) G. Benevento, W. Hu, and M. Raveri, Phys. Rev. D 101, 103517 (2020), arXiv:2002.11707 [astro-ph.CO] .
  * Alestas  _et al._ (2022) G. Alestas, D. Camarena, E. Di Valentino, L. Kazantzidis, V. Marra, S. Nesseris, and L. Perivolaropoulos, Phys. Rev. D 105, 063538 (2022), arXiv:2110.04336 [astro-ph.CO] .
  * Keeley  _et al._ (2019) R. E. Keeley, S. Joudaki, M. Kaplinghat, and D. Kirkby, JCAP 12, 035 (2019), arXiv:1905.10198 [astro-ph.CO] .
  * Zhou  _et al._ (2022) Z. Zhou, G. Liu, Y. Mu, and L. Xu, Mon. Not. Roy. Astron. Soc. 511, 595 (2022), arXiv:2105.04258 [astro-ph.CO] .
  * Alestas  _et al._ (2021a) G. Alestas, I. Antoniou, and L. Perivolaropoulos, Universe 7, 366 (2021a), arXiv:2104.14481 [astro-ph.CO] .
  * Bassett  _et al._ (2002) B. A. Bassett, M. Kunz, J. Silk, and C. Ungarelli, Mon. Not. Roy. Astron. Soc. 336, 1217 (2002), arXiv:astro-ph/0203383 .
  * Li and Shafieloo (2019) X. Li and A. Shafieloo, Astrophys. J. Lett. 883, L3 (2019), arXiv:1906.08275 [astro-ph.CO] .
  * Hernández-Almada  _et al._ (2020) A. Hernández-Almada, G. Leon, J. Magaña, M. A. García-Aspeitia, and V. Motta, Mon. Not. Roy. Astron. Soc. 497, 1590 (2020), arXiv:2002.12881 [astro-ph.CO] .
  * Efstathiou (2021) G. Efstathiou, Mon. Not. Roy. Astron. Soc. 505, 3866 (2021), arXiv:2103.08723 [astro-ph.CO] .
  * Vagnozzi (2021) S. Vagnozzi, Phys. Rev. D 104, 063524 (2021), arXiv:2105.10425 [astro-ph.CO] .
  * Cai  _et al._ (2022) R.-G. Cai, Z.-K. Guo, S.-J. Wang, W.-W. Yu, and Y. Zhou, (2022), arXiv:2202.12214 [astro-ph.CO] .
  * Bernal  _et al._ (2016) J. L. Bernal, L. Verde, and A. G. Riess, JCAP 10, 019 (2016), arXiv:1607.05617 [astro-ph.CO] .
  * Poulin  _et al._ (2019) V. Poulin, T. L. Smith, T. Karwal, and M. Kamionkowski, Phys. Rev. Lett. 122, 221301 (2019), arXiv:1811.04083 [astro-ph.CO] .
  * Smith  _et al._ (2020) T. L. Smith, V. Poulin, and M. A. Amin, Phys. Rev. D 101, 063523 (2020), arXiv:1908.06995 [astro-ph.CO] .
  * Moss  _et al._ (2021) A. Moss, E. Copeland, S. Bamford, and T. Clarke, (2021), arXiv:2109.14848 [astro-ph.CO] .
  * Berghaus and Karwal (2020) K. V. Berghaus and T. Karwal, Phys. Rev. D 101, 083537 (2020), arXiv:1911.06281 [astro-ph.CO] .
  * Hill and Baxter (2018) J. C. Hill and E. J. Baxter, JCAP 08, 037 (2018), arXiv:1803.07555 [astro-ph.CO] .
  * Murgia  _et al._ (2021) R. Murgia, G. F. Abellán, and V. Poulin, Phys. Rev. D 103, 063502 (2021), arXiv:2009.10733 [astro-ph.CO] .
  * Carneiro  _et al._ (2019) S. Carneiro, P. C. de Holanda, C. Pigozzo, and F. Sobreira, Phys. Rev. D 100, 023505 (2019), arXiv:1812.06064 [astro-ph.CO] .
  * Feng  _et al._ (2022) L. Feng, R.-Y. Guo, J.-F. Zhang, and X. Zhang, Phys. Lett. B 827, 136940 (2022), arXiv:2109.06111 [astro-ph.CO] .
  * Alestas  _et al._ (2021b) G. Alestas, L. Kazantzidis, and L. Perivolaropoulos, Phys. Rev. D 103, 083517 (2021b), arXiv:2012.13932 [astro-ph.CO] .
  * Marra and Perivolaropoulos (2021) V. Marra and L. Perivolaropoulos, Phys. Rev. D 104, L021303 (2021), arXiv:2102.06012 [astro-ph.CO] .
  * Pandey  _et al._ (2020) K. L. Pandey, T. Karwal, and S. Das, JCAP 07, 026 (2020), arXiv:1902.10636 [astro-ph.CO] .
  * Becker  _et al._ (2021) N. Becker, D. C. Hooper, F. Kahlhoefer, J. Lesgourgues, and N. Schöneberg, JCAP 02, 019 (2021), arXiv:2010.04074 [astro-ph.CO] .
  * Velten  _et al._ (2021) H. Velten, I. Costa, and W. Zimdahl, (2021), arXiv:2104.05352 [astro-ph.CO] .
  * Harvey  _et al._ (2015) D. Harvey, R. Massey, T. Kitching, A. Taylor, and E. Tittley, Science 347, 1462 (2015), arXiv:1503.07675 [astro-ph.CO] .
  * Krall  _et al._ (2017) R. Krall, F.-Y. Cyr-Racine, and C. Dvorkin, JCAP 09, 003 (2017), arXiv:1705.08894 [astro-ph.CO] .
  * Buen-Abad  _et al._ (2018) M. A. Buen-Abad, M. Schmaltz, J. Lesgourgues, and T. Brinckmann, JCAP 01, 008 (2018), arXiv:1708.09406 [astro-ph.CO] .
  * Pan  _et al._ (2018) Z. Pan, M. Kaplinghat, and L. Knox, Phys. Rev. D 97, 103531 (2018), arXiv:1801.07348 [astro-ph.CO] .
  * Escudero  _et al._ (2018) M. Escudero, L. Lopez-Honorez, O. Mena, S. Palomares-Ruiz, and P. Villanueva-Domingo, JCAP 06, 007 (2018), arXiv:1803.08427 [astro-ph.CO] .
  * Hill  _et al._ (2020) J. C. Hill, E. McDonough, M. W. Toomey, and S. Alexander, Phys. Rev. D 102, 043507 (2020), arXiv:2003.07355 [astro-ph.CO] .
  * Poulin  _et al._ (2018) V. Poulin, T. L. Smith, D. Grin, T. Karwal, and M. Kamionkowski, Phys. Rev. D 98, 083525 (2018), arXiv:1806.10608 [astro-ph.CO] .
  * Gogoi  _et al._ (2020) A. Gogoi, R. K. Sharma, P. Chanda, and S. Das, (2020), arXiv:2005.11889 [astro-ph.CO] .
  * Smith  _et al._ (2022) T. L. Smith, M. Lucca, V. Poulin, G. F. Abellan, L. Balkenhol, K. Benabed, S. Galli, and R. Murgia, (2022), arXiv:2202.09379 [astro-ph.CO] .
  * Aloni  _et al._ (2021) D. Aloni, A. Berlin, M. Joseph, M. Schmaltz, and N. Weiner, (2021), arXiv:2111.00014 [astro-ph.CO] .
  * Seto and Toda (2021) O. Seto and Y. Toda, Phys. Rev. D 103, 123501 (2021), arXiv:2101.03740 [astro-ph.CO] .
  * Dvorkin  _et al._ (2014) C. Dvorkin, M. Wyman, D. H. Rudd, and W. Hu, Phys. Rev. D 90, 083503 (2014), arXiv:1403.8049 [astro-ph.CO] .
  * Aboubrahim  _et al._ (2022) A. Aboubrahim, M. Klasen, and P. Nath, JCAP 04, 042 (2022), arXiv:2202.04453 [astro-ph.CO] .
  * Steigman (2013) G. Steigman, Phys. Rev. D 87, 103517 (2013), arXiv:1303.0049 [astro-ph.CO] .
  * Ghosh  _et al._ (2022) S. Ghosh, S. Kumar, and Y. Tsai, JCAP 05, 014 (2022), arXiv:2107.09076 [astro-ph.CO] .
  * Vagnozzi (2020) S. Vagnozzi, Phys. Rev. D 102, 023518 (2020), arXiv:1907.07569 [astro-ph.CO] .
  * Jedamzik  _et al._ (2021) K. Jedamzik, L. Pogosian, and G.-B. Zhao, Commun. in Phys. 4, 123 (2021), arXiv:2010.04158 [astro-ph.CO] .
  * D’Amico  _et al._ (2021) G. D’Amico, L. Senatore, P. Zhang, and H. Zheng, JCAP 05, 072 (2021), arXiv:2006.12420 [astro-ph.CO] .
  * Philcox  _et al._ (2021) O. H. E. Philcox, B. D. Sherwin, G. S. Farren, and E. J. Baxter, Phys. Rev. D 103, 023538 (2021), arXiv:2008.08084 [astro-ph.CO] .
  * Farren  _et al._ (2022) G. S. Farren, O. H. E. Philcox, and B. D. Sherwin, Phys. Rev. D 105, 063503 (2022), arXiv:2112.10749 [astro-ph.CO] .
  * Baxter and Sherwin (2021) E. J. Baxter and B. D. Sherwin, Mon. Not. Roy. Astron. Soc. 501, 1823 (2021), arXiv:2007.04007 [astro-ph.CO] .
  * Philcox  _et al._ (2022) O. H. E. Philcox, G. S. Farren, B. D. Sherwin, E. J. Baxter, and D. J. Brout, (2022), arXiv:2204.02984 [astro-ph.CO] .
  * Brout  _et al._ (2022) D. Brout _et al._ , (2022), arXiv:2202.04077 [astro-ph.CO] .
  * Wess and Zumino (1974) J. Wess and B. Zumino, Nucl. Phys. B 70, 39 (1974).
  * Beutler  _et al._ (2017) F. Beutler _et al._ (BOSS), Mon. Not. Roy. Astron. Soc. 464, 3409 (2017), arXiv:1607.03149 [astro-ph.CO] .
  * Alam  _et al._ (2015) S. Alam _et al._ (SDSS-III), Astrophys. J. Suppl. 219, 12 (2015), arXiv:1501.00963 [astro-ph.IM] .
  * Dawson  _et al._ (2013) K. S. Dawson _et al._ (BOSS), Astron. J. 145, 10 (2013), arXiv:1208.0022 [astro-ph.CO] .
  * Chudaykin  _et al._ (2020) A. Chudaykin, M. M. Ivanov, O. H. E. Philcox, and M. Simonović, Phys. Rev. D 102, 063533 (2020), arXiv:2004.10607 [astro-ph.CO] .
  * Ivanov  _et al._ (2020) M. M. Ivanov, M. Simonović, and M. Zaldarriaga, JCAP 05, 042 (2020), arXiv:1909.05277 [astro-ph.CO] .
  * van Uitert  _et al._ (2018) E. van Uitert _et al._ , Mon. Not. Roy. Astron. Soc. 476, 4662 (2018), arXiv:1706.05004 [astro-ph.CO] .
  * Hamana  _et al._ (2020) T. Hamana _et al._ , Publ. Astron. Soc. Jap. 72, Publications of the Astronomical Society of Japan, Volume 72, Issue 1, February 2020, 16, https://doi.org/10.1093/pasj/psz138 (2020), arXiv:1906.06041 [astro-ph.CO] .
  * Lesgourgues (2011) J. Lesgourgues, (2011), arXiv:1104.2932 [astro-ph.IM] .
  * Blas  _et al._ (2011) D. Blas, J. Lesgourgues, and T. Tram, JCAP 07, 034 (2011), arXiv:1104.2933 [astro-ph.CO] .
  * Mead  _et al._ (2015) A. Mead, J. Peacock, C. Heymans, S. Joudaki, and A. Heavens, Mon. Not. Roy. Astron. Soc. 454, 1958 (2015), arXiv:1505.07833 [astro-ph.CO] .
  * Mead (2017) A. Mead, Mon. Not. Roy. Astron. Soc. 464, 1282 (2017), arXiv:1606.05345 [astro-ph.CO] .
  * Mead  _et al._ (2020) A. Mead, S. Brieden, T. Tröster, and C. Heymans,  (2020), 10.1093/mnras/stab082, arXiv:2009.01858 [astro-ph.CO] .
  * Brinckmann and Lesgourgues (2019) T. Brinckmann and J. Lesgourgues, Phys. Dark Univ. 24, 100260 (2019), arXiv:1804.07261 [astro-ph.CO] .
  * Schöneberg  _et al._ (2021) N. Schöneberg, G. Franco Abellán, A. Pérez Sánchez, S. J. Witte, V. Poulin, and J. Lesgourgues, (2021), arXiv:2107.10291 [astro-ph.CO] .
  * Aghanim  _et al._ (2020b) N. Aghanim _et al._ (Planck), Astron. Astrophys. 641, A5 (2020b), arXiv:1907.12875 [astro-ph.CO] .
  * Beutler  _et al._ (2011) F. Beutler, C. Blake, M. Colless, D. H. Jones, L. Staveley-Smith, L. Campbell, Q. Parker, W. Saunders, and F. Watson, Mon. Not. Roy. Astron. Soc. 416, 3017 (2011), arXiv:1106.3366 [astro-ph.CO] .
  * Ross  _et al._ (2015) A. J. Ross, L. Samushia, C. Howlett, W. J. Percival, A. Burden, and M. Manera, Mon. Not. Roy. Astron. Soc. 449, 835 (2015), arXiv:1409.3242 [astro-ph.CO] .
  * Scolnic  _et al._ (2018) D. M. Scolnic _et al._ , Astrophys. J. 859, 101 (2018), arXiv:1710.00845 [astro-ph.CO] .
  * Camarena and Marra (2021) D. Camarena and V. Marra,  (2021), 10.1093/mnras/stab1200, arXiv:2101.08641 [astro-ph.CO] .


# PSR J2032+4127, the counterpart of TeV J2032+4130? Multiwavelength
Monitoring of the Approach to Periastron

, for the VERITAS Collaboration
University of California, Los Angeles
E-mail http://veritas.sao.arizona.edu

###### Abstract:

PSR J2032+4127 has recently been identified as being in a long period (45-50
years) binary in a highly eccentric orbit with the Be star MT91 213.
Periastron is due to occur in November 2017 and this rare occurrence has
prompted a multiwavelength monitoring campaign to determine if the system is a
$\gamma$-ray binary, and, if so, to study what would be only the second
$\gamma$-ray binary with a known compact object. In the same direction as TeV
J2032+4130, $\gamma$-ray emission from this binary system could be related to
the extended very high energy $\gamma$-ray emission from that region. As part
of this monitoring, observations are being conducted by Swift, Fermi-LAT and
VERITAS. We present the status of those observations, preliminary results and
the plan for continued monitoring through periastron.

##  1 Introduction

TeV J2032+4130 is an extended, very high energy (VHE, E $>$ 100 GeV) source
lying in the direction of the Cygnus OB2 region [1]. It was the first
unidentified VHE $\gamma$-ray source and after thirteen years of observations
it remains a mystery. The most promising potential association is with the
pulsar wind nebula of PSR J2032+4127 as has been suggested in [2].

PSR J2032+4127 has recently been identified as the compact object in a binary
system with the Be star MT91 213 [3]. This system has a highly eccentric orbit
with a period of 45-50 years [4]. Periastron is due to occur in November 2017
and a multiwavelength campaign is underway to monitor this rare event. These
observations will be used to test whether PSR J2032+4127/MT91 213 form a
$\gamma$-ray binary111Here we define $\gamma$-ray binaries as binary systems
containing a massive star and a compact object which have peak $\nu
F_{\nu}\;>\;$1 MeV as in [5] and potentially solve the mystery of the origin
of TeV J2032+4130.

PSR J2032+4127 was first identified by Fermi-LAT with a period of 142 ms [6].
Radio observations have since detected dramatic changes in its spin down rate,
at a level never seen in an isolated pulsar, from $\dot{\nu}$ =
-1$\times$10-12 s-1 to -2$\times$10-12 s-1 over the course of six years [3].
This has been attributed to Doppler shifting from the motion of the pulsar
within a binary system. The binary companion has been identified as the 15 M⊙
Be star MT91 213.

Figure 1: A schematic illustration of the orbital motion of PSR J2032+4127 and
its Be-star companion. Periastron is anticipated to take place in November
2017 [4].

Observations of the system have already shown a dramatic increase in the X-ray
flux, increasing by a factor of about 70 since 2002, and about 20 since 2010
[4]. This suggests that interactions between PSR J2032+4127 and the Be star
winds have already begun.

##  2 Comparison with Other Very High Energy $\gamma$-ray Binary Systems

$\gamma$-ray binaries are a select class of sources (only six very high energy
(VHE, E $>$ 100 GeV) $\gamma$-ray binaries have been detected [7]) and the
nature of the compact object is known in only one instance. Though the number
of known $\gamma$-ray binaries is limited, comparisons with other systems
containing Be stars shows the variety of potential outcomes of this periastron
passage. PSR B1259-63/LS 2883 is the only system where the compact object is
known (a radio pulsar) and is the most similar system to PSR J2032+4127/MT91
213. As well as VHE emission detected by H.E.S.S [8], Fermi-LAT has also
detected flaring activity about 30 days after periastron, with a peak flux
about 20-30 times pre-periastron levels, almost reaching the estimated total
spin down luminosity (fig. 2) [9]. Observations of H.E.S.S J0632+057 with
VERITAS and Swift have shown variable and correlated VHE and X-ray emission
[10, 11]. It has recently been detected in high energy (HE, 0.1 GeV $<$ E $<$
100 GeV) $\gamma$-rays using Fermi-LAT data by Malyshev and Chernyakova [12]
and is in the 3FHL catalog (3FHL J0632.7+0550) [13]. LS I +61 303 has been
shown to undergo significant orbital and superorbital variability in both the
HE and VHE bands [14]. Most notably, all three of these sources show
significant differences in their multiwavelength characteristics, in
particular the relationship between the HE and VHE $\gamma$-ray emission and
their dependencies on orbital phase. If PSR J2032+4127/MT91 213 is a
$\gamma$-ray binary it has the potential to provide significant insight into
this rare class of systems.

##  3 Swift Observations and Results

Swift is an X-ray telescope designed to detect and study $\gamma$-ray bursts.
As a result it is a flexible, rapid-response instrument capable of performing
a wider variety of observations. In this work we use the X-Ray Telescope (XRT)
which is sensitive to 0.2 – 10 keV X-rays (though we use a slightly narrower
0.3 – 10 keV band). Ho et al. [4] reported on Swift monitoring of PSR
J2032+4127 up to MJD 57633.16 (2016 September 2). Since then monitoring has
continued with observations taken every few days. Light curves have been
produced of the count rate binned by observation (fig. 2). They show that the
X-ray count rate is continuing to rise, reaching a factor of five higher than
the average level recorded in 2015 and and a factor of 35 higher than the
average level recorded in 2013 (not shown in fig. 2). There is also evidence
of variability in the light curve with peaks visible around MJD 57640 and
57860, the origin of this variability is uncertain.

Figure 2: Swift-XRT light curves of PSR J2032+4127 binned by observation [17,
18], a closer view of the data since Ho et al. [4] is shown in the inset.

##  4 Fermi-LAT Observations and Results

The Fermi Gamma-ray Space Telescope has been operating since its launch in
2008 with two $\gamma$-ray instruments: the Large Area Telescope (LAT) and the
Gamma-ray Burst Monitor (GBM). The LAT [15], Fermi’s primary instrument, is a
pair conversion $\gamma$-ray detector that is sensitive to high energy (HE)
$\gamma$-rays with energies from 20 MeV to greater than 500 GeV.

We have undertaken an analysis of Fermi-LAT Pass 8 data [16] using “SOURCE”
class photons. A joint likelihood analysis was performed with the events
separated into event types FRONT and BACK. Two different energy ranges were
tested, 300 MeV to 300 GeV and a higher energy threshold of 1 GeV to 300 GeV
to reduce the pulsar contribution. The region of interest was a square of
width 20° and the time interval covered one year, from MJD 57533 to 57897
(2016-05-24 to 2017-05-24). A base model derived from the 3FGL catalog was
used, with the _fermipy_ tool _find_sources_ used to identify new sources. Any
model components with a test statistic less than 25 were removed from the
analysis and spectral models were updated if necessary. The data were then
binned into twelve bins and a light curve was produced. Due to the variable
nature of Cygnus X-3, a known HE $\gamma$-ray emitter which lies 0.50° from
PSR J2032+4127 and was not detected in the analysis, a source was added in its
location. The data were then re-fit and another light curve produced to test
whether any variation detected in PSR J2032+4127 was a result of variation in
the Cygnus X-3 flux.

PSR J2032+4127 is detected strongly in both analyses, with test statistics of
1003 and 817, and fluxes of (6.00$\pm$0.40) $\times 10^{-9}$ cm-2 s-1 and
(2.00$\pm$0.12) $\times 10^{-9}$ for the analysis in the energy ranges 300 MeV
– 300 GeV and 1 GeV – 300 GeV respectively. The resultant light curves (fig.
3) show no marked increase, as is seen in the Swift data. They show minor
variation around the mean flux (taken from the fit with Cygnus X-3 in the
model) but, in comparison to a uniform flux model, have $\chi^{2}$ per degree
of freedom less than one. The introduction of Cygnus X-3 noticeably reduces
the flux for the 300 MeV – 1 GeV data point at MJD 57836 (2017-03-24). This is
likely due to flaring activity from Cygnus X-3 but will require further study
to confirm this.

A steady state analysis was also conducted using the whole Fermi-LAT dataset
up to MJD 57897. The measured fluxes from PSR J2032+4127 were (6.31$\pm$0.13)
$\times 10^{-9}$ cm-2 s-1 and (2.67$\pm$0.40) $\times 10^{-9}$ cm-2 s-1 over
the energy ranges 300 MeV – 300 GeV and 1 GeV – 300 GeV respectively. They
shows that the fluxes measured from MJD 57533 to 57897 have not increased
relative to the average fluxes up to MJD 57897.

(a) 300 MeV – 300 GeV

(b) 1 – 300 GeV

Figure 3: Fermi-LAT light curves covering the period MJD 57533 to 57897 for
PSR J2032+4127. Two background models were tested, one with Cygnus X-3
(orange) and (because it was not significantly detected in the analysis) one
without (blue). The average flux over the period of the light curve (black
dotted line) is also shown. In all cases the variation is not significant with
$\chi^{2}$/degrees of freedom of less than one.

##  5 VERITAS Observations and Results

The Very Energetic Radiation Imaging Telescope Array System (VERITAS) [19] is
an array of four IACTs, located at the Fred Lawrence Whipple Observatory in
southern Arizona (31°40’ N, 110°57’ W, 1.3 km a.s.l.). Each telescope is of
Davies-Cotton design, with a 12-m diameter reflector comprised of 345
hexagonal mirror facets. The focal length of each telescope is 12 m and each
telescope is equipped with a camera consisting of 499 close-packed
photomultiplier tube (PMT) “pixels” at the focus. The angular spacing between
PMTs is 0.15°, which yields a total field of view (FoV) of 3.5°. Full array
operations began in 2007 and, in the summer of 2009, the first telescope was
relocated to increase the sensitivity of the array [20]. Following a trigger
upgrade, in summer 2012, the cameras in each telescope were replaced with new,
high quantum efficiency PMTs which has resulted in a decrease of the array
energy threshold to about 85 GeV [21]. All the data presented in this work
were taken after the 2012 upgrade and using standard VERITAS operational
procedures. The results presented here were generated using one of the
standard VERITAS event reconstruction packages, similar to the scheme
described in [22] using a 0.1° integration radius (tuned for point sources).

PSR J2032+4127 lies at the edge of the extended source TeV J2032+4130
(9.5$\pm$1.2’ along the major axis and 4.0’$\pm$0.5’ minor axis). In this work
we did not study the emission from the extended source. Rather we have focused
on the emission from the direction of PSR J2032+4127, performing a point
source search at the location of PSR J2032+4127. We expect to have some level
of contamination from TeV J2032+4130 within the analysis, however, changes in
the flux will be indicative that at least part of the emission detected is
from the binary system, as it is highly unlikely that the we will detect
variability from TeV J2032+4130 due to its extended nature.

45.0 hours of data have been taken on PSR J2032+4127 since the 2012 upgrade,
giving 37.1 hours of quality-selected live time. Of this, 17.6 hours were
taken in fall 2012, 8.5 hours in fall 2016 and 11.0 hours taken in spring
2017. Further data are being taken in spring 2017 and data will also be taken
in fall 2017. Using these data, $\gamma$-ray emission is detected at the
4.82$\sigma$ level from the direction of PSR J2032+4127. The peak significance
is 5.69$\sigma$ (4.98$\sigma$ after accounting for the number of statistical
trials) at ($\alpha_{J2000}$, $\delta_{J2000}$) = ($20^{\mathrm{h}}\allowbreak
32^{\mathrm{m}}\allowbreak 17^{\mathrm{s}}$, $41\degree\allowbreak
28^{\prime}\allowbreak 30^{\prime\prime}$), 0.03°away from the location of PSR
J2032+4127 and well within the point spread function (PSF) which is 0.1° at 1
TeV. A significance sky map for the region is shown in fig. 4.

Figure 4: VERITAS significance sky map. The location of PSR J2032+4127 is
shown with black ”+”, Cygnus X-3 is shown with a cyan ”$\times$”. The
integration radius (which is equal to the PSF) is shown with the white circle
in the upper left corner.

To test for variation in the flux with time, the data was broken up into
yearly bins. The results are presented in table 1. None of the yearly datasets
reveal a more than a marginally significant detection. We also note that much
of the 2012 data were taken at a large offset from the location of PSR
J2032+4127. There is no change in the integral flux measured within the large
statistical errors on the measured values.

Table 1: The results of the VERITAS observations broken up by year of observation. All of the results are from the RR analysis. The significance (Sig) is given in unit of $\sigma$, the energy threshold (EThresh) is given in TeV and the Flux is integral flux in units of $10^{-13}$ cm-2 s-1 above the energy threshold.  Year | Live Time [hr] | Sig | Sig/$\sqrt{\mathrm{hr}}$ | EThresh | Flux
---|---|---|---|---|---
2012 | 17.62 | 2.59 | 0.78 | 0.45 | 2.9$\pm$3.8
2016 | 8.49 | 3.29 | 0.99 | 0.45 | 7.5$\pm$3.1
2017 | 11.0 | 4.17 | 1.25 | 0.45 | 7.0$\pm$2.9

Given the morphology reported in [2], we would expect to measure a flux that
is about 1/5 of the flux they report, assuming that it smoothly follows their
reported two dimensional Gaussian morphology. With a differential flux
normalization of (9.5$\pm$1.6stat$\pm$2.2sys)$\times 10^{-13}$ TeV${-1}$ cm-2
s-1 at 1 TeV and a spectral index of 2.10$\pm$0.14stat$\pm$0.21sys. This gives
a predicted integral flux above 0.45 TeV of about (4.4$\pm$1.0stat)$\times
10^{-13}$ for this analysis, below, but within statistical errors of, the flux
measured in this work. At present there is no evidence that the emission
observed is related due to the binary system, further observations are
required to investigate this in more detail.

##  6 Future Plans

With PSR J2032+4127 approaching periastron, the observation campaign is only
at the beginning. The X-ray flux shows significant variation and this suggests
that PSR J2032+4127 is beginning to interact with the stellar wind of MT91 213
and thus the possibility of detecting HE and VHE $\gamma$-rays from the binary
interaction is increasing. VERITAS will be conducting targeted observations in
the fall of 2017 to search for enhanced VHE emission.

As a survey instrument, Fermi-LAT exposure will continue, providing a
continuous monitor of PSR J2032+4127Ṫhey are also considering modifying the
observing strategy to increase the sensitivity to variability in the HE
waveband. A dedicated analysis is being conducted by members of the Fermi
collaboration which will search for enhanced emission from PSR J2032+4127 and
will include an off-pulse analysis search for changes in the continuous
emission from the region.

Swift observations are also ongoing, and we will continue to monitor the X-ray
flux through periastron. Swift also has an ultraviolet/optical telescope
(UVOT) which is co-aligned with the XRT. Data from this will also be analyzed
to search for variation in these wavelengths.

Using these three datasets (and other observations if available), a search
will be conducted for correlated emission between the wavebands in the hope
that this will identify a new $\gamma$-ray binary.

##  7 Conclusions

PSR J2032+4127 is a Fermi-LAT pulsar which is believed to be in a long period
binary system with the Be star MT91 213. Approaching periastron (expected in
November 2017) it is already showing significant variability in the X-ray
flux, with count rates reaching a factor of 35 times the level recorded in
2013 and also evidence of variability within the observed rise. A light curve
of one year of Fermi-LAT data shows no evidence for a change in flux using
monthly binning or by comparing the flux from the last year with that from the
whole mission. 37.1 hours of VERITAS data have been examined to search for
point source emission from the direction of PSR J2032+4127, with $\gamma$-ray
emission detected at the 4.82$\sigma$ level. Located within TeV J2032+4130, an
estimate was made on the expected flux from the extended source at
(4.4$\pm$1.0stat)$\times 10^{-12}$ cm-2 s-1, this is below the flux measured
in this work, but within statistical errors. There is no evidence of emission
from the binary system at this stage. Work is ongoing to analyze additional
archival data and future observations will be conducted to search for
variation in the VHE flux, which will be indicative of VHE emission from the
PSR J2032+4127/MT91 213 system.

##  8 Acknowledgements

This research is supported by grants from the U.S. Department of Energy Office
of Science, the U.S. National Science Foundation and the Smithsonian
Institution, and by NSERC in Canada. We acknowledge the excellent work of the
technical support staff at the Fred Lawrence Whipple Observatory and at the
collaborating institutions in the construction and operation of the
instrument. The VERITAS Collaboration is grateful to Trevor Weekes for his
seminal contributions and leadership in the field of VHE gamma-ray
astrophysics, which made this study possible.

This work made use of data supplied by the UK Swift Science Data Centre at the
University of Leicester.

## References

  *   * [1] Aharonian, F., et al. 2002, AAP, 393, L37
  * [2] Aliu, E., et al. 2014, ApJ, 783, 16
  * [3] Lyne, A. G., et al. 2015, MNRAS, 451, 581
  * [4] Ho, W. C. G. et al., 2017, MNRAS, 464, 1211
  * [5] Dubus, G. 2013, Astron Astrophys Rev, 21, 64
  * [6] Abdo, A. A., et al. 2009, Science, 325, 840
  * [7] http://tevcat2.uchicago.edu
  * [8] Aharonian, F., et al. 2005, A&A, 442, 1
  * [9] Abdo, A. A., et al. 2011, ApJL, 736, L11
  * [10] Acciari, V. A., et al. 2009, ApJL, 698, L94
  * [11] Aliu, E., et al. 2014, ApJ, 780, 168
  * [12] Malyshev, D. and Chernyakova, M., 2016, MNRAS, 463, 3074
  * [13] The Fermi-LAT Collaboration 2017. arXiv:1702.00664
  * [14] Aliu, E. et al. 2013, ApJ, 779, 88
  * [15] W. B. Atwood, et al. 2009, ApJ, 697, 1071
  * [16] Atwood, W. et al. 2013, in Proc. 4th Fermi Symposium
  * [17] Evans, P. A., 2007, A&A, 469, 379
  * [18] Evans, P. A., 2009, MNRAS, 397, 1177-1201
  * [19] https://veritas.sao.arizona.edu/
  * [20] A. N. Otte. 2009, arXiv:0907.4826
  * [21] N. Park. 2015, arXiv:1508.07070
  * [22] V. A. Acciari, et al. 2008, ApJ, 679, 1427