# Chiral symmetry restoration and the $\Delta$ matter formation in neutron
stars

Michał Marczenko<EMAIL_ADDRESS>Incubator of Scientific
Excellence - Centre for Simulations of Superdense Fluids, University of
Wrocław, plac Maksa Borna 9, PL-50204 Wrocław, Poland  Krzysztof Redlich
Institute of Theoretical Physics, University of Wrocław, plac Maksa Borna 9,
PL-50204 Wrocław, Poland  Chihiro Sasaki  Institute of Theoretical Physics,
University of Wrocław, plac Maksa Borna 9, PL-50204 Wrocław, Poland

###### Abstract

We analyze the effects of chiral symmetry restoration in hadronic matter,
including the lowest-lying baryonic resonance $\Delta$ based on the parity
doublet model. We study the role of $\Delta$ and its chiral partner on the
equation of state (EoS) of dense matter under neutron star (NS) conditions of
$\beta$-equilibrium and charge neutrality. We find that the softening of the
EoS driven by the early onset of $\Delta$ matter due to partial restoration of
chiral symmetry allows accommodating the modern multi-messenger astrophysical
constraints on the mass, radius, and tidal deformability. The softening above
the saturation density is accompanied by subsequent stiffening at high
densities. We also find that the matter composition in the NS cores may be
different upon variations of the repulsive interactions of $\Delta$ baryons in
hadronic matter.

###### pacs:

##  I Introduction

Neutron stars (NSs) are unique extraterrestrial laboratories to probe matter
under extreme conditions. Recently, the possibility of the formation of
baryons other than nucleons in the cores of NSs has become one of the central
issues in modern nuclear astrophysics Oertel  _et al._ (2017). Many extensive
studies address the appearance of hyperons, especially in the context of the
notorious hyperon puzzle (see, e.g., Glendenning (1985); Drago  _et al._
(2014a); Baym  _et al._ (2018); Lattimer and Prakash (2006); Raduta  _et al._
(2018); Fortin  _et al._ (2017); Oertel  _et al._ (2015); Gomes  _et al._
(2015); Katayama and Saito (2015); Bombaci (2017)). In contrast, much less
research has been conducted on the influence of $\Delta(1232)$ resonance on
the gross properties of the equation of state (EoS). In fact, for many years,
its possible contribution has been rather ignored because early results
suggested that $\Delta(1232)$ resonance appears at densities beyond the
typical central densities of compact objects Glendenning (1985). However,
recent studies based on the relativistic mean-field (RMF) model Drago  _et
al._ (2014b) and microscopic approaches Li and Sedrakian (2019) suggest an
early onset of $\Delta(1232)$ resonance, which, similarly to hyperons, leads
to a significant softening of the EoS. In the last decade, a number of studies
of $\Delta(1232)$ in dense matter were conducted within various approaches
Drago  _et al._ (2014b); Li  _et al._ (2018); Motta  _et al._ (2020); Li  _et
al._ (2020); Li and Sedrakian (2019); Cai  _et al._ (2015); Zhu  _et al._
(2016); Sahoo  _et al._ (2018); Sen (2021); Malfatti  _et al._ (2020); Ribes
_et al._ (2019); Schürhoff  _et al._ (2010); Drago  _et al._ (2014a);
Marczenko  _et al._ (2022). The consequences of the appearance of
$\Delta(1232)$ resonance are found to be crucial for the properties of the EoS
and the structure of NSs. Similar to the hyperon contribution, the early
appearance of $\Delta(1232)$ resonance creates a tension with the maximum NS
mass. Consequently, obtained maximal masses are far below the mass of the PSR
J0740+6620 pulsar Demorest  _et al._ (2010); Antoniadis  _et al._ (2013);
Fonseca  _et al._ (2016); Cromartie  _et al._ (2019); Fonseca  _et al._
(2021). On the other hand, the radius and tidal deformability of a
$1.4~{}M_{\odot}$ NS can be significantly reduced due to the softening of the
EoS linked to the early appearance of $\Delta(1232)$ resonance at low to
intermediate densities Li  _et al._ (2018); Li and Sedrakian (2019).

The recent advancements of multi-messenger astronomy on different dense-matter
astrophysical sources have led to further improvements in constraining the EoS
at low temperature and high density. The modern observatories for measuring
masses and radii of compact objects, the gravitational wave interferometers of
the LIGO/Virgo Collaboration (LVC) Abbott  _et al._ (2018, 2019), and the
X-ray observatory Neutron star Interior Composition Explorer (NICER) provide
new powerful constraints on the neutron-star mass-radius profile Riley  _et
al._ (2019); Miller  _et al._ (2019, 2021); Riley  _et al._ (2021). These
stringent constraints allow for a more systematic study of the influence of
the formation of various degrees of freedom inside the cores of NSs.
Therefore, it is crucial to explore the underlying role of $\Delta(1232)$
resonance in dense matter under extreme conditions.

Central densities of neutron stars lie up to a few times of normal nuclear
density so that it is to be expected that baryons, in particular nucleons and
$\Delta$, change their properties due to the restoration of chiral symmetry.
The recent lattice QCD (LQCD) results Aarts  _et al._ (2015, 2017, 2019)
exhibit a clear manifestation of the parity doubling structure for the low-
lying baryons around the chiral crossover. The observed behavior of parity
partners is likely an imprint of the chiral symmetry restoration in the
baryonic sector of QCD and is expected to occur also in cold dense matter,
including neutron-star conditions. Such properties of the chiral partners can
be described in the framework of the parity doublet model De Tar and Kunihiro
(1989); Jido  _et al._ (2000, 2001). The model has been applied to hot and
dense hadronic matter, neutron stars, as well as the vacuum phenomenology of
QCD Dexheimer  _et al._ (2008); Gallas  _et al._ (2010); Paeng  _et al._
(2012); Sasaki  _et al._ (2011); Gallas  _et al._ (2011); Zschiesche  _et al._
(2007); Benic  _et al._ (2015); Marczenko and Sasaki (2018); Marczenko  _et
al._ (2018, 2019); Marczenko (2020); Marczenko  _et al._ (2020, 2021, 2022);
Mukherjee  _et al._ (2017a, b); Dexheimer  _et al._ (2013); Steinheimer  _et
al._ (2011a); Weyrich  _et al._ (2015); Sasaki and Mishustin (2010); Yamazaki
and Harada (2019a, b); Ishikawa  _et al._ (2019); Steinheimer  _et al._
(2011b); Giacosa (2012); Motohiro  _et al._ (2015); Minamikawa  _et al._
(2021a, b).

In this work, we explore the implications of dynamical restoration of chiral
symmetry within the hadronic phase. To this end, we employ the parity doublet
model for the nucleon and $\Delta$ Jido  _et al._ (2000); Takeda  _et al._
(2018) and explore their implications on the structure of neutron stars. We
demonstrate how modern constraints on the masses, radii, and tidal
deformability of neutron stars can be achieved in this framework.

This paper is organized as follows. In Sec. II, we introduce the parity
doublet model for nucleon and $\Delta(1232)$ resonance. In Sec. III, we
discuss the obtained numerical results on the equation of state and chiral
structure in the isospin-symmetric matter and under neutron-star conditions.
In Sec. IV, we discuss the obtained neutron-star relations and confront the
results with recent observations. Finally, Sec. V is devoted to summary and
conclusions.

##  II Parity doublet model

$m_{N}~{}$[MeV] |  $m_{N}^{\star}~{}$[MeV] |  $m_{\Delta}~{}$[MeV] |  $m_{\Delta}^{\star}~{}$[MeV] |  $m_{\pi}~{}$[MeV] |  $m_{\omega}~{}$[MeV] |  $m_{\rho}~{}$[MeV] |  $f_{\pi}~{}$[MeV]
---|---|---|---|---|---|---|---
939 | 1500 | 1232 | 1700 | 140 | 783 | 775 | 93
Table 1: Physical vacuum inputs and the parity doublet model parameters used in this work. $n_{0}~{}[\rm fm^{-3}]$ | $E/A-m_{N}~{}[\rm MeV]$ |  $K~{}$[MeV] |  $E_{\rm sym}~{}$[MeV]
---|---|---|---
0.16 | -16 | 240 | 31
Table 2: Properties of the nuclear ground state at $\mu_{B}=923~{}$MeV and
$\mu_{Q}=0$, and the symmetry energy used in this work.

In this section, we briefly introduce the parity doublet model for nucleon and
$\Delta(1232)$ resonance capable of describing the chiral symmetry
restoration, following Ref. Takeda  _et al._ (2018). The model is composed of
the baryonic parity doublets for nucleon and $\Delta(1232)$ resonance, and
mesons as in the Walecka model Walecka (1974). The spontaneous chiral symmetry
breaking yields the mass splitting between the two baryonic parity partners,
in each parity doublet with given spin. In this work, we consider a system
with $N_{f}=2$. The baryonic degrees of freedom are coupled to the chiral
fields $\left(\sigma,\boldsymbol{\pi}\right)$, the vector-isoscalar field
($\omega_{\mu}$), and the vector-isovector field ($\boldsymbol{\rho}_{\mu}$).
The thermodynamic potential of the model in the mean-field approximation reads
Takeda  _et al._ (2018)

$\Omega=V_{\sigma}+V_{\omega}+V_{\rho}+\sum_{x=N,\Delta}\Omega_{x}\rm,$ (1)

with the index $N$ labeling collectively positive-parity and negative-parity
spin-$1/2$ nucleons, i.e., $N\in\\{p,n;p^{\star},n^{\star}\\}$, and spin-$3/2$
$\Delta$’s, i.e.,
$\Delta\in\\{\Delta_{++,+,0,-};\Delta^{\star}_{++,+,0,-}\\}$. The negative-
parity states are marked with asterisks. The mean-field potentials in Eq. (1)
read

$\displaystyle V_{\sigma}$
$\displaystyle=-\frac{\lambda_{2}}{2}\sigma^{2}+\frac{\lambda_{4}}{4}\sigma^{4}-\frac{\lambda_{6}}{6}\sigma^{6}-\epsilon\sigma\textrm{,}$
(2a) $\displaystyle V_{\omega}$
$\displaystyle=-\frac{m_{\omega}^{2}}{2}\omega^{2}\textrm{,}$ (2b)
$\displaystyle V_{\rho}$
$\displaystyle=-\frac{m_{\rho}^{2}}{2}\rho^{2}\textrm{,}$ (2c)

where $\lambda_{2}=\lambda_{4}f_{\pi}^{2}-\lambda_{6}f_{\pi}^{4}-m_{\pi}^{2}$,
and $\epsilon=m_{\pi}^{2}f_{\pi}$. $m_{\pi}$, $m_{\omega}$, and $m_{\rho}$ are
the $\pi$, $\omega$, and $\rho$ meson masses, respectively, and $f_{\pi}$ is
the pion decay constant. The parameters $\lambda_{4}$ and $\lambda_{6}$ are
fixed by the properties of the nuclear ground state. We note that the six-
point scalar interaction term in Eq. (2a) is essential in order to reproduce
the experimental value of the compressibility $K=240\pm 20~{}$MeV Shlomo  _et
al._ (2006); Motohiro  _et al._ (2015). In Eqs. (2b) and (2c), $\omega$ and
$\rho$ are the only non-vanishing expectation values of the $\omega_{\mu}$ and
$\boldsymbol{\rho}_{\mu}$ in the mean-field approximation, respectively. The
kinetic part of the thermodynamic potential, $\Omega_{x}$, reads

$\Omega_{x}=\gamma_{x}\int\frac{\mathrm{d}^{3}p}{(2\pi)^{3}}T\left(\ln\left(1-f_{x}\right)+\ln\left(1-\bar{f}_{x}\right)\right)\rm,$
(3)

where the factors $\gamma_{N}=2$ and $\gamma_{\Delta}=4$ denote the spin
degeneracy of both parity partners for nucleons and $\Delta$’s, respectively.
The particle (antiparticle) Fermi-Dirac distribution function reads

$\displaystyle f_{x}$ $\displaystyle=\frac{1}{1+e^{\beta(E_{x}-\mu_{x})}}\rm,$
(4a) $\displaystyle\bar{f}_{x}$
$\displaystyle=\frac{1}{1+e^{\beta(E_{x}+\mu_{x})}}\rm,$ (4b)

with $\beta$ being the inverse temperature, the dispersion relation
$E_{x}=\sqrt{\boldsymbol{p}^{2}+m_{x}^{2}}$ and $\mu_{x}$ is the effective
chemical potential.

The masses of the positive- and negative-parity chiral partners are given by

$m^{x}_{\pm}=\frac{1}{2}\left[\sqrt{\left(g_{1}^{x}+g_{2}^{x}\right)^{2}\sigma^{2}+4\left(m_{0}^{x}\right)^{2}}\mp\left(g^{x}_{1}-g^{x}_{2}\right)\sigma\right]\textrm{,}$
(5)

where $\pm$ sign denotes parity and $x=N,\Delta$ for nucleons and $\Delta$,
respectively. The positive-parity nucleons are identified as the positively
charged and neutral $N(938)$ states: proton ($p$) and neutron ($n$). Their
negative-parity counterparts, denoted as $p^{\star}$ and $n^{\star}$, are
identified as $N(1535)$ resonance Zyla  _et al._ (2020). The positive-parity
$\Delta$ states are identified with $\Delta(1232)$. Their negative-parity
chiral partners, $\Delta^{\star}$, are identified with $\Delta(1700)$ Zyla
_et al._ (2020). For given chirally invariant mass, $m_{0}^{x}$, the
parameters $g_{1}^{x}$ and $g_{2}^{x}$ are determined by the corresponding
vacuum masses (see Table 1).

The effective chemical potentials for nucleons and their chiral partners are
given by

$\displaystyle\mu_{p}$
$\displaystyle=\mu_{p^{\star}}=\mu_{B}+\mu_{Q}-g_{\omega}^{N}\omega-
g_{\rho}^{N}\rho\rm,$ (6a) $\displaystyle\mu_{n}$
$\displaystyle=\mu_{n^{\star}}=\mu_{B}-g_{\omega}^{N}\omega+g_{\rho}^{N}\rho\rm.$
(6b)

The effective chemical potentials for $\Delta$ and their chiral partners are
given by

$\displaystyle\mu_{\Delta_{++}}$
$\displaystyle=\mu_{\Delta^{\star}_{++}}=\mu_{B}+2\mu_{Q}-g_{\omega}^{\Delta}\omega-3g_{\rho}^{\Delta}\rho\rm,$
(7a) $\displaystyle\mu_{\Delta_{+}}$
$\displaystyle=\mu_{\Delta^{\star}_{+}}=\mu_{B}+\mu_{Q}-g_{\omega}^{\Delta}\omega-
g_{\rho}^{\Delta}\rho\rm,$ (7b) $\displaystyle\mu_{\Delta_{0}}$
$\displaystyle=\mu_{\Delta^{\star}_{0}}=\mu_{B}-g_{\omega}^{\Delta}\omega+g_{\rho}^{\Delta}\rho\rm,$
(7c) $\displaystyle\mu_{\Delta_{-}}$
$\displaystyle=\mu_{\Delta^{\star}_{-}}=\mu_{B}-\mu_{Q}-g_{\omega}^{\Delta}\omega+3g_{\rho}^{\Delta}\rho\rm.$
(7d)

The parameters, $g_{\omega}^{x}$ and $g_{\rho}^{x}$ control the coupling
strengths of baryons to $\omega$ and $\rho$ mesons, respectively Takeda  _et
al._ (2018).

In-medium profiles of the mean fields are obtained by extremizing the
thermodynamic potential in Eq. (1), leading to the following gap equations

$\displaystyle\frac{\partial\Omega}{\partial\sigma}$
$\displaystyle=\frac{\partial
V_{\sigma}}{\partial\sigma}+\sum_{x=N,\Delta}s_{x}\frac{\partial
m_{x}}{\partial\sigma}\textrm{,}$ (8a)
$\displaystyle\frac{\partial\Omega}{\partial\omega}$
$\displaystyle=\frac{\partial
V_{\omega}}{\partial\omega}+\sum_{x=N,\Delta}g^{x}_{\omega}n_{x}\textrm{,}$
(8b) $\displaystyle\frac{\partial\Omega}{\partial\rho}$
$\displaystyle=\frac{\partial
V_{\rho}}{\partial\rho}+g^{N}_{\rho}\left(n_{p}-n_{n}+n_{p^{\star}}-n_{n^{\star}}\right)$
(8c)
$\displaystyle+g_{\rho}^{\Delta}(3n_{\Delta_{++}}+n_{\Delta_{+}}-n_{\Delta_{0}}-3n_{\Delta_{-}}$
$\displaystyle\phantom{\;\;\;\;\;\;}+3n_{\Delta^{\star}_{++}}+n_{\Delta^{\star}_{+}}-n_{\Delta^{\star}_{0}}-3n_{\Delta^{\star}_{-}})\rm,$

where the scalar and vector densities are

$s_{x}=\gamma_{x}\int\frac{\mathrm{d}^{3}p}{\left(2\pi\right)^{3}}\frac{m_{x}}{E_{x}}\left(f_{x}+\bar{f}_{x}\right)$
(9)

and

$n_{x}=\gamma_{x}\int\frac{\mathrm{d}^{3}p}{\left(2\pi\right)^{3}}\left(f_{x}-\bar{f}_{x}\right)\textrm{,}$
(10)

respectively.

In the grand canonical ensemble, the thermodynamic pressure is obtained from
the thermodynamic potential as $P=-\Omega+\Omega_{0}$, where $\Omega_{0}$ is
the value of the thermodynamic potential in the vacuum.

The net-baryon number density and net-charge densities for a species $x$ are
defined as

$\displaystyle n^{x}_{B}$
$\displaystyle=-\frac{\partial\Omega_{x}(T,\mu_{B},\mu_{Q})}{\partial\mu_{B}}\textrm{,}$
(11a) $\displaystyle n^{x}_{Q}$
$\displaystyle=-\frac{\partial\Omega_{x}(T,\mu_{B},\mu_{Q})}{\partial\mu_{Q}}\textrm{,}$
(11b)

respectively, where $\Omega_{x}$ is the kinetic term in Eq. (3). The particle-
density fractions are defined as

$Y_{x}=\frac{n_{B}^{x}}{n_{B}}\textrm{.}$ (12)

The compressibility, symmetry energy, and its slope at saturation density are
given as

$\displaystyle K$ $\displaystyle=9n_{0}\frac{\partial\mu_{B}}{\partial
n_{B}}\Bigg{|}_{n_{B}=n_{0}}\textrm{,}$ (13a) $\displaystyle E_{\rm sym}$
$\displaystyle=\frac{1}{2}\frac{\partial^{2}\left(\epsilon/n_{B}\right)}{\partial\delta^{2}}\Bigg{|}_{\delta=0}\textrm{,}$
(13b) $\displaystyle L$ $\displaystyle=3n_{0}\frac{\partial E_{\rm
sym}}{\partial n_{B}}\Bigg{|}_{n_{B}=n_{0}}\textrm{,}$ (13c)

respectively. In Eq. (13b), $\epsilon$ is the energy density and the isospin
asymmetry parameter $\delta=\sum_{i}I_{i}n_{i}/n_{B}$, where $I_{i}$ is the
third component of the isospin operator of the i’th species#1#1#1We note that
this generalized definition reduces to the known definition of the asymmetry
parameter, $\delta=(n_{B}^{p}-n_{B}^{n})/n_{B}$, where the negative parity
states, $N^{\star}$ and $\Delta$ states are not populated, and their densities
vanish, which is the case at the saturation density.. We remark that the
obtained values of $L_{\rm sym}\approx 82~{}$MeV at saturation for
$m_{0}^{N}=550-700~{}$MeV agree with the commonly considered range of the
parameter (Oertel  _et al._ , 2017) and are also found in other parity-doublet
models (Motohiro  _et al._ , 2015). The most updated estimate, $L_{\rm
sym}=53^{+14}_{-15}~{}$MeV, is based on combined astrophysical data, PREX-II,
and recent effective chiral field theory results (Essick  _et al._ , 2021).
Furthermore, the recent analysis within density functional theory yields
$L_{\rm sym}=54\pm 8$ MeV (Reinhard  _et al._ , 2021). We note that the value
of the parameter $L_{\rm sym}$ can be decreased, e.g., by introducing an
$\omega-\rho$ interaction term to the Lagrangian Dexheimer  _et al._ (2019).

The strength of the $g^{N}_{\omega}$ coupling is fixed by the nuclear
saturation properties: the saturation density, $n_{0}$, and the binding
energy, $\epsilon/n_{B}-m_{N}$; the value of $g^{N}_{\rho}$ can be fixed by
fitting the value of symmetry energy. The properties of the nuclear ground
state and the symmetry energy are compiled in Table 2. The couplings of the
$\Delta$ resonance to the meson fields are poorly constrained due to limited
knowledge from experimental observations. The most advocated constraint was
obtained by an analysis of electromagnetic excitations of the $\Delta$ baryon
in the framework of relativistic mean-field model Wehrberger  _et al._ (1989).
It puts a constraint on the relative strength of the scalar and vector
couplings. Other phenomenological studies indicate an attractive $\Delta-N$
potential with no consensus on its actual size Horikawa  _et al._ (1980);
O’Connell and Sealock (1990); Lehr  _et al._ (2000); Nakamura  _et al._
(2010). We note that in the parity doublet model the values of the
$\Delta-\sigma$ couplings, $g_{1}^{\Delta}$ and $g_{2}^{\Delta}$, are uniquely
fixed by requiring the vacuum masses of $\Delta$ and $\Delta^{\star}$. On the
other hand, the nature of the repulsive interaction among $\Delta$ isobars and
their coupling to the $\omega$ and $\rho$ mean fields are still far from
consensus.

Figure 1: The $\sigma$ mean field in the parity doublet model for
$m_{0}^{N}=600~{}$MeV in isospin-symmetric matter at zero temperature as a
function of the baryon chemical potential. Figure 2: Thermodynamic pressure in
isospin-symmetric matter as a function of the net-baryon number density, in
units of the saturation density, $n_{0}=0.16~{}\rm fm^{-3}$.

It is customary to parametrize the $\Delta$-meson couplings in terms of the
nucleon-meson couplings:

$\displaystyle g_{\omega}^{\Delta}$
$\displaystyle=R_{\Delta}g_{\omega}^{N}\textrm{,}$ (14a) $\displaystyle
g_{\rho}^{\Delta}$ $\displaystyle=R_{\Delta}g_{\rho}^{N}\textrm{.}$ (14b)

For simplicity, in the present study, we fix $R_{\Delta}=1$. Detailed
discussion of the dependence of $R_{\Delta}$ is presented in Sec. IV.3. In
general, additional repulsion between $\Delta$’s would systematically shift
their onset in the stellar sequence to higher densities. This eventually would
prevent the neutron stars with $\Delta$ matter from existence in the
gravitationally stable branch of the sequence. We note that this effect is
similar to the case of repulsive interactions between quarks Marczenko  _et
al._ (2020).

Figure 3: Thermodynamic pressure for isospin-symmetric matter as a function of
the net-baryon number density, in units of the saturation density,
$n_{0}=0.16~{}\rm fm^{-3}$ (left panel), and under the NS conditions of
$\beta$-equilibrium and charge neutrality, as a function of the energy
density, $\epsilon$, at $T=0$ (right panel). In the left panel, the orange-
shaded region shows the flow constraint Danielewicz  _et al._ (2002). In the
right panel, the orange- and gray-shaded regions show the constraints obtained
by Annala  _et al._ (2020) and Abbott  _et al._ (2018), respectively.

In the present work, we take four representative values of $m_{0}^{N}=550$,
$600$, $650$, $700~{}$MeV. Because the onset of $\Delta$ matter depends on the
value of the chirally invariant mass $m_{0}^{\Delta}$ Takeda  _et al._ (2018),
we systematically study the influence of $\Delta$ on the EoS and compliance
with terrestrial constraints from heavy-ion collisions Danielewicz  _et al._
(2002), and the astrophysical constraints, i.e., the $2~{}M_{\odot}$ and the
tidal deformability constraints Hebeler  _et al._ (2013); Cromartie  _et al._
(2019); Fonseca  _et al._ (2021); Abbott  _et al._ (2018). In this work, we
put an additional constraint on the chirally invariant mass of $\Delta$.
Namely, we require that $m_{0}^{N}\leq m_{0}^{\Delta}$. Too low values of
$m_{0}^{\Delta}$ lead to the onset of $\Delta$ matter at subsaturation
densities; thus, it spoils the properties of the ground state Takeda  _et al._
(2018). Setting $m_{0}^{\Delta}=\infty$ suppresses the $\Delta$ states and the
EoS effectively corresponds to the purely nucleonic EoS. We note that by the
assumption, $L_{\rm sym}$ does not depend on the choice of $m_{0}^{\Delta}$.

The physical inputs used in this work are summarized in Tabs. 1 and 2.

##  III Chiral Structure and Equation of State

Figure 4: Sequences of masses for neutron stars vs. central net-baryon density
(left panel) and radius (right panel) as solutions of the TOV equations. The
dashed lines correspond to the purely nucleonic EoSs. The solid lines
correspond to the case $m_{0}^{N}=m_{0}^{\Delta}$. The region spanned between
the two lines mark the results obtained for $m_{0}^{N}<m_{0}^{\Delta}$ in each
case. The mass-density and mass-radius curves are plotted up to the maximal
mass solutions for given $m_{0}^{N}$. The orange band in the left panel shows
the PSR J0740+6620 mass constraint $M_{\rm max}=2.08\pm 0.07$ Fonseca  _et
al._ (2021). The inner (outer) orange band in the right panel shows the
$1\sigma$ credibility regions from the NICER analysis of observations of the
massive pulsar PSR J0740+6620 as dark orange Riley  _et al._ (2021) and light
orange Miller  _et al._ (2021) regions. The inner (outer) green and purple
bands show 50% (90%) credibility regions obtained from the recent GW170817
Abbott  _et al._ (2018) event for the low- and high-mass posteriors. Finally,
the inner (outer) gray region corresponds to the mass and radius constraint at
68.2% (95.4%) obtained for PSR J0030+0451 by the group analyzing NICER X-ray
data Miller  _et al._ (2021).

In this section, we discuss the influence of $\Delta$ matter on the EoS in the
parity doublet model in the isospin-symmetric matter, as well as under the NS
conditions of $\beta$-equilibrium and charge neutrality at $T=0$.

In Fig. 1, we show the baryon chemical potential dependence of the expectation
value of the $\sigma$ mean field for the isospin-symmetric matter. To
illustrate the effects of the onset of $\Delta$ matter on the chiral
structure, we fix $m_{0}^{N}=600~{}$MeV and study the dependence on the
parameter $m_{0}^{\Delta}$. In all cases shown in the figure, the results
exhibit similar behavior in the vicinity of the liquid-gas phase transition at
small values of the baryon chemical potential. Thus, the appearance of
$\Delta$ matter does not spoil the properties of the nuclear ground state. In
the case of the purely nucleonic EoS, the chiral transition is a smooth
crossover. The transition happens at $\mu_{B}\simeq 1400~{}$MeV and
corresponds to a sudden drop of $\sigma$ to about zero, causing the parity-
doublet nucleons ($N$, $N^{\star}$) to become nearly equally populated (cf.
Eq. (5)). In general, the appearance of $\Delta$ matter is seen as a change in
the stiffness of the EoS above the saturation density; thus, deviations from
the purely nucleonic EoS can be attributed to the onset of $\Delta$ matter.
For the case $m_{0}^{\Delta}=650~{}$MeV, $\sigma$ deviates smoothly from the
purely nucleonic result around $\mu_{B}=1100~{}$MeV, which corresponds to the
onset of positive-parity $\Delta$ matter. After this point, $\sigma$ changes
very slowly, and the population of chiral partners is seen only mildly around
$\mu_{B}=1400~{}$MeV and $\mu_{B}=1900~{}$MeV, for $N^{\star}$ and
$\Delta^{\star}$, respectively. For the case $m_{0}^{\Delta}=600~{}$MeV, the
swift decrease of $\sigma$ at $\mu_{B}=950$ MeV corresponds to the onset of
positive-parity $\Delta$ through a first-order phase transition. The chiral
symmetry gets restored around $\mu_{B}=1500~{}$MeV where the chiral partners
of the nucleon and $\Delta$ appear almost simultaneously. We note that for
smaller values of $m_{0}^{\Delta}$, $\Delta$ enters the matter at smaller
baryon chemical potentials before $N^{\star}$ is populated in the purely
nucleonic EoS. In both cases, it is associated with a drop of $\sigma$. At the
same time, the onset of $\Delta^{\star}$, and thus the full restoration of
chiral symmetry, is shifted to higher baryon chemical potentials.

In Fig. 2, we plot the thermodynamic pressure for the same set of EoSs. In
general, the chiral structure resembles the stiffness of the EoS. The
softening of the pressure in the purely nucleonic EoS around $4~{}n_{0}$
corresponds to the chiral symmetry restoration and the onset of $N^{\star}$.
Similar to the $\sigma$ expectation value, the onset of $\Delta$ matter is
signaled by deviations from the template nucleonic EoS. For
$m_{0}^{\Delta}=650~{}$MeV the onset is smooth and only softens the
thermodynamic pressure. However, for $m_{0}^{\Delta}=600~{}$MeV, the onset of
$\Delta$ through the first-order transition is pronounced in the density jump
of the order of $0.5~{}n_{0}$. In general, smaller values of $m_{0}^{\Delta}$
result in the earlier onset of $\Delta$ matter. However, the strength of the
transition depends on the choice of $m_{0}^{\Delta}$. Namely, higher values
yield weaker first-order transition, which turns into a second-order and
eventually becomes a smooth crossover, defined as a peak in
$\partial\sigma/\partial\mu_{B}$. When $m_{0}^{\Delta}\rightarrow\infty$, the
onset of $\Delta$ matter is shifted to higher densities and the EoS converges
to the purely nucleonic case. We note that the structure discussed for
$m_{0}^{N}=600~{}$MeV remains qualitatively the same for other values of the
parameter.

Figure 5: Onset densities of particle species for $m_{0}^{N}=600~{}$MeV as a
function of $m_{0}^{\Delta}$. The gray-shaded regions indicate densities
beyond the maximum mass configurations. The orange, solid line marks the
central density of $1.4~{}M_{\odot}$ neutron stars. We note that the onset
densities of the positive-parity nucleons, proton and neutron, are not shown
in the figure.

In the left panel of Fig. 3, shown are the EoSs for selected values of
$m_{0}^{N}$. To illustrate the effect of $\Delta$ matter on the EoS at
intermediate densities, we show results obtained for purely nucleonic EoS
(dashed line) together with the case $m_{0}^{\Delta}=m_{0}^{N}$ (solid line).
The regions bounded by the two results correspond to the range spanned by
solutions with $m_{0}^{N}<m_{0}^{\Delta}$ in each case. In the left panel of
the figure, we show the isospin-symmetric matter EoSs. The orange-shaded
region shows the proton flow constraint Danielewicz  _et al._ (2002). In
general, the low-density behavior in each case is similar, until the
deviations from the purely nucleonic EoSs are induced by the onset of $\Delta$
matter. The softening due to the onset of $\Delta$ resonance yields better
agreement with the proton flow constraint. The EoSs with $\Delta$ converge
back to the purely nucleonic EoS at high densities.

In the right panel of Fig. 3, we show the corresponding EoSs under the
neutron-star conditions. The grey- and orange-shaded envelopes show the
constraints derived in Abbott  _et al._ (2018) and Annala  _et al._ (2020),
respectively. Similar to the isospin-symmetric case, the low-density behavior
is comparable. For $m_{0}^{N}=550~{}$MeV, the EoSs with $m_{0}^{\Delta}\approx
m_{0}^{N}$ result in an appearance of $\Delta_{-}$ via a strong first-order
phase transition and underestimate the constraint at low densities. However,
the phase transition is followed by a subsequent stiffening as compared to the
purely nucleonic case and the EoS reaches the constraint at higher densities.
It suggests that the onset of $\Delta$ matter softens the EoS and stiffens it
at larger densities. This also resembles a higher speed of sound at
intermediate densities. We note that this effect is more readily pronounced
for smaller values of $m_{0}^{\Delta}$. For other parametrizations shown in
the figure, the EoSs fall into the region derived by the constraint. We note
that in the study of cold and dense QCD, commonly used are separate effective
models for the hadronic and quark matter phases (two-phase approaches) with a
priori assumed first-order phase transition, typically associated with
simultaneous chiral and deconfinement transitions Alvarez-Castillo  _et al._
(2016). As recently demonstrated in Marczenko  _et al._ (2022), such a strong
phase transition with large latent heat can occur also within hadronic matter
due to the onset of $\Delta$ matter being subject to chiral symmetry
restoration.

Figure 6: Particle fractions as a function of the net-baryon density in the
units of saturation density for $m_{0}^{N}=600~{}$MeV and different values of
$m_{0}^{\Delta}$. The top panel shows the results for the parity double model
without $\Delta$ matter. The gray-shaded regions indicate densities beyond the
maximum mass configurations. The solid (dashed) lines show the positive-parity
(negative-parity) particles. The color coding is the same as in Fig. 5. The
electrons and muons are shown as gray, dotted and dashed, lines, respectively.

##  IV Properties of neutron stars

###  IV.1 TOV solutions

Figure 7: The relation between tidal deformabilities $\Lambda_{1}$ and
$\Lambda_{2}$ of two compact stars that merged in the GW170817 event Abbott
_et al._ (2018). For comparison, also shown are the 50% and 90% fidelity
regions from the analysis of the GW signal by the LIGO-VIRGO Collaboration
Abbott  _et al._ (2018). Shown is only the physical region where
$\Lambda_{1}<\Lambda_{2}$.

In this section, we explore the impact of the emergence of the $\Delta$ matter
at supersaturation densities on the structure of neutron stars. The
composition of neutron-star matter requires $\beta$-equilibrium with leptons
($l$): electrons ($e$) and muons ($\mu$), included as free relativistic
particles, as well as the charge neutrality condition.

The EoS of dense matter plays a crucial role in determination of the structure
of neutron stars. Its microscopic details are linked to the specific values of
mass and radius, and therefore compactness. We use the EoSs introduced in the
previous section to solve the general-relativistic Tolman-Oppenheimer-Volkoff
(TOV) equations Tolman (1939); Oppenheimer and Volkoff (1939) for spherically
symmetric objects,

$\displaystyle\frac{\mathrm{d}P(r)}{\mathrm{d}r}$
$\displaystyle=-\frac{\left(\epsilon(r)+P(r)\right)\left(M(r)+4\pi
r^{3}P(r)\right)}{r\left(r-2M(r)\right)}\textrm{,}$ (15a)
$\displaystyle\frac{\mathrm{d}M(r)}{\mathrm{d}r}$ $\displaystyle=4\pi
r^{2}\epsilon(r)\textrm{,}$ (15b)

with the boundary conditions $P(r=R)=0$, $M=M(r=R)$, where $R$ and $M$ are the
radius and the mass of a neutron star, respectively. Once the initial
conditions are specified based on a given equation of state, namely the
central pressure $P_{c}$ and the central energy density $\epsilon_{c}$, the
internal profile of a neutron star can be calculated.

In general, there is one-to-one correspondence between the EoS and the mass-
radius relation calculated from Eqs. (15). In the left panel of Fig. 4, we
show the relationship of mass versus central net-baryon density, for the
calculated sequences of neutron stars for $m_{0}^{N}=550~{}$MeV,
$m_{0}^{N}=600~{}$MeV, $m_{0}^{N}=650~{}$MeV, $m_{0}^{N}=700~{}$MeV, together
with the state-of-the-art constraint on the maximum mass for the pulsar PSR
J0740+6620, $M=2.08\pm 0.07~{}M_{\odot}$ Fonseca  _et al._ (2021). Shown are
sequences for $m_{0}^{\Delta}=m_{0}^{N}$ (solid line) and pure nucleonic
sequence (dashed line). The color-shaded areas show the region spanned by
solutions obtained with $m_{0}^{\Delta}>m_{0}^{N}$. The relations are plotted
up to the maximally stable solutions. We point out that, similarly to the EoS,
the onset of positive-parity $\Delta$ leads to a softening of the mass-radius
sequence so that it is accompanied by a strong increase of the central
density. This is seen for $m_{0}^{N}=m_{0}^{\Delta}$ with plateaux at small NS
masses. We note that in general smaller values of $m_{0}^{N}$ yield larger
values of the maximum mass, owing to the stiffness of the EoS.

In the right panel of Fig. 4, we show the relationship of mass and radius.
Also shown are the state-of-the-art constraints: the high precision mass-
radius analysis of the massive pulsar PSR J0740+6620 by the NICER
collaboration Riley  _et al._ (2021); Miller  _et al._ (2021), constraint from
the recent GW170817 event Abbott  _et al._ (2018), and the constraint obtained
for PSR J0030+0451 by the group analyzing X-ray data Miller  _et al._ (2019).
The appearance of $\Delta$ matter affects the NS structure which is reflected
in the mass-radius relations. As discussed in Sec. III, the softening of the
EoS due to the appearance of $\Delta$ is followed by a subsequent stiffening
of the EoS. A substantial reduction of the radii is observed. For instance,
the radii of the $1.4~{}M_{\odot}$ NS obtained in the purely nucleonic EoS
reduce roughly up to $1~{}$km for $m_{0}^{N}=700~{}$MeV, to almost $2~{}$km
for $m_{0}^{N}=550~{}$MeV. On the other hand, the decrease of the star’s
maximum mass is seen only mildly. This, in turn, is a consequence of the
subsequent stiffening of the EoS at higher densities, which allows for massive
stars to sustain from the gravitational collapse. In Ref. Li  _et al._ (2018),
it was shown that the earlier onset of $\Delta$ yields a larger maximum mass.
We note that this stays in contrast to our work, where we find a mild decrease
of the maximum mass. This is due to additional softening of the EoS at high
densities that is provided by the onset of the negative-parity $N^{\star}$ and
$\Delta^{\star}$ due to the chiral symmetry restoration. This highlights the
importance of chiral symmetry restoration in the EoS and the properties of the
NS matter.

In Fig. 5 we study the onset densities of the particle species. As an example,
we show the results for a fixed value of $m_{0}^{N}=600~{}$MeV. We note that
the onset densities of the positive-parity nucleons, proton, and neutron, are
not shown in the figure. Their negative-parity counterparts, $n^{\star}$ and
$p^{\star}$, are always present for any choice of $m_{0}^{\Delta}$. They enter
the matter composition roughly between $3-5n_{0}$. Their appearance is always
sequential due to differences in their effective chemical potentials.
Conversely, for small values of $m_{0}^{\Delta}$, the positive-parity $\Delta$
resonances appear almost simultaneously through a first-order transition,
which triggers a sufficiently big density jump. At higher values of the
parameter, the differences in the onset density become more readily exposed.
The onset density of $\Delta_{-}$ changes only mildly within the shown range
of the parameter, while the densities for others increase and eventually the
$\Delta_{++}$, $\Delta_{+}$, and $\Delta_{0}$ are being populated in the
gravitationally unstable branch of the EoS. We find similar behavior for the
negative-parity $\Delta^{\star}$’s. They appear sequentially and no
$\Delta^{\star}$’s are present for $m_{0}^{\Delta}>670~{}$MeV. We stress that
similarly to matter constituents in symmetric matter, the appearance of the
chiral partners, reflects the partial chiral symmetry restoration. For
$m_{0}^{N}=600~{}$MeV the $\Lambda_{1.4}<580$ constraint is met for
$m_{0}^{\Delta}<650~{}$MeV (see Ref. Marczenko  _et al._ (2022)). Therefore,
the appearance of $\Delta$ matter is essential for softening the EoS at
intermediate densities. Although $\Delta^{\star}$’s are not populated at
central densities of $1.4~{}M_{\odot}$ NSs, the compliance with the
deformability constraint implies their onset at higher, gravitationally
stable, densities, and therefore, their presence in the cores of high-mass
neutron stars. This is a direct consequence of the chiral symmetry
restoration. Namely, the onset of $\Delta$ matter drastically decreases the
value of the $\sigma$ mean-field, and consequently, decreases the mass of
$\Delta^{\star}$’s towards the asymptotic value of $m_{0}^{\Delta}$, allowing
for their earlier onset. For sufficiently high $m_{0}^{\Delta}$, $\Delta$
matter is not populated in the gravitationally stable part of the mass-radius
sequence. The shift of the threshold of the onset density is qualitatively
similar as in the case of increasing the strength of repulsion between
$\Delta$’s, which may push their onset out of the gravitationally stable
branch of the stellar sequence. Thus, the corresponding EoSs would be
equivalent to the purely nucleonic EoS in this range. Lastly, we note that the
maximum-mass central density and the central density of a $1.4~{}M_{\odot}$ NS
decrease only mildly by approximately $0.5~{}n_{0}$ as the value of
$m_{0}^{\Delta}$ increases.

Figure 8: $\Delta$ matter fraction under the NS conditions of $\beta$
equilibrium and charge neutrality as a function of the net-baryon density in
the units of saturation density.

In Fig. 6, we show the neutron-star matter composition for
$m_{0}^{N}=600~{}$MeV, in the case without $\Delta$ (first panel),
$m_{0}^{\Delta}=750~{}$MeV (second panel), $m_{0}^{\Delta}=650~{}$MeV (third
panel), $m_{0}^{\Delta}=600~{}$MeV (fourth panel). The gray-shaded regions
mark the central densities that are gravitationally unstable. Because of
charge neutrality, the inclusion of $\Delta$ and $\Delta^{\star}$ influences
the composition of neutron stars. We observe that the onset of $\Delta^{-}$,
and consequently other $\Delta$ states, is shifted to higher values of
$m_{0}^{\Delta}$. Negatively charged $\Delta_{-}$ and $\Delta_{-}^{\star}$
partially replace the role of muons in compensating the charge of positively
charged baryons. For instance, the proton fraction drastically increases upon
the onset of $\Delta_{-}$. This may have consequences for the thermal
evolution of NSs and different cooling mechanisms, which are sensitive to the
proton fraction Klahn  _et al._ (2006). The direct Urca processes that include
the $\Delta$ resonance may be essential for rapid cooling of NSs as they can
become operative without the presence of exotic states nor the large proton
concentration ($11-14\%$), which is required for the canonical direct Urca
process Prakash  _et al._ (1992). Moreover, the threshold proton concentration
for the direct Urca process reduces by about 3% when the chiral symmetry is
restored, i.e., nucleons and their chiral partners are equally populated
Marczenko  _et al._ (2018).

Figure 9: Mass-radius sequences for different values of the parameter
$R_{\Delta}$ for $m_{0}^{N}=650~{}$MeV. Black, solid line shows the mass-
radius sequences obtained for purely nucleonic EoS.

###  IV.2 Tidal deformability

Figure 10: Various NS properties as functions of the parameter $R_{\Delta}$
for $m_{0}^{N}=650~{}$MeV and different values of $m_{0}^{\Delta}$. Gray-
shaded areas show regions excluded by the astrophysical constraints (see text
for details).

The dimensionless tidal deformability parameter $\Lambda$ can be computed
through its relation to the Love number $k_{2}$ Hinderer (2008); Damour and
Nagar (2009); Binnington and Poisson (2009); Yagi and Yunes (2013); Hinderer
_et al._ (2010),

$\Lambda=\frac{2}{3}k_{2}C^{-5}\textrm{,}$ (16)

where $C=M/R$ is the star compactness parameter, with $M$ and $R$ being the
total mass and radius of a star. The Love number $k_{2}$ reads

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

where $y=R\beta(R)/H(R)$. The functions $H(r)$ and $\beta(r)$ are given by the
following set of differential equations:

$\displaystyle\frac{\mathrm{d}\beta}{\mathrm{d}r}$ $\displaystyle=$
$\displaystyle 2\left(1-2\frac{M(r)}{r}\right)^{-1}$ (18) $\displaystyle
H\left\\{-2\pi\left[5\epsilon(r)+9P(r)+\frac{\mathrm{d}\epsilon}{\mathrm{d}P}(\epsilon(r)+P(r))\right]\phantom{\frac{3}{r^{2}}}\right.$
$\displaystyle\left.+\frac{3}{r^{2}}+2\left(1-2\frac{M(r)}{r}\right)^{-1}\left(\frac{M(r)}{r^{2}}+4\pi
rP(r)\right)^{2}\right\\}$
$\displaystyle+\frac{2\beta}{r}\left(1-2\frac{M(r)}{r}\right)^{-1}$
$\displaystyle\left\\{\frac{M(r)}{r}+2\pi
r^{2}(\epsilon(r)-P(r))-1\right\\}~{}\textrm{,}$
$\displaystyle\frac{\mathrm{d}H}{\mathrm{d}r}$ $\displaystyle=$
$\displaystyle\beta\textrm{.}$ (19)

The above equations have to be solved along with the TOV equations (15). The
initial conditions are $H(r\rightarrow 0)=c_{0}r^{2}$ and $\beta(r\rightarrow
0)=2c_{0}r$, where $c_{0}$ is a constant, which is irrelevant in the
expression for the Love number $k_{2}$.

In Fig. 7, we plot the tidal deformability parameters $\Lambda_{1}$ vs
$\Lambda_{2}$ of the high- and low-mass members of the binary merger together
with the 50% and 90% fidelity regions obtained by the LVC analysis of the
GW170817 event Abbott  _et al._ (2018). We note that the tidal deformability
parameter requires sufficiently soft EoS up to a few times the saturation
density. This is seen in the figure, where the smallest tidal deformability
and thus the best agreement with the constraint is obtained for smallest
values of $m_{0}^{\Delta}$, which corresponds to softest EoSs at low density.
On the other hand, the $2~{}M_{\odot}$ requires a sufficiently stiff equation
of state at higher densities. Inversely to the tidal deformability, the most
massive stars are obtained for the stiffest EoSs, namely for purely nucleonic
ones. Therefore, the two constraints are exclusive, which allows for their
precise determination. In Ref. Marczenko  _et al._ (2022), it was shown that
the chirally invariant masses can be roughly estimated to lie in the range
from 550 to 680 MeV. We remark that these values are at tension with the
results obtained in the LQCD simulations Aarts  _et al._ (2017) at vanishing
chemical potential and finite temperature where within the errors the masses
of the ground state nucleon and $\Delta$ do not deviate from their vacuum
masses. This might suggest additional medium dependence of the chirally
invariant masses, whereas the LQCD simulations were performed with large pion
mass, $m_{\pi}=400$ MeV. We also note that the interplay between the
$2~{}M_{\odot}$ and the tidal deformability constraints can be further used to
fix the allowed range of external model parameters. This would be of
particular use in a class of effective models in which the low- and high-
density regimes are not treated independently but rather combined in a
consistent unified framework.

$R_{\Delta}$ |  $m_{0}^{\Delta}~{}$[MeV] |  $R_{1.4}~{}$[km] | $n^{c}_{1.4}~{}[n_{0}]$ | $\Lambda_{1.4}$ | $Y^{1.4}_{\Delta+\Delta^{\star}}$ |  $R_{\rm max}~{}$[km] | $n^{c}_{\rm max}~{}[n_{0}]$ | $M_{\rm max}~{}[M_{\odot}]$ | $Y^{\rm max}_{\Delta+\Delta^{\star}}$
---|---|---|---|---|---|---|---|---|---
0.95 | 650 | 11.2 | 3.34 | 245 | 0.84 | 10.3 | 7.25 | 1.97 | 0.99
670 | 11.7 | 3.23 | 290 | 0.70 | 10.5 | 7.13 | 1.98 | 0.92
700 | 12.3 | 2.92 | 393 | 0.40 | 10.8 | 6.81 | 2.00 | 0.75
1.00 | 650 | 11.9 | 3.04 | 322 | 0.67 | 10.7 | 6.78 | 2.04 | 0.77
670 | 12.3 | 2.88 | 397 | 0.47 | 10.9 | 6.64 | 2.05 | 0.66
700 | 12.7 | 2.67 | 486 | 0.25 | 11.2 | 6.36 | 2.08 | 0.49
1.10 | 650 | 12.9 | 2.56 | 531 | 0.22 | 11.4 | 6.06 | 2.13 | 0.21
670 | 13.1 | 2.50 | 577 | 0.16 | 11.5 | 5.97 | 2.13 | 0.13
700 | 13.2 | 2.43 | 639 | 0.11 | 11.6 | 5.89 | 2.13 | 0.07
no $\Delta$ | 13.4 | 2.30 | 701 | 0.00 | 11.8 | 5.76 | 2.13 | 0.00
Table 3: Dependence of the NS properties on the parameter $R_{\Delta}$ for
$m_{0}^{N}=650~{}$MeV and different values of $m_{0}^{\Delta}$. The bottom row
contains values obtained in purely nucleonic EoS.

###  IV.3 $R_{\Delta}$ dependence

In this section, we explore the influence of the $\Delta$ couplings to vector
mesons on the properties of dense matter. We start by considering the isospin-
symmetric case, i.e., $\mu_{Q}=0$, which also implies that $\rho=0$. Firstly,
we analyze the asymptotic matter composition. Due to chiral symmetry
restoration at high density, baryon masses converge to constant values
$m_{0}^{x}$ (cf. Eq (5)). Because $m_{0}^{x}\ll\mu_{B}$, their contribution
can be neglected at large values of chemical potential, thus

$n_{B}^{x}\approx\frac{d_{x}}{6\pi^{2}}\mu_{x}^{3}\textrm{,}$ (20)

where $d_{x}$ and $\mu_{x}$ are the degeneracy factor and effective chemical
potential of the $x$’th species, respectively. The effective chemical
potential of the nucleon is $\mu_{N}=\mu_{B}-g_{\omega}^{N}\omega$. The
effective chemical potential for $\Delta$ can be written as follows

$\mu_{\Delta}=\mu_{B}-g_{\omega}^{\Delta}\omega=\mu_{N}-\left(R_{\Delta}-1\right)g_{\omega}^{N}\omega\textrm{,}$
(21)

thus, their difference

$\mu_{N}-\mu_{\Delta}=\left(R_{\Delta}-1\right)g_{\omega}^{N}\omega\textrm{,}$
(22)

is directly proportional to the value of $\omega$ mean field. In general, the
value of $\omega$ is expected to increase with density (cf. Eq. (8a)). Thus,
the difference given in Eq. (22) increases with density as well. For
$R_{\Delta}=1$, the chemical potentials are equal, i.e.,
$\mu_{N}=\mu_{\Delta}$. The densities of nucleons and $\Delta$ approach the
same value and the matter composition can be determined solely based on their
degeneracy factors, i.e., $Y_{x}=d_{x}/\sum_{i}d_{i}$. In this case, the
fraction of $\Delta$ matter amounts to $Y_{\Delta+\Delta^{\star}}=0.8$. For
$R_{\Delta}<1$, $\mu_{\Delta}>\mu_{N}$ and their difference grows with
density. The partial density $Y_{\Delta+\Delta^{\star}}$ will asymptotically
approach unity. Conversely, for $R_{\Delta}>1$, $\mu_{\Delta}<\mu_{N}$, thus
the system will be dominated by nucleons, and $Y_{\Delta+\Delta^{\star}}$ will
vanish at high density.

We note that the above analysis can be also performed for matter under NS
conditions of $\beta$ equilibrium and charge neutrality. This case is depicted
in Fig. 8, where shown are the results obtained for three values of the
parameter $R_{\Delta}=0.95,~{}1,~{}1.10$ (see Eq. (14)) for $m_{0}^{N}=650$
MeV and $m_{0}^{\Delta}=650,~{}670,~{}700~{}$MeV. The results converge to the
appropriate asymptotic values for different values of $R_{\Delta}$.

Densities in the cores of neutron stars can reach the values of the order of
$5-8~{}n_{0}$ (cf. Fig. 4). Already under such conditions, qualitative
differences in the abundance of $\Delta$ matter are vividly seen in Fig. 8,
even for small deviations of $R_{\Delta}$ from unity. For instance, in the
considered examples, for a canonical $1.4~{}M_{\odot}$ NS $\Delta$ matter
concentration can be up to four times smaller for $R_{\Delta}=1.10$ when
compared to $R_{\Delta}=0.95$. Interestingly, we find that for $R_{\Delta}>1$
the fraction of $\Delta$ matter exhibits a maximum at few times $n_{0}$, owing
to the asymptotic behavior, which requires that $\Delta$ matter should
disappear from the system at high densities. Similarly for $R_{\Delta}<1$, the
matter at high densities is composed solely of $\Delta$s, and the nucleons are
not present. We note that different studies suggested various ranges of the
repulsive forces (see, e.g, Maslov  _et al._ (2015); Drago  _et al._ (2014b)).
The differences in the matter composition due to the repulsive interactions
are also seen in the bulk NS properties. In Fig. 9, shown are the
corresponding mass-radius sequences. As expected, more compact solutions are
obtained for smaller values of $m_{0}^{\Delta}$ and $R_{\Delta}$. Contrary to
the previous finding that the decrease of $m_{0}^{\Delta}$ changes the maximum
mass only mildly, smaller values of $R_{\Delta}$ yield noticeably smaller
maximum masses. In Fig. 10, we plot various NS properties as functions of
$R_{\Delta}$. Besides the decrease in radius of the canonical
$1.4~{}M_{\odot}$ NS and the maximum mass for smaller values of $R_{\Delta}$,
we also observe a notable decrease in the tidal deformability parameter
$\Lambda_{1.4}$. Therefore, as far as the repulsive interactions are
concerned, the matter composition inside the cores of $2~{}M_{\odot}$ NSs may
be very different upon small variations in the strength of the repulsive
interactions among $\Delta$ isobars. Thus, determining the NS matter
composition has to be treated with particular care. For completeness, the NS
properties are listed in Table 3.

##  V Summary

We have explored the influence of the formation of $\Delta$ resonance in dense
matter on the bulk properties of neutron stars (NSs). To this end, we employed
the parity doublet model for nucleons and $\Delta$’s, which accounts for the
self-consistent treatment of the chiral symmetry breaking and its restoration
in the mesonic and baryonic sectors. We analyzed the equation of state (EoS)
under the NS conditions of $\beta$-equilibrium and charge neutrality.

We have shown how modern astrophysical constraints on the maximum mass Fonseca
_et al._ (2021), the tidal deformability from the binary merger GW170817
Abbott  _et al._ (2018), and recent simultaneous mass-radius constraints from
the NICER experiment Miller  _et al._ (2021, 2019); Riley  _et al._ (2019,
2021) allow us for a consistent determination of the EoS of dense matter. We
have demonstrated that the purely nucleonic EoSs obtained in the parity
doublet model are too stiff at intermediate densities and therefore they are
ruled out by the recently revised tidal deformability constraint. This is
caused by the apparent correlation between the properties of matter at
saturation and the compactness at a few times the saturation density.

We find that the early appearance of $\Delta(1232)$ resonance softens the EoS
at intermediate densities. As a result, the radius of canonical
$1.4~{}M_{\odot}$ reduces substantially providing better agreement with the
tidal deformability constraint from the GW170817 event. We observe a
subsequent stiffening of the EoS at higher densities. As a consequence, we
also find, in contrast to previous studies Li  _et al._ (2018), that the
maximum mass decreases compared to purely nucleonic results, and stays in
agreement with observational astrophysical data. This is traced back to the
restoration of chiral symmetry, which further softens the EoS at higher
densities and results in the onset of $N^{\star}$ and $\Delta^{\star}$
resonances in the high-mass part of the NS sequence.

We emphasize that an abrupt change in mass-radius profile in the high-mass
part of the sequence can be in general due to phase transition in nuclear
matter. However, as discussed above and also shown in Marczenko  _et al._
(2022); Somasundaram  _et al._ (2021), such a change is not necessarily linked
to a hadron-quark phase transition that implies the existence of quark matter
in the cores of NSs, as suggested in Alvarez-Castillo  _et al._ (2016); Li
_et al._ (2021); Christian and Schaffner-Bielich (2021).

To address in more detail the role of baryonic resonances in the presence of
chiral symmetry restoration, we have analyzed the repulsive $\Delta$-(vector
meson) interactions and studied their consequences for the phenomenological
description of NSs. We find that small repulsion yields more compact stellar
sequences with reduced maximum mass, and that variations of the strength of
these interactions lead to notable qualitative changes in the matter
composition at a few times the saturation density. Thus, the determination of
the structure of the NS matter requires detailed knowledge of the intricate
repulsive forces between individual hadronic constituents.

In the present study, we have assumed a hadronic EoS to conclude the
composition of neutron star interiors and the phase structure of dense matter.
A natural extension is to generalize the parity doublet model with $\Delta$
and include the quark degrees of freedom that may allow the existence of
hybrid quark-hadron stars through a deconfinement phase transition while
fulfilling the maximum mass constraint Marczenko  _et al._ (2020).
Furthermore, it is interesting to verify if signatures of the chiral symmetry
restoration related to a crossover or first-order phase transition in the
hadronic phase have an observable imprint on the GW emission from NS mergers,
and, thus, can be measurable in future GW detections.

Given the recent formulation of the three-flavor parity doubling Steinheimer
_et al._ (2011a); Sasaki (2018) and further lattice QCD studies Aarts  _et
al._ (2015, 2017), which indicate that the parity doubling occurs also in the
hyperon channels, it would be of great interest to extend this work and
include contributions of hyperons.

Furthermore, one anticipates that the $\Delta$ resonance may be essential for
rapid cooling of NSs Prakash  _et al._ (1992). Thus, it is alluring to examine
how the presence of the chiral partners of the nucleon and $\Delta$ affect the
NS thermal evolution. Work in these directions is in progress.

###### Acknowledgements.

This work is supported partly by the Polish National Science Centre (NCN)
under OPUS Grant No. 2018/31/B/ST2/01663 (K.R. and C.S.), Preludium Grant No.
2017/27/N/ST2/01973 (M.M.), and the program Excellence Initiative–Research
University of the University of Wroclaw of the Ministry of Education and
Science (M.M). K.R. also acknowledges the support of the Polish Ministry of
Science and Higher Education. M.M. acknowledges helpful discussions with David
E. Alvarez-Castillo. We also acknowledge fruitful discussions with David
Blaschke and Tobias Fischer.

## References

  * Oertel  _et al._ (2017) M. Oertel, M. Hempel, T. Klähn, and S. Typel, Rev. Mod. Phys. 89, 015007 (2017), arXiv:1610.03361 [astro-ph.HE] .
  * Glendenning (1985) N. K. Glendenning, Astrophys. J. 293, 470 (1985).
  * Drago  _et al._ (2014a) A. Drago, A. Lavagno, and G. Pagliara, Phys. Rev. D 89, 043014 (2014a), arXiv:1309.7263 [nucl-th] .
  * Baym  _et al._ (2018) G. Baym, T. Hatsuda, T. Kojo, P. D. Powell, Y. Song, and T. Takatsuka, Rept. Prog. Phys. 81, 056902 (2018), arXiv:1707.04966 [astro-ph.HE] .
  * Lattimer and Prakash (2006) J. M. Lattimer and M. Prakash, Nucl. Phys. A 777, 479 (2006).
  * Raduta  _et al._ (2018) A. R. Raduta, A. Sedrakian, and F. Weber, Mon. Not. Roy. Astron. Soc. 475, 4347 (2018), arXiv:1712.00584 [astro-ph.HE] .
  * Fortin  _et al._ (2017) M. Fortin, S. S. Avancini, C. Providência, and I. Vidaña, Phys. Rev. C 95, 065803 (2017), arXiv:1701.06373 [nucl-th] .
  * Oertel  _et al._ (2015) M. Oertel, C. Providência, F. Gulminelli, and A. R. Raduta, J. Phys. G 42, 075202 (2015), arXiv:1412.4545 [nucl-th] .
  * Gomes  _et al._ (2015) R. O. Gomes, V. Dexheimer, S. Schramm, and C. A. Z. Vasconcellos, Astrophys. J. 808, 8 (2015), arXiv:1411.4875 [astro-ph.SR] .
  * Katayama and Saito (2015) T. Katayama and K. Saito, Phys. Lett. B 747, 43 (2015), arXiv:1501.05419 [nucl-th] .
  * Bombaci (2017) I. Bombaci, JPS Conf. Proc. 17, 101002 (2017), arXiv:1601.05339 [nucl-th] .
  * Drago  _et al._ (2014b) A. Drago, A. Lavagno, G. Pagliara, and D. Pigato, Phys. Rev. C 90, 065809 (2014b), arXiv:1407.2843 [astro-ph.SR] .
  * Li and Sedrakian (2019) J. J. Li and A. Sedrakian, Astrophys. J. Lett. 874, L22 (2019), arXiv:1904.02006 [nucl-th] .
  * Li  _et al._ (2018) J. J. Li, A. Sedrakian, and F. Weber, Phys. Lett. B 783, 234 (2018), arXiv:1803.03661 [nucl-th] .
  * Motta  _et al._ (2020) T. F. Motta, A. W. Thomas, and P. A. M. Guichon, Phys. Lett. B 802, 135266 (2020), arXiv:1906.05459 [nucl-th] .
  * Li  _et al._ (2020) J. J. Li, A. Sedrakian, and F. Weber, Phys. Lett. B 810, 135812 (2020), arXiv:2010.02901 [astro-ph.HE] .
  * Cai  _et al._ (2015) B.-J. Cai, F. J. Fattoyev, B.-A. Li, and W. G. Newton, Phys. Rev. C 92, 015802 (2015), arXiv:1501.01680 [nucl-th] .
  * Zhu  _et al._ (2016) Z.-Y. Zhu, A. Li, J.-N. Hu, and H. Sagawa, Phys. Rev. C 94, 045803 (2016), [Addendum: Phys.Rev.C 102, 019905 (2020)], arXiv:1607.04007 [nucl-th] .
  * Sahoo  _et al._ (2018) H. S. Sahoo, G. Mitra, R. Mishra, P. K. Panda, and B.-A. Li, Phys. Rev. C 98, 045801 (2018), arXiv:1807.06453 [nucl-th] .
  * Sen (2021) D. Sen, Phys. Rev. C 103, 045804 (2021), arXiv:2103.14136 [nucl-th] .
  * Malfatti  _et al._ (2020) G. Malfatti, M. G. Orsaria, I. F. Ranea-Sandoval, G. A. Contrera, and F. Weber, Phys. Rev. D 102, 063008 (2020), arXiv:2008.06459 [astro-ph.HE] .
  * Ribes  _et al._ (2019) P. Ribes, A. Ramos, L. Tolos, C. Gonzalez-Boquera, and M. Centelles, Astrophys. J. 883, 168 (2019), arXiv:1907.08583 [astro-ph.HE] .
  * Schürhoff  _et al._ (2010) T. Schürhoff, S. Schramm, and V. Dexheimer, Astrophys. J. Lett. 724, L74 (2010), arXiv:1008.0957 [astro-ph.SR] .
  * Marczenko  _et al._ (2022) M. Marczenko, K. Redlich, and C. Sasaki, Astrophys. J. Lett. 925, L23 (2022), arXiv:2110.11056 [nucl-th] .
  * Demorest  _et al._ (2010) P. Demorest, T. Pennucci, S. Ransom, M. Roberts, and J. Hessels, Nature 467, 1081 (2010), arXiv:1010.5788 [astro-ph.HE] .
  * Antoniadis  _et al._ (2013) J. Antoniadis _et al._ , Science 340, 6131 (2013), arXiv:1304.6875 [astro-ph.HE] .
  * Fonseca  _et al._ (2016) E. Fonseca _et al._ , Astrophys. J. 832, 167 (2016), arXiv:1603.00545 [astro-ph.HE] .
  * Cromartie  _et al._ (2019) H. T. Cromartie _et al._ , Nat. Astron. 4, 72 (2019), arXiv:1904.06759 [astro-ph.HE] .
  * Fonseca  _et al._ (2021) E. Fonseca _et al._ , Astrophys. J. Lett. 915, L12 (2021), arXiv:2104.00880 [astro-ph.HE] .
  * Abbott  _et al._ (2018) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 121, 161101 (2018), arXiv:1805.11581 [gr-qc] .
  * Abbott  _et al._ (2019) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. X 9, 011001 (2019), arXiv:1805.11579 [gr-qc] .
  * Riley  _et al._ (2019) T. E. Riley _et al._ , Astrophys. J. 887, L21 (2019), arXiv:1912.05702 [astro-ph.HE] .
  * Miller  _et al._ (2019) M. C. Miller _et al._ , Astrophys. J. 887, L24 (2019), arXiv:1912.05705 [astro-ph.HE] .
  * Miller  _et al._ (2021) M. C. Miller _et al._ , Astrophys. J. Lett. 918, L28 (2021), arXiv:2105.06979 [astro-ph.HE] .
  * Riley  _et al._ (2021) T. E. Riley _et al._ , Astrophys. J. Lett. 918, L27 (2021), arXiv:2105.06980 [astro-ph.HE] .
  * Aarts  _et al._ (2015) G. Aarts, C. Allton, S. Hands, B. Jäger, C. Praki, and J.-I. Skullerud, Phys. Rev. D92, 014503 (2015), arXiv:1502.03603 [hep-lat] .
  * Aarts  _et al._ (2017) G. Aarts, C. Allton, D. De Boni, S. Hands, B. Jäger, C. Praki, and J.-I. Skullerud, JHEP 06, 034 (2017), arXiv:1703.09246 [hep-lat] .
  * Aarts  _et al._ (2019) G. Aarts, C. Allton, D. De Boni, and B. Jäger, Phys. Rev. D99, 074503 (2019), arXiv:1812.07393 [hep-lat] .
  * De Tar and Kunihiro (1989) C. E. De Tar and T. Kunihiro, Phys. Rev. D39, 2805 (1989).
  * Jido  _et al._ (2000) D. Jido, T. Hatsuda, and T. Kunihiro, Phys. Rev. Lett. 84, 3252 (2000), arXiv:hep-ph/9910375 [hep-ph] .
  * Jido  _et al._ (2001) D. Jido, M. Oka, and A. Hosaka, Prog. Theor. Phys. 106, 873 (2001), arXiv:hep-ph/0110005 [hep-ph] .
  * Dexheimer  _et al._ (2008) V. Dexheimer, S. Schramm, and D. Zschiesche, Phys. Rev. C77, 025803 (2008), arXiv:0710.4192 [nucl-th] .
  * Gallas  _et al._ (2010) S. Gallas, F. Giacosa, and D. H. Rischke, Phys. Rev. D82, 014004 (2010), arXiv:0907.5084 [hep-ph] .
  * Paeng  _et al._ (2012) W.-G. Paeng, H. K. Lee, M. Rho, and C. Sasaki, Phys. Rev. D85, 054022 (2012), arXiv:1109.5431 [hep-ph] .
  * Sasaki  _et al._ (2011) C. Sasaki, H. K. Lee, W.-G. Paeng, and M. Rho, Phys. Rev. D84, 034011 (2011), arXiv:1103.0184 [hep-ph] .
  * Gallas  _et al._ (2011) S. Gallas, F. Giacosa, and G. Pagliara, Nucl. Phys. A 872, 13 (2011), arXiv:1105.5003 [hep-ph] .
  * Zschiesche  _et al._ (2007) D. Zschiesche, L. Tolos, J. Schaffner-Bielich, and R. D. Pisarski, Phys. Rev. C75, 055202 (2007), arXiv:nucl-th/0608044 [nucl-th] .
  * Benic  _et al._ (2015) S. Benic, I. Mishustin, and C. Sasaki, Phys. Rev. D91, 125034 (2015), arXiv:1502.05969 [hep-ph] .
  * Marczenko and Sasaki (2018) M. Marczenko and C. Sasaki, Phys. Rev. D97, 036011 (2018), arXiv:1711.05521 [hep-ph] .
  * Marczenko  _et al._ (2018) M. Marczenko, D. Blaschke, K. Redlich, and C. Sasaki, Phys. Rev. D98, 103021 (2018), arXiv:1805.06886 .
  * Marczenko  _et al._ (2019) M. Marczenko, D. Blaschke, K. Redlich, and C. Sasaki, Universe 5, 180 (2019), arXiv:1905.04974 [nucl-th] .
  * Marczenko (2020) M. Marczenko, Eur. Phys. J. ST 229, 3651 (2020), arXiv:2005.14535 [nucl-th] .
  * Marczenko  _et al._ (2020) M. Marczenko, D. Blaschke, K. Redlich, and C. Sasaki, Astron. Astrophys. 643, A82 (2020), arXiv:2004.09566 [astro-ph.HE] .
  * Marczenko  _et al._ (2021) M. Marczenko, K. Redlich, and C. Sasaki, Phys. Rev. D 103, 054035 (2021), arXiv:2012.15535 [hep-ph] .
  * Mukherjee  _et al._ (2017a) A. Mukherjee, S. Schramm, J. Steinheimer, and V. Dexheimer, Astron. Astrophys. 608, A110 (2017a), arXiv:1706.09191 [nucl-th] .
  * Mukherjee  _et al._ (2017b) A. Mukherjee, J. Steinheimer, and S. Schramm, Phys. Rev. C96, 025205 (2017b), arXiv:1611.10144 [nucl-th] .
  * Dexheimer  _et al._ (2013) V. Dexheimer, J. Steinheimer, R. Negreiros, and S. Schramm, Phys. Rev. C87, 015804 (2013), arXiv:1206.3086 [astro-ph.HE] .
  * Steinheimer  _et al._ (2011a) J. Steinheimer, S. Schramm, and H. Stocker, Phys. Rev. C84, 045208 (2011a), arXiv:1108.2596 [hep-ph] .
  * Weyrich  _et al._ (2015) J. Weyrich, N. Strodthoff, and L. von Smekal, Phys. Rev. C92, 015214 (2015), arXiv:1504.02697 [nucl-th] .
  * Sasaki and Mishustin (2010) C. Sasaki and I. Mishustin, Phys. Rev. C82, 035204 (2010), arXiv:1005.4811 [hep-ph] .
  * Yamazaki and Harada (2019a) T. Yamazaki and M. Harada, Phys. Rev. D 99, 034012 (2019a), arXiv:1809.02359 [hep-ph] .
  * Yamazaki and Harada (2019b) T. Yamazaki and M. Harada, Phys. Rev. C100, 025205 (2019b), arXiv:1901.02167 [nucl-th] .
  * Ishikawa  _et al._ (2019) T. Ishikawa, K. Nakayama, and K. Suzuki, Phys. Rev. D99, 054010 (2019), arXiv:1812.10964 [hep-ph] .
  * Steinheimer  _et al._ (2011b) J. Steinheimer, S. Schramm, and H. Stocker, J. Phys. G38, 035001 (2011b), arXiv:1009.5239 [hep-ph] .
  * Giacosa (2012) F. Giacosa, Prog. Part. Nucl. Phys. 67, 332 (2012), arXiv:1111.4944 [hep-ph] .
  * Motohiro  _et al._ (2015) Y. Motohiro, Y. Kim, and M. Harada, Phys. Rev. C92, 025201 (2015), [Erratum: Phys. Rev.C95,no.5,059903(2017)], arXiv:1505.00988 [nucl-th] .
  * Minamikawa  _et al._ (2021a) T. Minamikawa, T. Kojo, and M. Harada, Phys. Rev. C 103, 045205 (2021a), arXiv:2011.13684 [nucl-th] .
  * Minamikawa  _et al._ (2021b) T. Minamikawa, T. Kojo, and M. Harada, Phys. Rev. C 104, 065201 (2021b), arXiv:2107.14545 [nucl-th] .
  * Takeda  _et al._ (2018) Y. Takeda, Y. Kim, and M. Harada, Phys. Rev. C97, 065202 (2018), arXiv:1704.04357 [nucl-th] .
  * Walecka (1974) J. D. Walecka, Annals Phys. 83, 491 (1974).
  * Shlomo  _et al._ (2006) S. Shlomo, V. M. Kolomietz, and G. Colò, Eur. Phys. J A30, 23 (2006).
  * Zyla  _et al._ (2020) P. A. Zyla _et al._ (Particle Data Group), PTEP 2020, 083C01 (2020).
  * Essick  _et al._ (2021) R. Essick, I. Tews, P. Landry, and A. Schwenk, Phys. Rev. Lett. 127, 192701 (2021), arXiv:2102.10074 [nucl-th] .
  * Reinhard  _et al._ (2021) P.-G. Reinhard, X. Roca-Maza, and W. Nazarewicz, Phys. Rev. Lett. 127, 232501 (2021), arXiv:2105.15050 [nucl-th] .
  * Dexheimer  _et al._ (2019) V. Dexheimer, R. de Oliveira Gomes, S. Schramm, and H. Pais, J. Phys. G 46, 034002 (2019), arXiv:1810.06109 [nucl-th] .
  * Wehrberger  _et al._ (1989) K. Wehrberger, C. Bedau, and F. Beck, Nucl. Phys. A 504, 797 (1989).
  * Horikawa  _et al._ (1980) Y. Horikawa, M. Thies, and F. Lenz, Nucl. Phys. A 345, 386 (1980).
  * O’Connell and Sealock (1990) J. S. O’Connell and R. M. Sealock, Phys. Rev. C 42, 2290 (1990).
  * Lehr  _et al._ (2000) J. Lehr, M. Effenberger, and U. Mosel, Nucl. Phys. A 671, 503 (2000), arXiv:nucl-th/9907091 .
  * Nakamura  _et al._ (2010) S. X. Nakamura, T. Sato, T. S. H. Lee, B. Szczerbinska, and K. Kubodera, Phys. Rev. C 81, 035502 (2010), arXiv:0910.1057 [nucl-th] .
  * Danielewicz  _et al._ (2002) P. Danielewicz, R. Lacey, and W. G. Lynch, Science 298, 1592 (2002), arXiv:nucl-th/0208016 [nucl-th] .
  * Annala  _et al._ (2020) E. Annala, T. Gorda, A. Kurkela, J. Nättilä, and A. Vuorinen, Nature Phys. 16, 907 (2020), arXiv:1903.09121 [astro-ph.HE] .
  * Hebeler  _et al._ (2013) K. Hebeler, J. M. Lattimer, C. J. Pethick, and A. Schwenk, Astrophys. J. 773, 11 (2013), arXiv:1303.4662 [astro-ph.SR] .
  * Alvarez-Castillo  _et al._ (2016) D. Alvarez-Castillo, S. Benic, D. Blaschke, S. Han, and S. Typel, Eur. Phys. J. A52, 232 (2016), arXiv:1608.02425 [nucl-th] .
  * Tolman (1939) R. C. Tolman, Phys. Rev. 55, 364 (1939).
  * Oppenheimer and Volkoff (1939) J. R. Oppenheimer and G. M. Volkoff, Phys. Rev. 55, 374 (1939).
  * Klahn  _et al._ (2006) T. Klahn _et al._ , Phys. Rev. C 74, 035802 (2006), arXiv:nucl-th/0602038 .
  * Prakash  _et al._ (1992) M. Prakash, M. Prakash, J. M. Lattimer, and C. J. Pethick, Astrophys. J. Lett. 390, L77 (1992).
  * Hinderer (2008) T. Hinderer, Astrophys. J. 677, 1216 (2008), arXiv:0711.2420 [astro-ph] .
  * Damour and Nagar (2009) T. Damour and A. Nagar, Phys. Rev. D80, 084035 (2009), arXiv:0906.0096 [gr-qc] .
  * Binnington and Poisson (2009) T. Binnington and E. Poisson, Phys. Rev. D80, 084018 (2009), arXiv:0906.1366 [gr-qc] .
  * Yagi and Yunes (2013) K. Yagi and N. Yunes, Phys. Rev. D88, 023009 (2013), arXiv:1303.1528 [gr-qc] .
  * Hinderer  _et al._ (2010) T. Hinderer, B. D. Lackey, R. N. Lang, and J. S. Read, Phys. Rev. D81, 123016 (2010), arXiv:0911.3535 [astro-ph.HE] .
  * Maslov  _et al._ (2015) K. A. Maslov, E. E. Kolomeitsev, and D. N. Voskresensky, Phys. Lett. B 748, 369 (2015), arXiv:1504.02915 [astro-ph.HE] .
  * Somasundaram  _et al._ (2021) R. Somasundaram, I. Tews, and J. Margueron, (2021), arXiv:2112.08157 [nucl-th] .
  * Li  _et al._ (2021) J. J. Li, A. Sedrakian, and M. Alford, Phys. Rev. D 104, L121302 (2021), arXiv:2108.13071 [astro-ph.HE] .
  * Christian and Schaffner-Bielich (2021) J.-E. Christian and J. Schaffner-Bielich, (2021), arXiv:2109.04191 [astro-ph.HE] .
  * Sasaki (2018) C. Sasaki, Nucl. Phys. A970, 388 (2018), arXiv:1707.05081 [hep-ph] .


# Calculating $B$-meson decay constants using domain-wall light quarks and
nonperturbatively tuned relativistic $b$-quarks

Center for Computational Science, Boston University,
3 Cummington Mall, Boston, MA 02215, USA
E-mail

###### Abstract:

We calculate $B$-physics quantities using the RBC/UKQCD 2+1 flavor domain-wall
plus Iwasaki lattices and the relativistic heavy quark action developed by
Christ, Li and Lin. After tuning these parameters nonperturbatively, we
present our preliminary results for the calculation of the decay constants
$f_{B}$ and $f_{B_{s}}$ analyzing data at two lattice spacings of $a\approx
0.11$ fm and $a\approx 0.08$ fm.

##  1 Introduction

Using lattice QCD to compute the nonperturbative contributions to $B$-physics
quantities allows us to strengthen constraints on the CKM unitarity triangle
and determine CKM matrix elements. In the game of over-constraining the apex
of the CKM unitarity triangle, $B$-physics may also reveal signs of physics
beyond the Standard Model. The standard global fit of the CKM unitarity
triangle uses $B$-meson decay constants and mixing as an input [1, 2,
3].111For the latest results see: http://ckmfitter.in2p3.fr/,
http://utfit.roma1.infn.it/ and www.latticesaverages.org. Experiments such as
BaBar, Belle, and LHCb measure $B_{q}-\overline{B}_{q}$ mixing in terms of
mass differences (oscillation frequencies) $\Delta M_{q}$ to subpercent
accuracy [4]. In the Standard Model, these are parameterized by [5]

$\displaystyle\Delta
M_{q}=\frac{G_{F}^{2}m^{2}_{W}}{6\pi^{2}}\eta_{B}S_{0}M_{B_{q}}{f_{B_{q}}^{2}B_{B_{q}}}\lvert
V_{tq}^{*}V_{tb}\rvert^{2},$ (1)

where the index $q$ denotes a $d$\- or a $s$-quark and $M_{B_{q}}$ the mass of
the $B_{q}$-meson. Perturbatively computed inputs are the QCD coefficient
$\eta_{B}$ [5] and the Inami-Lim function $S_{0}$ [6]; contributions to be
computed nonperturbatively are the leptonic $B$-meson decay constant
$f_{B_{q}}$ and the bag parameter $B_{B_{q}}$. Precise theoretical
determinations of $f^{2}_{B_{q}}B_{B_{q}}$ are needed for precise extraction
of the product of CKM matrix elements $V^{*}_{tq}V_{tb}$. It is in particular
advantageous to compute the ratio [7]

$\displaystyle\frac{\Delta M_{s}}{\Delta
M_{d}}=\frac{M_{B_{s}}}{M_{B_{d}}}\,{\xi^{2}}\,\frac{\lvert
V_{ts}\rvert^{2}}{\lvert V_{td}\rvert^{2}},$ (2)

because statistical and systematic uncertainties largely cancel and the
nonperturbative contribution is contained in the $SU(3)$ breaking ratio

$\displaystyle\xi$
$\displaystyle=\frac{f_{B_{s}}\sqrt{B_{B_{s}}}}{f_{B_{d}}\sqrt{B_{B_{d}}}}.$
(3)

Further, constraints of the CKM triangle from $BR(B\to\tau\nu)$ require
precise knowledge of the decay constant $f_{B}$ [8, 9, 10]. Currently, only
HPQCD and Fermilab/MILC have published results for the decay constants $f_{B}$
and $f_{B_{s}}$ as well as $\xi$ using 2+1 dynamical flavor gauge field
configurations and ensembles with different lattice spacing. However, both
groups use MILC’s Asqtad lattices [11, 12, 13, 14] and hence an independent
cross-check is desirable.

Simulating $B$-physics on the lattice provides a special challenge because an
additional scale given by the large $b$-quark mass needs to be accommodated.
Here we report on our project to compute $B$-physics quantities using the
nonperturbatively tuned relativistic heavy quark action for $b$-quarks and
domain-wall fermions for the light $u,\,d,\,s$-quarks.

##  2 Computational setup

We base our $B$-physics project on gauge-field configurations generated by the
RBC/UKQCD collaborations with 2+1 flavors of domain-wall fermions and the
Iwasaki gluon action [15, 16]. The ensembles and their parameters are listed
in Tab. 1. Using the domain-wall action [17, 18] we generate on each
configuration six light valence-quark propagators with quark masses
$am_{\text{val}}^{24}$ = 0.005, 0.010, 0.020, 0.030, 0.0343 and $0.040$ on the
coarser $24^{3}$ ensembles and $am_{\text{val}}^{32}$ = 0.004, 0.006, 0.008,
0.025, 0.0272 and 0.030 on the finer $32^{3}$ ensembles. On the two “coarse”
$24^{3}$ ensembles ($a\approx 0.11$fm; $a^{-1}=1.729$ GeV) with light sea-
quark mass $am_{l}=0.005$ and $0.01$ we place one time source per
configuration, while on the “fine” $32^{3}$ ensembles ($a\approx 0.086$ fm;
$a^{-1}=2.281$ GeV) with light sea-quark masses $am_{l}=0.004,\,0.006$ and
$0.008$ we place two time sources per configuration separated by half the
temporal extent of the lattice. The mass of the light quark propagators with
$am^{24}_{\text{val}}=0.0343$ and $am^{32}_{\text{val}}=0.0272$ is extremely
close to the value of the physical strange quark mass on the corresponding
ensembles: $am_{s}^{24}=0.0348(11)$ and $am_{s}^{32}=0.0273(11)$.

For the $b$ quarks we use the anisotropic Sheikholeslami-Wohlert (clover)
action with the relativistic heavy-quark (RHQ) interpretation [19, 20]. We
require two experimental inputs for tuning the three parameters, $m_{0}a$,
$c_{P}$, $\zeta$, nonperturbatively. In our notation the RHQ action is defined
by

$\displaystyle
S=\sum_{n,n^{\prime}}\bar{\Psi}_{n}\left\\{\\!m_{0}+\gamma_{0}D_{0}-\\!\frac{aD_{0}^{2}}{2}+\zeta\left[\vec{\gamma}\cdot\vec{D}-\frac{a\left(\vec{D}\right)^{2}}{2}\right]\\!-a\sum_{\mu\nu}\frac{ic_{P}}{4}\sigma_{\mu\nu}F_{\mu\nu}\\!\right\\}_{\\!\\!n,n^{\prime}}\\!\\!\\!\\!\Psi_{n^{\prime}}.$
(4)

Table 1: Lattice simulation parameters used in our $B$-physics program. The columns list the lattice volume, approximate lattice spacing, light ($m_{l}$) and strange ($m_{h}$) sea-quark masses, unitary pion mass, and number of configurations and time sources analyzed. |  |  |  |  |  | # time
---|---|---|---|---|---|---
$\left(L/a\right)^{3}\times\left(T/a\right)$ |  $\approx a$(fm) |  $am_{l}$ |  $am_{h}$ |  $M_{\pi}$(MeV) | # configs. | sources
$24^{3}\times 64$ | 0.11 | 0.005 | 0.040 | 329 | 1636 | 1
$24^{3}\times 64$ | 0.11 | 0.010 | 0.040 | 422 | 1419 | 1
$32^{3}\times 64$ | 0.086 | 0.004 | 0.030 | 289 | 628 | 2
$32^{3}\times 64$ | 0.086 | 0.006 | 0.030 | 345 | 889 | 2
$32^{3}\times 64$ | 0.086 | 0.008 | 0.030 | 394 | 544 | 2

We tune the parameters by probing seven points of the parameter space (see
Fig. 1) measuring each time the spin averaged $B_{s}$-meson mass
$\overline{M}$, the hyperfine splitting $\Delta_{M}$ and the ratio of the rest
mass over the kinetic mass $M_{1}/M_{2}$. Assuming linear dependence of these
measured quantities on the parameters for the seven points used, we obtain the
tuned parameters by matching $\overline{M}$ and $\Delta_{M}$ to experimental
results and demanding that $M_{1}/M_{2}$ equals one. Typically this requires a
few iterations to ensure that the tuned parameters are given by an
interpolation in a region where the measured parameter dependence is
sufficiently linear. Further details are presented in [21] where we also
predict bottomonium masses and mass splittings as test of our tuned
parameters.

Figure 1: The seven RHQ parameter sets are indicated by the six endpoints and the central point of the molecule-like structure in the $m_{0}a$, $c_{P}$ and $\zeta$ parameter space. Table 2: Tuned parameters $m_{0}a$, $c_{P}$, $\zeta$ of the RHQ action for simulating $b$-quarks on the ensembles listed in Table I. For details see [21]. $\approx a(fm)$ | $m_{sea}^{l}$ | $m_{0}a$ | $c_{P}$ | $\zeta$
---|---|---|---|---
0.11 | 0.005 | 8.43(7) | 5.7(2) | 3.11(9)
0.11 | 0.010 | 8.47(9) | 5.8(2) | 3.1(2)
0.086 | 0.004 | 4.07(6) | 3.7(1) | 1.86(8)
0.086 | 0.006 | 3.97(5) | 3.5(1) | 1.94(6)
0.086 | 0.008 | 3.95(6) | 3.6(1) | 1.99(8)

The computation of the $B$-meson decay constants presented here is performed
similarly to the tuning, i.e. we compute decay constants at the same seven
points in parameter space (see Fig. 1) and then interpolate these measurements
to the value of the tuned RHQ parameters given in Tab. 2. We measure the decay
constant $f_{B}$ by computing on the lattice the decay amplitude $\Phi_{B}$
which is proportional to the vacuum-to-meson matrix element of the heavy-light
axial vector current ${\cal A}_{\mu}=\bar{b}\gamma_{5}\gamma_{\mu}q$

$\displaystyle\langle 0|{\cal
A}_{\mu}|B(p)\rangle/\sqrt{M_{B}}=ip^{\mu}\Phi_{B}^{(0)}/M_{B}$ (5)

with $p_{\mu}$ the four-momentum of the $B$-meson. To reduce discretization
errors we use an $O(a)$-improved axial vector current on the lattice

$\displaystyle\Phi_{B}^{\text{imp}}=\Phi_{B}^{(0)}+c_{1}\Phi_{B}^{(1)},$ (6)

where the coefficient $c_{1}$ is computed at 1-loop with mean-field improved
lattice perturbation theory [22]. Multiplying the decay amplitude by the
renormalization factor $Z_{\Phi}$, the lattice spacing $a$ and the mass of the
$B_{q}$-meson $M_{B_{q}}$ we obtain the decay constant

$\displaystyle
f_{B_{q}}=Z_{\Phi}\Phi_{B}^{\text{imp}}a^{-3/2}/\sqrt{M_{B_{q}}}.$ (7)

In our final setup we intend to compute the renormalization factor $Z_{\Phi}$
following the mostly nonperturbative method outlined in [23], i.e. $Z_{\Phi}$
is split into a nonperturbative part containing the flavor-conserving factors
$Z_{V}^{ll}$ and $Z_{V}^{bb}$ and a perturbatively computed factor
$\varrho_{bl}$ which is expected to be close to one and to have a more
convergent series expansion in $\alpha_{s}$:

$\displaystyle Z_{\Phi}=\varrho_{bl}\sqrt{Z_{V}^{bb}Z_{V}^{ll}}.$ (8)

Here we use the nonperturbatively obtained $Z_{V}^{ll}$ [16], while we compute
$Z_{V}^{bb}$ and $\varrho_{bl}$ at 1-loop in mean-field improved lattice
perturbation theory [22]. Work is in progress to obtain $Z_{V}^{bb}$
nonperturbatively on all ensembles [24].

##  3 Preliminary results

In Fig. 2 we show our preliminary results for the renormalized,
$O(a)$-improved decay amplitude
$\Phi_{B}^{\text{ren}}=Z_{\Phi}\Phi_{B}^{\text{imp}}$. While a chiral
extrapolation to the physical $d$-quark mass is needed to obtain
$\Phi_{B_{d}}$, we may take advantage of the simulation data close to the
physical strange quark mass to read off a rough estimate of $\Phi_{B_{s}}$.
This indicates our result is in agreement with the existing literature and we
expect statistical uncertainties of a few MeV. Figure 3 shows our results for
the SU(3) breaking ratio of $\Phi_{B_{s}}/\Phi_{B_{q}}$. The data exhibit a
weak dependence on the light sea-quark mass as well as a mild dependence on
the lattice spacing.

##  4 Outlook

Currently work is in progress on the nonperturbative computation of
$Z_{V}^{bb}$ as well as checks on the perturbative computation of the
correction factor $\varrho_{bl}$ and the $O(a)$-improvement coefficient
$c_{1}$. We are also working on combined, correlated fits using expressions
based on SU(3) and SU(2) heavy-light meson chiral perturbation theory [26, 27,
28, 29]. With the individual parts being in good shape we hope to publish
results soon. We will move on to the computation of $B_{q}-\overline{B}_{q}$
mixing matrix elements and their ratio $\xi$.

Figure 2: Renormalized and $O(a)$-improved decay amplitude $\Phi_{B_{q}}$
computed on the full data set and shown in $a_{32}$ lattice units. For better
visibility data points for $m_{\text{sea}}^{l}=0.01$ and $0.008$ are plotted
with a small, positive, horizontal offset; data points for
$m_{\text{sea}}^{l}=0.004$ have a small, negative, horizontal offset. The
vertical gray lines with error-band indicate the physical values of the $d$\-
and $s$-quark masses [25, 15, 16]. Figure 3: Ratio of decay amplitudes
$\Phi_{B_{s}}/\Phi_{B_{q}}$ computed on the full data set and shown in
$a_{32}$ lattice units. For better visibility data points for
$m_{\text{sea}}^{l}=0.01$ and $0.008$ are plotted with a small, positive,
horizontal offset; data points for $m_{\text{sea}}^{l}=0.004$ have a small,
negative, horizontal offset. The vertical gray lines with error-band indicate
the physical values of the $d$\- and $s$-quark masses [25, 15, 16]. The open
black triangle (square) marks the value of the valence strange-quark mass in
the simulation at which the plotted ratio is one by construction.

## Acknowledgments

We thank our colleagues of the RBC and UKQCD collaborations for useful help
and discussions. Numerical computations for this work utilized USQCD resources
at Fermilab, in part funded by the Office of Science of the U.S. Department of
Energy, as well as computers at Brookhaven National Laboratroy and Columbia
University. O.W. acknowledges support at Boston University by the U.S. DOE
grant DE-FC02-06ER41440.

## References

  * [1] J. Charles _et al._ (CKMfitter Group), Eur.Phys.J. C41, 1 (2005), arXiv:hep-ph/0406184
  * [2] M. Bona _et al._ (UTfit), JHEP 0507, 028 (2005), arXiv:hep-ph/0501199
  * [3] J. Laiho, E. Lunghi, and R. S. Van de Water, Phys.Rev. D81, 034503 (2010), arXiv:0910.2928 [hep-ph]
  * [4] J. Beringer _et al._ (Particle Data Group), Phys.Rev. D86, 010001 (2012)
  * [5] A. J. Buras, M. Jamin, and P. H. Weisz, Nucl. Phys. B347, 491 (1990)
  * [6] T. Inami and C. S. Lim, Prog. Theor. Phys. 65, 297 (1981)
  * [7] C. W. Bernard, T. Blum, and A. Soni, Phys.Rev. D58, 014501 (1998), arXiv:hep-lat/9801039
  * [8] E. Lunghi and A. Soni, Phys.Rev.Lett. 104, 251802 (2010), arXiv:0912.0002 [hep-ph]
  * [9] M. Bona _et al._ (UTfit), Phys.Lett. B687, 61 (2010), arXiv:0908.3470 [hep-ph]
  * [10] A. Lenz, U. Nierste, and CKMfitter Group, Phys.Rev. D83, 036004 (2011), arXiv:1008.1593 [hep-ph]
  * [11] E. Gamiz _et al._ (HPQCD), Phys. Rev. D80, 014503 (2009), arXiv:0902.1815 [hep-lat]
  * [12] H. Na _et al._ (HPQCD), Phys.Rev. D86, 034506 (2012), arXiv:1202.4914 [hep-lat]
  * [13] A. Bazavov _et al._ (Fermilab Lattice and MILC), Phys.Rev. D85, 114506 (2012), arXiv:1112.3051 [hep-lat]
  * [14] A. Bazavov _et al._ (Fermilab Lattice and MILC), Phys.Rev. D86, 034503 (2012), arXiv:1205.7013 [hep-lat]
  * [15] C. Allton _et al._ (RBC and UKQCD), Phys. Rev. D78, 114509 (2008), arXiv:0804.0473 [hep-lat]
  * [16] Y. Aoki _et al._ (RBC and UKQCD), Phys.Rev. D83, 074508 (2011), arXiv:1011.0892 [hep-lat]
  * [17] D. B. Kaplan, Phys. Lett. B288, 342 (1992), arXiv:hep-lat/9206013
  * [18] Y. Shamir, Nucl. Phys. B406, 90 (1993), arXiv:hep-lat/9303005
  * [19] N. H. Christ, M. Li, and H.-W. Lin, Phys.Rev. D76, 074505 (2007), arXiv:hep-lat/0608006
  * [20] H.-W. Lin and N. Christ, Phys.Rev. D76, 074506 (2007), arXiv:hep-lat/0608005
  * [21] Y. Aoki _et al._ (RBC and UKQCD)(2012), arXiv:1206.2554 [hep-lat]
  * [22] C. Lehner, PoS Lattice2012, 126 (2012)
  * [23] A. X. El-Khadra _et al._ , Phys.Rev. D64, 014502 (2001), arXiv:hep-ph/0101023
  * [24] T. Kawanai, R. S. Van de Water, and O. Witzel, PoS Lattice2012, 109 (2012), arXiv:1211.0956 [hep-lat]
  * [25] T. Blum _et al._ , Phys.Rev. D82, 094508 (2010), arXiv:1006.1311 [hep-lat]
  * [26] J. L. Goity, Phys. Rev. D46, 3929 (1992), arXiv:hep-ph/9206230
  * [27] D. Arndt and C. D. Lin, Phys.Rev. D70, 014503 (2004), arXiv:hep-lat/0403012
  * [28] C. Aubin and C. Bernard, Phys. Rev. D73, 014515 (2006), arXiv:hep-lat/0510088
  * [29] C. Albertus _et al._ (RBC and UKQCD), Phys.Rev. D82, 014505 (2010), arXiv:1001.2023 [hep-lat]