# Relativistic hybrid stars with sequential first-order phase transitions
in light of multimessenger constraints

Jia Jie Li School of Physical Science and Technology, Southwest University,
Chongqing 400715, China  Armen Sedrakian Frankfurt Institute for Advanced
Studies, D-60438 Frankfurt am Main, Germany  Institute of Theoretical Physics,
University of Wroclaw, 50-204 Wroclaw, Poland  Mark Alford Department of
Physics, Washington University, St. Louis, MO 63130, USA<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
###### Abstract

In this work, we consider the properties of compact stars in which quark
matter has low- and high-density phases that are separated by a first-order
phase transition. Thus, unlike the commonly considered case of a single phase
transition from hadronic to quark matter, our models of hybrid stars contain
sequential phase transitions from hadronic matter to low- and then to high-
density quark matter phases. We extend our previous study of the parameter
space of hybrid stars with a single phase transition to those with sequential
phase transitions, taking into account the constraints on the mass and radius
of neutron stars from the NICER experiment, the experimental inferences of the
neutron skin thickness of the lead nucleus by the PREX-II experiment, and
constraints on the tidal deformability from the gravitational-wave event
GW170817. We determine the range of the masses for which both twin and triplet
configurations, i.e., identical-mass stars with two and three different values
of radii, arise.

Neutron stars (1108); Neutron star cores (1107); Nuclear astrophysics (1129);
High energy astrophysics (739); Gravitational waves astronomy (675)

##  1 Introduction

The possibility of phase transition from hadronic to quark matter and its
implications for the properties of compact stars (CSs) has been the focus of
researchers since early work on the subject several decades ago (Ivanenko &
Kurdgelaidze, 1965; Itoh, 1970; Collins & Perry, 1975); for recent reviews,
see Alford et al. (2008); Anglani et al. (2014); Pisarski et al. (2019). An
interesting special case is the possibility of a first-order phase transition
between the hadronic and quark phases, which occurs when mixed phases are
disfavored by surface tension and electrostatic energy costs (Alford et al.,
2001). It has been established since that the first-order phase transition
leads to a softening of the equation of state (EoS) and lower maximum masses
of the sequences of hybrid stars compared to their purely hadronic
counterparts. In a previous paper (Li et al., 2021b) we studied the
implications of the recent radius determination of PSR J0740+6620 by the NICER
experiment and the determination of the neutron skin by the PREX-II experiment
(which gives information about symmetry energy and its slope) for the
structure of hybrid stars with a strong first-order phase transition from
nucleonic to quark matter. It was argued that if the interpretation of the
PREX-II experiment implies a stiff EoS, and hence a large radius for the
nucleonic branch, an early first-order phase transition may relax the tensions
with the astrophysical inferences of radii for the relevant ranges of masses
of CSs. We also deduced the ranges of mass and radius of twin stars which,
however, were restricted to a relatively narrow domain of masses and radii.

The purpose of this work is to implement the idea of sequential first-order
phase transitions in a hybrid star (Alford & Sedrakian, 2017), in particular
the case where there are two such transitions, which can be realized in
Nambu–Jona-Lasinio (NJL) models (e.g., Blaschke et al., 2010; Bonanno &
Sedrakian, 2012; Klähn et al., 2013). Here, we will continue using the scheme
developed in our previous work (Li et al., 2020a, 2021b) which accounts for
the currently available multimessenger information coming from various
channels of astronomical observations and terrestrial experiments. Let us
give, at this point, a brief list of the constraints against which our models
developed below will be tested:

  * •

PSR J0030+0451: This is the first object with mass and radius inferred to a
precision of around 10% (Riley et al., 2019; Miller et al., 2019). These
quantities were extracted from the fits to the data coming from the NICER
observatory, which required modeling the soft X-ray pulses produced by the
stellar rotation and the hot spots on the star’s surface and fitting to the
waveforms. The two (independent) analyses predict ($68\%$ credible interval)
$M=1.34^{+0.15}_{-0.16}M_{\odot}$, $R=12.71^{+1.14}_{-1.19}$ km (Riley et al.,
2019) and $M=1.44^{+0.15}_{-0.14}M_{\odot}$, $R=13.02^{+1.24}_{{}^{-}1.06}$ km
(Miller et al., 2019).

  * •

PSR J0740+6620: The mass of this pulsar was initially measured (Cromartie et
al., 2020) using Shapiro delay to be $2.08^{+0.07}_{-0.07}\,M_{\odot}$. Its
mass and radius were determined using the NICER X-ray light curves with the
results for the radius $12.39^{+1.30}_{-0.98}$ (Riley et al., 2021) and
$13.71^{+2.61}_{-1.50}$ km (Miller et al., 2021) and corresponding mass
estimates $2.07^{+0.07}_{-0.07}\,M_{\odot}$ and
$2.08^{+0.09}_{-0.09}\,M_{\odot}$ ($68\%$ credible interval).

  * •

GW170817: The tidal deformability (TD) of a star of mass $\sim 1.4\,M_{\odot}$
in the GW170817 event by the LIGO-Virgo Collaboration (Abbott et al., 2017,
2018, 2019) was constrained to be below the (dimensionless) value
$\Lambda_{1.4}\leq 580$, which implies a soft EoS at the intermediate (a few
times nuclear saturation) density range. This observation is consistent with a
phase transition from hadronic to quark matter, which can improve the
agreement of theoretical models with the TD extracted from GW170817 event
(e.g., Annala et al., 2018; Paschalidis et al., 2018; Most et al., 2018; Tews
et al., 2018; Burgio et al., 2018; Alvarez-Castillo et al., 2019; Christian et
al., 2019; Montana et al., 2019; Sieniawska et al., 2019; Essick et al., 2020;
Li et al., 2020a; Miao et al., 2020; Li et al., 2021a; Malfatti et al., 2020;
Rodriguez et al., 2021; Tan et al., 2022). Furthermore, the phase transition
may lead to the emergence of new branches of stable CSs, which provides a new
diagnostics of phase transition in CSs, through observation of twin and
triplet stars, i.e., one or two hybrid stars having the same mass, but
different radii from, a purely hadronic star (e.g., Alford et al., 2013;
Alford & Sedrakian, 2017; Paschalidis et al., 2018; Alvarez-Castillo et al.,
2019; Christian et al., 2019; Montana et al., 2019; Li et al., 2020a;
Christian & Schaffner-Bielich, 2022).

Due to the high mass of PSR J0740+6620, it probes the region of densities that
is highly relevant to a possible phase transition to quark matter. Recent work
on the EoS has generally concluded, independent of the details of the EoS and
methods of comparison with the multimessenger data adopted, that the EoS of
the star must be moderately soft at intermediate densities and stiff enough at
high densities. The first requirement accounts for the small TD in GW170817,
while the second one accounts for the large mass of PSR J0740+6620 (e.g., Li
et al., 2021b; Tan et al., 2022; Biswas, 2022; Legred et al., 2021;
Raaijmakers et al., 2021; Huth et al., 2022; Zhang & Li, 2021; Tang et al.,
2021; Christian & Schaffner-Bielich, 2022; Jokela et al., 2022; Drischler et
al., 2022; Contrera et al., 2022).

The isospin dependence of the nuclear interaction is constrained by the
measurements of the neutron skin of nuclei. Most recently, the Lead Radius
Experiment Collaboration (PREX-II) measured the neutron skin thickness of the
lead nucleus, $R^{208}_{\rm skin}=0.283\pm 0.071$ fm (mean and $1\,\sigma$
deviation), in a parity-violating electron-scattering experiment (Adhikari et
al., 2021). The existing theoretical analyses (Reed et al., 2021; Reinhard et
al., 2021) do not converge to a consistent mean value of the symmetry energy,
$E_{\rm sym}$, and its slope, $L_{\rm sym}$, at nuclear saturation density.
While the first reference (Reed et al., 2021) infers $E_{\rm sym}=38.1\pm 4.7$
MeV and $L_{\rm sym}=106\pm 37$ MeV using a family of relativistic density
functionals (DFs), the second reference (Reinhard et al., 2021) finds $E_{\rm
sym}=32\pm 1$ MeV and $L_{\rm sym}=54\pm 8$ MeV using a larger number of
relativistic and nonrelativistic DFs. (These two analyses are still consistent
with each other at better than $2\sigma$ accuracy.) They also included
additional constraints from the experimental limits on the dipole
polarizability of 208Pb, which prefer DFs predicting a small value of $L_{\rm
sym}$ (Reinhard et al., 2021). While the second set of parameter values is
within the standard range (Lattimer & Lim, 2013; Danielewicz & Lee, 2014;
Oertel et al., 2017; Baldo & Burgio, 2016), the first set is not. In
particular, their large value of $L_{\rm sym}$ is in tension with the GW170817
deformability measurement if one assumes a purely nucleonic composition (Reed
et al., 2021). Other authors (Essick et al., 2021) have confirmed the tension
between the results of Reed et al. (2021) and astrophysical data using a
nonparametric EoS. Thus, following our previous work (Li et al., 2021b), we
will consider a broader range of $L_{\rm sym}$ to cover the possibilities
claimed in the work above.

The physics of dense quark matter allows for multiple quark phases with
distinct properties. One standard scenario is two-flavor color-superconducting
(2SC) at medium densities and the color-flavor-locked (CFL) phase at high
density. Thus, one of the aims of this work is to expand on the previous study
(Li et al., 2021b), which was restricted to twin configurations, by
considering triplets—three stars with the same mass but different radii—which
can arise if there are three disconnected branches in the mass-radius
(hereafter $M$-$R$) diagram (Alford & Sedrakian, 2017; Li et al., 2020a).

The paper is organized as follows. In Section 2 we briefly define the EoS that
we use to describe the hadronic and quark phases. In Section 3 we assess the
existence of twin or triplet configurations and confront the resultant EoS
models with the inferences from GW170817 and NICER observations of PSR
J0740+6620. Our conclusions are given in Section 4.

##  2 Construction of equation of state and basic picture

###  2.1 Nuclear matter equation of state

Figure 1:  $M$-$R\ $relation for nucleonic EoS within different pairs of
values of $Q_{\rm sat}$ and $L_{\rm sym}$ (in MeV). We show three ranges of
$M$-$R\ $curves, for $Q_{\text{sat}}=900\,\text{MeV}$ (gray),
$-100\,\text{MeV}$ (red) and $-500\,\text{MeV}$ (blue). For each, $L_{\rm
sym}$ is varied from $45\,\text{MeV}$ to $105\,\text{MeV}$. Thick solid lines
show combinations of ($Q_{\rm sat},L_{\rm sym}$) that could marginally meet
the GW170817 constraint, i.e., (900, 45), (-100, 65), and (-500, 85). In
addition, we show the TDs of CSs with mass $M/M_{\odot}=1.3,1.7$, and 2.0 for
two extreme models. Constraints at 90% credibility from multimessenger
astronomy are shown by shaded regions (Abbott et al., 2019; Riley et al.,
2019, 2021; Miller et al., 2019, 2021); see text for details.

We will adopt for the hadronic phase the same description as in the
predecessor paper (Li et al., 2021b). We provide below a brief account of our
approach for the sake of completeness. The hadronic matter is described within
a covariant density-functional (CDF) approach with density-dependent nucleon-
meson couplings (Lalazissis et al., 2005). The density dependence of the
coupling allows us to establish a one-to-one correspondence between our CDF
and the purely phenomenological expansion of the energy density of nuclear
matter (e.g., Margueron et al., 2018) in the vicinity of the saturation
density, $\rho_{\rm sat}$, with respect to the number density $\rho$ and
isospin asymmetry $\delta=(\rho_{\rm n}-\rho_{\rm p})/\rho$ where $\rho_{\rm
n(p)}$ is the neutron (proton) number density, as

$\displaystyle E(\chi,\delta)$ $\displaystyle\simeq$ $\displaystyle
E_{\rm{sat}}+\frac{1}{2!}K_{\rm{sat}}\chi^{2}+\frac{1}{3!}Q_{\rm{sat}}\chi^{3}$
(1)
$\displaystyle+\,E_{\rm{sym}}\delta^{2}+L_{\rm{sym}}\delta^{2}\chi+{\mathcal{O}}(\chi^{4},\chi^{2}\delta^{2}),$

where $\chi=(\rho-\rho_{\rm{sat}})/3\rho_{\rm{sat}}$. The coefficients of this
double expansion are referred to commonly as incompressibility, $K_{\rm sat}$,
the skewness, $Q_{\rm{sat}}$, the symmetry energy, $E_{\rm{sym}}$, and its
slope parameter, $L_{\rm{sym}}$. The mapping between the CDF and the
phenomenological expansion (1) allows us to express the gross properties of
CSs in terms of physically transparent quantities.

In this work, we use three sets of representative EoS models, taking three
values of $Q_{\rm sat}=-500,-100$, and 900 MeV, and exploring values of
$L_{\rm sym}$ ranging from 45 to $105$ MeV. Larger values of $L_{\rm sym}$
correspond to a stiffer EoS near nuclear saturation density, leading to larger
radii for a $1.4\,M_{\odot}$ star. Larger values of $Q_{\rm sat}$ mean that
the EoS is stiffer at high density, thereby increasing the maximum mass of a
static nucleonic CS (e.g., Li & Sedrakian, 2019; Zhang et al., 2018; Margueron
et al., 2018). For $Q_{\rm sat}=-500$ MeV the maximum mass is about
$2.1\,M_{\odot}$, which matches the mass measurement of PSR J0740+6620
(Cromartie et al., 2020; Fonseca et al., 2021); for $Q_{\rm sat}=-100$ MeV the
maximum mass is consistent with the (approximate) upper limit on the maximum
mass of static CSs $\sim 2.3\,M_{\odot}$ inferred from the analysis of the
GW170817 event (Rezzolla et al., 2018; Khadkikar et al., 2021); finally, for
$Q_{\rm sat}=900$ MeV the maximum mass is close to $2.5\,M_{\odot}$, which
would be compatible with the mass of the secondary in the GW190814 event
(Abbott et al., 2020b) and its interpretation as a nucleonic CS (Fattoyev et
al., 2020; Sedrakian et al., 2020; Li et al., 2020b). We choose the range of
$L_{\rm sym}$ between the central value and the lower limit of the 90%
credible interval (CI) of the PREX-II measurement (Adhikari et al., 2021; Reed
et al., 2021).

The $M$-$R\ $relations for our nucleonic EoS models are shown in Figure 1,
along with the current astrophysical observational constraints. These include
(i) the ellipses obtained by the two NICER modeling groups for PSR J0030+0451
and J0740+6620 (Riley et al., 2019, 2021; Miller et al., 2019, 2021); (ii) the
regions for each of the two CSs that merged in the gravitational-wave (GW)
event GW170817 (Abbott et al., 2019); and (iii) the mass of the secondary
component of GW190814 (Abbott et al., 2020b). All the regions/limits are given
at 90% CI. As seen from Figure 1, the softness of the EoS (as implied by the
GW170817 event) and stiffness at low and intermediate densities implied by the
large value of $L_{\rm sym}$ suggested by one of the analyses (Reed et al.,
2021) of the PREX-II experiment can be reconciled by an appropriate choice of
the parameters. Indeed, this could be accomplished by nucleonic EoS that trade
stiffness at high density for softness at low density, e.g., $L_{\rm
sym}\lesssim 85$ MeV for $Q_{\rm sat}\sim-500$ MeV and $L_{\rm sym}\lesssim
45$ MeV (which is at the lower end of the 90% CI of Reed et al. (2021)) for
$Q_{\rm sat}\sim 900$ MeV. Thus, we conclude that the nucleonic models are not
inconsistent with the current information available from multimessenger
astrophysics if fairly low values of $L_{\rm sym}$ are adopted. We next
explore how this situation changes when sequential first-order phase
transitions are allowed.

###  2.2 Quark matter equation of state

We will model below the EoS of the quark phase using a synthetic constant-
sound-speed (CSS) parameterization (Zdunik & Haensel, 2013; Alford et al.,
2013), which matches well with the predictions based on the NJL model
computations that include vector repulsion (e.g., Blaschke et al., 2010;
Bonanno & Sedrakian, 2012; Klähn et al., 2013). The extension of the CSS EoS
to the case of two sequential phase transitions is given by (Alford &
Sedrakian, 2017),

$\displaystyle
p(\varepsilon)=\left\\{\begin{array}[]{ll}p_{1},&\quad\varepsilon_{1}<\varepsilon<\varepsilon_{1}\\!+\\!\Delta\varepsilon_{1}\\\\[2.15277pt]
p_{1}+s_{1}\bigl{[}\varepsilon-(\varepsilon_{1}\\!+\\!\Delta\varepsilon_{1})\bigr{]},&\quad\varepsilon_{1}\\!+\\!\Delta\varepsilon_{1}<\varepsilon<\varepsilon_{2}\\\\[2.15277pt]
p_{2},&\quad\varepsilon_{2}<\varepsilon<\varepsilon_{2}\\!+\\!\Delta\varepsilon_{2}\\\\[2.15277pt]
p_{2}+s_{2}\bigl{[}\varepsilon-(\varepsilon_{2}\\!+\\!\Delta\varepsilon_{2})\bigr{]},&\quad\varepsilon>\varepsilon_{2}\\!+\\!\Delta\varepsilon_{2}\end{array}\right.$
(6)

where $p_{1,2}$ and $\varepsilon_{1,2}$ are the pressure and energy density at
which the transition from hadronic to quark matter and from low-density quark
phase (hereafter Q1) to high-density quark phase (hereafter Q2) takes place.
We recall that the last phase transition within the quark phase could be the
transition between the 2SC to CFL phases, but other options are not excluded.
Finally, $s_{1}$ and $s_{2}$ are squared sound speeds in phases Q1 and Q2,
which we quote below in units of the speed of light.

Figure 2:  Illustrative EoS and $M$-$R\ $relation for hybrid models. Schematic
plot showing the parameterizations of the EoS with single (light green) and
double (dark green) phase transitions that could predict the same maximum mass
in the quark branch. The emergence of triplet configurations that are
characterized by the minimal or maximum masses in each branch is shown in the
inset, with one example that obeys conditions (8) and one that does not.

Figure 2 provides a schematic picture of the double phase transition according
to equation (6). The six independent parameters which enter Equation (6) are

$\displaystyle\varepsilon_{1},\quad\Delta{\varepsilon_{1}},\quad\varepsilon_{2},\quad\Delta{\varepsilon_{2}},\quad
s_{1},\quad s_{2}\ .$ (7)

For a large enough jump in energy density, $\Delta\varepsilon_{1}$, when the
central pressure of a star rises above $p_{1}$ and Q1 quark matter appears in
the core, the star becomes unstable. 111Here and below we use the standard
stability criterion, which implies that CSs are unstable on the descending
branch of $M$-$R\ $diagram and stable on the ascending one. This is not always
the case for some boundary conditions on the interface between quark and
nuclear matter (see, e.g., Pereira et al., 2018; Gonçalves & Lazzari, 2022;
Rau & Sedrakian, 2022). However, it is also possible to regain stability at
higher central pressures: for a certain range of values of $\varepsilon_{1}$,
$\Delta\varepsilon_{1}$ and $s_{1}$, there can be a second stable branch or
“third family” of CSs. In this case, twin configurations may appear, i.e., two
stable CSs may have the same mass but different radii and, consequently,
deformabilities. If a second phase transition in the quark phase takes place,
then a third stable branch (or “fourth family”) of CSs containing Q2 quark
matter in the core can arise. In analogy to the above, for suitably chosen
parameters, triplet configurations can arise (Alford & Sedrakian, 2017; Li et
al., 2020a), in which case there are three stable CSs all having the same
masses but different radii as well as deformabilities.

For the convenience of subsequent discussion, let us define the maximal masses
of the purely hadronic star (“N star”), the star with a Q1 core (“Q1 star”),
and the star with a Q2 inner core and Q1 outer core (“Q2 star”), as

$M^{\rm{N}}_{\rm{max}},\quad M^{\rm{Q1}}_{\rm{max}},\quad
M^{\rm{Q2}}_{\rm{max}}\ .$

It is also useful to define the minimum values of the masses of the two
stellar branches which consist of stars with only Q1 in the core and with Q1
and Q2 in the core, respectively, as

$M^{\rm{Q1}}_{\rm{min}},\quad M^{\rm{Q2}}_{\rm{min}}\ .$

See Figure 2 for an illustration of these parameters. We show also in the
inset of this figure schematic $M$-$R\ $relations, which illustrate cases
where triplets of stars arise, with each new phase of matter introducing a new
family of CSs. The shaded regions in the inset define the range of masses for
which triplet configurations arise.

Figure 3:  $M$-$R\ $relations showing the stellar sequences have twin or
triplet configurations in our setup. The inset illustrates the case with a
fixed $s_{2}$ value but varying the values of $s_{1}$. For all curves the
values of $s_{1},s_{2}$ are as indicated in the plot. The dotted thin lines
indicate unstable configurations.

To study the possible role of quark phases in the context of the NICER results
for the masses, radii, and GW inferences for the TDs, we reduce the six-
dimensional parameter space (7) as follows:
(1) We work in terms of two physical parameters: the maximum mass on the
nucleonic branch, $M^{\rm{N}}_{\rm{max}}$ (which determines
$\varepsilon_{1}$), and the maximum mass on the hybrid branch,
$M^{\rm{Q2}}_{\rm{max}}$. (Note that in the case of a single phase transition
we adopt the convention of calling the quark phase “Q2”.)
(2) As in Li et al. (2020a), we impose the following conditions (discussed
below):

$\displaystyle M^{\rm{N}}_{\rm{max}}=M^{\rm{Q1}}_{\rm{max}},\quad
M^{\rm{Q1}}_{\rm{min}}=M^{\rm{Q2}}_{\rm{min}},\quad
M^{\rm{Q2}}_{\rm{max}}\geqslant M^{\rm{Q1}}_{\rm{max}}\ .$ (8)

(3) We require the Q2 branch to pass through both the PSR J0740+6620 and
J0030+0451 ellipses. This tightly constrains $s_{2}$: we find that its value
is highly correlated with $M^{\rm{Q2}}_{\rm{max}}$.
(4) We note that the value of $s_{1}$ has a negligible influence on the
$M$-$R\ $curve because the Q1 branch is very short.
Using the four constraints/correlations described above, we reduce the six-
parameter space (7) to a space with two physical parameters,
$M^{\rm{N}}_{\rm{max}}$ and $M^{\rm{Q2}}_{\rm{max}}$, which specify a nearly
unique $M$-$R\ $curve.

The conditions (8) ensure that, if in some range of masses twins of the
hadronic stars exist due to the phase transition to the Q1 phase, the second
phase transition leads to triplets in the same range of masses; see the inset
of Figure 2 for an illustration. In this inset, we show one example that obeys
conditions (8) and one that does not, where both models have the same values
for $M^{\rm{N}}_{\rm{max}}$ and $M^{\rm{Q2}}_{\rm{max}}$. The inset
illustrates the behavior that we typically find: if $M^{\rm Q1}_{\rm max}$
rises above $M^{\rm N}_{\rm max}$, achieved by increasing $\varepsilon_{2}$
(namely the width of the Q1 phase) then at fixed $M^{\rm Q2}_{\rm max}$ we
need to reduce the Q1-Q2 energy-density jump $\Delta\varepsilon_{2}$ so that
$M^{\rm Q2}_{\rm min}$ rises above $M^{\rm Q1}_{\rm min}$, which results in a
smaller mass range for triplets.

In conclusion, the constraint (8) helps to ensure that triplets exist over a
reasonable range of masses, $\Delta M_{\rm triplet}$, for a hybrid EoS model
with given values of $M^{\rm{N}}_{\rm{max}}$ and $M^{\rm{Q2}}_{\rm{max}}$.

Figure 3 illustrates how transitions to and in quark matter allow a nucleonic
EoS that is stiff at low densities to be consistent with astrophysical
constraints. In this example, the nucleonic EoS has $L_{\text{sym}}=105$ MeV
(stiff at low density) and $Q_{\text{sat}}=900$ MeV (stiff at high density)
and we require $M^{\rm Q2}_{\rm max}=2.1\,M_{\odot}$ and $M^{\rm N}_{\rm
max}=1.3\,M_{\odot}$, which corresponds to a first-order transition at
$\rho_{\text{tran}}=1.89\,\rho_{\rm sat}$. This creates a hybrid branch that
is compatible with the GW170817 and PSR J0740+6620 constraints. As an
illustration of the limited role of $s_{1}$, the insert of Figure 3 shows that
varying $s_{1}$ for a fixed $s_{2}$ changes the mass by $0.001\,M_{\odot}$.
Note that when we change $s_{1}$, the parameters $\Delta\varepsilon_{1}$,
$\varepsilon_{2}$ and $\Delta\varepsilon_{2}$ have to be accordingly adjusted,
in order to fulfill the first two conditions in Equation (8).

##  3 Existence of twin or triplet configurations

To assess the existence of twin/triplet configurations, we construct EoS from
the parameter space of our model that allow for mass twins or triplets that
are consistent with both NICER and GW measurements. Specifically, they yield
radii that are just above the 90% CI lower limit for PSR J0740+6620 from
NICER, and just below the 90% CI upper limit on the radius of a
$1.36\,M_{\odot}$ star. 222Note that $1.36\,M_{\odot}$ is the mass value
inferred for an equal-mass binary in the GW170817 event. In Section 3.1 we
will study the stiffer ones that require a transition to quark matter if they
are to obey the GW170817 constraint, and in Sections 3.2 and 3.3 we will move
on the softer EoS that do not require such a transition.

Figure 4:  Illustrative $M$-$R\ $relations for hybrid EoS models with stiff
nucleonic EoS with single (dashed lines) or double (solid lines) phase
translations at subcanonical masses of star sequences. We use nucleonic EoS
(with $Q_{\rm sat}=900,L_{\rm sym}=105$ MeV) and fix $M^{\rm N}_{\rm
max}=1.3\,M_{\odot}$ and $M^{\rm Q2}_{\rm max}=2.0,2.1\,M_{\odot}$. The
emergence of subcanonical-mass triplet configurations is shown in the inset.
For all the hybrid branches, the values of $s_{1},s_{2}$ are as indicated in
the plot.

###  3.1 Hybrid equation of state featuring sub-canonical-mass twins or
triplets

Nucleonic EoS models that are stiff at low and high density (with large
$L_{\rm sym}$, in the 100 MeV range, as suggested by Reed et al. (2021), and
positive values of $Q_{\rm sat}$) predict a large radius for a canonical-mass
CS that is in tension with the GW170817 inference; see Figure 1. We can
resolve this by positing a first-order phase transition at low density (i.e.,
in the hadronic branch ends below a mass of $1.4\,M_{\odot}$) leading to a
hybrid branch of more compact stars. It should be mentioned that the same
effect can be achieved in the EoS models where heavy baryons like $\Delta$
resonances appear early. This leads to softening of the EoS at intermediate
densities, reducing the radius and TD of the star (see, e.g., Li et al.,
2020a; Sedrakian et al., 2021, and references therein).

We have shown previously (Li et al., 2021b) that twin configurations can be
expected in the mass interval of $1.2$-$1.4\,M_{\odot}$. Having this in mind,
we fix $M^{\rm N}_{\rm max}=1.3\,M_{\odot}$, which corresponds to
$\rho_{\text{tran}}\thickapprox 2\,\rho_{\text{sat}}$, and study the general
characteristics of the resulting subcanonical-mass twins or triplets.

Figure 5:  Ranges of parameters characterizing twin and triplet configurations
for hybrid EoS models. Panel (a) shows the mass ranges of twins $\Delta M_{\rm
twin}$ and triplet $\Delta M_{\rm triplet}$; panel (b) shows the difference in
the radius $\Delta R_{1.3}$, and TD $\Delta\Lambda_{1.3}$ for twins (in green)
and triplets (in blue) with mass $M=1.3\,M_{\odot}$, as functions of the
parameters of the nucleonic EoS. In panel (a) those models with $\Delta
R_{1.3}\geq 1.0$ km for Q2-N pairs of stars are marked by bold lines. In panel
(b), for triplets, the column height shows the value for Q2-N pair of stars,
while the values for Q2-Q1 and Q1-N pairs can be read off by the lighter and
darker colored column segments, respectively. The colored numbers on the
x-axis mark those nucleonic models that meet GW170817’s inference.

Figure 4 shows exploratory examples based on the stiffest nucleonic EoS in our
collection with $Q_{\rm sat}=900$, $L_{\rm sym}=105$ MeV. We show $M$-$R$
curves with two values for the maximum mass on the Q2 branch, $M^{\rm Q2}_{\rm
max}=2.0$ and $2.1\,M_{\odot}$, and explore several values of $s_{1}$ and
$s_{2}$. It is seen that the NICER result for PSR J0740+6620 puts a strong
constraint on $M^{\rm Q2}_{\rm max}$ and $s_{2}$. For $M^{\rm
Q2}_{\text{max}}=2.1\,M_{\odot}$ it requires $s_{2}\approx 1.0$ and for
$M^{\rm Q2}_{\text{max}}=2.0\,M_{\odot}$ it requires $s_{2}\approx 0.7$. The
resultant hybrid stars are also compatible with the constraints from GW170817
and the NICER results for PSR J0030+0451.

Figure 5 summarizes the mass ranges of twins, $\Delta M_{\rm twin}$, and
triplets, $\Delta M_{\rm triplet}$, defined as the range between the maximum
value of the nucleonic star mass, $M^{\rm N}_{\rm max}$, and the common
minimum mass of the hybrid star branches, $M^{\rm Q1}_{\rm min}=M^{\rm
Q2}_{\rm min}$. It also gives the difference in the radius, $\Delta R_{1.3}$,
and TD, $\Delta\Lambda_{1.3}$, for stars with mass $M=1.3\,M_{\odot}$, as
functions of the parameters defining the hybrid EoS model.

From Figures 4 and 5, we observe the systematic features that were established
in the case of only twin stars in Li et al. (2021b).

First, increasing $M^{\rm Q2}_{\rm max}$ at fixed $s_{2}$ or decreasing
$s_{2}$ at fixed $M^{\rm Q2}_{\rm max}$ tends to reduce the mass range of
triplets and twins. At the same time, the instability region of mass (or
radius) between nucleonic and hybrid stars decreases as the value of $M^{\rm
Q2}_{\rm max}$ increases at fixed $s_{2}$, or as $s_{2}$ decreases at fixed
$M^{\rm Q2}_{\rm max}$. A smaller value of $s_{2}$, in the current setup,
allows a smaller $\Delta\varepsilon_{2}$ and for a resulting steep decrease in
the mass range for twins and/or triplets. For our parameter choice, in order
for the $M$-$R\ $curves to pass through the NICER 90% CI region for PSR
J0740+6620, it requires $s_{2}$ close to 1.0.

Second, for stiffer nucleonic EoS featuring larger values of $L_{\rm sym}$
and/or $Q_{\rm sat}$, the range of masses where twin/triplet stars appear is
larger. This is because for hybrid EoS with fixed $M^{\rm Q2}_{\rm max}$,
increased stiffness of the nucleonic EoS must be offset by a softer quark
phase, e.g., larger energy-density jumps $\Delta{\varepsilon_{1}}$,
$\Delta{\varepsilon_{2}}$, or width
$\varepsilon_{2}-\varepsilon_{1}-\Delta{\varepsilon_{1}}$. This makes the
hybrid branch more compact while the stiffer nucleonic EoS makes the hadronic
branch less compact, i.e., the difference in radius is larger. The mass ranges
of triplet stars $\Delta M_{\rm triplet}$ are typically $\sim
10^{-3}\,M_{\odot}$, which is one order smaller than the mass range for twin
stars, $\Delta M_{\rm twin}$.

Third, the differences in radii and deformabilities for triplet stars are
typically 40-50% smaller than for twin stars. The values of $\Delta R_{1.3}$
for N-Q2 pairs are, at most, about 1.0 km; the values for N-Q1 or Q1-Q2 pairs
are even smaller. Directly observing twins or triplets would therefore require
a radius measurement accuracy of better than 1 km. This corresponds to a
radius measurement accuracy of less than 10%, which is not yet available, but
NICER aims to achieve 5% accuracy in the future. 333NICER Home:
https://heasarc.gsfc.nasa.gov/docs/nicer/ The difference in the TDs of
triplets $\Delta\Lambda_{1.3}$ can be several hundred to one thousand, so,
given TD measurements with error bars in the few hundred range, there is a
possibility that TD measurements from the inspiral phase of CS mergers could
identify triplets.

Figure 6:  TDs of compact objects with single (in blue) and double (in red)
phase transitions for a fixed value of binary chirp mass
$\mathcal{M}=1.186\,M_{\odot}$ inferred from the GW170817 event, varying the
mass ratio. The hybrid models are constructed from nucleonic EoS by setting
$M^{\rm N}_{\rm max}=1.3\,M_{\odot}$ of the hadronic branch, $M^{\rm Q1}_{\rm
max}=1.3\,M_{\odot}$, $s_{1}=0.7$ of the Q1 branch, and $M^{\rm Q2}_{\rm
max}=2.1\,M_{\odot}$, $s_{2}=1.0$ of the Q2 branch. For models with mass
twins, the two types of pairs for stars with masses $M_{1}$ and $M_{2}$ are
Q2-N and Q2-Q1, while for models with mass triplets the three types of pairs
are Q2-N, Q2-Q1, and Q2-Q2. The shaded regions correspond to the 50% and 90%
CIs taken from the analysis of GW170817 within the PhenomPNRT model (Abbott et
al., 2018).

We now compare our theoretical TDs for hybrid star models with the
observational constraints for this quantity obtained from the analysis of the
GW170817 event (Abbott et al., 2018). We take the chirp mass as
$\mathcal{M}=1.186\,M_{\odot}$ inferred from this merger and carry out the
comparison only using the analysis which assumes the (more plausible) low-spin
case (Abbott et al., 2018). For this binary, the component masses are found to
be in the range $1.16$-$1.60\,{M}_{\odot}$ at 90% CI. In the scenarios that we
are exploring in this section, this implies that at least one of the
components should be a hybrid star.

Figure 6 displays the TDs $\Lambda_{1}$ and $\Lambda_{2}$ of the stars
involved in the binary with masses $M_{1}$ (the primary, which is defined as
the heavier of the pair) and $M_{2}$ (secondary) in the cases of single and
double phase transitions. We show curves for two models, one with a single Q2
phase and one with both Q1 and Q2 phases. The models are constructed from
soft-stiff nucleonic EoS ($Q_{\rm sat}=-100$ and $L_{\rm sym}=105$ MeV) by
setting $M^{\rm N}_{\rm max}=1.3\,M_{\odot}$ of the hadronic branch, and
combined quark EoS featuring $M^{\rm Q1}_{\rm max}=1.3\,M_{\odot}$,
$s_{1}=0.7$ of the Q1 branch, and $M^{\rm Q2}_{\rm max}=2.1\,M_{\odot}$,
$s_{2}=1.0$ of the Q2 branch. The diagonal line corresponds to the case of an
equal-mass binary with $M_{1,2}=1.362\,M_{\odot}$. The shaded areas correspond
to the 90% and 50% CIs, which are inferred from the analysis of the GW170817
event using the PhenomPNRT waveform model (Abbott et al., 2019).

The hypothesis of a stiff nucleonic EoS with a low-density transition to quark
matter is compatible with limits on the TD from GW170817. In Figure 6 we show
how this is possible: for the EoS models plotted there, both stars in the
GW170817 merger could be hybrid, either Q2-Q1 or Q2-Q2, or one could be a
heavy hybrid star and the other could be a nucleonic star near the top of the
nucleonic branch.

###  3.2 Hybrid equation of state featuring high-mass twins or triplets

Figure 7:  Same as in Figure 4, but for hybrid EoS models with an
intermediate-soft nucleonic EoS that allows single (dashed lines) or double
(solid lines) phase transitions to appear in high-mass stars. The results are
constructed from a nucleonic EoS ($Q_{\rm sat}=900,L_{\rm sym}=45$ MeV) by
fixing $M^{\rm N}_{\rm max}/M_{\odot}=1.7$ of the nucleonic branch, while
varying the maximum mass $M^{\rm Q2}_{\rm max}/M_{\odot}=2.0,2.1$ of the Q2
branch.

Consider next a class of less stiff nucleonic EoS models that could
(approximately) match GW170817’s inference without a phase transition to quark
matter. Figure 7 shows examples of the nucleonic EoS with $Q_{\rm sat}=900$,
$L_{\rm sym}=45$ MeV for a fixed value of $M^{\rm N}_{\rm
max}=1.70\,M_{\odot}$, which corresponds to
$\rho_{\text{tran}}=2.41\,\rho_{\rm sat}$. The values of the remaining
parameters — the sound speed squares, $s_{1}$ and $s_{2}$, and the maximum
mass, $M^{\rm Q2}_{\rm max}$ — are chosen such that the Q2 branch is located
close to the lower bound of the NICER’s 90% CI for the radius of PSR
J0740+6620. Within this setup, we find the acceptable models are those which
have the values of the two parameters $M^{\rm Q2}_{\rm max}/M_{\odot}$ and
$s_{2}$ defined as pairs, (2.1, 1.0) or (2.0, 0.8), which cover the
approximate mass ranges for twins and triplets.

In Figure 8 we show the ranges of parameters characterizing twin and triplet
configurations by varying continuously $L_{\rm sym}$ for some fixed values of
$Q_{\rm sat}$. Note that this figure includes for completeness also the
results for stiff nucleonic EoS that predict a $M$-$R\ $range outside of the
GW170817 ellipses.

Figure 8:  Same as in Figure 5, but for hybrid EoS models with $M^{\rm N}_{\rm
max}/M_{\odot}=1.7$. In panel (a) those models with $\Delta R_{1.7}\geq 1.0$
km for Q2-N pairs of stars are marked by bold lines. Notice that there are
almost no twin/triplet configurations for a nucleonic EoS models with $Q_{\rm
sat}=-500$ and $-100$ MeV.

The general features found for the case of hybrid models featuring
subcanonical-mass twin/triplet configurations above are replicated within this
class as well. However, we find (almost) no twin/triplet solutions for
nucleonic models with $Q_{\rm sat}=-500$, and -100 MeV; see Figure 8 (a). This
implies that the appearance of intermediate-mass twin/triplet configurations
requires a nucleonic EoS which is stiff in the entire relevant range, so that
the radius of a $1.7\,M_{\odot}$ star $R_{1.7}\gtrsim 13.5$ km. However, such
nucleonic EoS all predict $R_{1.4}\gtrsim 13.2$ km which is inconsistent with
GW170817’s constraint. As a result, the few valid hybrid EoS models in our
collection are those constructed from nucleonic EoS with $Q_{\rm sat}\sim 900$
and $L_{\rm sym}\sim 45$ MeV. Although we find twin/triplet configurations, as
seen in Figure 8 (b), the differences in their radii ($\Delta R\lesssim 1$ km)
and TDs ($\Delta\Lambda\lesssim 100$) are beyond current detection capability.

We now turn to the possible constraints placed by the TDs on our collection of
hybrid EoS models. To this end, we set the chirp mass
$\mathcal{M}=1.44\,M_{\odot}$ inferred from the GW190425 event (Abbott et al.,
2020a). For this binary, the 90% CIs for the component masses range from
$1.46$ to $1.87\,{M}_{\odot}$ if we restrict a low-spin prior (Abbott et al.,
2020a). The GW190425 observational analysis, however, does not provide
significantly novel information on the stellar matter EoS (Abbott et al.,
2020a). For instance, the estimation of the combined dimensionless TD
$\tilde{\Lambda}_{1.44}\lesssim 600$, which could be converted to
$\Lambda_{1.654}\lesssim 600$ for a star with mass $M=1.654\,M_{\odot}$. This
upper limit is consistent with the values of TDs predicted by our stiffest
nucleonic EoS ($Q_{\rm sat}=900,L_{\rm sym}=105$ MeV); see Figure 1 where the
values of TDs are quoted.

Figure 9:  TDs of compact objects with single (in blue) and double phase (in
red) transitions for a fixed value of binary chirp mass
$\mathcal{M}=1.44\,M_{\odot}$ inferred from the GW190425 event (Abbott et al.,
2020a). The hybrid EoS models are constructed from a nucleonic model with
$M^{\rm N}_{\rm max}=1.7\,M_{\odot}$ of the hadronic branch, $M^{\rm Q1}_{\rm
max}=1.7\,M_{\odot}$, $s_{1}=0.7$ of the Q1 branch, and $M^{\rm Q2}_{\rm
max}=2.1\,M_{\odot}$, $s_{2}=1.0$ of the Q2 branch. For models with mass
twins, the two types of pairs for stars with masses $M_{1}$ and $M_{2}$ are
Q2-N and N-N, while for models with mass triplet the three types of pairs are
Q2-N, Q2-Q1, and N-N.

Figure 9 displays the TDs $\Lambda_{1}$ and $\Lambda_{2}$ of the stars
involved in a binary with masses $M_{1}$ and $M_{2}$ in the cases of single
and double phase transitions for a selection of hybrid EoS. The models are
constructed from stiff-soft EoS with $Q_{\rm sat}=900$ and $L_{\rm sym}=45$
MeV by setting $M^{\rm N}_{\rm max}=1.7\,M_{\odot}$, $M^{\rm Q1}_{\rm
max}=1.7\,M_{\odot}$, $s_{1}=0.7$ and $M^{\rm Q2}_{\rm max}=2.1\,M_{\odot}$,
$s_{2}=1.0$. Let us recall that in the case of a single phase transition it is
assumed that the transition takes place directly from the hadronic phase to
the Q2 quark phase. The diagonal line corresponds to the case of an equal-mass
binary with $M_{1,2}=1.654\,M_{\odot}$. In this case the maximum mass of the
nucleonic branch of the sequence $M^{\rm N}_{\rm max}=1.70\,M_{\odot}$ is
slightly higher than the value of the equal-mass case $1.654\,M_{\odot}$
(Abbott et al., 2020a). The possible types of pairs differ from previous
results in Figure 6. In general, two types of pairs of CSs could be involved
in such a merger event, namely Q2-N and N-N for EoS models with a single phase
transition featuring mass twins; and three types of pairs of CSs, namely Q2-N,
Q2-Q1, and N-N for EoS models with two sequential phase transitions featuring
triplets. As far as observations are concerned, we note that the difference in
the TDs of triplets $\Delta\Lambda_{1.7}$ for N-Q2 pairs is, at most, $\sim$
100\. This makes it challenging to distinguish between nucleonic and hybrid
stars by analyzing the TDs.

###  3.3 Hybrid equation of state featuring massive twins or triplets

Figure 10:  Same as in Figure 4, but for hybrid EoS models with an
intermediate-soft nucleonic EoS that allows single (dashed lines) or double
(solid lines) phase translations to appear in high-mass stars. The results are
constructed from a nucleonic EoS ($Q_{\rm sat}=900,L_{\rm sym}=45$ MeV) by
fixing $M^{\rm N}_{\rm max}/M_{\odot}=2.0$ of the nucleonic branch, while
varying the maximum mass $M^{\rm Q2}_{\rm max}/M_{\odot}=1.9$, 2.0, and 2.1 of
the Q2 branch.

Finally, we consider an extreme class of models where the phase transitions
occur at high densities $\rho_{\text{tran}}\sim 3.0\,\rho_{\rm sat}$, i.e., a
strong phase transition takes place in CSs of the nucleonic branch with a mass
close to the value $\sim 2.0\,M_{\odot}$ of observed massive pulsars. In this
case, NICER’s 90% CI for PSR J0740+6620 does not provide significant
constraints on the possibility of phase transitions. This is because the
maximum mass of the nucleonic sequence satisfies the requirements set by the
NICER measurement and any hybrid branches of CSs are allowed to exist even
outside the 90% CI region.

Figure 10 shows model EoS for fixed $Q_{\rm sat}=900,L_{\rm sym}=45$ MeV,
$M^{\rm N}_{\rm max}/M_{\odot}=2.0$ (which corresponds to
$\rho_{\text{tran}}=2.76\,\rho_{\rm sat}$), and $(s_{1},s_{2})=(0.7,1.0)$ for
values of $M^{\rm Q2}_{\rm max}/M_{\odot}=1.9$, 2.0, and 2.1. For the model
with $M^{\rm Q2}_{\rm max}/M_{\odot}=1.9$, i.e., a case of $M^{\rm Q2}_{\rm
max}<M^{\rm N}_{\rm max}$, we drop the first condition in Equation (8) and use
instead $M^{\rm{Q2}}_{\rm{max}}=M^{\rm{Q1}}_{\rm{max}}$. It is seen from the
insets of Figure 10 that the most pronounced mass twins and/or mass triplets
appear in the sequences when $M^{\rm Q2}_{\rm max}/M_{\odot}\simeq 2.0$.
Figure 11 gives the ranges of parameters that characterize twin and triplet
configurations as a function of the model parameters defining the EoS. In this
class, again, we find no twin/triplet solutions for models with $Q_{\rm
sat}=-500$ MeV.

From Figures 10 and 11 (a) we conclude that the range of masses containing
twin/triplet configurations (if they are allowed by the parameters of the
nucleonic EoS) is somewhat extended compared to the $M^{\rm N}_{\rm
max}/M_{\odot}=1.3$ and 1.7 cases discussed in previous subsections. The
values of $\Delta M_{\rm twin}$ could be as large as $0.1\,M_{\odot}$, which
is in principle within the accuracy of mass measurements via the relativistic
Shapiro time delay. The radius difference between the radii of twins, $\Delta
R_{\rm twin}$, can reach up to 2-3 km, which is again within the accuracy that
is achieved routinely in NICER data analysis. However, in the case of triplet
configurations, the values of $\Delta M_{\rm triplet}$ are, at most,
$0.01\,M_{\odot}$, which is beyond the current detection capability. The
radius difference of triplets, $\Delta R_{\rm triplet}$, can be as large as 1
km for Q1-N or Q1-Q2 pairs, and 2 km for N-Q2 pairs; see Figure 11 (b). The
differences in the TDs of twins and triplets $\Delta\Lambda_{2.0}$, as
expected, are only several tens; see Figure 11 (b). This implies that TD
measurements are not useful for distinguishing very massive nucleonic stars
from their hybrid twins. We recall that we consider only the class of hybrid
EoS models for which the nucleonic EoS satisfies the GW170817 constraints.

Figure 11:  Same as in Figure 5, but for hybrid EoS models with $M^{\rm
N}_{\rm max}/M_{\odot}=2.0$. Notice that no twin/triplet configurations were
found for nucleonic EoS models with $Q_{\rm sat}=-500$ MeV.

###  3.4 Remarks

Let us close this section with some general remarks. Combining the results
presented above, it can be concluded that (i) if the nucleonic EoS is stiff
(i.e., there is tension with GW170817’s inference) and multiple stable
branches with twin or triplet configurations exist, then the transition from
the nucleonic branch likely happens at a low density,
$\rho_{\text{tran}}\lesssim 2.0\,\rho_{\rm sat}$, with a corresponding mass
$M^{\rm N}_{\rm max}\lesssim 1.3$-$1.4\,M_{\odot}$; and (ii) for a soft
nucleonic EoS, which is consistent with GW170817’s inference, the transition
density could be higher, $\rho_{\text{tran}}\sim 3.0\,\rho_{\rm sat}$, with a
corresponding mass $M^{\rm N}_{\rm max}\sim 2.0\,M_{\odot}$.

We found narrow ranges of masses and a modest radius difference for low-mass
and intermediate-mass twins/triplets, which are a consequence of the
constraint imposed by the large radius of PSR J0740+6620. It limits the
allowed range of the reduction of the radius of a hybrid star and makes it
challenging to distinguish between nucleonic and hybrid stars. The situation
is more optimistic if the nucleonic branch reaches the maximum value $\sim
2.0\,M_{\odot}$. Finally, note that in all the above discussed models, a
larger mass range for twin/triplet configurations exists for models that have
a large value of sound-speed square $s_{2}\leq 1$.

##  4 Conclusions

In this paper, we extended our analysis (Li et al., 2021b) of the two recent
observational/experimental results — the inference of the radius of PSR
J0740+6620 (Riley et al., 2021; Miller et al., 2021) and the neutron skin
thickness by analysis of the PREX-II experiment (Reed et al., 2021; Reinhard
et al., 2021) — to a broader class of models. These now include (i) double
phase transitions in the quark phase and (ii) transition to quark matter at
low, intermediate, and high densities. In doing so, we continued using a
density-functional approach to nucleonic matter which is calibrated to nuclear
phenomenology, and the CSS parameterization to describe high-density quark
matter.

We have assessed the existence of twin or triplet configurations for three
classes of nucleonic stars. For the class where nucleonic EoS models predict a
large radius for a canonical-mass ($M=1.4\,M_{\odot}$) star, we show that the
tension between GW170817’s inference for the radius and models with large
values of $L_{\rm sym}$ can be mitigated (if not completely removed) by a low-
density phase transition to quark matter and the formation of low-mass hybrid
stars. This leads always to the appearance of twins and/or triplets, in the
cases of single and double phase transitions, respectively. For the class of
less stiff nucleonic EoS models that are consistent with GW170817’s inference,
we find that for nucleonic branches extending up to $\sim 1.7\,M_{\odot}$,
twins or triplets only exist in very small mass ranges, i.e., $\Delta M_{\rm
twin}\lesssim 0.04\,M_{\odot}$ and $\Delta M_{\rm triplet}\lesssim
0.002\,M_{\odot}$.

If the nucleonic branch extends up to $\sim 2.0\,M_{\odot}$ (i.e., the radius
of PSR J0740+662 can be attributed to the nucleonic branch), twins and
triplets can exist in narrow mass ranges $\Delta M_{\rm twin}\lesssim
0.1M_{\odot}$ for twins and $\Delta M_{\rm triplet}\lesssim 0.01M_{\odot}$ for
triplets. As expected, the ranges of mass and radius (and TD) for the
existence of mass twins and mass triplets are larger for the stiff nucleonic
EoS model. The largest ranges of twins and triplets are generally supported by
models with $s_{1,2}\leq 1.0$, which allow for maximum masses $M^{\rm Q2}_{\rm
max}\sim 2.1\,M_{\odot}$. The most pronounced twins and/or triplet
configurations are found to appear in sequences with a mass of about
$2.0\,M_{\odot}$.

We also extended our previous analysis of the TD (Li et al., 2020a, 2021b) to
our current models, with a focus on sequences that contain twin and triplet
configurations. The resulting TD diagram can be used in future analysis of
binary neutron star merger events in a search for signatures of QCD phase
transition(s). In particular, we demonstrated that low-mass twins and triplets
differ quantitatively by their TDs, while massive twins and triplets differ
quantitatively by their radii; see Figures 5 and 11. This highlights the
prospects of confirmation of the existence of twin and triplet stars from
future measurements similar to those used in our analysis (i.e., X-ray
measurements of radii and GW measurements of TDs). Specifically, we expect a
significantly smaller radius (by about 1-2 km) for hybrid stars compared to
purely nucleonic stars. In merger events, finding stars with similar masses
but significantly different values of TD would be a clear indication of a
phase transition. In this case, the TDs will be drawn, respectively, from the
disconnected branches in the $\Lambda$-$\Lambda$ diagrams, one corresponding
to the nucleonic branch and another to the hybrid branch. There are other
observables that we have not discussed that can indicate a first-order phase
transition. For example, Bauswein et al. (2019) find that the dominant
postmerger GW frequency for hybrid and nucleonic stars differ significantly
from each other. Furthermore, the gravitational radiation from asymmetric
supernova explosions (Fischer, 2021; Bauswein et al., 2022) will carry the
imprints of separate phase transitions. Also, matter accretion onto a CS will
lead to its compression, two successive phase transitions and two separate
explosions with energy release associated with them (Zdunik et al., 2008;
Abdikamalov et al., 2009; Lin et al., 2011). Finally, we stress once again
that the present astrophysical and nuclear physics data do not prohibit the
existence of mass twin and/or mass triplet if strong first-order phase
transitions occur in dense matter.

J. L. acknowledges the support of the National Natural Science Foundation of
China (Grant No. 12105232), the Fundamental Research Funds for the Central
Universities (Grant No. SWU-020021), and by the Venture & Innovation Support
Program for Chongqing Overseas Returnees (Grant No. CX2021007). A. S. is
supported by the Deutsche Forschungsgemeinschaft Grant No. SE 1836/5-2 and the
Polish NCN Grant No. 2020/37/B/ST9/01937 at Wrocław University. M. A. is
supported by the U.S. Department of Energy, Office of Science, Office of
Nuclear Physics under Award No. DE-FG02-05ER41375.

## References

  * Abbott et al. (2017) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017, PhRvL, 119, 161101, doi: 10.1103/PhysRevLett.119.161101
  * Abbott et al. (2018) Abbott, B. P., Abbott, R., Abbot, T. D., et al. 2018, PhRvL, 121, 161101, doi: 10.1103/PhysRevLett.121.161101
  * Abbott et al. (2019) —. 2019, PhRvX, 9, 011001, doi: 10.1103/PhysRevX.9.011001
  * Abbott et al. (2020a) —. 2020a, ApJL, 892, L3, doi: 10.3847/2041-8213/ab75f5
  * Abbott et al. (2020b) Abbott, R., Abbott, T. D., Abraham, S., et al. 2020b, ApJL, 896, L44, doi: 10.3847/2041-8213/ab960f
  * Abdikamalov et al. (2009) Abdikamalov, E. B., Dimmelmeier, H., Rezzolla, L., & Miller, J. C. 2009, MNRAS, 392, 52, doi: 10.1111/j.1365-2966.2008.14056.x10.48550
  * Adhikari et al. (2021) Adhikari, D., Albataineh, H., Androic, D., et al. 2021, PhRvL, 126, 172502, doi: 10.1103/PhysRevLett.126.172502
  * Alford et al. (2013) Alford, M. G., Han, S., & Prakash, M. 2013, PhRvD, 88, 083013, doi: 10.1103/PhysRevD.88.083013
  * Alford et al. (2001) Alford, M. G., Rajagopal, K., Reddy, S., & Wilczek, F. 2001, PhRvD, 64, 074017, doi: 10.1103/PhysRevD.64.074017
  * Alford et al. (2008) Alford, M. G., Schmitt, A., Rajagopal, K., & Schäfer, T. 2008, RvMP, 80, 1455, doi: 10.1103/RevModPhys.80.1455
  * Alford & Sedrakian (2017) Alford, M. G., & Sedrakian, A. 2017, PhRvL, 119, 161104, doi: 10.1103/PhysRevLett.119.161104
  * Alvarez-Castillo et al. (2019) Alvarez-Castillo, D. E., Blaschke, D. B., Grunfeld, A. G., & Pagura, V. P. 2019, PhRvD, 99, 063010, doi: 10.1103/PhysRevD.99.063010
  * Anglani et al. (2014) Anglani, R., Casalbuoni, R., Ciminale, M., et al. 2014, RvMP, 86, 509, doi: 10.1103/RevModPhys.86.509
  * Annala et al. (2018) Annala, E., Gorda, T., Kurkela, A., & Vuorinen, A. 2018, PhRvL, 120, 172703, doi: 10.1103/PhysRevLett.120.172703
  * Baldo & Burgio (2016) Baldo, M., & Burgio, G. F. 2016, PrPNP, 91, 203, doi: 10.1016/j.ppnp.2016.06.006
  * Bauswein et al. (2019) Bauswein, A., Bastian, N.-U. F., Blaschke, D. B., et al. 2019, PhRvL, 122, 061102, doi: 10.1103/PhysRevLett.122.061102
  * Bauswein et al. (2022) Bauswein, A., Blaschke, D., & Fischer, T. 2022, arXiv e-prints, arXiv:2203.17188, doi: 10.48550/arXiv.2203.17188
  * Biswas (2022) Biswas, B. 2022, ApJ, 926, 75, doi: 10.3847/1538-4357/ac447b
  * Blaschke et al. (2010) Blaschke, D., Klahn, T., Lastowiecki, R., & Sandin, F. 2010, JPhG, 37, 094063, doi: 10.1088/0954-3899/37/9/094063
  * Bonanno & Sedrakian (2012) Bonanno, L., & Sedrakian, A. 2012, A&A, 539, A16, doi: 10.1051/0004-6361/201117832
  * Burgio et al. (2018) Burgio, G. F., Drago, A., Pagliara, G., Schulze, H. J., & Wei, J. B. 2018, ApJ, 860, 139, doi: 10.3847/1538-4357/aac6ee
  * Christian & Schaffner-Bielich (2022) Christian, J.-E., & Schaffner-Bielich, J. 2022, ApJ, 935, 122, doi: 10.3847/1538-4357/ac75cf
  * Christian et al. (2019) Christian, J.-E., Zacchi, A., & Schaffner-Bielich, J. 2019, PhRvD, 99, 023009, doi: 10.1103/PhysRevD.99.023009
  * Collins & Perry (1975) Collins, J. C., & Perry, M. J. 1975, PhRvL, 34, 1353, doi: 10.1103/PhysRevLett.34.1353
  * Contrera et al. (2022) Contrera, G. A., Blaschke, D., Carlomagno, J. P., Grunfeld, A. G., & Liebing, S. 2022, PhRvC, 105, 045808, doi: 10.1103/PhysRevC.105.045808
  * Cromartie et al. (2020) Cromartie, H. T., Fonseca, E., Ransom, S. M., et al. 2020, NatAs, 4, 72, doi: 10.1038/s41550-019-0880-2
  * Danielewicz & Lee (2014) Danielewicz, P., & Lee, J. 2014, NuPhA, 922, 1, doi: 10.1016/j.nuclphysa.2013.11.005
  * Drischler et al. (2022) Drischler, C., Han, S., & Reddy, S. 2022, PhRvC, 105, 035808, doi: 10.1103/PhysRevC.105.035808
  * Essick et al. (2020) Essick, R., Landry, P., & Holz, D. E. 2020, PhRvD, 101, 063007, doi: 10.1103/PhysRevD.101.063007
  * Essick et al. (2021) Essick, R., Tews, I., Landry, P., & Schwenk, A. 2021, PhRvL, 127, 192701, doi: 10.1103/PhysRevLett.127.192701
  * Fattoyev et al. (2020) Fattoyev, F. J., Horowitz, C. J., Piekarewicz, J., & Reed, B. 2020, PhRvC, 102, 065805, doi: 10.1103/PhysRevC.102.065805
  * Fischer (2021) Fischer, T. 2021, EPJA, 57, 270, doi: 10.1140/epja/s10050-021-00571-z10.48550
  * Fonseca et al. (2021) Fonseca, E., Cromartie, H. T., Pennucci, T. T., et al. 2021, ApJL, 915, L12, doi: 10.3847/2041-8213/ac03b8
  * Gonçalves & Lazzari (2022) Gonçalves, V. P., & Lazzari, L. 2022, EPJC, 82, 288, doi: 10.1140/epjc/s10052-022-10273-5
  * Huth et al. (2022) Huth, S., Pang, P. T. H., Tews, I., et al. 2022, Nature, 606, 276, doi: 10.1038/s41586-022-04750-w
  * Itoh (1970) Itoh, N. 1970, PThPh, 44, 291, doi: 10.1143/PTP.44.291
  * Ivanenko & Kurdgelaidze (1965) Ivanenko, D. D., & Kurdgelaidze, D. F. 1965, Astrophysics, 1, 251, doi: 10.1007/BF01042830
  * Jokela et al. (2022) Jokela, N., Järvinen, M., & Remes, J. 2022, PhRvD, 105, 086005, doi: 10.1103/PhysRevD.105.086005
  * Khadkikar et al. (2021) Khadkikar, S., Raduta, A. R., Oertel, M., & Sedrakian, A. 2021, PhRvC, 103, 055811, doi: 10.1103/PhysRevC.103.055811
  * Klähn et al. (2013) Klähn, T., Łastowiecki, R., & Blaschke, D. 2013, PhRvD, 88, 085001, doi: 10.1103/PhysRevD.88.085001
  * Lalazissis et al. (2005) Lalazissis, G. A., Niksic, T., Vretenar, D., & Ring, P. 2005, PhRvC, 71, 024312, doi: 10.1103/PhysRevC.71.024312
  * Lattimer & Lim (2013) Lattimer, J. M., & Lim, Y. 2013, ApJ, 771, 51, doi: 10.1088/0004-637X/771/1/51
  * Legred et al. (2021) Legred, I., Chatziioannou, K., Essick, R., Han, S., & Landry, P. 2021, PhRvD, 104, 063003, doi: 10.1103/PhysRevD.104.063003
  * Li et al. (2021a) Li, A., Miao, Z., Han, S., & Zhang, B. 2021a, ApJ, 913, 27, doi: 10.3847/1538-4357/abf355
  * Li & Sedrakian (2019) Li, J. J., & Sedrakian, A. 2019, PhRvC, 100, 015809, doi: 10.1103/PhysRevC.100.015809
  * Li et al. (2020a) Li, J. J., Sedrakian, A., & Alford, M. 2020a, PhRvD, 101, 063022, doi: 10.1103/PhysRevD.101.063022
  * Li et al. (2021b) —. 2021b, PhRvD, 104, L121302, doi: 10.1103/PhysRevD.104.L121302
  * Li et al. (2020b) Li, J. J., Sedrakian, A., & Weber, F. 2020b, PhLB, 810, 135812, doi: 10.1016/j.physletb.2020.135812
  * Lin et al. (2011) Lin, W., Li, B.-A., Xu, J., Ko, C. M., & Wen, D. H. 2011, PhRvC, 83, 045802, doi: 10.1103/PhysRevC.83.04580210.48550
  * Malfatti et al. (2020) Malfatti, G., Orsaria, M. G., Ranea-Sandoval, I. F., Contrera, G. A., & Weber, F. 2020, PhRvD, 102, 063008, doi: 10.1103/PhysRevD.102.063008
  * Margueron et al. (2018) Margueron, J., Hoffmann Casali, R., & Gulminelli, F. 2018, PhRvC, 97, 025806, doi: 10.1103/PhysRevC.97.025806
  * Miao et al. (2020) Miao, Z., Li, A., Zhu, Z., & Han, S. 2020, ApJ, 904, 103, doi: 10.3847/1538-4357/abbd41
  * Miller et al. (2019) Miller, M. C., Lamb, F. K., Dittmann, A. J., et al. 2019, ApJL, 887, L24, doi: 10.3847/2041-8213/ab50c5
  * Miller et al. (2021) —. 2021, ApJL, 918, L28, doi: 10.3847/2041-8213/ac089b
  * Montana et al. (2019) Montana, G., Tolos, L., Hanauske, M., & Rezzolla, L. 2019, PhRvD, 99, 103009, doi: 10.1103/PhysRevD.99.103009
  * Most et al. (2018) Most, E. R., Weih, L. R., Rezzolla, L., & Schaffner-Bielich, J. 2018, PhRvL, 120, 261103, doi: 10.1103/PhysRevLett.120.261103
  * Oertel et al. (2017) Oertel, M., Hempel, M., Klähn, T., & Typel, S. 2017, RvMP, 89, 015007, doi: 10.1103/RevModPhys.89.015007
  * Paschalidis et al. (2018) Paschalidis, V., Yagi, K., Alvarez-Castillo, D., Blaschke, D. B., & Sedrakian, A. 2018, PhRvD, 97, 084038, doi: 10.1103/PhysRevD.97.084038
  * Pereira et al. (2018) Pereira, J. P., Flores, C. V., & Lugones, G. 2018, ApJ, 860, 12, doi: 10.3847/1538-4357/aabfbf
  * Pisarski et al. (2019) Pisarski, R. D., Skokov, V. V., & Tsvelik, A. 2019, Univ, 5, 48, doi: 10.3390/universe5020048
  * Raaijmakers et al. (2021) Raaijmakers, G., Greif, S. K., Hebeler, K., et al. 2021, ApJL, 918, L29, doi: 10.3847/2041-8213/ac089a
  * Rau & Sedrakian (2022) Rau, P. B., & Sedrakian, A. 2022, arXiv e-prints, arXiv:2212.09828.  https://arxiv.org/abs/2212.09828
  * Reed et al. (2021) Reed, B. T., Fattoyev, F. J., Horowitz, C. J., & Piekarewicz, J. 2021, PhRvL, 126, 172503, doi: 10.1103/PhysRevLett.126.172503
  * Reinhard et al. (2021) Reinhard, P.-G., Roca-Maza, X., & Nazarewicz, W. 2021, PhRvL, 127, 232501, doi: 10.1103/PhysRevLett.127.232501
  * Rezzolla et al. (2018) Rezzolla, L., Most, E. R., & Weih, L. R. 2018, ApJL, 852, L25, doi: 10.3847/2041-8213/aaa401
  * Riley et al. (2019) Riley, T. E., Watts, A. L., Bogdanov, S., et al. 2019, ApJL, 887, L21, doi: 10.3847/2041-8213/ab481c
  * Riley et al. (2021) Riley, T. E., Watts, A. L., Ray, P. S., et al. 2021, ApJL, 918, L27, doi: 10.3847/2041-8213/ac0a81
  * Rodriguez et al. (2021) Rodriguez, M. C., Ranea-Sandoval, I. F., Mariani, M., et al. 2021, JCAP, 02, 009, doi: 10.1088/1475-7516/2021/02/009
  * Sedrakian et al. (2021) Sedrakian, A., Li, J.-J., & Weber, F. 2021, arXiv e-prints, arXiv:2105.14050.  https://arxiv.org/abs/2105.14050
  * Sedrakian et al. (2020) Sedrakian, A., Weber, F., & Li, J. J. 2020, PhRvD, 102, 041301, doi: 10.1103/PhysRevD.102.041301
  * Sieniawska et al. (2019) Sieniawska, M., Turczanski, W., Bejger, M., & Zdunik, J. L. 2019, A&A, 622, A174, doi: 10.1051/0004-6361/201833969
  * Tan et al. (2022) Tan, H., Dexheimer, V., Noronha-Hostler, J., & Yunes, N. 2022, PhRvL, 128, 161101, doi: 10.1103/PhysRevLett.128.161101
  * Tan et al. (2022) Tan, H., Dore, T., Dexheimer, V., Noronha-Hostler, J., & Yunes, N. 2022, PhRvD, 105, 023018, doi: 10.1103/PhysRevD.105.023018
  * Tang et al. (2021) Tang, S.-P., Jiang, J.-L., Han, M.-Z., Fan, Y.-Z., & Wei, D.-M. 2021, PhRvD, 104, 063032, doi: 10.1103/PhysRevD.104.063032
  * Tews et al. (2018) Tews, I., Margueron, J., & Reddy, S. 2018, PhRvC, 98, 045804, doi: 10.1103/PhysRevC.98.045804
  * Zdunik et al. (2008) Zdunik, J. L., Bejger, M., Haensel, P., & Gourgoulhon, E. 2008, A&A, 479, 515, doi: 10.1051/0004-6361:20078346
  * Zdunik & Haensel (2013) Zdunik, J. L., & Haensel, P. 2013, A&A, 551, A61, doi: 10.1051/0004-6361/201220697
  * Zhang & Li (2021) Zhang, N.-B., & Li, B.-A. 2021, ApJ, 921, 111, doi: 10.3847/1538-4357/ac1e8c
  * Zhang et al. (2018) Zhang, N.-B., Li, B.-A., & Xu, J. 2018, ApJ, 859, 90, doi: 10.3847/1538-4357/aac027


Sub-TeV hadronic interaction model differences and their impact on air showers

M. Schmelling1$\star$, Á. Pastor-Gutiérrez1, H. Schorlemmer2,3,1, R.D.
Parsons4,1

1 Max-Planck-Institut für Kernphysik, Heidelberg, Germany

2 IMAPP , Radboud University Nijmegen, Nijmegen, The Netherlands

3 Nationaal Instituut voor Kernfysica en Hoge Energie Fysica (NIKHEF), Science
Park, Amsterdam, The Netherlands

4 Institut für Physik, Humboldt-Universität zu Berlin, Berlin, Germany

*<EMAIL_ADDRESS>
51st International Symposium on Multiparticle Dynamics (ISMD2022)
Pitlochry, Scottish Highlands, 1-5 August 2022 10.21468/SciPostPhysProc.?

## Abstract

In the sub-TeV regime, the most widely used hadronic interaction models
disagree significantly in their predictions of particle spectra from cosmic
ray induced air showers. We investigate the nature and impact of model
uncertainties, focussing on air shower primaries with energies around the
transition between high and low energy hadronic interaction models, where the
dissimilarities are largest and which constitute the bulk of the interactions
in air showers.

##  1 Introduction

The description of air showers created by high energy cosmic ray primaries
hitting the atmosphere requires the modelling of hadronic interactions between
elementary particles and air nuclei over many orders of magnitude in energy.
Model differences in the description of high energy interactions are discussed
in [1], the low energy region has been studied in [2], some key results of
which are summarised below.

The scenario considered is sketched in fig. 1 and modelled by the CORSIKA
v7.64 [3] air shower simulation software package. A primary proton with a
total lab energy of 100 GeV and zero zenith angle interacts with a nitrogen
nucleus at an altitude of 17550 m. The observation level for the final state
is at 4100 m. Fixing the height of the first interaction removes geometric
effects caused by varying ground distances and puts the focus on differences
between the physics modelling. Here EPOS-LHC [4], QGSJetII-04 [5], SIBYLL 2.3c
[6] and UrQMD [7] are considered. UrQMD is designed for lab energies from less
than 100 MeV up to O(200) GeV, the other models can be applied for lab
energies above O(40) GeV. A full shower simulation starts with a high-energy
(HE) model for the initial part of the shower development and switches to a
low-energy (LE) model at a transition energy, which in CORSIKA is set to 80
GeV. With a 100 GeV primary all four models are suitable for the initial
interaction, the evolution below the transition energy is modelled by UrQMD.

Figure 1: Sketch of the scenario considered in this study.

##  2 Event classification

Ground level observables are affected by the interplay of what happens in the
first, highest energy, interaction and the subsequent lower energy processes.
Here the first interaction will be characterised by the inelasticity $\kappa$
of the event and the type of the leading particle, i.e. the secondary with the
highest energy. The inelasticity is defined as

$\kappa=1-\frac{E_{\rm LP}}{E_{\rm FI}}\approx 1-x_{F}^{\rm LP}\;,$ (1)

where $E_{\rm LP}$ is the energy of the leading particle and $E_{\rm FI}$ the
total energy of all final state particles. This definition ensures that
$\kappa$ is in the range $[0,1]$ and is insensitive to small violations of
energy momentum conservation that are observed in all models. The inelasticity
is related to Feyman’s scaling variable $x_{F}$ and provides a qualitative
measure for the amount of energy that goes into the production of new
particles.

Nucleons: | $p,n$
---|---
Muonic family: | $\mu,\pi,K$
EM component: | $\gamma,e$
Others: | $\Lambda,\Sigma,\Xi,\Omega\ldots$

Figure 2: Inelasticity distributions for events with different types of
leading particles in EPOS-LHC, and definition of the particle families used to
classify events. Antiparticles and different charge states are implied.

For the leading particle we differentiate between nucleons, the muonic family,
which contains particles that either directly or via decays contribute to the
ground level muon flux, the EM component, which leads to electromagnetic
showers and others, which are less important for the ground level observables.
Figure 2 shows the inelasticity distribution in EPOS-LHC for events with
different types of leading particles. For small inelasticities the subsequent
shower evolution is driven by leading nucleons, particles from the muonic
family become important at large inelasticities. In the following we consider
three regions: $\kappa_{1}\in[0,0.2]$, the elastic and diffractive region, the
transition region $\kappa_{2}\in[0.2,0.4]$, and $\kappa_{3}\in[0.4,1]$, the
highly inelastic regime.

##  3 Ground level observables

###  3.1 Lateral distributions of muons and electromagnetic energy flow

Figure 3 shows the lateral distributions of muons and the electromagnetic
energy flow at ground level when the leading particle of the first interaction
is a nucleon. As a function of the distance to the shower axis, the muon flux
is shown in units of muons per square meter, for the electromagnetic component
the energy deposit is given in GeV per square meter.

Figure 3: Lateral distributions electromagnetic energy flow (top row) and of
muon numbers (bottom row) for different inelasticity regions and different HE
interaction models.

The electromagnetic energy flow is dominated by events from the diffractive
region since in events with small inelasticity the bulk of the particle
production happens deeper in the atmosphere, with the consequence that more of
the electromagnetically interacting particles reach the ground. For the same
reason also a sizeable fraction of the muon flux comes from the diffractive
$\kappa_{1}$ region. In addition there is a large contribution from highly
inelastic events, where many particles are created that decay into muons. The
model with the largest number of muons close to the shower centre is
QGSJetII-04. The excess is even more pronounced for events with a leading
particle from the muonic family [2].

###  3.2 Muon flux at ground level

The overall numbers of muons at ground level are more similar since the area
at small distances from the shower centre is only small fraction of the total.
The average number of muons per event as a function of the inelasticity is
shown in fig. 4, together with the $\kappa$ distributions of those events that
dominate the muon flux at ground level, i.e. those with a leading nucleon or a
leading particle from the muonic family. In the muon numbers the HE models
show a discontinuity at $\kappa\approx 0.2$. Below that value the first two
interactions, above only the first interaction is modelled by one of EPOS-LHC,
QGSJetII-04 or SIBYLL 2.3c. The spread of the curves for $\kappa>0.2$ thus
shows the importance of the first interaction for ground level observables,
the jump at $\kappa\approx 0.2$ reflects the differences of the physics
modelling between the HE models and UrQMD at the transition energy. The
difference is largest for QGSJetII-04. Smaller but still significant jumps are
observed for EPOS-LHC and SIBYLL 2.3c.

Figure 4: Inelasticity distributions (left) of the first interaction for
events where the leading particle if a nucleon or belongs to the muonic
family, and (right) average number of ground-level muons per event as a
function of the inelasticity $\kappa$. The filled bands are for events with
leading nucleons in the first interaction, the open ones are for all events.

##  4 Conclusion

Studies of hadronic interaction models in the transition region between the HE
and LE regime reveal significant differences in the prediction of ground level
observables. Those can be traced to differences in the final states of the
first interaction generated by the HE models and a discontinuity in the
physics modelling when switching from a HE model to the LE model. To improve
the models, their predictions could be compared to data from accelerator
experiments, which are available at nucleon-nucleon centre-of-mass energies
$\sqrt{s_{NN}}$ of $O(10)$ GeV from SPS-data, $O(100)$ GeV from LHCb fixed
target data and $O(10)$ TeV from beam-beam collisions at the LHC.

##  Appendix A Energy-momentum violation in air shower models

Using the CRMC interface [8], the hadronic interaction models can be run
standalone in order to compare to accelerator data or for simple checks of the
kinematics. Taking for example an incoming proton with a lab momentum $q=100$
GeV/$c$ colliding with a nucleus at rest, momentum conservation requires that
the momentum sum of all final state particles is $100$ GeV/$c$. As shown in
fig. 5 there are deviations.

In SIBYLL-2.3c the effects are small, with slight offsets that are related to
the number of target nucleons involved in the interaction. This can be seen by
boosting the event to the nucleon-nucleon centre-of-mass, where the actual
event generation happens. For an interaction with $N>1$ target nucleons and a
boost that goes to the centre-of-mass for $N=1$, the boost generates a
momentum excess of $1-N$ times the centre-of-mass momentum of the incoming
proton. For SIBYLL-2.3c one finds delta-functions that are marginally
displaced from the expected values and correspond to the small shifts seen in
the lab system.

For EPOS-LHC the centre-of-mass distributions show a sizeable smearing around
the expected values, which correspond to violations of momentum conservation
of up to $10$ GeV/$c$ in the lab system. The large deviations can be avoided
by improving the numerical precision of the model, changing in line 17 of the
function epos-uti.f the value 0.5 to e.g.0.005 [9]:

if(iLHC.eq.1) errlim=max(0.00005,0.5/engy)

The boost to the nucleon-nucleon centre-of-mass allows one also to determine
the number of interacting target nucleons. Here, in SIBYLL-2.3c only up to 5
target nucleons interact whereas in EPOS-LHC the full range from 1–14 except
for N=9 and N=13 is covered.

Figure 5: Violation of momentum conservation in SIBYLL 2.3C and EPOS-LHC. The
left hand plots show the effect in the lab system, the right hand plots in the
nucleon-nucleon centre-of-mass system. The red dashed lines indicate the
expected values when $N=1,2,\ldots$ target nucleons are involved in the
interaction.

The importance of correctly modelling sub-TeV hadronic interactions for air
showers is underlined by fig. 6, which, for 1 PeV primaries, shows the average
number of interactions in the shower as a function of the particle energy. The
number of interactions rises exponentially as the particle energy drops, such
that the bulk of the particle production occurs at energies at or below the
transition energy. The plot also shows that the number of sub-TeV interactions
in SIBYLL 2.3c and QGSJetII-04 is up to 20% lower than for EPOS-LHC.

Figure 6: Average number of interactions as a function of the particle energy
in air showers initiated by 1 PeV protons hitting the atmosphere. The
histogram (green, right hand ordinate) is the prediction by EPOS-LHC, the
points (left hand ordinate) show the ratios to EPOS-LHC for SIBYLL 2.3c and
QGSJetII-04.

## References

  * [1] R.D. Parsons et al., Systematic differences due to high energy hadronic interaction models in air shower simulations in the 100 GeV–100 TeV range, Phys. Rev. D 100(2), 023010 (2019), https:// doi.org/10.1103/PhysRevD.100.023010.
  * [2] Á. Pastor-Gutiérrez et al., Sub-TeV hadronic interaction model differences and their impact on air showers, Eur.Phys.J. C 81 369 (2021),
https://doi.org/10.1140/epjc/s10052-021-09160-2.

  * [3] D. Heck et al.,CORSIKA : A Monte Carlo Code to Simulate Extensive Air Showers, Forschungszentrum Karlsruhe, Report FZK 6019 (1998),
https://web.iap.kit.edu/corsika/physics_description/corsika_phys.pdf.

  * [4] T. Pierog et al., EPOS-LHC: test of collective hadronization with data measured at the CERN large hadron collider Phys. Rev. C 92, 034906 (2015),
https://doi.org/10.1103/PhysRevC.92.034906.

  * [5] S. Ostapchenko, Monte Carlo treatment of hadronic interactions in enhanced Pomeron scheme: QGSJET-II model, Phys. Rev. D 83, 014018 (2011),
https://doi.org/10.1103/PhysRevD.83.014018.

  * [6] F. Riehn et al., The hadronic interaction model Sibyll 2.3c and Feynman scaling, 35th International Cosmic Ray Conference, PoS(ICRC2017), 301 (2017),
https://doi.org/10.22323/1.301.0301.

  * [7] S. Bass et al., Microscopic models for ultrarelativistic heavy ion collisions, Prog. Part. Nucl. Phys. 41 255 (1998), https://doi.org/10.1016/S0146-6410(98)00058-1.
  * [8] R. Ulrich et al., Cosmic Ray Monte Carlo Package, CRMC,
https://doi.org/10.5281/zenodo.4558706

  * [9] T. Pierog, private communication.