# Multiple Superconducting Phases and Unusual Enhancement of the Upper
Critical Field in UTe2

Dai Aoki1,2 E-mail<EMAIL_ADDRESS>Fuminori Honda1 Georg Knebel2 Daniel
Braithwaite2 Ai Nakamura1 DeXin Li1 Yoshiya Homma1 Yusei Shimizu1 Yoshiki J.
Sato1 Jean-Pascal Brison2 and Jacques Flouquet2 1IMR1IMR Tohoku University
Tohoku University Oarai Oarai Ibaraki Ibaraki 311-1313 311-1313 Japan
2Univ. Grenoble Alpes Japan
2Univ. Grenoble Alpes CEA CEA IRIG IRIG PHELIQS PHELIQS F-38000 Grenoble
F-38000 Grenoble France
France

###### Abstract

We performed AC calorimetry and magnetoresistance measurements under pressure
for $H\parallel a$-axis (easy-magnetization axis) in the novel heavy-fermion
superconductor UTe2. Thanks to the thermodynamic information, multiple
superconducting phases have been revealed under pressure and magnetic field.
The ($H,T$) phase diagram of superconductivity under pressure displays an
abrupt increase of the upper critical field ($H_{\rm c2}$) at low temperature
and in the high field region, and a strong convex curvature of $H_{\rm c2}$ at
high temperature. This behavior of $H_{\rm c2}$ and the multiple
superconducting phases require a state for the superconducting order parameter
more complex than a spin-triplet equal spin pairing. Above the superconducting
critical pressure, $P_{\rm c}$, we find strong indications that the possible
magnetic order is closer to antiferromagnetism than to ferromagnetism.

The recently discovered heavy-fermion superconductivity (SC) in UTe2 attracts
much attention [1, 2], because spin-triplet SC is most likely realized in this
system, lying at the proximity of ferromagnetic (FM) order. SC coexisting
microscopically with long-range FM order is already well studied [3, 4] in
three uranium compounds, UGe2 [5], URhGe [6] and UCoGe [7]. Because of the
strong internal field due to the FM moment, a spin-triplet state with equal-
spin pairing (ESP) is favored. In this case, SC can survive even under strong
internal exchange field. Furthermore SC can be even reinforced at high
magnetic field ($H$), pressure ($P$) and uniaxial stress by tuning the FM
fluctuations. When the field is applied along the intermediate hard-
magnetization axis ($b$-axis) in URhGe and UCoGe, the FM Curie temperature is
suppressed, and the FM fluctuations are remarkably enhanced. Then field-
reentrant (-reinforced) SC is observed at high fields, which highly exceeds
the so-called Pauli limit [8, 9].

A quite similar situation might be also realized in UTe2. UTe2 is a paramagnet
with a body-centered orthorhombic crystal structure (space group: $Immm$, #71,
$D_{2h}^{25}$). The large Sommerfeld coefficient, $\gamma\sim 120\,{\rm
mJ\,K^{-2}mol^{-1}}$, indicates strong correlations in the electronic states.
[10] The magnetization curves show a relatively large anisotropy between the
easy-magnetization axis ($a$-axis) and the hard-magnetization axes ($b$ and
$c$-axes). The $b$-axis is the hardest magnetization axis at low temperatures;
the magnetization curve shows a sharp 1st order metamagnetic transition at
$H_{\rm m}=35\,{\rm T}$, where the effective mass is strongly enhanced. [11,
12, 13] The value of $H_{\rm m}$ is well scaled by $T_{\chi_{\rm max}}$ ($\sim
35\,{\rm K}$), at which a broad maximum of susceptibility $\chi$ is observed
at low field.

The SC properties of UTe2 are spectacular. The SC transition occurs at $T_{\rm
sc}=1.6\,{\rm K}$ with the large specific heat jump. The residual
$\gamma$-value amounts to $\sim 40\,{\%}$ against the normal state
$\gamma$-value for the best quality sample, [14] and the entropy balance is
not satisfied, assuming a constant $\gamma$-value extrapolated from the normal
state above $T_{\rm sc}$. Strong electronic correlations are dominant in this
system, which is also indirectly confirmed by the failure of direct LDA band
structure calculations; they predict an Kondo insulating ground state [15, 2],
and large enough Coulomb repulsion ($U$) is required to turn it into a
metallic state [16, 17]. Hence, the upper critical field is large and
anisotropic, $H_{\rm c2}=6.5\,{\rm T}$ and $11\,{\rm T}$ for the $a$ and
$c$-axes, respectively. For the $b$-axis, field-reentrant SC appears at high
fields and survives up to $H_{\rm m}=35\,{\rm T}$. At $H_{\rm m}$, SC is
abruptly suppressed. Since the values of $H_{\rm c2}$ for all the field
directions highly exceed the Pauli-limit, a spin-triplet state seems
established. Furthermore, when the field direction is tilted from $b$ to
$c$-axis (hard- to hard-axis) by $\sim 30\,{\rm deg}$, SC reappears above
$H_{\rm m}\sim 40\,{\rm T}$. [18] From a microscopic point of view, the
development of FM fluctuations is detected by $\mu$SR [19] and NMR [20]
experiments. No long-range magnetic order was detected down to very low
temperatures by $\mu$SR and magnetization [21] measurements. The very small
decrease of the NMR Knight shift below $T_{\rm sc}$ supports a spin-triplet
scenario [22]. A point node gap structure was proposed by the specific heat,
thermal conductivity and penetration depth measurements [23, 24]. Furthermore,
topological SC is suggested by experiments and theory. [25, 26, 16]

Applying pressure also revealed surprising behaviors; $T_{\rm sc}$ initially
decreases slightly and then increases above $0.3\,{\rm GPa}$ up to $\sim
3\,{\rm K}$ at $1.2\,{\rm GPa}$. [27, 28] Further increasing pressure, SC is
suppressed at $P_{\rm c}\sim 1.5\,{\rm GPa}$, and a new ordered state which
probably corresponds to the magnetic order (MO) appears at $T_{\rm m}=3\,{\rm
K}$. In addition, remarkably, multiple SC phases are detected at zero field by
the AC calorimetry measurements above $0.2\,{\rm GPa}$ [27].

In order to study the multiple SC phases of UTe2 in details, we have performed
AC calorimetry measurements as a thermodynamic probe to determine ($H,T$)
phase diagrams. Magnetic field has been applied along the easy-magnetization
axis ($a$-axis) at different pressures. An unusual enhancement of $H_{\rm c2}$
at low temperature and in high field regions was found. Multiple SC phases are
found inside the SC domain under magnetic field and pressure. Our experiments
give here unique insights, which could not be observed in magnetoresistance
measurements restricted to the $H_{\rm c2}$ boundary between normal and SC
phases. [29]

Single crystal growth and experimental technique of AC calorimetry and
magnetoresistance measurements are described in Ref.Suppl.

Figure 1(a) shows the field dependence of the magnetoresistance for
$H\parallel a$-axis at $0.69\,{\rm GPa}$ at different temperatures. The
electrical current was applied along the $a$-axis. At the lowest temperature,
$0.05\,{\rm K}$, $H_{\rm c2}$ defined by the zero resistivity is enhanced up
to $9.5\,{\rm T}$ ($H_{\rm c2}\sim 6.5\,{\rm T}$ at ambient pressure). At
higher temperatures, a broad minimum appears at $H_{x}\sim 7.5\,{\rm T}$ in
the normal state. This anomaly may correspond to that detected at ambient
pressure ($\sim 6\mbox{--}10\,{\rm T}$) by magnetoresistance [12],
magnetization [11] and thermopower [31], which is most likely due to a
Lifshitz transition with Fermi surface reconstruction.

Figure 1(b) shows the ($H,T$) phase diagram at $0.69\,{\rm GPa}$ determined
from the results of magnetoresistance. The SC transition temperature is
$2.65\,{\rm K}$ at zero field and $H_{\rm c2}$ is $9.5\,{\rm T}$ at
$0.05\,{\rm K}$. The temperature dependence of $H_{\rm c2}$ plotted by solid
circles is quite unusual. Firstly, $H_{\rm c2}$ suddenly increases with a kink
at $0.6\,{\rm K}$ upon cooling. Secondly, the $H_{\rm c2}$ curve at high
temperature shows a strong convex (downward) curvature with a large initial
slope at $T_{\rm sc}$ ($-dH_{\rm c2}/dT=8.1\,{\rm T/K}$). This curvature could
be an indication of Pauli paramagnetic limitation. Note that a weak convex
curvature is also observed in specific heat measurements at ambient pressure
[24], and a similar behavior under pressure has already been observed on
magnetoresistance measurements [29]. The broad anomaly denoted by $H_{x}$ is
almost constant with slight increase with temperature and disappears above
$3.5\,{\rm K}$.

Figure 1: (Color online) (a) Magnetoresistance for $H\parallel a$-axis at
$0.69\,{\rm GPa}$ in UTe2 at different temperatures, 1.8, 1.7, 1.5, 1.2, 1.1,
0.94, 0.81, 0.56, 0.45, 0.41, 0.36, 0.32, 0.30, 0.27, 0.19 and 0.05 K. The
electrical current was applied along $a$-axis. (b) ($H,T$) phase diagram for
$H\parallel a$-axis at $0.69\,{\rm GPa}$ determined from the magnetoresistance
measurements in UTe2.

Figure 2 shows the results of AC calorimetry measurements at $0.70\,{\rm
GPa}$. In Fig. 2(a), the temperature dependence of AC calorimetry ($C_{\rm
AC}$) at zero field shows double transitions at $T_{\rm sc2}=2.8\,{\rm K}$ and
$T_{\rm sc3}=0.85\,{\rm K}$. The results are in good agreement with the
previous report [27], although $T_{\rm sc2}$ is slightly higher here. With
field, $T_{\rm sc2}$ shifts to lower temperature. On the other hand, $T_{\rm
sc3}$ slightly decreases with field and then increases again, as shown in Fig.
2(b). Here $C_{\rm AC}$ was obtained from the phase analysis, that is, $C_{\rm
AC}\propto 1/\tan(\phi)$, where $\phi$ is the phase calculated from the real
and imaginary components of the signal in the lock-in amplifier. At $5\,{\rm
T}$ the signal is rather small, but the anomaly is still visible. At higher
field, $7\,{\rm T}$, the anomaly becomes sharper again. The field scan of AC
calorimetry at different temperatures shows several anomalies, as shown in
Fig. 2(c). At $0.09\,{\rm K}$, three anomalies are detected at $9.3$, $5.9$
and $2.2\,{\rm T}$. The first anomaly corresponds to $H_{\rm c2}$ detected by
the magnetoresistance as described above. Note that the field scan of AC
calorimetry was measured at the constant frequency, thus the increase or
decrease of the signal at anomaly is arbitrary because of the field variation
of thermal conductivity, though the anomaly clearly displays the phase
boundary.

Figure 2: (Color online) (a) Temperature dependence of the AC calorimetry at
different fields at $0.70\,{\rm GPa}$ for $H\parallel a$-axis in UTe2. (b)
Temperature dependence of the AC calorimetry obtained by the phase analysis at
different fields. (c) Field dependence of the AC calorimetry at different
temperatures. The data are vertically shifted for clarity.

At $1.0\,{\rm GPa}$, three transitions are detected at $0.99$, $5.8$ and
$9.1\,{\rm T}$ in the field scan at the lowest temperature, $0.06\,{\rm K}$ as
shown in Fig. 3(a). The lower field transition is suppressed with increasing
temperature, and the higher two transitions merge at $0.62\,{\rm K}$. Note
that the transitions at $1.1$ and $1.8,{\rm K}$ lead to marked anomalies. In
the temperature scan at constant fields, $2$–$5\,{\rm T}$, no anomaly was
detected below $1\,{\rm K}$, in contrast to the case at $0.7\,{\rm GPa}$.

At $1.47\,{\rm GPa}$ which is presumably just below the SC critical pressure
$P_{\rm c}$, a single transition is detected at $6.7\,{\rm T}$ at $0.07\,{\rm
K}$ in the field scan, as shown in Fig. 3(b). The lower field transition
starts to appear above $0.8\,{\rm K}$ (at $4\,{\rm T}$) and the two
transitions merge at $2\,{\rm K}$. The temperature scan at low fields shows no
anomalies. The detected anomalies in the field scan seems a survival of SC, as
shown later in the ($H,T$) phase diagram.

Figure 3: (Color online) Field dependence of AC calorimetry at $1.0\,{\rm
GPa}$ (a) and $1.47\,{\rm GPa}$ (b) for $H\parallel a$-axis in UTe2. The data
at $2.8\,{\rm K}$ ($2.0\,{\rm K}$) at $1.0\,{\rm GPa}$ ($1.47\,{\rm GPa}$) are
scaled by factor, 0.2 (0.4). The data are vertically shifted for clarity.

At high pressures above $P_{\rm c}$ ($\sim 1.5\,{\rm GPa}$), SC is suppressed
and a new ordered phase, most likely magnetic, appears at $T_{\rm m}$. The
details are described in Ref.Suppl.

[tbh]

(Color online) (a) ($T,P$) phase diagram at zero field in UTe2. The solid
circles/squares are the present work. The open circles/squares are cited from
Ref.Bra19. (b)-(i) ($H,T$) phase diagrams for $H\parallel a$-axis at different
pressures. The lines are guides to the eyes. Figure 3(a) shows the ($T,P$)
phase diagram at zero field. The results are consistent with the previous ones
[27]. With increasing pressure, $T_{\rm sc1}$ initially decreases, then splits
above $0.2\,{\rm GPa}$. The lower transition temperature $T_{\rm sc3}$
decreases linearly with pressure and extrapolates to zero at $P_{\rm c}\sim
1.5\,{\rm GPa}$. The higher transition temperature $T_{\rm sc2}$ increases
under pressure and has a maximum near $1\,{\rm GPa}$. In the normal state,
$T_{\rm m}$ detected above $P_{\rm c}\sim 1.5\,{\rm GPa}$ is also in good
agreement with the previous results.

In Fig. 3(b)-(i), we show all the ($H,T$) phase diagrams at different
pressures. It is obvious that multiple SC phases exist at and above $0.4\,{\rm
GPa}$. Thus, UTe2 clearly belongs to the rare case of (spin-triplet)
superconductors with multiple SC states, like UPt3 [32], Th-doped UBe13 [33] ,
and of course, superfluid 3He [34].

The ($H,T$) phase diagram at $0.19\,{\rm GPa}$ is similar to that at ambient
pressure. At $0.40\,{\rm GPa}$, in zero field, two transitions $T_{\rm sc2}$
and $T_{\rm sc3}$ appear. The lower temperature $T_{\rm sc3}$ shows a
significant concave (upward) curvature with “S”-shape. On the other hand, the
higher temperature $T_{\rm sc2}$ shows a convex (downward) curvature. At
$0.54\,{\rm GPa}$, the “S”-shaped behavior of $T_{\rm sc3}$ and the convex
(downward) shape of $T_{\rm sc2}$ appear clearly. The phase boundary between
normal and superconducting state at $0.7\,{\rm GPa}$ nearly coincides with
that at $0.69\,{\rm GPa}$ detected by magnetoresistance measurements shown in
Fig. 1(b). The “S”-shape of $T_{\rm sc3}$ shows a vertical increase from 2 to
$6\,{\rm T}$. $T_{\rm sc3}$ seems to split further (at $0.4\,{\rm K}$ and
$1.3\,{\rm T}$), indicative of the evolution of the phase diagram at higher
pressure. At $1\,{\rm GPa}$, the “S”-shaped behavior is disconnected, and two
decoupled phases appear at high field and low field regions. The strong convex
(downward) curvature for $T_{\rm sc2}$ is very remarkable at pressures between
$0.4$ and $1\,{\rm GPa}$. At $1.47\,{\rm GPa}$ just below the critical
pressure $P_{\rm c}$ of SC, no anomaly is detected in the temperature scan at
zero field, while a clear double transition in the field scan associated with
SC is detected in the limited temperature range between $0.8$ and $1.4\,{\rm
K}$. At this pressure, $T_{\rm sc2}$ is rapidly suppressed as the first order
transition in the ($T,P$) phase diagram, as shown in Fig. 3(a) Therefore it is
quite difficult to detect $T_{\rm sc2}$ by a temperature scan. Note that the
magnetoresistance at $1.4\,{\rm GPa}$ for $H\parallel a$-axis reveals SC at
$1.9\,{\rm K}$ at zero field; nearly vertical initial slopes of $H_{\rm c2}$
are detected in the ($H,T$) phase diagram close to $P_{\rm c}$. [29]

At higher pressure above $P_{\rm c}$, anomalies are clearly detected again at
zero field, and the ($H,T$) phase diagrams are shown in Fig 3(h)(i). The
anomaly at $T_{\rm m}$ shifts to lower temperature with field, and it
collapses at $H_{\rm c}\sim 7$ and $10\,{\rm T}$ for $1.74$ and $1.82\,{\rm
GPa}$, respectively. For other field directions ($H\parallel b$ and $c$-axes)
[35, 36], $T_{\rm m}$ also decreases and is suppressed under magnetic field.
Thus the magnetic phase above $P_{\rm c}$ may be not a FM order. If the order
is FM, $T_{\rm m}$ should increase rapidly and change into a crossover when
the field is applied along the easy-magnetization axis. In the Ising FM case,
$T_{\rm m}$ would decrease for the field along hard-magnetization axis as
observed in URhGe and UCoGe. The decrease of $T_{\rm m}$ for all the field
directions in UTe2 is suggestive of an antiferromagnetic (AFM) order. The
phase diagram shown in Fig. 3(h)(i) is similar to the AFM one known in the
heavy fermion systems. $H_{\rm c}$ survives far above $T_{\rm m}$ up to $5$
and $7\,{\rm K}$ for $1.74$ and $1.82\,{\rm GPa}$, respectively. At
$1.82\,{\rm GPa}$, $H_{\rm c}$ slightly increases from $10$ to $11\,{\rm T}$
above $T_{\rm m}$. The shape of the anomaly at $5\,{\rm K}$ is clearly
different from those below $T_{\rm m}$, as shown in Ref. Suppl. These results
imply that the first order transition occurs at low temperatures, and that a
crossover occurs at high temperature above $T_{\rm m}$. Similar magnetic phase
diagrams are known in heavy fermion antiferromagnets, such as UPd2Al3 [37] and
CeRh2Si2 [38], where a crossover field is observed at high temperatures above
$T_{\rm N}$.

A very striking feature of our results is SC4 marked by an abrupt increase of
$H_{\rm c2}$ observed at low temperatures for pressures between $0.4$ and
$1\,{\rm GPa}$. As $H_{\rm c2}$ shows signs of Pauli limitation, an
explanation could be that it marks the appearence of a Fulde-Ferrell-Larkin-
Ovchinnikov (FFLO) state. Such a state has been recently revealed in the “two
dimensionnal” (2D) iron-based superconductors, KFe2As2 [39], and FeSe [40]. In
the FFLO state, the pair-breaking due to the Pauli paramagnetic effect is
reduced, because Cooper pairs between the two Zeeman split Fermi surfaces are
formed with a finite $\bm{q}$ for their center of mass. In UTe2, a FFLO
scenario could also be relevant, even for its spin-triplet state, precisely
because it could enhance the observed paramagnetic limit (arising from a
finite $\langle S_{Z}=0\rangle$ component of the superconducting order
parameter for that field direction). However, it is very unlikely that it
plays any role in the observed abrupt increase of $H_{\rm c2}$: first, because
UTe2 is a 3D superconductor, so the enhancement of $H_{\rm c2}$ should not be
larger than $6\%$; secondly, because the FFLO state never produces such a
sharp kink, only an inflection point (see e.g. Refs. Cho17,Kas20 ); and last,
because it is clear from our thermodynamic measurements, that the increase
emerges as a prolongation of the lower SC transition $T_{\rm sc3}$, existing
already at zero field, where no paramagnetic limitation exists.

Besides the identification of the symmetries of the different phases, and the
microscopic origin of such a complex SC phase diagram, a major question is the
observation of a strong convex (downward) curvature of $H_{\rm c2}(T)$ for
$H\parallel a$-axis. The main point, by analogy with the ferromagnetic
superconductors, is that it is difficult to imagine that the most favored
spin-triplet state has paramagnetic limitation along the field directions with
the largest susceptibility [10] and spin fluctuations [20, 19].

In order to solve this paradox, a phenomenological model has been proposed in
Refs. Kit20,Mac20, where the superconducting order parameter ($\bm{d}$-vector)
would have indeed no intrinsic paramagnetic limitation along the a-axis, but
where $T_{\rm sc}$ would be initially field-enhanced in this direction due to
the development of a strong magnetization. This may work for the weak effects
detected at zero pressure, but certainly cannot be applied to explain the very
strong curvature observed at and above $0.4\,{\rm GPa}$. From our pressure
measurement, it appears that the transition line at $T_{\rm sc3}$ is
reminiscent of $H_{\rm c2}(T)$ at $0.19\,{\rm GPa}$, where we do not observe
any paramagnetic limitations for $H\parallel a$-axis (keeping in mind that
$T_{\rm sc3}(H)$ is a transition between two superconducting phases, not an
$H_{\rm c2}$ line).

Therefore the strong paramagnetic limitation observed on $H_{\rm c2}$ ($T_{\rm
sc2}(H)$) should not be compared to the observation at zero pressure [24]. At
the opposite, it suggests that $T_{\rm sc2}$ emerges from an order parameter
different from that at zero pressure, originating from another, more
efficient, pressure-induced pairing mechanism. This would be coherent with the
observation that the pressure induced magnetic order above $P_{\rm c}$ is not
simple FM, so that magnetic fluctuations of a different type develop at the
proximity of $P_{\rm c}$. Most likely, owing to the very large $H_{\rm c2}$
observed along the $c$-axis in magnetoresistance measurements, with no sign of
paramagnetic limitation in this directions, it could be that the order
parameter is still the spin-triplet, but with a complex non-unitary
$\bm{d}$-vector having no component along the $c$-axis. The phase boundary
between SC3 and SC2 would be associated with a rotation of the
$\bm{d}$-vector, and SC4 would mark the restoration of the low pressure
dominant order parameter (at least up to $0.9\,{\rm GPa}$, the “easy-
magnetization” axis remains the $a$-axis.)

In UTe2, as in 4$f$ or 5$f$ heavy-fermion systems, pressure modifies the local
properties of U. This is well studied in Ce compounds, where the Kondo effect,
crystal field and volume drastically change with pressure. [42] As reported in
the recent paper, [29] the $P$ variation of the magnetocrystalline anisotropy
is a key player in the SC and PM properties. A fascinating new ingredient in
this system is the stabilization of a metallic state at $P=0$; pressure could
induce a critical interplay between hybridization and Coulomb repulsion,
leading to drastic changes with reinforcement of sharp electronic anomalies of
the density of states.

In recent reports with field scans along the $b$-axis, [29, 43] the link
between the pressure collapse of the metamagnetic transition at $H_{\rm m}$,
and the $H_{\rm c2}$ curves has been established. Contrary to the case of FM-
SC, the order parameter cannot be restricted to the spin-triplet state with
ESP. Thus, as shown here, field can select different order parameters with
well defined spin and orbital components, as it occurs for the clean case of
superfluid 3He with its A, A1 and B phases. A striking point is the strong
convex curvature of $H_{\rm c2}$ clearly observed here above $0.54\,{\rm
GPa}$. In our thermodynamic measurements, multiple SC phases of UTe2 have been
revealed in the ($H,T,P$) domain. Furthermore, even if the precise nature of
the MO phase above $P_{\rm c}$ remains to be solved, it appears much closer to
an AF phase than a FM one. So, this opens the possibility that FM fluctuations
collapse at $P_{\rm c}$ at the development of AF order. The relation with a
possible valence shift by pressure from U3+ to U4+ should be further explored,
as well as the pressure evolution of the band structure. The ($T,P$) phase
diagram of UTe2 remind us of the case of Tm chalcogenides, where band
structure, valence and magnetic ordered states are known to be in strong
interplay. [44, 45, 46, 47]

## Acknowledgements

We thank K. Machida, Y. Yanase, V. Mineev, S. Fujimoto, W. Knafo, K. Ishida
and Y. Tokunaga, for fruitful discussion. This work was supported by ERC
starting grant (NewHeavyFermion), KAKENHI (JP15H05884, JP15H05882, JP15K21732,
JP16H04006, JP15H05745, JP19H00646), and ICC-IMR.

## References

  * [1] S. Ran, C. Eckberg, Q.-P. Ding, Y. Furukawa, T. Metz, S. R. Saha, I.-L. Liu, M. Zic, H. Kim, J. Paglione, and N. P. Butch: Science 365, 684 (2019).
  * [2] D. Aoki, A. Nakamura, F. Honda, D. Li, Y. Homma, Y. Shimizu, Y. J. Sato, G. Knebel, J.-P. Brison, A. Pourret, D. Braithwaite, G. Lapertot, Q. Niu, M. Vališka, H. Harima, and J. Flouquet: J. Phys. Soc. Jpn. 88, 043702 (2019).
  * [3] D. Aoki and J. Flouquet: J. Phys. Soc. Jpn. 81, 011003 (2012).
  * [4] D. Aoki, K. Ishida, and J. Flouquet: J. Phys. Soc. Jpn. 88, 022001 (2019).
  * [5] S. S. Saxena, P. Agarwal, K. Ahilan, F. M. Grosche, R. K. W. Haselwimmer, M. J. Steiner, E. Pugh, I. R. Walker, S. R. Julian, P. Monthoux, G. G. Lonzarich, A. Huxley, I. Sheikin, D. Braithwaite, and J. Flouquet: Nature 406, 587 (2000).
  * [6] D. Aoki, A. Huxley, E. Ressouche, D. Braithwaite, J. Flouquet, J.-P. Brison, E. Lhotel, and C. Paulsen: Nature 413, 613 (2001).
  * [7] N. T. Huy, A. Gasparini, D. E. de Nijs, Y. Huang, J. C. P. Klaasse, T. Gortenmulder, A. de Visser, A. Hamann, T. Görlach, and H. v. Löhneysen: Phys. Rev. Lett. 99, 067006 (2007).
  * [8] F. Lévy, I. Sheikin, B. Grenier, and A. D. Huxley: Science 309, 1343 (2005).
  * [9] D. Aoki, T. D. Matsuda, V. Taufour, E. Hassinger, G. Knebel, and J. Flouquet: J. Phys. Soc. Jpn. 78, 113709 (2009).
  * [10] S. Ikeda, H. Sakai, D. Aoki, Y. Homma, E. Yamamoto, A. Nakamura, Y. Shiokawa, Y. Haga, and Y. Ōnuki: J. Phys. Soc. Jpn. Suppl 75, 116 (2006).
  * [11] A. Miyake, Y. Shimizu, Y. J. Sato, D. Li, A. Nakmura, Y. Homma, F. Honda, J. Flouquet, M. Tokunaga, and D. Aoki: J. Phys. Soc. Jpn. 88, 063706 (2019).
  * [12] W. Knafo, M. Vališka, D. Braithwaite, G. Lapertot, G. Knebel, A. Pourret, J.-P. Brison, J. Flouquet, and D. Aoki: J. Phys. Soc. Jpn. 88, 063705 (2019).
  * [13] S. Imajo, Y. Kohama, A. Miyake, C. Dong, J. Flouquet, K. Kindo, and D. Aoki: J. Phys. Soc. Jpn. 88, 083705 (2019).
  * [14] D. Aoki, A. Nakamura, F. Honda, D. Li, Y. Homma, Shimizu, Y. J. Sato, G. Knebel, J.-P. Brison, A. Pourret, D. Braithwaite, G. Lapertot, Q. Niu, M. Vališka, H. Harima, and J. Flouquet: JPS Conf. Proc. 30, 011065 (2020).
  * [15] H. Harima: JPS Conf. Proc. 29, 011006 (2020).
  * [16] J. Ishizuka, S. Sumita, A. Daido, and Y. Yanase: Phys. Rev. Lett. 123, 217001 (2019).
  * [17] Y. Xu, Y. Sheng, and Y. feng Yang: Phys. Rev. Lett. 123, 217002 (2019).
  * [18] S. Ran, I.-L. Liu, Y.S. Eo, D.J. Campbell, P. Neves, W.T. Fuhrman, S.R. Saha, C. Eckberg, H. Kim, J. Paglione, D. Graf, J. Singleton, and N.P. Butch: Nat. Phys. 15, 1250 (2019).
  * [19] S. Sundar, S. Gheidi, K. Akintola, A. M. Côté, S.R. Dunsiger, S. Ran, N. P. Butch, S. R. Saha, J. Paglione, and J. E. Sonier: Phys. Rev. B 100, 140502(R) (2019).
  * [20] Y. Tokunaga, H. Sakai, S. Kambe, T. Hattori, N. Higa, G. Nakamine, S. Kitagawa, K. Ishida, A. Nakamura, Y. Shimizu, Y. Homma, D. Li, F. Honda, and D. Aoki: J. Phys. Soc. Jpn. 88, 073701 (2019).
  * [21] C. Paulsen, G. Knebel, G. Lapertot, D. Braithwaite, A. Pourret, D. Aoki, F. Hardy, J. Flouquet, and J. P. Brison: arXiv:2002.12724.
  * [22] G. Nakamine, S. Kitagawa, K. Ishida, Y. Tokunaga, H. Sakai, S. Kambe, A. Nakamura, Y. Shimizu, Y. Homma, D. Li, F. Honda, and D. Aoki: J. Phys. Soc. Jpn. 88, 113703 (2019).
  * [23] T. Metz, S. Bae, S. Ran, I.-L. Liu, Y.S. Eo, W.T. Fuhrman, D. F. Agterberg, S. Anlage, N. P. Butch, and J. Paglione: Phys. Rev. B 100, 220504(R) (2019).
  * [24] S. Kittaka, Y. Shimizu, T. Sakakibara, A. Nakamura, D. Li, Y. Homma, F. Honda, D. Aoki, and K. Machida: arXiv:2002.06385.
  * [25] S. Bae, H. Kim, S. Ran, Y. S. Eo, I.-L. Liu, W. T. Fuhrman, J. Paglione, N. P. Butch, and S. M. Anlage: arXiv:1909.09032.
  * [26] L. Jiao, S. Howard, S. Ran, Z. Wang, J. O. Rodriguez, M. Sigrist, Z. Wang, N. Butch, and V. Madhavan: Nature 579, 523 (2020).
  * [27] D. Braithwaite, M. Vališka, G. Knebel, G. Lapertot, J. P. Brison, A. Pourret, M. E. Zhitomirsky, J. Flouquet, F. Honda, and D. Aoki: Commun. Phys. 2, 147 (2019).
  * [28] S. Ran, H. Kim, I.-L. Liu, S. Saha, I. Hayes, T. Metz, Y. S. Eo, J. Paglione, and N. P. Butch: arXiv:1909.06932.
  * [29] G. Knebel, M. Valǐska, M. Kimata, F. Honda, D. Li, D. Braithwaite, G. Lapertot, W. Knafo, A. Pourret, J.-P. Brison, J. Flouquet, and D. Aoki: arXiv:2003.08728.
  * [30] (Supplementalmaterial) Single crystal growth by the chemical vapor transport method and experimental technique under pressure are provided online. The temperature/field scans of AC calorimetry above the critical pressure are also provided online.
  * [31] Q. Niu, G. Knebel, D. Braithwaite, D. Aoki, G. Lapertot, G. Seyfarth, J.-P. Brison, J. Flouquet, and A. Pourret: Phys. Rev. Lett. 124, 086601 (2019).
  * [32] G. Bruls, D. Weber, B. Wolf, P. Thalmeier, B. Lüthi, A. de Visser, and A. Menovsky: Phys. Rev. Lett. 65, 2294 (1990).
  * [33] Y. Shimizu, S. Kittaka, S. Nakamura, T. Sakakibara, D. Aoki, Y. Homma, A. Nakamura, and K. Machida: Phys. Rev. B 96, 100505(R) (2017).
  * [34] A. J. Leggett: Rev. Mod. Phys. 47, 331 (1975).
  * [35] D. Braithwaite: private communication.
  * [36] W. Knafo: private communication.
  * [37] K. Sugiyama, T. Inoue, T. Ikeda, N. Sato, T. Komatsubara, A. Yamagishi, and M. Date: Physica B 186-188, 723 (1993).
  * [38] W. Knafo, R. Settai, D. Braithwaite, S. Kurahashi, D. Aoki, and J. Flouquet: Phys. Rev. B 95, 014411 (2017).
  * [39] C. woo Cho, J. H. Yang, N. F. Q. Yuan, J. Shen, T. Wolf, and R. Lortz: Phys. Rev. Lett. 119, 217002 (2017).
  * [40] S. Kasahara, Y. Sato, S. Licciardello, M. Čulo, S. Arsenijević, T. Ottenbros, T. Tominaga, J. Böker, I. Eremin, T. Shibauchi, N. E. H. J. Wosnitza, and Y. Matsuda: Phys. Rev. Lett. 124, 107001 (2020).
  * [41] K. Machida: J. Phys. Soc. Jpn. 89 033702 (2020).
  * [42] J. Flouquet and H. Harima: arXiv:0910.3110 (Kotai-Butsuri 47, 47 (2012) in Japanese).
  * [43] W.-C. Lin, D. J. Campbell, S. Ran, I.-L. Liu, H. Kim, A. H. Nevidomskyy, D. Graf, N. P. Butch, and J. Paglione: arXiv:2002.12885.
  * [44] V. N. Antonov, B. N. Harmon, and A. N. Yaresko: Physica B 312-313, 373 (2002).
  * [45] P. Link, I. N. Goncharenko, J. M. Mignot, T. Matsumura, and T. Suzuki: Phys. Rev. Lett. 80, 173 (1998).
  * [46] Y. Lassailly, C. Vettier, F. Holtzberg, J. Flouquet, C. M. E. Zeyen, and F. Lapierre: 28, 2880(R) (1983).
  * [47] J. Derr, G. Knebel, G. Lapertot, B. Salce, M. Méasson, and J. Flouquet: J. Phys.: Condens. Matter 18, 2089 (2006).

## Supplemental Material

High quality single crystals of UTe2 were grown using chemical vapor transport
(CVT) method with Iodine as the transport agent. [14] The magnetoresistance
under pressure was measured by the four-probe DC method using a MP35N piston
cylinder cell. The AC calorimetry measurements under pressure were performed
in another MP35N piston cylinder cell, using a Au/Fe-Au thermocouple and thin
Au wires as the thermometer and the heater, respectively. The thermocouple and
the heater were directly spot-welded on the small sample (approximately
$0.5\times 0.3\times 0.1\,{\rm mm}^{3}$ in dimension), and the signal was
detected by a lock-in amplifier with typical frequency of about $115\,{\rm
Hz}$. Daphne 7373 oil was used as pressure transmitting medium for both piston
cylinder cells. The pressure was determined by the SC transition temperature
of Pb. Both magnetoresistance and AC calorimetry measurements were done in a
top-loading dilution refrigerator at low temperatures down to $0.05\,{\rm K}$
and at high fields up to $15\,{\rm T}$.

Figure 4: (Color online) (a) Temperature dependence of the AC calorimetry at
different fields for $H\parallel a$-axis at $1.82\,{\rm GPa}$ in UTe2. (b)
Field dependence of the AC calorimetry at different temperatures. The data at
$5\,{\rm K}$ was scaled by factor 0.1 for comparison. The data are vertically
shifted for clarity.

Figure 4(a) shows the temperature dependence of the AC calorimetry at
$1.82\,{\rm GPa}$. At zero field, a jump corresponding to the magnetic order
is observed at $T_{\rm m}=3.5\,{\rm K}$. The jump becomes smaller under
magnetic fields and shifts slightly to lower temperatures. As shown in Fig.
4(b), the field scan reveals a clear anomaly at $10\,{\rm T}$ at the lowest
temperature. The anomaly survives even at high temperatures above $T_{\rm m}$
up to $7,{\rm K}$, but it becomes progressively a broad maximum, indicating a
crossover at high temperatures rather than a phase transition.


11institutetext: School of Physics and Optoelectronic, Yangtze University,
Jingzhou 434023, China; 22institutetext: Institute for Frontiers in Astronomy
and Astrophysics, Beijing Normal University, Beijing 102206, China;
33institutetext: Department of Astronomy, Beijing Normal University, Beijing
100875, China; 44institutetext: College of Physics, Hebei Normal University,
Shijiazhuang 050024, China; 55institutetext: School of Mechatronic
Engineering, Lishui Vocational $\&$ Technical College, Lishui 323000, China

# Revising the Hubble constant, spatial curvature and dark energy dynamics
with the latest observations of quasars

Tonghua Liu 1122 Shuo Cao 2233<EMAIL_ADDRESS>Xiaolei Li 44
<EMAIL_ADDRESS>Hao Zheng 55 Yuting Liu 2233 Wuzheng Guo 2233 Chenfa
Zheng 2233

(Received XXX/ Accepted XXX)

In this paper, we use a newly compiled sample of ultra-compact structure in
radio quasars and strong gravitational lensing systems with quasars acting as
background sources to constrain six spatially flat and non-flat cosmological
models ($\Lambda$CDM, PEDE and DGP). These two sets of quasar data (the time-
delay measurements of six strong lensing systems and 120 intermediate-
luminosity quasars calibrated as standard rulers) could break the degeneracy
between cosmological parameters ($H_{0}$, $\Omega_{m}$ and $\Omega_{k}$) and
therefore provide more stringent cosmological constraints for the six
cosmological models we study. A joint analysis of the quasar sample provides
model-independent estimations of the Hubble constant $H_{0}$, which is
strongly consistent with that derived from the local distance ladder by SH0ES
collaboration in the $\Lambda$CDM and PEDE model. However, in the framework of
a DGP cosmology (especially for the flat universe), the measured Hubble
constant is in good agreement with that derived from the the recent Planck
2018 results. In addition, our results show that zero spatial curvature is
supported by the current lensed and unlensed quasar observations and there is
no significant deviation from a flat universe. For most of cosmological model
we study (the flat $\Lambda$CDM, non-flat $\Lambda$CDM, flat PEDE, and non-
flat PEDE models), the derived matter density parameter is completely
consistent with $\Omega_{m}\sim 0.30$ in all the data sets, as expected by the
latest cosmological observations. Finally, according to the the statistical
criteria DIC, although the joint constraints provide substantial observational
support to the flat PEDE model, they do not rule out dark energy being a
cosmological constant and non-flat spatial hypersurfaces.

###### Key Words.:

Cosmological parameters – Gravitational lensing: strong – Quasars: general

##  1 Introduction

The $\rm{\Lambda}CDM$ model, i.e., the cosmological constant plus cold dark
matter model is so far the simplest and most natural model that could fit the
observations of type Ia supernovae (SNe Ia) (Riess et al., 2007; Alam et al.,
2017), Cosmic Microwave Background Radiation (CMB) (Planck Collaboration et
al., 2016), and strong gravitational lensing (Cao et al., 2012a; Cao & Zhu,
2014; Cao et al., 2015). However, it is still puzzled by several problems such
as the coincidence problem and the fine-tuning problem (Cao et al., 2011a).
Meanwhile, there exists some huge observational discrepancies if one tries to
motivate $\Lambda$ as a zero-point quantum vacuum energy, regarding the
estimations of different cosmological parameters in the framework of
$\rm{\Lambda}CDM$ model. One of the major issues is the inconsistency between
the Hubble constant ($H_{0}$) inferred from a $\Lambda$CDM fit to the CMB data
(temperature and polarization) from Planck satellite (Planck Collaboration et
al., 2020) and the local direct $H_{0}$ measurement by using so-called
standard candles (SNe Ia and Cepheid variables) from the Supernova $H_{0}$ for
the Equation of State collaboration (SH0ES) (Riess et al., 2019). Such tension
now has reached 4$\sigma$-6$\sigma$ with the accumulation of precise
astrophysical observations. Since there is no evidence for considerable
systematic errors in the Planck observation and the local measurements (Riess
et al., 2019; Di Valentino et al., 2018; Follin & Knox, 2018; Perivolaropoulos
& Skara, 2021), increasing attention is focusing on alternative cosmological
models beyond $\rm{\Lambda}CDM$, e.g., the Early Dark Energy models (Karwal &
Kamionkowski, 2016) and interacting dark energy models considering the
interaction between dark energy and dark matter (Caldera-Cabral et al., 2009;
Väliviita et al., 2010; Zheng et al., 2017; Li et al., 2017). Some recent
works suggested that the tension between the CMB and local determinations of
the Hubble constant could be greatly reduced, within the generalized Chaplygin
gas model (Yang et al., 2019). In the recent of Li & Shafieloo (2019), a
Phenomenologically Emergent Dark Energy (PEDE) model was introduced to shift
the constraints on $H_{0}$, which demonstrated its potential in addressing the
Hubble constant problem (Li & Shafieloo, 2019). In addition, modifications to
General Relativity (GR) theory (Jiménez Cruz & Escamilla-Rivera, 2021; Hashim
et al., 2021a, b; Briffa et al., 2020; Ren et al., 2021) or the well-known
Dvali-Gabadadze-Porrati (DGP) model physically motivated by possible
multidimensionality in the brane theory provides another way to deal with the
cosmological constant problem and alleviate the Hubble tension (Cao et al.,
2017; Xu & Wang, 2010; Xu, 2014; Giannantonio et al., 2008). Note that all of
these models could not only explain the late-time cosmic acceleration from
different mechanisms, but also describe the large-scale structure distribution
of the Universe (see Koyama (2016) for recent reviews). In this paper, we will
explore the validity of three spatially flat and non-flat cosmological models
($\Lambda$CDM, PEDE and DGP), focusing on the time delay measurements from
strongly lensed quasars and the angular size measurements of ultra-compact
structure in radio quasars. Specially, the Hubble constant, spatial curvature
and dark energy dynamics will be revisited with with such a newly compiled
quasar sample.

It is well known that the time-delays from strong gravitational lensing
systems provide an independent method to measure the Hubble constant (time-
delays are inversely proportional to the $H_{0}$). For a specific strong
gravitational lensing system, a distant active galactic nucleus (AGN), which
usually acts as the background source, is gravitationally lensed into multiple
images by a foreground early-type galaxy. Meanwhile, the light emitted from
background sources at the same time will arrive at the Earth at different
time. Due to the variable nature of quasars, the precise measurements of time-
delays between multiple images are realizable by monitoring flux variations of
the lens. Actually, the time-delays are directly related to the so-called
”time-delay distance”, a combination of three angular diameter distances of
the lensed quasar systems (from observer to lens, from observer to source, and
from lens to source). Such idea was recently realized by the $H_{0}$ Lenses in
COSMOGRAIL’s Wellspring (H0LiCOW) collaboration (Wong et al., 2020), which
presented the fits on the Hubble constant and other cosmological parameters
using a joint analysis of six gravitationally lensed quasars. In the framework
of six different cosmological models, their results showed that the derived
Hubble constant is in agreement with local distance ladder measurement in the
spatially flat $\Lambda$CDM model. However, the measured time-delays from
lensed quasars, which is only primarily sensitive to $H_{0}$, whereas
demonstrate their relatively weak constraints on other cosmological
parameters. For instance, the determined value of the matter density parameter
in the flat $\Lambda$CDM model, $\Omega_{m}=0.30^{+0.13}_{-0.13}$, would shift
to $\Omega_{m}=0.24^{+0.16}_{-0.13}$ in the non-flat $\Lambda$CDM model.
Moreover, the drawback of this method is that the fits on the Hubble constant
are strongly model dependent, i.e. the value of $H_{0}$ would shift to
$H_{0}=81.6^{+4.9}_{?5.3}\mathrm{\leavevmode\nobreak\ km\leavevmode\nobreak\
s^{-1}\leavevmode\nobreak\ Mpc^{-1}}$ when the dynamics of dark energy is
taken into consideration (Wong et al., 2020). For the discussions about model-
independent measurements of $H_{0}$, we refer the reader to see the following
works (Liao et al., 2019a, 2020; Lyu et al., 2020; Collett et al., 2019; Wei &
Melia, 2020; Qi et al., 2021; Taubenberger et al., 2019).

On the other hand, the angular-size/redshift relation of the ultra-compact
structures in unlensed radio quasars was also proposed for cosmological
applications (Kellermann, 1993). In the subsequent analysis, based on the
milliarcsecond angular size measurements from very-long-baseline
interferometry (VLBI) technique, Cao et al. (2017) demonstrated the
possibility of using intermediate-luminosity quasars as standard rulers for
cosmological inference, covering the redshift range of $0.462<z<2.73$. In the
framework of such a reliable cosmological standard ruler extended to higher
redshifts, great efforts have been made in the recent studies to set
observational limits on different cosmologies (Li et al., 2017; Zheng et al.,
2017), which shows that radio quasars could provide quite stringent
constraints on cosmological parameters. However, one issue which should be
discussed is the strong degeneracy between the Hubble constant $H_{0}$, the
matter density parameter $\Omega_{m}$ (Cao et al., 2017), the cosmic curvature
$\Omega_{k}$ (Qi et al., 2019b), and the equation of state of dark energy
$\omega$ (Di Valentino et al., 2020, 2021; Handley, 2021). Therefore, one may
expect that the combination of the latest observations of quasars, i.e., the
angular size of compact structure in radio quasars as standard rulers and the
time delays from gravitationally lensed quasars, would break the degeneracy
between the Hubble constant and other cosmological parameters, in the
framework of different cosmological models of interest. This paper is
organized as follows: In Sec. 2, we summarize the cosmological models to be
analyzed. In Sec. 3 we briefly describe the quasar data and the corresponding
analysis method. In Sec. 4 we report the results of constraints on the Hubble
constant, spatial curvature and dark energy dynamics with the latest quasar
data. Finally, we give our discussion and conclusions in Sec. 5.

##  2 Cosmological models

Let us now describe the models we are going to analyse in the next section
with the data sets. In this paper, we concentrate on three classes of
cosmological models in a spatially non-flat and flat universe, including the
standard $\Lambda$CDM model, the Phenomenologically Emergent Dark Energy
(PEDE) model, and the Dvali-Gabadadze-Porrati (DGP) model.

Assuming the Friedmann-Lematre-Robertson-Walker (FLRW) metric, with the non-
flat universe filled with ordinary pressureless matter (cold dark matter plus
baryons), dark energy, and negligible radiation, the Friedmann equation reads
as

$\displaystyle
H^{2}(z)=H_{0}^{2}[\Omega_{m}(1+z)^{3}+{\widetilde{\Omega}_{DE}}(z)+\Omega_{k}(1+z)^{2}],$
(1)

where $\Omega_{m}$ and $\Omega_{k}$ are the present values of the density
parameters of dust matter and spatial curvature, respectively. The dark energy
component of $\widetilde{\Omega}_{DE}(z)$ takes the following form of

$\widetilde{\Omega}_{DE}(z)=\Omega_{DE}\exp\big{\\{}3\int^{z}_{0}\frac{1+\omega(z^{\prime})}{1+z^{\prime}}dz^{\prime}\big{\\}},$
(2)

where $\Omega_{DE}$ is the present value of the dark energy density parameter,
and the equation of state of dark energy is defined as
$\omega(z)=p_{DE}/\rho_{DE}$, with $p_{DE}$ and $\rho_{DE}$ the pressure and
energy density of dark energy, respectively. When $\omega(z)=-1$ this so-
called XCDM parameterization reduces to the concordance $\Lambda$CDM model,
with $\widetilde{\Omega}_{DE}(z)=1-\Omega_{m}-\Omega_{k}$. Recently, another
kind of Phenomenologically Emergent Dark Energy (PEDE) model proposed in Li &
Shafieloo (2019, 2020) has attracted a lot of attention. In this cosmological
scenario, the density of dark energy, which has no effective presence in the
past and emerges in the later times, is written as

$\widetilde{\Omega}_{DE}(z)=\Omega_{DE}\times[1-\tanh(\log_{10}(1+z))].$ (3)

Note that compared with $\Lambda$CDM cosmology, the PEDE model has no extra
degree of freedom. For the third scenario, our idea of modifying the gravity
is based on the assumption that our universe is embedded in a higher
dimensional space-time, arising from the braneworld theory (Dvali et al.,
2000). In the DGP model, the cosmic acceleration is reproduced by the leak of
gravity into the bulk at large scales, which result in the accelerated
expansion of the Universe without the need of dark energy (Dvali et al.,
2001). In the framework of a non-flat DGP model, the Friedman equation is
modified as (Deffayet et al., 2002b, a; Multamäki et al., 2003)

$\displaystyle
H^{2}(z)=H_{0}^{2}[(\sqrt{\Omega_{m}(1+z)^{3}+\Omega_{r_{c}}}+\sqrt{\Omega_{r_{c}}})^{2}+\Omega_{k}(1+z)^{2}],$
(4)

where the density parameter $\Omega_{r_{c}}=1/(4r_{c}^{2}H_{0}^{2})$ is
associated with the length of $r_{c}$ where the leaking occurs. Based on the
normalization condition, $\Omega_{r_{c}}$ is also related to $\Omega_{m}$ and
$\Omega_{k}$ as
$(\sqrt{\Omega_{m}+\Omega_{r_{c}}}+\sqrt{\Omega_{r_{c}}})^{2}+\Omega_{k}=1$.

Summarizing, the Friedmann equations of the all the models presented in the
this section will be used to calculate the angular diameter distance

$\displaystyle
D_{A}(z;\mathbf{p})=\frac{c}{1+z}\int^{z}_{0}\frac{1}{H(z^{\prime})}dz^{\prime},$
(5)

where $\mathbf{p}$ denotes relevant cosmological model parameters, i.e.,
$\mathbf{p}\,=\,[H_{0},\Omega_{m}]$ for flat cosmological models and
$\mathbf{p}\,=\,[H_{0},\Omega_{m},\Omega_{k}]$ for non-flat cosmological
models.

##  3 Observational quasar data and Methodology

Here we work with large-scale data, selecting and combining complete samples
of quasars to investigate the late-time Universe. In this section, we describe
each data set and the methodology used for the cosmological analyses. These
are carried out using the probes separately, followed by the joint analysis.

Distance measurements from lensed quasars.— In strong lensing systems with
quasars acting as background sources, the time difference (time delay) between
two images of the source depends on the time-delay distance $D_{\mathrm{\Delta
t}}$ and the gravitational potential of the lensing galaxy by (Perlick, 1990a,
b; Treu & Marshall, 2016)

$\Delta t_{i,j}=\frac{D_{\mathrm{\Delta
t}}}{c}\Delta\phi_{i,j}(\mathbf{\xi}_{lens}).$ (6)

Here $c$ is the speed of light, and
$\Delta\phi_{i,j}=[(\mathbf{\theta}_{i}-\mathbf{\beta})^{2}/2-\psi(\mathbf{\theta}_{i})-(\mathbf{\theta}_{j}-\mathbf{\beta})^{2}/2+\psi(\mathbf{\theta}_{j})]$
represents the Fermat potential difference between the image $i$ and image
$j$, which is determined by the lens model parameters $\mathbf{\xi}_{lens}$
inferred from high resolution imaging of the host arcs. $\mathbf{\theta}_{i}$
and $\mathbf{\theta}_{j}$ are the angular positions of the image $i$ and $j$
in the lens plane. It is worth noting that the line-of-sight (LOS) mass
distribution to the lens could also affect the Fermat potential inference, the
contribution of which requires deep and wide field imaging of the area around
the lens system. The two-dimensional lensing potential at the image positions,
$\psi(\mathbf{\theta}_{i})$ and $\psi(\mathbf{\theta}_{j})$, and the unlensed
source position $\mathbf{\beta}$ can be determined by the lens mass model. The
time-delay distance $D_{\mathrm{\Delta t}}$ is a combination of three angular
diameter distances expressed as

$D_{\mathrm{\Delta
t}}\equiv(1+z_{d})\frac{D^{A}_{\mathrm{d}}(z;\mathbf{p})D^{A}_{\mathrm{s}}(z;\mathbf{p})}{D^{A}_{\mathrm{ds}}(z;\mathbf{p})},$
(7)

where the superscript ($A$) denotes the angular diameter distance, while the
subscripts ($d$ and $s$) represent the deflector (or lens) and the source,
respectively.

Moreover, the angular diameter distance to the deflector (or lens) can be
independently inferred from the kinematic modeling with additional information
of the lensing galaxy. The measured velocity dispersion provides the depth of
gravitational potential at the lensing position, while the time delay provides
the mass of the lensing galaxy enclosed within the position at which images
are formed. Therefore, the combination of the above two quantities will
generate the physical size of the system, on the base of which one could
obtain the measurement of $D_{d}$ at the lens position divided by the angular
separation of lensed images. More specifically, by choosing a suitable mass
density profile (such as the power-law lens distribution) and combine it with
the kinematic information of the lensing galaxy (the light distribution
function $\mathbf{\xi}_{light}$, the projected stellar velocity dispersion
$\sigma_{P}$, and the anisotropy distribution of the stellar orbits
$\beta_{ani}$), one can obtain the angular diameter distance to the lens

$D_{d}=\frac{1}{1+z_{d}}D_{\Delta
t}\frac{c^{2}H(\mathbf{\xi}_{lens},\mathbf{\xi}_{light},\beta_{ani})}{\sigma^{2}_{P}},$
(8)

where the function $H$ captures all of the model components calculated from
the sky angle (from the imaging data) and the anisotropy distribution of the
stellar orbit (from the spectroscopy). Here, we summarise the crucial points
required by the present work and one could refer to Birrer et al. (2016, 2019)
for more details. Note that the cosmological constraints obtained from the
$D_{d}$ sample are generally weaker than those from the $D_{\Delta t}$ sample.
However, the previous analysis has also demonstrated its potential in breaking
the possible degeneracies among cosmological parameters, particularly those
between cosmic curvature and the redshift-varying equation of state of dark
energy in some non-flat dark energy models (Jee et al., 2016).

The latest sample of strong-lensing systems with time-delay observations,
recently released by the H0LiCOW collaboration consist of six lensed quasars
covering the redshift range of $0.654<z_{s}<1.789$ (Wong et al., 2020):
B1608+656 (Suyu et al., 2010; Jee et al., 2019), RXJ1131-1231 (Suyu et al.,
2013, 2014; Chen et al., 2019), HE0435-1223 (Wong et al., 2017; Chen et al.,
2019), SDSS1206+4332 (Birrer et al., 2019), WFI2033-4723 (Rusu et al., 2020),
and PG1115+080 (Chen et al., 2019). The redshifts of both lens and source, the
time-delay distances, and the angular diameter distance to the lenses for
these lensed quasar systems are summarized in Table 2 of Wong et al. (2020).
Note that for the lens of B1608+656, its $D_{\Delta t}$ measurement is given
in the form of skewed log-normal distribution (due to the absence of blind
analysis of relevant cosmological quantities), while the derived $D_{\Delta
t}$ for other five lenses are given in the form of Monte Carlo Markov chains
(MCMC). As for the measurements of the angular diameter distance to the the
lens $D_{d}$, only four strong lensing systems (B1608+656, RXJ1131-1231,
SDSS1206+4332, and PG1115+080) are used in our statistical analysis, which are
provided in the form of MCMC. We remark here that a kernel density estimator
is used to compute the posterior distributions of $\mathcal{L}_{(D_{\Delta
t},D_{d})}$ or $\mathcal{L}_{D_{d}}$ from chains, which allows to account for
any correlations between $D_{\Delta t}$ and $D_{d}$ in
$\mathcal{L}_{(D_{\Delta t},D_{d})}$. The distances posterior distributions
for six time delay distances (denoted as ”$6D_{\Delta t}$” for simplicity) and
four angular diameter distances to the lenses (denoted as ”$4D_{d}$” for
simplicity) are available on the website of H0LiCOW 111http://www.h0licow.org.
For more works in cosmology by using strong lensing time delays, we refer the
reader to see the literature (Ding et al., 2021; Bag et al., 2022; Sonnenfeld,
2021; Liao, 2021; Rathna Kumar et al., 2015).

Distance measurements from radio quasars.— Being the brightest sources in the
universe, quasars have considerable potential to be used as useful
cosmological probes (Liu et al., 2020a, 2021a), despite of the extreme
variability in their luminosity and high observed dispersion. For instance,
Risaliti & Lusso (2019) attempted to use quasars as standard candles through
the non-linear relation between their intrinsic UV emission from an accretion
disk and the X-ray emission from hot corona, through the analysis and
refinement of the quasar sample with well-measured X-ray and UV fluxes. In
this work, we focus on the ”angular size-distance” relation of ultra-compact
structure in radio quasars that can be observed up to high redshifts, with
milliarcsecond angular sizes measured by the very-long-baseline interferometry
(VLBI). Specially, with the the signal received at multiple radio telescopes
across Earth’s surface, together with the registered correlated intensities
considering the different arrival times at various facilities, the
characteristic angular size of a distant radio quasar is defined as

$\theta={2\sqrt{-\ln\Gamma\ln 2}\over\pi B}$ (9)

where $B$ is the interferometer baseline measured in wavelengths and the
visibility modulus $\Gamma=S_{c}/S_{t}$ is the ratio between the total and
correlated flux densities. With gradually refined selection technique and the
elimination of possible systematic errors, Cao et al. (2018) compiled a sample
of 120 intermediate-luminosity radio quasars ($10^{27}W/Hz<L<10^{28}W/Hz$)
with reliable measurements of the angular size of the compact structure from
updated VLBI observations. It is now understood that the dispersion in linear
size is greatly mitigated (Cao et al., 2019), i.e., the linear sizes of these
standard rulers show negligible dependence on both redshifts and intrinsic
luminosity. Our quasar data come from a newly compiled sample of these
standard rulers from observations of 120 intermediate-luminosity quasars with
angular sizes $\theta(z)$ and redshifts $z$ listed in Table 1 of Cao et al.
(2017), which extend the Hubble diagram to a redshift range $0.46<z<2.76$
currently inaccessible to the traditional resorts.

The corresponding theoretical predictions for the angular sizes of the compact
structure can be written as

$\theta(z)=\frac{l_{m}}{D^{A}(z;\mathbf{p})},$ (10)

where $D^{A}(z)$ is the angular diameter distance at redshift $z$, which is
related to different combinations of cosmological parameters $\mathbf{p}$
(Hubble constant and the dimensionless expansion rate expansion rate). The
intrinsic length $l_{m}$ needs to be calibrated with external indicators such
SNe Ia. In this analysis, we adopt the calibration results of such quantity
$l_{m}=11.03\pm 0.25$ pc through a new cosmology-independent technique, by
using the Gaussian Process to reconstruct the expansion history of the
Universe from 24 cosmic chronometer measurements (Cao et al., 2018). The data
of ultra-compact structures in radio quasars, which is denoted as ”QSO” in
this work, has been extensively used for cosmological applications in the
literature (Li et al., 2017; Cao et al., 2018; Qi et al., 2017; Ma et al.,
2017; Cao et al., 2019).

Now the posterior likelihood $\mathcal{L}_{QSO}$ for can be constructed by

$\displaystyle\mathcal{L}_{QSO}=\prod_{i=1}^{120}\frac{1}{\sqrt{2\pi(\sigma_{i}^{2}+\sigma_{sys}^{2})}}\times\exp\bigg{\\{}-\frac{\left[\theta(z_{i};\mathbf{p})-\theta_{obs,i}\right]^{2}}{2(\sigma_{sta}^{2}+\sigma_{sys}^{2})}\bigg{\\}},$
(11)

where $\theta_{obs,i}$ is the observed angular size for the ith data point in
the sample, $\sigma_{sta}$ is the observational statistical uncertainty for
the ith quasar. According to the error strategy proposed in Cao et al. (2017),
an additional $10\%$ systematic uncertainty ($\sigma_{sys}=0.1\theta_{obs}$)
in the observed angular sizes is also considered, accounting for the intrinsic
spread in the linear size (Cao et al., 2019; Qi et al., 2021). Such strategy
has been extensively applied in the subsequent cosmological studies with such
standard ruler data, which extended our understanding of the evolution of the
Universe to $z\sim 3$ (Ryan, 2021; Vavryčuk & Kroupa, 2020; Melia, 2018).

Summarizing, we use different combinations of lensed and unlensed quasars by
adding the log-likelihood of cosmological parameters $H_{0}$, $\Omega_{m}$ and
$\Omega_{k}$ (if available) for the posterior distributions of six lens time-
delay distances $\mathcal{L}_{D_{\Delta t}}$, four angular diameter distances
to the lenses $\mathcal{L}_{D_{d}}$, and 120 angular diameter distances to the
radio quasars $\mathcal{L}_{QSO}$. The final log-likelihood

$\mathrm{ln}{\cal L}=\mathrm{ln}({\cal L}_{\rm{QSO}})+\mathrm{ln}({\cal
L}_{(D_{\Delta t},D_{d})}),$ (12)

is sampled in the framework of Python MCMC module EMCEE (Foreman-Mackey et
al., 2013). Meanwhile, in this analysis we introduce the DIC to evaluate which
model is more consistent with the observational data, focusing on the Deviance
Information Criterion (DIC) to compare the goodness of fit on models with
different numbers of parameters (Spiegelhalter et al., 2002). The DIC is
defined as

$DIC=D(\bar{\theta})+2p_{D}=\overline{D(\theta)}+p_{D},$ (13)

where $D(\theta)=-2\mathrm{ln}{\cal L(\theta)}+C$, $C$ is a normalized
constant depending only on the data that disappears from ant derived quantity,
and $p_{D}=\overline{D(\theta)}-D(\bar{\theta})$ is the effective number of
model parameters, with the deviance of the likelihood $D$. Specially, when we
use $\chi^{2}=-2\mathrm{ln}{\cal L(\theta)}$ to describe the $p_{D}$, it can
be rewritten as

$p_{D}=\overline{\chi^{2}(\theta)}-\chi^{2}(\bar{\theta}).$ (14)

Compared with the widely-used Akaike Information Criterion (AIC) or Bayesian
Information Criterion (BIC), the advantage of DIC lies in the fact that it is
determined by the qualities which can be easily obtained from Monte Carlo
posterior samples. Moreover, parameters that are unconstrained by the data
would also be appropriately treated in the framework of DIC (Liddle, 2007).

Table 1: The best-fit values and 1$\sigma$ uncertainties for the parameters ($H_{0}$, $\Omega_{m}$, $\Omega_{r_{c}}$ (if available), and $\Omega_{k}$) in each cosmological model and quasar data set. The $\chi^{2}$ and DIC values for all models are also added for comparison. Model | Data set | $H_{0}(\mathrm{\leavevmode\nobreak\ km\leavevmode\nobreak\ s^{-1}\leavevmode\nobreak\ Mpc^{-1}})$ | $\Omega_{m}$ | $\Omega_{k}$ | $\Omega_{r_{c}}$ | $\chi^{2}$ | DIC
---|---|---|---|---|---|---|---
Flat $\Lambda$CDM | 6$D_{\Delta t}$ | $73.26^{+1.74}_{-1.84}$ | $0.297^{+0.133}_{-0.127}$ | $-$ | $-$ | 39.98 | 42.64
| 6$D_{\Delta t}$+QSO | $73.51^{+1.73}_{-1.82}$ | $0.281^{+0.050}_{-0.041}$ | $-$ | $-$ | 358.99 | 364.21
| 6$D_{\Delta t}$+$4D_{d}$+QSO | $73.42^{+1.74}_{-1.84}$ | $0.285^{+0.049}_{-0.041}$ | $-$ | $-$ | 385.38 | $391.33$
Flat PEDE | 6$D_{\Delta t}$ | $74.94^{+1.89}_{-2.03}$ | $0.259^{+0.173}_{-0.137}$ | $-$ | $-$ | 39.57 | 42.76
| 6$D_{\Delta t}$+QSO | $75.13^{+1.87}_{-2.07}$ | $0.294^{+0.046}_{-0.040}$ | $-$ | $-$ | 358.34 | 363.81
| 6$D_{\Delta t}$+$4D_{d}$+QSO | $75.08^{+1.83}_{-2.05}$ | $0.297^{+0.046}_{-0.038}$ | $-$ | $-$ | 384.53 | 390.68
Flat DGP | 6$D_{\Delta t}$ | $67.75^{+2.43}_{-3.95}$ | $0.251^{+0.184}_{-0.137}$ | $-$ | $0.140^{+0.056}_{-0.059}$ | 40.43 | 45.74
| 6$D_{\Delta t}$+QSO | $67.84^{+1.12}_{-1.19}$ | $0.243^{+0.043}_{-0.037}$ | $-$ | $0.143^{+0.015}_{-0.016}$ | 360.02 | 367.72
| 6$D_{\Delta t}$+$4D_{d}$+QSO | $67.84^{+1.10}_{-1.12}$ | $0.246^{+0.043}_{-0.037}$ | $-$ | $0.142^{+0.014}_{-0.016}$ | 386.42 | 393.67
Non-flat $\Lambda$CDM | 6$D_{\Delta t}$ | $74.38^{+2.12}_{-2.34}$ | $0.242^{+0.164}_{-0.129}$ | $0.258^{+0.167}_{-0.253}$ | $-$ | 39.80 | 44.67
| 6$D_{\Delta t}$+QSO | $73.39^{+2.06}_{-2.11}$ | $0.274^{+0.083}_{-0.080}$ | $0.036^{+0.223}_{-0.214}$ | $-$ | 358.31 | 366.89
| 6$D_{\Delta t}$+$4D_{d}$+QSO | $73.78^{+1.99}_{-2.17}$ | $0.254^{+0.083}_{-0.074}$ | $0.100^{+0.214}_{-0.213}$ | $-$ | 384.69 | 393.81
Non-flat PEDE | 6$D_{\Delta t}$ | $74.84^{+15.81}_{-8.72}$ | $0.247^{+0.180}_{-0.131}$ | $-0.027^{+0.340}_{-0.306}$ | $-$ | 39.49 | 44.23
| 6$D_{\Delta t}$+QSO | $71.67^{+16.75}_{-7.38}$ | $0.315^{+0.094}_{-0.106}$ | $-0.089^{+0.365}_{-0.276}$ | $-$ | 357.82 | 366.39
| 6$D_{\Delta t}$+$4D_{d}$+QSO | $75.05^{+2.01}_{-2.19}$ | $0.275^{+0.081}_{-0.077}$ | $0.076^{+0.222}_{-0.201}$ | $-$ | 383.81 | 393.45
Non-flat DGP | 6$D_{\Delta t}$ | $68.31^{+15.37}_{-10.45}$ | $0.251^{+0.186}_{-0.142}$ | $0.023^{+0.316}_{-0.345}$ | $0.135^{+0.143}_{-0.116}$ | 40.31 | 46.24
| 6$D_{\Delta t}$+QSO | $66.26^{+13.57}_{-8.04}$ | $0.253^{+0.089}_{-0.081}$ | $-0.050^{+0.328}_{-0.305}$ | $0.151^{+0.107}_{-0.101}$ | 359.47 | 367.72
| 6$D_{\Delta t}$+$4D_{d}$+QSO | $66.35^{+14.92}_{-8.10}$ | $0.253^{+0.086}_{-0.083}$ | $-0.047^{+0.348}_{-0.311}$ | $0.150^{+0.110}_{-0.104}$ | 385.19 | 395.27

##  4 Results and discussion

Figure 1: The 2-D plots and 1-D marginalized distributions with 1-$\sigma$ and
2-$\sigma$ contours of cosmological parameters ($H_{0}$ and $\Omega_{m}$) in
the framework of flat $\Lambda$CDM (upper), flat PEDE (lower left) and flat
DGP (lower right) models with lensed quasars ($6D_{\Delta t}$, $4D_{d}$) and
radio quasars (QSO).

Figure 2: The 2-D plots and 1-D marginalized distributions with 1-$\sigma$ and
2-$\sigma$ contours of cosmological parameters ($H_{0}$, $\Omega_{m}$ and
$\Omega_{k}$) in the framework of non-flat $\Lambda$CDM (upper), non-flat PEDE
(lower left) and non-flat DGP (lower right) models with lensed quasars
($6D_{\Delta t}$, $4D_{d}$) and unlensed quasars (QSO).

In order to demonstrate the constraining power of the latest observations of
quasars, we use different data combinations (6$D_{\Delta t}$, 6$D_{\Delta
t}$+QSO, and 6$D_{\Delta t}$+4$D_{d}$+QSO) to place constraints on
cosmological parameters in $\Lambda$CDM, PEDE and DGP models both in flat and
non-flat cases. The posterior one-dimensional (1-D) probability distributions
and two-dimensional (2-D) confidence regions of the cosmological parameters
for the six flat and non-flat models are shown in Figs. 1-2. We list the
marginalized best-fitting parameters and $1\sigma$ uncertainties for all
models and data combinations in Table 1. The corresponding $\chi^{2}$ and DIC
values are also listed in Table 1.

Flat cases.— The constraint results for three flat cosmological models are
presented in Fig. 1, in which we show the 2-D confidence regions (with
1$\sigma$ and 2$\sigma$ limits) as well as 1-D marginalized distributions from
different data combinations. Our findings show that the constraints we obtain
from the combination of the latest observations of quasars are more reliable
than the those derived from independent quasar sample. Fortunately, the
results in Fig. 2 show that the radio quasar data place remarkable constraints
on all parameters of the three cosmological models, and the degeneracy among
different model parameters are well broken. When combined with 120 unlensed
quasar data, the 6 time-delay lensed quasars would produce tighter constraints
on the matter density parameter in all of the three cosmological scenarios.
The measured $\Omega_{m}$ range from a low value of $0.243^{+0.043}_{-0.037}$
(flat DGP) to a high value of $0.294^{+0.046}_{-0.040}$ (flat PEDE). In
particular, for flat $\Lambda$CDM the derived matter density parameter from
120 unlensed quasars and 6 time-delay lensed quasars
$\Omega_{m}=0.297^{+0.133}_{-0.127}$ shows a perfect agreement with the TT,
TE, EE+lowE+lensing results of Planck Collaboration ($\Omega_{m}=0.3103\pm
0.0057$) (Planck Collaboration et al., 2020). Such findings are different from
the those obtained from the latest compilation of X-ray+UV quasars acting as
standard candles, which favors a larger value of the matter density parameter
at higher redshifts (Risaliti & Lusso, 2019). See also Lian et al. (2021) for
more discussion about this issue. In spite of the low mean values of
$\Omega_{m}$ in the flat DGP model, the constraints on $\Omega_{m}$ obtained
with 6$D_{\Delta t}$, 6$D_{\Delta t}$+QSO, and 6$D_{\Delta t}$+QSO+4$D_{d}$
data are mutually consistent with those obtained from other astrophysical
probes (Xu & Wang, 2010; Xu, 2014; Giannantonio et al., 2008). For comparison,
the corresponding fits on the parameter of $\Omega_{r_{c}}$ in flat DGP model
are also displayed in Table 1. Finally, our analysis demonstrates that the
matter density parameter plays an important role in the determination of the
Hubble constant, which can be clearly seen from the anti-correlation between
$\Omega_{m}$ and $H_{0}$ in Fig. 1.

The constraints on the Hubble constant are between
$H_{0}=67.75^{+2.43}_{-3.95}$ $\mathrm{\leavevmode\nobreak\
km\leavevmode\nobreak\ s^{-1}\leavevmode\nobreak\ Mpc^{-1}}$ (DGP) and
$H_{0}=74.94^{+1.89}_{-2.03}$ $\mathrm{\leavevmode\nobreak\
km\leavevmode\nobreak\ s^{-1}\leavevmode\nobreak\ Mpc^{-1}}$ (PEDE) with 6
time-delay lensed quasars, which shift to $H_{0}=67.84^{+1.12}_{-1.19}$
$\mathrm{\leavevmode\nobreak\ km\leavevmode\nobreak\
s^{-1}\leavevmode\nobreak\ Mpc^{-1}}$ (DGP) and $H_{0}=75.13^{+1.87}_{-2.07}$
$\mathrm{\leavevmode\nobreak\ km\leavevmode\nobreak\
s^{-1}\leavevmode\nobreak\ Mpc^{-1}}$ (PEDE) with the combined 6$D_{\Delta
t}$+QSO data. Specially, for the flat $\Lambda$CDM and PEDE model, the mean
values of $H_{0}$ obtained with 6$D_{\Delta t}$+4$D_{d}$+QSO data are more
consistent with the recent determinations of $H_{0}$ from the Supernovae H0
for the SH0ES collaboration (Riess et al., 2019). However, in the framework of
flat DGP model, the measured value of $H_{0}$ with 1$\sigma$ uncertainty
($67.84^{+1.12}_{-1.19}$ $\mathrm{\leavevmode\nobreak\ km\leavevmode\nobreak\
s^{-1}\leavevmode\nobreak\ Mpc^{-1}}$), which is 3.6$\sigma$ lower than the
statistical estimates of the SH0ES results, demonstrates a perfect agreement
with that derived from the recent Plank CMB observations (Planck Collaboration
et al., 2020). In addition, relative to the 6$D_{\Delta t}$+QSO constraints,
the Hubble constant derived from the combined 6$D_{\Delta t}$+4$D_{L}$+QSO
data are a little higher than those values measured from the 6$D_{\Delta
t}$+QSO case. However, these differences are not statistically significant
given the error bars, as can be seen from the numerical results summarized in
Table. 1.

Non-flat cases.— As was revealed in the recent studies of Di Valentino et al.
(2020, 2021); Handley (2021), the discrepancy between the Hubble constant and
cosmic curvature measured locally and inferred from Planck highlights the
importance of considering non-flat cosmological models in this work. The
constraint results for three flat cosmological models are presented in Fig. 2,
in which we show the 2-D confidence regions(with 1$\sigma$ and 2$\sigma$
limits) as well as 1-D marginalized distributions from different data
combinations. The numerical results are also summarized in Table 1. One can
see from the upper panel of Fig. 2 that, in the concordance $\rm{\Lambda}CDM$
cosmology a stringent constraint on the Hubble constant could be obtained from
6 time-delay quasar data $D_{\Delta t}$ ($H_{0}=74.38^{+2.12}_{-2.34}$
$\mathrm{\leavevmode\nobreak\ km\leavevmode\nobreak\
s^{-1}\leavevmode\nobreak\ Mpc^{-1}}$). However, the MCMC chains failed to
converge for the other two model parameters, i.e., matter density and cosmic
curvature parameters ($\Omega_{m}$, $\Omega_{k}$). Such issue could be
appropriately addressed with stringent constraints produced by the combination
of 120 QSO sample and 6 time-delay lensed quasars, with the best-fitting
values with 68.3% confidence level for the three parameters:
$H_{0}=73.39^{+2.06}_{-2.11}$ $\mathrm{\leavevmode\nobreak\
km\leavevmode\nobreak\ s^{-1}\leavevmode\nobreak\ Mpc^{-1}}$,
$\Omega_{m}=0.274^{+0.083}_{-0.080}$, and
$\Omega_{k}=0.036^{+0.223}_{-0.214}$. With the combined data sets 6$D_{\Delta
t}$+QSO+4$D_{d}$, we also get stringent constraints on the model parameters
$H_{0}=73.78^{+1.99}_{-2.17}$ $\mathrm{\leavevmode\nobreak\
km\leavevmode\nobreak\ s^{-1}\leavevmode\nobreak\ Mpc^{-1}}$,
$\Omega_{m}=0.254^{+0.083}_{-0.074}$, and
$\Omega_{k}=0.100^{+0.214}_{-0.213}$. Compared with the previous results
obtained in other model-independent methods (Qi et al., 2019b), our analysis
results also demonstrate the strong degeneracy between the Hubble constant,
the matter density parameter and cosmic curvature, which would be effectively
broken by the combination of the latest observations of quasars, i.e., the
angular size of compact structure in radio quasars as standard rulers and the
time delays from gravitationally lensed quasars. The combination of the quasar
data sets, justified by their consistency within $1\sigma$, retains the same
correlation between $H_{0}$ and $\Omega_{k}$ as distinct samples of quasars.
For the determination of Hubble constant, our constraint in the framework of
non-flat $\Lambda$CDM cosmology is well consistent with the local Hubble
constant measurement from SH0ES collaboration (Riess et al., 2019). The
determination of $\Omega_{k}$ suggests no significant deviation from flat
spatial hypersurfaces, although favoring a somewhat positive value in the non-
flat $\Lambda$CDM case.

In the case of non-flat PEDE model, it can be clearly seen from the comparison
plots presented in Fig. 2 that there is a consistency between 6$D_{\Delta t}$,
6$D_{\Delta t}$+QSO, and 6$D_{\Delta t}$+4$D_{d}$+QSO data sets. Our results
confirm that the combination of compact structure in radio quasars and angular
diameter distances to the lenses (4$D_{d}$) could break the degeneracy between
cosmological parameters and lead to a more stringent constraints on all of the
cosmological parameters, which is the most unambiguous result of the current
lensed + unlensed quasar data set. Interestingly, the 6$D_{\Delta t}$+QSO data
generate a higher matter density parameter
$\Omega_{m}=0.315^{+0.094}_{-0.106}$ compared with other quasar samples
($\Omega_{m}=0.247^{+0.180}_{-0.131}$ for 6$D_{\Delta t}$ and
$\Omega_{m}=0.275^{+0.081}_{-0.077}$ for 6$D_{\Delta t}$+4$D_{d}$+QSO). For
the Hubble constant inferred from 6$D_{\Delta t}$, 6$D_{\Delta t}$+QSO, and
6$D_{\Delta t}$+4$D_{d}$+QSO data sets, it’s obvious that the combination of
6$D_{\Delta t}$ with 120 unlensed quasar data will result in a lower $H_{0}$,
but adding 4$D_{d}$ data to the combination would increase the median value of
$H_{0}$ comparing to the $H_{0}$ values obtained from 6$D_{\Delta t}$ alone.
The estimated values of the Hubble constant in the non-flat PEDE model are
between $H_{0}=71.67^{+16.75}_{-7.38}$ $\mathrm{\leavevmode\nobreak\
km\leavevmode\nobreak\ s^{-1}\leavevmode\nobreak\ Mpc^{-1}}$ and
$H_{0}=75.05^{+2.01}_{-2.19}$ $\mathrm{\leavevmode\nobreak\
km\leavevmode\nobreak\ s^{-1}\leavevmode\nobreak\ Mpc^{-1}}$ with 6$D_{\Delta
t}$ and 6$D_{\Delta t}$+4$D_{d}$+QSO data, which in broad terms agree very
well with the standard ones reported by the SH0ES collaboration (Riess et al.,
2019). For the determination of cosmic curvature in the non-flat PEDE model,
our results show that there is no significant evidence indicating its
deviation from zero (spatially flat geometry). More specifically, the two
quasar samples of 6$D_{\Delta t}$ and 6$D_{\Delta t}$+QSO favor closed
geometry ($\Omega_{k}=-0.027^{+0.340}_{-0.306}$,
$\Omega_{k}=-0.089^{+0.365}_{-0.276}$), while an open universe is favoured by
6$D_{\Delta t}$+4$D_{d}$+QSO data sets with
$\Omega_{k}=0.076^{+0.222}_{-0.201}$.

Finally, we perform a comparative analysis of the current lensed+unlensed
quasar data set in the non-flat DGP model. The $1\sigma$ and $2\sigma$
confidence level contours for parameter estimations are shown in Fig. 2, with
the marginalized best-fitting parameters and $1\sigma$ uncertainties
summarized in Table 1. Let us note that for such modified gravity model
gravity arising from the braneworld theory, the best-fitting matter density
parameter will be considerably shifted to a lower value. Our final assessments
of the matter density with the corresponding 1$\sigma$ uncertainties
($\Omega_{m}=0.251^{+0.186}_{-0.142}$, $\Omega_{m}=0.253^{+0.089}_{-0.081}$,
and $\Omega_{m}=0.253^{+0.086}_{-0.083}$ with 6$D_{\Delta t}$, 6$D_{\Delta
t}$+QSO, and 6$D_{\Delta t}$+4$D_{d}$+QSO data sets) are consistent with the
standard ones reported by other astrophysical probes, such as growth factors
combined with CMB+BAO+SNe Ia observations (Xia, 2009), Hubble parameter
combined with CMB+SNe Ia data (Fang et al., 2008), statistical analysis of
strong gravitational lensing systems (Cao et al., 2011b, 2012b; Ma et al.,
2019; Liu et al., 2020b), and the observations of compact structure in
intermediate-luminosity radio quasars (Cao et al., 2017; Liu et al., 2019,
2021b). For comparison, the corresponding fits on the parameter of
$\Omega_{r_{c}}$ are also displayed in Table 1. We remark here that the DGP
model will reduce to the concordance $\Lambda$CDM model when $\Omega_{rc}=0$.
Considering its non-vanishing value revealed in this analysis and previous
works ($\Omega_{r_{c}}\sim 0.14$), our fitting result shows the DGP model
fails to recover and is only a bit worse than the $\Lambda$CDM under the
current observational tests (see also Fang et al. (2008); Cao et al. (2011b)).
We also investigate how sensitive our results on $H_{0}$ and $\Omega_{k}$ are
on the choice of this cosmological model. On the one hand, the measured value
of $H_{0}$ with 1$\sigma$ uncertainty derived from the latest quasar sample,
$H_{0}=68.31^{+15.37}_{-10.45}$ $\mathrm{\leavevmode\nobreak\
km\leavevmode\nobreak\ s^{-1}\leavevmode\nobreak\ Mpc^{-1}}$ (6$D_{\Delta t}$)
and $H_{0}=66.26^{+13.57}_{-8.04}$ $\mathrm{\leavevmode\nobreak\
km\leavevmode\nobreak\ s^{-1}\leavevmode\nobreak\ Mpc^{-1}}$ (6$D_{\Delta
t}$+QSO) shows a perfect consistency with that derived from the recent Planck
CMB observation (Planck Collaboration et al., 2020). On the other hand, our
final assessments of the cosmic curvature with corresponding 1$\sigma$
uncertainty are $\Omega_{k}=0.023^{+0.316}_{-0.345}$ and
$\Omega_{k}=-0.050^{+0.328}_{-0.305}$ with the two quasar samples (6$D_{\Delta
t}$ and 6$D_{\Delta t}$+QSO), which are more consistent with flat spatial
hypersurfaces than what Qi et al. (2019b) found. From the joint analyses with
lensed+unlensed quasar data, we find that the Hubble constant and the spatial
curvature parameter are constrained to be $H_{0}=66.35^{+14.92}_{-8.10}$
$\mathrm{\leavevmode\nobreak\ km\leavevmode\nobreak\
s^{-1}\leavevmode\nobreak\ Mpc^{-1}}$ and
$\Omega_{k}=-0.047^{+0.348}_{-0.311}$, which furthermore confirms the above
conclusions.

Figure 3: Determination of Hubble constant from the combined data of lensed
quasars ($6D_{\Delta t}$ and $4D_{d}$) and unlensed radio quasars (QSO), in
the framework of six spatially flat (F) and non-flat (NF) cosmological models
($\Lambda$CDM, PEDE and DGP).  Figure 4: Graphical representation of the DIC
values for six spatially flat (F) and non-flat (NF) cosmological models
($\Lambda$CDM, PEDE and DGP), based on the combined data of lensed quasars
($6D_{\Delta t}$ and $4D_{d}$) and unlensed radio quasars (QSO). Figure 5:
Cosmological constraints on the flat PEDE model from the combined quasar data,
with different systematical uncertainties in the angular size measurements of
unlensed radio quasars (QSO).

Now let us remark on the $H_{0}$ and $\Omega_{k}$ measurements from the newly
compiled sample of ultra-compact structure in radio quasars and strong
gravitational lensing systems with quasars acting as background sources. The
$H_{0}$ determination from six spatially flat and non-flat cosmological models
($\Lambda$CDM, PEDE and DGP) are displayed in Fig. 3. A joint analysis of the
quasar sample (the time-delay measurements of six strong lensing systems, four
angular diameter distances to the lenses and 120 intermediate-luminosity
quasars calibrated as standard rulers) provides model-independent estimations
of the Hubble constant $H_{0}$, which is strongly consistent with that derived
from the local distance ladder by SH0ES collaboration in the $\Lambda$CDM and
PEDE model. However, in the framework of a DGP cosmology (especially for the
flat universe), the measured Hubble constant is in good agreement with that
derived from the the recent Planck 2018 results. Meanwhile, our results also
demonstrate that zero spatial curvature is supported by the current lensed and
unlensed quasar observations and there is no significant deviation from a flat
universe. Finally, our results have provided independent evidence for the
accelerated expansion of the Universe, with the existence of dominant dark
energy density ($\Omega_{m}\sim 0.30$) in the framework of six cosmologies
classified into different categories. These are the most unambiguous result of
the current quasar data sets.

Model comparison.— The values of DIC for all models are reported in Table 1.
According to the number of model parameters, these six cosmological models can
be divided into two classes: the two-parameter models including the flat
$\Lambda$CDM, flat PEDE and flat DGP; and the three-parameter models including
non-flat $\Lambda$CDM, non-flat PEDE and non-flat DGP. For the two-parameter
models, the flat PEDE provides the smallest information criterion result
(DIC=390.68) among all of the flat cosmological models. However, we note that
the difference of DIC between the flat PEDE and flat $\Lambda$CDM model is
only $\Delta$DIC=0.65, which means these two models are comparable to each
other according to this criterion. The flat DGP model is the worst one to
explain the current lensed and unlensed quasar observations, since the DIC
value it yields is the largest among the two-parameter models
($\Delta$DIC=2.99). As for the three-parameter models, the DIC results show
that the non-flat $\Lambda$CDM model is still not the best one. The non-flat
PEDE model, which is a little bit better than non-flat $\Lambda$CDM, performs
best in explaining the quasar data (DIC=393.45), with positive evidence
against non-flat DGP ($\Delta$DIC$=1.82$).

We also provide a graphical representation of the DIC results in Fig. 4 which
directly shows the results in the IC test for each model. Out of all the
candidate models, it is obvious that the flat PEDE and flat $\Lambda$CDM are
the two most favored models in the data combination of lensed+unlensed
quasars, Following them are the flat DGP, non-flat $\Lambda$CDM, and non-flat
PEDE that give comparable fits to the data. According to the DIC results, the
most disfavored model is non-flat DGP, with strong evidence against such
cosmological scenario among the six models we study here.

##  5 Conclusions

In this paper, we analyze six spatially flat and non-flat cosmological models
($\Lambda$CDM, PEDE and DGP) using a newly compiled sample of ultra-compact
structure in radio quasars and strong gravitational lensing systems with
quasars acting as background sources. This study is strongly motivated by the
need for revisiting the Hubble constant, spatial curvature and dark energy
dynamics in the framework of different cosmological models of interests, and
searching for implications for the non-flat Universe and extensions of the
standard cosmological model (the spatially flat $\Lambda$CDM). The inclusion
of such a newly compiled quasar sample in the cosmological analysis is crucial
to this aim as it will extend the Hubble diagram to a high-redshift range, in
which predictions from different cosmological models can be distinguished
(Capozziello et al., 2006). From the constraints derived using the updated
observations of quasars, we can identify some relatively model-independent
features.

In all cosmological models, the cosmological parameters obtained from distinct
quasar samples are consistent and the combination of the latest observations
of quasars, i.e., the time-delay measurements of six strong lensing systems,
four angular diameter distances to the lenses and 120 intermediate-luminosity
quasars calibrated as standard rulers, would break the degeneracy between the
Hubble constant and other cosmological parameters. The lensed quasar
($6D_{\Delta t}$ and $4D_{d}$) and unlensed radio quasar (QSO) data
combination produces the most reliable constraints. In particular, for most of
cosmological model we study (the flat $\Lambda$CDM, non-flat $\Lambda$CDM,
flat PEDE, and non-flat PEDE models), the derived matter density parameter is
completely consistent with $\Omega_{m}\sim 0.30$ in all the data sets, as
expected by the latest cosmological observations (Wong et al., 2020; Cao et
al., 2017; Liu et al., 2019; Li et al., 2017; Cao et al., 2020; Lian et al.,
2021) and Planck Collaboration results ($\Omega_{m}=0.3103\pm 0.0057$) (Planck
Collaboration et al., 2020). Nevertheless, the DGP model in both flat and non-
flat cases shows a deviation from this prediction, with statistical lower
values of $\Omega_{m}=0.246^{+0.043}_{-0.037}$ and
$\Omega_{m}=0.253^{+0.086}_{-0.083}$ for the combined sample of 6$D_{\Delta
t}$+4$D_{d}$+QSO. A joint analysis of the quasar sample provides model-
independent estimations of the Hubble constant $H_{0}$, which is strongly
consistent with that derived from the local distance ladder by SH0ES
collaboration (Riess et al., 2019) in the $\Lambda$CDM and PEDE model.
However, in the framework of a DGP cosmology (especially for the flat
universe), the measured value of $H_{0}$ with 1$\sigma$ uncertainty, which is
3.6$\sigma$ lower than the statistical estimates of the SH0ES results,
demonstrates a perfect agreement with that derived from the recent Plank CMB
observations (Planck Collaboration et al., 2020). Our findings also confirm
the flatness of our universe (Collett et al., 2019; Wei & Melia, 2020; Qi et
al., 2021), i.e., the most unambiguous result of the current lensed and
unlensed quasar observations, although there is some room for a little spatial
curvature energy density in the non-flat $\Lambda$CDM, non-flat PEDE and non-
flat DGP cases. Finally, we statistically evaluate which model is more
consistent with the observational quasar data. Concerning the ranking of
competing dark energy models, the flat PEDE is the most favored model out of
all the candidate models, while the non-flat DGP is substantially penalized by
the DIC criteria. However, our analysis still does not rule out dark energy
being a cosmological constant and non-flat spatial hypersurfaces.

Considering the possible controversy around the systematics of the observed
angular sizes of compact radio quasars (Kellermann, 1993), the other
reasonable strategy to quantify the such effect of systematics is taking
$\sigma_{sys}$ as an additional free parameter, which should be fitted
simultaneously with the cosmological parameters of interest. Such strategy has
been extensively applied in the derivation of quasar distances based on the
non-linear relation between their UV and X-ray fluxes, based on the largest
quasar sample consists of 12,000 objects with both X-ray and UV observations
(Risaliti & Lusso, 2019). In this paper, we perform a sensitivity analysis by
introducing an overall 5% and 20% systematical uncertainty to the angular size
measurements of compact radio quasars, in order to investigate how the
cosmological constraints on flat PEDE is altered by different choices of
systematics. In the framework of such error strategy in the construction of
posterior likelihood $\mathcal{L}_{QSO}$, the matter density parameter and
Hubble constant respectively change to $\Omega_{m}=0.300^{+0.040}_{-0.034}$,
$H_{0}=74.82^{+1.81}_{-2.04}$ $\mathrm{\leavevmode\nobreak\
km\leavevmode\nobreak\ s^{-1}\leavevmode\nobreak\ Mpc^{-1}}$ (5% systematical
uncertainty), and $\Omega_{m}=0.331^{+0.048}_{-0.040}$,
$H_{0}=74.93^{+1.83}_{-2.05}$ $\mathrm{\leavevmode\nobreak\
km\leavevmode\nobreak\ s^{-1}\leavevmode\nobreak\ Mpc^{-1}}$ (20% systematical
uncertainty). The comparison of the resulting constraints on $\Omega_{m}$ and
$H_{0}$ based on different systematical uncertainties is shown in Fig. 5. In
general, one can easily check that the derived value of $\Omega_{m}$ is more
sensitive to the systematical uncertainties of angular size measurements,
i.e., a larger systematical uncertainties will shift the matter density
parameter to a relatively higher value. This illustrates the importance of a
larger quasar sample from future VLBI observations based on better uv-coverage
(Pushkarev & Kovalev, 2015). As a final remark, from the observational point
of view, one can see that the 120 unlensed radio quasars have perfect coverage
of source redshifts in six strong lensing systems ($z\sim 3$). Given the
usefulness of compact radio quasars and strongly lensed quasars acting as
standard rulers at high redshifts, we pin our hopes on a large amount of
intermediate-luminosity radio quasars detected by future VLBI surveys at
different frequencies (Pushkarev & Kovalev, 2015), and strongly lensed quasars
with well-measured time delays discovered by future surveys of Large Synoptic
Survey Telescope (Collett, 2015).

###### Acknowledgements.

This work was supported by the National Natural Science Foundation of China
under Grants Nos. 12203009, 12021003, 11690023, 11920101003 and A2020205002;
the Strategic Priority Research Program of the Chinese Academy of Sciences,
Grant No. XDB23000000; the Interdiscipline Research Funds of Beijing Normal
University; the China Manned Space Project (Nos. CMS-CSST-2021-B01 and CMS-
CSST-2021-A01); and the CAS Project for Young Scientists in Basic Research
under Grant No. YSBR-006. X. Li was supported by National Natural Science
Foundation of China under Grants No. 1200300, Hebei NSF under Grant No.
A2020205002 and the fund of Hebei Normal University under Grants No. L2020B02.

## References

  * Alam et al. (2017) Alam, U., Bag, S., & Sahni, V. 2017, Phys. Rev. D, 95, 023524. doi:10.1103/PhysRevD.95.023524
  * Bag et al. (2022) Bag, S., Shafieloo, A., Liao, K., et al. 2022, ApJ, 927, 191. doi:10.3847/1538-4357/ac51cb
  * Birrer et al. (2016) Birrer, S., Amara, A., & Refregier, A. 2016, J. Cosmology Astropart. Phys., 2016, 020. doi:10.1088/1475-7516/2016/08/020
  * Birrer et al. (2019) Birrer, S., Treu, T., Rusu, C. E., et al. 2019, MNRAS, 484, 4726. doi:10.1093/mnras/stz200
  * Briffa et al. (2020) Briffa, R., Capozziello, S., Said, J. L., et al. 2020, Classical and Quantum Gravity, 38, 055007. doi:10.1088/1361-6382/abd4f5
  * Caldera-Cabral et al. (2009) Caldera-Cabral, G., Maartens, R., & Ureña-López, L. A. 2009, Phys. Rev. D, 79, 063518. doi:10.1103/PhysRevD.79.063518
  * Cao et al. (2011a) Cao, S., Liang, N., & Zhu, Z.-H. 2011a, MNRAS, 416, 1099. doi:10.1111/j.1365-2966.2011.19105.x
  * Cao et al. (2011b) Cao, S., Zhu, Z.-H., & Zhao, R. 2011b, Phys. Rev. D, 84, 023005. doi:10.1103/PhysRevD.84.023005
  * Cao et al. (2012a) Cao, S., Pan, Y., Biesiada, M., et al. 2012a, J. Cosmology Astropart. Phys., 2012a, 016. doi:10.1088/1475-7516/2012/03/016
  * Cao et al. (2012b) Cao, S., Covone, G., & Zhu, Z.-H. 2012b, ApJ, 755, 31. doi:10.1088/0004-637X/755/1/31
  * Cao & Zhu (2014) Cao, S. & Zhu, Z.-H. 2014, Phys. Rev. D, 90, 083006. doi:10.1103/PhysRevD.90.083006
  * Cao et al. (2015) Cao, S., Biesiada, M., Gavazzi, R., et al. 2015, ApJ, 806, 185. doi:10.1088/0004-637X/806/2/185
  * Cao et al. (2017) Cao, S., Zheng, X., Biesiada, M., et al. 2017, A&A, 606, A15. doi:10.1051/0004-6361/201730551
  * Cao et al. (2018) Cao, S., Biesiada, M., Qi, J., et al. 2018, European Physical Journal C, 78, 749. doi:10.1140/epjc/s10052-018-6197-y
  * Cao et al. (2019) Cao, S., Qi, J., Biesiada, M., et al. 2019, Physics of the Dark Universe, 24, 100274. doi:10.1016/j.dark.2019.100274
  * Cao et al. (2020) Cao, S., Ryan, J., & Ratra, B. 2020, MNRAS, 497, 3191. doi:10.1093/mnras/staa2190
  * Capozziello et al. (2006) Capozziello, S., Cardone, V. F., & Troisi, A. 2006, J. Cosmology Astropart. Phys., 2006, 001. doi:10.1088/1475-7516/2006/08/001
  * Chen et al. (2019) Chen, G. C.-F., Fassnacht, C. D., Suyu, S. H., et al. 2019, MNRAS, 490, 1743. doi:10.1093/mnras/stz2547
  * Collett (2015) Collett, T. E. 2015, ApJ, 811, 20. doi:10.1088/0004-637X/811/1/20
  * Collett et al. (2019) Collett, T., Montanari, F., & Räsänen, S. 2019, Phys. Rev. Lett., 123, 231101. doi:10.1103/PhysRevLett.123.231101
  * Deffayet et al. (2002a) Deffayet, C., Dvali, G., & Gabadadze, G. 2002a, Phys. Rev. D, 65, 044023. doi:10.1103/PhysRevD.65.044023
  * Deffayet et al. (2002b) Deffayet, C., Landau, S. J., Raux, J., et al. 2002b, Phys. Rev. D, 66, 024019. doi:10.1103/PhysRevD.66.024019
  * Di Valentino et al. (2018) Di Valentino, E., Melchiorri, A., Fantaye, Y., et al. 2018, Phys. Rev. D, 98, 063508. doi:10.1103/PhysRevD.98.063508
  * Di Valentino et al. (2020) Di Valentino, E., Melchiorri, A., & Silk, J. 2020, Nature Astronomy, 4, 196. doi:10.1038/s41550-019-0906-9
  * Di Valentino et al. (2021) Di Valentino, E., Melchiorri, A., & Silk, J. 2021, ApJ, 908, L9. doi:10.3847/2041-8213/abe1c4
  * Ding et al. (2021) Ding, X., Treu, T., Birrer, S., et al. 2021, MNRAS, 503, 1096. doi:10.1093/mnras/stab484
  * Dvali et al. (2000) Dvali, G., Gabadadze, G., & Porrati, M. 2000, Physics Letters B, 485, 208. doi:10.1016/S0370-2693(00)00669-9
  * Dvali et al. (2001) Dvali, G., Gabadadze, G., & Shifman, M. 2001, Physics Letters B, 497, 271. doi:10.1016/S0370-2693(00)01329-0
  * Fang et al. (2008) Fang, W., Wang, S., Hu, W., et al. 2008, Phys. Rev. D, 78, 103509. doi:10.1103/PhysRevD.78.103509
  * Follin & Knox (2018) Follin, B. & Knox, L. 2018, MNRAS, 477, 4534. doi:10.1093/mnras/sty720
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., et al. 2013, PASP, 125, 306. doi:10.1086/670067
  * Giannantonio et al. (2008) Giannantonio, T., Song, Y.-S., & Koyama, K. 2008, Phys. Rev. D, 78, 044017. doi:10.1103/PhysRevD.78.044017
  * Handley (2021) Handley, W. 2021, Phys. Rev. D, 103, L041301. doi:10.1103/PhysRevD.103.L041301
  * Hashim et al. (2021a) Hashim, M., El Hanafy, W., Golovnev, A., et al. 2021a, J. Cosmology Astropart. Phys., 2021, 052. doi:10.1088/1475-7516/2021/07/052
  * Hashim et al. (2021b) Hashim, M., El-Zant, A. A., El Hanafy, W., et al. 2021b, J. Cosmology Astropart. Phys., 2021, 053. doi:10.1088/1475-7516/2021/07/053
  * Jee et al. (2016) Jee, I., Komatsu, E., Suyu, S. H., et al. 2016, J. Cosmology Astropart. Phys., 2016, 031. doi:10.1088/1475-7516/2016/04/031
  * Jee et al. (2019) Jee, I., Suyu, S. H., Komatsu, E., et al. 2019, Science, 365, 1134. doi:10.1126/science.aat7371
  * Jiménez Cruz & Escamilla-Rivera (2021) Jiménez Cruz, N. M. & Escamilla-Rivera, C. 2021, European Physical Journal Plus, 136, 51. doi:10.1140/epjp/s13360-020-00950-3
  * Koyama (2016) Koyama, K. 2016, Reports on Progress in Physics, 79, 046902. doi:10.1088/0034-4885/79/4/046902
  * Karwal & Kamionkowski (2016) Karwal, T. & Kamionkowski, M. 2016, Phys. Rev. D, 94, 103523. doi:10.1103/PhysRevD.94.103523
  * Kellermann (1993) Kellermann, K. I. 1993, Nature, 361, 134. doi:10.1038/361134a0
  * Li et al. (2017) Li, X., Cao, S., Zheng, X., et al. 2017, arXiv:1708.08867
  * Li & Shafieloo (2019) Li, X. & Shafieloo, A. 2019, ApJ, 883, L3. doi:10.3847/2041-8213/ab3e09
  * Li & Shafieloo (2020) Li, X. & Shafieloo, A. 2020, ApJ, 902, 58. doi:10.3847/1538-4357/abb3d0
  * Li et al. (2019) Li, X., Shafieloo, A., Sahni, V., et al. 2019, ApJ, 887, 153. doi:10.3847/1538-4357/ab535d
  * Liao et al. (2019a) Liao, K., Shafieloo, A., Keeley, R. E., et al. 2019, ApJ, 886, L23. doi:10.3847/2041-8213/ab5308
  * Liao et al. (2020) Liao, K., Shafieloo, A., Keeley, R. E., et al. 2020, ApJ, 895, L29. doi:10.3847/2041-8213/ab8dbb
  * Liao (2021) Liao, K. 2021, ApJ, 906, 26. doi:10.3847/1538-4357/abc876
  * Lian et al. (2021) Lian, Y., Cao, S., Biesiada, M., et al. 2021, MNRAS, 505, 2111. doi:10.1093/mnras/stab1373
  * Liu et al. (2019) Liu, T., Cao, S., Zhang, J., et al. 2019, ApJ, 886, 94. doi:10.3847/1538-4357/ab4bc3
  * Liu et al. (2020a) Liu, T., Cao, S., Biesiada, M., et al. 2020a, ApJ, 899, 71. doi:10.3847/1538-4357/aba0b6
  * Liu et al. (2020b) Liu, T., Cao, S., Zhang, J., et al. 2020b, MNRAS, 496, 708. doi:10.1093/mnras/staa1539
  * Liu et al. (2021a) Liu, T., Cao, S., Zhang, S., et al. 2021a, European Physical Journal C, 81, 903. doi:10.1140/epjc/s10052-021-09713-5
  * Liu et al. (2021b) Liu, T., Cao, S., Biesiada, M., et al. 2021b, MNRAS, 506, 2181. doi:10.1093/mnras/stab1868
  * Liddle (2007) Liddle, A. R. 2007, MNRAS, 377, L74. doi:10.1111/j.1745-3933.2007.00306.x
  * Lyu et al. (2020) Lyu, M.-Z., Haridasu, B. S., Viel, M., et al. 2020, ApJ, 900, 160. doi:10.3847/1538-4357/aba756
  * Ma et al. (2017) Ma, Y., Zhang, J., Cao, S., et al. 2017, European Physical Journal C, 77, 891. doi:10.1140/epjc/s10052-017-5461-x
  * Ma et al. (2019) Ma, Y.-B., Cao, S., Zhang, J., et al. 2019, European Physical Journal C, 79, 121. doi:10.1140/epjc/s10052-019-6630-x
  * Melia (2018) Melia, F. 2018, EPL (Europhysics Letters), 123, 39001. doi:10.1209/0295-5075/123/39001
  * Multamäki et al. (2003) Multamäki, T., Gaztañaga, E., & Manera, M. 2003, MNRAS, 344, 761. doi:10.1046/j.1365-8711.2003.06880.x
  * Perlick (1990a) Perlick, V. 1990a, Classical and Quantum Gravity, 7, 1319. doi:10.1088/0264-9381/7/8/011
  * Perlick (1990b) Perlick, V. 1990b, Classical and Quantum Gravity, 7, 1849. doi:10.1088/0264-9381/7/10/016
  * Perivolaropoulos & Skara (2021) Perivolaropoulos, L. & Skara, F. 2021, Phys. Rev. D, 104, 123511. doi:10.1103/PhysRevD.104.123511
  * Planck Collaboration et al. (2016) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, A&A, 594, A13. doi:10.1051/0004-6361/201525830
  * Planck Collaboration et al. (2020) Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020, A&A, 641, A6. doi:10.1051/0004-6361/201833910
  * Pushkarev & Kovalev (2015) Pushkarev, A. B. & Kovalev, Y. Y. 2015, MNRAS, 452, 4274. doi:10.1093/mnras/stv1539
  * Qi et al. (2017) Qi, J.-Z., Cao, S., Biesiada, M., et al. 2017, European Physical Journal C, 77, 502. doi:10.1140/epjc/s10052-017-5069-1
  * Qi et al. (2019a) Qi, J., Cao, S., Biesiada, M., et al. 2019a, Phys. Rev. D, 100, 023530. doi:10.1103/PhysRevD.100.023530
  * Qi et al. (2019b) Qi, J.-Z., Cao, S., Zhang, S., et al. 2019b, MNRAS, 483, 1104. doi:10.1093/mnras/sty3175
  * Qi et al. (2021) Qi, J.-Z., Zhao, J.-W., Cao, S., et al. 2021, MNRAS, 503, 2179. doi:10.1093/mnras/stab638
  * Rathna Kumar et al. (2015) Rathna Kumar, S., Stalin, C. S., & Prabhu, T. P. 2015, A&A, 580, A38. doi:10.1051/0004-6361/201423977
  * Ren et al. (2021) Ren, X., Wong, T. H. T., Cai, Y.-F., et al. 2021, Physics of the Dark Universe, 32, 100812. doi:10.1016/j.dark.2021.100812
  * Riess et al. (2007) Riess, A. G., Strolger, L.-G., Casertano, S., et al. 2007, ApJ, 659, 98. doi:10.1086/510378
  * Riess et al. (2019) Riess, A. G., Casertano, S., Yuan, W., et al. 2019, ApJ, 876, 85. doi:10.3847/1538-4357/ab1422
  * Risaliti & Lusso (2019) Risaliti, G. & Lusso, E. 2019, Nature Astronomy, 3, 272. doi:10.1038/s41550-018-0657-z
  * Ryan (2021) Ryan, J. 2021, J. Cosmology Astropart. Phys., 2021, 051. doi:10.1088/1475-7516/2021/08/051
  * Rusu et al. (2020) Rusu, C. E., Wong, K. C., Bonvin, V., et al. 2020, MNRAS, 498, 1440. doi:10.1093/mnras/stz3451
  * Sonnenfeld (2021) Sonnenfeld, A. 2021, A&A, 656, A153. doi:10.1051/0004-6361/202142062
  * Spiegelhalter et al. (2002) Spiegelhalter, D. J., Best, N. G., Carlin, B. P. & Van Der Linde, A. 2020, Journal of the royal statistical society: Series b (statistical methodology), 64, 583. doi:10.1111/1467-9868.00353
  * Suyu et al. (2010) Suyu, S. H., Marshall, P. J., Auger, M. W., et al. 2010, ApJ, 711, 201. doi:10.1088/0004-637X/711/1/201
  * Suyu et al. (2013) Suyu, S. H., Auger, M. W., Hilbert, S., et al. 2013, ApJ, 766, 70. doi:10.1088/0004-637X/766/2/70
  * Suyu et al. (2014) Suyu, S. H., Treu, T., Hilbert, S., et al. 2014, ApJ, 788, L35. doi:10.1088/2041-8205/788/2/L35
  * Taubenberger et al. (2019) Taubenberger, S., Suyu, S. H., Komatsu, E., et al. 2019, A&A, 628, L7. doi:10.1051/0004-6361/201935980
  * Treu & Marshall (2016) Treu, T. & Marshall, P. J. 2016, A&A Rev., 24, 11. doi:10.1007/s00159-016-0096-8
  * Väliviita et al. (2010) Väliviita, J., Maartens, R., & Majerotto, E. 2010, MNRAS, 402, 2355. doi:10.1111/j.1365-2966.2009.16115.x
  * Vavryčuk & Kroupa (2020) Vavryčuk, V. & Kroupa, P. 2020, MNRAS, 497, 378. doi:10.1093/mnras/staa1936
  * Wang et al. (2020) Wang, B., Qi, J.-Z., Zhang, J.-F., et al. 2020, ApJ, 898, 100. doi:10.3847/1538-4357/ab9b22
  * Wei & Melia (2020) Wei, J.-J. & Melia, F. 2020, ApJ, 897, 127. doi:10.3847/1538-4357/ab959b
  * Wong et al. (2017) Wong, K. C., Suyu, S. H., Auger, M. W., et al. 2017, MNRAS, 465, 4895. doi:10.1093/mnras/stw3077
  * Wong et al. (2020) Wong, K. C., Suyu, S. H., Chen, G. C.-F., et al. 2020, MNRAS, 498, 1420. doi:10.1093/mnras/stz3094
  * Xia (2009) Xia, J.-Q. 2009, Phys. Rev. D, 79, 103527. doi:10.1103/PhysRevD.79.103527
  * Xu (2014) Xu, L. 2014, J. Cosmology Astropart. Phys., 2014, 048. doi:10.1088/1475-7516/2014/02/048
  * Xu & Wang (2010) Xu, L. & Wang, Y. 2010, Phys. Rev. D, 82, 043503. doi:10.1103/PhysRevD.82.043503
  * Yang et al. (2019) Yang, W., Pan, S., Vagnozzi, S., et al. 2019, J. Cosmology Astropart. Phys., 2019, 044. doi:10.1088/1475-7516/2019/11/044
  * Zheng et al. (2017) Zheng, X., Biesiada, M., Cao, S., et al. 2017, J. Cosmology Astropart. Phys., 2017, 030. doi:10.1088/1475-7516/2017/10/030