# QHC21 equation of state of neutron star matter – in light of 2021 NICER data

Toru Kojo  Key Laboratory of Quark and Lepton Physics (MOE) and Institute of
Particle Physics, Central China Normal University, Wuhan 430079, China  Gordon
Baym  Illinois Center for Advanced Studies of the Universe  Department of
Physics, University of Illinois at Urbana-Champaign, 1110 W. Green Street,
Urbana, IL 61801, USA  RIKEN iTHEMS Program, Wako, Saitama 351-0198, Japan
Tetsuo Hatsuda  RIKEN iTHEMS Program, Wako, Saitama 351-0198, Japan

###### Abstract

The recent NICER measurement of the radius of the neutron star PSR J0740+6620,
and the inferred small variation in neutron star radii from $1.4M_{\odot}$ to
$2.1M_{\odot}$, suggest that the neutron star equation of state remains
relatively stiff up to baryon densities $n\sim$ 2-4 times nuclear saturation
density, $n_{0}$ – the region where we expect hadronic matter to be undergoing
transformation into quark matter. To delineate the physics from the nuclear to
the quark matter regimes we use the quark-hadron-crossover (QHC) template to
construct an updated equation of state, QHC21. We include nuclear matter
results primarily based on chiral effective field theory, but also note
results of using nuclear matter variational calculations based on empirical
nuclear forces, thus covering the range of uncertainties in the nuclear
equation of state. To allow for a possible early transition to quark degrees
of freedom we begin the crossover regime from nucleons to quarks at
$1.5n_{0}$. The resulting equations of state are stiffer than our earlier
QHC19 at $\lesssim 2n_{0}$, predicting larger radii in substantial agreement
with the NICER data, with accompanying peaks in sound velocity at 2-4$n_{0}$.
We discuss possible microscopic mechanisms underlying stiffening of the
equation of state.

###### Subject headings:

dense matter, equation of state, stars: neutron, quarks

††slugcomment: Submitted to ApJ

##  1\. Introduction

Recent data from the NICER observatory has yielded the important inference
that the equatorial radius of the most massive neutron star known, PSR
J0740+6620, with mass $M/M_{\odot}=2.08\pm 0.07$ (Fonseca et al., 2021), is
13.7${}^{+2.6}_{-1.5}$ km, as analyzed by Miller et al. (2021) combining NICER
and XMM data for PSR J0740+6620. With nuclear physics constraints at low
density, and gravitational radiation data from GW170817 added in, their
inferred radius drops to $12.35\pm 0.75$ km for a 2.08 $M_{\odot}$ neutron
star. The analysis of Riley et al. (2021) using the NICER and XMM data for PSR
J0740+6620 gives $12.39^{+1.30}_{-0.98}$ km. These analyses lead to an updated
radius of 1.4$M_{\odot}$ neutron stars,111Here $R_{M/M_{\odot}}$ denotes the
radius of a star of mass $M$. $R_{1.4}=12.45\pm 0.65$ km (Miller et al., 2021)
and $12.33^{+0.76}_{-0.81}$ km (Riley et al., 2021) to $12.18^{+0.56}_{-0.79}$
km (Raaijmakers et al., 2021).

Our aim in this paper is to understand the implications of the new NICER data
for the equation of state of neutron star matter within the Quark-Hadron
Crossover (QHC) framework (Baym et al., 2018, 2019), which includes nucleonic
degrees of freedom at lower densities, and allow for their transitioning to
quark degrees of freedom at higher densities. QHC provides a useful baseline
to describe the available neutron star data, in which variations such as weak
first order phase transitions can be treated as perturbations. This approach
lets one relate neutron star observables to nuclear physics as well as to
quark physics, such as hadron spectroscopy and the quantum chromodynamics
(QCD) phase diagram.222QHC employs the picture of quark-hadron continuity – a
smooth crossover from baryonic to quark degrees of freedom with increasing
density, suggested in the context of the color superconductivity by Schäfer &
Wilczek (1999); Hatsuda et al. (2006) and of quantum percolation by Fukushima
et al. (2020). Neutron star phenomenology in the QHC framework was originally
introduced in (Masuda et al., 2013a, b) to study the stiffening of the
equation of state at high density in the presence of the strongly interacting
quark matter. We develop and show here the results of our new equation of
state, QHC21, which well accommodates the NICER data.

The inferred NICER radii provide clues to the properties of dense matter at
baryon densities beyond nuclear saturation density $n_{0}$, and indicate the
extent to which the central density of PSR J0740+6620 is high enough to
contain quark matter. The relatively large radii, compared with earlier
equations of state based on variational calculations of nuclear matter
starting from empirical nuclear forces (Akmal et al., 1998; Togashi et al.,
2017; Baym et al., 2019), imply that neutron star matter is relatively stiff
up to $\sim$ 2-4 $n_{0}$, a result confirmed in the equations of state of
Miller et al. (2021) used to fit the data. The radii of $1.4M_{\odot}$ neutron
stars, for which the central baryon density is $\sim$ 2-3$n_{0}$, contain
information on nuclear descriptions for 1-3$n_{0}$. Furthermore, the lack of
significant variation of the radii from $1.4M_{\odot}$ to $\sim 2.1M_{\odot}$
illuminates the properties of dense matter at densities $\gtrsim 3n_{0}$, and
rules out substantial softening of matter (as would result from a strong first
order phase transition) between 2-3$n_{0}$ and 4-5$n_{0}$ (Drischler et al.,
2021a); the radical shrinkage in the radius such softening would induce is
incompatible with the new NICER data.

The radii inferred from NICER are about 0.7-0.8 km larger than those predicted
by our earlier QHC19 equation of state (Baym et al., 2019). There, for
example, $R_{1.4}\simeq 11.6$ km, a value essentially determined by the
nuclear equation of state at densities in the range of 2$n_{0}$. Similarly,
$2.1M_{\odot}$ neutron stars were predicted in QHC19 to have radii $\simeq$
11.3-11.5 km. The challenge is how QHC can allow stiffer equations of state in
order to yield such larger radii. The ways in which matter, up to $\sim$ 2-4
$n_{0}$, can be stiffer than in the QHC19 description include stiffening in
the nucleonic regime and possible early onset of the transition to quark
degrees of freedom. Indeed, the NICER data provide the first observational
test of the relatively soft equation of states for neutron star matter
obtained by microscopic variational calculations (Akmal et al., 1998; Togashi
et al., 2017).

The NICER analysis of Miller et al. (2021) assumes the Togashi equation of
state (Togashi et al., 2017) for the crust up to baryon densities $n\lesssim
0.5n_{0}$, and parameterized equations of state – a polytrope model, a model
with a parameterized speed of sound, and a Gaussian model – beyond that
density. A more microscopically motivated approach is to use nuclear equations
of state at all densities below $\sim 1.5n_{0}$, where they are reasonably
constrained by controlled uncertainty estimates. The chiral effective field
theory ($\chi$EFT) formalism (Lonardoni et al., 2020; Drischler et al., 2021a,
b), an expansion in the effective momentum transfers, provides such a control.
As discussed by Drischler et al. (2021a) (see their Fig. 1), the uncertainty
in the pressure within $\chi$EFT is small at nuclear matter saturation
density, and grows to $\sim$ 25% at $\sim 1.4n_{0}$; the uncertainties at
2$n_{0}$ are too large to continue to use this approach at such density. The
analyses of Riley et al. (2021); Raaijmakers et al. (2021) use $\chi$EFT and
polytrope or constant speed of sound models.333 Recent analyses similar to
NICER but with nuclear constraints from $\chi$EFT (Pang et al., 2021) yield
somewhat lower estimates for the radii $R_{1.4}=11.94^{+0.76}_{-0.87}$ km and
$R_{2.08}=11.96^{+0.80}_{-0.75}$ km than the NICER analyses. Previously the
same group (Dietrich et al., 2020) gave $R_{1.4}=11.75^{+0.86}_{-0.81}$ km
without the data of PSR J0740+6620; the new NICER data leads to an increase in
the estimate by $\simeq 0.2$ km, illustrating how data for $2M_{\odot}$
neutron stars constrains neutron stars of lower mass and the equation of state
in the core.

QHC19 used the Togashi nuclear equation of state at $n\leq 2n_{0}$, quark
matter equations of state calculated within the Nambu-Jona-Lasinio (NJL) model
for $n\gtrsim 5n_{0}$, and a smooth interpolation between these equations of
state. The Togashi equation of state, is based on precision two-body nuclear
forces and empirical three-body nuclear forces, computed using variational
techniques for nuclear many-body problems, as in the Akmal-Pandharipande-
Ravenhall (APR) equation of state (Akmal et al., 1998); the same nuclear
forces are used consistently in the nuclear liquid and crust equations of
state.

Figure 1.—  Nuclear equations of state for neutron star matter in beta
equilibrium based on $\chi$EFT (collected in Raaijmakers et al. (2021)) and
variational methods based on phenomenological nuclear forces (Akmal et al.,
1998; Togashi et al., 2017). The “central” $\chi$EFT curve corresponds to the
average of the data in Drischler et al. (2021a).

Figure 1 shows the range of the pressure in nuclear equations of state vs.
baryon density in $\chi$EFT approaches (labeled Lynn, Hebeler, Tews and
Drischler) and variational approaches (labeled APR and Togashi). To allow for
stiffer matter, we focus on the nuclear equation of states constructed by
Drischler et al. (2021a) based on beyond-leading-order (N3LO) $\chi$EFT, and
the Togashi equation of state, thus covering the uncertainty range of existing
nuclear calculations at $\lesssim 1.5n_{0}$. Specifically we take the central
value of the N3LO $\chi$EFT pressure at given density. This pressure is larger
than Togashi by $\sim$ 30-40%, leading to larger radii.

We explore a greater range of physics above nuclear matter density, permitting
possibly stiffer interpolated equations of state, by using explicit nuclear
degrees of freedom only up to $n\sim 1.5n_{0}$. At the same time we also
explore quark equations of state at lower densities, $n\sim$ 3-4$n_{0}$, than
in QHC19 ($\gtrsim 5n_{0}$). In the crossover regime between the nucleonic and
quark regimes we carry out a smooth interpolation. By lowering the density for
onset of quark matter can be stiffer at lower densities in the interpolation
region.

Although theoretical calculations in the crossover region are difficult,
possibly involving baryon-like and quark-like degrees of freedom at the same
time (McLerran & Reddy, 2019; Kojo, 2021a), the physics in this regime is
highly constrained by having to match the nuclear and quark matter equations
of state at the boundaries of the crossover region. The constraints that the
sound velocity not exceed the speed of light, and that the matter in the
interpolation region be thermodynamically stable further constrain possible
interpolations. The interpolating curves reveal the general trends (Masuda et
al., 2013a, b, 2016; Kojo et al., 2015; Annala et al., 2020).

##  2\. Physics of QHC21

We detail here the construction of QHC21. We specify the thermodynamics in
terms of the pressure $P$ as a function of baryon chemical potential
$\mu_{{}_{B}}$. Nuclear equations of state are used up to a density $n^{L}$,
and quark equations of state above a density $n^{U}$ ($>n^{L}$). In the
crossover between these regimes, we smoothly interpolate the pressure, the
baryon density $n=\partial P/\partial\mu_{{}_{B}}$, and the susceptibility
$\chi_{{}_{B}}=\partial^{2}P/\partial\mu_{{}_{B}}^{2}$ at the matching
boundaries. Thermodynamic stability imposes the requirement that the
interpolated pressure leads to non-negative $\chi_{{}_{B}}$, and causality
that the adiabatic speed of sound $c_{s}=\sqrt{\partial
P/\partial\varepsilon\,}$ does not exceed the speed of light, $c$. Here
$\varepsilon=-P+\mu_{{}_{B}}n$ is the energy density.

In QHC21 we consider primarily the $\chi$EFT equation of state in the
nucleonic regime, but also note results using the Togashi equation of state
there (labeled QHC21T). To describe quark matter we employ the Nambu–Jona-
Lasinio (NJL) model. Although the NJL model does not include confining effects
it is a good framework in which to capture dynamical effects, e.g., chiral
symmetry breaking, concisely (Hatsuda & Kunihiro, 1994), and in addition it
reproduces quantities that are not very sensitive to confinement such as
effective quark masses and low energy constants in the chiral effective
Lagrangian for pions. QHC thus incorporates general features of quark matter
by exploring a wide range of the NJL model parameters consistent with neutron
star observables. The coupling constant $G$ and the ultraviolet cutoff
$\Lambda$ of the NJL model, specified in Baym et al. (2018), are the relevant
parameters describing chiral symmetry breaking. In addition, we include the
short range density-density repulsion in QCD (Kunihiro, 1991; Song et al.,
2019), of strength measured by a coupling constant $g_{{}_{\rm V}}$, and an
attraction between quarks, proportional to a coupling $H$, which governs the
color superconductivity pairing gap.444For simplicity we assume that the
effective coupling constants $g_{{}_{\rm V}}$ and $H$ are independent of
density. However, one expects from QCD that $g_{{}_{\rm V}}$ falls with
increasing density (Song et al., 2019); see Appendix A. The quark pairing
interaction, $H$, should behave similarly.

The quark equations of state, parametrized by given $(g_{{}_{\rm V}},H)$, and
the nuclear equation of state are first interpolated in a given range
$(n^{L},n^{U})$ by using fifth order interpolating polynomials in
$\mu_{{}_{B}}$, for which the coefficients are uniquely fixed by six boundary
conditions. We then determine the ranges of $(g_{{}_{\rm V}},H)$ that are
consistent with stability, causality, and the existence of the known high mass
neutron stars. This procedure predicts $M$ vs. $R$, thermodynamic quantities
such as baryon density as a function of $\mu_{{}_{B}}$, as well as more
microscopic quantities such as effective quark masses and diquark paring gaps.

Table 1 The interpolation ranges $(n^{L},n^{U})$ in units of $n_{0}$, and the interaction parameter sets $(g_{{}_{\rm V}},H)$ for QHC21 Aχ-Dχ, and QHC21T AT-DT. Tabulated for each parameter set are the maximum mass, $M_{\rm max}$, in units of $M_{\odot}$, radii $R$, in km, and central baryon densities, $n^{c}$, in units of $n_{0}$, of neutron stars of masses 1.4, 2.08 $M_{\odot}$, and $M_{\rm max}$.  | QHC21 | QHC21T
---|---|---
| Aχ | Bχ | Cχ | Dχ | AT | BT | CT | DT
$(n^{L},n^{U})$ | (1.5, 3.5) $n_{0}$ | (1.5, 3.5) $n_{0}$
$g_{{}_{\rm V}}/G$ | 1.0 | 1.1 | 1.2 | 1.3 | 0.9 | 1.0 | 1.1 | 1.2
$H/G$ | 1.50 | 1.52 | 1.54 | 1.56 | 1.50 | 1.52 | 1.55 | 1.57
$M_{\rm max}$ | 2.19 | 2.25 | 2.31 | 2.37 | 2.13 | 2.20 | 2.26 | 2.32
$R_{1.4}$ | 12.4 | 12.4 | 12.4 | 12.3 | 11.8 | 11.9 | 11.7 | 11.8
$R_{2.08}$ | 12.0 | 12.2 | 12.4 | 12.5 | 11.5 | 11.8 | 11.9 | 12.0
$R_{\rm M_{max}}$ | 11.7 | 11.5 | 11.4 | 11.2 | 10.9 | 11.1 | 11.1 | 11.3
$n^{c}|_{1.4}$ | 2.6 | 2.5 | 2.5 | 2.5 | 2.8 | 2.7 | 2.8 | 2.8
$n^{c}|_{2.08}$ | 4.3 | 4.0 | 3.5 | 3.3 | 4.9 | 4.2 | 4.0 | 3.6
$n^{c}|_{M_{\rm max}}$ | 6.0 | 5.8 | 5.6 | 5.4 | 6.4 | 6.1 | 6.0 | 5.8

##  3\.  Neutron Star Structure

Figures 2-4 display, for $(n^{L},n^{U})/n_{0}=(1.5,3.5)$, the results of QHC21
in the $(g_{{}_{\rm V}},H)$ plane for the maximum mass and the corresponding
central density, as well as the radius and central density of 2.08 $M_{\odot}$
and 1.4 $M_{\odot}$ stars. We first determine the domain of $(g_{{}_{\rm
V}},H)$ that leads to equations of state satisfying the causality condition
and that are consistent with the lower bound for the maximum mass, $M_{\rm
max}/M_{\odot}\geq 2.08\pm 0.07$; these allowed regions are shown as colored
in Figs. 2-4. As a guide we choose several characteristic values in the
$(g_{{}_{\rm V}},H)$ plane, focussing on $(g_{{}_{\rm V}},H)/G=(1.0,1.50)$
(set Aχ), $(1.1,1.52)$ (set Bχ), $(1.2,1.54)$ (set Cχ), and $(1.3,1.56)$ (set
Dχ). We denote these parameter sets with a subscript $\chi$ to distinguish
them from the sets A, B, C, and D in QHC19. In addition Table 1 summarizes the
parameters and results of both QHC21 and QHC21T (with parameter sets denoted
by a subscript T). In Appendix B we vary the range of interpolation, and find
that the domains of $(g_{{}_{\rm V}},H)$ can vary by 10-20%, with no
qualitative changes.

Figure 2.—  The maximum mass $M_{\rm max}$ (upper panel) and the central
density $n^{c}$ (lower panel) with the QHC21 equation of state. The matching
densities are $(n^{L},n^{U})/n_{0}=(1.5,3.5)$.

The upper panel of Fig. 2 shows the maximum neutron star mass $M_{\rm max}$,
in the color scale on the right, as a function of ($g_{{}_{\rm V}},H$)
measured in units of the NJL coupling $G$. Smaller repulsion $g_{{}_{\rm V}}$
leads to smaller maximum masses; the requirement $M_{\rm max}\geq
2.08M_{\odot}$ is satisfied only for $g_{{}_{\rm V}}/G\gtrsim 0.84$, while
$g_{{}_{\rm V}}/G\lesssim 0.74$ leads to $M_{\rm max}/M_{\odot}\lesssim 2.01$,
The largest allowed555This mass is in practice a lower bound on the absolute
maximum, since we use only a minimal interpolation between the nuclear and
quark regimes. As noted in Appendix B, narrowing the interpolation window can
lead to larger maximum mass, owing to small differences of the resultant
$P(\mu_{{}_{B}})$ from that with a simple interpolation.$M_{\rm
max}/M_{\odot}$ is $\simeq 2.43$ with $(g_{{}_{\rm V}},H)\simeq(1.42,1.59)$.

The lower panel of Fig. 2 shows the central density, $n^{c}$, at $M_{\rm
max}$. A larger $M_{\rm max}$ allows a smaller $n^{c}$, consistent with the
general rule that the stiffer the equation of state the larger the maximum
mass and the smaller the central density at the maximum. The restriction
$M_{\rm max}/M_{\odot}\geq 2.08$ imposes the condition $n^{c}/n_{0}\lesssim
6.4$. At $M_{\rm max}/M_{\odot}\simeq 2.40$, $n^{c}/n_{0}\simeq 5.3$.

Figure 3.—  The radius (upper panel) and the central density (lower panel) at
$2.08M_{\odot}$, in the ($g_{{}_{\rm V}},H$) plane with the QHC21 equation of
state.

Figure 3 shows the dependence of the radius and the central density of
$2.08M_{\odot}$ neutron stars on $g_{{}_{\rm V}}$ and $H$. Within the allowed
band of $(g_{{}_{\rm V}},H)$, the radii range from $\simeq 10.8$ km to $\simeq
12.9$ km. Radii $\lesssim 11$ km, significantly smaller than $R_{1.4}\simeq
12$ km, appear for equations of state whose maximum mass barely reaches
$2.08M_{\odot}$, since the radius in general shrinks radically near the
maximum mass.

Apropos of whether the central density of PSR J0740+6620 with a mass
$2.08M_{\odot}$ is high enough to accommodate quark matter, we see in the
lower panel of Fig. 3 that the central density is $\lesssim 5n_{0}$, slightly
smaller than the density where baryons of radius $\sim 0.8$ fm begin to
overlap, but considerably larger than the density where pure hadronic
calculations are applicable.

The central density, and thus the presence of quark matter, is sensitive to
how close the mass is to the maximum, $M_{\rm max}$, allowed by the equation
of state, since near $M_{\rm max}$ the radius shrinks rapidly, by 0.5-1km (see
upper panel of Fig. 5). The radius of PSR J0740+6620 being comparable to
$R_{1.4}$ indicates that it mass is not near the maximum; we expect that
$M_{\rm max}$ is larger than $2.08M_{\odot}$ by at least 0.5$M_{\odot}$. For
$(g_{{}_{\rm V}},H)$ satisfying this condition, $n^{c}\lesssim 5n_{0}$.

Figure 4.—  The radius (upper panel) and the central density (lower panel) at
$1.4M_{\odot}$ for the same equations of state as in Figs. 2 and 3.

The range of radii of $1.4M_{\odot}$ neutron stars, $\simeq 11.8$-$12.7$ km,
is consistent with the NICER result $R_{1.4}=12.4\pm 0.6$ km. The
corresponding central density is 2.1-3.2$n_{0}$. These radii are $\sim
0.3$-$1.2$ km larger than in QHC19 ($R_{1.4}\simeq 11.5$ km, with
$(n^{L},n^{U})/n_{0}=(2.0,5.0)$), a consequence primarily of using the stiffer
$\chi$EFT nuclear equation of state. The Togashi equation of state, even if
one allows more freedom for the interpolated pressure in the density range
1.5-2.0 $n_{0}$, leads to radii $\simeq$ 11.6-12.2 km, on the smaller side of
the NICER band.

Having described the overall trend of the masses, radii, and central densities
in the parameter space, we look more closely at results for the parameter sets
Aχ-Dχ for QHC21. The $M$-$R$ relation for these parameters is shown in the
upper panel of Fig. 5; the lower panel shows the corresponding results for
QHC21T for the parameter sets AT-DT. The numbers next to the heavy dots
indicate the central density of the star, in units of $n_{0}$, for the given
mass. As we have matched interpolating functions with nuclear equations of
state up to second derivatives at $n=1.5n_{0}$, the deviation of the $M$-$R$
curves from what the pure, i.e., without quark matter, $\chi$EFT or Togashi
equations of state would give appears only after $n$ increases considerably
beyond the matching point $n=1.5n_{0}$. (This trend can also be seen in the
speed of sound in Fig. 7.) The central $\chi$EFT equation of state in Fig. 1,
which is available only for $n\lesssim 2.2n_{0}$, predicts overall larger
radii for low core densities, leading to better agreement of QHC21 with the
NICER data. As we see in Fig. 6, in QHC21 the equations of state become
stiffer than QHC21T and pure Togashi at $\sim 2n_{0}$, producing
$1.4M_{\odot}$ neutron stars at a lower central density than with Togashi
alone.

An important feature of the QHC21 equation of state is the behavior of the
sound velocity. As seen in Fig. 7, $c_{s}^{2}$ exceeds the conformal limit
$c_{s}^{2}=c^{2}/3$ at $n\simeq$ 2.1-2.3$n_{0}$, a considerably lower density
than $\simeq 2.8n_{0}$ with Togashi alone, and also less than the $\sim
2.6n_{0}$ in QHC19. The peaks in $c_{s}^{2}$ for the QHC21 set Aχ-Dχ or AT-DT
are at $\simeq 2.1$-$3.0n_{0}$, a considerably density than in QHC19 they are
in the range $\simeq 3.0$-$4.0n_{0}$.

The peak in the sound velocity is a novel feature of dense QCD matter at
finite baryon density, as originally pointed out phenomenologically by Masuda
et al. (2013b) and recently analyzed theoretically by Kojo & Suenaga (2021).
It does not exist in a hot QCD plasma (Bazavov et al., 2014), where the speed
of sound develops a dip instead of a peak in the crossover region. The finite
temperature transition from a hadron resonance gas to a quark gluon plasma is
largely driven by non-relativistic resonances which are energetically
disfavored but are important due to their large entropies. Such entropic
effects are absent in cold dense matter.

Finally, we compare the pressure in QHC21 with constraints deduced by the
NICER analyses. Figure 8 shows the pressure in QHC21 as a function of the
baryon mass density, $\rho=m_{{}_{N}}n$, where $m_{{}_{N}}=1.67\times
10^{-24}$ g is the nucleon mass. Then Fig. 9 compares the pressure as a
function of baryon chemical potential $\mu_{{}_{B}}$, with the constraint “All
Measurements (Gaussian Process)” in Miller et al. (2021) at a 95%CL.666 The
error bands in Fig. 8 are primarily at larger $P$ at given $n$ than the QHC21
curves, while in Fig. 9 they occur for smaller $P$ at given $\mu_{{}_{B}}$.
This reversal is a consequence of thermodynamics. For example, if the energy
density, $\varepsilon$, is scaled by a factor $1+\eta$ at fixed $n$, then $P$
and $\mu_{{}_{B}}$ are similarly scaled. But at fixed $\mu_{{}_{B}}$, the
scaling of $P$ has the opposite sign: $\delta
P|_{\mu_{{}_{B}}}=-\eta\varepsilon$. The QHC21 equations of state are quite
consistent with the NICER inferences.

Figure 5.—  Upper panel: Mass vs. radius for QHC21 (Aχ-D${}_{\chi})$ and the
“central” $\chi$EFT equations of state Drischler et al. (2021a). Lower panel:
Mass vs. radius for QHC21T (AT-DT) and pure Togashi equations of state. In
both figures, the heavy dots indicate central densities
$n^{c}=1.3,1.5,1.7,1.9,2.1$ and 3.5 $n_{0}$. The bands (95% CL) correspond to
“PSR+GW+J0030+J0740” in Legred et al. (2021) (yellow) and “All Measurements
(Gaussian Process)” in Miller et al. (2021) (green).  Figure 6.—  Mass vs.
central density for QHC21 Aχ-Dχ (bold lines), QHC21T AT-DT (thin lines), and
the pure Togashi equations of state.  Figure 7.—  The square of the sound
speed, $c_{s}^{2}$, in units of $c^{2}$, as function of baryon density $n$ for
the QHC21 Aχ-Dχ (bold lines), QHC21T AT-DT (thin lines), and pure Togashi
equations of state. The conformal limit $c_{s}^{2}=c^{2}/3$ is shown as a
guide. The band (95%CL) shows the constraint from “PSR+GW+J0030+J0740” deduced
by Legred et al. (2021)).  Figure 8.—  Pressure in QHC21(Aχ,T, Dχ,T) compared
with the limiting constraints deduced by the NICER analyses, as a function of
baryon rest mass density $\rho=m_{N}n$. The vertical lines, $\rho$ = 0.5, 1.0,
2.0, and 6.0 $m_{N}n_{0}$, are shown as guides. The bands correspond to
“PSR+GW+J0030+J0740” in Legred et al. (2021) (blue) and “All Measurements
(Gaussian Process)” in Miller et al. (2021) (yellow). Both bands are 95% CL.
The band appears narrow in the region $\rho\sim 10^{15}$ gm/cm3 only as a
consequence of plot being logarithmic.  Figure 9.—  Pressure as a function of
baryon chemical potential $\mu_{{}_{B}}$. The heavy dots on the QHC21 curves
indicate densities $n=1,2$ and $6n_{0}$. The band (95% CL) shows the equation
of state bound “All Measurements (Gaussian Process)” in Miller et al. (2021) .

##  4\. Discussion

The NICER analyses strengthen the evidence for a peak in the sound velocity,
reflecting the stiffening of the equation of state beyond the nucleonic
regime. In the quark regime, matter can be stiffened by several physical
effects (as discussed in greater detail in Appendix A). First, the vector
repulsion between quarks, which allows quark matter to support neutron stars
above $2M_{\odot}$, stiffens the matter in general, and produces sound
velocities above the conformal limit, $c/\sqrt{3}$. Furthermore, quark pairing
correlations leading to color superconductivity (Alford et al., 2005)777Recent
work suggests that an early onset of quark matter (Annala et al., 2020) or
quark Pauli blocking effects in the transition regime (McLerran & Reddy, 2019;
Kojo, 2021a) can stiffen matter at densities lower than predicted by typical
hadronic models (see Fig. 2 in Kojo (2021b) for possible scenarios), and that
this stiffening at low density is consistent with observational constraints
(Tews et al., 2018).  or baryonic correlations in quarkyonic matter (McLerran
& Reddy, 2019) can also stiffen the equation of state.888Underlying diquark
pairing is the color magnetic interaction in QCD, which plays a crucial role
as well in hadron spectroscopy (De Rujula et al., 1975) to baryon-baryon
interactions at short distances (Oka & Yazaki, 1980; Park et al., 2020).

Within QHC21, the central density of PSR J0740+6620, $\lesssim 5n_{0}$, is
well above densities where pure hadronic calculations are valid, and is
entering the transition regime to strongly interacting quark matter. The cores
of higher mass neutron stars, possibly to be discovered in the future, could
reach the transition regime and beyond. Whether fully developed quark matter
exists in the cores will have to await fully microscopic calculations of
matter undergoing the transition from nucleonic to quark degrees of freedom.
Nonetheless the properties of quark matter strongly influence physics in the
transition, possibly even slightly above saturation density.

Tables of the equations of state QHC21 for parameter sets Aχ-Dχ and AT-DT are
available on the CompOSE archive at https://compose.obspm.fr.

## Acknowledgments

G.B. would like to thank Cole Miller and Fred Lamb for informative discussions
of the recent NICER results, including at the Aspen Center for Physics
(supported by supported by National Science Foundation Grant No. PHY-1607611).
We thank Christian Drischler, Cole Miller, Geert Raajimakers, and Isaac Legred
for kindly providing us with their data and instructions on using it. T.H. is
supported by JSPS Grant No. JP18H05236 and JST CREST Grant No. JPMJCR19T2, and
T.K. by NSFC Grant No. 11875144.

##  Appendix A A. The conformal limit for the sound velocity

In QCD at ultrahigh densities in the absence of significant length scales, the
velocity approaches the conformal limit, $c_{s}=c/\sqrt{3}$, from below; this
limiting value is the conformal bound. We focus in this Appendix on two
questions: The first is how the equation of state determines whether the sound
velocity can exceed the conformal bound.999Tews et al. (2018), also Greif et
al. (2019), assuming the validity of nuclear calculations up to $n\sim
1.5n_{0}$ where $c_{s}^{2}$ is only $\sim$ 0.1 - 0.2, point out that to have
sufficient stiffening at higher density in order to reach neutron stars over
two solar masses, it is necessary that $c_{s}^{2}$ exceed 1/3 there. On the
other hand, Annala et al. (2020) construct stiff equations of state by
allowing $c_{s}^{2}$ to reach $\sim 1/3$ from below at $n=1.1$-$1.5n_{0}$, a
velocity well above the nuclear results; such early stiffening results in
stiff equations of state that do not violate the conformal limit. In
considering quark matter in neutron stars we are particularly interested in
the domain where quark descriptions are natural but the matter is not dense
enough for perturbative treatments (valid at $\gtrsim 40n_{0}$). While such a
domain, above 5$n_{0}$, may not show up in the cores of neutron stars,
understanding $c_{s}$ there should give important constraints on the behavior
of equations of state in the range of baryon density $\sim 2$-$5n_{0}$ which
is only weakly constrained by perturbative results alone. The second question
we discuss is the behavior of the sound velocity in QED, a theory that is not
asymptotic free. As we show, in the perturbative treatment of the massless
electron gas, the sound velocity exceeds $c/\sqrt{3}$. We set $c=1$ here for
simplicity.

The sound velocity involves a microscopic interplay between the kinetic energy
density, $\varepsilon_{\rm kin}$ and the interaction energy, $\varepsilon_{\rm
int}$. We first show schematically how the conformal bound can be violated in
relativistic matter. The quark number density, $n_{{}_{\rm Q}}$, is given in
terms of the quark Fermi momentum, $p_{{}_{\rm F}}$, for equal mass quarks, by
$n_{{}_{\rm Q}}=N_{f}p_{{}_{\rm F}}^{3}/\pi^{2}$, with $N_{f}$ the number of
quark flavors present; the kinetic energy density for massless quarks, is thus
$\sim p_{{}_{\rm F}}^{4}\sim n_{{}_{\rm Q}}^{4/3}$. Calculationally it is
simplest to write the total energy density $\varepsilon(n_{{}_{\rm Q}})$ in
terms of $p_{{}_{\rm F}}$ as

$\displaystyle\varepsilon=ap_{{}_{\rm F}}^{4}\left(1+h(p_{{}_{\rm
F}})\right),$ (A1)

where $a=3N_{f}/4\pi^{2}$ and $h(p_{{}_{\rm F}})=\varepsilon_{\rm
int}/\varepsilon_{\rm kin}$ is the interaction energy relative to the kinetic
energy. Then the pressure is

$\displaystyle P=\frac{a}{3}p_{{}_{\rm F}}^{4}(1+h+p_{{}_{\rm F}}h^{\prime}),$
(A2)

where primes denote derivatives with respect to $p_{{}_{\rm F}}$. Equations.
(A1) and (A2) imply

$\displaystyle\frac{\partial\varepsilon}{\partial p_{{}_{\rm F}}}=ap_{{}_{\rm
F}}^{3}\left(4(1+h)+p_{{}_{\rm F}}h^{\prime}\right),\quad\ \ \frac{\partial
P}{\partial p_{{}_{\rm F}}}=\frac{a}{3}p_{{}_{\rm
F}}^{3}\left(4(1+h)+6p_{{}_{\rm F}}h^{\prime}+p_{{}_{\rm
F}}^{2}h^{\prime\prime}\right),$ (A3)

so that

$\displaystyle c_{s}^{2}=\frac{1}{3}\left[1+\frac{5p_{{}_{\rm
F}}h^{\prime}+p_{{}_{\rm F}}^{2}h^{\prime\prime}}{4(1+h)+p_{{}_{\rm
F}}h^{\prime}}\right].$ (A4)

In the absence of $h^{\prime}$ and $h^{\prime\prime}$, the $h$ term does not
affect the sound velocity.

Let us first consider $h$ having a power-law (PL) behavior, $h=bp_{{}_{\rm
F}}^{\zeta}$, where $b$ is a constant. Then,

$\displaystyle c_{s}^{2}({\rm PL})=\frac{1}{3}\left[1+\frac{(4+\zeta)\zeta
h}{4+(4+\zeta)h}\right].$ (A5)

We see that a repulsive interaction, $h>0$, with $\zeta>0$ increases the sound
velocity above the conformal bound. An attractive interaction, $h<0$, can also
stiffen the equation of state if $\zeta<0$. Such density dependence is typical
in correlation effects near the Fermi surface. Pairing correlations in color
superconductivity (Alford et al., 2005) or baryonic correlations in quarkyonic
matter (McLerran & Reddy, 2019) are in a regime with $\zeta\sim 2$ to within
the density dependence of the gap (which is not known well in the non-
perturbative regime). For example, for the repulsive vector interaction
between quarks in the NJL model, $\sim g_{{}_{V}}n_{Q}^{2}$, one has
$\zeta=2$, and in the three-flavor NJL model employed in the text, we have
$h_{{}_{\rm NJL}}=(4/\pi^{2})g_{{}_{\rm V}}p_{{}_{\rm F}}^{2}$. Then

$\displaystyle c_{s}^{2}({\rm NJL})=\frac{1}{3}\left[1+\frac{6h_{{}_{\rm
NJL}}}{2+3h_{{}_{\rm NJL}}}\right],$ (A6)

which always exceeds 1/3 and approaches to 1 in the high density limit.

On the other hand, in the perturbative QCD (pQCD) regime at ultrahigh
densities, beyond the range one actually encounters in neutron stars,
$\zeta=0$, and an additional physical effect comes into play, the strong
coupling constant $\alpha_{s}$ in dense matter runs with $p_{{}_{\rm F}}$ (or
equivalently the chemical potential) as $\sim 1/\ln(p_{{}_{\rm
F}}/\Lambda_{\rm QCD})$ with $\Lambda_{\rm QCD}$ the QCD scale parameter
($\simeq$ 340 MeV for $N_{f}=3$). To first order in $\alpha_{s}$, this
dependence leads explicitly to

$\displaystyle\varepsilon_{{}_{\rm pQCD}}$ $\displaystyle=$
$\displaystyle\frac{3N_{f}}{4\pi^{2}}p_{{}_{\rm
F}}^{4}\left(1+\frac{2\alpha_{s}}{3\pi}\right);$ (A7)

thus $h_{{}_{\rm pQCD}}=2\alpha_{s}/3\pi$. To lowest order

$\displaystyle p_{{}_{\rm F}}\frac{\partial\alpha_{s}}{\partial p_{{}_{\rm
F}}}=-\frac{33-2N_{f}}{24\pi}\alpha_{s}^{2},$ (A8)

while $h^{\prime\prime}$ is of order $\alpha_{s}^{3}$. From Eq. (A4) we then
find for three flavors,

$\displaystyle c_{s}^{2}({\rm
pQCD})\simeq\frac{1}{3}\left(1+\frac{5}{4}p_{{}_{\rm
F}}h^{\prime}\right)=\frac{1}{3}\,\left(1-\frac{5(33-2N_{f})}{144\pi^{2}}\alpha_{s}^{2}\right)\stackrel{{\scriptstyle
N_{f}=3}}{{\longrightarrow}}\frac{1}{3}\,\left(1-\frac{15}{16}\left(\frac{\alpha_{s}}{\pi}\right)^{2}\right),$
(A9)

approaching the asymptotic value 1/3 from below. We see here that the
asymptotic approach to the conformal limit is the result of a competition
between the gluonic factor 33, which produces an approach from below, and the
fermionic term $2N_{f}$, which produces an approach from above. In QCD the
gluonic term wins, and the approach is from below.

Equation (A6), valid at intermediate densities relevant to neutron star
interiors, and Eq. (A9) for asymptotically high density, are smoothly joined
by a density dependent vector interaction, $g_{{}_{\rm V}}^{*}$, as introduced
in Song et al. (2019),

$\displaystyle g_{{}_{\rm V}}^{*}(p_{{}_{\rm
F}})=\frac{4\pi}{3}\frac{\alpha_{s}(p_{{}_{\rm F}})}{9m_{g}^{2}+8p_{{}_{\rm
F}}^{2}}.$ (A10)

Here $m_{g}\simeq 400$ MeV is a non-perturbative scale proportional to
$\Lambda_{\rm QCD}$, and the running coupling for $p_{{}_{\rm F}}<m_{g}$ is
frozen at $\alpha_{s}(m_{g})$ with $g_{{}_{\rm
V}}\equiv{4\pi\alpha_{s}(m_{g})}/27m_{g}^{2}$. With this density dependent
vector coupling, $h=(4/\pi^{2})g_{{}_{\rm V}}^{*}p_{{}_{\rm F}}^{2}$, which
approaches $h_{{}_{\rm NJL}}$ at $p_{{}_{\rm F}}\ll m_{g}$ and $h_{{}_{\rm
pQCD}}$ at $p_{{}_{\rm F}}\gg m_{g}$. The Fermi momentum at which
${c_{s}^{*}}^{2}-1/3$ changes sign from positive to negative is $\simeq
1.5\Lambda_{\rm QCD}$, corresponding to a baryon density $n\sim 10n_{0}$,

Quantum electrodynamics, in contrast to QCD, has no “charged” gluons, so that
it is not asymptotically free and the sound velocity exceeds 1/3. To see this
in detail, we note that the energy density in the relativistic limit to first
order in $\alpha_{e}$ is

$\displaystyle\varepsilon_{{}_{\rm QED}}$ $\displaystyle=$
$\displaystyle\frac{1}{4\pi^{2}}p_{{}_{\rm
F}}^{4}\left(1+\frac{\alpha_{e}}{2\pi}\right),$ (A11)

so that $h=\alpha_{e}/2\pi$. The evolution of $\alpha_{e}$ in QED is
equivalent to dropping the 33, letting $N_{f}=1$ and $\Lambda_{{}_{{}_{\rm
QCD}}}\rightarrow m_{e}$ in $\alpha_{s}$, and multiplying by a factor 2 which
arises from the algebraic difference of the interactions in QED and QCD. Then,
to lowest order, we have $p_{{}_{\rm F}}\partial\alpha_{e}/\partial p_{{}_{\rm
F}}=\alpha_{e}^{2}/6\pi$. The resulting sound speed, to order
$\alpha_{e}^{2}$, is given by

$\displaystyle c_{s}^{2}({\rm
QED})=\frac{1}{3}\,\left(1+\frac{5}{48}\left(\frac{\alpha_{e}}{\pi}\right)^{2}\right),$
(A12)

which exceeds 1/3 (and grows as the Fermi momentum increases). In the low
density non-relativistic limit, the sound speed of the electron gas is
$c_{s}^{2}({\rm QED})=(p_{{}_{\rm F}}/m_{e})^{2}/3$, which is smaller than
1/3.

##  Appendix B B. Exploring other interpolation ranges

While we take the interpolation range $n^{L}=1.5n_{0}$ and $n^{U}=3.5n_{0}$
for our main results in the text, we examine in this Appendix the dependence
of the NJL parameters and the resulting mass-radius relations on the matching
densities $n^{L}$ and $n^{U}$. We consider $n^{L}=(1.3,1.7)n_{0}$ and
$n^{U}=(3.5,4.5)n_{0}$ in the four possible combinations, continuing to use
the $\chi$EFT equation of state in the nuclear regime. Figure 10 shows the
allowed ranges of $(g_{{}_{\rm V}},H)$, with the parameter sets Aχ-Dχ in the
main text marked as guides. As we can see, changing the interpolation range
modifies the domains of $(g_{{}_{\rm V}},H)$ by only 10-20%.

Increasing the size of interpolating interval should in principle
monotonically extend the window for $(g_{{}_{\rm V}},H)$, as well as allow a
larger mass range. For example, the results of $(n^{L},n^{U})=(1.3,4.5)n_{0}$
should contain the results of $(1.7,3.5)n_{0}$ as special cases. However in
our modeling where we use simple polynomials which cannot cover all allowable
curves, this is not the case. We use simple structureless polynomials so as
not to introduce exotic structures by hand, and yet the interpolation results
in a highly nontrivial structure in the equation of state, e.g., the peak in
sound velocity. Owing to our not exploring all possible interpolation curves,
the maximum mass we find in a given interpolation window should be understood
as a lower bound of the maximum mass.

Shown in Fig. 11 are $c_{s}^{2}$ vs. $n$ (upper panels), $M$ vs. $R$ (middle
panels), and $M$ vs. $n$ (lower panels), for the matching densities
$n^{L}=1.5n_{0}$ and $n^{U}=(3.5,4.0,4.5)n_{0}$. For given $g_{{}_{\rm
V}}=1.1$, and $1.2G$, the lower (bold lines) and upper (thin lines) bounds on
$H$ are as shown in Fig. 10. The plots for intermediate values of $H$ lie
between these lines.

There are clear correlations between the locations of the peaks in the sound
velocity, the radii of $1.4M_{\odot}$ neutron stars, and the growth of the
neutron star mass as a function of the central density. A peak in sound
velocity at lower density indicates rapid stiffening which results in a larger
radius. In the upper panels of Fig. 11, the locations of the peaks change by
$\sim 1$-$1.5n_{0}$ as $H$ is varied at given $g_{{}_{\rm V}}$, and the radii
change by $\sim$ 0.3-0.7 km, as shown in the middle panels. The central value
of $H$ leads to $R_{1.4}\simeq$ 12.2-12.4 km, and with the variation of $H$,
the $M$-$R$ relations remain largely consistent with the NICER data.

Figure 10.—  The maximum mass $M_{\rm max}$ as a function of $(g_{{}_{\rm
V}},H)$ for the four combinations of matching densities $n^{L}$= (1.3, 1.7)
$n_{0}$ and $n^{U}$ = (3.5, 4.5) $n_{0}$. The parameter sets
$A_{\chi}$-$D_{\chi}$ in the main text are indicated only as guides.

Figure 11.—  The relations among the sound velocity, mass-radius, and core
density of neutron stars; $c_{s}^{2}$ vs. $n$ (upper panels), $M$ vs $R$
(middle panels), and $M$ vs. $n$ (lower panels), for matching densities
$n^{L}$ =1.5 $n_{0}$ and $n^{U}$=3.5, 4.0, and 4.5 $n_{0}$. For given
$g_{{}_{\rm V}}=1.1$, and 1.2 $G$, the coupling $H$ has a lower (shown as
bold) and an upper (shown as thin) bound (see Fig. 10). The plots for
intermediate $H$ lie basically between the bold and thin lines. The error
bands are as those shown in Figs. 5 and 7.

## References

  * Akmal et al. (1998) Akmal, A., Pandharipande, V., & Ravenhall, D. 1998, Phys. Rev. C, 58, 1804, doi: 10.1103/PhysRevC.58.1804
  * Alford et al. (2005) Alford, M., Braby, M., Paris, M. W., & Reddy, S. 2005, Astrophys. J., 629, 969, doi: 10.1086/430902
  * Annala et al. (2020) Annala, E., Gorda, T., Kurkela, A., Nattila, J., & Vuorinen, A. 2020, Nature Phys., doi: 10.1038/s41567-020-0914-9
  * Baym et al. (2019) Baym, G., Furusawa, S., Hatsuda, T., Kojo, T., & Togashi, H. 2019, Astrophys. J., 885, 42, doi: 10.3847/1538-4357/ab441e
  * Baym et al. (2018) Baym, G., Hatsuda, T., Kojo, T., et al. 2018, Rept. Prog. Phys., 81, 056902, doi: 10.1088/1361-6633/aaae14
  * Bazavov et al. (2014) Bazavov, A., et al. 2014, Phys. Rev. D, 90, 094503, doi: 10.1103/PhysRevD.90.094503
  * De Rujula et al. (1975) De Rujula, A., Georgi, H., & Glashow, S. L. 1975, Phys. Rev. D, 12, 147, doi: 10.1103/PhysRevD.12.147
  * Dietrich et al. (2020) Dietrich, T., Coughlin, M. W., Pang, P. T. H., et al. 2020, Science, 370, 1450, doi: 10.1126/science.abb4317
  * Drischler et al. (2021a) Drischler, C., Han, S., Lattimer, J. M., et al. 2021a, Phys. Rev. C, 103, 045808, doi: 10.1103/PhysRevC.103.045808
  * Drischler et al. (2021b) Drischler, C., Holt, J. W., & Wellenhofer, C. 2021b, Ann. Rev. Nucl. Part. Sci., 71, 403, doi: 10.1146/annurev-nucl-102419-041903
  * Fonseca et al. (2021) Fonseca, E., et al. 2021, Astrophys. J. Lett., 915, L12, doi: 10.3847/2041-8213/ac03b8
  * Fukushima et al. (2020) Fukushima, K., Kojo, T., & Weise, W. 2020, Phys. Rev. D, 102, 096017, doi: 10.1103/PhysRevD.102.096017
  * Greif et al. (2019) Greif, S. K., Raaijmakers, G., Hebeler, K., Schwenk, A., & Watts, A. L. 2019, Mon. Not. Roy. Astron. Soc., 485, 5363, doi: 10.1093/mnras/stz654
  * Hatsuda & Kunihiro (1994) Hatsuda, T., & Kunihiro, T. 1994, Phys. Rept., 247, 221, doi: 10.1016/0370-1573(94)90022-1
  * Hatsuda et al. (2006) Hatsuda, T., Tachibana, M., Yamamoto, N., & Baym, G. 2006, Phys. Rev. Lett., 97, 122001, doi: 10.1103/PhysRevLett.97.122001
  * Kojo (2021a) Kojo, T. 2021a, Phys. Rev. D, 104, 074005, doi: 10.1103/PhysRevD.104.074005
  * Kojo (2021b) —. 2021b, AAPPS Bull., 31, 11, doi: 10.1007/s43673-021-00011-6
  * Kojo et al. (2015) Kojo, T., Powell, P. D., Song, Y., & Baym, G. 2015, Phys. Rev. D, 91, 045003, doi: 10.1103/PhysRevD.91.045003
  * Kojo & Suenaga (2021) Kojo, T., & Suenaga, D. 2021.  https://arxiv.org/abs/2110.02100
  * Kunihiro (1991) Kunihiro, T. 1991, Phys. Lett. B, 271, 395, doi: 10.1016/0370-2693(91)90107-2
  * Legred et al. (2021) Legred, I., Chatziioannou, K., Essick, R., Han, S., & Landry, P. 2021, Phys. Rev. D, 104, 063003, doi: 10.1103/PhysRevD.104.063003
  * Lonardoni et al. (2020) Lonardoni, D., Tews, I., Gandolfi, S., & Carlson, J. 2020, Phys. Rev. Res., 2, 022033, doi: 10.1103/PhysRevResearch.2.022033
  * Masuda et al. (2013a) Masuda, K., Hatsuda, T., & Takatsuka, T. 2013a, Astrophys. J., 764, 12, doi: 10.1088/0004-637X/764/1/12
  * Masuda et al. (2013b) —. 2013b, PTEP, 2013, 073D01, doi: 10.1093/ptep/ptt045
  * Masuda et al. (2016) —. 2016, Eur. Phys. J. A, 52, 65, doi: 10.1140/epja/i2016-16065-6
  * McLerran & Reddy (2019) McLerran, L., & Reddy, S. 2019, Phys. Rev. Lett., 122, 122701, doi: 10.1103/PhysRevLett.122.122701
  * Miller et al. (2021) Miller, M. C., et al. 2021, Astrophys. J. Lett., 918, L28, doi: 10.3847/2041-8213/ac089b
  * Oka & Yazaki (1980) Oka, M., & Yazaki, K. 1980, Phys. Lett. B, 90, 41, doi: 10.1016/0370-2693(80)90046-5
  * Pang et al. (2021) Pang, P. T. H., Tews, I., Coughlin, M. W., et al. 2021.  https://arxiv.org/abs/2105.08688
  * Park et al. (2020) Park, A., Lee, S. H., Inoue, T., & Hatsuda, T. 2020, Eur. Phys. J. A, 56, 93, doi: 10.1140/epja/s10050-020-00078-z
  * Raaijmakers et al. (2021) Raaijmakers, G., Greif, S. K., Hebeler, K., et al. 2021, Astrophys. J. Lett., 918, L29, doi: 10.3847/2041-8213/ac089a
  * Riley et al. (2021) Riley, T. E., et al. 2021, Astrophys. J. Lett., 918, L27, doi: 10.3847/2041-8213/ac0a81
  * Schäfer & Wilczek (1999) Schäfer, T., & Wilczek, F. 1999, Phys. Rev. Lett., 82, 3956, doi: 10.1103/PhysRevLett.82.3956
  * Song et al. (2019) Song, Y., Baym, G., Hatsuda, T., & Kojo, T. 2019, Phys. Rev. D, 100, 034018, doi: 10.1103/PhysRevD.100.034018
  * Tews et al. (2018) Tews, I., Carlson, J., Gandolfi, S., & Reddy, S. 2018, Astrophys. J., 860, 149, doi: 10.3847/1538-4357/aac267
  * Togashi et al. (2017) Togashi, H., Nakazato, K., Takehara, Y., et al. 2017, Nucl. Phys. A, 961, 78, doi: 10.1016/j.nuclphysa.2017.02.010


GenEx – Exclusive Meson Generator

Maciej Trzebiński111This was partially supported by Polish Ministry of Science
and Higher Education under the Mobility Plus programme (1285/MOB/IV/2015/0)
and by Polish National Science Centre under the grant number
UMO-2015/17/D/ST2/03530.

Institute of Nuclear Physics Polish Academy of Sciences, Radzikowskiego 152
st., Cracow, Poland

> Exclusive light meson production processes possible to be measured at RHIC
> and LHC accelerators are briefly described. This includes Pomeron and photon
> induced continuum of pions and kaons as well as $f_{0}$ and $\rho^{0}$
> resonant production. Next, GenEx, a new Monte Carlo generator recently
> developed in Cracow, is presented. Its purpose is to generate resonant and
> non-resonant pion and kaon production and diffractive Bremsstrahlung. Its
> future development will take into account spin (polarization) effects and
> simulate corrections due to absorption and re-scattering.

> PRESENTED AT
>
>
>
>
> Presented at EDS Blois 2017, Prague,
> Czech Republic, June 26-30, 2017

##  1 Introduction

Collisions at hadron accelerators are dominated by soft processes. Absence of
a hard scale in these events prevents one from using perturbation theory.
Instead, in order to calculate the properties of the produced particles (e.g.
energy or angular distributions), one has to use approximative methods.

In about a half of collisions at LHC or RHIC one or both outgoing protons stay
intact. Such final state is possible due to the exchange of a colourless
object: photon (electromagnetic) or Pomeron (strong interaction). Such an
exchange results in a presence of the rapidity gap between the centrally
produced system and scattered protons. Thus, such events are of diffractive
nature.

Diffractive studies are one of the important parts of the physics programme
for the LHC experiments. A particularly interesting classes are exclusive
processes, where all centrally-produced particles are detected. An example is
the process of pion pair production $pp\rightarrow p\pi^{+}\pi^{-}p$, which
can be used for studies of low-mass resonances, including searches for
glueballs.

In order to include detector effects (acceptance, efficiency) in theory-to-
data comparison Monte Carlo generators are needed. There are few MC generators
which are specialized in simulation of the soft exclusive processes, e.g.
SuperCHIC [1] or DIME [2]. In this paper GenEx [3, 4], a tool complementary to
the existing ones in terms of implemented processes and calculation methods,
will be briefly described.

##  2 Soft Exclusive Production at the LHC

Diffractive protons are usually scattered at very small angles (fraction of
milliradian). In order to measure them, special devices that allow detector
presence in close vicinity of the beam are required. At the LHC Roman pots are
commonly used. ATLAS [5] has two systems of such detectors: ALFA [6, 7] and
AFP [8]. Similar set is installed around CMS/TOTEM Interaction Point [9].
Roman pots are usually located far away (few hundreds meters) from the
Interaction Point (IP). It means that there are several LHC magnets between
them and the IP. The settings of these magnets, commonly called machine
optics, play a key role in the exclusive analysis. The detailed description of
the properties of optics sets used at the LHC can be found in [10]. In
addition, in this reference a geometric acceptance of AFP and ALFA detectors
is presented.

###  2.1 Non-resonant Pion Pair Production

Diagram for the non-resonant pion pair production (also called continuum) is
shown in Fig. 1 (left); a Pomeron is “emitted” from each proton resulting in a
four particles present in final state: scattered protons and (central) pions.
The object exchanged in the $t$-channel is an off-shell pion. The model
implemented in GenEx is based on calculations done in Ref. [11]. Feasibility
studies for the ATLAS detector are shown in [12].

Figure 1: Diagrams of: non-resonant (left) and $f_{0}$ resonance (right) pion
pair production.

###  2.2 Pion Pair Production via $f_{0}$ Resonance

Diagram for the exclusive pion pair production via the $f_{0}$ resonance is
shown in Fig. 1 (right). In order to calculate it, one must know the Pomeron-
Pomeron-meson vertices. These elements have been discussed in detail in Ref.
[13]. It was shown there that different couplings structures with
corresponding orbital angular momentum and spin of two Pomerons are possible.

###  2.3 Photon-induced Continuum and $\rho^{0}$ Photoproduction

The dominant diagram of the exclusive pion pair continuum production is a
Pomeron-induced one. However, the production of a photon-induced continuum
(see Fig. 2 (left)) is also possible. It should be taken into account in the
Monte Carlo generators as it contributes up to the several percent of the
exclusive pion pair production, see Ref. [14].

Figure 2: Diagrams of: non-resonant (left) and $\rho^{0}$ (right) pion pair
photo-production.

Another interesting class of photon-induced processes is a resonant $\rho^{0}$
photoproduction. Feynman diagram is shown Fig. 2 (right) whereas details are
discussed in [14].

###  2.4 Diffractive Bremsstrahlung

Diffractive Bremsstrahlung is typically an electromagnetic process (see Fig. 3
(left)). However, as discussed in [15], high energy photons can be radiated in
the elastic proton-proton scattering (Fig. 3 (right)). This idea was extended
in [16] by introducing the proton form-factor into the calculations and by
considering also other mechanisms leading to the $pp\gamma$ final state, e.g.
a virtual photon re-scattering. Recently, diffractive Bremsstrahlung was
implemented into the GenEx Monte Carlo generator [17].

Figure 3: Diagrams of electromagnetic (left) and diffractive (right)
Bremsstrahlung. The blobs represent various mechanisms of photon emission.

##  3 GenEx Generator

GenEx [3] is a C++ class structure for the construction of a Monte Carlo event
generators which can produce unweighted events within the relativistic phase
space. Generator is self-adapting to the provided matrix element and
acceptance cuts. Current version contain resonant and non-resonant exclusive
meson production processes with vector and tensor Pomeron as well as
diffractive Bremsstrahlung production [18]. In near future spin (polarization)
effects and absorption/re-scattering corrections will be added.

Current version of generator is based on several core classes:

  * •

TFoam – class of adaptive Monte Carlo simulator,

  * •

TDensity – class with integrand function (calculated accordingly to a given
event),

  * •

TEventMaker2toN – generates two leading particles and a central blob which
then is decayed by TDecay into $N-2$ remaining particles.

Recently, a significant effort is taken in order to simplify the
implementation of new processes. For reference see [19].

##  4 Summary

Exclusive light meson production should be measurable at LHC and RHIC. A Monte
Carlo generator containing continuum and resonances would be useful to
understand recorded data. Such tool will be especially needed in order to
include detector effects (acceptance, efficiency) in theory-data comparison.

One of such generators is GenEx. Its first version is already available for
tests (see [18]) and contains: the generator structure, pion/kaon continuum
and diffractive Bremsstrahlung processes. The developments are ongoing – in
the near future a new version of GenEx would allow to: generate resonant
production ($f_{0}(500)$, $f_{0}(980)$, $f_{0}(1370)$, $f_{0}(1500)$,
$f_{2}(1270)$, $f_{2}^{\prime}(1520)$ and $\rho_{0}$), take into account spin
(polarization) effects and simulate corrections due to absorption and re-
scattering.

## References

  * [1] L. Harland-Lang et al., SuperCHIC: Monte Carlo Event Generator for Central Exclusive Production, https://superchic.hepforge.org/.
  * [2] L. Harland-Lang et al., Dime Monte Carlo: event generator for exclusive meson pair production via double Pomeron exchange, http://dimemc.hepforge.org/.
  * [3] R. A. Kycia, J. Chwastowski, R. Staszewski, J. Turnau, GenEx: A simple generator structure for exclusive processes in high energy collisions, arXiv:1411.6035.
  * [4] P. Erland, R. Staszewski, M. Trzebinski, R. Kycia, Elastic Hadron Scattering in Various Pomeron Models, Acta Phys. Pol. B 48 (2017) 981.
  * [5] ATLAS Collaboration, The ATLAS Experiment at the CERN Large Hadron Collider, J. Instrum. 3 (2008) S08003.
  * [6] ATLAS Luminosity and Forward Physics Community, Technical Design Report, CERN-LHCC-2008-004, https://cds.cern.ch/record/1095847.
  * [7] S. Abdel Khalek et al., The ALFA Roman Pot Detectors of ATLAS, JINST 11 (2016) P11013.
  * [8] ATLAS Collaboration, Letter of Intent for the Phase-I Upgrade of the ATLAS Experiment, CERN-LHCC-2011-012, https://cds.cern.ch/record/1402470.
  * [9] TOTEM Collaboration, Technical Design Report, CERN-LHCC-2004-002,
TOTEM Collaboration, TOTEM Upgrade Proposal, CERN-LHCC-2013-009.

  * [10] M. Trzebinski, Machine Optics Studies for the LHC Measurements, in proceedings of XXXIV-th SPIE Joint Symposium Wilga vol. 9290 26. Proc. SPIE Int. Soc. Opt. Eng., 0277, 786; arXiv:1408.1836
  * [11] P. Lebiedowicz, O. Nachtmann, A. Szczurek, Central exclusive diffractive production of $\pi^{+}\pi^{-}$ continuum, scalar and tensor resonances in $pp$ and $p\overline{p}$ scattering within tensor pomeron approach, Phys. Rev. D 93 (2016) 054015.
  * [12] R. Staszewski, P. Lebiedowicz, M. Trzebinski, J. Chwastowski, A. Szczurek, Exclusive $\pi^{+}\pi^{-}$ Production at the LHC with Forward Proton Tagging, Acta Physica Polonica B 42 (2011) 1861.
  * [13] P. Lebiedowicz, O. Nachtmann, A. Szczurek, Exclusive central diffractive production of scalar and pseudoscalar mesons tensorial vs. vectorial pomeron, Annals Phys. 344 (2014) 301-339.
  * [14] P. Lebiedowicz, O. Nachtmann, A. Szczurek, The $\rho^{0}$ and Drell-Soding contributions to central exclusive production of $\pi^{+}\pi^{-}$ pairs in proton-proton collisions at high energies, Phys. Rev. D 91 (2015) 074023.
  * [15] V. A. Khoze, J. W. Lamsa, R. Orava, M. G. Ryskin, Forward Physics at the LHC: Detecting Elastic pp Scattering by Radiative Photons, JINST 6 (2011) P01005.
  * [16] P. Lebiedowicz and A. Szczurek, Exclusive diffractive photon bremsstrahlung at the LHC, Phys. Rev. D 87 (2013) 114013.
  * [17] J. J. Chwastowski, S. Czekierda, R. Staszewski, M. Trzebinski, Diffractive Bremsstrahlung at High-$\beta^{*}$ LHC Case Study, Eur. Phys. J. C 77 (2017) 216.
  * [18] https://github.com/rkycia/GenEx.
  * [19] P. Buglewicz, K. Cieśla, Calculations in tensor pomeron model – an object-oriented implementation, talk at Challenges in Photon Induced Interactions, Cracow, 5-8 Sep. 2017.