# Hadron-quark Pasta Phase in Massive Neutron Stars

Min Ju School of Physics, Nankai University, Tianjin 300071, China  Jinniu Hu
School of Physics, Nankai University, Tianjin 300071, China;
<EMAIL_ADDRESS>Hong Shen School of Physics, Nankai University,
Tianjin 300071, China<EMAIL_ADDRESS>
###### Abstract

The structured hadron-quark mixed phase, known as the pasta phase, is expected
to appear in the core of massive neutron stars. Motivated by the recent
advances in astrophysical observations, we explore the possibility of the
appearance of quarks inside neutron stars and check its compatibility with
current constraints. We investigate the properties of the hadron-quark pasta
phases and their influences on the equation of state (EOS) for neutron stars.
In this work, we extend the energy minimization (EM) method to describe the
hadron-quark pasta phase, where the surface and Coulomb contributions are
included in the minimization procedure. By allowing different electron
densities in the hadronic and quark matter phases, the total electron chemical
potential with the electric potential remains constant, and local $\beta$
equilibrium is achieved inside the Wigner–Seitz cell. The mixed phase
described in the EM method shows the features lying between the Gibbs and
Maxwell constructions, which is helpful for understanding the transition from
the Gibbs construction (GC) to the Maxwell construction (MC) with increasing
surface tension. We employ the relativistic mean-field model to describe the
hadronic matter, while the quark matter is described by the MIT bag model with
vector interactions. It is found that the vector interactions among quarks can
significantly stiffen the EOS at high densities and help enhance the maximum
mass of neutron stars. Other parameters like the bag constant can also affect
the deconfinement phase transition in neutron stars. Our results show that
hadron-quark pasta phases may appear in the core of massive neutron stars that
can be compatible with current observational constraints.

Neutron stars — Nuclear astrophysics — neutron star cores — gravitational
waves

##  1 Introduction

The appearance of deconfined quark matter, which is expected to appear in the
core of massive neutron stars, has received increasing attention recently
because of its relevance to astrophysical observations (Lattimer & Prakash,
2016; Baym et al., 2018; Annala et al., 2020). In the last decade, several
breakthrough discoveries in astronomy provided valuable information and
constraints on the properties of neutron stars. The precise mass measurements
of PSR J1614-2230 ($1.908\pm 0.016M_{\odot}$; Demorest et al. 2010; Fonseca et
al. 2016; Arzoumanian et al. 2018), PSR J0348+0432 ($2.01\pm 0.04M_{\odot}$;
Antoniadis et al. 2013), and PSR J0740+6620 ($2.08\pm 0.07M_{\odot}$;
Cromartie et al. 2020; Fonseca et al. 2020) constrain the maximum neutron-star
mass $M_{\rm{max}}$ to be larger than about $2M_{\odot}$, which poses a
challenge to our understanding of the equation of state (EOS) of superdense
matter. The recent observations by the Neutron Star Interior Composition
Explorer (NICER) provided a simultaneous measurement of the mass and radius
for PSR J0030+0451, which was reported to have a mass of
$1.44_{-0.14}^{+0.15}M_{\odot}$ with a radius of $13.02_{-1.06}^{+1.24}$ km
(Miller et al., 2019) and a mass of $1.34_{-0.16}^{+0.15}M_{\odot}$ with a
radius of $12.71_{-1.19}^{+1.14}$ km (Riley et al., 2019) by two independent
groups. The new measurements by NICER for the most massive known neutron star,
PSR J0740+6620 ($2.08\pm 0.07M_{\odot}$), showed that it has a radius of
$13.7_{-1.5}^{+2.6}$ km by Miller et al. (2021) and $12.39_{-0.98}^{+1.30}$ km
by Riley et al. (2021). In particular, the discovery of gravitational waves
from a binary neutron-star merger event GW170817 has opened a new era of
multimessenger astronomy (Abbott et al., 2017). Based on the observations of
GW170817, the tidal deformability of a canonical $1.4M_{\odot}$ neutron star
was estimated to be $70<\Lambda_{1.4}<580$ and the corresponding radius was
inferred to be $10.5<R_{1.4}<13.3$ km (Abbott et al., 2018). More recently,
the gravitational-wave events, GW190425 (Abbott et al., 2020a) and GW190814
(Abbott et al., 2020b), were reported by LIGO and Virgo Collaborations. The
total mass of the GW190425 system is as large as $3.4^{+0.3}_{-0.1}M_{\odot}$,
which is even more massive than any neutron-star binary observed so far; hence
its gravitational-wave analyses may offer valuable information for the EOS at
high densities and possible phase transitions inside neutron stars. The
GW190814 event was detected from a compact binary coalescence involving a
22.2–24.3$M_{\odot}$ black hole and a 2.50–2.67$M_{\odot}$ compact object that
could be either the heaviest neutron star or the lightest black hole ever
observed (Abbott et al., 2020b). The discovery of GW190814 has triggered many
theoretical efforts exploring the nature of the secondary object and its
implications for high-density EOS (Fattoyev et al., 2020; Huang et al., 2020;
Tews et al., 2021). The gravitational-wave analyses (Essick & Landry, 2020;
Tews et al., 2021) suggested that the secondary in GW190814 is more likely to
be a black hole, but its possibility as a neutron star cannot be ruled out.
Several studies have suggested that the secondary in GW190814 may be a rapidly
rotating neutron star (Li et al., 2020; Most et al., 2020; Tsokaros et al.,
2020; Zhang & Li, 2020). Moreover, it could be also considered as a heavy
neutron star containing deconfined quark matter, where the inclusion of quarks
has significant impact on the observations (Tan et al., 2020; Demircik et al.,
2021; Dexheimer et al., 2021). In light of these recent developments in
astronomy, it would be interesting and informative to explore possible
structures of the hadron-quark mixed phase in massive neutron stars.

For the description of hadron-quark phase transition in neutron stars, both
Gibbs and Maxwell constructions are often used depending on the surface
tension at the interface (Bhattacharyya et al., 2010). In the limit of zero
surface tension, the Gibbs equilibrium conditions are satisfied between the
two coexisting phases, and only global charge neutrality is imposed in the
mixed phase (Glendenning, 1992; Schertler et al., 2000; Yang & Shen, 2008; Xu
et al., 2010; Wu & Shen, 2017). On the other hand, the Maxwell construction
(MC) is valid for sufficiently large surface tension, where local charge
neutrality is enforced and the transition takes place at constant pressure
(Bhattacharyya et al., 2010; Han et al., 2019; Wu & Shen, 2019). It is
noteworthy that only bulk contributions are involved in the Gibbs and Maxwell
constructions, whereas the finite-size effects like surface and Coulomb
energies are neglected. In a more realistic case, due to the competition
between surface and Coulomb energies, some geometric structures may be formed
in the hadron-quark mixed phase, known as hadron-quark pasta phases
(Heiselberg et al., 1993; Endo et al., 2006; Maruyama et al., 2007; Yasutake
et al., 2014; Spinella et al., 2016; Weber et al., 2019; Wu & Shen, 2019).
This structured mixed phase is analogous to nuclear pasta phase in the inner
crust of neutron stars. Several methods have been developed to study the
properties of hadron-quark pasta phases. In the coexisting phases (CP) method
(Wu & Shen, 2019), the hadronic and quark phases are assumed to satisfy the
Gibbs conditions for phase equilibrium, while the surface and Coulomb energies
are taken into account perturbatively. A more realistic description of the
pasta phase has been developed in a series of works (Endo et al., 2006;
Maruyama et al., 2007; Yasutake et al., 2014), where the Thomas–Fermi
approximation was used to describe the density profiles of hadrons and quarks
in the Wigner–Seitz cell. For simplicity, the particle densities in the two
coexisting phases are generally assumed to be spatially constant, and the
charge screening effect is neglected. In our previous works (Wu & Shen, 2019;
Ju et al., 2021), we proposed an energy minimization (EM) method for improving
the treatment of surface and Coulomb energies, which play a key role in
determining the structure of the pasta phases. By incorporating the surface
and Coulomb contributions in the EM procedure, one can derive the equilibrium
conditions for coexisting phases that are different from the Gibbs conditions
used in the CP method. In the present work, we further extend the EM method
for exploring the hadron-quark pasta phases in massive neutron stars.

In the EM method, the hadron-quark pasta phase is described within the
Wigner–Seitz approximation, where the whole space is divided into equivalent
cells with a geometric symmetry. The hadronic and quark phases in a charge-
neutral cell are assumed to be separated by a sharp interface. In the more
realistic description of Tatsumi et al. (2003), there is a thin boundary layer
between two separate bulk phases, which can lead to a difference in the
electric potential between the two phases. As a result, the electron densities
in the hadronic and quark phases are allowed to be different from each other,
whereas the electron chemical potentials are equal. According to this
argument, the MC could be stable although the electron densities in two
separate bulk phases are different. The difference in the electron densities
between two separate phases can be understood as a result of the charge
screening effect. In the present study, we would like to evaluate how
important the charge screening effect is for the hadron-quark pasta phases.
Meanwhile, the results in the EM method will be compared to those obtained
with a uniform electron gas in the CP method.

To describe the hadronic phase, we employ the relativistic mean-field (RMF)
model and choose the recently proposed BigApple parameterization (Fattoyev et
al., 2020). In the past decades, the RMF approach based on various energy
density functionals has been successfully applied in the description of finite
nuclei and infinite nuclear matter. Several popular RMF models, such as NL3
(Lalazissis et al., 1997), TM1 (Sugahara & Toki, 1994), and IUFSU (Fattoyev et
al., 2010), have been widely used in astrophysical applications, since they
can not only reproduce experimental data of finite nuclei but also predict
large enough neutron-star masses. The BigApple model was proposed by Fattoyev
et al. (2020) after the discovery of GW190814, in which they considered
various constraints from astrophysical observations as well as the ground-
state properties of finite nuclei. The BigApple model predicts a maximum
neutron-star mass of $M_{\rm{max}}=2.6M_{\odot}$, while the resulting radius
and tidal deformability of neutron stars are consistent with GW170817 and
NICER observations. It is well known that nuclear symmetry energy and its
slope play an important role in understanding the properties of neutron stars
(Oertel et al., 2017; Ji et al., 2019). There exists a positive correlation
between the symmetry energy slope $L$ and the neutron-star radius (Alam et
al., 2016). The NL3 and TM1 models have a rather large slope parameter, which
results in too large a radius and tidal deformability for a canonical
$1.4M_{\odot}$ neutron star as compared to the estimations from astrophysical
observations (Ji et al., 2019). In the present work, we prefer to employ the
BigApple model with a small slope $L$, which is more consistent with the
analysis of the GW170817 event. Generally, the appearance of hyperons at high
densities would considerably soften the EOS and reduce the maximum neutron-
star mass (Oertel et al., 2017). However, the influence of hyperons may be
suppressed by incorporating quark degrees of freedom or introducing additional
repulsion for hyperons (Lonardoni et al., 2015; Gomes et al., 2019).
Currently, there are large uncertainties in the hyperon–nucleon and
hyperon–hyperon interactions due to limited experimental data. For simplicity,
we do not include hyperons in the present calculations and focus on the
influence of deconfined quarks in the core of neutron stars.

We utilize a modified MIT bag model with vector interactions, often referred
to as the vMIT model (Gomes et al., 2019; Han et al., 2019), to describe the
quark phase. It has been shown in the literature that including vector
interactions among quarks can significantly stiffen the EOS at high densities
and help to produce massive neutron stars in agreement with the astronomical
observations (Klähn & Fischer, 2015). The vector interaction in the vMIT model
is introduced via the exchange of a vector meson that is analogous to the
$\omega$ meson in the Walecka model (Lopes et al., 2021). The vector coupling
constant is usually treated as a free parameter and a universal coupling is
assumed for all quark flavors. In Han et al. (2019) and Gomes et al. (2019),
the influences of the vector interaction on the properties of compact stars
were studied using the Gibbs and Maxwell constructions to describe the hadron-
quark mixed phase without considering possible geometric structures. In the
present work, we intend to investigate how the hadron-quark pasta phases are
affected by the vector interactions among quarks within the EM method.

We have two aims in this paper. The first is to explore the possibility of the
appearance of hadron-quark pasta phases in the core of massive neutron stars.
By comparing with current constraints, we investigate the compatibility
between the hadron-quark phase transition and astrophysical observations, as
well as the influence of the quark vector interactions. The second aim is to
extend the EM method for describing the hadron-quark pasta phases with the
charge screening effect. By allowing different electron densities in separate
bulk phases, the local $\beta$ equilibrium is achieved inside the Wigner–Seitz
cell. Although in more precise Thomas–Fermi calculations the electron density
should depend on the position, a simplified treatment of the charge screening
effect with different constant electron densities in the two bulk phases is
helpful for understanding the transition from the GC to the MC. Comparing with
the GC, the EM method used in this work involves the effects of finite size
and charge screening, which are essential for determining the structure of the
hadron-quark mixed phase.

This paper is arranged as follows. In Section 2, we briefly describe the
framework for describing the hadronic and quark phases. The EM method for the
hadron-quark pasta phases is presented in Section 3. In Section 4, we discuss
numerical results of the hadron-quark pasta phases and neutron-star
properties. Section 5 is devoted to a summary.

##  2 Models for hadronic and quark phases

In this section, we briefly describe the RMF model for the hadronic matter and
the vMIT model for the quark matter. In addition, we explain the model
parameters used in our calculations.

###  2.1 Hadronic Phase

We employ the RMF model with the BigApple parameterization to describe the
hadronic matter, where nucleons interact through the exchange of various
mesons including the isoscalar–scalar $\sigma$ meson, the isoscalar–vector
$\omega$ meson, and the isovector–vector $\rho$ meson. The Lagrangian density
for hadronic matter consisting of nucleons ($p$ and $n$) and leptons ($e$ and
$\mu$) is written as

$\displaystyle\mathcal{L}$ $\displaystyle=$
$\displaystyle\sum_{i=p,n}\bar{\psi}_{i}\left\\{i\gamma_{\mu}\partial^{\mu}-M-g_{\sigma}\sigma-\gamma_{\mu}\left[g_{\omega}\omega^{\mu}+\frac{g_{\rho}}{2}\tau_{a}\rho^{a\mu}\right]\right\\}\psi_{i}$
(1)
$\displaystyle+\frac{1}{2}\partial_{\mu}\sigma\partial^{\mu}\sigma-\frac{1}{2}m^{2}_{\sigma}\sigma^{2}-\frac{1}{3}g_{2}\sigma^{3}-\frac{1}{4}g_{3}\sigma^{4}$
$\displaystyle-\frac{1}{4}W_{\mu\nu}W^{\mu\nu}+\frac{1}{2}m^{2}_{\omega}\omega_{\mu}\omega^{\mu}+\frac{1}{4}c_{3}\left(\omega_{\mu}\omega^{\mu}\right)^{2}$
$\displaystyle-\frac{1}{4}R^{a}_{\mu\nu}R^{a\mu\nu}+\frac{1}{2}m^{2}_{\rho}\rho^{a}_{\mu}\rho^{a\mu}+\Lambda_{\rm{v}}\left(g_{\omega}^{2}\omega_{\mu}\omega^{\mu}\right)\left(g_{\rho}^{2}\rho^{a}_{\mu}\rho^{a\mu}\right)$
$\displaystyle+\sum_{l=e,\mu}\bar{\psi}_{l}\left(i\gamma_{\mu}\partial^{\mu}-m_{l}\right)\psi_{l},$

where $W^{\mu\nu}$ and $R^{a\mu\nu}$ represent the antisymmetric field tensors
for $\omega^{\mu}$ and $\rho^{a\mu}$, respectively. Under the mean-field
approximation, the meson fields are treated as classical fields, which are
denoted by $\sigma=\left\langle\sigma\right\rangle$,
$\omega=\left\langle\omega^{0}\right\rangle$, and
$\rho=\left\langle\rho^{30}\right\rangle$. These mean fields can be obtained
by solving a set of coupled equations in the RMF model.

For the hadronic matter in $\beta$ equilibrium, the chemical potentials
satisfy the relations $\mu_{p}=\mu_{n}-\mu_{e}$ and $\mu_{\mu}=\mu_{e}$. At
zero temperature, the chemical potentials of nucleons and leptons are
expressed as

$\displaystyle\mu_{i}$ $\displaystyle=$
$\displaystyle\sqrt{{k_{F}^{i}}^{2}+{M^{\ast}}^{2}}+g_{\omega}\omega+\frac{g_{\rho}}{2}\tau_{3}\rho,\hskip
14.22636pti=p,n,$ (2) $\displaystyle\mu_{l}$ $\displaystyle=$
$\displaystyle\sqrt{{k_{F}^{l}}^{2}+m_{l}^{2}},\hskip 85.35826ptl=e,\mu,$ (3)

with $\tau_{3}=+1$ and $-1$ for protons and neutrons, respectively. The
effective nucleon mass is defined as $M^{\ast}=M+g_{\sigma}{\sigma}$. The
total energy density and pressure in hadronic matter are written as

$\displaystyle\varepsilon_{\rm{HP}}$ $\displaystyle=$
$\displaystyle\sum_{i=p,n}\varepsilon^{i}_{\rm{FG}}+\sum_{l=e,\mu}\varepsilon^{l}_{\rm{FG}}$
(4)
$\displaystyle+\frac{1}{2}m^{2}_{\sigma}{\sigma}^{2}+\frac{1}{3}{g_{2}}{\sigma}^{3}+\frac{1}{4}{g_{3}}{\sigma}^{4}+\frac{1}{2}m^{2}_{\omega}{\omega}^{2}$
$\displaystyle+\frac{3}{4}{c_{3}}{\omega}^{4}+\frac{1}{2}m^{2}_{\rho}{\rho}^{2}+3{\Lambda}_{\rm{v}}\left(g^{2}_{\omega}{\omega}^{2}\right)\left(g^{2}_{\rho}{\rho}^{2}\right),$
$\displaystyle P_{\rm{HP}}$ $\displaystyle=$
$\displaystyle\sum_{i=p,n}P^{i}_{\rm{FG}}+\sum_{l=e,\mu}P^{l}_{\rm{FG}}$ (5)
$\displaystyle-\frac{1}{2}m^{2}_{\sigma}{\sigma}^{2}-\frac{1}{3}{g_{2}}{\sigma}^{3}-\frac{1}{4}{g_{3}}{\sigma}^{4}+\frac{1}{2}m^{2}_{\omega}{\omega}^{2}$
$\displaystyle+\frac{1}{4}{c_{3}}{\omega}^{4}+\frac{1}{2}m^{2}_{\rho}{\rho}^{2}+\Lambda_{\rm{v}}\left(g^{2}_{\omega}{\omega}^{2}\right)\left(g^{2}_{\rho}{\rho}^{2}\right),$

where $\varepsilon^{i}_{\rm{FG}}$ and $P^{i}_{\rm{FG}}$ denote the Fermi gas
contributions of species $i$ with a mass $m_{i}$ and degeneracy $N_{i}$,

$\displaystyle\varepsilon^{i}_{\rm{FG}}$ $\displaystyle=$ $\displaystyle
N_{i}\int_{0}^{k^{i}_{F}}\frac{d^{3}k}{(2\pi)^{3}}\sqrt{k^{2}+m_{i}^{2}},$ (6)
$\displaystyle P^{i}_{\rm{FG}}$ $\displaystyle=$
$\displaystyle\frac{N_{i}}{3}\int_{0}^{k^{i}_{F}}\frac{d^{3}k}{(2\pi)^{3}}\frac{k^{2}}{\sqrt{k^{2}+m_{i}^{2}}}.$
(7)

For nucleons, the effective mass $m_{i}=M^{\ast}$ is used and $N_{i}=2$
represents the spin degeneracy.

In the present calculations, we employ the BigApple parameterization given in
Fattoyev et al. (2020), which could provide an accurate description of ground-
state properties for finite nuclei across the nuclear chart. Its predictions
for infinite nuclear matter at the saturation density
$n_{0}=0.155\,\rm{fm}^{-3}$ are $E_{0}=-16.34\,\rm{MeV}$ for energy per
nucleon, $K=227\,\rm{MeV}$ for incompressibility,
$E_{\rm{sym}}=31.3\,\rm{MeV}$ for symmetry energy, and $L=39.8\,\rm{MeV}$ for
the slope of symmetry energy. The small slope parameter $L$ in the BigApple
model leads to acceptable radius and tidal deformability for a canonical
$1.4M_{\odot}$ neutron star, as compared to the estimations from astrophysical
observations. Moreover, the BigApple model predicts a maximum neutron-star
mass of $2.6M_{\odot}$, which is sufficiently large for exploring possible
phase transitions in neutron stars.

###  2.2 Quark Phase

We use a modified MIT bag model with vector interactions (vMIT) to describe
the quark phase. The inclusion of repulsive vector interactions among quarks
plays a crucial role in obtaining a stiff high-density EOS that is required by
the observations of $\sim 2M_{\odot}$ neutron stars. In the vMIT model, the
vector interaction is introduced via the exchange of a vector meson with the
mass $m_{V}$. We consider the quark matter consisting of three flavor quarks
($u$, $d$, and $s$) and leptons ($e$ and $\mu$). The Lagrangian density of the
vMIT model in the mean-field approximation is written as

$\displaystyle\mathcal{L}$ $\displaystyle=$
$\displaystyle\sum_{i=u,d,s}\left[\bar{\psi}_{i}\left(i\gamma_{\mu}\partial^{\mu}-m_{i}-g_{V}\gamma_{\mu}V^{\mu}\right)\psi_{i}-B\right.$
(8)
$\displaystyle\left.+\frac{1}{2}m_{V}^{2}V_{\mu}V^{\mu}\right]\Theta+\sum_{l=e,\mu}\bar{\psi_{l}}\left(i\gamma_{\mu}\partial^{\mu}-m_{l}\right)\psi_{l},$

where $B$ denotes the bag constant and $\Theta$ is the Heaviside step function
representing the confinement of quarks inside the bag. The nonzero mean field
$V_{0}$ is calculated from the equation of motion for the vector meson,

$\displaystyle m_{V}^{2}V_{0}$ $\displaystyle=$ $\displaystyle
g_{V}\sum_{i=u,d,s}n_{i},$ (9)

with $n_{i}$ being the number density of the quark flavor $i$. The quark
chemical potential is then given by

$\displaystyle\mu_{i}=\sqrt{{k^{i}_{F}}^{2}+m_{i}^{2}}+g_{V}V_{0},$ (10)

which is clearly enhanced by the vector potential.

The total energy density and pressure in quark matter are written as

$\displaystyle\varepsilon_{\rm{QP}}$ $\displaystyle=$
$\displaystyle\sum_{i=u,d,s}\varepsilon^{i}_{\rm{FG}}+\sum_{l=e,\mu}\varepsilon^{l}_{\rm{FG}}$
(11)
$\displaystyle+\frac{1}{2}\left(\frac{g_{V}}{m_{V}}\right)^{2}\left(n_{u}+n_{d}+n_{s}\right)^{2}+B,$
$\displaystyle P_{\rm{QP}}$ $\displaystyle=$
$\displaystyle\sum_{i=u,d,s}P^{i}_{\rm{FG}}+\sum_{l=e,\mu}P^{l}_{\rm{FG}}$
(12)
$\displaystyle+\frac{1}{2}\left(\frac{g_{V}}{m_{V}}\right)^{2}\left(n_{u}+n_{d}+n_{s}\right)^{2}-B,$

where $\varepsilon^{i}_{\rm{FG}}$ and $P^{i}_{\rm{FG}}$ denote the Fermi gas
contributions of species $i$, as given by Equations (6) and (7). For the quark
flavor $i$, the degeneracy $N_{i}=6$ arises from the spin and color degrees of
freedom, while $m_{i}$ represents the current quark mass. The vector
interactions among quarks are crucial for high-density EOS. In practice, we
vary the vector coupling $G_{V}=\left(g_{V}/m_{V}\right)^{2}$ in the range of
0–0.3 fm2 in order to examine the impact of vector interactions. We adopt the
current quark masses $m_{u}=m_{d}=5.5$ MeV and $m_{s}=95$ MeV in our
calculations. As for the bag constant, we mainly use the value of
$B^{1/4}=180$ MeV. It is well known that the bag constant could significantly
affect the EOS of quark matter and consequently influence the hadron-quark
phase transition. We will compare the results with different choices of $B$ in
the vMIT model for quark matter, to evaluate the influence of the bag constant
on the hadron-quark pasta phases in neutron stars.

##  3 Hadron-quark pasta phases

To describe the hadron-quark pasta phases, we employ the (EM) method within
the Wigner–Seitz approximation, where the whole space is divided into
equivalent cells with a geometric symmetry. The coexisting hadronic and quark
phases in a charge-neutral cell are assumed to be separated by a sharp
interface, while the particle densities in each phase are taken to be uniform
for simplicity. This is analogous to the compressible liquid-drop model used
in the study of nuclear liquid-gas phase transition at subnuclear densities
(Lattimer & Swesty, 1991; Bao et al., 2014). It is well known that the surface
tension at the interface plays a crucial role in determining the structure of
the mixed phase. The Gibbs and Maxwell constructions, respectively, correspond
to the two extreme cases of zero and large surface tension. The hadron-quark
mixed phase described in the EM method can be understood as an intermediate
state lying between the Gibbs and Maxwell constructions. For the EM method
used in our previous works (Wu & Shen, 2019; Ju et al., 2021), the electrons
are assumed to be uniformly distributed throughout the whole cell. This is
consistent with the GC, but contradicts the MC where the requirement of local
charge neutrality leads to different electron densities in the two phases. In
order to understand the transition from the GC to the MC with increasing
surface tension, we extend the EM method by allowing different electron
densities in the hadronic and quark phases in the present study. This may be
caused by the difference in the electric potential between the two phases,
considered as a result of the charge screening effect.

In the EM method, the total energy density of the mixed phase is expressed as

$\displaystyle\varepsilon_{\rm{MP}}$ $\displaystyle=$
$\displaystyle\chi\varepsilon_{\rm{QP}}+\left(1-\chi\right)\varepsilon_{\rm{HP}}+\varepsilon_{\rm{surf}}+\varepsilon_{\rm{Coul}},$
(13)

where $\chi=V_{\rm{QP}}/(V_{\rm{QP}}+V_{\rm{HP}})$ denotes the volume fraction
of the quark phase. The energy densities, $\varepsilon_{\rm{HP}}$ and
$\varepsilon_{\rm{QP}}$, are given by Equations (4) and (11), respectively.
The first two terms of Equation (13) represent the bulk contributions, while
the last two terms come from the finite-size effects. The surface and Coulomb
energy densities are calculated from

$\displaystyle{\varepsilon}_{\rm{surf}}$ $\displaystyle=$
$\displaystyle\frac{D\sigma\chi_{\rm{in}}}{r_{D}},$ (14)
$\displaystyle{\varepsilon}_{\rm{Coul}}$ $\displaystyle=$
$\displaystyle\frac{e^{2}}{2}\left(\delta
n_{c}\right)^{2}r_{D}^{2}\chi_{\rm{in}}\Phi\left(\chi_{\rm{in}}\right),$ (15)

with

$\displaystyle\Phi\left(\chi_{\rm{in}}\right)=\left\\{\begin{array}[]{ll}\frac{1}{D+2}\left(\frac{2-D\chi_{\rm{in}}^{1-2/D}}{D-2}+\chi_{\rm{in}}\right),&D=1,3,\\\
\frac{\chi_{\rm{in}}-1-\ln{\chi_{\rm{in}}}}{D+2},&D=2,\\\ \end{array}\right.$
(18)

where $D=1,2,3$ denotes the geometric dimension of the cell, and $r_{D}$
represents the size of the inner phase. $\chi_{\rm{in}}$ is the volume
fraction of the inner phase, i.e., $\chi_{\rm{in}}=\chi$ for droplet, rod, and
slab configurations, and $\chi_{\rm{in}}=1-\chi$ for tube and bubble
configurations. The charge density difference between the hadronic and quark
phases is defined as

$\displaystyle\delta n_{c}=n_{c}^{\rm{HP}}-n_{c}^{\rm{QP}},$ (19)

with

$\displaystyle n_{c}^{\rm{HP}}$ $\displaystyle=$ $\displaystyle
n_{p}-n_{e}^{\rm{HP}}-n_{\mu}^{\rm{HP}},$ (20) $\displaystyle n_{c}^{\rm{QP}}$
$\displaystyle=$
$\displaystyle\frac{2}{3}n_{u}-\frac{1}{3}n_{d}-\frac{1}{3}n_{s}-n_{e}^{\rm{QP}}-n_{\mu}^{\rm{QP}}.$
(21)

The surface tension $\sigma$ can significantly affect the structure of the
mixed phase (Endo et al., 2006; Yasutake et al., 2014; Wu & Shen, 2019). At
present, the value of $\sigma$ is poorly known, so it is usually taken as a
free parameter. In this study, we use the value of $\sigma=40$ MeV fm-2; which
is close to the prediction of the MIT bag model using the multiple reflection
expansion method (Ju et al., 2021).

The energy density of the mixed phase, given in Equation (13), is calculated
as a function of the following variables: $n_{p}$, $n_{n}$, $n_{u}$, $n_{d}$,
$n_{s}$, $n_{e}^{\rm{HP}}$, $n_{\mu}^{\rm{HP}}$, $n_{e}^{\rm{QP}}$,
$n_{\mu}^{\rm{QP}}$, $\chi$, and $r_{D}$. The values of these variables are
determined by solving a set of equilibrium equations between the hadronic and
quark phases at a given baryon density $n_{b}$. In the EM method, the
equilibrium conditions for coexisting two phases in the cell are derived by
minimizing the total energy density (13) under the constraints of baryon
number conservation and globe charge neutrality, which are written as

$\displaystyle\frac{\chi}{3}\left(n_{u}+n_{d}+n_{s}\right)+\left(1-\chi\right)\left(n_{p}+n_{n}\right)$
$\displaystyle=$ $\displaystyle n_{b},$ (22) $\displaystyle\chi
n_{c}^{\rm{QP}}+\left(1-\chi\right)n_{c}^{\rm{HP}}$ $\displaystyle=$
$\displaystyle 0.$ (23)

We introduce the Lagrange multipliers $\mu_{n}$ and $\mu_{e}$ for the
constraints, and then perform the minimization for the function

$\displaystyle w$ $\displaystyle=$
$\displaystyle\varepsilon_{\rm{MP}}-\mu_{n}\left[\frac{\chi}{3}\left(n_{u}+n_{d}+n_{s}\right)+\left(1-\chi\right)\left(n_{p}+n_{n}\right)\right]$
(24) $\displaystyle+\mu_{e}\left[\chi
n_{c}^{\rm{QP}}+\left(1-\chi\right)n_{c}^{\rm{HP}}\right].$

According to the definition of the chemical potential,
$\mu_{n}=\partial\varepsilon_{\rm{MP}}/\partial n_{n}^{\rm{MP}}$ and
$\mu_{e}=\partial\varepsilon_{\rm{MP}}/\partial n_{e}^{\rm{MP}}$ correspond to
the chemical potentials of neutrons and electrons in the mixed phase,
respectively. By minimizing $w$ with respect to the particle densities, we
obtain the following equilibrium conditions for chemical potentials:

$\displaystyle\mu_{e}$ $\displaystyle=$
$\displaystyle\mu_{e}^{\rm{HP}}-\frac{2\varepsilon_{\rm{Coul}}}{(1-\chi)\delta
n_{c}}=\mu_{e}^{\rm{QP}}+\frac{2\varepsilon_{\rm{Coul}}}{\chi\delta n_{c}},$
(25) $\displaystyle\mu_{p}$ $\displaystyle=$
$\displaystyle\mu_{n}-\mu_{e}^{\rm{HP}},$ (26) $\displaystyle\mu_{u}$
$\displaystyle=$
$\displaystyle\frac{1}{3}\mu_{n}-\frac{2}{3}\mu_{e}^{\rm{QP}},$ (27)
$\displaystyle\mu_{d}$ $\displaystyle=$
$\displaystyle\mu_{s}=\frac{1}{3}\mu_{n}+\frac{1}{3}\mu_{e}^{\rm{QP}}.$ (28)

It is necessary to clarify the definition of the chemical potential,
especially for charged particles. When the electric potential is taken into
account, the chemical potential of a charged particle is gauge dependent as
discussed in Tatsumi et al. (2003). In the mixed phase studied here, the
hadronic and quark phases are charged separately, which lead to nonzero
electric potentials in the two phases. We can understand the terms associated
with $\varepsilon_{\rm{Coul}}$ in Equation (25) as the contribution from the
electric potential. The chemical potentials of nucleons and quarks in the
above equations are defined by Equations (2) and (10), where the electric
potential is not taken into account. For the electrons, Equation (25) implies
that the total chemical potential including the contribution from the electric
potential remains constant throughout the cell, although the electron
densities in the hadronic and quark phases are different from each other.
Moreover, Equations (26)–(28) imply that local $\beta$ equilibrium should be
reached, which is also satisfied in the Gibbs and Maxwell constructions.

The equilibrium condition for the pressures is derived by minimizing $w$ with
respect to the volume fraction $\chi$, which is expressed as

$\displaystyle P_{\rm{HP}}$ $\displaystyle=$ $\displaystyle
P_{\rm{QP}}-\frac{2\varepsilon_{\rm{Coul}}}{\delta
n_{c}}\left[\frac{n_{c}^{\rm{QP}}}{\chi}+\frac{n_{c}^{\rm{HP}}}{1-\chi}\right]$
(29)
$\displaystyle\mp\frac{\varepsilon_{\rm{Coul}}}{\chi_{\rm{in}}}\left(3+\chi_{\rm{in}}\frac{{\Phi}^{{}^{\prime}}}{\Phi}\right),$

where the sign of the last term is “$-$” for droplet, rod, and slab
configurations, while it is “$+$” for tube and bubble configurations. The
pressure of the mixed phase is calculated from the thermodynamic relation,
$P_{\rm{MP}}=\mu_{n}n_{b}-\varepsilon_{\rm{MP}}$, which is actually equal to
$-w$ given in Equation (24). Comparing with the Gibbs condition of equal
pressures, the additional terms in Equation (29) are caused by the finite-size
effects. When the surface and Coulomb energies are neglected by taking the
limit $\sigma\rightarrow 0$, all the equilibrium equations derived above would
reduce to the Gibbs conditions.

By minimizing $w$ with respect to the size $r_{D}$, we obtain the well-known
relation ${\varepsilon}_{\rm{surf}}=2{\varepsilon}_{\rm{Coul}}$, which leads
to the formula for the size of the inner phase,

$\displaystyle r_{D}$ $\displaystyle=$
$\displaystyle\left[\frac{\sigma{D}}{e^{2}\left(\delta
n_{c}\right)^{2}\Phi}\right]^{1/3},$ (30)

whereas the size of the Wigner–Seitz cell is given by

$\displaystyle r_{C}$ $\displaystyle=$
$\displaystyle\chi_{\rm{in}}^{-1/D}r_{D}.$ (31)

In practice, we solve the equilibrium equations at a given baryon density
$n_{b}$ for all pasta configurations, and then determine the thermodynamically
stable state that has the lowest energy density. The hadron-quark mixed phase
exists only in the density range where its energy density is lower than that
of both hadronic matter and quark matter.

|
---|---
Figure 1: Pressure $P$ as a function of the neutron chemical potential $\mu_{n}$ for the BigApple model and the vMIT model with different vector coupling $G_{V}$ (a) and different bag constant $B$ (b). The filled circles indicate the crossing points of the hadronic and quark curves corresponding to the phase transition in the Maxwell construction. Figure 2: Energy densities of the mixed phase obtained using the EM method with $B^{1/4}=180$ MeV and $\sigma=40$ MeV fm-2 relative to those of the GC ($\sigma=0$). The filled circles indicate the transition points between different pasta phases. The results of the Maxwell construction (MC) are shown by green dotted lines. Figure 3: Particle density profiles in the Wigner–Seitz cell for a slab configuration at $n_{b}=0.6\,\rm{fm}^{-3}$. The results are obtained with $B^{1/4}=180$ MeV and $G_{V}=0.2\,\rm{fm}^{2}$. The hadron-quark interface and the cell boundary are indicated by the vertical lines. The results of the EM method (red thick lines) are compared with those of the CP method (blue thin lines) to show the charge screening effect. Figure 4: Size of the Wigner–Seitz cell ($r_{C}$) and that of the inner phase ($r_{D}$) as a function of the baryon density $n_{b}$ with $B^{1/4}=180$ MeV. The results of the EM and CP methods are shown by thick and thin lines, respectively. Figure 5: Charge densities of hadronic and quark phases, $n_{c}^{\rm{HP}}$ and $n_{c}^{\rm{QP}}$, as a function of the baryon density $n_{b}$ with $B^{1/4}=180$ MeV. The results of the EM and CP methods are shown by solid and dashed lines, respectively. Table 1: Onset Densities of the Hadron-Quark Pasta Phases and Pure Quark Matter Obtained Using Different Model Parameters and Methods. Model Parameter | Method | Onset Density (fm-3)
---|---|---
$B^{1/4}$ (MeV) |  $G_{V}$(fm2) |  | Droplet | Rod | Slab | Tube | Bubble | Quark
180 | 0.1 | EM | 0.378 | 0.460 | 0.514 | 0.615 | 0.651 | 0.688
180 | 0.2 | EM | 0.451 | 0.535 | 0.594 | 0.700 | 0.736 | 0.770
180 | 0.3 | EM | 0.544 | 0.640 | 0.710 | 0.832 | 0.872 | 0.907
170 | 0.2 | EM | 0.382 | 0.442 | 0.488 | 0.571 | 0.594 | 0.609
170 | 0.2 | CP | 0.414 | 0.419 | 0.484 | – | – | 0.580
180 | 0.2 | CP | 0.483 | 0.503 | 0.593 | – | – | 0.732
Table 2: Properties of Neutron Stars with the Maximum Mass $M_{\mathrm{max}}$.

Model Parameter | Method | $M_{\mathrm{max}}$ | $n_{b}^{c}$ | $R_{\mathrm{QP}}$ | $R_{\mathrm{MP}}$ | $R$
---|---|---|---|---|---|---
$B^{1/4}$ (MeV) |  $G_{V}$ $(\rm{fm}^{2})$ |  | $(M_{\odot})$ | $(\rm{fm}^{-3})$ | (km) | (km) | (km)
$180$ | 0.1 | EM | 1.93 | 0.68 | – | 5.23 | 13.20
$180$ | 0.2 | EM | 2.32 | 0.69 | – | 5.27 | 13.10
$180$ | 0.3 | EM | 2.52 | 0.73 | – | 4.67 | 12.83
$170$ | 0.2 | EM | 2.08 | 0.76 | 4.54 | 7.00 | 12.65
$170$ | 0.2 | CP | 2.15 | 0.71 | 3.17 | 5.88 | 12.98
$180$ | 0.2 | CP | 2.36 | 0.70 | – | 4.46 | 13.05

Note: The central baryon density is denoted by $n_{b}^{c}$, while
$R_{\mathrm{QP}}$, $R_{\mathrm{MP}}$, and $R$ correspond to the radii of pure
quark phase, mixed phase, and whole star, respectively.

##  4 Results and Discussion

In this section, we present numerical results for the hadron-quark pasta
phases, which are likely to appear in the core of massive neutron stars. We
employ the RMF model with the BigApple parameterization to describe the
hadronic matter. For the quark matter, we employ the MIT bag model with vector
interactions (vMIT), while the effects of the vector coupling $G_{V}$ and the
bag constant $B$ are discussed. The hadron-quark mixed phases are computed by
the EM method described in the previous section, and the properties of massive
neutron stars are calculated by using the EOS with quarks.

###  4.1 Hadron-quark mixed phases

In neutron-star matter, with increasing density the structured hadron-quark
mixed phase is expected to appear, which depends on the models used to
describe the hadronic and quark phases. In the present study, we employ the
RMF model with the BigApple parameterization for the hadronic phase, while the
quark phase is described by the vMIT model. To analyze the influence of model
parameters on the deconfinement phase transition, it is convenient to use the
Maxwell construction in which the phase transition appears at the crossing of
the hadronic EOS with the quark EOS in the pressure and chemical potential
plane. This is due to that two coexisting phases in the Maxwell construction
must have the same pressure and baryon chemical potential. We show in Figure 1
the pressure $P$ as a function of the neutron chemical potential $\mu_{n}$ for
the BigApple model and the vMIT model with different vector coupling $G_{V}$
(left panel) and different bag constant $B$ (right panel). According to the
Maxwell equilibrium conditions, the phase transition takes place at the
crossing of the hadronic and quark EOS curves. One can see from the left panel
that a larger $G_{V}$ in the vMIT model corresponds to a higher transition
pressure, which implies that the deconfinement phase transition is delayed
accordingly. In the right panel, we see that the transition pressure increases
as the bag constant $B$ increases. Therefore, it is reasonable to expect that
the formation of hadron-quark pasta phases may be delayed with increasing
vector coupling $G_{V}$ and bag constant $B$ in the vMIT model.

We compute the properties of a structured hadron-quark mixed phase in the EM
method, where the surface and Coulomb energies are included in the
minimization procedure. Compared to the Gibbs calculation without finite-size
effects, we obtain a higher energy for the pasta phase, since the surface and
Coulomb terms are always positive. In Figure 2, we present the energy
densities of pasta phases obtained in the EM method with $B^{1/4}=180$ MeV and
$\sigma=40$ MeV fm-2 relative to those of the GC (i.e., $\sigma=0$). The
filled circles represent the transition points between different pasta phases.
It is shown that the pasta configuration changes from droplet to rod, slab,
tube, and bubble with increasing baryon density $n_{b}$. For comparison, the
energy densities of pure hadronic matter and pure quark matter are displayed
by black dotted-dashed and dashed lines, respectively, whereas the results of
the Maxwell construction are shown by green dotted lines. It is seen that the
energy density of the Maxwell construction is higher than that of the pasta
phase, which implies that the structured mixed phase is more stable than the
Maxwell construction in our present calculations. However, for large enough
surface tension $\sigma$, the Maxwell construction may have lower energy than
the pasta phase, which has been discussed in the literature (Maslov et al.,
2019). To examine the influence of the quark vector interactions, we present
the results with $G_{V}=0.1$, 0.2, and 0.3 fm2 in the left, middle, and right
panels, respectively. One can see that as $G_{V}$ increases, the mixed phase
is shifted toward higher densities with a wider range. This tread is in
agreement with the behavior observed in the Maxwell construction (see the left
panel of Figure 1). In Table 1, we present the onset densities of the hadron-
quark pasta phases and pure quark matter obtained using different model
parameters and methods. From this table, one can see the effect of $G_{V}$ and
the difference between the EM and CP methods.

It is interesting to study the geometric structure of the hadron-quark mixed
phase, which may be attributed to the competition between surface and Coulomb
energies. For checking the influence of the method used, we compare the
results using the EM method developed in the present work with those obtained
in the simple CP method (Wu & Shen, 2019). We emphasize that in the CP method
the Gibbs conditions for phase equilibrium are enforced, while the surface and
Coulomb energies are taken into account perturbatively, and therefore the
charge screening effect is disregarded in this case. In contrast, the EM
method incorporates the surface and Coulomb contributions self-consistently in
deriving the equilibrium conditions, which leads to the rearrangement of
charged particles known as the charge screening effect. In order to see how
large the charge screening effect is, we display in Figure 3 the particle
density profiles in the Wigner–Seitz cell for a slab configuration at
$n_{b}=0.6\,\rm{fm}^{-3}$. The calculations are carried out with $B^{1/4}=180$
MeV and $G_{V}=0.2\,\rm{fm}^{2}$. The hadron-quark interface and the cell
boundary are indicated by the vertical lines. One can see that the structure
size of the EM method is larger than that of the CP method. This is because
the charge screening effect tends to reduce the net charge density in each
phase in order to lower the Coulomb energy. In the negatively charged quark
matter and positively charged hadronic matter, the particle densities $n_{d}$,
$n_{s}$, and $n_{p}$ are reduced in the EM method, whereas $n_{u}$ is enhanced
in comparison to that of the CP method. The electron densities in the quark
and hadronic matter are, respectively, reduced and enhanced due to the same
reason. It is expected that as the surface tension $\sigma$ increases, the
charge screening effect becomes more pronounced, and finally the local charge
neutrality is reached in the Maxwell construction.

In Figure 4, the size of the Wigner–Seitz cell ($r_{C}$) and that of the inner
phase ($r_{D}$) are plotted as a function of the baryon density $n_{b}$, where
the vMIT model with $B^{1/4}=180$ MeV are adopted for quark matter. The
results with $G_{V}=0$ and 0.2 fm2 are shown in the left and right panels,
respectively. We compare the results obtained from the EM method (thick lines)
with those from the CP method (thin lines) in order to check the charge
screening effect. It is found that both $r_{D}$ and $r_{C}$ of the EM method
are larger than that of the CP method, which can be explained by the
rearrangement of charged particles as shown in Figure 3. According to Equation
(30), a small charge density difference $\delta n_{c}$ leads to a large pasta
size $r_{D}$. The value $\delta n_{c}=n_{c}^{\rm{HP}}-n_{c}^{\rm{QP}}$ can be
observed in Figure 5, where the charge densities $n_{c}^{\rm{HP}}$ and
$n_{c}^{\rm{QP}}$ are shown as a function of the baryon density $n_{b}$. One
can see that the magnitudes of $n_{c}^{\rm{HP}}$ and $n_{c}^{\rm{QP}}$ in the
EM method are significantly reduced due to the charge screening effect, and
hence their difference $\delta n_{c}$ is also reduced in comparison to the CP
results. The reduction of $\delta n_{c}$ caused by the charge rearrangement in
the EM method leads to the increase of $r_{D}$ and $r_{C}$ in Figures 3 and 4.
In addition, a comparison between the EM and CP methods can be found in Table
1, where the onset and configuration of pasta phases are somewhat different.
As shown in the last two lines of Table 1, due to relatively large Coulomb
energies of the CP method, the tube and bubble configurations are
energetically disfavored and will not appear before the transition to pure
quark matter. A similar analysis is available for understanding the impact of
the vector coupling $G_{V}$ on the pasta size. Compared to the case of
$G_{V}=0$ in the left panel of Figure 4, the values of $r_{D}$ and $r_{C}$
obtained with $G_{V}=0.2\,\rm{fm}^{2}$ in the right panel are relatively
small. This can be understood from the behavior of $\delta n_{c}$ in Figure 5,
where the values of $\delta n_{c}=n_{c}^{\rm{HP}}-n_{c}^{\rm{QP}}$ with
$G_{V}=0.2\,\rm{fm}^{2}$ are clearly larger than those with $G_{V}=0$, which
lead to smaller $r_{D}$ and $r_{C}$ in the right panel of Figure 4.
Furthermore, we can see that as the density increases, $r_{D}$ in the droplet,
rod, and slab configurations increases, whereas it decreases in the tube and
bubble phases. These trends are related to a monotonic increase of the quark
volume fraction $\chi$ during the phase transition.

Figure 6: Pressure $P$ as a function of the baryon density $n_{b}$ for hadronic, mixed, and quark phases. The results of pasta phases obtained using the EM method with $B^{1/4}=180$ MeV are compared to those of the GC and MC. Figure 7: Particle fraction $Y_{i}$ as a function of the baryon density $n_{b}$ for hadronic, mixed, and quark phases. The shaded area denotes the mixed phase region, where the results of pasta phases are obtained using the EM method with $B^{1/4}=180$ MeV and $G_{V}=0.2\,\rm{fm}^{2}$. |
---|---
Figure 8: Mass–radius relations of neutron stars for different model
parameters. The results of the pure hadronic EOS are compared to those
including the hadron-quark mixed phase described in the EM method for
different vector coupling $G_{V}$ (a) and different bag constant $B$ (b). The
filled circles indicate the onset of the star containing hadron-quark pasta
phases, while the filled squares denote the star with the maximum mass. The
red horizontal line with arrows at both ends represents the constraint on
$R_{1.4}$ inferred from GW170817 (Abbott et al., 2018). The green and purple
shaded areas correspond to simultaneous measurements of the mass and radius
from NICER for PSR J0030+0451 (Miller et al., 2019) and PSR J0740+6620 (Miller
et al., 2021), respectively. The mass constraint from GW190814 (Abbott et al.,
2020b) is depicted by the pink horizontal bar. Figure 9: Dimensionless tidal
deformability $\Lambda$ as a function of the neutron-star mass $M$. The
results of the pure hadronic EOS are compared to those including the hadron-
quark mixed phase described in the EM method with $B^{1/4}=180$ MeV and
$G_{V}=0.2$ fm2. The constraint $70\leq\Lambda_{1.4}\leq 580$ is inferred from
the analysis of GW170817 (Abbott et al., 2018), while
$458\leq\Lambda_{1.4}\leq 889$ is achieved based on the assumption that
GW190814’s secondary is a neutron star (Abbott et al., 2020b).

###  4.2 Properties of Neutron Stars

In Figure 6, we show the pressure $P$ as a function of the baryon density
$n_{b}$ for hadronic, mixed, and quark phases. The calculations of hadron-
quark pasta phases are performed in the EM method, where the hadronic matter
is described by the BigApple model and the quark matter by the vMIT model with
$B^{1/4}=180$ MeV. For comparison, the results of the Gibbs and Maxwell
constructions are displayed by red solid and green dotted lines, respectively.
The pressures with the Maxwell construction remain constant during the phase
transition, whereas the pressures with the GC increase with $n_{b}$ over a
broad range. It is shown that the results of pasta phases lie between the
Gibbs and Maxwell constructions. As the vector coupling $G_{V}$ increases from
left to right panels, we see that the hadron-quark mixed phases appear at
higher densities and pressures. Especially, the ending of the mixed phase
shows a more significant $G_{V}$ dependence than the beginning, which may be
attributed to the increasing quark fraction during the phase transition.

In Figure 7, we display the particle fraction $Y_{i}=n_{i}/n_{b}$ as a
function of the baryon density $n_{b}$. The shaded area denotes the mixed
phase region, where the results of hadron-quark pasta phases are obtained in
the EM method. We employ the BigApple model for hadronic matter and the vMIT
model with $B^{1/4}=180$ MeV and $G_{V}=0.2\,\rm{fm}^{2}$ for quark matter. At
low densities, the matter consists of neutrons, protons, and electrons,
whereas the muons appear at about $0.11\,\rm{fm}^{-3}$. When the deconfined
quarks are present in the mixed phase, the quark fractions $Y_{u}$, $Y_{d}$,
and $Y_{s}$ increase rapidly together with a decrease of neutron fraction
$Y_{n}$. Meanwhile, $Y_{e}$ and $Y_{\mu}$ decrease significantly because the
quark matter is negatively charged that takes the role of electrons for
satisfying the constraint of global charge neutrality. On the other hand, the
hadronic matter in the mixed phase is positively charged, so the proton
fraction $Y_{p}$ increases slightly at the beginning of the mixed phase. As
the baryon density $n_{b}$ increases, $Y_{p}$ and $Y_{n}$ decrease to very low
values due to the increase of the quark volume fraction $\chi$ in the phase
transition. At sufficiently high densities, the matter turns into a pure quark
phase, where $Y_{u}\approx Y_{d}\approx Y_{s}\approx 1/3$ is nearly achieved
under the conditions of charge neutrality and chemical equilibrium. The direct
Urca (dUrca) process, known as the most efficient mechanism for neutron-star
cooling, is closely related to the particle fractions shown in Figure 7. In
the hadronic matter, the dUrca process, i.e., the electron capture by a proton
and the $\beta$-decay of a neutron, is mainly determined by the proton
fraction $Y_{p}$, which must exceed a critical value to allow simultaneous
energy and momentum conservation (Fantina et al., 2013). In the case of simple
$npe$ matter, the dUrca process can occur only for $Y_{p}\geq 1/9$. When muons
are present under the equilibrium condition $\mu_{e}=\mu_{\mu}$, the critical
$Y_{p}$ for the dUrca process is in the range of $11.1\%-14.8\%$ (Lattimer et
al., 1991). According to observations of thermal radiation from neutron stars,
the dUrca process is unlikely to occur in neutron stars with masses below
$1.5M_{\odot}$, since it would lead to an unacceptably fast cooling in
disagreement with the observations (Fantina et al., 2013). In our present
calculations, a $1.5M_{\odot}$ neutron star remains in the pure hadronic phase
with a central density of about $0.34\,\rm{fm}^{-3}$, where the proton
fraction is lower than the critical value for the dUrca process. We note that
the BigApple model predicts a rather large threshold density of
$0.64\,\rm{fm}^{-3}$ for the dUrca process, which is related to its small
symmetry energy slope $L$ (see the discussion in Ji et al. 2019). This result
is compatible with the cooling observations of neutron stars.

The properties of static neutron stars are obtained by solving the well-known
Tolman–Oppenheimer–Volkoff (TOV) equation with the EOS over a wide range of
densities. In the present calculations, we use the Baym–Pethick–Sutherland EOS
(Baym et al., 1971) for the outer crust below the neutron drip density, while
the inner crust EOS is based on a Thomas–Fermi calculation using the TM1e
parameterization of the RMF model (Shen et al., 2020). The crust EOS is
matched to the EOS of uniform neutron-star matter at the crossing of the two
segments. At very high densities, the quark degrees of freedom are taken into
account using the EM method for the hadron-quark pasta phases. In Figure 8, we
display the predicted mass-radius relations of neutron stars, together with
several constraints from astrophysical observations. Compared to the results
using pure hadronic EOS (solid lines), the inclusion of quarks leads to an
obvious reduction of the maximum neutron-star mass $M_{\mathrm{max}}$. As
shown in the left panel, the results of massive stars are strongly dependent
on the vector coupling $G_{V}$. It is found that both the onset of hadron-
quark pasta phases (filled circles) and the values of $M_{\mathrm{max}}$
(filled squares) increase with increasing $G_{V}$.

For a canonical $1.4M_{\odot}$ neutron star, only the pure hadronic phase is
present. The prediction of the BigApple model for the radius $R_{1.4}$ is
compatible with the constraints inferred from the GW170817 event and NICER for
PSR J0030+0451.

In order to examine the impact of the bag constant $B$, we compare in the
right panel the results between $B^{1/4}=180$ MeV (blue dashed line) and
$B^{1/4}=170$ MeV (green dashed line). The reduction of $M_{\mathrm{max}}$
with $B^{1/4}=170$ MeV is more pronounced than with $B^{1/4}=180$ MeV. This is
because the formation of hadron-quark pasta phases appears at lower densities
for $B^{1/4}=170$ MeV (see Table 1), which leads to a stronger softening of
the EOS. To analyze the effects of the hadron-quark phase transition in more
detail, we show in Table 2 the resulting properties of neutron stars with the
maximum mass. It is found that both $G_{V}$ and $B$ can significantly affect
the structure and maximum mass of neutron stars, whereas the influence of the
charge screening (i.e., the difference between the EM and CP methods) is
relatively small. In most cases, a structured hadron-quark mixed phase can be
formed in massive neutron stars inside the radius $R_{\mathrm{MP}}$. The pure
quark matter appears in the interior of neutron stars only for the cases with
small $B$ where the central density $n_{b}^{c}$ is larger than the onset of
the pure quark phase. We find that for all cases listed in Table 2, the
central density of a canonical $1.4M_{\odot}$ neutron star is not high enough
to form hadron-quark pasta phases, so it remains in a pure hadronic phase.
These results are qualitatively consistent with the arguments of Annala et al.
(2020), which suggested that the matter in the core of a $1.4M_{\odot}$
neutron star is compatible with nuclear model calculations, whereas the matter
in the interior of a $2M_{\odot}$ neutron star exhibits characteristics of the
deconfined quark phase. In Figure 9, we plot the dimensionless tidal
deformability $\Lambda$ as a function of the neutron-star mass $M$. The
results of the pure hadronic EOS are compared to those including the hadron-
quark mixed phase described in the EM method with $B^{1/4}=180$ MeV and
$G_{V}=0.2$ fm2. It is shown that the inclusion of quarks leads to small
differences in $\Lambda$ for massive neutron stars, which is related to the
reduction of the radius as shown in Figure 8. Considering the constraints on
$\Lambda_{1.4}$ inferred from the analysis of GW170817 and GW190814, the
BigApple model provides acceptable results for the tidal deformability of
neutron stars. The effects of the hadron-quark phase transition can be
observed only for very massive stars as shown in the inset of Figure 9.

##  5 Summary

Motivated by the recent advances in astrophysical observations, we studied the
properties of the hadron-quark pasta phases, which may appear in the interior
of massive neutron stars. The structured mixed phase is described within the
Wigner–Seitz approximation, where the whole space is divided into equivalent
cells with a geometric symmetry. The coexisting hadronic and quark phases
inside the cell are assumed to be separated by a sharp interface with constant
densities in each phase. We extended the EM method for describing the hadron-
quark pasta phases by allowing different electron densities in the two
coexisting phases, which is helpful for understanding the transition from the
GC to the MC. In the EM method, the surface and Coulomb energies are included
in the minimization procedure, which results in different equilibrium
equations from the Gibbs conditions. Comparing to the simple CP method, the
treatment of surface and Coulomb energies in the EM method leads to the
rearrangement of charged particles, known as the charge screening effect,
which can significantly affect the structure of the hadron-quark mixed phase.
The resulting EOS obtained from the EM method was found to lie between those
of the GC and MC.

In the present study, we have employed the RMF model with the BigApple
parameterization to describe the hadronic matter, while the quark matter is
described by a modified MIT bag model with vector interactions (vMIT). The
BigApple model could provide a good description of finite nuclei across the
nuclear chart, while its prediction for the maximum neutron-star mass is as
large as $2.6M_{\odot}$. In addition, the BigApple model predicts acceptable
radius and tidal deformability for a canonical $1.4M_{\odot}$ neutron star, as
compared to the estimations from astrophysical observations. For the quark
matter in the vMIT model, the vector interactions among quarks could
significantly stiffen the EOS at high densities and help to enhance the
maximum mass of neutron stars. We found that as the vector coupling $G_{V}$
increases, the hadron-quark pasta phases appear at higher densities, and
especially the ending of the mixed phase shows a more significant $G_{V}$
dependence than the beginning. In addition, a larger $G_{V}$ leads to a
smaller pasta size, which is caused by relatively large charge density.
Meanwhile, the pasta size is also affected by the charge screening effect in
the EM method, where the rearrangement of charged particles can lower the
Coulomb energy and enhance the pasta size.

We investigated the properties of massive neutron stars containing the hadron-
quark pasta phases. It was found that the inclusion of quarks could
considerably soften the EOS and reduce the maximum mass of neutron stars. The
results of massive stars clearly depend on the vector coupling $G_{V}$ and the
bag constant $B$ used. A large $G_{V}$ corresponds to a weak reduction of
$M_{\mathrm{max}}$ relative to that of pure hadronic stars. In most cases
considered in this study, a structured hadron-quark mixed phase can be formed
in the interior of massive neutron stars, whereas a canonical $1.4M_{\odot}$
neutron star remains in a pure hadronic phase. Our results showed that the
presence of quarks inside neutron stars could be compatible with current
constraints inferred from NICER data and gravitational-wave observations.

This work was supported in part by the National Natural Science Foundation of
China (grants Nos. 12175109 and 11775119).

## References

  * Abbott et al. (2017) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017, Phys. Rev. Lett., 119, 161101
  * Abbott et al. (2018) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2018, Phys. Rev. Lett., 121, 161101
  * Abbott et al. (2020a) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2020, ApJ, 892, L3
  * Abbott et al. (2020b) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2020, ApJ, 896, L44
  * Alam et al. (2016) Alam, N., Agrawal, B. K., Fortin, M., Pais, H., Providência, C., Raduta, Ad. R., & Sulaksono, A. 2016, Phys. Rev. C, 94, 052801(R)
  * Annala et al. (2020) Annala, E., Gorda, T., Kurkela, A., Nättilä, J., & Vuorinen, A. 2020, NatPh, 16, 907
  * Antoniadis et al. (2013) Antoniadis, J., Freire, P. C. C., Wex, N., et al. 2013, Sci, 340, 448
  * Arzoumanian et al. (2018) Arzoumanian, Z., et al. 2018, ApJS, 235, 37
  * Bao et al. (2014) Bao, S. S., Hu, J. N., Zhang, Z. W., & Shen, H. 2014, Phys. Rev. C, 90, 045802
  * Baym et al. (1971) Baym, G., Bethe, H. A., & Pethick, C. J. 1971, Nucl. Phys. A, 175, 225
  * Baym et al. (2018) Baym, G., Hatsuda, T., Kojo, T., Powell, P. D., Song, Y., & Takatsuka, T. 2018, RPPh, 81, 056902
  * Bhattacharyya et al. (2010) Bhattacharyya, A., Mishustin, I. N., & Greiner, W. 2010, JPhG, 37, 025201
  * Cromartie et al. (2020) Cromartie, H. T., Fonseca, E., Ransom, S. M., et al. 2020, NaAs, 4, 72
  * Demircik et al. (2021) Demircik, T., Ecker, C., & Järvinen, M. 2021, ApJ, 907, L37
  * Demorest et al. (2010) Demorest, P. B., Pennucci, T., Ransom, S. M., Roberts, M. S. E., & Hessels, J. W. T. 2010, Natur, 467, 1081
  * Dexheimer et al. (2021) Dexheimer, V., Gomes, R., Klähn, T., Han, S., & Salinas, M. 2021, Phys. Rev. C, 103, 025808
  * Endo et al. (2006) Endo, T., Maruyama, T., Chiba, S., & Tatsumi, T. 2006, PThPh, 115, 337
  * Essick & Landry (2020) Essick, R., & Landry, P. 2020, ApJ, 904, 80
  * Fantina et al. (2013) Fantina, A. F., Chamel, N., Pearson, J. M., & Goriely, S. 2013, A&A, 559, A128
  * Fattoyev et al. (2020) Fattoyev, F. J., Horowitz, C. J., Piekarewicz, J., & Reed, B. 2020, Phys. Rev. C, 102, 065805
  * Fattoyev et al. (2010) Fattoyev, F. J., Horowitz, C. J., Piekarewicz, J., & Shen, G. 2010, Phys. Rev. C, 82, 055802
  * Fonseca et al. (2016) Fonseca, E., et al. 2016, ApJ, 832, 167
  * Fonseca et al. (2020) Fonseca, E., et al. 2021, ApJ, 915, L12
  * Glendenning (1992) Glendenning, N. K. 1992, Phys. Rev. D, 46, 1274
  * Gomes et al. (2019) Gomes, R. O., Char, P., & Schramm, S. 2019, ApJ, 877, 139
  * Han et al. (2019) Han, S., Mamun, M. A. A., Lalit, S., Constantinou, C., & Prakash, M. 2019, Phys. Rev. D, 100, 103022
  * Heiselberg et al. (1993) Heiselberg, H., Pethick, C. J., & Staubo, E. F. 1993, Phys. Rev. Lett., 70, 1355
  * Huang et al. (2020) Huang, K. X., Hu, J. N., Zhang, Y., & Shen, H. 2020, ApJ, 904, 39
  * Ji et al. (2019) Ji, F., Hu, J. N., Bao, S. S., & Shen, H. 2019, Phys. Rev. C, 100, 045801
  * Ju et al. (2021) Ju, M., Wu, X. H., Ji, F., Hu, J. N., & Shen, H. 2021, Phys. Rev. C, 103, 025809
  * Klähn & Fischer (2015) Klähn, T. & Fischer, T. 2015, ApJ, 810, 134
  * Lalazissis et al. (1997) Lalazissis, G. A., Köning, J., & Ring, P. 1997, Phys. Rev. C, 55, 540
  * Lattimer & Swesty (1991) Lattimer, J. M., & Swesty, F. D. 1991, Nucl. Phys. A, 535, 331
  * Lattimer et al. (1991) Lattimer, J. M., Pethick, C. J., Prakash, M., & Haensel, P. 1991, Phys. Rev. Lett., 66, 2701
  * Lattimer & Prakash (2016) Lattimer, J. M., & Prakash, M. 2016, PhR, 621, 127
  * Li et al. (2020) Li, J. J., Sedrakian, A., & Weber, F. 2020, PhLB, 810, 135812
  * Lonardoni et al. (2015) Lonardoni, D., Lovato, A., Gandolfi, S., & Pederiva, F. 2015, Phys. Rev. Lett., 114, 092301
  * Lopes et al. (2021) Lopes, L., L., Biesdorf, C., & Menezes, D. P. 2021, PhyS, 96, 065303
  * Maruyama et al. (2007) Maruyama, T., Chiba, S., Schulze, H.-J., & Tatsumi, T. 2007, Phys. Rev. D, 76, 123015
  * Maslov et al. (2019) Maslov, K., Yasutake, N., Blaschke, D., Ayriyan, A., Grigorian, H., Maruyama, T., Tatsumi, T., & Voskresensky, D. N. 2019, Phys. Rev. C, 100, 025802
  * Miller et al. (2019) Miller, M. C., Lamb, F. K., Dittmann, A. J., et al. 2019, ApJ, 887, L24
  * Miller et al. (2021) Miller, M. C., Lamb, F. K., Dittmann, A. J., et al. 2021, ApJ, 918, L28
  * Most et al. (2020) Most, E. R., Papenfort, L. J., Weih, L. R., & Rezzolla, L. 2020, MNRAS Lett, 499, L82
  * Oertel et al. (2017) Oertel, M., Hempel, M., Klähn, T., & Typel, S. 2017, RvMP, 89, 015007
  * Riley et al. (2019) Riley, T. E., Watts, A. L., Bogdanov, S., et al. 2019, ApJ, 887, L21
  * Riley et al. (2021) Riley, T. E., Watts, A. L., Ray, P. S., et al. 2021, ApJL, 918, L27
  * Schertler et al. (2000) Schertler, K., Greiner, C., Schaffner-Bielich, J., & Thoma, M. H. 2000, Nucl. Phys. A, 677, 463
  * Shen et al. (2020) Shen, H., Fan, J., Hu, J. N., & Sumiyoshi, K. 2020, ApJ, 891, 148
  * Spinella et al. (2016) Spinella, W. M., Weber, F., Contrera, G. A., & Orsaria, M. G. 2016, EPJA, 52, 61
  * Sugahara & Toki (1994) Sugahara, Y., & Toki, H. 1994, Nucl. Phys. A, 579, 557
  * Tan et al. (2020) Tan, H., Noronha-Hostler, J., & Yunes, N. 2020, Phys. Rev. Lett., 125, 261104
  * Tatsumi et al. (2003) Tatsumi, T., Yasuhira, M., & Voskresensky, D. 2003, Nucl. Phys. A, 718, 359
  * Tews et al. (2021) Tews, I., Pang, P. T. H., Dietrich, T., Coughlin, M. W., Antier, S., Bulla, M., Heinzel, J., & Issa, L. 2021, ApJ, 908, L1
  * Tsokaros et al. (2020) Tsokaros, A., Ruiz, M., & Shapiro, S. L. 2020, ApJ, 905, 48
  * Weber et al. (2019) Weber, F., Farrell, D., Spinella, W. M., Malfatti, G., Orsaria, M. G., Contrera, G. A., & Maloney, I. 2019, Univ, 5, 169
  * Wu & Shen (2017) Wu, X. H. & Shen, H. 2017, Phys. Rev. C, 96, 025802
  * Wu & Shen (2019) Wu, X. H. & Shen, H. 2019, Phys. Rev. C, 99, 065802
  * Xu et al. (2010) Xu, J., Chen, L. W., Ko, C. M., & Li, B. A. 2010, Phys. Rev. C, 81, 055803
  * Yang & Shen (2008) Yang, F. & Shen, H. 2008, Phys. Rev. C, 77, 025801
  * Yasutake et al. (2014) Yasutake, N., Łastowiecki, R., Benić, S., Blaschke, D., Maruyama, T., & Tatsumi, T. 2014, Phys. Rev. C, 89, 065803
  * Zhang & Li (2020) Zhang, N. B., & Li, B. A. 2020, ApJ, 902, 38