# Exploring universal characteristics of neutron star matter with relativistic
ab initio equations of state

Sibo Wang  Department of Physics, Chongqing University, Chongqing 401331,
China  Chencan Wang  School of Physics, Nankai University, Tianjin 300071,
China  Hui Tong<EMAIL_ADDRESS>College of Physics and Materials Science,
Tianjin Normal University, Tianjin 300387, China  Strangeness Nuclear Physics
Laboratory, RIKEN Nishina Center, Wako, 351-0198, Japan

###### Abstract

Starting from the relativistic realistic nucleon-nucleon ($NN$) interactions,
the relativistic Brueckner-Hartree-Fock (RBHF) theory in the full Dirac space
is employed to study neutron star properties. First, the one-to-one
correspondence relation for gravitational redshift and mass is established and
used to infer the masses of isolated neutron stars by combining gravitational
redshift measurements. Next, the ratio of the moment of inertia $I$ to mass
times radius squared $MR^{2}$ as a function of the compactness $M/R$ is
obtained, and is consistent with the universal relations in the literature.
The moment of inertia for $1.338M_{\odot}$ pulsar PSR J0737-3039A
$I_{1.338M_{\odot}}$ is predicted to be 1.356$\times 10^{45}$, 1.381$\times
10^{45}$, and $1.407\times 10^{45}\ \mathrm{g~{}cm^{2}}$ by the RBHF theory in
the full Dirac space with $NN$ interactions Bonn A, B, and C, respectively.
Finally, the quadrupole moment of neutron star is calculated under the slow-
rotation and small-tidal-deformation approximation. The equations of state
constructed by the RBHF theory in the full Dirac space, together with those by
the projection method and momentum-independence approximation, conform to
universal $I$-Love-$Q$ relations as well. By combing the tidal deformability
from GW170817 and the universal relations from relativistic ab initio methods,
the moment of inertia of a neutron star with 1.4 solar mass is also deduced as
$I_{1.4M_{\odot}}=1.22^{+0.40}_{-0.25}\times 10^{45}\mathrm{g\ cm^{2}}$.

##  I Introduction

Neutron stars are one of the most compact objects in the universe: their
central densities can reach as high as 5 to 10 times the saturation density of
nuclear matter, $\rho_{0}\approx 0.16\ \text{fm}^{-3}$ Lattimer and Prakash
(2004), which is far beyond what can be achieved in terrestrial laboratories.
Therefore, neutron stars are ideal laboratories for studying ultradense
matter, and have established close connections among nuclear physics, particle
physics, and astrophysics.

The astrophysical observations of the global properties of neutron stars
provide important constraints for the equation of state (EOS) of dense matter
Lattimer and Prakash (2007); Özel and Freire (2016); Lattimer and Prakash
(2016); Burgio  _et al._ (2021), which is the only ingredient needed to unveil
the structure of neutron stars theoretically. The high-precision mass
measurements of massive neutron stars constitute nowadays one of the most
stringent astrophysical constraints on the nuclear EOS; such measurements
include PSR J1614-2230 ($1.928\pm 0.017M_{\odot}$) Demorest  _et al._ (2010);
Fonseca  _et al._ (2016), PSR J0348+0432 ($2.01\pm 0.04M_{\odot}$) Antoniadis
_et al._ (2013), and PSR J0740+6620 ($2.08\pm 0.07M_{\odot}$) Cromartie  _et
al._ (2020); Fonseca  _et al._ (2021). Recently, the Neutron star Interior
Composition Explorer (NICER) mission reported two independent Bayesian
parameter estimations of the mass and equatorial radius of the millisecond
pulsar PSR J0030+0451: $1.34_{-0.16}^{+0.15}\ M_{\odot}$ and
$12.71_{-1.19}^{+1.14}\ \text{km}$ Riley  _et al._ (2019) as well as
$1.44_{-0.14}^{+0.15}\ M_{\odot}$ and $13.02_{-1.06}^{+1.24}\ \text{km}$
Miller  _et al._ (2019). In combination with constraints from radio timing,
gravitational wave (GW) observations, and nuclear physics experiments, these
posterior distributions have been used to infer the properties of the dense
matter EOS (see Ref. Riley  _et al._ (2021) and references therein). Moreover,
two independent Bayesian estimations of the radius for the massive millisecond
pulsar PSR J0740+6620 have also been reported Riley  _et al._ (2021); Miller
_et al._ (2021).

Another unique probe for studying the properties of dense matter was extracted
from the recent observation of GW signals emitted from a binary neutron star
merger, i.e., GW170817 Abbott and et al. (2017). The tidal deformability,
which denotes the mass quadrupole moment response of a neutron star to the
strong external gravitational field induced by its companion Damour  _et al._
(1992); Hinderer (2008); Flanagan and Hinderer (2008); Damour and Nagar
(2009); Postnikov  _et al._ (2010), can be inferred from the GW signals. The
limits on the tidal deformability have been widely used to constrain the
neutron star radius Fattoyev  _et al._ (2018); Annala  _et al._ (2018); Most
_et al._ (2018); Tews  _et al._ (2018), the asymmetric nuclear matter EOS
Malik  _et al._ (2018); Zhang  _et al._ (2018); Tong  _et al._ (2020), and
hence the neutron skin thickness of 208Pb Fattoyev  _et al._ (2018).

Besides, as rotating objects, the internal structures of neutron stars are
strongly constrained by the moment of inertia, which can be determined from
the measurements of spin-orbit coupling in double pulsar systems Lyne  _et
al._ (2004). Such a measurement of the moment of inertia for neutron star
would have crucial implications for delimiting the EOS significantly Lattimer
and Schutz (2005) and can be used to distinguish neutron stars from quark
stars Yagi and Yunes (2013a). Special attention has been attracted by the
system PSR J0737-3039 Burgay  _et al._ (2003); Lyne  _et al._ (2004); Kramer
and Wex (2009), which is the only currently known double pulsar system. It is
hoped that the moment of inertia of the 1.338$M_{\odot}$ primary component in
this system, i.e., PSR J0737-3039A, will be measured eventually to within 10%
Kramer  _et al._ (2021), and could be used to impose new constraints on the
EOS Lattimer (2021).

The global properties of neutron stars, like masses, radii, tidal
deformabilities, and moments of inertia are highly sensitive to the EOS for
neutron star matter Lim  _et al._ (2019); Li and Sedrakian (2019); Tong  _et
al._ (2022). Nevertheless, it has been shown Yagi and Yunes (2013a, b) that,
for slowly rotating neutron stars, there exist universal relations between the
moment of inertia $I$, the tidal deformability $\Lambda$ (or Love number), and
the quadrupole moment $Q$ of neutron stars, i.e., the so-called $I$-Love-$Q$
relations, which are approximately independent of the internal composition and
the EOS for neutron star matter. Wide attention has been attracted by these
universal relations (see Ref. Yagi and Yunes (2017) for a review). Although so
far the reasons for these universal behaviors are not well understood Yagi
_et al._ (2014); Jiang and Yagi (2020), attempts have been made to combine the
universal relations with GW detections to infer neutron star properties Landry
and Kumar (2018); Kumar and Landry (2019).

The robustness of the $I$-Love-$Q$ relations has been extensively studied with
EOSs constructed from a variety of nuclear models (see Refs. Yagi and Yunes
(2017); Wei  _et al._ (2019) and references therein), including the
relativistic Brueckner-Hartree-Fock (RBHF) theory Müther  _et al._ (1987);
Brockmann and Machleidt (1990); Li  _et al._ (1992); Engvik  _et al._ (1996);
Gross-Boelting  _et al._ (1999); Katayama and Saito (2015). Since the 1980s
the RBHF theory has played an important role in understanding the properties
of dense nuclear matter from realistic nucleon-nucleon ($NN$) interactions
Anastasio  _et al._ (1980, 1981). In the RBHF theory, the single-particle
motion of the nucleon in nuclear matter is described with the Dirac equation,
where the medium effects are absorbed into the single-particle potentials. In
principle, the scalar and the vector components of the single-particle
potentials should be determined in the full Dirac space Nuppenau  _et al._
(1989), i.e., by considering the positive-energy states (PESs) and negative-
energy states (NESs) simultaneously. However, to avoid the difficulties
induced by NESs, the RBHF calculations are primarily performed in the Dirac
space without NESs Müther  _et al._ (1987); Brockmann and Machleidt (1990); Li
_et al._ (1992); Engvik  _et al._ (1996); Horowitz and Serot (1984); Gross-
Boelting  _et al._ (1999).

Recently, a self-consistent RBHF calculation in the full Dirac space was
achieved for symmetric nuclear matter (SNM) Wang  _et al._ (2021, 2022a) and
asymmetric nuclear matter (ANM) Wang  _et al._ (2022b). By decomposing the
matrix elements of single-particle potential operator in the full Dirac space,
the momentum-dependent scalar and vector components of the single-particle
potentials are determined uniquely Wang  _et al._ (2021). The long-standing
controversy about the isospin dependence of the effective Dirac mass in
relativistic ab initio calculations of ANM is also clarified Wang  _et al._
(2022b).

The RBHF theory in the full Dirac space has been applied to neutron stars Wang
_et al._ (2022b); Tong  _et al._ (2022), where the mass, radius, and tidal
deformability are calculated with realistic $NN$ interactions Bonn A, B, and C
Machleidt (1989). The maximum mass of a neutron star is found less than
$2.4M_{\odot}$ and the neutron star radius for $1.4M_{\odot}$ is predicted
about 12 km, which are consistent with the astrophysical observations of
massive neutron stars and simultaneous mass-radius estimations by NICER Miller
_et al._ (2019). The tidal deformabilities for a $1.4M_{\odot}$ neutron star
are predicted as 376, 473, and 459 for the three parametrizations of $NN$
interactions respectively, and all lie in the region
$\Lambda_{1.4M_{\odot}}=190^{+390}_{-120}$ inferred from the revised analysis
by LIGO and Virgo Collaborations Abbott and et al. (2018).

In this work, we employ the RBHF theory in the full Dirac space to study other
global properties of neutron stars, including the gravitational redshift,
moment of inertia, and quadrupole moment under the slow-rotation and small-
tidal-deformation approximation. The main focus will be the relation between
the moment of inertia and the compactness parameter, as well as the universal
$I$-Love-$Q$ relations. This paper is organized as follows. In Sec. II, the
theoretical framework of the RBHF theory and structure equations for neutron
star properties are briefly described. The obtained results and discussions
are presented in Sec. III. Finally, a summary is given in Sec. IV.

##  II Theoretical framework

###  II.1 The relativistic Brueckner-Hartree-Fock theory

In the RBHF calculations, one of the most important procedure is the self-
consistent determination of the single-particle potential operator
$\mathcal{U}$ of the nucleons, which is generally divided into scalar and
vector components Serot and Walecka (1986)

$\mathcal{U}(\bm{p})=U_{S}(p)+\gamma^{0}U_{0}(p)+\bm{\gamma\cdot\hat{p}}\,U_{V}(p).$
(1)

Here $\hat{\bm{p}}=\bm{p}/p$ is the unit vector parallel to the momentum
$\bm{p}$. The quantities $U_{S}(p)$, $U_{0}(p)$, and $U_{V}(p)$ are the scalar
potential, the timelike part and the spacelike part of the vector potential.

In principle, the scalar and the vector components of the single-particle
potentials can only be determined uniquely in the full Dirac space. However,
to avoid the numerical difficulties in the full Dirac space, different
approximations are proposed to extract the single-particle potentials in the
Dirac space without NESs. The momentum-independence approximation Brockmann
and Machleidt (1990) assumes that the single-particle potentials are
independent of the momentum, and the spacelike component of the vector
potential, $U_{V}$, is negligible. The scalar potential $U_{S}$ and the
timelike part of the vector potential, $U_{0}$, are then extracted from the
single-particle potential energies at two selected momenta. In the projection
method Gross-Boelting  _et al._ (1999), the effective $NN$ interaction $G$
matrix, which is obtained by solving the in-medium scattering equation, is
projected onto a complete set of five Lorentz invariant amplitudes, from which
the single-particle potentials are calculated analytically. However, the
choice of these Lorentz invariant amplitudes is not unique.

Only by decomposing the matrix elements of $\mathcal{U}$ in the full Dirac
space, can the Lorentz structure and momentum dependence of single-particle
potentials be uniquely determined Nuppenau  _et al._ (1989). The theoretical
framework for the RBHF theory in the full Dirac space has been described in
detail in Ref. Wang  _et al._ (2021) for SNM and Ref. Tong  _et al._ (2022)
for ANM. In this work this method is used to construct the EOS of neutron star
matter, which is regarded as beta equilibrium nuclear matter consisting of
protons, neutrons, electrons, and muons Krastev and Sammarruca (2006). Using
the relativistic Bonn A potential Machleidt (1989), the RBHF theory in the
full Dirac space for nuclear matter is applicable for density in the range
0.08-0.57 $\text{fm}^{-3}$. For lower density in the crust of a neutron star,
the EOS introduced with the Baym-Bethe-Pethick (BBP) Baym  _et al._ (1971) and
Baym-Pethick-Sutherland (BPS) model Baym  _et al._ (1971) is used. For higher
density, we follow the strategy proposed in Ref. Rhoades and Ruffini (1974)
and applied in Refs. Gandolfi  _et al._ (2012); Wang  _et al._ (2022b), where
the neutron-star matter EOS above a critical density $\rho_{c}=0.57\
\text{fm}^{-3}$ is replaced with the maximally stiff or causal one, which
predicts the most rapid increase of pressure with energy density without
violating the causality limit.

###  II.2 Mass, radius, gravitational redshift, and tidal deformability

The stable configurations of a cold, spherically symmetric, and nonrotating
neutron star can be obtained from the Tolman-Oppenheimer-Volkov (TOV)
equations Oppenheimer and Volkoff (1939); Tolman (1939). Adopting natural
units $G=c=1$, the TOV equations are given by

$\displaystyle\frac{dP(r)}{dr}=$ $\displaystyle\
-\frac{[P(r)+\mathcal{E}(r)][M(r)+4\pi r^{3}P(r)]}{r[r-2M(r)]},$ (2a)
$\displaystyle\frac{dM(r)}{dr}=$ $\displaystyle\ 4\pi r^{2}\mathcal{E}(r),$
(2b)

where $P(r)$ is the pressure at neutron star radius $r$, $M(r)$ is the total
neutron star mass inside a sphere of radius $r$, and $\mathcal{E}(r)$ is the
total energy density. These differential equations can be solved numerically
with a given central pressure $P_{c}$ and $M(0)=0$. The quantity $R$ for
$P(R)=0$ denotes the radius of the neutron star, and $M(R)$ is its mass. The
gravitational redshift which relates the mass of the neutron star to its
radius is defined as

$z=\left(1-2M/R\right)^{-1/2}-1.$ (3)

Since the radius of the neutron star is harder to observe relative to its
mass, the simultaneous measurements of the mass and the gravitational redshift
would provide a clear radius determination.

The tidal deformability is defined as

$\Lambda=\frac{2}{3}k_{2}C^{-5}.$ (4)

where $C=M/R$ is the compactness parameter. The second Love number $k_{2}$
Hinderer (2008); Hinderer  _et al._ (2010) is calculated by

$\begin{split}k_{2}=&\
\frac{8C^{5}}{5}(1-2C)^{2}[2-y_{R}+2C(y_{R}-1)]\times\left\\{6C[2-y_{R}+C(5y_{R}-8)]\right.\\\
&\ +4C^{3}[13-11y_{R}+C(3y_{R}-2)+2C^{2}(1+y_{R})]\\\ &\
\left.+3(1-2C)^{2}[2-y_{R}+2C(y_{R}-1)]\ln(1-2C)\right\\}^{-1},\end{split}$
(5)

where $y_{R}=y(R)$ is the solution of the following nonlinear, first-order
differential equation

$r\frac{dy(r)}{dr}+y^{2}(r)+F(r)y(r)+r^{2}Q(r)=0.$ (6)

Here the two functions $F(r)$ and $Q(r)$ depend on the known mass, radius,
pressure, and energy density profiles of the star:

$\displaystyle F(r)=$ $\displaystyle\
\left[1-\frac{2M(r)}{r}\right]^{-1}\left\\{1-4\pi
r^{2}[\mathcal{E}(r)-P(r)]\right\\},$ (7a) $\displaystyle Q(r)=$
$\displaystyle\
\left\\{4\pi\left[5\mathcal{E}(r)+9P(r)+\frac{\mathcal{E}(r)+P(r)}{\partial
P/\partial\mathcal{E}}\right]-\frac{6}{r^{2}}\right\\}\times\left[1-\frac{2M(r)}{r}\right]^{-1}$
$\displaystyle\ -\left[\frac{2M(r)}{r^{2}}+8\pi
rP(r)\right]^{2}\times\left[1-\frac{2M(r)}{r}\right]^{-2}$ (7b)

The differential equation (6) for $k_{2}$ can be solved together with the TOV
equations and the initial condition $y(0)=2$.

###  II.3 The moment of inertia

The moment of inertia is calculated under the slow-rotation approximation
pioneered by Hartle and Thorne Hartle (1967); Hartle and Thorne (1968), where
the frequency $\Omega$ of a uniformly rotating neutron star is far smaller
than the Kepler frequency at the equator:

$\Omega\ll\Omega_{\text{max}}\simeq\sqrt{M/R^{3}}.$ (8)

In the slow-rotation approximation the moment of inertia of a uniformly
rotating, axially symmetric neutron star is given by the following expression
Fattoyev and Piekarewicz (2010):

$I=\frac{8\pi}{3}\int^{R}_{0}r^{4}e^{-\nu(r)}\frac{\bar{\omega}(r)}{\Omega}\frac{\mathcal{E}(r)+P(r)}{\sqrt{1-2M(r)/r}}dr.$
(9)

The quantity $\nu(r)$ is a radially dependent metric function and is defined
as

$\nu(r)=\frac{1}{2}\ln\left(1-\frac{2M}{R}\right)-\int^{R}_{r}\frac{M(x)+4\pi
x^{3}P(x)}{x^{2}[1-2M(x)/x]}dx.$ (10)

The frame-dragging angular velocity $\bar{\omega}$ is usually obtained by the
dimensionless relative frequency $\tilde{w}\equiv\bar{\omega}/\Omega$, which
satisfies the following second-order differential equation:

$\frac{d}{dr}\left[r^{4}j(r)\frac{d\tilde{\omega}(r)}{dr}\right]+4r^{3}\frac{dj(r)}{dr}\tilde{\omega}(r)=0,$
(11)

where $j(r)=e^{-\nu(r)}\sqrt{1-2M(r)/r}$ for $r\leq R$. The relative frequency
$\tilde{\omega}(r)$ is subject to the following two boundary conditions

$\displaystyle\tilde{\omega}^{\prime}(0)=$ $\displaystyle\ 0,$ (12a)
$\displaystyle\tilde{\omega}(R)+\frac{R}{3}\tilde{\omega}^{\prime}(R)=$
$\displaystyle\ 1.$ (12b)

It should be noted that under the slow-rotation approximation the moment of
inertia does not depend on the stellar frequency $\Omega$.

###  II.4 The quadrupole moment

It has been shown Yagi and Yunes (2013a, b) that there exist universal
relations between the moment of inertia, the Love number, and the quadrupole
moment of neutron stars. Physically, the moment of inertia quantifies how fast
a neutron star can spin for a fixed angular momentum, the quadrupole moment
describes how much a neutron star is deformed away from sphericity due to
rotation, and the Love number characterizes how easily a neutron star can be
deformed due to an external tidal field. These quantities can be computed by
numerically solving for the interior and exterior gravitational field of a
neutron star in a slow-rotation Hartle (1967); Hartle and Thorne (1968)
approximation and in a small-tidal-deformation approximation Hinderer (2008);
Hinderer  _et al._ (2010). In this work the quadrupole moment is calculated by
following the detailed instructions described in Ref. Yagi and Yunes (2013b).
In order to investigate the universal $I$-Love-$Q$ relations, the following
dimensionless quantities are introduced:

$\bar{I}\equiv\frac{I}{M^{3}},\qquad\bar{Q}\equiv-\frac{QM}{(I\Omega)^{2}}.$
(13)

##  III Results and discussions

Figure 1: (Color online) The speed of sound squared for neutron star matter as
a function of density (left) and the neutron star mass-radius relation (right)
obtained by the RBHF theory with the potential Bonn A in the full Dirac space
(red solid line), in comparison with those obtained by the projection method
(green dashed line) and the momentum-independence approximation (gray dotted
line). The dark (light) blue and purple regions indicate the 68% (95%)
confidence intervals constrained by the NICER analysis of PSR J0030+0451
Miller  _et al._ (2019) and PSR J0740+6620 Miller  _et al._ (2021). The red
dotted line corresponds to the maximally stiff EOS.

In the left panel of Fig. 1, the speed of sound squared for neutron star
matter from the RBHF theory with potential Bonn A is depicted as a function of
density. The results calculated with the projection method and momentum-
independence approximation show slowly increasing tendencies with the increase
of density. In the full Dirac space, the speed of sound squared increases
quickly with the density, and reaches 0.77 at $0.57\ \text{fm}^{-3}$. Above
this density, the amplitude of scalar potential $U_{S}$ exceeds the nucleon
rest mass, and the RBHF iteration in the full Dirac space is very difficult to
achieve. This fact might be related to the extension of the Bonn potential to
the full Dirac space, which was determined with PESs only. When the central
density of a neutron star is fixed at $\rho_{c}=0.57\ \text{fm}^{-3}$, our
calculations could support a neutron star with mass equal to $1.9M_{\odot}$.
To further explore the maximum mass of a neutron star from the RBHF theory in
the full Dirac space, we continue with an EOS where the speed of sound is
equal to the speed of light Rhoades and Ruffini (1974); Gandolfi  _et al._
(2012); Wang  _et al._ (2022b) (red dotted line). This would provide an upper
bound on the maximum mass of a neutron star.

Based on the EOS from the RBHF theory, the mass-radius relations of a neutron
star can be calculated from the TOV equations, which are shown in the right
panel of Fig. 1. The 68% and 95% contours of the joint probability density
distribution of the mass and radius of PSR J0030+0451 Miller  _et al._ (2019)
and PSR J0740+6620 Miller  _et al._ (2021) from the NICER analysis are also
shown. It can be found that the results obtained in the full Dirac space, with
the projection method, and with the momentum-independence approximation are
consistent with the recent constraints by NICER. The maximum masses
$M_{\text{max}}$ predicted by the three methods are $2.43M_{\odot}$,
$2.31M_{\odot}$, and $2.18M_{\odot}$ respectively, which are consistent with
the available astrophysical constraints from massive neutron star
observations, such as PSR J1614-2230 ($1.928\pm 0.017M_{\odot}$) Demorest  _et
al._ (2010); Fonseca  _et al._ (2016), PSR J0348+0432 ($2.01\pm
0.04M_{\odot}$) Antoniadis  _et al._ (2013), and PSR J0740+6620 ($2.08\pm
0.07M_{\odot}$) Cromartie  _et al._ (2020); Fonseca  _et al._ (2021). The
radii $R_{1.4M_{\odot}}$ of a $1.4M_{\odot}$ neutron star from the three
methods are $11.97\ \text{km}$, $12.38\ \text{km}$, and $12.35\ \text{km}$,
respectively.

Table 1: Neutron star properties and nuclear matter properties at saturation density calculated by the RBHF theory in the full Dirac space with potential Bonn A, in comparison with the results obtained by the RBHF calculation with the projection method and the momentum-independence approximation. Model | $M_{\text{max}}$ | $R_{1.4M_{\odot}}$ | $\rho_{1.4M_{\odot}}$ | $\rho_{0}$ | $E/A$ | $E_{\mathrm{sym}}$ | $L$
---|---|---|---|---|---|---|---
($M_{\odot}$) | (km) | (fm-3) | (fm-3) | (MeV) | (MeV) | (MeV)
Full Dirac space | 2.43 | 11.97 | 0.46 | 0.188 | -15.40 | 33.1 | 65.2
Projection method | 2.31 | 12.38 | 0.42 | 0.179 | -16.15 | 34.7 | 68.8
Mom.-ind. app. | 2.18 | 12.35 | 0.43 | 0.178 | -15.36 | 33.2 | 67.3

In Table 1 we summarize the maximum mass of a neutron star and the radius and
the central density for a $1.4M_{\odot}$ neutron star obtained by the RBHF
theory in the full Dirac space, together with those obtained with the
projection method and momentum-independence approximation. The smallest value
for $R_{1.4M_{\odot}}$ found in the full Dirac space corresponds to the
softest EOS for neutron star matter below a density of about 0.4 fm-3. This
fact can be further related to nuclear matter properties, which are also
summarized in Tab. 1, including the binding energy per nucleon $E/A$, the
symmetry energy $E_{\text{sym}}$, and its slope $L$. It is found that the
symmetry energy and its slope at the saturation density obtained in the full
Dirac space are the smallest among the three methods. This explains the reason
for smallest $R_{1.4M_{\odot}}$ and shows how important it is to take both the
PESs and NESs into account.

Figure 2: (Color online) The relation between the gravitational redshift $z$
of a neutron star with its mass $M$ (left) and radius $R$ (right) obtained by
the RBHF theory in the full Dirac space with the potentials Bonn A, B, and C.
Astrophysical observation for the gravitational redshift of isolated neutron
stars RBS 1223 Hambaryan  _et al._ (2014), RX J1856.5-3754 Hambaryan  _et al._
(2014), and RX J0720.4-3125 Hambaryan, V.  _et al._ (2017) are also shown.

Figure 2 shows the relations between the gravitational redshift $z$ of a
neutron star with its mass $M$ and radius $R$ obtained by the RBHF theory in
the full Dirac space with the potentials Bonn A, B, and C. It can be seen that
the $z$-$M$ relation is not strongly affected by the $NN$ interactions. With
the increase of mass, the gravitational redshift shows a monotonically
increasing behavior. A backbending phenomenon is found for the gravitational
redshift with the decrease of radius, similar to the case for $M$-$R$
relations as shown in Fig. 1.

Table 2: The masses of isolated neutron star RBS 1223 Hambaryan  _et al._ (2014), RX J1856.5-3754 Hambaryan  _et al._ (2014), and RX J0720.4-3125 Hambaryan, V.  _et al._ (2017) predicted by combining gravitational redshift measurements and the RBHF theory in the full Dirac space with potential Bonn A. The radii are also given. The uncertainties in the last two columns are all from the gravitational redshift measurements. System | Gravitational redshift $z$ | Mass $M[M_{\odot}]$ | Radius $R[\text{km}]$
---|---|---|---
RBS 1223 Hambaryan  _et al._ (2014) | $0.16^{+0.03}_{-0.02}$ | $1.03^{+0.15}_{-0.11}$ | $11.85^{+0.05}_{-0.04}$
RX J1856.5-3754 Hambaryan  _et al._ (2014) | $0.22^{+0.06}_{-0.12}$ | $1.33^{+0.25}_{-0.64}$ | $11.94^{+0.03}_{-0.18}$
RX J0720.4-3125 Hambaryan, V.  _et al._ (2017) | $0.205^{+0.006}_{-0.003}$ | $1.258^{+0.028}_{-0.014}$ | $11.922^{+0.008}_{-0.004}$

The clear one-to-one correspondence relation for gravitational redshift and
mass established in the left panel in Fig. 2 can be used to infer the mass of
a isolated neutron star, when the observation of the gravitational redshift is
provided. In Fig. 2, astrophysical observations of gravitational redshift for
isolated neutron stars RX J0720.4-3125 Hambaryan, V.  _et al._ (2017), RBS
1223 Hambaryan  _et al._ (2014), and RX J1856.5-3754 Hambaryan  _et al._
(2014) are shown as shadow bands. The predicted masses by combining these
observations and the theoretical calculations from the RBHF theory in the full
Dirac space with Bonn A potential are listed in the third column of Table 2.
The uncertainties are from the gravitational redshift measurements. In Table 2
the corresponding radii are also shown in the last column.

Figure 3: (Color online) The ratio of the moment of inertia $I$ to $MR^{2}$ as
a function of the compactness parameter $M/R$ obtained with the RBHF theory in
the full Dirac space with $NN$ interactions Bonn A, B, and C. The probability
distribution (orange) from the Bayesian analysis together with the fitted
(orange) line in Ref. Lim  _et al._ (2019) are shown for comparison. The
empirical (purple) band from Ref. Lattimer and Schutz (2005) is also shown.
See the text for details.

In Fig. 3 we display the ratio of the moment of inertia $I$ to $MR^{2}$ as a
function of the compactness parameter $M/R$ obtained by the RBHF theory in the
full Dirac space with $NN$ interactions Bonn A, B, and C. Lattimer et al.
Lattimer and Schutz (2005) showed that, in the absence of phase transition and
other effects that strongly soften the EOS at supranuclear densities, there is
a relatively unique relation between the quantity $I/MR^{2}$ and $M/R$:

$I/MR^{2}\simeq(0.237\pm 0.008)(1+2.844C+18.91C^{4}).$ (14)

This relation is shown as the purple band in Fig. 3. It is found that our
results are consistent with the universal relations obtained in Ref. Lattimer
and Schutz (2005) for the range where $M/R>0.08\ M_{\odot}/\text{km}$. The
derivation for smaller compactness is unimportant, since the observational
evidence for neutron star masses and radii lie in the ranges of $1.2\
M_{\odot}\leq M\leq 2.2\ M_{\odot}$ and $9\ \text{km}\leq R\leq 15\ \text{km}$
respectively, which leads to the range of compactness as $0.08\
M_{\odot}/\text{km}<M/R<0.24\ M_{\odot}/\text{km}$. Lim and collaborators Lim
_et al._ (2019) have investigated neutron star moments of inertia from
Bayesian posterior probability distributions of the nuclear EOSs that
incorporate information from microscopic many-body theory and empirical data
of finite nuclei. The probability distribution for $I/MR^{2}$ is shown as in
Fig. 3. They found that over the entire range of neutron star compactness
$0<C\leq 0.25\ M_{\odot}/\text{km}$, their results can be well fitted with the
following formula:

$I/MR^{2}=\frac{C+27.178C^{4}}{0.0871+2.183C}.$ (15)

This result is also shown as the orange dashed line in Fig. 3. It can be seen
that our results are very close to that obtained by Lim et al., especially for
the neutron stars with small compactness.

Table 3: The momenta of inertia for $1.338\ M_{\odot}$ pular PSR J0737-3039A predicted by the RBHF theory in the full Dirac space with potentials Bonn A, B, and C. The results obtained with the projection method Gross-Boelting  _et al._ (1999); van Dalen  _et al._ (2004) and momentum-independence approximation (Mom.-ind. approx.) Brockmann and Machleidt (1990) are also shown. For comparison, the values obtained with nonrelativistic ab initio calculation Morrison  _et al._ (2004), predicted with the tidal deformability in GW170817 combing universal relations Landry and Kumar (2018) as well as that inferred from Bayesian analysis (95% credibility) Lim  _et al._ (2019) are also shown. The last column is the corresponding radius. Model | Potential | $I_{1.338M_{\odot}}~{}[10^{45}\ \mathrm{g~{}cm^{2}}]$ |  $R_{1.338M_{\odot}}$ [km]
---|---|---|---
| A | 1.356 | 11.94
Full Dirac space | B | 1.381 | 12.11
| C | 1.407 | 12.26
| A | 1.440 | 12.34
Projection method | B | 1.465 | 12.49
| C | 1.487 | 12.60
| A | 1.431 | 12.32
Mom.-ind. approx. | B | 1.452 | 12.46
| C | 1.471 | 12.57
Variational calculation (APR) Morrison  _et al._ (2004) | AV18 $+\ \delta v$ $+$ UIX* | 1.24 | 11.56
GW170817 $+$ universal relations Landry and Kumar (2018) |  | $1.15^{+0.38}_{-0.24}$ |
Bayesian analysis Lim  _et al._ (2019) |  | $1.36^{+0.15}_{-0.32}$ | $12.2^{+0.7}_{-1.9}$

Although the ratio of the neutron star moment of inertia $I$ to $MR^{2}$ has a
universal function of the compactness parameter $M/R$, the moment of inertia
itself depends sensitively on the neutron star’s internal structure. It has
been suggested Lattimer and Schutz (2005) that a measurement accuracy of 10%
for $I$ is sufficient to place strong constraints on the EOS.

In Table 3, we show the momenta of inertia $I_{1.338M_{\odot}}$ and radius
$R_{1.338M_{\odot}}$ for PSR J0737-3039A predicted by the RBHF theory in the
full Dirac space with three parametrizations for $NN$ interactions. The
results obtained by the projection method Gross-Boelting  _et al._ (1999); van
Dalen  _et al._ (2004) and momentum-independence approximation Brockmann and
Machleidt (1990) are also shown. The RBHF theory in the full Dirac space leads
to minimum values compared to the approximations in the Dirac space without
NESs. This is understandable since the RBHF theory in the full Dirac space
gives the minimum radius of a neutron star for the fixed canonical mass, as
shown in Table 1.

The moment of inertia for PSR J0737-3039A predicted by the RBHF theory in the
full Dirac space with Bonn A is $1.356\times 10^{45}\ \mathrm{g~{}cm^{2}}$,
which is very close to the most probable value $1.36\times 10^{45}\
\mathrm{g~{}cm^{2}}$ obtained from Bayesian analysis (95% credibility) Lim
_et al._ (2019). The result from the nonrelativistic ab initio variational
calculations Akmal  _et al._ (1998) is also shown in Table 3, where the
Argonne $v18$ interaction (AV18) Wiringa  _et al._ (1995) is used, together
with the relativistic boost corrections to the two-nucleon interaction as well
as three-nucleon interactions modeled with the Urbana force Pudliner  _et al._
(1995). The nonrelativistic ab initio calculation leads to a moment of inertia
smaller that what we obtain, similarly to the case for radius. In Ref. Landry
and Kumar (2018), by using well-known universal relations among neutron star
observables, the reported 90% credible bound on the tidal deformability
$\Lambda_{1.4M_{\odot}}=190^{+390}_{-120}$ from GW170817 Abbott and et al.
(2018) has been translated into a direct constraint on the moment of inertia
of PSR J0737-3039A, giving $I_{1.338M_{\odot}}=1.15^{+0.38}_{-0.24}\times
10^{45}\ \mathrm{g~{}cm^{2}}$. It can be seen that the results with the three
methods for RBHF theory are consistent with this constraint.

Figure 4: (Color online) Top: The $I$-Love (left) and $Q$-Love (right)
relations for slowly-rotating neutron stars with EOSs obtained by the RBHF
theory in the full Dirac space, projection method, and momentum-independence
approximations with potentials Bonn A, B, and C. The solid curves show the fit
by using Eq. (16). Bottom: Absolute fractional difference from the fit.

Let us now compare the EOSs obtained by the relativistic ab initio
calculations, i.e., the RBHF theory in the full Dirac space, projection
method, and momentum-independence approximation with Bonn potentials, to the
universal $I$-Love-$Q$ relations. The $I$-Love as well as $Q$-Love relations
and $I$-$Q$ relations are shown in the top panels of Figs. 4 and 5,
respectively. The dimensionless moment of inertia $\bar{I}$ and dimensionless
quadrupole moment $\bar{Q}$ are defined in Eq. (15). A single parameter along
the curve is the mass or compactness, which increases to the left of the
plots. Similarly to Ref. Yagi and Yunes (2017), we only show data with the
mass of an isolated, nonrotating configuration in the range
$1M_{\odot}<M<M_{\text{max}}$ with $M_{\text{max}}$ representing the maximum
mass for such a configuration. One observes that the universal relations hold
very well. Since the relations are insensitive to EOS, one can construct a
single fit (black solid curves) given by Yagi and Yunes (2013b, 2017)

$\ln y_{i}=a_{i}+b_{i}\ln x_{i}+c_{i}(\ln x_{i})^{2}+d_{i}(\ln
x_{i})^{3}+e_{i}(\ln x_{i})^{4},$ (16)

where coefficients are listed in Table 4. These coefficients are very close to
that in Ref. Yagi and Yunes (2017), where a large number of EOSs are
considered. The bottom panels of Figs. 4 and 5 show the absolute fractional
difference between all the data and the fit, which is less than $1\%$ in the
whole range.

The universal relations between $\bar{I}$ and $\Lambda$ allows one to extract
the momenta of inertia of a neutron star with $1.4$ solar mass,
$\bar{I}_{1.4M_{\odot}}$, from the tidal deformability
$\Lambda_{1.4M_{\odot}}$ from GW170817. The revised analysis from LIGO and
Virgo Collaborations, $\Lambda_{1.4M_{\odot}}=190^{+390}_{-120}$ Abbott and et
al. (2018), corresponds to $\bar{I}_{1.4M_{\odot}}=10.30^{+3.39}_{-2.10}$ as
shown in the left panel of Fig. 4. From $\bar{I}_{1.4M_{\odot}}$ and the
relation $\bar{I}=I/M^{3}$ we obtain
$I_{1.4M_{\odot}}=1.22^{+0.40}_{-0.25}\times 10^{45}\mathrm{g\ cm^{2}}$. These
values are consistent with the results
$\bar{I}_{1.4M_{\odot}}=11.10^{+3.64}_{-2.28}$ and
$I_{1.4M_{\odot}}=1.15^{+0.38}_{-0.24}\times 10^{45}\mathrm{g\ cm^{2}}$ in
Ref. Landry and Kumar (2018), where the $I$-Love relation is obtained by using
a large set of candidate neutron star EOSs based on relativistic mean-field
and Skyrme-Hartree-Fock theory.

Figure 5: (Color online) The same as the left panel in Fig. 4, but for $I$-$Q$ relations. Table 4: Numerical coefficients for the fitting formula of the $I$-Love, $I$-$Q$, and $Q$-Love relations given in Eq. (16).  $y_{i}$ | $x_{i}$ | $a_{i}$ | $b_{i}$ | $c_{i}$ | $d_{i}$ | $e_{i}$
---|---|---|---|---|---|---
$\bar{I}$ | $\Lambda$ | 1.493 | 0.06409 | 0.02104 | $-5.381\times 10^{-4}$ | $1.957\times 10^{-6}$
$\bar{I}$ | $\bar{Q}$ | 1.387 | 0.5722 | 0.01043 | $0.02342$ | $6.245\times 10^{-4}$
$\bar{Q}$ | $\Lambda$ | 0.1899 | 0.09937 | 0.04380 | $-3.430\times 10^{-3}$ | $7.054\times 10^{-5}$

##  IV Summary

In summary, the RBHF theory in the full Dirac space has been employed to study
the gravitational redshift, moment of inertia, and quadrupole moment of
neutron star under the slow-rotation and small-tidal-deformation
approximation. The one-to-one correspondence relation for gravitational
redshift and mass is established and used to infer the masses of isolated
neutron stars by combining gravitational redshift measurements. The ratio of
the moment of inertia $I$ to $MR^{2}$ as a function of the compactness $M/R$
is obtained, and is consistent with the universal relations shown by Lattimer
et al. Lattimer and Schutz (2005) and that from Bayesian posterior probability
distributions by Lim et al. Lim  _et al._ (2019). Using $NN$ interactions Bonn
A, B, and C, the moment of inertia for $1.338M_{\odot}$ pulsar PSR J0737-3039A
is predicted to be $1.356\times 10^{45}$, $1.381\times 10^{45}$, and
$1.407\times 10^{45}\ \mathrm{g~{}cm^{2}}$, which are consistent with the
constraint translated from the tidal deformability deduced from GW170817 with
universal relations among neutron star observables. The EOSs constructed by
the RBHF theory in the full Dirac space, together with those by the projection
method and momentum-independence approximation, are compared successfully to
universal $I$-Love-$Q$ relations. By combing the tidal deformability
$\Lambda_{1.4M_{\odot}}$ from GW170817 and the numerical fitting for these
universal relations from relativistic ab initio EOSs, the moment of inertia of
neutron star with 1.4 solar mass is deduced as
$I_{1.4M_{\odot}}=1.22^{+0.40}_{-0.25}\times 10^{45}\mathrm{g\ cm^{2}}$.

###### Acknowledgements.

We thank Yeunhwan Lim for providing the probability distribution data for the
moment of inertia from the Bayesian analysis and Armen Sedrakian for carefully
reading the manuscript. This work was partly supported by the National Natural
Science Foundation of China (NSFC) under Grant No. 12147102, the Fundamental
Research Funds for the Central Universities under Grants No. 2020CDJQY-Z003
and No. 2021CDJZYJH-003, and the MOST-RIKEN Joint Project “Ab initio
investigation in nuclear physics.” Part of this work was achieved by using the
supercomputer OCTOPUS at the Cybermedia Center, Osaka University under the
support of the Research Center for Nuclear Physics of Osaka University.

## References

  * Lattimer and Prakash (2004) J. M. Lattimer and M. Prakash, Science 304, 536 (2004).
  * Lattimer and Prakash (2007) J. M. Lattimer and M. Prakash, Physics Reports 442, 109 (2007).
  * Özel and Freire (2016) F. Özel and P. Freire, Annual Review of Astronomy and Astrophysics 54, 401 (2016).
  * Lattimer and Prakash (2016) J. M. Lattimer and M. Prakash, Physics Reports 621, 127 (2016).
  * Burgio  _et al._ (2021) G. Burgio, H.-J. Schulze, I. Vidaña, and J.-B. Wei, Progress in Particle and Nuclear Physics 120, 103879 (2021).
  * Demorest  _et al._ (2010) P. B. Demorest, T. Pennucci, S. M. Ransom, M. S. E. Roberts, and J. W. T. Hessels, Nature 467, 1081 (2010).
  * Fonseca  _et al._ (2016) E. Fonseca _et al._ , ApJ 832, 167 (2016).
  * Antoniadis  _et al._ (2013) J. Antoniadis _et al._ , Science 340, 1233232 (2013).
  * Cromartie  _et al._ (2020) H. T. Cromartie, E. Fonseca, S. M. Ransom, P. B. Demorest, Z. Arzoumanian, H. Blumer, P. R. Brook, M. E. DeCesar, T. Dolch, J. A. Ellis, R. D. Ferdman, E. C. Ferrara, N. Garver-Daniels, P. A. Gentile, M. L. Jones, M. T. Lam, D. R. Lorimer, R. S. Lynch, M. A. McLaughlin, C. Ng, D. J. Nice, T. T. Pennucci, R. Spiewak, I. H. Stairs, K. Stovall, J. K. Swiggum, and W. W. Zhu, Nature Astronomy 4, 72 (2020).
  * Fonseca  _et al._ (2021) E. Fonseca, H. T. Cromartie, T. T. Pennucci, P. S. Ray, A. Y. Kirichenko, S. M. Ransom, P. B. Demorest, I. H. Stairs, Z. Arzoumanian, L. Guillemot, A. Parthasarathy, M. Kerr, I. Cognard, P. T. Baker, H. Blumer, P. R. Brook, M. DeCesar, T. Dolch, F. A. Dong, E. C. Ferrara, W. Fiore, N. Garver-Daniels, D. C. Good, R. Jennings, M. L. Jones, V. M. Kaspi, M. T. Lam, D. R. Lorimer, J. Luo, A. McEwen, J. W. McKee, M. A. McLaughlin, N. McMann, B. W. Meyers, A. Naidu, C. Ng, D. J. Nice, N. Pol, H. A. Radovan, B. Shapiro-Albert, C. M. Tan, S. P. Tendulkar, J. K. Swiggum, H. M. Wahl, and W. W. Zhu, The Astrophysical Journal Letters 915, L12 (2021).
  * Riley  _et al._ (2019) T. E. Riley, A. L. Watts, S. Bogdanov, P. S. Ray, R. M. Ludlam, S. Guillot, Z. Arzoumanian, C. L. Baker, A. V. Bilous, D. Chakrabarty, K. C. Gendreau, A. K. Harding, W. C. G. Ho, J. M. Lattimer, S. M. Morsink, and T. E. Strohmayer, The Astrophysical Journal 887, L21 (2019).
  * Miller  _et al._ (2019) M. C. Miller, F. K. Lamb, A. J. Dittmann, S. Bogdanov, Z. Arzoumanian, K. C. Gendreau, S. Guillot, A. K. Harding, W. C. G. Ho, J. M. Lattimer, R. M. Ludlam, S. Mahmoodifar, S. M. Morsink, P. S. Ray, T. E. Strohmayer, K. S. Wood, T. Enoto, R. Foster, T. Okajima, G. Prigozhin, and Y. Soong, The Astrophysical Journal 887, L24 (2019).
  * Riley  _et al._ (2021) T. E. Riley, A. L. Watts, P. S. Ray, S. Bogdanov, S. Guillot, S. M. Morsink, A. V. Bilous, Z. Arzoumanian, D. Choudhury, J. S. Deneva, K. C. Gendreau, A. K. Harding, W. C. G. Ho, J. M. Lattimer, M. Loewenstein, R. M. Ludlam, C. B. Markwardt, T. Okajima, C. Prescod-Weinstein, R. A. Remillard, M. T. Wolff, E. Fonseca, H. T. Cromartie, M. Kerr, T. T. Pennucci, A. Parthasarathy, S. Ransom, I. Stairs, L. Guillemot, and I. Cognard, The Astrophysical Journal Letters 918, L27 (2021).
  * Miller  _et al._ (2021) M. C. Miller, F. K. Lamb, A. J. Dittmann, S. Bogdanov, Z. Arzoumanian, K. C. Gendreau, S. Guillot, W. C. G. Ho, J. M. Lattimer, M. Loewenstein, S. M. Morsink, P. S. Ray, M. T. Wolff, C. L. Baker, T. Cazeau, S. Manthripragada, C. B. Markwardt, T. Okajima, S. Pollard, I. Cognard, H. T. Cromartie, E. Fonseca, L. Guillemot, M. Kerr, A. Parthasarathy, T. T. Pennucci, S. Ransom, and I. Stairs, The Astrophysical Journal Letters 918, L28 (2021).
  * Abbott and et al. (2017) B. P. Abbott and et al. (LIGO Scientific Collaboration and VIRGO Collaboration), Phys. Rev. Lett. 119, 161101 (2017).
  * Damour  _et al._ (1992) T. Damour, M. Soffel, and C. Xu, Phys. Rev. D 45, 1017 (1992).
  * Hinderer (2008) T. Hinderer, Astrophys. J. 677, 1216 (2008).
  * Flanagan and Hinderer (2008) E. E. Flanagan and T. Hinderer, Phys. Rev. D 77, 021502 (2008).
  * Damour and Nagar (2009) T. Damour and A. Nagar, Phys. Rev. D 80, 084035 (2009).
  * Postnikov  _et al._ (2010) S. Postnikov, M. Prakash, and J. M. Lattimer, Phys. Rev. D 82, 024016 (2010).
  * Fattoyev  _et al._ (2018) F. J. Fattoyev, J. Piekarewicz, and C. J. Horowitz, Phys. Rev. Lett. 120, 172702 (2018).
  * Annala  _et al._ (2018) E. Annala, T. Gorda, A. Kurkela, and A. Vuorinen, Phys. Rev. Lett. 120, 172703 (2018).
  * Most  _et al._ (2018) E. R. Most, L. R. Weih, L. Rezzolla, and J. Schaffner-Bielich, Phys. Rev. Lett. 120, 261103 (2018).
  * Tews  _et al._ (2018) I. Tews, J. Margueron, and S. Reddy, Phys. Rev. C 98, 045804 (2018).
  * Malik  _et al._ (2018) T. Malik, N. Alam, M. Fortin, C. Providência, B. K. Agrawal, T. K. Jha, B. Kumar, and S. K. Patra, Phys. Rev. C 98, 035804 (2018).
  * Zhang  _et al._ (2018) N.-B. Zhang, B.-A. Li, and J. Xu, The Astrophysical Journal 859, 90 (2018).
  * Tong  _et al._ (2020) H. Tong, P. Zhao, and J. Meng, Phys. Rev. C 101, 035802 (2020).
  * Lyne  _et al._ (2004) A. G. Lyne, M. Burgay, M. Kramer, A. Possenti, R. Manchester, F. Camilo, M. A. McLaughlin, D. R. Lorimer, N. D’Amico, B. C. Joshi, J. Reynolds, and P. C. C. Freire, Science 303, 1153 (2004).
  * Lattimer and Schutz (2005) J. M. Lattimer and B. F. Schutz, The Astrophysical Journal 629, 979 (2005).
  * Yagi and Yunes (2013a) K. Yagi and N. Yunes, Science 341, 365 (2013a).
  * Burgay  _et al._ (2003) M. Burgay, N. D’Amico, A. Possenti, R. N. Manchester, A. G. Lyne, B. C. Joshi, M. A. McLaughlin, M. Kramer, J. M. Sarkissian, F. Camilo, V. Kalogera, C. Kim, and D. R. Lorimer, Nature 426, 531 (2003).
  * Kramer and Wex (2009) M. Kramer and N. Wex, Classical and Quantum Gravity 26, 073001 (2009).
  * Kramer  _et al._ (2021) M. Kramer, I. H. Stairs, R. N. Manchester, N. Wex, A. T. Deller, W. A. Coles, M. Ali, M. Burgay, F. Camilo, I. Cognard, T. Damour, G. Desvignes, R. D. Ferdman, P. C. C. Freire, S. Grondin, L. Guillemot, G. B. Hobbs, G. Janssen, R. Karuppusamy, D. R. Lorimer, A. G. Lyne, J. W. McKee, M. McLaughlin, L. E. Münch, B. B. P. Perera, N. Pol, A. Possenti, J. Sarkissian, B. W. Stappers, and G. Theureau, Phys. Rev. X 11, 041050 (2021).
  * Lattimer (2021) J. Lattimer, Annual Review of Nuclear and Particle Science 71, 433 (2021).
  * Lim  _et al._ (2019) Y. Lim, J. W. Holt, and R. J. Stahulak, Phys. Rev. C 100, 035802 (2019).
  * Li and Sedrakian (2019) J. J. Li and A. Sedrakian, Phys. Rev. C 100, 015809 (2019).
  * Tong  _et al._ (2022) H. Tong, C. Wang, and S. Wang, The Astrophysical Journal 930, 137 (2022).
  * Yagi and Yunes (2013b) K. Yagi and N. Yunes, Phys. Rev. D 88, 023009 (2013b).
  * Yagi and Yunes (2017) K. Yagi and N. Yunes, Physics Reports 681, 1 (2017).
  * Yagi  _et al._ (2014) K. Yagi, L. C. Stein, G. Pappas, N. Yunes, and T. A. Apostolatos, Phys. Rev. D 90, 063010 (2014).
  * Jiang and Yagi (2020) N. Jiang and K. Yagi, Phys. Rev. D 101, 124006 (2020).
  * Landry and Kumar (2018) P. Landry and B. Kumar, The Astrophysical Journal 868, L22 (2018).
  * Kumar and Landry (2019) B. Kumar and P. Landry, Phys. Rev. D 99, 123026 (2019).
  * Wei  _et al._ (2019) J.-B. Wei, A. Figura, G. F. Burgio, H. Chen, and H.-J. Schulze, Journal of Physics G: Nuclear and Particle Physics 46, 034001 (2019).
  * Müther  _et al._ (1987) H. Müther, M. Prakash, and T. Ainsworth, Phys. Lett. B 199, 469 (1987).
  * Brockmann and Machleidt (1990) R. Brockmann and R. Machleidt, Phys. Rev. C 42, 1965 (1990).
  * Li  _et al._ (1992) G. Q. Li, R. Machleidt, and R. Brockmann, Phys. Rev. C 45, 2782 (1992).
  * Engvik  _et al._ (1996) L. Engvik, E. Osnes, M. Hjorth-Jensen, G. Bao, and E. Ostgaard, The Astrophysical Journal 469, 794 (1996).
  * Gross-Boelting  _et al._ (1999) T. Gross-Boelting, C. Fuchs, and A. Faessler, Nucl. Phys. A 648, 105 (1999).
  * Katayama and Saito (2015) T. Katayama and K. Saito, Phys. Lett. B 747, 43 (2015).
  * Anastasio  _et al._ (1980) M. R. Anastasio, L. S. Celenza, and C. M. Shakin, Phys. Rev. Lett. 45, 2096 (1980).
  * Anastasio  _et al._ (1981) M. R. Anastasio, L. S. Celenza, and C. M. Shakin, Phys. Rev. C 23, 2273 (1981).
  * Nuppenau  _et al._ (1989) C. Nuppenau, Y. Lee, and A. MacKellar, Nucl. Phys. A 504, 839 (1989).
  * Horowitz and Serot (1984) C. Horowitz and B. D. Serot, Phys. Lett. B 137, 287 (1984).
  * Wang  _et al._ (2021) S. Wang, Q. Zhao, P. Ring, and J. Meng, Phys. Rev. C 103, 054319 (2021).
  * Wang  _et al._ (2022a) S. Wang, H. Tong, and C. Wang, Phys. Rev. C 105, 054309 (2022a).
  * Wang  _et al._ (2022b) S. Wang, H. Tong, Q. Zhao, C. Wang, P. Ring, and J. Meng, Phys. Rev. C 106, L021305 (2022b).
  * Machleidt (1989) R. Machleidt, Adv. Nucl. Phys. 19, 189 (1989).
  * Abbott and et al. (2018) B. P. Abbott and et al. (The LIGO Scientific Collaboration and the Virgo Collaboration), Phys. Rev. Lett. 121, 161101 (2018).
  * Serot and Walecka (1986) B. Serot and J. Walecka, Advance in Nuclear Physics 16, 1 (1986).
  * Krastev and Sammarruca (2006) P. G. Krastev and F. Sammarruca, Phys. Rev. C 74, 025808 (2006).
  * Baym  _et al._ (1971) G. Baym, H. A. Bethe, and C. J. Pethick, Nuclear Physics A 175, 225 (1971).
  * Baym  _et al._ (1971) G. Baym, C. Pethick, and P. Sutherland, Astrophys. J. 170, 299 (1971).
  * Rhoades and Ruffini (1974) C. E. Rhoades and R. Ruffini, Phys. Rev. Lett. 32, 324 (1974).
  * Gandolfi  _et al._ (2012) S. Gandolfi, J. Carlson, and S. Reddy, Phys. Rev. C 85, 032801(R) (2012).
  * Oppenheimer and Volkoff (1939) J. R. Oppenheimer and G. M. Volkoff, Phys. Rev. 55, 374 (1939).
  * Tolman (1939) R. C. Tolman, Phys. Rev. 55, 364 (1939).
  * Hinderer  _et al._ (2010) T. Hinderer, B. D. Lackey, R. N. Lang, and J. S. Read, Phys. Rev. D 81, 123016 (2010).
  * Hartle (1967) J. B. Hartle, Astrophys. J. 150, 1005 (1967).
  * Hartle and Thorne (1968) J. B. Hartle and K. S. Thorne, The Astrophysical Journal 153, 807 (1968).
  * Fattoyev and Piekarewicz (2010) F. J. Fattoyev and J. Piekarewicz, Phys. Rev. C 82, 025810 (2010).
  * Hambaryan  _et al._ (2014) V. Hambaryan, R. Neuhäuser, V. Suleimanov, and K. Werner, Journal of Physics: Conference Series 496, 012015 (2014).
  * Hambaryan, V.  _et al._ (2017) Hambaryan, V., Suleimanov, V., Haberl, F., Schwope, A. D., Neuhäuser, R., Hohle, M., and Werner, K., A&A 601, A108 (2017).
  * van Dalen  _et al._ (2004) E. N. E. van Dalen, C. Fuchs, and A. Faessler, Nucl. Phys. A 744, 227 (2004).
  * Morrison  _et al._ (2004) I. A. Morrison, T. W. Baumgarte, S. L. Shapiro, and V. R. Pandharipande, The Astrophysical Journal 617, L135 (2004).
  * Akmal  _et al._ (1998) A. Akmal, V. R. Pandharipande, and D. G. Ravenhall, Phys. Rev. C 58, 1804 (1998).
  * Wiringa  _et al._ (1995) R. B. Wiringa, V. G. J. Stoks, and R. Schiavilla, Phys. Rev. C 51, 38 (1995).
  * Pudliner  _et al._ (1995) B. S. Pudliner, V. R. Pandharipande, J. Carlson, and R. B. Wiringa, Phys. Rev. Lett. 74, 4396 (1995).


∎

11institutetext: J. Zhang 22institutetext: Z. Li 33institutetext: L. Ma
44institutetext: G. Zheng 55institutetext: L. Chen(🖂) 66institutetext: School
of Physics and Information Technology, Shaanxi Normal University, Xi’an
710062, P. R. China
66email<EMAIL_ADDRESS>77institutetext: J. Zhang 88institutetext: School
of Physics and Electronic-Electrical Engineering, Ningxia University, Yinchuan
750021, P. R. China

# Oscillatory cooperation prevalence emerges from misperception

Jing Zhang  Zhao Li  Jiqiang Zhang  Lin Ma  Guozhong Zheng  Li Chen

(Received: date / Accepted: date)

###### Abstract

Oscillatory behaviors are ubiquitous in nature and the human society. However,
most previous works fail to reproduce them in the two-strategy game-
theoretical models. Here we show that oscillatory behaviors naturally emerge
if incomplete information is incorporated into the cooperation evolution of a
non-Markov model. Specifically, we consider a population playing prisoner’s
dilemma game, where each individual can only probabilistically get access to
their neighbors’ payoff information and store them within their memory with a
given length. They make their decisions based upon these memories.
Interestingly, we find that the level of cooperation generally cannot
stabilize but render quasi-periodic oscillation, and this observation is
strengthened for a longer memory and a smaller information acquisition
probability. The mechanism uncovered shows that there are misperceived payoffs
about the player’s neighborhood, facilitating the growth of cooperators and
defectors at different stages that leads to oscillatory behaviors as a result.
Our findings are robust to the underlying structure of the population. Given
the omnipresence of incomplete information, our findings may provide a
plausible explanation for the phenomenon of oscillatory behaviors in the real
world.

###### Keywords:

Cooperation Evolutionary game theory Prisoner’s dilemma Oscillation

##  1 Introduction

Human civilization is the result of large-scale cooperation in every aspect of
socioeconomic activities, ranging from raising offsprings to the division of
labor Axelrod and Hamilton (1981); Smith and Szathmary (1997). Altruistic
cooperation behaviors are also widely observed in different ecosystems in the
wild Clutton-Brock (2009). However, according to Darwin’s _The Origin of
Species_ , “survival of the fittest” implies the innate selfish nature of
individuals. This means that cooperative behaviors are actually
counterintuitive, because cooperators helping others is at a certain cost,
thus put themselves in an inferior competitive position to those free-riders
Rider  _et al._ (1984); Smith (1982); Gintis  _et al._ (2000); Nowak (2006a).
As a result, deciphering the mechanism of cooperation becomes a fundamental
challenge, which has attracted many researchers’ attention from different
fields and now is a highly interdisciplinary field Gintis  _et al._ (2000).

To proceed, the evolutionary game theory Nowak and Sigmund (2004) has been
introduced and is proved to be a useful theoretical framework. By analysing
stylized social dilemmas (e.g. the prisoner’s dilemma, the snowdrift game, and
the public goods game), several mechanisms have been proposed Nowak (2006a),
such as direct Trivers (1971) and indirect reciprocity Zahavi and Zahavi
(1999), kin Hamilton (1964) and group selection Keller (1999); Queller (1964),
reward and punishment Nowak (2006b), social diversity Sigmund  _et al._ (2001)
and social hierarchy Liang  _et al._ (2021), spatial or network reciprocity
Nowak and May (1992). Particularly, theoretically considering the fact that
human populations are structured with fixed neighbors can support cooperation.
The argument behind is that cooperator clusters are more likely to form in the
structured population, which potentially are able to defeat the invasion of
defectors, compared to the well-mixed scenario. However, human behavioral
experiments do not support this theoretic prediction in general Traulsen  _et
al._ (2010). This unsatisfactory situation suggests that the evolution of
cooperation in realistic population could be far more complex than what
current game-theoretic models captured, and the experiment-driven modeling
approach is needed. It’s worthy noting that the dynamical reciprocity is
recently also proposed as a counterpart mechanism Liang  _et al._ (2022),
showing that the game-game interactions could also improve the cooperation
prevalence, though remains to be validated in the experiments.

While most of previous works aim to understand the emergence of cooperation
where the prevalence finally settles down, there are quite some systems where
the oscillatory behaviors are often seen. The most prominent example is the
population oscillation of different species in many ecological systems, which
is vital to preserving biodiversity Kerr  _et al._ (2002); Reichenbach  _et
al._ (2007). An early work Nowak and Sigmund (1993) shows that in a population
consisting of some competing strategies, strategic mutation can lead to
persistent periodic or even chaotic oscillations in the frequencies of the
strategies and the level of cooperation. Although the model is formulated
within the iterated prisoner’s dilemma framework, the strategies are defined
by a quadruple of parameters, where there are as much as 16 species are
engaged in the evolution. Later, a mainstream explanation for oscillation is
via the cyclic dominance Sinervo and Lively (1996); Szabó and Hauert (2002);
Kerr  _et al._ (2002); Reichenbach  _et al._ (2007); Imhof  _et al._ (2005);
Szolnoki  _et al._ (2014). A popular example is the Rock-Paper-Scissors game
Zhou (2016), where the fractions of three strategies naturally oscillate given
the population is spatially structured Szolnoki  _et al._ (2010) or they are
capable of mobility Szabó and Hauert (2002); Kerr  _et al._ (2002);
Reichenbach  _et al._ (2007). However, this sort of models require three or
more species and thus cannot explain how the oscillatory behaviors emerge in
the 2-specie populations, such as the classic 10-year population cycle of
snowshoe hares and Canada lynx in the boreal forests of North America Elton
and Nicholson (1942); Moran (1953); Krebs  _et al._ (2001). Recently, Yang et
al. Yang and Li (2022) observed a tide-like burst in a probabilisitic
migration model of prisoner’s dilemma, where migration is driven by conformity
and self-centered inequity aversion norms. It’s also worth noting that
recently Zhang et al. Zhang  _et al._ (2020) introduced the reinforcement
learning framework, where periodic burst-like oscillation is also seen in
pairwise games.

Here we provide a different explanation for oscillations in the pairwise game
within the classic framework, where we incorporate two commonly seen
ingredients: the non-Markov effect and incomplete information. While the
former considers the historical impact, the latter is due to the limited
information acquisition capability, because their neighbors’ payoff
information is not always accessible. In fact, the argument behind the non-
markov effect is well-acknowledged that the decision-making for individuals is
based on the knowledge of the past records rather than a single round Szabó
and Fath (2007). Plenty of theoretical works show that memory can considerably
improve the cooperation level Wang  _et al._ (2006); Qin  _et al._ (2008);
Yong-Kui  _et al._ (2010); Jiang  _et al._ (2009), and several related factors
such as the length of memory, the temptation of defection, and the structure
of networks have been studied Luo  _et al._ (2016); Yong-Kui  _et al._ (2010);
Yang  _et al._ (2012); Stewart and Plotkin (2016); Liu  _et al._ (2020);
Horvath  _et al._ (2012); Shu  _et al._ (2019); Ye  _et al._ (2017); Wang  _et
al._ (2020). The impact of incomplete information scenario, where players do
not possess full information about their opponents, has also been investigated
separately in previous studies. Ref. Vuolevi and Van Lange (2012) shows
experimentally that incompleteness of information about the opponent’s past
behaviors undermines both expectations about another person’s cooperation as
well as one’s own cooperation. Similar observations are made that selfish
behaviors are the best way of protecting oneself against non-cooperative
behaviors in the presence of incomplete information Tazdaït  _et al._ (2008).
However, in a finitely repeated PD, incomplete information about their
partners’ preference is found to trigger the possibility of coopepration
Dijkstra and van Assen (2017).

In this work, we find that when the two ingredients are incorporated
simultaneously, oscillatory behaviors emerge, where the cooperation prevalence
rises and falls in a quasi-periodic way. Specifically, we consider a
structured population of individuals where they can only probabilistically get
access to their neighbors’ payoffs and store this information in their memory;
by comparing the perceived payoffs to their own values they decide whether to
imitate their neighbours’ cooperation propensity or not. We find oscillation
is likely to occur when the memory is long and the payoff information is hard
to obtain. Further analysis shows that the oscillation is an intrinsic
property of the system, the presence of incomplete information leads to a
misperception of neighbors’ payoff that causes cooperation explosion and
decline. The longer memory and stronger incomplete information are, the
stronger perception errors are expected.

This paper is organized as follows. We introduce our model in Sec. 2.
Numerical simulation results on the $2d$ square lattice are shown in Sec. 3.
The mechanism is analysed in Sec. 4. Finally, we conclude and discuss our work
in Sec. 5.

Figure 1:  Typical time series of cooperation level $f_{C}$ for four different
combination of the memory length $M$ and information acquisition probability
$p$: (a) $M=1$, $p=1.0$; (b) $M=100$, $p=1.0$; (c) $M=20$, $p=0.2$; and (d)
$M=100$, $p=0.2$. Parameter: $b=1.2$.

##  2 Model

We study the prisoner’s dilemma (PD) game with $N$ individuals that are
located on an $L\times L$ square lattice with a periodic boundary condition.
PD is a typical pairwise game for many social dilemmas, with the strategy
being either cooperation (C) or defection (D). Mutual cooperation brings the
reward $R$, mutual defection yields the punishment $P$ for each, and mixed
encounter gives the cooperator the sucker’s payoff $S$, yet the temptation $T$
for the defector. The conditions of $T>R>P>S$ and $2R>T+S$ are required for
PD. It’s easy to see that mutual defection is the Nash equilibrium, even
though the mutual cooperation is optimal for their collective profits. We
adopt the weak version of PD Nowak and May (1992) in this work, $R=1$,
$P=S=0$, and $T=b$, where $1.0\leq b\leq 2.0$.

To incorporate the history, each player (e.g. player $i$) has $k_{i}$ memories
$\mathcal{M}_{i}(j)=\\{\Pi_{i}^{1}(j),...,\Pi_{i}^{M}(j)\\}$ of length $M$ to
record the payoff information for each neighbor $j\in\Omega_{i}$, where
$\Omega_{i}$ denote the neighborhood of player $i$, and $k_{i}=|\Omega_{i}|$
the degree of node $i$ and is four in the 2d square lattice. Due to the
uncertainties in realistic surroundings, this information is not always
accessible in each round. We assume that the payoff information in the
neighborhood is obtained with an _acquisition probability_ $p$. The larger the
value of $p$, the more likely to obtain their neighbors’ payoff information,
and vice versa. The case of $p=1$ recovers to the complete information
scenario that many previous work assumed Szabó and Fath (2007); Liu  _et al._
(2010).

Our model follows Monte Carlo (MC) simulation procedure. At the very start,
each player’s strategy is randomly assigned with either $C$ or $D$. The
implementation of an elementary MC step is as follows. Firstly, we randomly
select a player $i$ and one of its neighbor player $j$, and calculate their
total payoffs $\Pi_{i,j}$ by playing with all their neighbors. Secondly, the
obtained payoff of player $j$ is then recorded in the memory of player $i$
with the probability $p$; if succeed $\Pi_{i}^{1}(j)=\Pi_{j}$, while the rest
are shifted by one site with the earliest payoff information $\Pi_{i}^{M}(j)$
being removed accordingly; otherwise the memory $\mathcal{M}_{i}(j)$ will not
be updated.

Next, the probability for player $i$ to adjust its strategy $s_{i}$ is given
by the Fermi rule Szabó and Tőke (1998)

$\displaystyle W(\rho_{j}^{C}\rightarrow
s_{i})=\dfrac{1}{1+\exp[(\overline{\Pi}_{i}-\overline{\Pi}_{j}^{\prime})/K]},$
(1)

where $\overline{\Pi}_{i}$ is the averaged payoffs for player $i$ in the
latest $M$ rounds since player $i$ knows exactly her/his own payoffs, and
$\overline{\Pi}_{j}^{\prime}=\sum_{m=1}^{M}\Pi_{i}^{m}(j)/M$ is the average
payoffs of $j$ according to the memory $\mathcal{M}_{i}(j)$. Their subtraction
in Eq. (1) can then be interpreted as the perceived payoff difference by
player $i$. $K$ qualifies the surrounding noise during the imitation process,
which is usually interpreted as the bounded rationality. $K\to 0$ indicates
that players are absolutely rational and the imitation is deterministic;
$K\to\infty$, on the other hand, means that the decision-making is completely
random and is not affected by the neighbors. Here $K=0.1$ throughout the work
Szabó and Tőke (1998). Instead of directly copying $s_{j}$ for imitation,
player $i$ adopts $j$’s cooperation propensity defined as
${\rho}_{j}^{C}=N^{C}_{j}/M$, where $N^{C}_{j}$ is the number of times acting
as a cooperator for player $j$ in its latest $M$ rounds. Different from the
payoff information, the neighbours’ strategy $s_{j}$ is always accessible for
players $i$ because the two interact directly. Finally, player $i$ behaves as
a cooperator with the probability ${\rho}_{j}^{C}$ once the imitation is
successful, otherwise no strategy change is made.

Notice that the above setup follows a typical asynchronous updating procedure.
A complete Monte Carlo step (MCS) consists of $N$ elementary steps, meaning
that every player updates its state exactly once on average. We compute the
cooperation prevalence $f_{C}=\frac{1}{N}\sum_{i=1}^{N}s_{i}$ as our primary
order parameter, measuring the overall preference in cooperation of the
population. If not stated otherwise, $L=100$ throughout our study, though the
impact of size will be also discussed in the latter part.

Figure 2:  Color-coded phase diagram in the $M-p$ space. (a) The difference
between the average maximal and the minimal cooperation level $\langle
f_{C}^{max}\rangle-\langle f_{C}^{min}\rangle$; (b) The period corresponding
to the principal frequency ($T=1/\omega_{max}$) by Fourier analysis of the
time series for the given parameter combination of $M$ and $p$. Other
parameter: $b=1.2$.

##  3 Numerical results

Figure 1 reports typical time series for several combinations of information
acquisition probability $p$ and memory length $M$ by fixing the temptation
$b=1.2$. Let’s first see the extreme case where the memory is absent ($M=1$
and $p=1$, Fig. 1(a)), which recovers the traditional Markov model Szabó  _et
al._ (2005). As can be seen, no cooperation is seen for the given parameter
$b$, in line with the previous result. When the non-Markov effect is
incorporated but with full information acquisition ($M=100$ and $p=1$),
cooperation arises and the fraction is above $80\%$, see Fig. 1(b). This means
that by taking the history into account in the decision-making, the
cooperation prevalence can be promoted, which is also reported in previous
studies Liu  _et al._ (2010). But once the incomplete information is
incorporated, the cooperation prevalence becomes unstable and its average
value is slightly decreased, see Fig. 1(c,d). When the memory is short ($M=20$
and $p=0.2$, Fig. 1(c)), the cooperation prevalence $f_{C}$ shows noise-like
fluctuations, while a longer memory ($M=100$ and $p=0.2$) gives rise to
oscillatory behaviors, as shown in Fig. 1(d). A closer lookup shows this
oscillation is quasi-periodic, confirmed by the pronounced peaks by doing
Fourier analysis (see Sec. I in SM SM ).

To further investigate the properties of the cooperation oscillation, we
compute the difference between the average maximal prevalence at the peak
$\langle f_{C}^{max}\rangle$ and the minimal prevalence at the valley $\langle
f_{C}^{min}\rangle$ as a function of the memory length $M$ and the acquisition
probability $p$. The obtained phase diagram Fig. 2(a) confirms that
oscillatory behaviors are more likely to occur for a longer memory and a small
value of $p$, corresponding to strong non-Markov effect and severe incomplete
information scenarios. Fig. 2(b) plots the phase diagram of corresponding
period $T=1/\omega_{max}$, where $\omega_{max}$ is the principal frequency in
the power spectrum. It shows that the strongly non-Markov effect and severely
incomplete information scenarios give rise to a long period of oscillation.

To develop some intuition, typical spatial patterns are provided in Fig. 3
with the same parameters used as in Fig. 1(d). We can see that in Fig. 3(a),
where the cooperation level $f_{c}$ is at its minimum in time, some cooperator
clusters emerge and start to grow. Afterwards, these clusters are merged with
each other [Fig. 3(b)], reaching its maximum at some time point [Fig. 3(c)].
The later snapshots show that defectors increase but in a scattered manner
[Fig. 3(d)], and the cooperation prevalence $f_{C}$ declines to its minimum
[Fig. 3(e)]. Afterwards, the cooperator clusters re-emerge and grow [Fig.
3(f)], the above process repeats again and again.

For completeness, we also provide the corresponding patterns for the
cooperation propensity, the average cooperation level over the last $M$
rounds, see Fig. 4. Compared to the scattered patterns shown in Fig. 3, these
plots show that the cooperation propensity varies smoothly across the domain,
though the two figures are highly correlated. Detailed examination of the
spatial patterns by compactness analysis shows that the cooperators are much
closely clustered than the defectors, both of their compactness levels
oscillate with the same frequency of $f_{c}$, but is not in perfect synchrony
(see Sec. II in SM SM ).

Figure 3:  Cooperation state patterns at different times, where yellow and
blue sites denote cooperators and defectors, respectively. From (a) to (c),
cooperation is in the explosive stage, followed by the cooperation decline
(d,e), and then rising again (f). Parameters: $M=100$, $p=0.2$, and $b=1.2$.

##  4 Mechanism analysis

To understand the formation of cooperation oscillation, we now turn to its
dynamical mechanism analysis. Without loss of generality, we focus on the case
with the parameters used in Fig. 1(d) in this section.

Figure 4:  The corresponding cooperation propensity patterns defined as
${\rho}_{j}^{C}=N^{C}_{j}/M$. Other settings are exactly the same as Fig. (3).
Figure 5:  The time evolution of the average effective payoff difference for
four different strategy combinations for $i-j$ pair, where $i$ is the focal
node and $j$ is one of $i$’s neighbors. The inset shows an enlarged evolution
section illustrating two phases: I – cooperation explosion and II –
cooperation decline. Parameters: $M=100$, $b=1.2$, and $p=0.2$.

Specifically, we monitor the effective payoff differences for different pairs
to understand how they are correlated to the strategy update, and thus also
the evolution of cooperation prevalence $f_{C}$, see Fig. 5. For example,
$\Delta\Pi_{CC}=\overline{\Pi}_{i}-\overline{\Pi}_{j}^{\prime}$ denotes the
difference between the average payoff of the focal player $i$ and the
perceived payoff of its neighbor $j$, where both players are cooperators. By
combination, there are four scenarios in total. A critical observation in Fig.
5 is that all four $\Delta\Pi$ are oscillating. According to the Eq. (1), this
implies that the probabilities to imitate are time-varying, their profiles are
found to synchronously fluctuate with the derivatives of the cooperation
prevalence $f^{\prime}$ (see the inset in Fig. 5). Note that, different from
the traditional case of directly imitating their neighbor strategy where C-C
or D-D pairs remain unchanged by the mutual imitation, the updating in our
case could still lead to the strategy change due to the probabilistic nature
in the imitation. That’s why we also need to monitor the cases of C-C and D-D
pairs.

To see how the cooperation oscillates, we divide the evolution into two
phases, see the inset in Fig. 5.

_Phase I — cooperation explosion._ This stage corresponds to Fig. 3(a-c) and
Fig. 4(a-c), where the cooperation prevalence is rising. The key feature is
that the values of the four $\Delta\Pi$s are larger than their time averages,
this is especially true for $\Delta\Pi_{CC}$. This means that the strategy
updating in this stage is less frequent since the probability
$W(\rho_{j}^{C}\rightarrow s_{i})$ is relatively small. This is important for
cooperation clusters to keep their high cooperation prevalence $\rho_{C}$.
When cooperator clusters are formed, they are at an advantageous position over
defectors, for the same logic behind the network reciprocity Nowak and May
(1992); therefore, the cooperation clusters expand and the low cooperation
region shrinks as shown from Fig. 3(a) to 3(c) and Fig. 4(a) to 4(c).

_Phase II — cooperation decline_. Once the domain is dominated by the high
cooperation prevalence [i.e. Fig. 3(c) and Fig. 4(c)], shifting to defection
would be a profitable move for individuals in such an environment. That’s what
happens from Fig. 3(c)-3(e) and Fig. 4(c)-4(e), where the cooperation
prevalence declines. In this phase, all four $\Delta\Pi$ are lower than their
time averages, meaning that the updating events occur very frequently.

At the end of phase II, the domain is now dominated by defectors, defectors
are now at the disadvantage when they meet up with cooperation clusters. As a
result, the cooperation decline stops, the regions with low $\rho_{C}$ are
invaded again by cooperation clusters, and the dynamical process of Phase I
restarts, and repeats again and again.

In both phases, the presence of incomplete information (i.e. the acquisition
probability $p<1$) facilitates the cooperation explosion and decline,
respectively. In Phase I, a small $p$ inhibits the spread of defection within
the cooperation clusters. Consider the scenario where a player (say $j$)
defects within a highly cooperative neighborhood, this brings $j$ a high
payoff as a result, but due to the presence memory and the acquisition
probability, this payoff information is not immediately reflected in the
payoff difference perceived by its neighbors. The payoff advantage of
defection needs a longer time to spread within the cooperation clusters; this
leaves a time window for cooperation cluster’s growth till a relatively high
level of cooperation. In Phase II, our scenario facilitates the spread of
defection on the contrary. Consider the invasion of defection (say player $j$)
into regions with a relatively high cooperation [say player $i$, c.f. Fig.
4(d)], the invasion is not possible in the absence of incomplete information
because defectors have lower payoffs when they are clustered, showing no
advantage in this scenario. But what player $j$ perceived at this time point
is higher $\overline{\Pi}_{j}^{\prime}>\overline{\Pi}_{j}$ because
$\overline{\Pi}_{j}^{\prime}$ includes early decent payoffs when $j$ was also
a cooperator, this helps its invasion and leads to the extinction of
relatively high cooperation regions. In both phases, the degree of cooperation
explosion and decline varies due to the stochastic nature, the amplitude of
oscillation is fluctuating, as we seen in Fig. 1(d).

In brief, due to the presence of incomplete information ($p<1$), there is a
perception mismatch that inhibits the invasion of defection in the first phase
and facilitates its invasion in the second phase. That leads to a cooperation
explosion and followed by a cooperation decline, and the process repeats,
oscillation is thus formed. The longer the memory, the stronger mismatch we
expected. That explains the oscillatory behaviors are only observed in the
case of a large $M$ and a small $p$. As a comparison, when the incomplete
information is removed from the model ($p=1$), the memory effect alone is
insufficient to trigger the cooperation oscillation, the payoff difference is
relatively stable without perception mismatch, see Fig. 6. Their spatial
patterns show no cooperation explosion and decline after transient evolution,
see see Sec. III in SM SM .

Figure 6:  The time evolution of the average effective payoff difference for
four different strategy combinations for $i-j$ pair, where $i$ is the focal
node and $j$ is one of $i$’s neighbors. Parameters: $M=100$, $b=1.2$, and
$p=1.0$.

##  5 Discussion and Conclusion

In summary, our work extends the previous game-theoretic model of cooperation,
by allowing for both the memory effect and the incomplete information. While
the former factor is well-motivated and has been extensively studied
previously Wang  _et al._ (2006); Qin  _et al._ (2008); Yong-Kui  _et al._
(2010); Jiang  _et al._ (2009); Szabó and Fath (2007), the latter considers
the fact that players in the real world are not always aware of their
neighbors’ payoff information, not the complete information scenario assumed
as in most previous work Szabó and Fath (2007). The inclusion of memory leads
to the promotion of cooperation prevalence as also revealed previously, we
uncover here that further inclusion of the incomplete information de-
stabilizes the lifted prevalence and yields a quasi-periodic oscillation. The
oscillation is strengthened for a long memory and a strong incomplete
information situation. Further analysis shows that the oscillation is an
intrinsic property when both ingredients are included. They together lead to a
misperception of individual neighborhood’s payoff, erroneously making
decisions. The consequence is that the invasion of defection is inhibited and
facilitated in the phase of cooperation explosion and decline, respectively,
causing cooperation oscillation as a result.

Since the oscillation is an intrinsic property of the system, and no synchrony
is found in the prevalence for different locations, we expect that the
oscillation amplitude of the cooperation prevalence $f_{C}$ becomes smaller as
the system size becomes larger. This is confirmed in our further experiments
(see Sec. IV in SM SM ). From the perspective of statistical physics, the
observation of oscillation in $f_{C}$ is a finite-size effect, and cannot be
expected in an infinitely large system. However, if we focus on a local region
of the population, oscillation is still observable. Since the populations in
the real world are always of finite size, the oscillatory behaviors should
still be observed. When the population is structured in complex networks,
numerical experiments show qualitatively the same phenomena (data not shown).
However, once the game is replaced with the snowdrift game, no oscillation is
seen, this is because the failure of spatial reciprocity in this game, where
cooperators are not well-clustered for their changes Hauert and Doebeli
(2004).

Compared to the existing game-theoretic models for explaining oscillatory
behaviors, our pairwise game model does not require three or more species
Szolnoki  _et al._ (2014), or seek for other complex mechanisms such as
mobility or conformity Reichenbach  _et al._ (2007); Yang and Li (2022). In
this sense, our model provides a relatively simple framework to understand the
emergence of oscillation. Given the omnipresence of incomplete information
scenario and the non-Markov process in the realistic games, our model may
provide a plausible explanation for a range of oscillatory phenomenon.

###### Acknowledgements.

We are supported by the Natural Science Foundation of China under Grants Nos.
12075144 and 61703257. ZJQ is supported by the Natural Science Foundation of
China under Grants No. 12165014.

### Data availability

The datasets generated during the current study are available from the
corresponding author on reasonable request.

## Declarations

### Conflict of interest

The authors declare that they have no conflict of interest.

## References

  * Axelrod and Hamilton (1981) R. Axelrod and W. D. Hamilton, Science 211, 1390 (1981).
  * Smith and Szathmary (1997) J. M. Smith and E. Szathmary,  _The major transitions in evolution_ (OUP Oxford, 1997).
  * Clutton-Brock (2009) T. Clutton-Brock, Nature 462, 51 (2009).
  * Rider  _et al._ (1984) R. Rider _et al._ , Journal of Economic Behavior & Organization 5, 406 (1984).
  * Smith (1982) J. M. Smith,  _Evolution and the Theory of Games_ (Cambridge university press, 1982).
  * Gintis  _et al._ (2000) H. Gintis _et al._ ,  _Game theory evolving: A problem-centered introduction to modeling strategic behavior_ (Princeton university press, 2000).
  * Nowak (2006a) M. A. Nowak,  _Evolutionary dynamics: exploring the equations of life_ (Harvard university press, 2006).
  * Nowak and Sigmund (2004) M. A. Nowak and K. Sigmund, Science 303, 793 (2004).
  * Trivers (1971) R. L. Trivers, The Quarterly review of biology 46, 35 (1971).
  * Zahavi and Zahavi (1999) A. Zahavi and A. Zahavi,  _The handicap principle: A missing piece of Darwin’s puzzle_ (Oxford University Press, 1999).
  * Hamilton (1964) W. D. Hamilton, Journal of Theoretical Biology 7, 17 (1964).
  * Keller (1999) L. Keller,  _Levels of selection in evolution_ , Vol. 21 (Princeton University Press, 1999).
  * Queller (1964) D. C. Queller, Nature 201, 1145 (1964).
  * Nowak (2006b) M. A. Nowak, Science 314, 1560 (2006b).
  * Sigmund  _et al._ (2001) K. Sigmund, C. Hauert, and M. A. Nowak, Proceedings of the National Academy of Sciences 98, 10757 (2001).
  * Liang  _et al._ (2021) R. Liang, J. Zhang, G. Zheng, and L. Chen, Physica A: Statistical Mechanics and its Applications 567, 125726 (2021).
  * Nowak and May (1992) M. A. Nowak and R. M. May, Nature 359, 826 (1992).
  * Traulsen  _et al._ (2010) A. Traulsen, D. Semmann, R. D. Sommerfeld, H.-J. Krambeck, and M. Milinski, Proceedings of the National Academy of Sciences 107, 2962 (2010).
  * Liang  _et al._ (2022) R. Liang, Q. Wang, J. Zhang, G. Zheng, L. Ma, and L. Chen, Physical Review E 105, 054302 (2022).
  * Kerr  _et al._ (2002) B. Kerr, M. A. Riley, M. W. Feldman, and B. J. Bohannan, Nature 418, 171 (2002).
  * Reichenbach  _et al._ (2007) T. Reichenbach, M. Mobilia, and E. Frey, Nature 448, 1046 (2007).
  * Nowak and Sigmund (1993) M. Nowak and K. Sigmund, Proceedings of the National Academy of Sciences 90, 5091 (1993).
  * Sinervo and Lively (1996) B. Sinervo and C. M. Lively, Nature 380, 240 (1996).
  * Szabó and Hauert (2002) G. Szabó and C. Hauert, Physical Review Letters 89, 118101 (2002).
  * Imhof  _et al._ (2005) L. A. Imhof, D. Fudenberg, and M. A. Nowak, Proceedings of the National Academy of Sciences 102, 10797 (2005).
  * Szolnoki  _et al._ (2014) A. Szolnoki, M. Mobilia, L.-L. Jiang, B. Szczesny, A. M. Rucklidge, and M. Perc, Journal of the Royal Society Interface 11, 20140735 (2014).
  * Zhou (2016) H.-J. Zhou, Contemporary Physics 57, 151 (2016).
  * Szolnoki  _et al._ (2010) A. Szolnoki, Z. Wang, J. Wang, and X. Zhu, Physical Review E 82, 036110 (2010).
  * Elton and Nicholson (1942) C. Elton and M. Nicholson, The Journal of Animal Ecology , 215 (1942).
  * Moran (1953) P. A. Moran, Australian Journal of Zoology 1, 291 (1953).
  * Krebs  _et al._ (2001) C. J. Krebs, R. Boonstra, S. Boutin, and A. R. Sinclair, BioScience 51, 25 (2001).
  * Yang and Li (2022) Z. Yang and Z. Li, Nonlinear Dynamics , 1 (2022).
  * Zhang  _et al._ (2020) S.-P. Zhang, J.-Q. Zhang, L. Chen, and X.-D. Liu, Nonlinear Dynamics 99, 3301 (2020).
  * Szabó and Fath (2007) G. Szabó and G. Fath, Physics Reports 446, 97 (2007).
  * Wang  _et al._ (2006) W.-X. Wang, J. Ren, G. Chen, and B.-H. Wang, Physical Review E 74, 056113 (2006).
  * Qin  _et al._ (2008) S.-M. Qin, Y. Chen, X.-Y. Zhao, and J. Shi, Physical Review E 78, 041129 (2008).
  * Yong-Kui  _et al._ (2010) L. Yong-Kui, L. Zhi, C. Xiao-Jie, and W. Long, Chinese Physics B 19, 090203 (2010).
  * Jiang  _et al._ (2009) L. L. Jiang, M. Zhao, H. X. Yang, J. Wakeling, B. H. Wang, and T. Zhou, Physical Review E 80, 31144 (2009).
  * Luo  _et al._ (2016) C. Luo, X. Zhang, H. Liu, and R. Shao, Physica A: Statistical Mechanics and its Applications 450, 560 (2016).
  * Yang  _et al._ (2012) H. X. Yang, Z. Rong, P. M. Lu, and Y. Z. Zeng, Physica A Statistical Mechanics and Its Applications 391, 4043 (2012).
  * Stewart and Plotkin (2016) A. J. Stewart and J. B. Plotkin, Scientific reports 6, 1 (2016).
  * Liu  _et al._ (2020) C. Liu, J. Wang, X. Li, and C. Xia, Physica A: Statistical Mechanics and its Applications 554, 124691 (2020).
  * Horvath  _et al._ (2012) G. Horvath, J. Kovářík, and F. Mengel, Journal of Theoretical Biology 300, 193 (2012).
  * Shu  _et al._ (2019) F. Shu, M. Li, and X. Liu, Chaos, Solitons & Fractals 120, 17 (2019).
  * Ye  _et al._ (2017) W. Ye, W. Feng, C. Lü, and S. Fan, Applied Mathematics and Computation 307, 31 (2017).
  * Wang  _et al._ (2020) J. Wang, R. Wang, F. Yu, Z. Wang, and Q. Li, Applied Mathematics and Computation 370, 124887 (2020).
  * Vuolevi and Van Lange (2012) J. H. Vuolevi and P. A. Van Lange, Acta Psychologica 141, 67 (2012).
  * Tazdaït  _et al._ (2008) T. Tazdaït, A. Caparros, and J.-C. Péreau, Rationality and Society (2008).
  * Dijkstra and van Assen (2017) J. Dijkstra and M. A. van Assen, The Journal of Mathematical Sociology 41, 1 (2017).
  * Liu  _et al._ (2010) Y. Liu, Z. Li, X. Chen, and L. Wang, Physica A: Statistical Mechanics and its Applications 389, 2390 (2010).
  * Szabó and Tőke (1998) G. Szabó and C. Tőke, Physical Review E 58, 69 (1998).
  * Szabó  _et al._ (2005) G. Szabó, J. Vukov, and A. Szolnoki, Physical Review E 72, 047107 (2005).
  * (53) See Supplemental Material for additional numerical results.
  * Hauert and Doebeli (2004) C. Hauert and M. Doebeli, Nature 428, 643 (2004).