Phase transitions in neutron stars and their links to gravitational waves

TL;DR

This paper reviews how gravitational wave observations, especially from events like GW170817, can reveal phase transitions in neutron star cores, potentially indicating quark matter presence through g-mode signals.

Contribution

It provides a comprehensive review of neutron star matter, the equation of state, and how multimessenger observations constrain phase transitions within neutron stars.

Findings

Detection of g-modes could signal hadron-quark phase transitions.

Gravitational wave data constrains neutron star core composition.

Observations suggest possible existence of quark matter in neutron stars.

Abstract

The recent direct observation of gravitational wave event $GW170817$ and its $GRB170817A$ signal has opened up a new window to study neutron stars and heralds a new era of Astronomy referred to as the Multimessenger Astronomy. Both gravitational and electromagnetic waves from a single astrophysical source have been detected for the first time. This combined detection offers an unprecedented opportunity to place constraints on the neutron star matter equation of state. The existence of a possible hadron-quark phase transition in the central regions of neutron stars is associated with the appearance of g-modes, which are extremely important as they could signal the presence of a pure quark matter core in the centers of neutron stars. Observations of g-modes with frequencies between 1 kHz and 1.5 kHz could be interpreted as evidence of a sharp hadron-quark phase transition in the cores of neutron stars. In this article, we shall review the description of the dense matter composing neutron stars, the determination of the equation of state of such matter, and the constraints imposed by astrophysical observations of these fascinating compact objects.

Tables

Table 1. Table 1: Properties of nuclear matter at saturation density for the hadronic DDRMF parametrizations, GM1L [ 82 ] and DD2 [ 83 ] .

Saturation Property	GM1L	DD2
$n_0$ (fm-3)	0.153	0.149
$E_0$ (MeV)	$-16.30$	$-16.02$
$K_0$ (MeV)	300.0	242.7
$m^{\ast }/m_N$	0.70	0.56
$J$ (MeV)	32.5	32.8
$L_0$ (MeV)	55.0	55.3

Table 2. Table 2: The minimum estimated energy (in units of M ⊙ c 2 ) M_{\odot}c^{2}) required in each mode in order to lead to a detection with signal-to-noise ratio of S/N = 8 from a pulsating NS at different distances (local event or at the distance of the Virgo cluster). We use typical order of magnitude frequencies and damping times for each mode.

Detector	Distance	$f$-mode	$p_1$-mode	$g$-mode
LIGO/Virgo	10 kpc	$8.9\times {10}^{-8}$	$3.2\times {10}^{-6}$	$8.0\times {10}^{-9}$
LIGO/Virgo	15 Mpc	0.2	7.2	$1.8\times {10}^{-2}$
Einstein	10 kpc	$8.9\times {10}^{-10}$	$7.2\times {10}^{-8}$	$1.8\times {10}^{-10}$
Einstein	15 Mpc	$2\times {10}^{-3}$	0.16	$4\times {10}^{-4}$

Equations

$${\mathcal{L}}_B$$

$$g_{\rho B}(n)=g_{\rho B}(n_0)\exp [-a_{\rho }(x-1)],$$

$$\text{vspace}2.0cm{g}_{iB}(n)=g_{iB}(n_0)a_i\frac{1+b_i(x+d_i{)}^2}{1+c_i(x+d_i{)}^2},\mathrm{for}i=\sigma ,\omega ,$$

$$\mathcal{L}$$

$$j_a^S(x)$$

$$j_a^P(x)$$

$$\mathrm{d}{s}^2$$

$${\xi }^r$$

$${\xi }^{\theta }$$

$${\xi }^{\varphi }$$

$$\frac{E_{\mathrm{gw}}}{M_{\odot }{c}^2}=3.47\times 10^{36}{\left(\frac{S}{N}\right)}^2\frac{1+4Q^2}{4Q^2}{\left(\frac{D}{10\mathrm{kpc}}\right)}^2{\left(\frac{f}{1\mathrm{kHz}}\right)}^2{\left(\frac{S_n}{1{\mathrm{Hz}}^{-1}}\right)}^2,$$

Full text

1Grupo de Gravitación, Astrofísica y Cosmología, Facultad de Ciencias Astronómicas y Geofísicas, Universidad Nacional de La Plata, Paseo del Bosque S/N, 1900, La Plata, Argentina.

2 CONICET, Godoy Cruz 2290, 1425, CABA, Argentina.

3 IFLP, UNLP, CONICET, Facultad de Ciencias Exactas, Diagonal 113 e/63 y 64, 1900, La Plata, Argentina.

4 AIM, CEA, CNRS, Université Paris-Saclay, Université Paris Diderot, Sorbonne Paris Cité, F-91191 Gif-sur-Yvette, France.

5 Instituto Argentino de Radioastronomía (CCT-La Plata, CONICET; CICPBA), C.C. No. 5, 1894 Villa Elisa, Argentina.

6 Department of Sciences, Wentworth Institute of Technology, 550 Huntington Avenue, Boston, MA 02115, USA.

7 Universidade Federal do ABC, Centro de Ciências Naturais e Humanas, Avenida dos Estados 5001- Bangú, CEP 09210-580, Santo André, SP, Brazil.

8 Department of Physics, San Diego State University, San Diego, CA 92182, USA.

9 University of California at San Diego, La Jolla, CA 92093, USA.

1 Introduction

One of the principal aims of Hadron Physics is to provide a comprehensive theoretical description of hadrons at the quark level, which constitute the fundamental building blocks of matter. The hadrons are composite objects made up of three quarks (baryons) or of quark-antiquark pairs (mesons). Light baryons are formed by different combination of three light quarks, $u$ (up), $d$ (down), and $s$ (strange)111Despite the fact that heavy baryons and mesons composed of $c$ (charm) and $b$ (bottom) quarks have been produced [1], we will ignore them. Quark chemical potentials inside neutron stars are at most on the order of 500 MeV, too small to create a population of these states.. The gauge field theory with the local symmetry group SU(3) which describes the strong interactions of colored quarks and gluons is known as Quantum Chromodynamics (QCD). In the framework of QCD, neutrons and protons, which are the building blocks of ordinary atomic nuclei and of nuclear matter, are the fundamental degrees of freedom of strongly interacting matter at low temperatures and densities. At high temperatures and/or densities, QCD predicts a phase transition of ordinary baryonic matter, where quarks are confined inside of hadrons, to a new state of matter called quark-gluon plasma (QGP), where the hadronic boundaries have disappeared and quarks and gluons roam about freely.

It is understood that the universe was filled with such a plasma immediately after the Big Bang. Moreover, it is being speculated that quark matter may play a role in the demise of massive stars in supernova explosions, could exist temporarily in the cores of (hot) proto-neutron stars, and may form a permanent component of matter in the cores of (cold) neutron stars (hereafter, NSs). It is therefore extremely important to explore possible phase transitions in super-dense matter, understand the thermodynamic behavior of the QGP, and to delineate realistic models for the equation of state (EoS) of such matter.

Due to the property of asymptotic freedom the interactions between quarks become increasingly weaker with decreasing separation. This is the regime of high-momentum transfer (high densities and/or high temperatures), where perturbative techniques may be used to determine the properties of quark matter. However, in the regime of low momentum transfer (low densities and low temperatures, increasing separation), the strong coupling constant becomes large and the theory is highly non-perturbative. In this way, it is theoretically well-established that confinement and chiral symmetry breaking cannot be obtained in a simple perturbation-theoretical analysis of QCD. Thus, to describe hadron masses or hadron decay constants it is essential to develop formalisms that allows one to study the theory in the non-perturbative region. A standard method to solve non-perturbative QCD analytically does not exist. As an alternative, lattice QCD (LQCD) [2] has been developed, which explores the dynamics of quarks and gluons numerically on a discretized space-time grid. A very important consequence of the discretization is that the measurement of the functional integrals are products of a large but countable number of differentials.

Due to ever increasing computing power it has become possible to decrease the spacing of the space-time grid and increase the size of the system so as to approach both the continuum limit as well as the thermodynamic limit. At the same time it has become possible to sample a sufficiently large number of configurations in such numerical simulations so that the statistical errors decrease significantly. However, the Monte Carlo-type methods commonly used to evaluate the partition function are reliable only in the case where the chemical potential is zero. The extension of the calculations to finite chemical potentials present serious difficulties (see Refs. [3, 4], and references therein). A path to treating QCD at finite chemical potential is to resort to effective models of QCD, that is, models based on Lagrangians that are simplified compared to the full QCD Lagrangian while maintaining the fundamental properties of the full theory.

Studies of the QCD phase diagram predict that the crossover shown by LQCD calculations at vanishing chemical potential becomes a first-order phase transition at intermediate temperatures and high baryon chemical potentials [5]. In particular, this suggests that the matter in the ultra-dense cores of NSs, which is compressed to densities several times that of nuclear matter ($n_{0}\sim$ 0.16 fm*-3*), may undergo such a transition.

NSs are dense, neutron-packed remnants of massive stars that blew apart in supernova explosions [6]. They are typically about 20 kilometers across and spin rapidly, often making many hundred rotations per second. The record holder among the rapidly spinning neutron stars (pulsars) is currently PSR J1748-2446ad which rotates at 716 Hz [7].

Many NSs are radio pulsars, emitting radio waves that appear from the Earth to pulse on and off like a lighthouse as the star rotates at very high speeds. NSs in X-ray binaries accrete material from a companion star and flare to life with tremendous bursts of X-rays. Pulsars can also be powered by extreme magnetic fields. Some NSs, known as magnetars, have been shown to have surface magnetic fields of 1014-1015 Gauss [8, 9, 10, 11]. It is believed that the magnetic fields in their cores are even stronger by several orders of magnitude. From model calculations it is known that the cores of NSs contain most of their gravitational masses. The nuclear EoS associated with these regimes, however, is still largely unknown since the description of the matter in the ultra-dense cores becomes increasingly uncertain with density. This concerns both the many-body techniques that are used to model such matter as well as open issues regarding the compositional structure of the core, where aside from neutrons and protons, hyperons, the delta resonance, meson condensates, and/or deconfined quarks may exist [12]. These features make NSs unique astrophysical laboratories for a wide range of physical studies, which range from the exploration of nuclear processes on the surfaces and inner crusts, to the exploration of novel states of matter in the cores of such objects.

Determining the NS EoS, that is, the functional relationship between pressure and energy density, allows one to obtain the relationship between the mass and radius of a NS as well as other key properties such as the tidal deformability and the moment of inertia. Strong constraints on the appearance of additional degrees of freedom in neutron star matter, and hence on the EoS used to describe such matter, follow from the discovery of NSs with masses of around 2M*⊙* [13, 14, 15, 16]. Moreover, important additional constraints on the EoS of NSs follow from the recent direct detection of gravitational waves (GWs) emitted from a merger of two NSs–the gravitational-wave event GW170817–whose electromagnetic counterpart (GRB 170817A) has been observed as well [17, 18, 19, 20, 21, 22, 23]. This event serves, for instance, to explore the tidal deformability and radii of neutron stars and to impose further constraints on the nuclear EoS [24, 25, 26, 27]. Last but not least we mention non-radial oscillation modes of NSs. Such modes may also lead to the emission of GWs (see Refs. [28, 29, 30] and references therein), which could shed additional light on the behavior of the matter in the interiors of NSs.

In this article, we shall review the composition of cold catalyzed neutron star matter, the determination of the equation of state of such matter, and the constraints imposed by astrophysical observations of NSs. We shall focus on the nature and properties of matter at supra-nuclear densities and low temperatures $T\leq 1$ MeV ($\simeq 10^{10}$ K). The article is organized as follows. In section 2 we discuss the QCD phase diagram, effective models (and their extensions) for the description of hadronic and quark matter, and the constraints on quark matter imposed by LQCD studies. In section 3 we discuss the nature of the hadron-quarks phase transition. In section 4 we discuss constraints on the nuclear EoS from radioastronomical observations of radio pulsars and X-ray observations of isolated NSs and X-ray binaries. In section  5 we briefly describe the post GW170817 stage and discuss the NSs merger numerical simulation and the constraints to the EoS from the gravitational-wave event GW170817. In section 6 we present some basic aspects of NSs oscillations and their link to GW emission from isolated NSs and theoretical results obtained from a representative hybrid EoS. We discuss the relationship between a possible phase transition in the inner core of a NS and GWs. We also provide an estimate of the intensity of a GW signal linked to non-radial oscillation modes. A summary and outlook is given in section 7.

2 Probing the dense matter EoS I. Terrestrial laboratories constraints

2.1 Probing dense matter beyond nuclei

In the present review the focus is on the zone of the QCD phase diagram of high chemical potential and low temperatures. To study the behavior of matter subjected to such extreme conditions, as it could be found inside NSs, we should take into account the degrees of freedom relevant to low and high densities. From the theoretical point of view, it was originally believed that the QCD phase diagram consisted of only two phases: i) the hadronic phase of confined quarks and gluons with broken chiral symmetry, and ii) the deconfined phase with restored chiral symmetry, known as QGP. Only at the end of the 90’s the possibility of new intermediate phases in the QCD phase diagram, like color superconductivity in which quark Cooper pairs are formed, began to be accepted [31] and later studied in detail in compact stars environments [32]. It has also been suggested that there may be non-homogeneous phases [33, 34, 35, 36], and that there may be a phase of confined quarks with restored chiral symmetry (quarkyonic phase) [37, 38, 39]. A schematic version of the conjectured phase diagram of QCD is presented in figure 1.

According to LQCD simulations, a discontinuous chiral transition is expected at some critical end point (CEP) in the phase diagram of QCD, immediately after the crossover line in figure 1. The location and even the existence of the CEP are still open problems [40]. Recent advances have been made in understanding the phase diagram of QCD in the presence of strong magnetic fields. The intensity of the magnetic field can affect the CEP location. Different model calculations and LQCD simulations predict magnetic catalysis at zero temperature (the quark condensate increases as a function of the magnetic field) [41]. The so-called inverse magnetic catalysis in 2$+$1 flavors, in which the magnetic field suppresses the light $u$ and $d$ quark condensates close to the chiral restoration temperature, has been obtained using LQCD data [42]. The phenomenon of inverse magnetic catalysis found in LQCD calculations was also verified theoretically considering two flavors of quarks [43]. The nature of strong interactions under extreme external magnetic fields and its effect on the QCD phase diagram remains an open issue.

As mentioned in the introduction, the calculations of the LQCD are limited to zero chemical potential for realistic quark masses as a consequence of the limitations of Monte Carlo techniques to evaluate the partition function. Due to the sign problem222The interpretation of the canonical density operator as an evolution operator in imaginary time is useful to study QCD at finite density, making the effective action a complex quantity. This fact prevents comparing the probabilities associated with different configurations., it is only possible to investigate the QCD phase diagram at low chemical potentials developing theoretical alternative methods and using accurately extrapolation of the QCD data. Using a Taylor expansion [44, 45, 46] or analytical continuation from imaginary chemical potential [47, 48, 49, 50], the behavior of the transition temperature at finite density in the low density region of the QCD phase diagram has been studied (see Ref. [3] for a complete review of this issue). At zero chemical potential some observables (chiral condensate and its susceptibility) that are related to the spontaneous chiral symmetry breaking give a critical temperature around $T_{c}(\mu=0)\simeq 155$ MeV, obtained from extrapolations to the continuum limit for $2+1$ quark flavors [51, 52, 53]. In addition, some experimental data from the Beam Energy Scan at RHIC by the STAR collaboration and by the ALICE Collaboration at the LHC show some signals of deconfinement phase transition from the determination of the chemical freeze-out temperatures and baryochemical potentials [54, 55, 56, 57, 58].

Although the low density regime of the QCD phase diagram can be probed by the LQCD calculations through the approaches mentioned before, the formulation of LQCD at higher values of the chemical potential fails. Thus, phenomenological approaches are necessary to interpret the recent experimental results at NICA-FAIR energies [59]. These approaches basically consist on the assumption the production of new degrees of freedom in the early stage of nucleus-nucleus collisions is a statistical process. The broad set of experimental data obtained at NICA-FAIR energies agrees with the interpretation of the results using statistical models based on such statistical process. This enables to extract information from the transition zone between the hadronic and quark-gluon phases at the region of the QCD phase diagram covered by these experiments.

Due the separation of the energy scales in the treatment of strong interactions (low and high energy regimes), powerful tools to study the dense matter behavior are effective theories that keep some of the features of the QCD. The most popular effective models used for hadron interactions are the Nambu–Jona Lasinio (NJL) models [60, 61, 62, 63]. Other models, like the Chiral Sigma Model [64, 65, 66, 67], the Polyakov–Quark-Meson Model [68, 69, 70, 71], or the Field Correlator Method [72, 73, 74, 75, 76] have also been used to describe the quark matter phase.

So far, theoretical physicists have not managed to find an effective QCD model capable of reproducing both the hadron and quark regimes in a single EoS. Matter becomes simpler to describe at very high densities due to the asymptotic freedom property, but very complicated in the opposite limit due to the confinement. What happens in the middle is speculation that comes from combining the data extracted from heavy ion collisions and experimental physics with different theoretical models.

NSs can provide a unique window to explore strong interactions under extreme conditions. Thus, to model the matter in the inner and outer cores of NSs, where the bulk of their mass is concentrated, theoretical astrophysicists consider two distinct phases composed of different types of particles. In this way, a hybrid EoS will be a multicomponent system composed of two independent phases: the hadronic phase at low densities, and the quark phase at high densities.

2.2 Hadronic phase at low densities

To describe the interaction between baryons, we use Quantum Hadrodynamics, the relativistic quantum field theory of interacting many-body systems based on hadronic degrees of freedom. Such description can be achieved by modeling the interaction between two baryons as a coupling of the baryons to meson fields, which reproduce the fundamental attractive and repulsive nature of the strong force. This approach is similar to the one-boson-exchange modeling of the nuclear interaction that very successfully describes the wealth of nucleon-nucleon scattering data. However, instead of fixing the meson-nucleon coupling constants to scattering data, they are fit to the properties of bulk nuclear matter at the nuclear saturation density extracted from recent experimental or theoretical results from nuclear physics.

A popular model of Quantum Hadrodynamics is known as the Walecka model [77, 78, 79] in which all the baryon states, $B$, (nucleons, hyperons which possess strange quarks and delta resonances) interact via the exchange of scalar, vector and isovector mesons ($\sigma$, $\omega$, $\rho$, respectively). In this model, the meson field strengths are taken to be equal to their mean values in the so called relativistic mean-field (RMF) approximation, and the baryon-baryon interactions are described in terms of meson fields. Thus, the scalar meson $\sigma$ describes attraction between baryons, the vector meson $\omega$ describes repulsion, and the isovector meson $\rho$ is necessary for describing baryon-baryon interactions in isospin asymmetric systems. The EoS of isospin asymmetric nuclear matter is crucial in understanding the structure of NSs, specially for their cooling history [80].

The Lagrangian of nuclear matter constructed on the basis of the original Walecka model and extended by the inclusion of nonlinear scalar self-interaction terms is given by [81],

[TABLE]

where $g_{\sigma B}(n)$, $g_{\omega B}(n)$ and $g_{\rho B}(n)$ are the (density dependent) meson-baryon coupling constants, $n$ is the baryon number density, ${\bm{\tau}}$ are the Pauli isospin matrices, and $\gamma^{\mu}$ are the Dirac matrices.

In the standard RMF the meson-nucleon coupling constants, $g_{iN}(n)=g_{iN}(n_{0})$, $i\in\{\sigma,\omega,\rho\}$ and the coefficients of the nonlinear scalar self interactions $b_{\sigma}$ and $c_{\sigma}$, are set so as to reproduce the properties of isospin symmetric nuclear matter. Symmetric nuclear matter, $n_{p}=n_{n}$, is a reasonable approximation for the interior of heavy atomic nuclei and can be parameterized to reproduce empirically determined properties of nuclear matter at the nuclear saturation density, $n_{0}$, thereby constraining the EoS. These saturation properties are the binding energy per nucleon ($E_{0}$), the nuclear incompressibility ($K_{0}$), the isospin asymmetry energy ($J$) and its slope ($L_{0}$), and the effective mass ($m^{*}/m_{N}$), all fixed to their $n_{0}$ values. However, the standard RMF approximation is not parameterized to fix $L_{0}$ at the saturation density, a quantity that has become more tightly constrained in recent years and may have a significant effect on the composition and properties of NSs (see Refs. [84, 85] and references therein). One of the ways to set $L_{0}$ is to consider medium effects of the effective interaction, choosing density functionals for the meson-baryon couplings. More specifically, the density dependence is extracted from Dirac-Brueckner calculations in which nucleon interactions in the nuclear medium are obtained from nucleon-nucleon potentials consistent with scattering experiments [86]. Being based on first principles calculations, density dependent coupling constants give more empirical robustness to the density-dependent relativistic mean-field (DDRMF) than in the standard RMF. The functional forms for the density dependent isovector- , scalar- and vector-meson-baryon coupling constants are given by [87, 88]

[TABLE]

where $x=n/n_{0}$. This choice of parametrization accounts for nuclear medium effects by making the meson-baryon coupling constants dependent on the local baryon number density [89]. The parameters $a_{i}$, $b_{i}$, $c_{i}$ and $d_{i}$ are adjusted by comparing model predictions to experimental data such as binding energies, charge and diffraction radii, spin-orbit splittings and neutron skin thickness of finite nuclei [90]. The density dependence of the meson-baryon couplings also eliminates the need for the non-linear self-interactions in the DDRMF Lagrangian. Additionally, as a consequence of the density dependent coupling, a rearrangement contribution term in the chemical potential and pressure is mandatory to maintain thermodynamical consistency [91].

The determination of the meson-hyperon coupling constants at the saturation density is more complicated. Since hyperons are not present in the atomic nuclei, it is not possible to adapt the meson-hyperon coupling constants to any property of the nuclear matter as in the case of nucleons. To scape this problem, the simplest proposal is to suppose that the coupling constants have the same value for all baryons. This is the so-called universal coupling. Another alternative is to determine the meson-hyperon coupling constants based on the spin-isospin flavor symmetry SU(3) [82].

The saturation properties associated with the DDRMF parameterizations are given in table 1. Results using EoS calculated employing the GM1L and DD2 parametrizations are shown and discussed in sections 5 and 6.

2.3 Quark phase at low, intermediate and high densities

As it has been mentioned in 2.1, to study the physical processes among quarks in the low momentum range of QCD, it is possible to use effective models that are described by a Lagrangian that accounts for the main features of QCD at low energies. For example, NJL-like models describe interactions between constituent quarks and give a simple scheme to study spontaneous chiral symmetry breaking, a key feature of QCD in its low density/temperature phase, and its manifestations in hadronic physics, such as dynamical quark mass generation, the appearance of a quark pair condensate and the role of pions as Goldstone bosons. In its simplest form, the NJL model uses only a local scalar-isoscalar and pseudoscalar-isovectorial interaction between fermions. The problem with such local interactions is, however, that the model should be regularized in order to avoid divergences in quantities such as the quark self-energies and the meson masses. In the 90’s, it was suggested to use non-local interactions rather than local interactions [92]. A natural way to introduce non-locality in the quark-quark interaction is to use a phenomenologically modified (effective) gluon propagator. These modified propagator leads to a non-local regulator for quark currents that have the advantage to decay at large momentum transfer, thus enabling to numerically perform the integrals to infinity without using a momentum cut-off. In general, these models lack confinement, which is one of the key properties of QCD. The Polyakov loop dynamics is a suitable framework to achieve quark confinement in these kind of models, in particular at finite temperatures and vanishing densities.

Polyakov-NJL (PNJL) models have an effective potential ${\cal{U}}(\Phi,T)$ defined for the (complex) $\Phi$ field, which is conveniently chosen to reproduce, at the mean field level, results obtained in LQCD calculations. Different ansatz for $\cal{U}$ are available in the literature (see, for example, Refs. [93, 94, 95]). The effective potential exhibits the feature of a phase transition from color confinement ($T<T_{0}$, the minimum of the effective potential being at $\Phi=0$) to color deconfinement ($T>T_{0}$, the minimum of the effective potential occurring at $\Phi\neq 0$), where the parameter $T_{0}$ is the critical temperature for the deconfinement phase transition. This parameter is the only free parameter of the Polyakov loop once the effective potential is chosen. If a nonlocal extension of the SU(3) PNJL model is considered, the Lagrangian can be written as [96, 97, 98, 99, 100, 101, 102]:

[TABLE]

where $j_{a}^{S,P}(x)$ are the scalar and pseudoscalar currents given by

[TABLE]

with the regulator $\tilde{Rg}(z)$ written into momentum space, $Rg(p)=\mathrm{exp}(-p^{2}/\Lambda^{2})$, and $\Lambda$ a cutoff adjusted to set the range of the non-locality. In this model, the quark current masses and coupling constants have be chosen so as to reproduce the phenomenological values of pion decay constant $f_{\pi}$, and the meson masses $m_{\pi}$, $m_{\rm K}$ and $m_{\eta^{\prime}}$. Setting the strange quark mass to the updated value of $m_{s}=95$ MeV [1], the parameters $\Lambda=1071.38$ MeV, $G_{s}\Lambda^{2}=10.78$, $H\Lambda^{5}=-353.29$ and $m_{u}=m_{d}=3.63$ MeV [101] are obtained.

To extend the model at finite temperature and study different thermodynamical quantities, the Matsubara imaginary time formalism is used [103]. Once the grand canonical potential is obtained for the mean field approximation after the bosonization of equation 2, the phase diagram of figure 2 can be calculated, where the traced (blue) line is a crossover phase transition defined by the critical points having a peak in their chiral or baryonic susceptibility, the solid (black) line indicates the first order phase transition defined by the critical points where the two minims of the grand canonical potential are equal, and the solid dot indicates the location of the critical end point (CEP) where the first order phase transition line ends and a second order phase transition occurs. Freeze-out and lattice data points mentioned previously are shown in figure 2. These data points have also been compared to several analytical extensions of LQCD in [3]. It is worthly to remark that in both local and nonlocal extensions of NJL models there are absence of a hadronic phase of confined quarks. The overlayed chemical freeze-out data in the phase diagram of the non-local model are presented only by completeness. Chemical freeze-out parameters are analyzed from particle ratios in heavy-ion collision experiments using the statistical model of the hadron resonance gas [4]. The non-local PNJL model at zero temperature has been used to describe quark matter in compact stars interiors in several previous works [104, 105, 106, 107].

Recently, the Field Correlator Method, a non-perturbative approximation of QCD parametrized through the gluonic condensate and the quark-antiquark static potential for long distances (string confinement), has also been used as an alternative effective model to study quark matter [75, 76, 108, 109, 110].

Some thermodynamical functions can be determined from the standard derivation of the grand canonical potential. In figure 3 the results obtained with the non-local SU(3) PNJL model are compared to the continuum extrapolation for 2+1 flavors in LQCD obtained with the stout action by the Wuppertal-Budapest Collaboration [111] and the results of HISQ Action obtained by HotQCD Collaboration [112]. For completeness the calculations have been performed using several model parameterizations which consider different ranges of strange quark masses, for local (Hatsuda and Kunihiro (HK) [62] and Lutz, Klimt, and Weise (LKW) [113]) and non-local (for $m_{s}=140.7$ MeV [104, 105] and $m_{s}=95$ MeV [101]) PNJL models. The results obtained from the FCM has also been included.

Although the existence of color superconductivity had been suggested in the mid-70’s [114], the possibility of this alternative phase of quark matter at low temperatures and high densities was mostly accepted twenty years later. Then, it was argued that its ground state consists of quark pairs which form a color superconductor where the attractive interaction between quarks emerges naturally through the exchange of gluons [31]. As seen in the Bardeen-Cooper-Schrieffer (BCS) theory of ordinary superconductivity, quarks with equal but opposite momenta at their respective Fermi surfaces form Cooper pairs to lower the energy of the system. Color superconductivity in the interiors of compact stars has been studied by many authors over the last 15 years [71, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124].

In particular, color superconducting NJL 2SC+s quark matter cores in HSs has been studied in Ref. [106]. In the NJL 2SC+s phase, green and red up and down quarks form pairs and are embedded in a gas of free strange quarks. Results using this high density EoS are shown and discussed in section 5.

3 Dense hybrid-matter in NSs structure

Different EoSs define different types of stars: NSs, HSs, and strange stars [125]. Allowing for the possibility of quark deconfinement in NSs, a hadron-quark phase transition may commence when the pressure of the quark phase equals that of the hadronic phase. The nature of the phase transition depends on the hadron-quark surface tension, $\sigma_{HQ}$, between the two phases that is still quite uncertain. Recent works have typically placed $\sigma_{HQ}\lesssim 30$ MeV/fm2, though there are suggestions that values greater than 100 MeV/fm2 could be possible too [126, 127, 128, 129]. If $\sigma_{HQ}\gtrsim 70$ MeV/fm2 the quark-hadron phase transition will be one of constant pressure with an equation of state that is discontinuous in energy density [130, 131]. The result is a sharp interface between phases of pure hadronic matter and pure quark matter at a given NS radius, as shown schematically in figure 4(a).

The phase transition is achieved by applying the Maxwell construction (at zero temperature, $T$) and imposing $P_{H}$ ($\mu_{Bt},T=0$) = $P_{Q}$ ($\mu_{Bt},T=0$), where $\mu_{Bt}$ is the baryon chemical potential at the transition. $P_{H}$ and $P_{Q}$ represent the pressure of the hadronic and the quark phase respectively. The condition for a stable phase transition from hadron to quark matter is fulfilled if the Gibbs free energy of quark phase is less than the Gibbs free energy of the hadron phase, for pressures above the transition pressure. If $\sigma_{HQ}\lesssim 70$ MeV/fm2 the phase transition results in the formation of a stable coexistent (mixed) phase, gradually converting NS matter from hadronic matter to deconfined quark matter with increasing density, as shown in figure 4(b). This phase transition satisfies the Gibbs condition for phase equilibrium, $P_{H}$($\mu_{n}$,$\mu_{e},T=0$) = $P_{Q}$($\mu_{n}$,$\mu_{e},T=0$), where $\mu_{e}$ is the electron chemical potential. The isospin restoring force favors a positively charged hadronic phase as an increased number of protons reduces isospin asymmetry. This necessarily results in a negatively charged quark phase, electric charge neutrality being achieved globally.

The presence of quark matter enables the hadronic regions of the mixed phase to arrange themselves to be more isospin symmetric than in the purely hadronic phase by transferring charge to the quark phase in equilibrium with it. The symmetry energy will be lowered thereby at only a small cost in rearranging the quark Fermi surfaces. The electrons play only a minor role when neutrality can be achieved among the baryon-charge carrying particles. The stellar implication of this charge rearrangement is that the mixed phase region of the star will have positively charged regions of nuclear matter and negatively charged regions of quark matter. Because of the competition between the Coulomb and the surface energies associated with the nuclear and quark regions, the mixed phase may develop geometrical structures, schematically illustrated in figure 5, as is similarly expected of the sub-nuclear liquid-gas phase transition [132, 133]. The quantity $\chi\equiv V_{q}/V_{TOTAL}$ in such figure denotes the volume proportion of quark matter, $V_{q}$, in the mixed phase. By definition, $\chi$ therefore varies between 0 and 1, depending on how much hadronic matter has been converted to quark matter.

Competition between the Coulomb and surface energies establishes the shapes, sizes, and spacings of the rare phase in the background of the other in order to minimize the lattice energy [134]. The change in energy accompanied by developing such geometrical structures is likely to be very small in comparison with the volume energy [135] and, thus, may not affect the global properties of a NS. However, the geometrical structure of the mixed phase may be very important for irregularities like pulsar glitches in the timing structure of pulsar spin-down as well as for the thermal and transport properties of NSs [136]. The effect of a crystalline quark-hadron mixed phase on the neutrino emissivity from the cores of NSs has been recently studied in Refs. [136, 137].

Summarizing, the phase transition that might occur between the inner and outer cores of a NS could be described within different scenarios, depending on the value of $\sigma_{HQ}$. There are two limiting cases, the Maxwell construction, which assumes $\sigma_{HQ}\to\infty$ and the bulk Gibbs construction, which assumes $\sigma_{HQ}\to 0$. For intermediate values of $\sigma_{HQ}$, the full Gibbs formalism, i.e. a full pasta phase calculation, must be used taking into account the formation of phase structures in the mixed phase [138]. Recently, several works have included finite-size effects at the hadron-quark phase transition in a phenomenological way, mimicking the mixed phase by a function which depends on the broadening of the phase transition [30, 131, 139]. Therefore, the Gibbs construction for the mixed phase can be simulated through a continuous interpolation between the hadronic and quark phases to mimic the mixing or percolation [30]. Results using a sharp phase transition are shown in section 5.

4 Probing the dense matter EoS II. Astronomical constraints pre GW17081

NSs were first observed serendipitously in the radio wavelengths as pulsars: extremely stable periodic signals that spin-down slowly with high regularity as they loose rotational energy through dipolar radiation. Soon after their discovery, radio pulsars were interpreted as rapidly-rotating highly-magnetized NSs. Currently, $\sim$2000 pulsars are known333See http://www.atnf.csiro.au/research/pulsar/psrcat. Despite their precise temporal behavior, more than 500 glitches in about 180 pulsars have been detected [140]. Those glitches are characterized by sudden changes in the spin period, followed by a gradual relaxation phase. They could be originated during the fast coupling of the spin rates of the crust and the core of the NS. Continuous pulsar monitoring at radio wavelengths can give clues about the NS interior composition through timing analysis of glitching activity [141, 142].

NSs are also observed in the X-ray band. Most of them are found in X-ray binaries (XRBs) where the X-ray radiation is produced when matter from the companion stars is accreted by the NSs. Isolated NSs have also been detected in X-rays, either young NSs in supernova remnants, such as the Crab and Vela pulsars, but also middle-aged close-by NSs, such as the famous Magnificent Seven [143]. Some NSs have also been detected in the UV-optical bands and more recently in the near infrared [144].

The information gathered from pulsar timing of glitches and X-ray spectral analysis of NSs showing thermal cooling suggests that a hadronic crust must be present in at least the vast majority of the NSs observed. Those crusts are essentially composed by neutron-rich nuclei forming an ordered lattice, and a gas of free neutrons and electrons in beta equilibrium. NSs in X-ray binaries are able to accumulate new, light material on their hard surfaces that burns through unstable thermonuclear explosions known as X-ray bursts. Detailed temporal and spectral studies of the X-ray emission of NSs can lead to an estimation of their masses, radii, moments of inertia and internal composition. Several more involved reviews focused on the observational achievements and issues encountered can be found in Refs. [145, 146, 147, 148, 149]. We now describe some of the most important astrophysical results previous to the discovery of GW170817.

Besides the two NSs involved in GW170817, there are about 40 NSs with estimated masses444See https://stellarcollapse.org/nsmasses. Some of them are extremely precisely determined from their orbital parameters. NSs masses tend to concentrate around the canonical 1.4 M⊙ but the most stringent observational constraints on the NS EoS arise from two pulsating NSs in binary systems with white dwarfs companions: PSR J1614–2230 [13] with $M=1.908\pm 0.016$ M⊙ [16], and PSR J0348+0432, with $M=2.01\pm 0.04$ M⊙ [15]. These highly precise measurements based on timing analysis, challenge the high-density EoS models, but also impose limits on the EoS at lower densities (see figure 6 in section 5).

On the other hand, NSs spin can be precisely determined either from pulsar radio or X-ray timing observations. A theoretical NS maximum spin frequency is given by the Kepler limit. Thus, spin measurements might also be used to constrain the EoS. Unfortunately, the maximum observed spin frequency of 716 Hz for a binary radio pulsar in Terzan 5 [150] is not enough to stablish strong constraints on the NSs EoS [151], but this situation would change if a faster-spinning NS is discovered.

In other eclipsing binary systems containing pulsars with low-mass stellar companions, higher masses have been estimated, but with significantly larger systematics errors. In these systems, pulsars strongly irradiate their stellar companions, leading to less reliable mass determinations. For instance, a mass of $M=2.40\pm 0.12$ M⊙ has been estimated for PSR B1757+20 [152] and, more recently, a mass of $M=2.27\pm 0.17$ M⊙ for PSR J2215+5135 [153]. However, uncertainties in the inclination, make the NS mass also compatible with lower values. Independent measurements of the orbital inclination of these systems could lead to more precise mass determinations [154].

Although maximum-mass measurements alone can provide constraints on the nucleon interactions [155], the current limits at $\sim 2$ M⊙ do not rule out the appearance of exotic particles in the core of NSs [156, 157]. Radius determinations are needed to provide stronger constraints on their EoS. However, a very precise simultaneous determination of the mass and radius of the same NS is still pending. Several approaches based on X-ray observations were developed for this purpose, from the study of thermal emission of isolated NSs to accreting NSs in quiescent state and thermonuclear X-ray bursts in low-mass XRBs. In general, radii observations are challenging and remain strongly model dependent. Measurements are affected by external uncertainties on the interstellar extinction and source distance, but also by intrinsic causes as anisotropies in the emission, magnetic fields, rotation and specially atmosphere composition and residual accretion in XRBs (see Ref. [149] for a thorough analysis of the different issues). Depending on the type of system, different (dis)advantages are found.

4.1 Isolated NSs

Isolated NSs, by definition, do not experience accretion. The X-ray thermal emission coming from these sources originates in their internal cooling. In these cases, radius measurements are extremely difficult due to distance indetermination and poorly known atmosphere compositions. Moreover, in these sources it is quite hard to properly estimate their masses. NSs are born in violent supernovae explosions with very high initial temperatures. Isolated NSs experience long-term cooling throughout their lives. This cooling directly depends on their internal composition, being mainly governed by neutrino emissivities and heat-transport coefficients. If present, superfluidity can strongly modify the cooling characteristic timescales (see Ref. [158] for a detailed review). NS cooling calculations present many uncertainties due to several not-fully-understood processes, but the observation of an enhanced cooling for about ten years in the young NS in Cas A [159], may give direct information on its internal composition. Despite this result having been questioned [160, 161], if this rapid cooling is real, it could be explained by a transition to a neutron superfluidity phase in the NS core [162], which would become the first evidence of superfluidity and superconductivity at supra-nuclear densities within a NS. In a recent paper, the authors of Ref. [163] considered the cooling of several isolated NSs and found that, for a nuclear EoS, the threshold density for direct Urca processes is reached in NSs with masses in the 1.7–2.0 M⊙ range. However, 2 M⊙ NSs are less likely to cool down by the direct Urca process due to neutron-triplet superfluidity suppression.

4.2 Quiescent emission from NS cooling

NSs in low-mass XRBs tend to be discovered by wide-field X-ray monitors during transient episodes of strong accretion (i.e. the outburst state). In those episodes, their luminosity is highly increased and the NS crust significantly heated. When accretion stops, or strongly reduces, the NS enters in the so-called quiescent state. During quiescence, the emission is dominated by thermal radiation originating at the NS surface, which can be studied using pointed X-ray telescopes like Chandra, XMM-Newton and NuSTAR. If the distance to the source is known, a direct estimation of the NS radius can be achieved by fitting a NS atmosphere model to the time-evolving spectra. By means of a systematic study of five quiescent sources, it was shown in Ref. [164] that, depending on the assumptions made on the NS atmosphere composition, the EoS of dense matter is either incompatible with a strong phase transition, or this possibility is allowed. They also predict that the most probable radius for a canonical NS of 1.4 M⊙ should be in the 10.5–12.7 km range. Discrepancies arising from chemical composition uncertainties are well illustrated by two studies on NGC 6397: for a hydrogen atmosphere, the authors of Ref. [165] estimated $R=9.4\pm 1.2$ km, while for helium, a radius larger than 2 km was found [166]. More recently, the authors of Ref. [167] added several corrections like those resulting from the inclusion of NS rotation, and, based on a comprehensive analysis of 12 sources, they estimate the radius of a 1.5 M⊙ NS to be in the range 10.1–11.1 Km.

4.3 Thermonuclear X-ray bursts in accreting binary systems

NSs in low-mass XRBs show episodes of unstable thermonuclear burning of fresh hydrogen and helium deposited on their surfaces by the accretion processes [168, 169], leading to bright flares called Type I X-ray bursts. Some of these bursts are strong enough to reach the Eddington luminosity. In those cases, the radiation pressure triggers the expansion of the outer layers of the NS atmosphere, leading to a Photospheric Radius Expansion (PRE) burst. During these events, the luminosity increases by a factor of $\sim$10 in about 0.5–5 s. After that quick rise, the flux decreases exponentially for 10–100 s down to the persistent level. After such fast expansion, the atmosphere contracts again until it reaches back the NS surface, at touchdown, starting the so-called cooling phase. If the cooling process is homogeneous throughout the NS surface, the radius of the NS can be estimated by fitting the time-resolved spectra during the cooling phase of the bursts. Moreover, assuming that the Eddington luminosity is reached at the touchdown point, the mass and the distance to the NS can also be estimated if the chemical composition is known [170].

Although the method seems very powerful, the strength of the attainable constraints is limited by several uncertainties associated with some assumptions: distance and atmospheric composition have a strong impact on the results [168, 164], and it is not fully clear that touchdown is completely reached in all PRE bursts neither if that the cooling phase strictly follows the expected behavior [171, 172]. Moreover, the presence of residual accretion affecting the spectral fits cannot be fully discarded. Following this approach and considering a carefully selected set of sources, the authors of Ref. [173] found that the radius of an 1.4 M⊙ star is favored in the 9.8–11.0 km range. However, despite statistically significant constraints are often found using this method, results may suffer of systematics that should be evaluated with caution [149].

A specific spectral model for the cooling phase of NS X-ray bursts was developed to constrain the mass and the radius of the NS in 4U 1724–307 assuming different chemical compositions [174]. The method was applied later to three sources showing X-ray bursts at low mass-accretion rates, leading to radii estimations in the 10.5–12.8 km range for a canonical 1.4 M⊙ NS [175]. More recently, the authors of Ref. [176] fitted the same spectra using atmospheric burst models and obtained preferred radii and masses of 12.4$\pm$0.4 km and 1.9$\pm$0.3 M⊙, respectively. This method allows to simultaneously fit constraints to the atmosphere composition.

4.4 Emission lines

The detection of an atomic-line feature at the surface of a NS would serve to directly estimate its gravitational redshift, and thus the ratio between it mass and it radius. The authors of Ref. [177] have reported a narrow absorption feature in EXO 0748–676, but more recent observations failed to confirm the detection. The spinning of NSs should widen emission lines, making them much harder to detect, but recent studies suggest that line profiles from rotating NSs might actually be narrower than the initially predicted [178].

4.5 Glitch activity

In two recent works the glitch population has been used to probe the masses of NSs. The standard glitch scenario, assuming that the entire excess of angular momentum in the core is transferred to the crust during the largest glitch detected in a given pulsar, the authors of Ref. [179] placed limits on the NS mass. For example, for Vela and PSR J0537–6910, they found masses of 1.35$\pm$0.08 M⊙ and 1.25$\pm$0.06 M⊙, respectively. A different approach was suggested by the authors of Ref. [180], who performed a detailed thermal evolution model using a nuclear EoS with superfluidity. For the same two pulsars, the later authors obtained 1.51$\pm$0.04 M⊙ and 1.83$\pm$0.04 M⊙.

In turn, the highly-magnetized AXP 1E 2259+586 has experienced a so-called anti-glitch, or a sudden spin-down [181]. This glitch was accompanied by an X-ray burst. As its behavior cannot be accounted by the standard glitch theory, different models assuming an internal [182, 183, 184] or external [185] origin have been proposed. More recently, new anti-glitches have been reported in the accreting pulsar NGC 300 ULX-1 [186]. Despite the present uncertainties, a precise detection of a glitch in a pulsar with known mass in a binary system would substantially increase our understanding of NS interiors.

5 Post GW170817 era and restrictions to the NS equation of state

Every time a new electromagnetic (EM) window has been opened a revolution has occurred. Seeing in infrared, radio, x-rays and $\gamma$-rays has decisively changed our understanding of the Universe. Three years ago a completely different window has opened up: we started to see in gravitational waves. This new ability allows us to study the most energetic events in the Universe from a radically different point of view and to detect objects we had not even imagined. Gravitational waves are not absorbed or reflected by matter like electromagnetic radiation, for this reason through their detection it is possible to see objects without the EM problems Astronomy has. Moreover, as the Universe become transparent to GWs much earlier than to EM waves, their detection and study would help us to understand the first instants after the Big Bang.

Until 2015, GWs had been only indirectly detected and the most important of these indirect observations was associated to the Hulse-Taylor pulsar [187, 188, 189]. On September 14th, 2015 at 9:50:45 UTC, the LIGO detectors in Hanford, Washington and Livingston, Louisiana, detected GWs emitted during a binary black hole merger [190]. This detection marked the beginning of a new era: the gravitational-wave astronomy era. Analysis of the signal of GW150914 that lasted $\sim 0.2$s has confirmed the existence of both gravitational waves predicted by Einstein and of solar-mass black holes. For the first time in history, General Relativity was tested in the strong-field regime, and no deviation from theoretical predictions was found [191].

After decades of development, LIGO and Virgo collaborations had achieved the direct detection of GWs emitted from ten binary black hole mergers: GW150914, GW151012, GW151226, GW170104, GW170608, GW170814, GW170729, GW170809, GW170818 and GW170823 (see Ref. [192] and references therein). By getting some adjustments and improvements in the detectors sensitivity, the strongest evidence to confirm the hypothesis that mergers of binary stars are the cause of short gamma-ray bursts would appear in 2017: the detection of the GW170817 and GRB 170817A counterpart from a merger of two NSs.

5.1 The GW170817 and GRB 170817A connection

The LIGO/Virgo Collaboration detected for the first time GWs emitted from a binary NS merger GW170817 [193]. On August 17, 2017, at 12: 41: 04.4 UTC the fusion of the two NSs was registered. The first electromagnetic signal detected was GRB 170817A, a short gamma-ray burst, detected by Fermi and INTEGRAL 1.74 $\pm$ 0.05 seconds after the fusion time and which lasted approximately 2 seconds. The compact objects that merged have been determined to be compatible with low-mass NSs with masses in the range $\sim 1.2$ to $\sim 1.6M_{\odot}$. The nature of the $\sim 2.7M_{\odot}$ remnant is not clear; it could be either an extremely massive NS or a low-mass BH [194].

The combined data obtained from GW170817 and GRB 170817A unambiguously associates binary NS mergers with short gamma-ray bursts (sGRB) and helped us to understand the nature and properties of these systems, constituting central engines of sGRB [195]. The NS-NS merger took place at a distance of $\sim 40$Mpc form the Earth. The analysis of the data from GW170817 and GRB 170817A sheds some light on some fundamental physical aspects. It allows us to constrain the difference between the speed of light in vacuum, $c$, and the speed of gravity, $c_{g}$, a difference determined to be between $-3\times 10^{-15}<(1-c_{g}/c)<7\times 10^{-16}$ [196]. Moreover, it served to put limits on the violation of Lorentz invariance and test the equivalence principle.

5.2 Mergers of NSs and restriction to the EoS from GWs

In a NS-NS merger simulation , linearized gravity can not be used to describe the system properly. In such events, the dynamical evolution is critically influenced by the backreaction of the GWs emitted. This kind of situation is studied in the “post-Newtonian” (PN) approximation. Within this scenario, backreaction is considered when the 2.5PN order is analyzed (for details, see Ref. [197]). PN expansions do not converge when the velocities involved are a considerable fraction of the speed of light in vacuum. To overcome this problem different proposals of resummation of series have been developed. One commonly used in the literature is the so called Effective One Body model, first developed for binary BH mergers [198] and extended later to calculate waveforms produced in binary NS mergers using full 3D numerical relativity [199]. Simulations solving the dynamics of mergers in full 3D general relativity spread rapidly also thanks to free access codes like Einstein Toolkit [200] and LORENE (see, https://lorene.obspm.fr/). Final stages of the merger and the post-merger phase must be studied using numerical relativity. A review of the initial data of binary NSs mergers and advanced simulations with approximate treatments of gravity can be found in Ref. [201].

Binary NS merger simulations in full 3D general relativity were performed for the first time in the early 2000’s [202, 203]. In these pioneering works, mergers of polytropic NSs with equal masses were studied. Since then, huge advances have been made in the area of numerical modeling of binary NS mergers. Many different EoS for NS in merger simulations have been used and a range of mass and spins of the components of the binary system has been covered to produce a comprehensive catalogue of waveforms that can be compared with observable data (for details, see Ref. [204] and references therein).

GW waveforms produced during the early stages of a NS merger are characterized almost unambiguously by the tidal deformability $\Lambda$, which is a measure of how a NS is deformed due to the gravitational field of the companion compact object (see, for example, Ref. [205]). A qualitative study of the astrophysical implications of measuring the pre-merger waveform in NS-NS or NS-BH mergers is performed in Ref. [206]. In this paper, the authors conclude that, as waveforms deviate from the pointlike particle approximation, this information would be a powerful tool to help unravel the nature of compact objects EoS, helping us understand the behavior of matter when it is compress to extreme NS-like pressures.

Using linearized theory, the tidal deformability can be written in terms of the dimensionless tidal Love number $k_{2}$ (see, Ref. [207], and references therein). Recently, the authors of Ref. [208] have performed calculations using hybrid EoS within the “Constant-Speed-of-Sound” (CSS) parametrization, a generic parametrization of the quark matter EoS with a density-independent speed of sound [209], to show how tidal deformability behaves in the presence of sharp phase transitions.

Analysis of data from GW170817 put an upper limit on the effective tidal deformability $\tilde{\Lambda}$ (i.e. the mass weighted-average of the individual tidal deformabilities) as well as the chirp mass of the NS-NS binary system. The $\tilde{\Lambda}<800$ limit allowed to conclude, with $90\%$ confidence, that the radius of a NS with a mass of $1.4M_{\odot}$, $R_{1.4}$, cannot exceed $\sim 13.6$ km [210]. Other studies (see, for example, Refs. [24, 27, 211, 212, 213]) arrive to similar conclusions. This result can be used to discard several hadronic EoS, like the famous NL3 and GM3 parametrizations of the Walecka RMF as shown in figure 6. In addition, it has been argued that the radius of the non-rotating maximum-mass configuration must be larger than $\sim 9.6$ km [214]. The improved NS EoSs obtained with DDRMF parametrizations are compatible with the constraint of $R_{1.4}$ derived from the GW170817 data (for details, see figure 7).

After a detailed analysis of data from GW170817 and its EM counterpart, additional restrictions to the EoS have been found. According to Ref. [215], the pressure at two times saturation density must be in the range $P(2\,n_{0})=3.5^{+2.7}_{-1.7}\times 10^{34}{\rm dyn/cm}^{2}$, and at six times saturation density in the range $P(6\,n_{0})=9.0^{+7.9}_{-2.6}\times 10^{35}{\rm dyn/cm}^{2}$ at the $90\%$ confidence level. The first restriction is not very stringent, but the second one serves to discard very soft hadronic EoSs.

Multiwavelength (X-ray, UV, optical, infrared, and radio bands) analyses of the EM counterpart of GW170817 (called AT2017gfo) tentatively constraint $\tilde{\Lambda}>400$ [216]. These constraints serve to exclude hadronic EoSs that are either too stiff or too soft. Some criticisms have been raised up regarding the results presented in Ref. [216]. The authors of Ref. [217] suggest that $\tilde{\Lambda}>400$ restriction might have been an overestimation. The authors of Ref. [218] have studied the role of vector-isovector meson cross interactions and their effect on a series of hadronic EoSs and on macroscopic quantities, in particular the M-R relationship and the tidal deformability. As a general result, they found the additional interaction reduces the radius of the canonical NS, allowing all the presented models to satisfy the restrictions imposed by GW170817.

Several model-dependent studies have been performed analyzing the existence of HS. Pure quark matter cores in HS were analyzed using the Maxwell construction for the hadron-quark phase transition in Ref. [219] Alternatively, several works indicate that HS might exist but without pure quark matter in their cores, rather a hadron-quark mixed-phase core (see, for example, Refs. [220, 221]). Moreover, the implications of the radial constraints on a $1.4M_{\odot}$ NS suggest that HS with quark cores can exist only if the hadron-quark phase transition occurs at an energy density similar to that of nuclear saturation [222].

Another interesting series of works are those dedicated to finding GWs from the post-merger phase. In Ref. [223] the relevance of the information of detecting GWs from the post-merger phase is emphasized, specially the connection of $f_{peak}$, the dominant post-merger frequency which contains information about the structure and composition of the compact remnant [224, 225]. In addition, it has been found that a strong discontinuity in $f_{peak}$ as a function of the total mass would be an indicator of the coexistence of two families of compact stars: hadronic and quark matter stars [226]. On the other hand, the authors of Refs.[227, 228] have analyzed a large number of accurate numerical-relativity simulations of binaries with nuclear EoSs finding an universal correlation between the low-frequency peak related to the merger process and the total compactness of the stars in the binary system. This also provides a powerful tool to set constraints on the EoS.

A search of GWs emitted in the post-merger phase of GW170817 has been performed, but no evidence of emission has been found yet [229].

6 GWs from isolated NSs

The most promising astrophysical scenario in which isolated NSs could radiate gravitational waves is probably related to the birth of proto-NSs in core-collapse supernovae due to the rate of events and the energy involved (see Ref. [230] and references therein). Additionally, isolated NSs could radiate GWs if the stars are deformed and rotating. Other circumstances that could lead to emission of GWs from isolated NSs are pulsar glitches, magneto-elastic oscillations in magnetars, and “quasi-normal” mode oscillations excited by the rearrangement of a compact star after a hadron-quark phase transition or a star quake.

It is expected that future LIGO/Virgo observational runs might be able to detect not only GWs emitted during mergers of compact objects but also those emitted by isolated NS.

The EoS of NS matter could be constrained by the analysis of the GW emission of stellar pulsation modes. Such oscillations could be excited in newly born proto-NSs [231] due to starquakes and glitches [232], by accretion in a binary system, or by the rearrangement of the star following the conversion of a hadronic star into a quark star [233, 234]. It can also occur in the post-merger phase of a binary NS system if the remnant does not collapse into a black hole [226, 29, 235].

In this Section we present some basic aspects of NSs oscillations and their link to GW emission from isolated NSs. The focus will be put on non-radial oscillations, in particular on the $g-$mode whose detection could be a strong indicator of the occurrence of a sharp phase transition from hadronic matter to quark matter in the interior of NSs. We also briefly discuss the role of dynamical and secular instabilities in NSs.

6.1 Quasi-normal modes of NSs

In a pulsating star, energy is lost by the emission of gravitational radiation; such vibrations are not described by normal modes but by “quasi-normal” modes. The pulsation equations are obtained by introducing small perturbations in the Einstein equations [236, 29]. A general linear perturbation of a scalar quantity (e.g. $\delta\epsilon$ in the energy density or $\delta p$ in the pressure) can be written as a sum of quasi-normal modes that are characterized by the indices $(\ell,m)$ of the spherical harmonic functions $Y_{\ell}^{m}$ and a time dependence of the form $e^{i\omega_{\ell}t}$, where $\omega_{l}$ is a complex frequency (the real part represents the actual oscillation frequency and the imaginary part the damping produced by the emission of GWs). The perturbations in the fluid’s four-velocity $\delta u_{\mu}$ can be expressed in terms of vector harmonics, and the metric perturbation $\delta g_{\mu\nu}$ in terms of spherical, vector, and tensor harmonics. These perturbations can be classified as polar or axial, according to their parity (which is defined as the change in sign under a combination of a reflection in the equatorial plane and rotation by $\pi$). Polar perturbations have parity $(-1)^{\ell}$ and axial perturbations have parity $(-1)^{\ell+1}$. The important point here is that, for non-rotating relativistic stars, it can be shown that the polar and axial perturbations are completely decoupled (see Ref. [237] and references therein). Polar modes produce spheroidal fluid deformations, and axial modes toroidal deformations. Modes are said to be fluid if they involve a significant fluid motion, and spacetime if they only involve vibrations of spacetime. Polar fluid modes are analogous to the Newtonian fluid pulsations. The most well known in the literature are the fundamental ($f$), the pressure ($p$), and the gravity ($g$) modes, but there are other relevant modes such as the shear mode in a solid crust or the Alfvén modes in the presence of magnetic fields. In the case of axial modes, there are torsional modes and rotational modes (Rossby and inertial modes). There are also polar and axial spacetime modes, which only involve vibrations of spacetime (the so called $w$ modes).

It is important to note at this point that the global properties of NSs, such as masses and radii, are only marginally modified by rotation as long as the stars’ spin frequencies are well below the Kepler frequency. With the exception of pulsar J1748-244ad, which rotates at 716 Hz, this is the case for all NSs known to date. Therefore, restricting ourselves to non-rotating compact objects, polar non-radial perturbations introduce changes in the Schwarszchild line element that can be written as

[TABLE]

The Lagrangian displacements that perturb the fluid can be written as

[TABLE]

The functions $H_{1}$, $W$, $K$ and $X$ (a linear combination of $V$, $W$ and $H_{0}$) satisfy a system of first order differential equations. Outside the compact object, fluid perturbations vanish and the perturbations are ruled by the Zerilli equation. For a given $\ell$, the quasi-normal modes of the star are solutions of the pulsation equations that behave as pure outgoing waves at infinity, are regular at the NS center, and verify that the interior solution matches continuously with the exterior perturbation at the NS surface (see Ref. [238] and references therein).

6.2 Exotic stars, HSs and oscillation modes

Stellar pulsation modes may emit gravitational radiation that could be detected by the LIGO/Virgo observatories as well as by a new generation of planned GW detectors such as KAGRA, the Einstein Telescope and IndIGO (see, for example, Ref. [239] and references therein). Particularly, the identification of features in the $f$ and $g$ modes may signal unequivocally the appearance of deconfined quark matter in NSs or reveal the existence of exotic stars.

The energy of the fundamental fluid mode ($f$) is quite effectively channeled into GW emission. Its properties have been extensively analyzed for several EoSs. An interesting and observationally useful property of the $f$-mode is that it can be described by universal fitting formula that are argued to be quite EoS-independent. In fact, it has been shown that the oscillation frequency $f$ of the fundamental mode has a reasonably linear dependence on the square root of the average density, and the damping time $\tau$ can be fit with a simple formula involving the stellar mass and radius [240, 28, 241, 242]. These results suggest that the mass and radius of a compact object could be inferred if $f$ and $\tau$ were detected by the new generation of gravitational wave detectors. In a recent work, the properties of the fundamental mode have been systematically studied using a set of EoS connecting state-of-the-art calculations at low and high densities [238]. Specifically, a low density model based on the chiral effective field theory (EFT) and high density results based on perturbative Quantum Chromodynamics (QCD) are matched through different interpolating polytropes fulfilling thermodynamic stability and subluminality of the speed of sound, together with the additional requirement that the equations of state support a two solar mass NS. Using this family of EoSs, the frequency and the damping time of the $f$-mode are constrained within narrow, quite model-independent windows (see figure 8) which are in good agreement with many previous works that used a variety of phenomenological EoSs. In principle, GWs arising from the $f$-mode of non-rotating compact stars should fall within the window presented in figure 8. If not, the emitting object could be some kind of exotic star (see, for example, Ref. [243]).

The previous results hold for non-rotating objects. However, it is expected that all NSs should have at least some degree of rotation. Moreover, promising scenarios for detecting the $f-$mode such as NS mergers are likely to involve rapidly spinning NSs. In spite of its relevance, a generalization of GW asteroseismology to rotating objects has only been carried out recently within the Cowling approximation (where the spacetime metric is unperturbed and only the fluid perturbations are considered). The $f$-mode frequency as measured in the stellar rotating frame can be written as a quadratic polynomial in terms of the angular frequency $\Omega$ [244, 245]. For objects spinning at a significant fraction of the mass shedding angular frequency, the $f$-mode oscillation frequency is very different from the non-rotating value. Moreover, when rotation is sufficiently fast, NSs can be destabilized by the so-called Chandrasekhar-Friedman-Schutz instability [246, 247], and some non-axisymmetric modes may grow in amplitude rather than decaying under the emission of GWs (see Section 6.4).

Another important family of fluid “quasi-normal” modes are the gravity modes ($g$) which arise because gravity tends to smooth out material inhomogeneities along equipotential level-surfaces. They appear only in the presence of a significant temperature or composition gradient or as a consequence of sharp discontinuities due to first order phase transitions [29, 248, 249, 250]. Unlike pressure ($p$) modes where the radial displacement is much bigger than the tangential one, for $g$-modes the tangential displacement is much bigger than the radial one.

If density discontinuities actually exist inside NSs due to the existence of a deconfined quark matter core, the associated $g$-mode would be a valuable indicator of its presence. Moreover, as shown in Refs. [251, 252, 253, 250, 30] the $g$-mode frequency strongly depends on the amplitude of the discontinuity, as it is shown in figure 9. Thus, it could be used to infer the value of the density jump and shed some light on the properties of the hadron-quark transition.

NJL-like models typically yield equations of state that feature a first-order transition to quark matter. These kinds of models have been extensively used to describe the matter in the inner cores of HSs [254, 104, 105, 107, 106]. The first order deconfinement transition of hadronic matter at low temperature but high density can be modeled through a hybrid EoS, which describes NS matter at low densities in terms of hadrons and of deconfined quarks at high densities, assuming a sharp first-order phase transition between hadronic and quark matter and considering the CSS parametrization. Thus, supposing the hadron-quark surface tension is high enough to disfavor mixed phases, and restricting EoSs to those that allow NS to reach 2$M_{\odot}$, it was found that the appearance of the quark matter core either destabilizes the star immediately (this is typical for non-local NJL models) or leads to a very short HS branch in the mass-radius relation (this is typical for local NJL models). By using the CSS parametrization it has been shown that the reason for the near-absence of HSs is that the transition pressure is fairly high and the transition is strongly first-order [106]. This situation changes if the possibility of formation of a color superconducting phase is considered in the description of dense matter. Color superconductivity lowers the transition pressure, making it possible to construct stable hydrid stars containing a diquark matter core. According to the results presented in Ref. [30], the maximum mass star corresponding to the different families constructed through the hybrid EoS with diquarks would be composed of 60-75 $\%$ quark matter in a superconducting phase.

In figure 9 the frequencies of the $g$-modes have been obtained for the different hybrid EoSs as a function of the quotient between the quark-hadron energy gap, $\Delta\epsilon$, and energy density at the phase transition, $\epsilon_{\rm trans}$, one of the parameters of the CSS parametrization. This result suggests that the obtained relationship is of universal nature and could be used to better our understanding of the behavior of the dense matter in NSs. After a detection of a $g$-mode, it may be possible to infer the existence of a sharp phase transition inside NSs and to deduce a value for one of the CSS parameters, constraining the EoS.

Figure 10 shows the frequencies of the $g$-modes as a function of gravitational mass for compact objects constructed using the hybrid EoSs presented in Ref. [30].

As it has been mentioned in subsection 5.2, some of the features of the postmerger evolution could reveal important information about the properties of the EoS of compact objects. In particular, the study of the postmerger phase could provide evidence of the existence of a sharp phase transition in the interiors of NSs. One of the key issues would be the presence of a discontinuity in the dominant post merger GW frequency, $f_{peak}$. The imprint on GWs emitted in a NS merger due to a discontinuity in the EoS has also been analyzed in Ref. [255]. Using a particular hybrid, temperature-dependent EoS [256], the authors found that the frequency of the dominant post-merger mode suffers a significant modification when compared to the waveforms emitted when normal NSs merge. Performing simulations using full general relativity, observational signatures of the appearance of quark matter in the post-merger high-temperature phase were found. When it occurs, this phase transition leads to a ringdown GW signal that can be distinguished from the one produced if a purely hadronic EoS is used [257]. In this way, the GWs emitted as a consequence of a binary NS merger can be used to probe the behavior of matter at high-temperature and density and to study, from a different point of view, properties of the QCD phase diagram.

6.3 Intensity of GW signals

If the frequencies of the mentioned oscillation modes were detectable with the Advanced LIGO/Virgo detector or Einstein telescope, we could probe the internal composition of NSs. To set constraints on the different theoretical possibilities of EoSs, we need better observation data, which probably comes with the improvement in sensitivity of the third generation detectors.

From the data, and assuming certain parameters, for example that detection requires certain signal-to-noise ratio and the distance to the source, it is possible to estimate the minimum energy required in each mode to be detectable with Advanced LIGO/Virgo or Einstein telescope, considering the respective characteristic frequency. Useful order of magnitude calculations that allow us to estimate the amount of energy that must be channeled in a particular GW mode in order to be detectable with a given detector can be written [29]

[TABLE]

where $S_{n}$ is the noise power spectral density of the detector, that for a frequency of 1kHz is $S_{n}^{1/2}\sim 2\times 10^{-23}{\rm Hz}^{-1/2}$ for Advanced LIGO/Virgo and is supposed to be, in the same band, $S_{n}^{1/2}\sim 10^{-24}{\rm Hz}^{-1/2}$ for the Einstein telescope. $D$ is the distance to the source, $f$ the frequency of the mode, $Q=\pi f\tau$ the quality factor (which introduces a slight dependence with the damping time) and $S/N$ the desired signal-to-noise ratio.

In table 2 we present the minimum energy that should be channeled into a specific oscillation mode in order to be detected. To obtain this estimate, we have used typical frequencies of 2 kHz, 7 kHz and 1 kHz for the excited modes $f$, $p_{1}$ and $g$, respectively [250, 30]. The obtained energetics are not unrealistic; numerical core-collapse simulations estimate that the energy emitted in GWs in such an event should be of the order of $10^{-9}-10^{-8}{\rm M_{\odot}{c}^{2}}$ [258]. A conservative upper limit might be taken from simulations of the collapse of a rapidly rotating NS to form a black hole. Such an event would radiate GWs with a total energy of $\sim 10^{-6}{\rm M_{\odot}{c}^{2}}$ [259].

6.4 Other sources of GWs: dynamical and secular instabilities

The above considerations were focused on non rotating stars. However, rotating compact objects may be affected by nonaxisymmetric dynamical and secular instabilities. The most studied type of rotational dynamical instability is the so called bar-mode instability. The $m=2$ bar-mode instability would be excited when $\beta=T/|W|$ is larger than 0.24 [260, 261], where $T$ the rotational kinetic energy and $|W|$ the gravitational energy. However, differential rotation may change the scenario significantly. NSs with a high degree of differential rotation may be dynamically unstable for $\beta\gtrsim 0.01$ [262, 263]. Additionally, an $m=1$ one-armed spiral instability may become unstable if differential rotation is sufficiently strong [264, 265].

Nevertheless, in practice only very young proto-NSs are affected by differential rotation. In fact, it is expected that a few minutes after a NS is born in a core collapse supernova, rigid rotation sets in due to the presence of viscosity or a sufficiently strong magnetic field [266]. Thus, these dynamical instabilities are not thought to have an impact in most cold catalyzed NSs.

Another class of nonaxisymmetric instabilities are secular instabilities, which require the presence of dissipation due, for example, to viscosity and/or GWs. In this context the so-called $r-$modes have attracted considerable attention. They are associated with large-scale currents in NSs that couple to gravitational radiation and remove energy and angular momentum from the star in the form of GWs [267, 268, 269, 270]. In the absence of viscous dissipation, they are unstable at all rotation frequencies leading to an exponential rise of the $r-$mode amplitude. However, when viscous damping is taken into account the star is stable at low frequencies and instability windows are expected only at high enough frequencies [270, 271]. In addition, the fact that fast-spinning compact stars are observed suggests that nonlinear damping mechanisms are present and limit the exponential grow of $r-$modes that would destroy the NS [272]. If this instability is stopped at a large amplitude, $r-$modes may be a strong and continuous source of GWs in some objects which can be used to probe the NS interior structure.

7 Summary and Future Prospects

In this work we have presented a general overview of some of the hot topics and open issues concerning the behavior of the matter in the dense cores of NSs. Particular emphasis has been put on the EoS associated with such matter and its connection to NS observables constrained by the gravitational wave event GW170817. Models for the EoS have been presented which are based on the relativistic (density dependent) mean-field approximation to model hadronic matter and a nonlocal extension of the three-flavor Nambu-Jona-Lasinio model for quark matter. Assuming that the latter is in the 2SC+s phase, we have presented hybrid EoSs with sharp quark-hadron interface. Stable hybrid stars are found to be composed of 60-75$\%$ of diquark matter. Motivated by the gravitational-wave event GW170817, we show the results obtained for $g$-mode frequencies of oscillating NSs [30]. Since these modes depend strongly on the amplitude of the energy density discontinuity at the quark-hadron interface, $g$-modes may be used to learn about the presence and size of the energy density jump that characterized the quark-hadron interface. An analysis of the $f$ and $g$ modes suggest that they can be described by fits which are rather independent of the nuclear EoS. Aside from constraints derived from $f$ and $g$ modes, GW from NSs mergers impose restrictions on the masses and radii of NSs, which, in turn, sheds light on their interior compositions. In particular, such detections will also be very helpful to determine the sought-after upper mass limit of NSs.

Although LIGO’s second observing Run O2 ended on august of 2017, it is estimated that the third observing Run O3 will start at the beginning of 2019555LIGO/Virgo O3 run started 2019 April 1 at 15:00 UTC. In this new era, data is made publicly available to allow for rapid follow-up of events by the electromagnetic astrophysical community (https://gracedb.ligo.org/).. These new searches for gravitational-waves from stellar-mass coalescing compact binaries will help astrophysicists to provide new constraints on the EoS of dense nuclear matter composing NSs.

In addition, observational efforts are emphasized through several projects to allow for precise determinations of NS masses and radii. Large effective areas combined with high spectral and timing resolutions of the next generation of X-ray instruments such as the already working Neutron Star Interior Composition ExploreR (NICER666https://www.nasa.gov/nicer) and future planned ATHENA777http://sci.esa.int/athena/ will substantially help to improve our current picture. NICER was designed to reach about 5% accuracy on the determination of the mass and the radius of the closely PSR J0437–4715 by directly fitting the phase-resolved shape of its X-ray pulses. In a recent paper, the authors use NICER observations to update the X-ray timing of radio- and $\gamma$-quiet pulsars PSR J1412+7922 and PSR J1849–0001 [273]. Moreover, five X-ray bursts on 4U 1820–30 were already detected with NICER. An analysis of the first of them has been recently published [274]. The high sensitivity of ATHENA will improve radius estimates of isolated NSs, and cooling tails of X-ray bursts, and allow for searches of spectral features in NS atmospheres. At radio wavelengths, the Square Kilometer Array (SKA888https://www.skatelescope.org/) will become the largest radio facility, beginning science operations very soon, allowing for the discovery of a significant number of new pulsars, detect hundreds of glitches, and to perform new NS mass determinations.

Acnowledgements

This research was supported by Universidad Nacional de La Plata and CONICET (Argentina), Grants G140, G157, X824 and PIP-0714, and the National Science Foundation (USA), grant number PHY-1714068. Germán Lugones acknowledges CNPq Brazil for partial financial support.