# Quantum sensing to suppress systematic errors with measurements in the Zeno
regime

Alisa Shimada  Research Center for Emerging Computing Technologies, National
institute of Advanced Industrial Science and Technology (AIST), Central2,
1-1-1 Umezono, Tsukuba, Ibaraki 305-8568, Japan  Hideaki Hakoshima  Research
Center for Emerging Computing Technologies, National institute of Advanced
Industrial Science and Technology (AIST), Central2, 1-1-1 Umezono, Tsukuba,
Ibaraki 305-8568, Japan  Suguru Endo  NTT Secure Platform Laboratories, NTT
Corporation, Musashino 180-8585, Japan  Kaoru Yamamoto  NTT Secure Platform
Laboratories, NTT Corporation, Musashino 180-8585, Japan  Yuichiro Matsuzaki
<EMAIL_ADDRESS>Research Center for Emerging Computing
Technologies, National institute of Advanced Industrial Science and Technology
(AIST), Central2, 1-1-1 Umezono, Tsukuba, Ibaraki 305-8568, Japan

###### Abstract

Quantum magnetic field sensing is an important technology for material science
and biology. Although experimental imperfections affect the sensitivity,
repetitions of the measurements decrease the estimation uncertainty by a
square root of the total number of measurements if there are only statistical
errors. However, it is difficult to precisely characterize the coherence time
of the system because it fluctuates in time in realistic conditions, which
induces systematic errors. In this case, due to residual bias of the measured
values, estimation uncertainty cannot be lowered than a finite value even in
the limit of the infinite number of measurements. On the basis of the fact
that the decoherence dynamics in the so-called Zeno regime are not significant
compared to other regimes, we propose a novel but very simple protocol to use
measurements in the Zeno regime for reducing systematic errors. Our scheme
allows the estimation uncertainty $\delta^{2}\omega$ to scale as $L^{1/4}$
where $L$ denotes the number of measurements even when we cannot precisely
characterize the coherence time.

##  I Introduction

Measurements of small magnetic fields is an essential technique in several
fields such as biology and material science. Recently, quantum systems have
been used for sensitive magnetic field sensing. External magnetic fields shift
the resonant frequency of the quantum system, and measurements of a phase
shift allow us to detect magnetic fields via the frequency shift Degen  _et
al._ (2017). Nitrogen vacancies in diamond Taylor  _et al._ (2008); Maze  _et
al._ (2008); Balasubramanian  _et al._ (2008), superconducting flux qubits Bal
_et al._ (2012); Toida  _et al._ (2019); Budoyo  _et al._ (2020), and
optically pumped atoms Kominis  _et al._ (2003) have been used for such
quantum magnetic field sensing.

Experimental imperfections such as decoherence typically deteriorate the
sensitivity Degen  _et al._ (2017). For the detection of the small magnetic
fields, it is essential to decrease the estimation uncertainty. When
experimental parameters such as g factor and coherence time of the qubits are
known, we can estimate the magnetic fields from the experimental results
without bias. In this case, the estimation uncertainty decreases as the total
number of measurements increases, because there are only statistical errors.

However, in the actual experiments, the coherence time of the system generally
fluctuates in time Yan  _et al._ (2016); Abdurakhimov  _et al._ (2019), which
makes it difficult to know the accurate value of the coherence time. In such a
case, the lack of knowledge of the coherence time induces systematic errors in
the estimation. Under the effect of the systematic errors, the estimation
uncertainty is lower bounded by a specific value even in the limit of the
large number of measurements Wolf  _et al._ (2015); Budoyo  _et al._ (2020).
This could be a problem to measure small magnetic fields.

On the other hand, quantum mechanics predicts that frequent measurements with
short time intervals can suppress the evolution of quantum states Misra and
Sudarshan (1977); Facchi and Pascazio (2008). In the so-called Zeno-regime,
‘survival probability’ that the state remains in the initial state after a
short time $t$ scales as $P_{\rm{suv}}\simeq 1-\Gamma^{2}t^{2}$, where
$\Gamma$ denotes a decay rate Schulman  _et al._ (1994); Nakazato  _et al._
(1996) Erez  _et al._ (2004); Facchi and Pascazio (2008); Nakazato  _et al._
(1996); Koshino and Shimizu (2005). This phenomenon has been experimentally
demonstrated in several systems Itano  _et al._ (1990); Kwiat  _et al._
(1995); Fischer  _et al._ (2001); Streed  _et al._ (2006); Helmer  _et al._
(2009); Wolters  _et al._ (2013); Kakuyanagi  _et al._ (2015); Kondo  _et al._
(2016); Kalb  _et al._ (2016).

Here, we propose to use measurements in the Zeno regime for quantum sensing to
suppress the systematic errors when we cannot characterize the precise value
of the coherence time. In the Zeno regime, the system coupled with the
environment shows a quadratic decay, and the decay dynamics can be much slower
for a short time region than the other regions. This means that the qubit
state is not significantly affected by the decoherence in the short time
regime. We show that, by performing the measurements in the Zeno regime, we
can suppress the effect of the systematic errors induced by the lack of
knowledge of the coherence time. Moreover, in our scheme, the estimation
uncertainty can be arbitrarily small by increasing the number of measurements,
even when we do not know the precise value of the coherence time. Although the
measurements in the Zeno regime was discussed in entanglement enhanced sensing
Matsuzaki  _et al._ (2011); Chin  _et al._ (2012); Macieszczak (2015), we
firstly utilize this concept to suppress the systematic errors.

##  II quantum sensing

Let us review the standard quantum sensing, focusing on the case of a qubit.
We explain two cases, without and with systematic errors.

###  II.1 Quantum sensing without systematic errors

First, we review the sensing without systematic errors. Here, for simplicity,
we ignore any experimental imperfection such as decoherence. We consider the
Hamiltonian as follows

$\displaystyle H=\sum_{j=1}^{L}\frac{\omega}{2}\hat{\sigma}^{(j)}_{z},$ (1)

where $\omega$ denotes a frequency of the qubit and $\hat{\sigma}_{z}$ denote
a Pauli matrix spanned by $|1\rangle$ and $|-1\rangle$. We assume that
$\omega$ has a linear dependence on applied magnetic fields $B$. This means
that, if we know the value of the $\omega$, we can determine the value of the
magnetic fields. The sensing protocol can be performed as following. First, we
prepare the state of
$\bigotimes_{j=1}^{L}|+\rangle_{j}=\frac{1}{\sqrt{2}}(|1\rangle_{j}+|-1\rangle_{j})$.
Second, let the state evolve by the Hamiltonian to obtain
$\bigotimes_{j=1}^{L}|\psi(t)\rangle_{j}$ where
$|\psi(t)\rangle_{j}=\frac{1}{\sqrt{2}}(|1\rangle_{j}+e^{i\omega
t}|-1\rangle_{j})$. Third, perform a projective measurement with
$\hat{\mathcal{P}}_{j}=(\openone+\hat{\sigma}^{(j)}_{y})/2$ for $j=1,2,\cdots
L$ where the probability to have this projection is $P^{\prime}=(1+\sin\omega
t)/2\simeq(1+\omega t)/2$ for weak magnetic fields with $|\omega|t\ll 1$.
Fourth, repeat these three steps $N$ times. Finally, we obtain $L$ measurement
results $\\{s_{m}\\}_{m=1}^{L}$ from $L$ qubits, and estimate the frequency
$\omega$ based on the measurement results. More specifically, we obtain the
estimate value as $\omega_{\rm{est}}=(2(\sum_{m=1}^{L}\frac{s_{m}}{LN})-1)/t$.
We can calculate the uncertainty of the estimation as follows.

$\displaystyle\delta^{2}\omega$ $\displaystyle=$
$\displaystyle\overline{(\omega-\omega_{\rm{est}})^{2}}$ (2) $\displaystyle=$
$\displaystyle\frac{P^{\prime}(1-P^{\prime})}{|\frac{dP^{\prime}}{d\omega}|^{2}LN}$
$\displaystyle=$ $\displaystyle\frac{1}{t^{2}LN}$

where the overline denotes the statistical average. In our paper, we basically
consider a case of $N=1$. The uncertainty scales as
$\delta^{2}\omega=\Theta(L^{-2})$, and so we can decrease the uncertainty by
increasing the number of qubits

###  II.2 Quantum sensing with systematic errors

We review the sensing with systematic errors. Suppose that there are some
unknown experimental imperfections such as systematic errors, and the actual
probability $P^{\prime}$ to have a projection of
$\hat{\mathcal{P}}_{j}=(\openone+\hat{\sigma}^{(j)}_{y})/2$ for $j=1,2,\cdots
L$ is different from what an experimentalist expected to be true. Let us
define such a wrong (true) probability as $P$ ($P^{\prime}$). We assume that
these probabilities have a linear dependence on $\omega$ such as $P=x+y\omega$
and $P^{\prime}=x^{\prime}+y^{\prime}\omega$. In this case, the uncertainty of
the estimation for a small $\omega$ can be calculated as follows Sugiyama
(2015); Takeuchi  _et al._ (2019).

$\displaystyle\delta^{2}\omega=\frac{1}{y^{2}}\left[\frac{P^{\prime}(1-P^{\prime})}{L}+(x-x^{\prime})^{2}\right]$
(3)

It is worth mentioning that, even in a limit of large $L$, the uncertainty has
a finite non-zero value as follows:

$\displaystyle\lim_{L\rightarrow\infty}\delta^{2}\omega=\frac{(x-x^{\prime})^{2}}{y^{2}}.$
(4)

This finite value comes from systematic errors due to the lack of knowledge of
the coherence time. The goal of this paper is to suppress such systematic
errors.

##  III Quantum sensing with measurements in Zeno regime

In this section, we explain our scheme to suppress the systematic errors with
measurements in Zeno regime.

###  III.1 Setup

Let us explain the setup of our scheme. We consider a three level system (or a
spin-1 system) , which describes an NV center in diamond Taylor  _et al._
(2008); Maze  _et al._ (2008); Balasubramanian  _et al._ (2008) or a
capacitively shunted flux qubit You  _et al._ (2007); Yan  _et al._ (2016).
Although there are three levels such as $|1\rangle$, $|-1\rangle$, and
$|0\rangle$, we mainly use $|1\rangle$ and $|-1\rangle$ for magnetic field
sensing. However, due to an energy relaxation, we have unwanted population of
the state $|0\rangle$. Throughout this paper, although this system is not a
qubit because it has a small population of the third level, we call this a
qubit because we mainly use $|1\rangle$ and $|-1\rangle$. In the magnetic
field sensing scheme, there is a state preparation step to have an initial
state of
$\bigotimes_{j}^{L}|+\rangle_{j}=\bigotimes_{j}^{L}\frac{1}{\sqrt{2}}(|1\rangle_{j}+|-1\rangle_{j})$,
an exposure step to interact the spin-1 systems with target magnetic fields,
and a readout step to measure the state. We assume that the state preparation
time and readout time is much longer than the exposure time. After interacting
the $L$ qubits with the magnetic fields for sensing, we have the following
state.

$\displaystyle\rho$ $\displaystyle=$
$\displaystyle\bigotimes_{j=1}^{L}\rho_{j}$ $\displaystyle\rho_{j}$
$\displaystyle=$
$\displaystyle(1-\epsilon)|\psi(t)\rangle_{j}\langle\psi(t)|_{j}+\epsilon|0\rangle_{j}\langle
0|_{j}$ (5) $\displaystyle|\psi(t)\rangle_{j}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2}}(|1\rangle_{j}+e^{i\omega t}|-1\rangle_{j}$
(6)

where $\epsilon$ denotes an error rate due to the energy relaxation. Since
every unstable system shows a quadratic decay for a short time scale Schulman
_et al._ (1994); Nakazato  _et al._ (1996), the error can be approximately
described by $1-\epsilon\simeq 1-(t/T^{\prime}_{1})^{2}$, where
$T^{\prime}_{1}$ denotes a energy relaxation time as long as we are interested
in a short time scale. We assume that the experimentalists do not know the
precise value of $T^{\prime}_{1}$. In the experiment, the energy relaxation
time can fluctuate in time Yan  _et al._ (2016); Abdurakhimov  _et al._
(2019). In this case, we could not estimate the exact value of
$T^{\prime}_{1}$ because it changes before a precise estimation of
$T^{\prime}_{1}$. In our paper, we assume that, $T^{\prime}_{1}$ is constant
during the interaction between the spin-1 system and magnetic fields. However,
it is worth mentioning that $T^{\prime}_{1}$ could change if we consider a
longer time scale such as an initialize time, a readout time, or a repetition
time (when we have $N\geq 2$), which makes it difficult to estimate the exact
value of the $T_{1}^{\prime}$.

###  III.2 Uncertainty of the estimation for a Gaussian decay

Let us calculate the uncertainty of the estimation when we have a Gaussian
decay such as $\epsilon=1-e^{-(t/T^{\prime}_{1})^{2}}$. We perform projective
measurement with $\hat{\mathcal{P}}_{j}=(\openone+\hat{\sigma}^{(j)}_{y})/2$
for $j=1,2,\cdots L$ on the state described by Eq. (5). Here,
$\hat{\sigma}^{(j)}_{y}$ is the Pauli matrix spanned by $|1\rangle$ and
$|-1\rangle$. The probability to obtain this projection is as follows.

$\displaystyle P^{\prime}\simeq\frac{e^{-(t/T^{\prime}_{1})^{2}}}{2}(1+\omega
t)$ (7)

where we assume $|\omega t|\ll 1$. Since we do not know the precise value of
$T^{\prime}_{1}$, we use $T_{1}(\neq T_{1}^{\prime})$ to estimate the
projection probability. This means that we consider the projection probability
as $P\simeq\frac{e^{-(t/T_{1})^{2}}}{2}(1+\omega t)$, which is different from
the true probability $P^{\prime}$. In this case, the imperfect knowledge of
the energy relaxation time induce the systematic errors. We can calculate the
estimation uncertainty from Eq. (3) as follows.

$\displaystyle\delta^{2}\omega$ $\displaystyle\simeq$
$\displaystyle\frac{e^{2(t/T_{1})^{2}}}{t^{2}}\Bigl{[}\frac{e^{-(t/T^{\prime}_{1})^{2}}(2-e^{-(t/T^{\prime}_{1})^{2}})}{L}\Bigr{.}$
(8) $\displaystyle+$
$\displaystyle\Bigl{.}(e^{-(t/T_{1})^{2}}-e^{-(t/T^{\prime}_{1})^{2}})^{2}\Bigr{]},$

where the first term comes from the statistical error while the second term
comes from the systematic errors due to the lack of the knowledge of the
coherence time.

Here, we try to minimize the uncertainty by optimizing the interaction time
$t$. In the conventional strategy, since we believe that $T_{1}\simeq
T_{1}^{\prime}$, we obtain the following ’wrong’ uncertainty by substituting
$T_{1}^{\prime}=T_{1}$ in Eq. (8):

$\displaystyle\delta^{2}\omega_{\rm{est}}=\frac{2e^{(t/T_{1})^{2}}-1}{(t/T_{1})^{2}LT_{1}^{2}}.$
(9)

Actually, we can minimize $\delta^{2}\omega_{\rm{est}}$ with $t=cT_{1}$ and
$c\simeq 0.876$. However, since the actual coherence time $T_{1}^{\prime}$ is
different from $T_{1}$, the optimal $t$ obtained above induces the bias. To
check this, we calculate the ’actual’ uncertainty under $T_{1}^{\prime}\neq
T_{1}$ for $t=cT_{1}$ using Eq. (8) and obtain

$\displaystyle\delta^{2}\omega_{\rm{conv}}$ $\displaystyle\simeq$
$\displaystyle\frac{e^{2c^{2}}}{t^{2}}\Bigl{[}\frac{e^{-(cT_{1}/T_{1}^{\prime})^{2}}(2-e^{-(cT_{1}/T_{1}^{\prime})^{2}})}{L}\Bigr{.}$
(10) $\displaystyle+$
$\displaystyle\Bigl{.}(e^{-c^{2}}-e^{-(cT_{1}/T_{1}^{\prime})^{2}})^{2}\Bigr{]}.$

We can confirm that even in the limit of large $L$, the uncertainty approaches
to
$\lim_{L\rightarrow\infty}\delta^{2}\omega_{\rm{conv}}=(e^{-c^{2}}-e^{-(cT_{1}/T_{1}^{\prime})^{2}})^{2})$,
and so the uncertainty is lower bounded due to the systematic errors. This
could be a serious problem to use the quantum sensor for practical purposes.

Figure 1: The uncertainty of the estimation against the number of qubits with
systematic errors. We set $T_{1}=1$ $s$, $T_{1}^{\prime}=1.4$ $s$, and
$t=2L^{-1/4}$ $s$. The estimation uncertainty of our scheme becomes smaller
than that of the conventional scheme for $L\geqq 26$.

On the other hand, in our strategy, we set the interaction time $t$ in the
Zeno regime to suppress the systematic errors. We consider a case with
$t/T^{\prime}_{1}\ll 1$. To satisfy this condition, we remark that it is not
necessary to know the exact value of $T_{1}^{\prime}$, but just necessary to
know the order of $T_{1}^{\prime}$. In this case, we can perform a Taylor
expansion on the estimation uncertainty for $t/T^{\prime}_{1}\ll 1$, and we
obtain the following.

$\displaystyle\delta^{2}\omega$ $\displaystyle\simeq$
$\displaystyle\frac{1}{Lt^{2}}+\frac{1}{t^{2}}\left[\left(\frac{t}{T_{1}}\right)^{2}-\left(\frac{t}{T^{\prime}_{1}}\right)^{2}\right]^{2}$
(11) $\displaystyle\geq$
$\displaystyle\frac{2}{\sqrt{L}}\left|\left(\frac{1}{T_{1}}\right)^{2}-\left(\frac{1}{T^{\prime}_{1}}\right)^{2}\right|$

where we use an inequality of arithmetic and geometric means. The interaction
time that minimize $\delta^{2}\omega$ is given as follows:

$\displaystyle\frac{1}{Lt^{2}}$ $\displaystyle=$
$\displaystyle\frac{1}{t^{2}}\left[\left(\frac{t}{T_{1}}\right)^{2}-\left(\frac{t}{T^{\prime}_{1}}\right)^{2}\right]^{2}$
$\displaystyle\Leftrightarrow t$ $\displaystyle=$ $\displaystyle
t_{\rm{opt}}=L^{-1/4}\left|\left(\frac{1}{T_{1}}\right)^{2}-\left(\frac{1}{T^{\prime}_{1}}\right)^{2}\right|^{-1/2}.$
(12)

It is noteworthy that, since we do not know the precise value of
$T^{\prime}_{1}$, we cannot choose the optimal time
$t_{\rm{opt}}=L^{-1/4}|(1/T_{1})^{2}-(1/T^{\prime}_{1})^{2}|^{-1/2}$. Instead,
we can take $t=\tau L^{-1/4}$ where $\tau$ can be determined by a rough
estimation of $T_{1}^{\prime}$. Therefore, by taking $t=\Theta(L^{-1/4})$, we
obtain $\delta^{2}\omega=\Theta(L^{-1/2})$ regardless of the choice of $\tau$.
This means that, by increasing the number of qubits, we can decrease the
estimation uncertainty as small as we want, which is in stark contrast to the
conventional scheme to have a lower bound of the estimation uncertainty. Note
that, although such measurements in the Zeno regime were utilized in
entanglement enhanced sensing Matsuzaki  _et al._ (2011); Chin  _et al._
(2012); Macieszczak (2015), we firstly utilize this concept to suppress the
systematic errors of the quantum sensing.

We plot the uncertainty of the estimation with the systematic errors to
compare the uncertainty of our scheme with that of the conventional scheme, as
shown in Fig. 1. The uncertainty of the conventional scheme approaches to a
finite non-zero value, while the uncertainty of our scheme decreases as we
increase the number of qubits. In our scheme, the statistical error is larger
than that of the conventional scheme due to the short interaction time between
the qubits and magnetic fields. So, for the small number of the qubits where
the statistical error is more dominant than the systematic errors, the
uncertainty of our scheme is larger than that of the conventional scheme.
However, as we increase the number of qubits, the systematic errors become
more dominant than the statistical error, and our scheme shows a better
performance than the conventional scheme.

###  III.3 Uncertainty of the estimation for a realistic decay

Here, we calculate the estimation uncertainty for a more realistic decay. We
consider a spin-boson model with a Lorentzian form factor, and this is a
typical noise model when the spin-1 system is coupled with a leaky cavity
Koshino and Shimizu (2005). (See the details in Appendix A). The error rate
can be described as follows:

$\displaystyle\epsilon$ $\displaystyle=$ $\displaystyle 1-|f(t)|^{2}$
$\displaystyle f(t)$ $\displaystyle=$
$\displaystyle\frac{1+\sqrt{1-2\gamma/\Delta}}{2e^{i\lambda_{1}t}\sqrt{1-2\gamma/\Delta}}-\frac{1-\sqrt{1-2\gamma/\Delta}}{2e^{i\lambda_{2}t}\sqrt{1-2\gamma/\Delta}}\
\ \ \ $ (13)

where $\gamma$ denotes the strength of the noise, and $\Delta$ denotes the
linewidth of the form factor. Also, $\lambda_{1}$ and $\lambda_{2}$ are the
solutions of the following equation:

$\displaystyle\lambda(\lambda+i\Delta)-\gamma\Delta/2=0.$ (14)

We can find the coherence time of the system as
$\tilde{T}^{\prime}_{1}=\sqrt{2/(\gamma\Delta)}$ by approximating the error
rate as $\epsilon\simeq 1-\frac{\gamma\Delta t^{2}}{2}$. The $f(t)$ depends on
$\gamma$ and $\tilde{T}^{\prime}_{1}$. We assume that, while we know the value
of $\gamma$, we do not know the exact value of $\tilde{T}^{\prime}_{1}$, and
we have $\tilde{T}_{1}$ as an expected coherence time where
$\tilde{T}^{\prime}_{1}\neq\tilde{T}_{1}$. By substituting
$\Delta=\frac{2}{\gamma(\tilde{T_{1}})^{2}}$ with Eq. (13), we obtain
$f_{\rm{est}}(t)$ that we expect, while the true value is $f_{\rm{true}}(t)$
where $\Delta=\frac{2}{\gamma(\tilde{T^{\prime}_{1}})^{2}}$ is substituted
with Eq. (13).

Let us calculate the uncertainty of the estimation when we have the decay such
as $\epsilon=1-|f(t)|^{2}$. We perform projective measurement with
$\hat{\mathcal{P}}_{j}=(\openone+\hat{\sigma}^{(j)}_{y})/2$ for $j=1,2,\cdots
L$ on the state described by Eq. (5). The probability to obtain this
projection is as follows:

$\displaystyle P\simeq\frac{1-|f(t)|^{2}}{2}(1+\omega t),$ (15)

where we assume $|\omega t|\ll 1$. If the true value of the coherence time is
given, we have $P^{\prime}\simeq\frac{1-|f_{\rm{true}}(t)|^{2}}{2}(1+\omega
t)$, while we have $P\simeq\frac{1-|f_{\rm{est}}(t)|^{2}}{2}(1+\omega t)$ when
we have $f_{\rm{est}}(t)$. By using Eq. (3), we obtain

$\displaystyle\delta^{2}\omega$ $\displaystyle=$
$\displaystyle\frac{1}{t^{2}(1-|f_{\rm{est}}(t)|^{2})^{2}}\Bigl{[}\frac{(1-|f_{\rm{true}}(t)|^{2})(1+|f_{\rm{true}}(t)|^{2})}{L}\Bigr{.}$
(16) $\displaystyle+$
$\displaystyle\Bigl{.}(|f_{\rm{est}}(t)|^{2}-|f_{\rm{true}}(t)|^{2})^{2}\Bigr{]}$

As long as we are interested in a short time region, we can obtain an
approximate form of the error rate such as $1-|f_{\rm{true}}(t)|^{2}\simeq
1-t^{2}/(\tilde{T}^{\prime}_{1})^{2}\simeq
e^{-(t/\tilde{T}^{\prime}_{1})^{2}}$ and $1-|f_{\rm{est}}(t)|^{2}\simeq
e^{-(t/\tilde{T}_{1})^{2}}$. In this case, the uncertainty of Eq. (16) has the
same form as Eq. (8), and we can adopt the same strategy as the case of the
Gaussian decay. Therefore, we can take $t=\Theta(L^{-1/4})$ to obtain the
scaling of $\delta^{2}\omega=\Theta(L^{-1/2})$.

Figure 2: The uncertainty of the estimation against the number of qubits with
systematic errors. We set $\tilde{T}_{1}=2$ $s$, $\tilde{T}^{\prime}_{1}=2.4$
$s$, $\gamma=1$ $s^{-1}$, and $t=2L^{-1/4}$ $s$. The estimation uncertainty of
our scheme becomes smaller than that of the conventional scheme for $L\geqq
279$.

To quantify the performance of our scheme, we plot the uncertainty of the
estimation of our scheme, and also compare it with that of the conventional
scheme in Fig. 2. In the conventional scheme, similar to the Gaussian decay
case, we substitute $\tilde{T}^{\prime}_{1}=\tilde{T}_{1}$ in Eq. (16), and
obtain $\delta^{2}\omega_{\rm{conv}}$. After we fix the values of $\gamma$,
$\tilde{T}_{1}$, and $\tilde{T}^{\prime}_{1}$, we minimize
$\delta^{2}\omega_{\rm{conv}}$ by choosing the optimal $t$, and use this $t$
for the plot. On the other hand, in our scheme, we take $t=\tau L^{-1/4}$.
Similar to the case of the Gaussian decay, the statistical error is relevant
for the small number of the qubits, and the uncertainty of our scheme is
larger than that of the conventional scheme. However, the uncertainty of the
estimation of the conventional scheme has a lower bound even in the large $L$
due to the systematic errors. On the other hand, in our scheme, the
uncertainty decreases as we increase the number of qubits. Therefore, for the
large number of the qubits, our scheme shows better performance than the
conventional one.

##  IV Conclusion

In conclusion, we propose a scheme to use measurements in the Zeno regime for
the suppression of systematic errors in quantum magnetic field sensing. When
we do not know a precise value of the coherence time of the quantum system,
such a lack of knowledge induces the systematic errors for quantum sensing. In
this case, the uncertainty of the estimation of the target magnetic field is
lower bounded even in the limit of the infinite number of the repetitions of
the measurements. To suppress such systematic errors, we utilized the fact
that every decoherence process shows a quadratic decay in a short time region
in the Zeno regime where the effect of the decoherence does not significantly
affect the dynamics. We show that, by taking the interaction time (between the
system and magnetic fields) as $\Theta(L^{-1/4})$, we obtain a scaling of
$\delta^{2}\omega=\Theta(L^{-1/2})$ where $L$ is the number of measurements.
Therefore, by increasing the number of the repetition, we can decrease the
uncertainty as small as we want even under the effect of the systematic errors
due to the imperfect knowledge of the coherence time. Our results pave the way
for the practically useful quantum magnetic fields sensors.

###### Acknowledgements.

We are grateful to Shiro Kawabata for useful discussions. This work was
supported by Leading Initiative for Excellent Young Researchers MEXT Japan and
JST presto (Grant No. JPMJPR1919) Japan.

##  Appendix A Decay dynamics under the effect of the decoherence

Here, we discuss the decay dynamics of the three level system under the effect
of the decoherence. Especially, we consider an NV center for the three level
system.

The Hamiltonian of the NV center coupled with an environment is described as
follows.

$\displaystyle H$ $\displaystyle=$ $\displaystyle H_{nv}+H_{I}+H_{E}$
$\displaystyle H_{nv}$ $\displaystyle=$ $\displaystyle
D_{0}\hat{S}_{z}^{2}+g\mu_{b}B\hat{S}_{z}$ $\displaystyle H_{I}$
$\displaystyle=$
$\displaystyle\hat{S}_{x}\int_{\mu}g_{\mu}(\hat{b}_{\mu}+\hat{b}_{\mu}^{\dagger})d\mu$
$\displaystyle H_{E}$ $\displaystyle=$
$\displaystyle\int_{\mu}\omega_{\mu}\hat{b}^{\dagger}_{\mu}\hat{b}_{\mu}d\mu$
(17)

where $D_{0}$ denotes a zero field splitting, $g$ denotes a g factor,
$\mu_{b}$ denotes a Bohr magneton, $g_{\mu}$ denotes the coupling strength
between the NV center and environment, $\omega_{\mu}$ denotes the frequency of
the environment, $\hat{S}_{x}$ ($\hat{S}_{y}$) denote the spin $1$ operator
for $x$ ($y$), and $\hat{b}_{\mu}$ denotes an annihilation (creation) operator
of the environmental mode. In the interaction picture, we have

$\displaystyle H_{I}(t)=\int_{\mu}g_{\mu}(|B\rangle\langle
0|\hat{b}_{\mu}e^{-i\Delta\omega_{\mu}}+|0\rangle\langle
B|(\hat{b}_{\mu}^{\dagger})e^{i\Delta\omega_{\mu}})d\mu$

where $\Delta\omega_{\mu}=\omega_{\mu}-D_{0}-g\mu_{b}B$ and
$|B\rangle=\frac{1}{\sqrt{2}}(|1\rangle+|-1\rangle)$. Here, we can consider a
subspace spanned by $|B\rangle$ and $|0\rangle$ because the dark state
$|D\rangle=\frac{1}{\sqrt{2}}(|1\rangle-|-1\rangle)$ is not involved in the
dynamics, In this case, we can adopt the standard results of the open quantum
system for a qubit. We consider the following Lorentzian form factor Koshino
and Shimizu (2005).

$\displaystyle|g_{\mu}|^{2}=\frac{\gamma}{2\pi}\frac{\Delta^{2}}{(\mu-\mu_{0})^{2}+\Delta^{2}}$
(18)

where $\gamma$ denotes the strength of the noise, $\mu_{0}$ denotes a central
frequency of the environment, and $\Delta$ denotes the linewidth of the form
factor. In our paper, we set $\mu_{0}=D_{0}\simeq D_{0}+g\mu_{b}B$ where
$g\mu_{b}B$ is assumed to be much smaller than any other parameters. We
consider an initial state as $|B\rangle$, and we can calculate the survival
probability as follows Koshino and Shimizu (2005):

$\displaystyle P_{\rm{suv}}$ $\displaystyle=$ $\displaystyle|f(t)|^{2}$
$\displaystyle f(t)$ $\displaystyle=$
$\displaystyle\frac{1+\sqrt{1-2\gamma/\Delta}}{2\sqrt{1-2\gamma/\Delta}}e^{-i\lambda_{1}t}-\frac{1-\sqrt{1-2\gamma/\Delta}}{2\sqrt{1-2\gamma/\Delta}}e^{-i\lambda_{2}t}$

Also, $\lambda_{1}$ and $\lambda_{2}$ are the solutions of the following
equation:

$\displaystyle\lambda(\lambda+i\Delta)-\gamma\Delta/2=0.$ (20)

The density matrix in the interaction picture after a time $t$ is described as
follows:

$\displaystyle\rho_{I}(t)=P_{\rm{suv}}|B\rangle\langle
B|+(1-P_{\rm{suv}})|0\rangle\langle 0|.$ (21)

Therefore, by going back to the Schrödinger picture, we have

$\displaystyle\rho(t)=P_{\rm{suv}}|\psi(t)\rangle\langle\psi(t)|+(1-P_{\rm{suv}})|0\rangle\langle
0|,$ (22)

where $|\psi(t)\rangle=\frac{1}{\sqrt{2}}(|+1\rangle+e^{i\omega t}|-1\rangle)$
and $\omega=2g\mu_{b}B$. We can thus define an error rate as
$\epsilon=1-P_{\rm{suv}}$.

## References

  * Degen  _et al._ (2017) C. L. Degen, F. Reinhard, and P. Cappellaro, Reviews of modern physics 89, 035002 (2017).
  * Taylor  _et al._ (2008) J. Taylor, P. Cappellaro, L. Childress, L. Jiang, D. Budker, P. Hemmer, A. Yacoby, R. Walsworth, and M. Lukin, Nature Physics 4, 810 (2008).
  * Maze  _et al._ (2008) J. R. Maze, P. L. Stanwix, J. S. Hodges, S. Hong, J. M. Taylor, P. Cappellaro, L. Jiang, M. G. Dutt, E. Togan, A. Zibrov, _et al._ , Nature 455, 644 (2008).
  * Balasubramanian  _et al._ (2008) G. Balasubramanian, I. Chan, R. Kolesov, M. Al-Hmoud, J. Tisler, C. Shin, C. Kim, A. Wojcik, P. R. Hemmer, A. Krueger, _et al._ , Nature 455, 648 (2008).
  * Bal  _et al._ (2012) M. Bal, C. Deng, J.-L. Orgiazzi, F. Ong, and A. Lupascu, Nature communications 3, 1 (2012).
  * Toida  _et al._ (2019) H. Toida, Y. Matsuzaki, K. Kakuyanagi, X. Zhu, W. J. Munro, H. Yamaguchi, and S. Saito, Communications Physics 2, 1 (2019).
  * Budoyo  _et al._ (2020) R. P. Budoyo, K. Kakuyanagi, H. Toida, Y. Matsuzaki, and S. Saito, Applied Physics Letters 116, 194001 (2020).
  * Kominis  _et al._ (2003) I. Kominis, T. Kornack, J. Allred, and M. V. Romalis, Nature 422, 596 (2003).
  * Yan  _et al._ (2016) F. Yan, S. Gustavsson, A. Kamal, J. Birenbaum, A. P. Sears, D. Hover, T. J. Gudmundsen, D. Rosenberg, G. Samach, S. Weber, _et al._ , Nature communications 7, 1 (2016).
  * Abdurakhimov  _et al._ (2019) L. V. Abdurakhimov, I. Mahboob, H. Toida, K. Kakuyanagi, and S. Saito, Applied Physics Letters 115, 262601 (2019).
  * Wolf  _et al._ (2015) T. Wolf, P. Neumann, K. Nakamura, H. Sumiya, T. Ohshima, J. Isoya, and J. Wrachtrup, Physical Review X 5, 041001 (2015).
  * Misra and Sudarshan (1977) B. Misra and E. G. Sudarshan, Journal of Mathematical Physics 18, 756 (1977).
  * Facchi and Pascazio (2008) P. Facchi and S. Pascazio, Journal of Physics A: Mathematical and Theoretical 41, 493001 (2008).
  * Schulman  _et al._ (1994) L. Schulman, A. Ranfagni, and D. Mugnai, Physica Scripta 49, 536 (1994).
  * Nakazato  _et al._ (1996) H. Nakazato, M. Namiki, and S. Pascazio, International Journal of Modern Physics B 10, 247 (1996).
  * Erez  _et al._ (2004) N. Erez, Y. Aharonov, B. Reznik, and L. Vaidman, Physical Review A 69, 062315 (2004).
  * Koshino and Shimizu (2005) K. Koshino and A. Shimizu, Physics reports 412, 191 (2005).
  * Itano  _et al._ (1990) W. M. Itano, D. J. Heinzen, J. Bollinger, and D. Wineland, Physical Review A 41, 2295 (1990).
  * Kwiat  _et al._ (1995) P. Kwiat, H. Weinfurter, T. Herzog, A. Zeilinger, and M. A. Kasevich, Physical Review Letters 74, 4763 (1995).
  * Fischer  _et al._ (2001) M. Fischer, B. Gutiérrez-Medina, and M. Raizen, Physical review letters 87, 040402 (2001).
  * Streed  _et al._ (2006) E. W. Streed, J. Mun, M. Boyd, G. K. Campbell, P. Medley, W. Ketterle, and D. E. Pritchard, Physical review letters 97, 260402 (2006).
  * Helmer  _et al._ (2009) F. Helmer, M. Mariantoni, E. Solano, and F. Marquardt, Physical Review A 79, 052115 (2009).
  * Wolters  _et al._ (2013) J. Wolters, M. Strauß, R. S. Schoenfeld, and O. Benson, Physical Review A 88, 020101 (2013).
  * Kakuyanagi  _et al._ (2015) K. Kakuyanagi, T. Baba, Y. Matsuzaki, H. Nakano, S. Saito, and K. Semba, New Journal of Physics 17, 063035 (2015).
  * Kondo  _et al._ (2016) Y. Kondo, Y. Matsuzaki, K. Matsushima, and J. G. Filgueiras, New Journal of Physics 18, 013033 (2016).
  * Kalb  _et al._ (2016) N. Kalb, J. Cramer, D. J. Twitchen, M. Markham, R. Hanson, and T. Taminiau, Nature communications 7, 1 (2016).
  * Matsuzaki  _et al._ (2011) Y. Matsuzaki, S. C. Benjamin, and J. Fitzsimons, Phys. Rev. A 84, 012103 (2011).
  * Chin  _et al._ (2012) A. W. Chin, S. F. Huelga, and M. B. Plenio, Phys. Rev. Lett. 109, 233601 (2012).
  * Macieszczak (2015) K. Macieszczak, Phys. Rev. A 92, 010102 (2015).
  * Sugiyama (2015) T. Sugiyama, Physical Review A 91, 042126 (2015).
  * Takeuchi  _et al._ (2019) Y. Takeuchi, Y. Matsuzaki, K. Miyanishi, T. Sugiyama, and W. J. Munro, Physical Review A 99, 022325 (2019).
  * You  _et al._ (2007) J. You, X. Hu, S. Ashhab, and F. Nori, Physical Review B 75, 140515 (2007).


# Core States of Neutron Stars from Anatomizing their Scaled Structure
Equations

Bao-Jun<EMAIL_ADDRESS>Quantum Machine Learning Laboratory, Shadow
Creator Inc., Shanghai 201208, China  Bao-An<EMAIL_ADDRESS>Department of Physics and Astronomy, Texas A$\&$M University-Commerce,
Commerce, TX 75429-3011, USA  Zhen<EMAIL_ADDRESS>Sino-
French Institute of Nuclear Engineering and Technology, Sun Yat-Sen
University, Zhuhai 519082, China

###### ABSTRACT

Given an Equation of State (EOS) for neutron star (NS) matter, there is a
unique mass-radius sequence characterized by a maximum mass
$M_{\textrm{NS}}^{\max}$ at radius $R_{\max}$. We first show analytically that
the $M_{\textrm{NS}}^{\max}$ and $R_{\max}$ scale linearly with two different
combinations of NS central pressure $P_{\textrm{c}}$ and energy density
$\varepsilon_{\textrm{c}}$ by dissecting perturbatively the dimensionless
Tolman-Oppenheimer-Volkoff (TOV) equations governing NS internal variables.
The scaling relations are then verified via 87 widely used and rather diverse
phenomenological as well as 17 microscopic NS EOSs with/without considering
hadron-quark phase transitions and hyperons by solving numerically the
original TOV equations. The EOS of densest NS matter allowed before it
collapses into a black hole (BH) is then obtained. Using the universal
$M_{\textrm{NS}}^{\max}$ and $R_{\max}$ scalings and NICER (Neutron Star
Interior Composition Explorer) and XMM-Newton mass-radius observational data
for PSR J0740+6620, a very narrow constraining band on the NS central EOS is
extracted directly from the data for the first time without using any specific
input EOS model.

## 1\. Introduction

Understanding the nature and Equation of State (EOS) of supradense matter in
neutron stars (NSs) (Walecka, 1974; Collins & Perry, 1975; Chin, 1977;
Freedman & McLerran, 1977a, b, c; Wiringa et al., 1988; Akmal et al., 1998)
has been a longstanding and shared goal of both astrophysics and nuclear
physics (Shuryak, 1980; Lattimer & Prakash, 2001; Danielewicz et al., 2002;
Alford et al., 2008; Li et al., 2008; Özel & Freire, 2016; Watts et al., 2016;
Oertel et al., 2017; Baym et al., 2018; Vidaña, 2018; Baiotti, 2019; Orsaria
et al., 2019; Dexheimer et al., 2021; Lattimer, 2021; Lovato et al., 2022;
Sorensen et al., 2023). In particular, novel phases and new degrees of freedom
are expected to appear in NS cores where the density is the highest. Finding
the EOS of densest visible matter in the Universe is an ultimate goal of
astrophysics in the era of high-precision multi-messenger astronomy
(Sathyaprakash et al., 2019). However, despite of much effort using various
data and models over the last few decades information about the NS core EOS
remains elusive and ambiguous (De et al., 2018; Radice et al., 2018; Tews et
al., 2018; Bauswein et al., 2019; Baym et al., 2019; Montaña et al., 2019;
Most et al., 2019; Xie & Li, 2019, 2021; Zhang & Li, 2019, 2019, 2021;
Drischler et al., 2020; Weih et al., 2020; Zhao & Lattimer, 2020; Miao et al.,
2021; Pang et al., 2021; Raaijmakers et al., 2021; Breschi et al., 2022; Ecker
& Rezzolla, 2022; Huth et al., 2022; Tan et al., 2022a, b; Perego et al.,
2022) partially because of the strong degeneracy between the core EOSs and the
still poorly understood crust. Can one extract directly the core EOS from some
NS observables without using any specific model EOS?

In this work, we show first analytically and then numerically that the maximum
mass $M_{\textrm{NS}}^{\max}$ and the corresponding radius $R_{\max}$ on the
mass-radius sequence for any EOS scales universally with
$\Gamma_{\textrm{c}}\equiv\varepsilon_{\textrm{c}}^{-1/2}[\widehat{P}_{\textrm{c}}/(1+3\widehat{P}_{\textrm{c}}^{2}+4\widehat{P}_{\textrm{c}})]^{3/2}$
and
$\nu_{\textrm{c}}\equiv\varepsilon_{\textrm{c}}^{-1/2}[\widehat{P}_{\textrm{c}}/(1+3\widehat{P}_{\textrm{c}}^{2}+4\widehat{P}_{\textrm{c}})]^{1/2}$,
respectively, where $\widehat{P}_{\textrm{c}}$ is the ratio of pressure
$P_{\textrm{c}}$ over energy density $\varepsilon_{\textrm{c}}$ at NS centers.
The correlation between $M_{\textrm{NS}}^{\max}$ and $\Gamma_{\textrm{c}}$, in
particular, gives the EOS of densest matter allowed in NSs before they
collapse into black holes (BHs). Applying these universal scalings to
NICER/XMM-Newton’s mass-radius observational data (abbreviated as NICER since
then) for PSR J0740+6620 (Fonseca et al., 2021; Riley et al., 2021; Miller et
al., 2021; Salmi et al., 2022), a very narrow constraining band on the NS core
EOS is directly extracted without using any model EOS.

## 2\. NS Radius/Mass Scalings by Perturbatively Solving the Scaled Structural
Equations

Information about the EOS in the entire NS is normally obtained by comparing
observations with predictions of a given EOS model by integrating the Tolman-
Oppenheimer-Volkoff (TOV) equations (Tolman, 1939; Oppenheimer & Volkoff,
1939) from the NS center to its surface. In the following, we present a novel
approach to access directly the NS core EOS without the hinderance of its
still unconstrained EOSs in other regions especially the crust.

Using the energy density scale $W=G^{-1}(4\pi
G\varepsilon_{\textrm{c}})^{-1/2}$ and length scale $Q=(4\pi
G\varepsilon_{\textrm{c}})^{-1/2}$, the reduced mass $\widehat{M}\equiv M/W$,
radius $\widehat{r}\equiv r/Q$, pressure $\widehat{P}\equiv
P/\varepsilon_{\textrm{c}}$ and energy density
$\widehat{\varepsilon}\equiv\varepsilon/\varepsilon_{\textrm{c}}$ (with
$\varepsilon_{\textrm{c}}$ being the central energy density) satisfy the
dimensionless TOV equations (adopting $c=1$),

$\frac{\mathrm{d}}{\mathrm{d}\widehat{r}}\widehat{P}=-\frac{(\widehat{P}+\widehat{\varepsilon})(\widehat{r}^{3}\widehat{P}+\widehat{M})}{\widehat{r}^{2}-2\widehat{M}\widehat{r}},~{}~{}\frac{\mathrm{d}}{\mathrm{d}\widehat{r}}\widehat{M}=\widehat{r}^{2}\widehat{\varepsilon}.$
(1)

The smallness of $\widehat{P}$ enables solving perturbatively the above
coupled evolution equations of $\widehat{P}$ and $\widehat{M}$ given an EOS
$P(\varepsilon)$. Starting from the NS center at $\widehat{r}=0$, the
(reduced) energy density $\widehat{\varepsilon}$, pressure $\widehat{P}$ and
mass $\widehat{M}$ can be expanded as $\widehat{\varepsilon}\approx
1+a_{1}\widehat{r}+a_{2}\widehat{r}^{2}+\cdots$,
$\widehat{P}\approx\widehat{P}_{\textrm{c}}+b_{1}\widehat{r}+b_{2}\widehat{r}^{2}+\cdots$,
and $\widehat{M}\approx c_{1}\widehat{r}+c_{2}\widehat{r}^{2}+\cdots$,
respectively. Moreover, the reduced EOS $\widehat{P}(\widehat{\varepsilon})$
can be written as $\widehat{P}=\sum_{k=1}d_{k}\widehat{\varepsilon}^{k}$ with
coefficients $\\{d_{k}\\}$. The inner-relations among the coefficients
$\\{a_{k}\\}$, $\\{b_{k}\\}$ and $\\{c_{k}\\}$ are determined by the NS EOS
and the TOV equations. Particularly, as shown in the APPENDIX, all $a_{k}$’s
and $b_{k}$’s with odd $k$ are zero and the first nonzero coefficient is
$b_{2}=-6^{-1}(1+3\widehat{P}_{\textrm{c}}^{2}+4\widehat{P}_{\textrm{c}})$,
$c_{1}=c_{2}=0$ and $c_{k}=a_{k-3}/k$ for $k\geq 3$.

Within the perturbative approach outlined above, the reduced radius
$\widehat{R}$ defined via the vanishing of pressure $\widehat{P}$ is
determined to order $b_{2}$ by
$\widehat{P}_{\textrm{c}}+b_{2}\widehat{R}^{2}=0$. The resulting
$\widehat{R}=[6\widehat{P}_{\textrm{c}}/(1+3\widehat{P}_{\textrm{c}}^{2}+4\widehat{P}_{\textrm{c}})]^{1/2}$
multiplied by $Q\sim 1/\sqrt{\varepsilon_{\textrm{c}}}$ makes the physical
radius $R=Q\widehat{R}$ to scale as

$R\sim\nu_{\textrm{c}}\equiv\frac{\widehat{P}_{\textrm{c}}^{1/2}}{\sqrt{\varepsilon_{\textrm{c}}}}\cdot\left(\frac{1}{1+3\widehat{P}_{\textrm{c}}^{2}+4\widehat{P}_{\textrm{c}}}\right)^{1/2}.$
(2)

Since the reduced NS mass is $\widehat{M}(\widehat{R})\approx
c_{3}\widehat{R}^{3}=\widehat{R}^{3}/3$, its physical mass
$M_{\textrm{NS}}\equiv M\sim\widehat{R}^{3}/\sqrt{\varepsilon_{\textrm{c}}}$
then scales as

$M_{\textrm{NS}}\sim\Gamma_{\textrm{c}}\equiv\frac{\widehat{P}_{\textrm{c}}^{3/2}}{\sqrt{\varepsilon_{\textrm{c}}}}\cdot\left(\frac{1}{1+3\widehat{P}_{\textrm{c}}^{2}+4\widehat{P}_{\textrm{c}}}\right)^{3/2}.$
(3)

The above scalings relate directly the radius and mass observables of a NS
with its central EOS
${P}_{\textrm{c}}(\varepsilon_{\textrm{c}})=P(\varepsilon_{\textrm{c}})$.
These scalings are general since they are direct consequences of the TOV
equations themselves without assuming any particular structure, composition
and EOS for the NS. Once the $\nu_{\textrm{c}}$ or $\Gamma_{\textrm{c}}$ is
constrained within a certain range by observational data, it can be used to
determine the NS central EOS. Moreover, the scaling for NS mass is expected to
be more accurate than that for the radius which is affected more by the
uncertain low-density EOS for the NS crust (Baym et al., 1971), see FIG. 1.
Actually the low-density EOS may affect the crust thickness and consequently
the NS radius sizably while has little influence on its mass (Xu et al.,
2009).

Extensive studies on scaling relations among NS observables exist in the
literature, see, e.g., Refs. Rhoades & Ruffini (1974); Hartle (1978); Kalogera
& Baym (1996); Lattimer & Prakash (2005); Tsui & Leung (2005); Yagi & Yunes
(2013); Wen et al. (2019); Jiang & Yagi (2020) and references therein. These
EOS-independent universial scalings are mostly for NS bulk properties, such as
the moment of inertia ($I$), tidal Love number ($L$), (spin-induced)
quadrupole moment (${Q}$), compactness ($C$), and frequencies ($f$ or
$\omega$) of various oscillation modes of NSs. Knowing one or more
observables, these scalings (e.g., $I$-$L$-${Q}$, $I$-$L$-$C$ or
$M\omega$-$C$) enable the finding of other observables or bulk properties that
have not been or hard to be measured. They are completely different from the
EOS-independent scalings between the radius and mass of most massive NSs and
the combinations of the reduced core pressure and energy density that we have
found and used here for the first time. Among the earlier scalings closest to
ours is the one showing $M_{\textrm{NS}}^{\max}\sim
D_{\textrm{M}}\varepsilon_{\textrm{c}}^{-1/2}$ (see, e.g., Refs. Hartle
(1978); Kalogera & Baym (1996); Lattimer & Prakash (2005)) with the
coefficient $D_{\textrm{M}}$ estimated using either unrealistic or certain
selected microscopic dense matter EOSs. In addition, a numerical power law was
given as $R\sim D_{\textrm{R}}P^{\delta}$ with $\delta\approx
0.23\mbox{$\sim$}0.26$ (Lattimer & Prakash, 2001), where $P$ is the pressure
at about 1-2 times the nuclear saturation density $\rho_{\textrm{sat}}$. Since
the coefficients $D_{\textrm{M}}$ and $D_{\textrm{R}}$ are obtained quite
empirically, the corresponding predictions either have sizable EOS model
dependence or lack solid/clear physical origins. Considering our method and
most importantly, it is the scaling of the observables
$M_{\textrm{NS}}^{\max}$ and radius $R_{\max}$ with the internal variables
(different combinations of the core pressure and energy density via
$\nu_{\textrm{c}}$ and $\Gamma_{\textrm{c}}$) that enable us to extract the
relationship (core EOS)
${P}_{\textrm{c}}(\varepsilon_{\textrm{c}})=P(\varepsilon_{\textrm{c}})$ once
either one of the observables is constrained observationally, instead of just
the individual value of $\varepsilon_{\textrm{c}}$ or the pressure around
$(1\mbox{$\sim$}2)\rho_{\textrm{sat}}$ in the previous studies of scalings.
Thus, as outlined in the previous section and discussed in detail in the
APPENDIX, our method and results are completely novel with respect to the
exisiting work in the literature.

## 3\. Validations of NS Radius/mass Scalings by Solving the Original
Structural Equations

The maximum mass $M_{\textrm{NS}}^{\max}$ and radius $R_{\max}$ characterize
the densest stable NS configuration on the mass-radius curve before it
collapses to a BH for a given NS EOS. The $M_{\textrm{NS}}^{\max}$ is required
to be at least as large as the currently observed most massive NS. Using the
signals from GW170817 (Abbott et al., 2018), the $M_{\textrm{NS}}^{\max}$
(also known as the $M_{\textrm{TOV}}$ for non-rotating NSs) is currently
predicted to be about $2.01$$\sim$$2.16$ solar mass $M_{\odot}$ (Rezzolla et
al., 2018), similarly Margalit & Metzger (2017) showed it should be $\lesssim
2.17M_{\odot}$ by including the electromagnetic information). In addition, a
consistent prediction of $M_{\textrm{NS}}^{\max}\lesssim
2.16\mbox{$\sim$}2.28M_{\odot}$ through general relativistic
magnetohydrodynamics simulations was also given (Ruiz et al., 2018). We notice
that these predictions are all somewhat EOS-model dependent.

To test our mass and radius scalings, shown in FIG. 1 are the
$R_{\max}$-$\nu_{\textrm{c}}$ (upper panel) and
$M_{\textrm{NS}}^{\max}$-$\Gamma_{\textrm{c}}$ (lower panel) correlations by
using 87 phenomenological and 17 extra microscopic NS EOSs with and/or without
considering hadron-quark phase transitions and hyperons by solving numerically
the original TOV equations. The phenomenological EOSs are from three rather
different classes of EOS models (Iida & Sato, 1997): the covariant
relativistic mean-field (RMF) models (Serot & Walecka, 1986), the non-
relativistic energy density functionals (EDF) of the Gogny-like type (Cai &
Li, 2022) as well as the conventional Skyrme-type (Vautherin & Brink, 1972;
Stone & Reinhard, 2007; Zhang & Chen, 2016). Within each EOS class, different
EOS models are generated by perturbing the model coupling constants within
their empirical uncertain ranges. In addition, shown with the lightblue “+”
symbols are results from using 17 extra microscopic EOS models. These include
(1) the APR EOS (Akmal et al., 1998) based on microscopic nuclear many-body
theories using 2-body and 3-body nucleon forces; (2) the Dirac-Brueckner-
Hartree-Fock method (sets MPA1 (Müther et al., 1987) and ENG (Engvik et al.,
1996)); (3) the quark mean field model (Li et al., 2020); (4) the RMF
calculations with hyperons (Glendenning, 1985) (as well as sets DD2 (Typel et
al., 2010) and H4 (Lackey et al., 2006)); (5) the chiral mean field model
with/without hyperons (Dexheimer & Schramm, 2008); (6) the hybrid ALF2 EOS
(Alford et al., 2005) with mixtures of nucleonic and quark matter, and (7) a
few quark matter EOSs (Prakash et al., 1995; Chu et al., 2019). The observed
strong linear correlations demonstrate vividly that the
$R_{\max}$-$\nu_{\textrm{c}}$ and
$M_{\textrm{NS}}^{\max}$-$\Gamma_{\textrm{c}}$ scalings are nearly universal.
It is particularly interesting to notice that EOSs allowing phase transitions
and/or hyperon formations predict consistently the same scalings. By
performing linear fits of the results obtained from the EOS samples, the
scaling relations are determined to be $R_{\textrm{max}}/\textrm{km}\approx
1.05_{-0.03}^{+0.03}\times 10^{3}\nu_{\textrm{c}}+0.64_{-0.29}^{+0.29}$ with
its Pearson’s coefficient about 0.958 and
$M_{\textrm{NS}}^{\textrm{max}}/M_{\odot}\approx 0.173_{-0.003}^{+0.003}\times
10^{4}\Gamma_{\textrm{c}}-0.106_{-0.035}^{+0.035}$ with its Pearson’s
coefficient about 0.986, respectively, here $\nu_{\textrm{c}}$ and
$\Gamma_{\textrm{c}}$ are measured in $\textrm{fm}^{3/2}/\textrm{MeV}^{1/2}$
(as $\varepsilon_{\textrm{c}}$ is measured in $\textrm{MeV}/\textrm{fm}^{3}$).
Moreover, the standard errors for the radius and mass fittings are about 0.39
km and 0.05$\,M_{\odot}$ respectively for the EOS samples considered. In FIG.
1, the condition $M^{\max}_{\textrm{NS}}\gtrsim 1.2M_{\odot}$ used is
necessary to mitigate influences of uncertainties in modeling the crust EOS
(Baym et al., 1971) for low-mass NSs. For the heavier NSs studied here, it is
reassuring to see that although the above 104 EOSs predicted quite different
crust properties, they all fall closely around the same scaling lines
consistently, especially for the
$M_{\textrm{NS}}^{\max}$-$\Gamma_{\textrm{c}}$ relation.

FIG. 1:  The $R_{\max}$-$\nu_{\textrm{c}}$ (upper panel) and
$M_{\textrm{NS}}^{\max}$-$\Gamma_{\textrm{c}}$ (lower panel) correlations with
104 EOS samples (colored symbols), the constraints on the mass (Fonseca et
al., 2021) and radius (Riley et al., 2021) of PSR J0740+6620 are shown by the
pink hatched bands. The two orange arrows and captions nearby in each panel
indicate the corresponding ranges for $\nu_{\textrm{c}}$ and
$\Gamma_{\textrm{c}}$ defined in Eq. (2) and Eq. (3), respectively. The
scaling relations for $R_{\max}/\textrm{km}$ and
$M_{\textrm{NS}}^{\max}/M_{\odot}$ are also shown (with their uncertainties
given in the text).

Using the universal relation between $M_{\textrm{NS}}^{\max}$ and
$\Gamma_{\textrm{c}}$ given and shown in panel (b) of FIG. 1, the EOS of
densest NS matter is found to be
$P_{\textrm{c}}(\varepsilon_{\textrm{c}})\approx
f^{2/3}\varepsilon_{\textrm{c}}^{4/3}\cdot[1+4f^{2/3}\varepsilon_{\textrm{c}}^{1/3}+19f^{4/3}\varepsilon_{\textrm{c}}^{2/3}+\cdots]$,
which depends on the NS mass through
$f\equiv(M_{\textrm{NS}}^{\max}/M_{\odot}+0.106)/1730$. It implies that a
heavier NS has a smaller upper $\varepsilon_{\textrm{c}}$ (as the ratio
$P_{\textrm{c}}/\varepsilon_{\textrm{c}}$ is bounded from above) and therefore
is easier to collapse into a BH (Hawking & Ellis, 1973).

In addition, besides the correlation between $R_{\max}$ and $\nu_{\textrm{c}}$
we also found a strong positive correlation between $R$ (for a given dense
matter EOS) and $\nu_{\textrm{c}}$ along its mass-radius (M-R) curve near the
maximum-mass configuration. Considering now two points on different sides of
the maximum-mass configuration with identical mass, since these two
configurations would have different $\nu_{\textrm{c}}$’s (as the latter
positively correlates with $R$), their $\Gamma_{\textrm{c}}$’s would be
essentially different, implying that the correlation between mass and
$\Gamma_{\textrm{c}}$ starts to break down when going away from the maximum-
mass point (besides the NS crust effects mentioned earlier). Similarly, one
needs to adopt the assumption that the NS is in (near) the maximum-mass
configuration (roughly with $M_{\textrm{NS}}^{\max}\gtrsim 1.2M_{\odot}$) when
the observational mass data is used to infer the central EOS from the mass
scaling (e.g., the constraints shown in FIG. 2). Thus, our results are more
accurate for massive NSs and the constraints on EOS essentially hold (via the
factor $f\equiv(M_{\textrm{NS}}^{\max}/M_{\odot}+0.106)/1730$) if stable NSs
even heavier than PSR J0740+6620 are identified in the future. Moreover, since
the NS maximum-mass configuration on the M-R curve is the closest point to BHs
and some properties of the latter are known to be nearly universal, it is
reasonable to expect that the scalings discussed in this work are most
accurate for most massive NSs.

Extracting precisely the radius is much harder than the mass from observations
of NSs. NICER’s radius measurements for PSR J0740+6620 with its mass $\approx
2.08_{-0.07}^{+0.07}M_{\odot}$ (Fonseca et al., 2021) have been analyzed
independently using different approaches, leading to slightly different
results for its radius. We report here results of our analyses using the
observational radius bands extracted from three studies by the NICER
Collaboration (Riley et al., 2021; Miller et al., 2021; Salmi et al., 2022).
The observational radius band of $12.39_{-0.98}^{+1.30}\,\textrm{km}$ at a 68%
confidence level from Riley et al. (2021) and the mass band given above are
shown (using pink hatched bands) in the upper and lower panels of FIG. 1,
respectively. Their intersections with the two scaling lines indicate that
$10.3\lesssim
10^{3}\nu_{\textrm{c}}/[\textrm{fm}^{3/2}/\textrm{MeV}^{1/2}]\lesssim 12.4$
and
$12.2\lesssim{10^{4}\Gamma_{\textrm{c}}}/[{\textrm{fm}^{3/2}/\textrm{MeV}^{1/2}]}\lesssim
13.1$, respectively, in PSR J0740+6620 (indicated by the orange arrows).
Different observational constraints on the radius may lead to different ranges
for $\nu_{\textrm{c}}$. We found that the radius constraint of
$13.7_{-1.5}^{+2.6}\,\textrm{km}$ (Miller et al., 2021) gives $11.0\lesssim
10^{3}\nu_{\textrm{c}}/[\textrm{fm}^{3/2}/\textrm{MeV}^{1/2}]\lesssim 14.9$.
In addition, the very recent analysis including NICER’s background estimates
(Salmi et al., 2022) found a radius of $12.90_{-0.97}^{+1.25}\,\textrm{km}$,
leading to $10.8\lesssim
10^{3}\nu_{\textrm{c}}/[\textrm{fm}^{3/2}/\textrm{MeV}^{1/2}]\lesssim 12.9$.
The core EOSs extracted using these three analyses of PSR J0740+6620 are
compared in a.Fig. 2 in the APPENDIX. In the following presentation, we focus
on results obtained using the radius band from Riley et al. (2021).

For a comparison, we notice that the hadron-quark hybrid EOS ALF2 predicts its
$M_{\textrm{NS}}^{\max}\approx 2.09M_{\odot}$ (Alford et al., 2005) that is
very close to the mass of PSR J0740+6620. It also predicts
$P_{\textrm{c}}\approx 310\,\textrm{MeV}/\textrm{fm}^{3}$ and
$\varepsilon_{\textrm{c}}\approx 1100\,\textrm{MeV}/\textrm{fm}^{3}$,
therefore $\Gamma_{\textrm{c}}\approx 12.4\times
10^{-4}\,\textrm{MeV}^{1/2}/\textrm{fm}^{3/2}$, which is consistent with the
$\Gamma_{\textrm{c}}$ extracted for PSR J0740+6620 from the mass scaling.

## 4\. Central EOS via NICER’s Mass/radius Measurements and Limits for Speed
of Sound

As discussed earlier, once the $\Gamma_{\textrm{c}}$ or $\nu_{\textrm{c}}$ is
determined within certain range by astrophysical observations, the NS central
EOS $P_{\textrm{c}}(\varepsilon_{\textrm{c}})$ is subsequently constrained.
The $\Gamma_{\textrm{c}}$ and $\nu_{\textrm{c}}$ extracted above from NICER’s
observation of PSR J0740+6620 provide two independent constraints on the NS
central EOS as shown in FIG. 2 by the blue and chocolate bands, respectively.
Since the uncertainty of the radius is currently about 3 times larger than
that of the mass measurement, the resulting constraining band of the central
EOS from the radius measurement is wider. Consequently, the central EOS is
essentially given by the mass measurement via $\Gamma_{\textrm{c}}$ while the
radius effectively sets the upper- and lower-limits for
$\varepsilon_{\textrm{c}}$. In the overlapping area of the two constraining
bands, the most probable NS central EOS can be written as
$P_{\textrm{c}}(\varepsilon_{\textrm{c}})\approx
0.012\varepsilon_{\textrm{c}}^{4/3}\cdot[1+0.047\varepsilon_{\textrm{c}}^{1/3}+0.0026\varepsilon_{\textrm{c}}^{2/3}+\cdots]$
where $\varepsilon_{\textrm{c}}$ varies within about
640$\sim$1210$\,\textrm{MeV}/\textrm{fm}^{3}$.

FIG. 2: Constraints on the central EOS
$P_{\textrm{c}}(\varepsilon_{\textrm{c}})$ independently from the NICER
measurements on the mass of $2.08_{-0.07}^{+0.07}M_{\odot}$ (Fonseca et al.,
2021) (blue band) and the radius of $12.39_{-0.98}^{+1.30}\,\textrm{km}$
(Riley et al., 2021) (chocolate band) for PSR J0740+6620, respectively. The
limit $P_{\textrm{c}}/\varepsilon_{\textrm{c}}\lesssim 0.374$ from causality
is also shown (black dotted), see Eq. (4). The inset plots our central EOS
(magenta crossed band) and predictions/constraints from a few other approaches
(the energy density varies from 500 to 1500$\,\textrm{MeV}/\textrm{fm}^{3}$),
see text for details.

To illustrate some impact of the constraining band on the NS central EOS
inferred above, the inset of FIG. 2 compares our result (shown as magenta
crossed band) with predictions/constraints on the dense NS matter EOS from
several other approaches: (1) (grey shaded) combination (Annala et al., 2018)
of low-density nuclear, high-density perturbative QCD (pQCD) theories and the
tidal deformabilitiies from GW170817 (Abbott et al., 2018); (2) (yellow solid)
pQCD with constraint of the speed of sound at very high densities (Komoltsev &
Kurkela, 2022); (3) (red dash-dotted) similar physical inputs as (1) combining
Gaussian-process (GP) inference algorithm (Gorda et al., 2022); (4) (pink
vertical hatched) consistent construction of EOS ensembles consistent with
nuclear theories, pQCD and astronomical observations and continuous speed of
sound (Altiparmak et al., 2022); (5) (black dashed) nonparametric inference
via GP adopting similar astronomical observations together with (6) (green
horizontal hatched) the prediction on the central pressure-density of PSR
J0740+6620 (Legred et al., 2021); and (7) (blue solid) the constraint from an
earlier model analysis of the NICER data (Miller et al., 2021). Our
constraining band on the central NS EOS directly from the NICER data using the
mass and radius scalings without using any specific input EOS model is mostly
consistent with but much narrower than the NS EOS bands from the previous
studies listed above.

Observational constraints on both $\nu_{\textrm{c}}$ and $\Gamma_{\textrm{c}}$
enable us to determine simultaneously the NS central pressure $P_{\textrm{c}}$
and energy density $\varepsilon_{\textrm{c}}$. More specifically, we obtain
from the relation
$\Gamma_{\textrm{c}}/\nu_{\textrm{c}}=\widehat{P}_{\textrm{c}}/(1+3\widehat{P}_{\textrm{c}}^{2}+4\widehat{P}_{\textrm{c}})$
that $\widehat{P}_{\textrm{c}}\approx 0.24_{-0.07}^{+0.05}$ for PSR
J0740+6620. Equivalently, the most probable values of the central energy
density and pressure are $\varepsilon_{\textrm{c}}\approx
901_{-287}^{+214}\,\textrm{MeV}/\textrm{fm}^{3}$ and $P_{\textrm{c}}\approx
218_{-125}^{+93}\,\textrm{MeV}/\textrm{fm}^{3}$ (black “+” in FIG. 2),
respectively. This value of central energy density is consistent with that
(about 900$\sim$1200 $\textrm{MeV}/\textrm{fm}^{3}$) from a joint Bayesian
inference using NS mass and radius constraints from GW170817 and several other
astronomical observations (Al-Mamun et al., 2021). Using similar observational
data, Brandes et al. (2023) found very recently for a $2.1M_{\odot}$ NS (with
the inferred radius about 11.6 km) that $P_{\textrm{c}}\approx
312_{-169}^{+226}\,\textrm{MeV}/\textrm{fm}^{3}$ and
$\varepsilon_{\textrm{c}}\approx
904_{-327}^{+329}\,\textrm{MeV}/\textrm{fm}^{3}$. Our results are consistent
with theirs but have much smaller uncertainties.

The stiffness of NS EOS can be measured effectively by the speed of sound
squared $s^{2}\equiv\mathrm{d}P/\mathrm{d}\varepsilon$ (Bedaque & Steiner,
2015; Fujimoto et al., 2022). Using the heuristic relation
$M_{\textrm{NS}}\sim\widehat{R}^{3}/\sqrt{\varepsilon_{\textrm{c}}}$ and the
correlation in panel (b) of FIG. 1, we obtain the central speed of sound
squared for the maximum mass configuration,

$s_{\textrm{c}}^{2}=\widehat{P}_{\textrm{c}}\left(1+\frac{1}{3}\frac{1+3\widehat{P}_{\textrm{c}}^{2}+4\widehat{P}_{\textrm{c}}}{1-3\widehat{P}_{\textrm{c}}^{2}}\right).$
(4)

Essentially, $s_{\textrm{c}}^{2}$ grows (minishes) as
$\widehat{P}_{\textrm{c}}$ increases (decreases). For PSR J0740+6620, using
the constraints for $\widehat{P}_{\textrm{c}}$ given above we obtain
$s_{\textrm{c}}^{2}\approx 0.45_{-0.18}^{+0.14}$. The condition
$s_{\textrm{c}}^{2}\leq 1$ leads in addition to an upper bound
$\widehat{P}_{\textrm{c}}\lesssim 0.374$. Therefore although the causality
condition requires apparently $\widehat{P}_{\textrm{c}}\leq 1$, the superdense
nature of core NS matter (indicated by the nonlinear dependence of
$s_{\textrm{c}}^{2}$ on $\widehat{P}_{\textrm{c}}$) renders it to be much
smaller.

Finally, we notice here that slightly different core EOSs are extracted by
using the radii extracted by Miller et al. (2021) and Salmi et al. (2022)
compared to the results presented above using the radius from Riley et al.
(2021). This is what one expects and it verifies the sensitivity of our
approach. Interestingly, as shown numerically in the APPENDIX, the
constraining bands on the core EOSs in PSR J0740+6620 from using the three
radii (Riley et al., 2021; Miller et al., 2021; Salmi et al., 2022) are all
very narrow and they largely overlap. The remaining small widths of the
constraining bands reflect mainly the observational uncertainties as no model
EOS is used in extracting the NS core EOS. As the radius uncertainty becomes
smaller with the ongoing and planned new measurements using more advanced
high-precision observatories (Sathyaprakash et al., 2019; Watts, 2019), our
method can be further applied to refine the core EOS of massive NSs.

## 5\. Conclusion

The mass and radius of the maximum mass configuration of any dense matter EOS
scale linearly with two different combinations of their central pressure and
energy density. The EOS of densest NS matter allowed before it collapses into
a BH is obtained without using any specific model EOS. Using NICER’s mass-
radius observational data for PSR J0740+6620, a very narrow constraining band
on its central EOS is extracted directly from the observational data itself
without using any model EOS for the first time.

Acknowledgements: We thank Peng-Cheng Chu for helpful discussions and
providing us a few quark matter EOSs. This work is supported in part by the
U.S. Department of Energy, Office of Science, under Award No. DE-SC0013702,
the CUSTIPEN (China-U.S. Theory Institute for Physics with Exotic Nuclei)
under US Department of Energy Grant No. DE-SC0009971, the National Natural
Science Foundation of China under Grant No.12235010.

## A. Perturbative Treatments of the Dimensionless TOV Equations

The dimensionless TOV equations are given in Eq. (1), where the reduced
pressure $\widehat{P}$, reduced energy density $\widehat{\varepsilon}$ and the
reduced mass $\widehat{M}$ are all functions of the reduced radius
$\widehat{r}$. Without loosing generality, we expand them in polynomials,

$\displaystyle\widehat{\varepsilon}=$ $\displaystyle
1+\sum_{k=1}a_{k}\widehat{r}^{k}\approx
1+a_{1}\widehat{r}+a_{2}\widehat{r}^{2}+a_{3}\widehat{r}^{3}+\cdots,$ (a1)
$\displaystyle\widehat{P}=$
$\displaystyle\widehat{P}_{\textrm{c}}+\sum_{k=1}b_{k}\widehat{r}^{k}\approx\widehat{P}_{\textrm{c}}+b_{1}\widehat{r}+b_{2}\widehat{r}^{2}+b_{3}\widehat{r}^{3}+\cdots,$
(a2) $\displaystyle\widehat{M}=$
$\displaystyle\sum_{k=1}c_{k}\widehat{r}^{k}\approx
c_{1}\widehat{r}+c_{2}\widehat{r}^{2}+c_{3}\widehat{r}^{3}+\cdots,$ (a3)

here $\widehat{P}_{\textrm{c}}={P}_{\textrm{c}}/{\varepsilon}_{\textrm{c}}$
and
$\widehat{\varepsilon}_{\textrm{c}}=\varepsilon_{\textrm{c}}/\varepsilon_{\textrm{c}}=1$.
By putting these polynomials into the dimensionless TOV equations, and
expanding them to the same order of $\widehat{r}$, we obtain,

$b_{1}=0,$ (a4)

and

$c_{1}=c_{2}=0,~{}~{}c_{3}=\frac{1}{3},~{}~{}c_{k}=\frac{a_{k-3}}{k},~{}~{}\textrm{for}\;k\geq
4.$ (a5)

Moreover,

$\displaystyle b_{2}=$
$\displaystyle-\frac{1}{6}\left(1+3\widehat{P}_{\textrm{c}}^{2}+4\widehat{P}_{\textrm{c}}\right),$
(a6) $\displaystyle b_{3}=$
$\displaystyle-\frac{a_{1}}{36}\left(7+15\widehat{P}_{\textrm{c}}\right),$
(a7) $\displaystyle b_{4}=$
$\displaystyle\left(-\frac{3a_{2}}{10}-\frac{5}{36}+\frac{1}{4}\widehat{P}_{\textrm{c}}^{2}+\frac{1}{3}\widehat{P}_{\textrm{c}}\right)\widehat{P}_{\textrm{c}}-\frac{2a_{2}}{15}+\frac{2}{9}\widehat{P}_{\textrm{c}}-\frac{a_{1}^{2}}{16}.$
(a8)

Next, the pressure could be expanded in terms of the energy density as,

$\widehat{P}=\sum_{k=1}d_{k}\widehat{\varepsilon}^{k}\approx
d_{1}\widehat{\varepsilon}+d_{2}\widehat{\varepsilon}^{2}+d_{3}\widehat{\varepsilon}^{3}+\cdots$
(a9)

which could be rewritten as,

$\displaystyle\widehat{P}\approx$ $\displaystyle
d_{1}+d_{2}+d_{3}+\cdots+a_{1}(d_{1}+2d_{2}+3d_{3}+\cdots)\widehat{r}$
$\displaystyle+\left[(d_{2}+3d_{3}+\cdots)a_{1}^{2}+(d_{1}+2d_{2}+3d_{3}+\cdots)a_{2}\right]\widehat{r}^{2}$
$\displaystyle+\left[\left(a_{1}^{2}d_{3}+2a_{2}d_{2}+6a_{2}d_{3}+\cdots\right)a_{1}+\left(d_{1}+2d_{2}+3d_{3}+\cdots\right)a_{3}\right]\widehat{r}^{3}.$
(a10)

The zeroth-order term gives the sum rule for $\\{d_{k}\\}$ as,

$\widehat{P}_{\textrm{c}}=d_{1}+d_{2}+d_{3}+\cdots=\sum_{k=1}d_{k}.$ (a11)

Moreover, since $b_{1}=0$, we obtain

$a_{1}(d_{1}+2d_{2}+3d_{3}+\cdots)=0.$ (a12)

The speed of sound squared at the center gives another sum rule for
$\\{d_{k}\\}$,

$s_{\textrm{c}}^{2}=\left.\frac{\mathrm{d}\widehat{P}}{\mathrm{d}\widehat{\varepsilon}}\right|_{\widehat{r}=0\leftrightarrow\widehat{\varepsilon}_{\textrm{c}}=1}=d_{1}+2d_{2}+3d_{3}+\cdots=\sum_{k=1}kd_{k},$
(a13)

which is generally nonzero, thus

$a_{1}=0.$ (a14)

Another argument is since $\widehat{P}$ and $\widehat{\varepsilon}$ are both
even functions of $\widehat{r}$, only even terms of $\widehat{r}$ survive, one
then naturally has $b_{1}=0$ and $a_{1}=0$. According to the expression for
$b_{3}\propto a_{1}$, one then obtains $b_{3}=0$ (and immediately $a_{3}=0$
could be obtained considering the coefficients in front of $\widehat{r}^{3}$
of Eq. (A. Perturbative Treatments of the Dimensionless TOV Equations)).

The next-order term is $b_{4}$, i.e.,

$b_{4}=\frac{\widehat{P}_{\textrm{c}}}{12}\left(1+3\widehat{P}_{\textrm{c}}^{2}+4\widehat{P}_{\textrm{c}}\right)-\frac{a_{2}}{30}\left(4+9\widehat{P}_{\textrm{c}}\right).$
(a15)

Here the coefficient $a_{2}$ could be obtained as follows: the speed of sound
squared $s^{2}$ is

$s^{2}=\frac{\mathrm{d}\widehat{P}}{\mathrm{d}\widehat{\varepsilon}}=\frac{\mathrm{d}\widehat{P}}{\mathrm{d}\widehat{r}}\cdot\frac{\mathrm{d}\widehat{r}}{\mathrm{d}\widehat{\varepsilon}}=\frac{2b_{2}\widehat{r}+4b_{4}\widehat{r}^{3}+\cdots}{2a_{2}\widehat{r}+4a_{4}\widehat{r}^{3}+\cdots}=\frac{b_{2}+2b_{4}\widehat{r}^{2}+\cdots}{a_{2}+2a_{4}\widehat{r}^{2}+\cdots}.$
(a16)

Taking it at the center gives $s_{\textrm{c}}^{2}=b_{2}/a_{2}$, see Eq. (4)
for the expression of $s_{\textrm{c}}^{2}$ (for the maximum-mass
configuration). Consequently, $a_{2}=b_{2}/s_{\textrm{c}}^{2}$. The $b_{6}$
coefficient is given by,

$b_{6}=-\frac{1}{216}\left(1+9\widehat{P}_{\textrm{c}}^{2}\right)\left(1+3\widehat{P}_{\textrm{c}}^{2}+4\widehat{P}_{\textrm{c}}\right)-\frac{a_{2}^{2}}{30}+\left(\frac{2}{15}\widehat{P}_{\textrm{c}}^{2}+\frac{1}{45}\widehat{P}_{\textrm{c}}-\frac{1}{54}\right)a_{2}-\frac{5+12\widehat{P}_{\textrm{c}}}{63}a_{4},$
(a17)

here besides the coupling to $a_{2}$, it also couples to $a_{4}$ (which has no
simple relation of $\widehat{P}_{\textrm{c}}$ as $a_{2}$ does).

a.FIG. 1: Correction $Y$ as a function of $\widehat{P}_{\textrm{c}}$.

It is necessary to point out that if we truncate the expansion of the pressure
to order $b_{2}$, then only the combination $\widehat{P}_{\textrm{c}}$ comes
into play since
$b_{2}=-6^{-1}(1+3\widehat{P}_{\textrm{c}}^{2}+4\widehat{P}_{\textrm{c}})$,
i.e., it is independent of the coefficients $\\{a_{k}\\}$’s. Considering high-
order contributions from like $b_{4}$, certain model dependence on the dense
matter EOS may emerge. For instance, the small coefficients like $1/216,1/63$,
etc., in the expression of $b_{6}$ reduce the effects of these high-order
terms on the scaling of $\widehat{R}$ given in the main text. As an example,
we consider effects of the $b_{4}$-term on the dependence of $R$ versus
$\nu_{\textrm{c}}$. We write out

$\widehat{P}_{\textrm{c}}+b_{2}\widehat{R}^{2}+b_{4}\widehat{R}^{4}=0.$ (a18)

Using the expression for $\widehat{R}^{2}\approx
6^{-1}\widehat{P}_{\textrm{c}}(1+3\widehat{P}_{\textrm{c}}^{2}+4\widehat{P}_{\textrm{c}})$
from the $b_{2}$-order equation, one has

$\widehat{P}_{\textrm{c}}+\underbrace{b_{2}\widehat{R}^{2}+\underbrace{\left[\frac{\widehat{P}_{\textrm{c}}}{12}\left(1+3\widehat{P}_{\textrm{c}}^{2}+4\widehat{P}_{\textrm{c}}\right)-\frac{a_{2}}{30}\left(4+9\widehat{P}_{\textrm{c}}\right)\right]}_{b_{4}}\underbrace{\left[\frac{1}{6}\widehat{P}_{\textrm{c}}\left(1+3\widehat{P}_{\textrm{c}}^{2}+4\widehat{P}_{\textrm{c}}\right)\right]}_{\approx\widehat{R}^{2}}\widehat{R}^{2}}_{\overline{b}_{2}\widehat{R}^{2}}=0,$
(a19)

where the effective coefficient $\overline{b}_{2}$ is given by,

$\overline{b}_{2}=\underbrace{-\frac{1}{6}\left(1+3\widehat{P}_{\textrm{c}}^{2}+4\widehat{P}_{\textrm{c}}\right)}_{b_{2}}\times\underbrace{\left[1+\frac{a_{2}}{30}\widehat{P}_{\textrm{c}}\left(4+9\widehat{P}_{\textrm{c}}\right)-\frac{\widehat{P}_{\textrm{c}}^{2}}{12}\left(1+3\widehat{P}_{\textrm{c}}^{2}+4\widehat{P}_{\textrm{c}}\right)\right]}_{Y}=Yb_{2},$
(a20)

where $a_{2}=b_{2}/s_{\textrm{c}}^{2}$ with $s_{\textrm{c}}^{2}$ given by Eq.
(4) of the main text, a.FIG. 1 gives the correction $Y$ as a function of
$\widehat{P}_{\textrm{c}}$, from which one finds that $0.92\lesssim Y\lesssim
1$ for $0\lesssim\widehat{P}_{\textrm{c}}\lesssim 0.4$, e.g., $Y\approx 0.96$
for $\widehat{P}_{\textrm{c}}\approx 0.2$. In this sense, the $b_{4}$ term
does not essentially influence the scaling form for $R$ given by truncating
the equation of pressure to order $b_{2}$. The effects of the $b_{6}$-term
could be estimated similarly, and the effects on the dependence of $R$ versus
$\nu_{\textrm{c}}$ are small.

## B. Derivations of the Core Speed of Sound $s_{\textrm{c}}^{2}$

Using the radius scaling,

$\widehat{R}\sim\left(\frac{\widehat{P}_{\textrm{c}}}{1+3\widehat{P}_{\textrm{c}}^{2}+4\widehat{P}_{\textrm{c}}}\right)^{1/2},$
(a21)

we can write the NS mass $M_{\textrm{NS}}$ in the heuristic form
$M_{\textrm{NS}}\sim\widehat{R}^{3}/\sqrt{\varepsilon_{\textrm{c}}}\sim
F(\varepsilon_{\textrm{c}})$ as,

$F(\varepsilon_{\textrm{c}})=\frac{\Phi}{\sqrt{\varepsilon_{\textrm{c}}}}\left(\frac{\widehat{P}_{\textrm{c}}}{1+3\widehat{P}_{\textrm{c}}^{2}+4\widehat{P}_{\textrm{c}}}\right)^{3/2},$
(a22)

where $\Phi$ is a positive coefficient. Taking derivative of $F$ with respect
to $\varepsilon_{\textrm{c}}$ gives,

$\displaystyle\frac{\mathrm{d}F}{\mathrm{d}\varepsilon_{\textrm{c}}}=\frac{1}{2}\frac{F(\varepsilon_{\textrm{c}})}{\varepsilon_{\textrm{c}}}\left[3\left(\frac{\varepsilon_{\textrm{c}}}{P_{\textrm{c}}}\frac{\mathrm{d}{P}_{\textrm{c}}}{\mathrm{d}\varepsilon_{\textrm{c}}}-1\right)\frac{1-3\widehat{P}_{\textrm{c}}^{2}}{1+3\widehat{P}_{\textrm{c}}^{2}+4\widehat{P}_{\textrm{c}}}-1\right],$
(a23)

where one uses the expression for the derivative of $\widehat{R}$ with respect
to $\widehat{P}_{\textrm{c}}$,

$\frac{\mathrm{d}\widehat{R}}{\mathrm{d}\widehat{P}_{\textrm{c}}}=\widehat{R}\left(\frac{1-3\widehat{P}_{\textrm{c}}^{2}}{2\widehat{P}_{\textrm{c}}+8\widehat{P}_{\textrm{c}}^{2}+6\widehat{P}_{\textrm{c}}^{3}}\right),$
(a24)

and the basic relation,

$\frac{\mathrm{d}\widehat{P}_{\textrm{c}}}{\mathrm{d}\varepsilon_{\textrm{c}}}=\frac{1}{\varepsilon_{\textrm{c}}}\left(\frac{\mathrm{d}{P}_{\textrm{c}}}{\mathrm{d}\varepsilon_{\textrm{c}}}-\frac{P_{\textrm{c}}}{\varepsilon_{\textrm{c}}}\right).$
(a25)

Making $\mathrm{d}F/\mathrm{d}\varepsilon_{\textrm{c}}=0$ gives the speed of
sound
$s_{\textrm{c}}^{2}\equiv\mathrm{d}P_{\textrm{c}}/\mathrm{d}\varepsilon_{\textrm{c}}$
at the center of a NS of Eq. (4) in the main text. After obtaining the
$s_{\textrm{c}}^{2}$, we have,

$\frac{\partial
R}{\partial\varepsilon_{\textrm{c}}}\sim\frac{\partial}{\partial\varepsilon_{\textrm{c}}}\left(\frac{\widehat{R}}{\sqrt{\varepsilon_{\textrm{c}}}}\right)=\left(\frac{s_{\textrm{c}}^{2}}{\widehat{P}_{\textrm{c}}}-1\right)\frac{1-3\widehat{P}_{\textrm{c}}^{2}}{1+3\widehat{P}_{\textrm{c}}^{2}+4\widehat{P}_{\textrm{c}}}-1=-\frac{2}{3},$
(a26)

i.e., as $\varepsilon_{\textrm{c}}$ increases, the radius $R$ decreases, as
expected.

## C. Analyses using Different Radius Constraints on PSR J0740+6620

In Table 1, we show the radius constraints on PSR J0740+6620 using different
inference/analysis methods (Riley et al., 2021; Miller et al., 2021; Salmi et
al., 2022), as well as the resulted factor $\nu_{\textrm{c}}$, central energy
density $\varepsilon_{\textrm{c}}$, central pressure $P_{\textrm{c}}$ and the
central speed of sound squared $s_{\textrm{c}}^{2}$. The central EOS is shown
in a.FIG. 2, from which one can find that the constraining bands for
$P_{\textrm{c}}(\varepsilon_{\textrm{c}})$ in all three cases are narrow. The
uncertainties on $P_{\textrm{c}}(\varepsilon_{\textrm{c}})$ mainly reflect the
observational uncertainty on the radius of PSR J0740+6620.

Reference | radius (km) | $10^{3}\nu_{\textrm{c}}$ | $\varepsilon_{\textrm{c}}$ | $P_{\textrm{c}}$ | $s_{\textrm{c}}^{2}$
---|---|---|---|---|---
Riley et al. (2021) | $12.39_{-0.98}^{+1.30}$ | $11.2_{-0.9}^{+1.2}$ | $902_{-287}^{+214}$ | $218_{-125}^{+93}$ | $0.45_{-0.18}^{+0.14}$
Miller et al. (2021) | $13.7_{-1.5}^{+2.6}$ | $12.4_{-1.4}^{+2.5}$ | $656_{-339}^{+187}$ | $124_{-99}^{+53}$ | $0.32_{-0.14}^{+0.08}$
Salmi et al. (2022) | $12.90_{-0.97}^{+1.25}$ | $11.7_{-0.9}^{+1.2}$ | $794_{-235}^{+181}$ | $173_{-89}^{+69}$ | $0.39_{-0.13}^{+0.09}$
Table 1: Most probable values for the central energy density
$\varepsilon_{\textrm{c}}$, central pressure $P_{\textrm{c}}$ and the central
speed of sound squared $s_{\textrm{c}}^{2}$ for PSR J0740+6620, here its mass
about $2.08_{-0.07}^{+0.07}M_{\odot}$ and the corresponding
$12.2\lesssim{10^{4}\Gamma_{\textrm{c}}}/[{\textrm{fm}^{3/2}/\textrm{MeV}^{1/2}]}\lesssim
13.1$ are used for all three cases. Here $10^{3}\nu_{\textrm{c}}$ is measured
in unit $\textrm{fm}^{3/2}/\textrm{MeV}^{1/2}$, $\varepsilon_{\textrm{c}}$ and
$P_{\textrm{c}}$ are measured in $\textrm{MeV}/\textrm{fm}^{3}$.  a.FIG. 2:
Constraining bands for the central EOS of PSR J0740+6620 using different
radius constraints (Riley et al., 2021; Miller et al., 2021; Salmi et al.,
2022).

## References

  * Abbott et al. (2018) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2018, Phys. Rev. Lett., 121, 161101, doi: 10.1103/PhysRevLett.121.161101
  * Akmal et al. (1998) Akmal, A., Pandharipande, V. R., & Ravenhall, D. G. 1998, Phys. Rev. C, 58, 1804, doi: 10.1103/PhysRevC.58.1804
  * Al-Mamun et al. (2021) Al-Mamun, M., Steiner, A. W., Nättilä, J., et al. 2021, Phys. Rev. Lett., 126, 061101, doi: 10.1103/PhysRevLett.126.061101
  * Alford et al. (2005) Alford, M., Braby, M., Paris, M., & Reddy, S. 2005, The Astrophysical Journal, 629, 969, doi: 10.1086/430902
  * Alford et al. (2008) Alford, M. G., Schmitt, A., Rajagopal, K., & Schäfer, T. 2008, Rev. Mod. Phys., 80, 1455, doi: 10.1103/RevModPhys.80.1455
  * Altiparmak et al. (2022) Altiparmak, S., Ecker, C., & Rezzolla, L. 2022, The Astrophysical Journal Letters, 939, L34, doi: 10.3847/2041-8213/ac9b2a
  * Annala et al. (2018) Annala, E., Gorda, T., Kurkela, A., & Vuorinen, A. 2018, Phys. Rev. Lett., 120, 172703, doi: 10.1103/PhysRevLett.120.172703
  * Baiotti (2019) Baiotti, L. 2019, Progress in Particle and Nuclear Physics, 109, 103714, doi: https://doi.org/10.1016/j.ppnp.2019.103714
  * Bauswein et al. (2019) Bauswein, A., Bastian, N.-U. F., Blaschke, D. B., et al. 2019, Phys. Rev. Lett., 122, 061102, doi: 10.1103/PhysRevLett.122.061102
  * Baym et al. (2019) Baym, G., Furusawa, S., Hatsuda, T., Kojo, T., & Togashi, H. 2019, The Astrophysical Journal, 885, 42, doi: 10.3847/1538-4357/ab441e
  * Baym et al. (2018) Baym, G., Hatsuda, T., Kojo, T., et al. 2018, Reports on Progress in Physics, 81, 056902, doi: 10.1088/1361-6633/aaae14
  * Baym et al. (1971) Baym, G., Pethick, C., & Sutherland, P. 1971, ApJ, 170, 299, doi: 10.1086/151216
  * Bedaque & Steiner (2015) Bedaque, P., & Steiner, A. W. 2015, Phys. Rev. Lett., 114, 031103, doi: 10.1103/PhysRevLett.114.031103
  * Brandes et al. (2023) Brandes, L., Weise, W., & Kaiser, N. 2023, Phys. Rev. D, 107, 014011, doi: 10.1103/PhysRevD.107.014011
  * Breschi et al. (2022) Breschi, M., Bernuzzi, S., Godzieba, D., Perego, A., & Radice, D. 2022, Phys. Rev. Lett., 128, 161102, doi: 10.1103/PhysRevLett.128.161102
  * Cai & Li (2022) Cai, B.-J., & Li, B.-A. 2022, Nuclear Equation of State and Single-nucleon Potential from Gogny-like Energy Density Functionals Encapsulating Effects of Nucleon-nucleon Short-range Correlations.  https://arxiv.org/abs/2210.10924
  * Chin (1977) Chin, S. 1977, Annals of Physics, 108, 301, doi: https://doi.org/10.1016/0003-4916(77)90016-1
  * Chu et al. (2019) Chu, P.-C., Zhou, Y., Qi, X., et al. 2019, Phys. Rev. C, 99, 035802, doi: 10.1103/PhysRevC.99.035802
  * Collins & Perry (1975) Collins, J. C., & Perry, M. J. 1975, Phys. Rev. Lett., 34, 1353, doi: 10.1103/PhysRevLett.34.1353
  * Danielewicz et al. (2002) Danielewicz, P., Lacey, R., & Lynch, W. G. 2002, Science, 298, 1592, doi: 10.1126/science.1078070
  * De et al. (2018) De, S., Finstad, D., Lattimer, J. M., et al. 2018, Phys. Rev. Lett., 121, 091102, doi: 10.1103/PhysRevLett.121.091102
  * Dexheimer et al. (2021) Dexheimer, V., Noronha, J., Noronha-Hostler, J., Yunes, N., & Ratti, C. 2021, Journal of Physics G: Nuclear and Particle Physics, 48, 073001, doi: 10.1088/1361-6471/abe104
  * Dexheimer & Schramm (2008) Dexheimer, V., & Schramm, S. 2008, The Astrophysical Journal, 683, 943, doi: 10.1086/589735
  * Drischler et al. (2020) Drischler, C., Furnstahl, R. J., Melendez, J. A., & Phillips, D. R. 2020, Phys. Rev. Lett., 125, 202702, doi: 10.1103/PhysRevLett.125.202702
  * Ecker & Rezzolla (2022) Ecker, C., & Rezzolla, L. 2022, The Astrophysical Journal Letters, 939, L35, doi: 10.3847/2041-8213/ac8674
  * Engvik et al. (1996) Engvik, L., Osnes, E., Hjorth-Jensen, M., Bao, G., & Ostgaard, E. 1996, ApJ, 469, 794, doi: 10.1086/177827
  * Fonseca et al. (2021) Fonseca, E., Cromartie, H. T., Pennucci, T. T., et al. 2021, The Astrophysical Journal Letters, 915, L12, doi: 10.3847/2041-8213/ac03b8
  * Freedman & McLerran (1977a) Freedman, B. A., & McLerran, L. D. 1977a, Phys. Rev. D, 16, 1130, doi: 10.1103/PhysRevD.16.1130
  * Freedman & McLerran (1977b) —. 1977b, Phys. Rev. D, 16, 1147, doi: 10.1103/PhysRevD.16.1147
  * Freedman & McLerran (1977c) —. 1977c, Phys. Rev. D, 16, 1169, doi: 10.1103/PhysRevD.16.1169
  * Fujimoto et al. (2022) Fujimoto, Y., Fukushima, K., McLerran, L. D., & Praszałowicz, M. 2022, Phys. Rev. Lett., 129, 252702, doi: 10.1103/PhysRevLett.129.252702
  * Glendenning (1985) Glendenning, N. K. 1985, ApJ, 293, 470, doi: 10.1086/163253
  * Gorda et al. (2022) Gorda, T., Komoltsev, O., & Kurkela, A. 2022, Ab-initio QCD calculations impact the inference of the neutron-star-matter equation of state.  https://arxiv.org/abs/2204.11877
  * Hartle (1978) Hartle, J. B. 1978, Physics Reports, 46, 201, doi: https://doi.org/10.1016/0370-1573(78)90140-0
  * Hawking & Ellis (1973) Hawking, S., & Ellis, G. 1973, The Large Scale Structure of Space-time, Cambridge University Press, Section 9.1
  * Huth et al. (2022) Huth, S., Pang, P. T. H., Tews, I., et al. 2022, Nature, 606, 276, doi: 10.1038/s41586-022-04750-w
  * Iida & Sato (1997) Iida, K., & Sato, K. 1997, The Astrophysical Journal, 477, 294, doi: 10.1086/303685
  * Jiang & Yagi (2020) Jiang, N., & Yagi, K. 2020, Phys. Rev. D, 101, 124006, doi: 10.1103/PhysRevD.101.124006
  * Kalogera & Baym (1996) Kalogera, V., & Baym, G. 1996, The Astrophysical Journal, 470, L61, doi: 10.1086/310296
  * Komoltsev & Kurkela (2022) Komoltsev, O., & Kurkela, A. 2022, Phys. Rev. Lett., 128, 202701, doi: 10.1103/PhysRevLett.128.202701
  * Lackey et al. (2006) Lackey, B. D., Nayyar, M., & Owen, B. J. 2006, Phys. Rev. D, 73, 024021, doi: 10.1103/PhysRevD.73.024021
  * Lattimer (2021) Lattimer, J. M. 2021, Annual Review of Nuclear and Particle Science, 71, 433, doi: 10.1146/annurev-nucl-102419-124827
  * Lattimer & Prakash (2001) Lattimer, J. M., & Prakash, M. 2001, The Astrophysical Journal, 550, 426, doi: 10.1086/319702
  * Lattimer & Prakash (2005) —. 2005, Phys. Rev. Lett., 94, 111101, doi: 10.1103/PhysRevLett.94.111101
  * Legred et al. (2021) Legred, I., Chatziioannou, K., Essick, R., Han, S., & Landry, P. 2021, Phys. Rev. D, 104, 063003, doi: 10.1103/PhysRevD.104.063003
  * Li et al. (2020) Li, A., Zhu, Z.-Y., Zhou, E.-P., et al. 2020, Journal of High Energy Astrophysics, 28, 19, doi: https://doi.org/10.1016/j.jheap.2020.07.001
  * Li et al. (2008) Li, B.-A., Chen, L.-W., & Ko, C. M. 2008, Physics Reports, 464, 113, doi: https://doi.org/10.1016/j.physrep.2008.04.005
  * Lovato et al. (2022) Lovato, A., Dore, T., Pisarski, R. D., et al. 2022, Long Range Plan: Dense matter theory for heavy-ion collisions and neutron stars.  https://arxiv.org/abs/2211.02224
  * Margalit & Metzger (2017) Margalit, B., & Metzger, B. D. 2017, The Astrophysical Journal Letters, 850, L19, doi: 10.3847/2041-8213/aa991c
  * Miao et al. (2021) Miao, Z., Jiang, J.-L., Li, A., & Chen, L.-W. 2021, The Astrophysical Journal Letters, 917, L22, doi: 10.3847/2041-8213/ac194d
  * Miller et al. (2021) Miller, M. C., Lamb, F. K., Dittmann, A. J., et al. 2021, The Astrophysical Journal Letters, 918, L28, doi: 10.3847/2041-8213/ac089b
  * Montaña et al. (2019) Montaña, G., Tolós, L., Hanauske, M., & Rezzolla, L. 2019, Phys. Rev. D, 99, 103009, doi: 10.1103/PhysRevD.99.103009
  * Most et al. (2019) Most, E. R., Papenfort, L. J., Dexheimer, V., et al. 2019, Phys. Rev. Lett., 122, 061101, doi: 10.1103/PhysRevLett.122.061101
  * Müther et al. (1987) Müther, H., Prakash, M., & Ainsworth, T. 1987, Physics Letters B, 199, 469, doi: https://doi.org/10.1016/0370-2693(87)91611-X
  * Oertel et al. (2017) Oertel, M., Hempel, M., Klähn, T., & Typel, S. 2017, Rev. Mod. Phys., 89, 015007, doi: 10.1103/RevModPhys.89.015007
  * Oppenheimer & Volkoff (1939) Oppenheimer, J. R., & Volkoff, G. M. 1939, Phys. Rev., 55, 374, doi: 10.1103/PhysRev.55.374
  * Orsaria et al. (2019) Orsaria, M. G., Malfatti, G., Mariani, M., et al. 2019, Journal of Physics G: Nuclear and Particle Physics, 46, 073002, doi: 10.1088/1361-6471/ab1d81
  * Özel & Freire (2016) Özel, F., & Freire, P. 2016, Annual Review of Astronomy and Astrophysics, 54, 401, doi: 10.1146/annurev-astro-081915-023322
  * Pang et al. (2021) Pang, P. T. H., Tews, I., Coughlin, M. W., et al. 2021, The Astrophysical Journal, 922, 14, doi: 10.3847/1538-4357/ac19ab
  * Perego et al. (2022) Perego, A., Logoteta, D., Radice, D., et al. 2022, Phys. Rev. Lett., 129, 032701, doi: 10.1103/PhysRevLett.129.032701
  * Prakash et al. (1995) Prakash, M., Cooke, J. R., & Lattimer, J. M. 1995, Phys. Rev. D, 52, 661, doi: 10.1103/PhysRevD.52.661
  * Raaijmakers et al. (2021) Raaijmakers, G., Greif, S. K., Hebeler, K., et al. 2021, The Astrophysical Journal Letters, 918, L29, doi: 10.3847/2041-8213/ac089a
  * Radice et al. (2018) Radice, D., Perego, A., Zappa, F., & Bernuzzi, S. 2018, The Astrophysical Journal Letters, 852, L29, doi: 10.3847/2041-8213/aaa402
  * Rezzolla et al. (2018) Rezzolla, L., Most, E. R., & Weih, L. R. 2018, The Astrophysical Journal Letters, 852, L25, doi: 10.3847/2041-8213/aaa401
  * Rhoades & Ruffini (1974) Rhoades, C. E., & Ruffini, R. 1974, Phys. Rev. Lett., 32, 324, doi: 10.1103/PhysRevLett.32.324
  * Riley et al. (2021) Riley, T. E., Watts, A. L., Ray, P. S., et al. 2021, The Astrophysical Journal Letters, 918, L27, doi: 10.3847/2041-8213/ac0a81
  * Ruiz et al. (2018) Ruiz, M., Shapiro, S. L., & Tsokaros, A. 2018, Phys. Rev. D, 97, 021501, doi: 10.1103/PhysRevD.97.021501
  * Salmi et al. (2022) Salmi, T., Vinciguerra, S., Choudhury, D., et al. 2022, The Astrophysical Journal, 941, 150, doi: 10.3847/1538-4357/ac983d
  * Sathyaprakash et al. (2019) Sathyaprakash, B. S., Buonanno, A., Lehner, L., et al. 2019, Extreme Gravity and Fundamental Physics.  https://arxiv.org/abs/1903.09221
  * Serot & Walecka (1986) Serot, B. D., & Walecka, J. D. 1986, Adv. Nucl. Phys., 16, 1
  * Shuryak (1980) Shuryak, E. V. 1980, Physics Reports, 61, 71, doi: https://doi.org/10.1016/0370-1573(80)90105-2
  * Sorensen et al. (2023) Sorensen, A., Agarwal, K., Brown, K. W., et al. 2023, Dense Nuclear Matter Equation of State from Heavy-Ion Collisions.  https://arxiv.org/abs/2301.13253
  * Stone & Reinhard (2007) Stone, J., & Reinhard, P.-G. 2007, Progress in Particle and Nuclear Physics, 58, 587, doi: https://doi.org/10.1016/j.ppnp.2006.07.001
  * Tan et al. (2022a) Tan, H., Dexheimer, V., Noronha-Hostler, J., & Yunes, N. 2022a, Phys. Rev. Lett., 128, 161101, doi: 10.1103/PhysRevLett.128.161101
  * Tan et al. (2022b) Tan, H., Dore, T., Dexheimer, V., Noronha-Hostler, J., & Yunes, N. 2022b, Phys. Rev. D, 105, 023018, doi: 10.1103/PhysRevD.105.023018
  * Tews et al. (2018) Tews, I., Carlson, J., Gandolfi, S., & Reddy, S. 2018, The Astrophysical Journal, 860, 149, doi: 10.3847/1538-4357/aac267
  * Tolman (1939) Tolman, R. C. 1939, Phys. Rev., 55, 364, doi: 10.1103/PhysRev.55.364
  * Tsui & Leung (2005) Tsui, L. K., & Leung, P. T. 2005, Monthly Notices of the Royal Astronomical Society, 357, 1029, doi: 10.1111/j.1365-2966.2005.08710.x
  * Typel et al. (2010) Typel, S., Röpke, G., Klähn, T., Blaschke, D., & Wolter, H. H. 2010, Phys. Rev. C, 81, 015803, doi: 10.1103/PhysRevC.81.015803
  * Vautherin & Brink (1972) Vautherin, D., & Brink, D. M. 1972, Phys. Rev. C, 5, 626, doi: 10.1103/PhysRevC.5.626
  * Vidaña (2018) Vidaña, I. 2018, Proceedings of the Royal Society of London Series A, 474, 20180145, doi: 10.1098/rspa.2018.0145
  * Walecka (1974) Walecka, J. 1974, Annals of Physics, 83, 491, doi: https://doi.org/10.1016/0003-4916(74)90208-5
  * Watts (2019) Watts, A. L. 2019, in American Institute of Physics Conference Series, Vol. 2127, Xiamen-CUSTIPEN Workshop on the Equation of State of Dense Neutron-Rich Matter in the Era of Gravitational Wave Astronomy, 020008, doi: 10.1063/1.5117798
  * Watts et al. (2016) Watts, A. L., Andersson, N., Chakrabarty, D., et al. 2016, Rev. Mod. Phys., 88, 021001, doi: 10.1103/RevModPhys.88.021001
  * Weih et al. (2020) Weih, L. R., Hanauske, M., & Rezzolla, L. 2020, Phys. Rev. Lett., 124, 171103, doi: 10.1103/PhysRevLett.124.171103
  * Wen et al. (2019) Wen, D.-H., Li, B.-A., Chen, H.-Y., & Zhang, N.-B. 2019, Phys. Rev. C, 99, 045806, doi: 10.1103/PhysRevC.99.045806
  * Wiringa et al. (1988) Wiringa, R. B., Fiks, V., & Fabrocini, A. 1988, Phys. Rev. C, 38, 1010, doi: 10.1103/PhysRevC.38.1010
  * Xie & Li (2019) Xie, W.-J., & Li, B.-A. 2019, The Astrophysical Journal, 883, 174, doi: 10.3847/1538-4357/ab3f37
  * Xie & Li (2021) —. 2021, Phys. Rev. C, 103, 035802, doi: 10.1103/PhysRevC.103.035802
  * Xu et al. (2009) Xu, J., Chen, L.-W., Li, B.-A., & Ma, H.-R. 2009, The Astrophysical Journal, 697, 1549, doi: 10.1088/0004-637X/697/2/1549
  * Yagi & Yunes (2013) Yagi, K., & Yunes, N. 2013, Science, 341, 365, doi: 10.1126/science.1236462
  * Zhang & Li (2019) Zhang, N.-B., & Li, B.-A. 2019, European Physical Journal A, 55, 39, doi: 10.1140/epja/i2019-12700-0
  * Zhang & Li (2019) Zhang, N.-B., & Li, B.-A. 2019, The Astrophysical Journal, 879, 99, doi: 10.3847/1538-4357/ab24cb
  * Zhang & Li (2021) —. 2021, The Astrophysical Journal, 921, 111, doi: 10.3847/1538-4357/ac1e8c
  * Zhang & Chen (2016) Zhang, Z., & Chen, L.-W. 2016, Phys. Rev. C, 94, 064326, doi: 10.1103/PhysRevC.94.064326
  * Zhao & Lattimer (2020) Zhao, T., & Lattimer, J. M. 2020, Phys. Rev. D, 102, 023021, doi: 10.1103/PhysRevD.102.023021