Phases of Hadron-Quark Matter in (Proto) Neutron Stars

In the first part of this paper, we investigate the possible existence of a structured hadron-quark mixed phase in the cores of neutron stars. This phase, referred to as the hadron-quark pasta phase, consists of spherical blob, rod, and slab rare phase geometries. Particular emphasis is given to modeling the size othis phase in rotating neutron stars. We use the relativistic mean-field theory to model hadronic matter and the non-local three-flavor Nambu-Jona-Lasinio model to describe quark matter. Based on these models, the hadron-quark pasta phase exists only in very massive neutron stars, whose rotational frequencies are less than around 300 Hz. All other stars are not dense enough to trigger quark deconfinement in their cores. Part two of the paper deals with the quark-hadron composition of hot (proto) neutron star matter. To this end we use a local three-flavor Polyakov-Nambu-Jona-Lasinio model which includes the 't Hooft (quark flavor mixing) term. It is found that this term leads to non-negligible changes in the particle composition of (proto) neutron stars made of hadron-quark matter.


Introduction
Neutron stars are extremely compact objects. They are typically about 20 to 30 km across and spin rapidly, often making many hundreds of rotations per second (pulsar PSR J1748-2446ad is the fastest spinning neutron star known to date, rotating at 716 Hz.) Depending on mass and rotational frequency, gravity compresses the matter in the core regions of neutron stars to densities that are several times higher than the density of atomic nuclei. At such enormous densities, atomic nuclei are squeezed so tightly together that new types of particles and novel states of matter may be created. These include the creation of hyperons and delta isobars and the formation of a mixed phase of hadrons and quarks deep in the central regions of neutron stars. As shown by Glendenning [1,2], the quark-hadro mixed phase region could segregate phases by net charge to minimize the total energy of the system, leading to the formation of a hadron-quark Coulomb lattice. Since this lattice is qualitatively similar to the hypothesized "nuclear pasta" structures in the crustal regions of neutron stars [1][2][3][4][5], it is referred to as hadron-quark pasta phase. 1 All these possibilities make neutron stars and the events creating them (e.g., GW 170817) superb astrophysical laboratories for a wide range of physical studies. With observational data accumulating rapidly from both orbiting (HST, ATHENA, NICER) and ground based (skA, FAST, LIGO, VIRGO) observatories there has never been a more exciting time than today to study neutron stars.

Hadronic Composition of Neutron Star Matter at Finite Temperature
A neutron star is born as a result of the gravitational collapse of the iron core of a massive evolved progenitor star (M ∼ 8-25 M ) in a type-II supernova explosion. The iron core of such a star collapses when its mass reaches the Chandrasekhar limit, which depends on the mass of the progenitor star. Due to the Fermi pressure of the nucleons, the collapse stops when nuclear matter density is reached and the core bounces back. Shortly after core bounce (some 10 ms) a hot, lepton rich neutron star, called a proto-neutron star, is formed [8][9][10][11][12][13]. This star consists of a shocked envelope with an entropy per baryon of s ∼ 4-10 (in units of the Boltzmann constant) and an unshocked core with s ∼ 1. The envelope and the core contain nearly the same mass of about 0.6-0.8 M . During the so-called Kelvin-Helmholtz cooling phase the lepton number decreases in the proton-neutron star due to the loss of neutrinos and the proto-neutron star evolves in about 10 to 30 s to a hot, lepton poor neutron star with an entropy per baryon of s ∼ 1-2. After several minutes, this hot neutron star cools to a cold neutron star with temperatures T < 1 MeV, which then continues to slowly cool via neutrino and photon emission until its thermal radiation becomes too weak to be detectable with X-ray telescopes, after about 10 7 years.
To study hot hadronic matter at conditions that exist in the cores of (proto) neutron stars, we use the relativistic nuclear mean-field (RMF) approximation and choose the GM1L, HV, and DD2 nuclear parametrizations listed in Table 1 that satisfy current constraints from studies of symmetric nuclear matter and from observations of neutron stars. The GM1L parametrization is based on the standard RMF model and was proposed in Refs. [1,14] but has been constructed to satisfy recent constraints from symmetric nuclear matter at saturation density while simultaneously being able to satisfy the 2.01 M mass constraint set by PSR J0348+0432 with the inclusion of hyperons in SU(3) flavor symmetry. The DD2 parametrization is also based on the standard RMF approach but utilizes a density dependent meson-baryon coupling scheme, with the parameters of the density dependence set such that the model reproduces the properties of finite nuclei [15][16][17][18]. The lagrangian of the density-dependent relativistic mean-field model is given by [17][18][19][20] where g σB (n), g ωB (n) and g ρB (n) are density dependent meson-baryon coupling constants and n = ∑ B=N,Y,∆ n B is the total baryon number density. In this work we consider all members of the spin-1/2 baryon octet, which include nucleons N = (n, p) and hyperons Y = (Λ, Σ − , Σ 0 , Σ + , Ξ − , Ξ 0 ), as well as the spin 3/2 delta isobar quartet (∆ − , ∆ 0 , ∆ + , ∆ ++ ). The density dependent coupling constants are given by [19] for i = σ, ω, and 1 This phase has nothing to do with the crystalline quark matter phases that have been introduced in the context of discussing color superconducting quark matter phases [6,7].
for the ρ meson. This approach is parametrized to reproduce the properties of symmetric nuclear matter at saturation density n 0 . The parametrizations used in this work (GM1L, HV, and DD2) are shown in Table 1. The properties listed in this table are the energy per nucleon (E 0 ), nuclear incompressibility (K 0 ), isospin asymmetry energy (J), and effective nucleon mass (m * /m N ). The isovector meson-baryon coupling constant is treated density-dependent. Its numerical value has been fit to the slope of the asymmetry energy (L 0 ) at n 0 . The parameters a i , b i , c i , d i , and a ρ in Equations (2) and (3) are fixed by the binding energies, charge and diffraction radii, spin-orbit splittings, and the neutron skin thickness of finite nuclei. The density dependence of the meson-baryon couplings of the DD2 parametrization eliminates the need for non-linear self-interactions of the σ-meson field. The non-linear terms [23,24] of the σ-meson field in Equation (1) therefore only contribute to the GM1L and HV models. The meson-hyperon coupling constants have been determined following the Nijmegen extended soft core (ESC08) model [25]. The relative isovector meson-hyperon coupling constants were scaled with the hyperon isospin. For the ∆-isobar we used x σ∆ = x ω∆ = 1.1 and x ρ∆ = 1.0, where x i∆ = g i∆ /g iN (for details, see Ref. [14]).
The meson mean-field equations following from Equation (1) are given by where I 3B is the 3-component of isospin and n s B and n B are the scalar and particle number densities of each baryon B, given by The quantity f B∓ denotes the Fermi-Dirac distribution function, and E * B stands for the effective single-baryon energy, 3 of 15 The quantity γ B = 2J B + 1 accounts for the spin degeneracy, and m * B = m B − g σB (n)σ denotes the effective baryon mass. The effective baryon chemical potential, µ * B , is defined as where R is the so-called rearrangement term given by which is important for achieving thermodynamic consistency. The total baryonic pressure of the matter follows from while the energy density of the system follows from the standard Gibbs relation (P, T, n) = −P + To complete the description of the matter in the cores of (proto) neutron stars, leptons (L) need to be included in the treatment. Leptons can be treated as free, relativistic Fermions whose total thermodynamic potential is given by (γ L = 2) with The sum in Equation (13) is over electrons and muons (with masses m L ) as well as over massless neutrinos, ν e . Neutrinos do not contribute to cold neutron stars matter, but are present in the cores of newly formed proto-neutron stars.

Electric Charge Neutrality, Chemical Equilibrium and Conserved Charges
The composition of the matter in (proto-) neutron stars is constrained by the conditions of electric charge neutrality, chemical equilibrium, and baryon number conservation. Electric charge neutrality and baryon number conservation leads to where q B and q L denote the electric charges of baryons and leptons, respectively, and n B and n L are the respective number densities of these particles. The condition of chemical equilibrium among the different particle species present at a given density reads where µ n , µ e , and µ ν e are the chemical potentials of neutrons, electrons, and neutrinos, respectively. The lepton chemical potentials obey Neutrinos are trapped in the core during the very early stages in the life of a neutron star. Accepted lepton fractions at that stage of stellar evolution are [8] During this stage the system is characterized by three independent chemical potentials, µ n , µ e , and µ ν e , which, together with the hadronic field equations and the conditions of electric charge neutrality and baryon number conservation, determine the composition of hot (proto-) neutron star matter. The condition Y Lµ = 0 has to be imposed in addition since no muons are present when neutrinos are trapped. When a proto-neutron star cools down, the matter in the core becomes transparent to neutrinos quickly, leading to µ ν e = µν µ = 0 for cold neutron stars. In this case, the number of independent chemical potentials reduces from three to two (µ n and µ e ). Figure 1 shows the equation of state (EoS) of neutron star matter at three different temperatures. The matter is composed of neutrons, protons, hyperons, delta isobars, and leptons in chemical equilibrium with each other. Since the number of hyperons and delta isobars in the matter increases drastically with temperature, the pressure at a given density drops for increasing temperatures. In Figures 2 and 3 we show the particle compositions of neutron stars at zero as well as finite temperature. As can be seen, the complexity of the compositions intensifies quickly with increasing temperature, lepton content, and assumptions about the neutrino population.

Quark Matter Modeled with the Non-Local NJL Model
A model widely used to describe the properties of deconfined 3-flavor quark matter is the Nambu-Jona-Lasinio model [26]. In the following, we present a non-local extension of this model (denoted n3NJL). The effective action of this model reads (for details, see Refs. [27,28]) where ψ ≡ (u, d, s) T andm denotes the current quark mass matrix give bym = diag(m u , m d , m s ). The currents j S,P a (x) and j µ V (x) depend on form factors which describe non-local interactions among  the quarks. The quantity T abc in the 't Hooft term accounts for quark-flavor mixing. The current quark massm of up and down quarks and the coupling constants G S and H in Equation (20) are fitted to the pion decay constant, f π , and meson masses m π , m K , m η [29,30]. This leads tom = 6.2 MeV, Λ = 706.0 MeV, G S Λ 2 = 15.04, and HΛ 5 = −337.71. The current mass of the strange quark is treated as a free parameter. Its value is set to m s = 140.7 MeV. It is customary to express the vector interaction, G V , in terms of the strong coupling constant, G S , according to ξ V ≡ G V /G S . The value of ξ V is then typically varied from 0 to around 0.09. The expression of the thermodynamic potential that follows from S E of Equation (20) is given, for the mean-field approximation, by The quantitiesσ f , ω f , andS f denote the scalar, vector, and auxiliary mean-fields of quarks, respectively. The quantity E f denotes the energy-momentum relation of free quarks given by are Gaussian form factors, which describe the nonlocal nature of the interactions among quarks [28]. The quantitiesS f denote auxiliary mean fields [31]. The scalar and vector quark mean-fields are obtained by minimizing the thermodynamic potential, ∂Ω/∂σ f = 0 and ∂Ω/∂ω f = 0. The number densities of quarks are obtained from the thermodynamic potential as n f = ∂Ω/∂µ f .
To determine the equation of state of the system is given in terms of the quark mean-fieldsσ f ,ω f and the neutron and electron chemical potentials µ n and µ e . The pressure p Q and the energy density Q of the system at zero temperature are given by and where Ω 0 is chosen such that P Q vanishes at T = µ = 0 [27,28].

Pasta Structure of the Hadron-Quark Phase
In this section, we provide a brief discussion of the mathematical treatment of the mixed phase region of hadrons and quarks in neutron stars. The starting point is the Gibbs condition for phase equilibrium between hadronic (H) matter and quark (Q) matter [1,2], where P H and P Q denote the pressures of hadronic matter and quark matter, respectively [1,2]. The quantities {φ H } and {φ Q } in Equation (24) collectively denote the field variables of hadronic (σ, ω,ρ) and quark (σ f ,ω f ,S f ) matter as well as the Fermi momenta (k B , k λ ) that characterize a solution of the equations of hadronic matter described in Section 2. The quantity χ ≡ V Q /V describes the volume proportion of quark matter, V Q , in the unknown volume V. χ, therefore, varies between 0 and 1, where 0 means that 0% quark matter is present and 1 means that the matter consist of 100% quark matter. Equation (24) is solved subject to the conditions of global baryon number conservation and global electric charge conservation. The global conservation of baryon charge is expressed as [1,2] where n Q and n H are the baryon number densities of the quark and hadronic phases, respectively. Global electric charge neutrality is given by [1,2] where q Q and q H denoting the electric charge densities of the quark phase and hadronic phase, respectively. As demonstrated by Glendenning [1,2], a mixed phase of hadron-quark matter may arrange itself in such a way that the total energy of the system is minimized. For matter described by global electric charge neutrality, this is the same as minimizing the contributions to the total energy due to phase segregation, which includes the surface and Coulomb energy contributions. Mathematically, the relations for the Coulomb ( C ) and the surface ( S ) energy densities are given by [1,2,32] where q H denotes the charge density of the hadronic phase and q Q stands for the charge density of the quark phase. The symbol r is the radius of the rare phase structures (i.e., quark blobs embedded in hadronic matter at the onset of phase equilibrium). The surface tension between hadronic matter and quark matter is denoted by α(χ) where 0 ≤ χ ≤ 1. The quantity f D (x) in Equation (27) is defined as where D is the dimensionality of the lattice and x ≡ min(χ, 1 − χ). The phase rearrangements will result in the formation of geometrical structures of the rare phase distributed in a lattice that is immersed in the dominant phase. For the sake of simplicity, the rare phase structures are approximated by spherical blobs, rods, and slabs. The spherical blobs populate sites in a three dimensional (D = 3) body centered cubic (BCC) lattice, the rods in a two dimensional (D = 2) triangular lattice, and the slabs in a simple one dimensional (D = 1) lattice. At a volume ratio of χ = 0.5 both hadronic and quark matter exist as slabs in the same proportion, and at χ > 0.5 the hadronic phase becomes the rare phase with its geometry evolving in reverse order (from slabs to rods to blobs) with increasing density [1,2,32]. The determination of the surface tension, α, of the quark-hadro interface is hampered by the fact that a single theory that can accurately describe both hadronic matter and quark matter does not exist. We therefore make use of an approximation proposed by Gibbs where the surface tension is taken to be proportional to the difference in the energy densities of the interacting phases, The quantity L is proportional to the surface thickness which should be on the order of the range of the strong interaction (1 fm). The quantity η denotes a proportionality constant. For the results shown in this work we maintained the energy density proportionality but set the parameter η = 0.08 so that the surface tension falls below 70 MeV fm −2 , which is a reasonable upper limit for the existence of a hadron-quark mixed phase [33]. The surface tension as a function of χ is given in Figure 4 for the nuclear GM1L parametrization introduced in Section 2.  The size of the rare phase structures is given by the radius (r) and is determined by minimizing the sum of the Coulomb and surface energies, ∂(E C + E S )/∂r, and solving for r Rare phase structures are centered in the primitive cell of the lattice, taken to be a Wigner-Seitz cell of the same geometry as the rare phase structure. The Wigner-Seitz cell radius R is set so that the primitive cell is charge neutral, i.e., R(χ) = rx −1/D . In Figure 5 we show r and R as a function of the quark fraction in the mixed phase. Both r and R increase with an increase in the baryonic degrees of freedom, particularly when χ < ∼ 0.5 and the vector interaction is included. The transport properties of a hypothetical quark-hadro Coulomb  [14]. A similar figure for the GM1L parametrization can be found in Ref. [14]. lattice in the core of a neutron star and the neutrino emissivity in this phase have been investigated in Refs. [32,[34][35][36].

Hadron-Quark Lattices in the Cores of Rotating Neutron Stars
In this section we use the equations of state discussed in Sections 2 and 4 to explore the possible existence of hadron-quark lattices in the cores of (cold) rotating neutron stars. Depending on rotational frequency, such stars are highly deformed and experience significant density changes (up to 50%) in their cores which drastically alters their core compositions. The line element that is used to model rotating neutron stars has the form [37] ds 2 = −e 2ν dt 2 + e 2ψ (dφ − ωdt) 2 + e 2µ dθ 2 + e 2λ dr 2 , where ν, ψ, µ and λ denote metric functions and ω stand for the angular velocity of the local inertial frames. Each quantity depends on the radial coordinate r, the polar angle θ, and on the rotating star's angular velocity Ω. The metric functions and the frame dragging frequency are to be computed from Einstein's field equation, where T κσ = T κσ ( , P( )) denotes the energy momentum tensor of the stellar matter, whose equation of state is given by P( ). The other quantities in Equation (33) are the Ricci tensor R κσ , the curvature scalar R, and the metric tensor, g κσ . A strict limit on rapid rotation is given by the Kepler frequency, Ω K , at which mass shedding from the equator terminates stable rotation. This frequency is of great theoretical value because it sets an absolute limit on rapid rotation, but it may not be accessible to neutron stars as other instabilities set tighter constraints on stable rotation [38,39]. As a case in point, even the most rapidly spinning neutron star observed to date, PSR J1748-244ad, spinning at 716 Hz (i.e., 1.4 milliseconds) [40], rotates only at a fraction of the Kepler frequency. The general relativistic expression for the Kepler frequency is obtained by evaluating δ ds 2 = 0 in the star's equatorial plane, which leads to [37,38] where r ≡ ∂/∂r. This equation is to be computed self-consistently in combination with Einstein's field equation at the equator of a rotating neutron star. In Figure 6 we show the gravitational mass versus central energy density relationships for both neutron stars spinning at their Kepler frequencies, as well as non-rotating neutron stars. As one moves along either of these two lines, the baryon mass M B of the stars changes, increasing with larger central density. Isolated neutron stars spinning down due to the loss of energy would follow a path of constant baryon mass, from the Kepler frequency curve down to the non-rotating one. Four such evolutionary paths (marked with arrows) are depicted in Figure 6. One sees that neutron stars with baryon masses less than 2.20 M do not become dense enough during spin down to produce hadron-quark phases in their cores. This is only possible for the most massive neutron stars of our study, which contains a Coulomb lattice made of quark blobs and quark rods if the star's rotational frequency has dropped to 300 Hz (25% of Ω K ). Other quark phases (e.g., quark slabs) would not be present in stable neutron stars.

Quark Matter at Finite Temperature Modeled with the Local PNJL Model
The non-local 3-flavor NJL model has been introduced in Section 4. This model has then be used, together with the model for hadronic matter introduced in Section 2, to determine the possible possible pasta structure of hadrons and quarks in the hadron-quark mixed phase of neutron star matter. In this section, we introduce a local version of the Polyakov-Nambu-Jona-Lasinio (PNJL) model which will be used to study quark matter at zero as well as at finite temperatures [13,[41][42][43]. The Polyakov loop effective potential, which relates quark confinement to color confinement, is included, so that from here on we refer to this model as PNJL. The general expression for the mean-field thermodynamic potential of the model can be written as where Ω 0 is defined by the condition that Ω vanishes at zero temperature and chemical potential, T = µ = 0. Allowing for three independent chemical potentials µ f for the three quark flavors f = (u, d, s), the SU(3) regularized thermodynamic potential in the presence of the quark condensates and the corresponding free term is given by where N c = 3 is the number of quark colors and E f = p 2 + M 2 f is the energy of a quark of flavor f . The color background fields due the coupling to the Polyakov loop are φ c = cφ 3 , i.e., φ r = −φ g = φ 3 and φ b = 0. The sums over the flavor and color indices are over f = (u, d, s) and c = (r, g, b), respectively. The constituent quark masses M f are given by with f , j, k = u, d, s indicating cyclic permutations of the quark flavors. For our calculations we have used the Polyakov loop effective potential as given in Ref. [44], where the coefficients a 1 (T, T 0 ) and a 2 (T, T 0 ) are fixed by Lattice QCD simulations of gluon dynamics and the traced Polyakov loop Φ = [2 cos(φ 3 /T) + 1]/3 [30,45]. The parameter T 0 denotes the critical temperature of the deconfinement phase transition. The effective potential is key to describe the phase transition from the color confined state (T < T 0 , with the minimum of the effective potential at Φ = 1) to the color deconfined state (T > T 0 , with the minima of the effective potential at Φ = 0). T 0 is the only free parameter of the Polyakov loop once the effective potential has been chosen. In this work we will set T 0 = 195.0 MeV [46]. The scalar coupling constant, G s , the 't Hooft coupling constant, H, the quark masses, and the three-momentum ultraviolet cutoff, Λ, are model parameters. The respective values of these quantities are m u = m d = 5.5 MeV, m s = 140.7 MeV, Λ = 602.3 MeV, G s Λ 2 = 3.67 and HΛ 5 = −12.36 [47]. The grand canonical thermodynamic potential Ω v due to vector interactions among quarks is given by where the number density of a quark of flavor f in the mean field approximation is given by As already mentioned in Section 4, the vector coupling constant G v is treated as a free parameter and is generally expressed in terms of the strong coupling constant G s , that is, The minimization of the thermodynamic potential with respect to the quark condensates and the Polyakov-loop color field φ 3 leads to a system of coupled, non-linear equations that are to be solved numerically together with the conditions of electric charge neutrality, baryon number conservation, and chemical equilibrium, similar to the ones discussed in Section 3 for hadronic matter. Once the thermodynamic potential has been computed, the system's pressure, P = −Ω, and the quark number density, n q = ∑ f n f , can readily be determined. The EoS of quark matter then follows as (P, T, n f ) = −P + TS + ∑ f µ f n f , where S = ∂P ∂T and n f = ∂P To construct the hadron-quark phase transition we, as described in Section 5, the Gibbs condition G H (P, T, {φ H }) = G Q (P, T, {φ Q }), where G H respectively G Q denote the Gibbs free energy per baryon of the hadronic (H) and the quark (Q) phase at a given pressure and transition temperature. The Gibbs energy of each phase (i = H, Q) is given by where the sum over j is over all the particles present in each phase. Figure 7 shows the Gibbs energy of hadronic matter and quark matter as a function of pressure. A sharp (Maxwell-like) first-order phase transition is assumed. Values of ξ v = 0 and 0.1 have been chosen for the ratio of the vector-to-scalar quark coupling constant to demonstrate the dependence of the Gibbs energy on the repulsion among quarks. The points where the curves cross each other mark the location of phase equilibrium between hadronic and quark matter. Figure 8 shows the composition of cold quark matter in chemical equilibrium computed for ξ v = 0 and 0.1. It turns out that the 't Hooft term responsible for quark flavor mixing alters the quark population significantly, which, in turn, severely changes the populations of leptons and muons.
In Figure 9, we show the hadron quark composition of hot neutron star matter. Neutrinos are included. This composition, therefore, corresponds to matter that may exist in the cores of proto-neutron stars. The shaded areas indicate a sharp (Maxwell) hadron-quark phase transition, computed for the DD2 nuclear parametrization. Three-flavor quark matter is described by the local PNJL model described above. As can be seen, quark flavor mixing shifts the hadron-quark phase transition toward a smaller baron number density. This trend is to be expected from the results for the particle densities shown in Figure 8, which shows that flavor mixing shifts the presence of quarks to lower densities.

Summary and Future Prospects
The purpose of this paper is to provide a short overview of the role of quark matter for neutron star matter at different temperatures. As shown, the hadron-quark composition of such matter depends critically on temperature. For proto-neutron star matter neutrinos are to be included in the treatment too, which also impacts the various particle fractions present at a given temperature and density. If hadron-quark matter exists in cold neutron star, it may form a Coulomb lattice. As demonstrated here (for a non-local version of the NJL model), this lattice would be made almost entirely of spherical quark blobs. The size of the lattice depends on the central density of a neutron star, which, most interestingly, links the size of the lattice to the spin-frequency of a neutron star. We find that the lattice could be produced during spin-down of massive pulsars, once their spin frequency has dropped below around 300 Hz (∼ 25% of the Kepler frequency). The study of quark matter with a local PNJL shows that the 't Hooft quark flavor mixing term leads to non-negligible changes in the particle composition of (proto-) neutron stars made of hadron-quark matter.
We have analyzed the role of the 't Hooft mixing term, associated with the axial U(1) symmetry anomaly of QCD, for quark matter that is locally electrically charge neutral, color neutral, and in chemical equilibrium (i.e., matter under neutron star matter conditions). This symmetry is expected to be restored (not necessarily together with chiral symmetry) at high temperature, which implies that the mixing coupling constant H → 0. For high-density matter as existing in the cores of neutron stars, the mixing (six-point interaction) term is less relevant than the four-point scalar and vector interaction terms (see Equations (20) and (21)). At the mean-field level, the axial U(1) symmetry anomaly of QCD does not change the thermodynamics of the quarks [48].
Our results show that the inclusion of the mixing term lowers the pressure at which the quark-hadro phase transition occurs. It will be interesting to explore how the mixing term modifies the possible formation of a pasta phase in the hadron-quark mixed phase if color super-conductivity is taken into account, since the mixing term can destabilize a color superconductor by preventing the formation of diquark condensates [49]. The existence of a solid hadron-quark pasta phase in the core of a neutron star may be linked to sudden pulsar spin-ups (glitches) and the subsequent healing of the pulsar period. A very prominent example of a pulsar that displays regular glitch activity is PSR J0537-6910. The glitch rate of this pulsar is 3.2 yr −1 [50]. A link between pulsar glitches and solid neutron-matter cores inside of neutron stars has been pointed out may years ago already [51][52][53]. This suggests that one should investigate the connection between a solid hadron-quark pasta phase and the pulsar glitch phenomenon. In particular, it needs to be studied whether or not the sudden release of elastic energy stored in a solid hadron-quark pasta phase could explain macro-glitches (which are related to sudden changes in the rotational frequency of a pulsar on the order of ∆Ω/Ω ∼ 10 −6 ) at the rate (months to years) at which pulsar glitches are typically observed. Another issue that needs to be investigated concerns the role of neutron superfluidity. Since neutron stars with and without hadron-quark matter cores contain nuclear crusts with identical compositions (but varying thickness), both types of stars possess nuclear crusts containing 1 S 0 neutron superfluids. The situation is less clear for the neutrons in the cores of both types of stars. It needs to be studied if both types of stars may possess 3 S 1 neutrons superfluids in their cores, or whether this phase is restricted to neutron stars without quark matter cores.
Gravitational-wave detectors such as LIGO and VIRGO have opened up a new observational window on the inner workings of neutron stars. This is particularly the case for the neutron stars in the elliptic galaxy NGC 4993 that produced the gravitational-wave event GW170817 observed with LIGO and VIRGO in 2017. The detection of the gravitational and electromagnetic waves produced by the coalescence of these neutron stars has offered an unprecedented opportunity for imposing constraints on the quark-hadro composition of neutron stars and the equation of state associated with such objects (see Refs. [54][55][56], and references therein).
Author Contributions: The authors contributed equally to the theoretical and numerical aspects of the work presented in this paper.

Conflicts of Interest:
The authors declare no conflict of interest.