Neutron star cooling and GW170817 constraint within quark-meson coupling models

TL;DR

This study evaluates five quark-meson coupling models to analyze neutron star cooling and GW170817 constraints, finding correlations between the symmetry energy slope and cooling rates, but no single model fits all observations.

Contribution

It compares different QMC models with various confinement potentials and assesses their ability to satisfy GW170817 constraints and neutron star cooling data.

Findings

No model simultaneously satisfies GW170817 and cooling observations.

Cooling rate correlates with the slope of the symmetry energy.

Pasta phase has negligible impact on results.

Abstract

In the present work we used five different versions of the quark-meson coupling (QMC) model to compute astrophysical quantities related to the GW170817 event and to neutron star cooling process. Two of the models are based on the original bag potential structure and three versions consider a harmonic oscillator potential to confine the quarks. The bag-like models also incorporate the pasta phase used to describe the inner crust of neutron stars. We show that the pasta phase always play a minor or negligible role in all studies. Moreover, while no clear correlation between the models that satisfy the GW170817 constraints and the slope of the symmetry energy is found, a clear correlation is observed between the slope and the fact that the cooling is fast or slow, i.e., fast (slow) cooling is related to higher (lower) values of the slope. We did not find one unique model that can describe, at the same time, GW170817 constraints and give a perfect description of the possible cooling processes.

Tables

Table 1. Table 1: Nuclear matter and stellar properties obtained with the QMC and QMC ω ​ ρ 𝜔 𝜌 \omega\rho models.

Model	$B/A$	$n_0$	${\Lambda }_v$	$g_{\rho }$	$M_N^{\ast }/M_N$	$J$	$L_0$	$K_0$	$M_{max}$	$R_{M_{max}}$	$R_{M=1.4M_{\odot }}$
	(MeV)	(fm-3)				(MeV)	(MeV)	(MeV)	(M⊙)	(km)	(km)
QMC	-16.4	0.15	0.00	8.6510	0.77	34.50	90.00	295	2.14	11.51	13.55
QMC$\omega \rho$	-16.4	0.15	0.03	9.0078	0.77	30.92	69.17	295	2.07	10.96	12.83

Table 2. Table 2: Nuclear matter and stellar properties obtained from the parametrizations of the MQMC model. The free parameters are also given. G σ q ​ 2 superscript subscript 𝐺 𝜎 𝑞 2 G_{\sigma}^{q2} , G ω 2 superscript subscript 𝐺 𝜔 2 G_{\omega}^{2} and G ρ 2 superscript subscript 𝐺 𝜌 2 G_{\rho}^{2} are given in 10 − 5 superscript 10 5 10^{-5} MeV -2 . For all parametrizations, one has B / A = − 16.4 𝐵 𝐴 16.4 B/A=-16.4  MeV and n 0 = 0.15 subscript 𝑛 0 0.15 n_{0}=0.15  fm -3 .

Model	$m_q$	$K_0$	$M_N^{\ast }/M_N$	$J$	$L_0$	$a$	$V_0$	$G_{\sigma }^{q2}$	$G_{\omega }^2$	$G_{\rho }^2$	$M_{max}$	$R_{M_{max}}$	$R_{M=1.4M_{\odot }}$
	(MeV)	(MeV)	-	(MeV)	(MeV)	(fm-3)	(MeV)				($M_{\odot }$)	(km)	(km)
MQMC1	$210.61$	$295$	$0.84$	$34.50$	$93.20$	$0.95$	$-92.27$	$5.13$	$8.33$	$14.67$	$1.97$	$11.43$	$13.55$
MQMC2	$210.61$	$295$	$0.84$	$30.92$	$82.46$	$0.95$	$-92.27$	$5.13$	$8.33$	$12.18$	$1.97$	$11.34$	$13.32$
MQMC3	$210.61$	$295$	$0.84$	$25.00$	$64.70$	$0.95$	$-92.27$	$5.13$	$8.33$	$8.08$	$1.97$	$11.18$	$12.94$

Equations

$$[ito0.0pt/\text{hss}\partial$$

$${\psi }_{q_N}(\mathbf{r},t)$$

$${ϵ}_{q_N}={\Omega }_{q_N}+R_N\left(g_{\omega }^q\omega +\frac{1}{2}{g}_{\rho }^q{\tau }_z{\rho }_{03}\right),$$

$${\beta }_{q_N}=\sqrt{\frac{{\Omega }_{q_N}-R_N{m}_q^{\ast }}{{\Omega }_{q_N}+R_N{m}_q^{\ast }}},$$

$${\cal N}_{q_{N}}^{-2}=2R_{N}^{3}j_{0}^{2}(x_{q})\left[\Omega_{q}(\Omega_{q}-1)+R_{N}m_{q}^{*}/2\right]\Big{/}x_{q}^{2}~{},$$

$$j_{0_N}(x_{q_N})={\beta }_{q_N}{j}_{1_N}(x_{q_N}).$$

$$E_N^{\mathrm{bag}}=\sum _q{n}_q\frac{{\Omega }_{q_N}}{R_N}-\frac{Z_N}{R_N}+\frac{4}{3}\pi R_N^3{B}_N,$$

$$\frac{d{M}_N^{\ast }}{d{R}_N^{\ast }}=0,N=p,n.$$

$$\epsilon =\frac{1}{2}{m}_{\sigma }^2\sigma +\frac{1}{2}{m}_{\omega }^2{\omega }_0^2+\frac{1}{2}{m}_{\rho }^2{\rho }_{03}^2+3{\Lambda }_v{g}_{\omega }^2{g}_{\rho }^2{\omega }_0^2{\rho }_{03}^2$$

$$+\sum _N\frac{1}{{\pi }^2}{\int }_0^{k_N}{k}^{2}dk[k^2+M_N^{\ast 2}{]}^{1/2},$$

$$p=-\frac{1}{2}{m}_{\sigma }^2\sigma +\frac{1}{2}{m}_{\omega }^2{\omega }_0^2+\frac{1}{2}{m}_{\rho }^2{\rho }_{03}^2+{\Lambda }_v{g}_{\omega }^2{g}_{\rho }^2{\omega }_0^2{\rho }_{03}^2$$

$$+\sum _N\frac{1}{{\pi }^2}{\int }_0^{k_N}{k}^{4}dk/[k^2+M_N^{\ast 2}{]}^{1/2}.$$

$${\omega }_0=\frac{g_{\omega }(n_p+n_n)}{m_{\omega }^{{\ast }^2}},{\rho }_{03}=\frac{g_{\rho }(n_p-n_n)}{2m_{\rho }^{{\ast }^2}},$$

$$n_B=n_p+n_n=\sum _N\frac{2k_N^3}{3{\pi }^2},N=p,n.$$

$$m_{\omega }^{{\ast }^2}=m_{\omega }^2+2{\Lambda }_v{g}_{\omega }^2{g}_{\rho }^2{\rho }_{03}^2$$

$$m_{\rho }^{{\ast }^2}=m_{\rho }^2+2{\Lambda }_v{g}_{\omega }^2{g}_{\rho }^2{\omega }_0^2.$$

$$\frac{\partial \epsilon }{\partial \sigma }=0.$$

$${\mu }_p={\mu }_n-{\mu }_e,{\mu }_e={\mu }_{\mu }.$$

$$n_p=n_e+n_{\mu },$$

$$P^I=P^{II},$$

$${\mu }_p^I={\mu }_p^{II},$$

$${\mu }_n^I={\mu }_n^{II},$$

$$f(n_p^I-n_e^I)+(1-f)(n_p^{II}-n_e^{II})=0.$$

$$f=\frac{n_B-n_B^{II}}{n_B^I-n_B^{II}}.$$

$${\epsilon }_{matter}=f{\epsilon }^I+(1-f){\epsilon }^{II}+{\epsilon }_e.$$

$$\epsilon ={\epsilon }_{matter}+{\epsilon }_{surf}+{\epsilon }_{Coul}.$$

$${\epsilon }_{Coul}=\frac{2\alpha }{4^{2/3}}(e^2\pi \Phi {)}^{1/3}{\left[\sigma D(n_p^I-n_p^{II})\right]}^{2/3},$$

$$\displaystyle\Phi=\left\{\begin{array}[]{c}\left(\frac{2-D\alpha^{1-2/D}}{D-2}+\alpha\right)\frac{1}{D+2},\quad D=1,3\\ \frac{\alpha-1-\ln\alpha}{D+2},\quad D=2\end{array}\right.$$

$${\mathcal{L}}_{\text{mbox}MQMC}$$

$$\left\{\vec{\alpha }\cdot \vec{k}+{\gamma }^0\left[m_q-V_{\sigma }+\frac{V(r)}{2}\right]+\frac{V(r)}{2}+V_{\omega }+\frac{{\tau }_z{V}_{\rho }}{2}\right\}{\psi }_q$$

$$={\epsilon }_q{\psi }_q$$

$$\displaystyle\vec{\alpha}\cdot\vec{k}=\left(\begin{array}[]{cc}0&\vec{\sigma}\cdot\vec{k}\\ \vec{\sigma}\cdot\vec{k}&0\end{array}\right),\quad\psi_{q}=\left(\begin{array}[]{c}\varphi\\ \chi\end{array}\right),$$

Full text

I Introduction

The observation of the binary neutron star system GW170817 by the LIGO-VIRGO scientific collaboration and also in the X-ray, ultraviolet, optical, infrared, and radio bands gave rise to the new era of multimessenger astronomy. All these joint observations gave support to the idea that GW170817 was produced when two neutron stars merged in NGC4993. The observation of the electromagnetic counterpart of GW170817 by the Fermi Gamma-ray observatory corroborated that binary neutron star mergers are associated to short gamma ray bursts (sGRB) and kilonova emissions probably powered by the radioactive decay of r-process nuclei synthesized in the ejecta Abbott_2017 . Therefore multimessenger observations are an excellent tool to extract precious information of compact stars under extreme conditions.

In the next years more signal detections will be possible by further upgrades of dubbed A+ and LIGO Voyager LIGOthecnicalReport and by the joint collaboration of KAGRA and LIGO India Aasi2015 ; Acernese_2014 ; PhysRevD.88.043007 . In addition, planned third generation observatories such as Einstein Telescope (ET) and Cosmic Explorer (CE) may bring us the detection of binary neutron stars at cosmological distances PhysRevD.96.084004 ; PhysRevD.96.084039 ; Punturo_2010 .

As can be seen the observation of binary systems in different channels is of outstanding importance in order to establish stronger constraints on neutron star physics. As a matter of fact neutron stars have been, for many decades, objects of intense astrophysical research. Electromagnetic observations were used to establish some limits on the mass and radius of neutron stars. In light of the recent detections, the gravitational wave channel has become a new way to observe NS and also provide new neutron star observables, like the tidal polarizability arising just before the merging.

The dimensionless tidal deformability or tidal polarizability (TP) is related to the induced deformation that a neutron star undergoes by the influence of the tidal field of its neutron star companion in the binary system 1995MNRAS.275..301K . That influence is expected to be detected in the low frequency range of the inspiral stage, where the effect is a small correction in the waveform phase. As each different neutron star composition has a characteristic response to the tidal field, the TP can be used to discriminate between different equations of state (EOS). Therefore, this new and very interesting physics offers an exciting possibility to investigate the neutron star composition in a very clear form.

In addition, thermal evolution studies are a complementary way of probing properties of the neutron star. The investigation of the cooling of compact stars has been proved a promising method for exploring the properties of neutron stars, as it is strongly connected to both micro and macroscopic realms (TSURUTA1965, ; Maxwell1979, ; Page2006, ; Weber2007, ; Negreiros2010, ; 2018ApJ…863..104N, ). There is a wealth of literature on the cooling of neutron stars and related phenomena such as magnetic field (Aguilera2008, ; Pons2009, ; Niebergal2010a, ; Negreiros2017b, ), deconfined quark matter (Horvath1991, ; Blaschke2000a, ; Shovkovy2002, ; Grigorian2005, ; Alford2005a, ), superfluidity (Levenfish1994, ; Schaab1997, ; Alford2005, ; Page2009, ), rotation (Negreiros2012, ; Negreiros2013, ; Negreiros2017, ), among others (Weber2005, ; Alford2005, ; Gusakov2005, ; Negreiros2010, ).

In the present work we present five different versions of the quark-meson coupling (QMC) model equations of state chosen to describe neutron star (NS) matter. The first two are respectively the original QMC model with an underlying bag structure and its counterpart with the inclusion of an interaction between the meson $\omega$ and $\rho$ fields. In both cases, the pasta phase is also considered in the description of the NS inner crust. The other three versions are modified QMC (MQMC), where the parameters are adjusted so that the constituent quarks are confined to a flavor-independent harmonic oscillator potential jpg ; npa .

The above mentioned models are confronted with recent astrophysical constraints, including the ones predicted by the GW170817. We also perform cooling simulations for all models and compare the results with observed data. This study will allow us (in conjunction with the aforementioned studies) a better evaluation of the quality of the underlying microscopic models adopted.

The formalism used for the EOS is shown in section II, the necessary equations to estimate the quantities related to the GW constraints are given in section III and for the cooling process in section IV. We present and discuss our results in section V and in the last section we make our final remarks.

II Equations of State

In this section we present the formalism of the original QMC model, its counterpart with the inclusion of the $\omega\rho$ interaction and in the sequel, the modified QMC model.

II.1 QMC and QMC$\omega\rho$ models

In the QMC model, the nucleon in nuclear medium is assumed to be a static spherical MIT bag in which quarks interact with the scalar ($\sigma$) and vector ($\omega$, $\rho$) fields, and those are treated as classical fields in the mean field approximation (MFA) guichon . The quark field, $\psi_{q_{N}}$, inside the bag then satisfies the equation of motion:

[TABLE]

where $m_{q}$ is the current quark mass, and $g_{\sigma}^{q}$, $g_{\omega}^{q}$ and $g_{\rho}^{q}$ denote the quark-meson coupling constants. The normalized ground state for a quark in the bag is given by

[TABLE]

where

[TABLE]

and,

[TABLE]

with the normalization factor given by

[TABLE]

where $\Omega_{q_{N}}\equiv\sqrt{x_{q_{N}}^{2}+(R_{N}\,m_{q}^{*})^{2}}$, $m_{q}^{*}=m_{q}-g_{\sigma}^{q}\,\sigma$, $R_{N}$ is the bag radius of nucleon $N$ and $\chi_{q}$ is the quark spinor. The bag eigenvalue for nucleon $N$, $x_{q_{N}}$, is determined by the boundary condition at the bag surface

[TABLE]

The energy of a static bag describing nucleon $N$ consisting of three quarks in ground state is expressed as

[TABLE]

where $Z_{N}$ is a parameter which accounts for zero-point motion of nucleon $N$ and $B_{N}$ is the bag constant. The set of parameters used in the present work is determined by enforcing stability of the nucleon (here, the “bag”), much like in alex09 , so there is a single value for proton and neutron masses. The effective mass of a nucleon bag at rest is taken to be $M_{N}^{*}=E_{N}^{\rm bag}.$

The equilibrium condition for the bag is obtained by minimizing the effective mass, $M_{N}^{*}$ with respect to the bag radius

[TABLE]

By fixing the bag radius $R_{N}=0.6$ fm and the bare nucleon mass $M=939$ MeV the unknowns $Z_{N}=4.0050668$ and $B^{1/4}_{N}=210.85$MeV are then obtained. Furthermore, the desired values for the binding energy and saturation density (Table 1), are achieved by setting $g_{\sigma}^{q}=5.9810$, $g_{\omega}=8.9817$, where $g_{\omega}=3g^{q}_{\omega}$ and $g_{\rho}=g^{q}_{\rho}$. The meson masses are $m_{\sigma}=550$ MeV, $m_{\omega}=783$ MeV and $m_{\rho}=770$ MeV. With this parameterization, some of the bulk properties at saturation density, namely, the compressibility, the symmetry energy and the slope of the symmetry energy, are calculated and the values can be seen in Table 1. Other parameter sets are possible, as discussed in Gramspasta ; PandaKrein ; Santos-09 ; panda . The values given in the first line of Table 1 are very close to the most accepted values (see dutra14 ; oertel16 , for instance) and $J$ and $L_{0}$ can be easily controlled by the inclusion of a $\omega-\rho$ interaction, as discussed in EPJA_2014 ; Prafulla_2012 ; Cavagnoli_2011 . The larger the value of this interaction, the lower the values of the symmetry energy and its slope. In the present work, we follow the calculations mentioned in Gramspasta , where an $\omega-\rho$ interaction strength that results in a symmetry energy equal to 22 MeV at 0.1 fm*-3* is included, with a consequent change in the $g_{\rho}$ coupling constant. The new values of the symmetry energy and its slope at saturation are also given in Table 1.

In a relativistic mean field (RMF) approximation, the total energy density of the nuclear matter reads

[TABLE]

and the pressure is,

[TABLE]

The vector mean field $\omega_{0}$ and $\rho_{03}$ are determined through

[TABLE]

where

[TABLE]

is the baryon density, and the effective masses of the meson fields are

[TABLE]

and

[TABLE]

Finally, the mean field $\sigma$ is fixed by imposing that

[TABLE]

Our interest lies on stellar matter in $\beta$-equilibrium conditions, i.e.,

[TABLE]

Charge neutrality requires that

[TABLE]

and these two conditions imply that a free gas of leptons (electrons and muons) have to be added to the energy density and pressure of the system.

• •

Pasta phases

We construct the pasta phases within the QMC and QMC$\omega\rho$ models using the coexisting phases method Maruyama-05 ; Avancini-08 . For a given total density $n_{B}$ the pasta structures are built with different geometrical forms, usually called sphere (bubble), cylinder (tube), and slab, in three, two, and one dimensions, respectively. This is achieved from the Gibbs conditions, that impose that both phases have the same pressure, proton and neutron chemical potentials. For stellar matter, the following equations must be solved simultaneously

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $I$ ($II$) represents the high (low) density region, $n_{p}$ is the global proton density and $f$ is the volume fraction of the phase $I$, that reads

[TABLE]

The hadronic matter energy reads:

[TABLE]

Adding the surface and Coulomb terms to Eq. (25) results in the total energy density:

[TABLE]

Minimizing $\varepsilon_{surf}+\varepsilon_{Coul}$ with respect to the size of the droplet/bubble, cylinder/tube or slabs, we obtain maruyama $\varepsilon_{surf}=2\varepsilon_{Coul}$ where

[TABLE]

with $\alpha=f$ for droplets, tubes and slaps, and $\alpha=1-f$ for tubes and bubbles. $\Phi$ is given by

[TABLE]

$\sigma$ is the surface tension, which measures the energy per area necessary to create a planar interface between the two regions and is calculated using an adapted geometric approach Gramspasta .

Notice that the pasta phase is only present at the low density regions of the neutron stars and in this region muons are not present, although they are present in the EOS that describes the homogeneous region.

II.2 MQMC model

The Lagrangian density of the modified quark-meson coupling (MQMC) model, extended in comparison with its standard form jpg ; npa in order to take into account asymmetric nuclear matter rnm ; hss ; prc16 , is given by

[TABLE]

The antisymmetric field tensors $F_{\mu\nu}$ and $\vec{B}_{\mu\nu}$ are given by $F_{\mu\nu}=\partial_{\nu}\omega_{\mu}-\partial_{\mu}\omega_{\nu}$ and $\vec{B}_{\mu\nu}=\partial_{\nu}\vec{\rho}_{\mu}-\partial_{\mu}\vec{\rho}_{\nu}$.

In the original QMC model guichon introduced previously, quarks interact each other through an MIT-like bag potential and the nucleon is described as a bag composed of three quarks. Each set of three quarks interact with another one through meson exchanges, as already mentioned. The same kind of nucleon interaction is also present in the MQMC model. However, the interaction between quarks inside the nucleon is taken into account via a confining harmonic oscillator potential given by $U(r)=\frac{1}{2}(1+\gamma^{0})V(r)$, with $V(r)=ar^{2}+V_{0}$, instead of the bag-like treatment. The potential intensity and depth are related, respectively, to the $a$ and $V_{0}$ constants.

From the Euler-Lagrange equations applied to the Lagrangian density in Eq. (31), one can write the Dirac equation for the quarks as

[TABLE]

with

[TABLE]

and

[TABLE]

where $\sigma$, $\omega_{0}$ and $\rho_{03}$ are the classical meson fields in the mean-field approximation.

If we write the small component ($\chi$) of the Dirac field in terms of the larger one ($\varphi$), and replace it back into one of the coupled equations for $\chi$ and $\varphi$, we will find

[TABLE]

with

[TABLE]

and

[TABLE]

The Schrödinger equation of the tridimensional harmonic oscillator is recognized in Eq. (39), with the lowest order energy identified in the right-hand side. Such an identification gives rise to

[TABLE]

by taking into account that $\frac{3}{2}\omega=\frac{\varepsilon_{q}^{*}-m^{*}_{q}-V_{0}}{2}$ with $\omega=\sqrt{\frac{a}{\varepsilon_{q}^{*}+m^{*}_{q}}}$, in units of $\hbar=c=1$.

Since the center of mass motion of the composite nucleon in the MQMC model is also bound by the harmonic potential between the quarks, corrections in the nucleon wave function must to be taken into account if we consider the composite state as a translationally invariant one. Here, we follow the procedure used in Ref. npa in order to extract center of mass effects from nucleon observables. Firstly, we obtain the center of mass energy, given by npa

[TABLE]

with

[TABLE]

and

[TABLE]

The effective nucleon mass in the medium as the center of mass corrected energy of the three independent quarks is taken into account is then expressed as npa

[TABLE]

and the mean squared nucleon radius, also corrected for center of mass effects, is written as npa

[TABLE]

The harmonic oscillator parameters $a$ and $V_{0}$ are determined by imposing the vacuum values for $M_{N}^{*}$ and $\left<r_{N}^{2}\right>$. Here we adopt $M_{N}^{*}(n=0)=939$ MeV and $\left<r_{N}^{2}\right>(n=0)=0.8^{2}$ fm2.

The equations of state (EoS) and field equations of the MQMC model are given as in the QMC one, by taking $\Lambda_{v}=0$. More specifically, the energy density and pressure of the MQMC model are given by Eqs. (11) and (12), respectively. The mean fields $\omega_{0}$, $\rho_{03}$ and $\sigma$ are obtained as indicated in Eqs. (13) and (17), all of them with the restriction that the $\Lambda_{v}$ parameter is set equal to zero.

The free parameters $G_{\sigma}^{q2}\equiv(g_{\sigma}^{q}/m_{\sigma})^{2}$ and $G_{\omega}^{2}\equiv(g_{\omega}/m_{\omega})^{2}$ are found by imposing the nuclear matter saturation at $n=n_{0}=0.15$ fm*-3* with a binding energy of $B/A=16$ MeV. Finally, $G_{\rho}^{2}\equiv(g_{\rho}/m_{\rho})^{2}$ is determined by fixing a particular value for $J\equiv\mathcal{S}(n=n_{0})$, with the symmetry energy given by hss ; prc16

[TABLE]

with $k_{F}$ being the Fermi momentum. The input free parameter $m_{q}$ is used to control the incompressibility at the saturation density $K_{0}=K(n=n_{0})$, with $K=9\partial p/\partial n$. Here we restrict the MQMC model to present the same $J$ and $K_{0}$ values as those from the QMC models in Sec. II, see Table 1. Such parametrizations are named as MQMC1 and MQMC2. We also generate a third one, namely, MQMC3 in which one has $J=25$ MeV. Notice that all these parametrizations present $J$ and $K_{0}$ values inside the ranges of $25\,\mbox{MeV}\leqslant J\leqslant 35\,\mbox{MeV}$ dutra14 and $250\,\mbox{MeV}\leqslant K_{0}\leqslant 315\,\mbox{MeV}$ stone , respectively.

For the stellar matter calculations, we proceed as described in Sec. II concerning the $\beta$-equilibrium conditions on the chemical potentials and densities. In particular, the nucleon chemical potentials in the MQMC model are given by,

[TABLE]

with the upper (lower) sign for protons (neutrons). The nuclear matter and stellar properties, along with the free parameters obtained from the parametrizations of the MQMC model, are given in Table. 2.

III GW170817 constraints

The gravitational Love number depends directly on the detailed structure of the neutron star (NS). Therefore the observation of Love numbers can offer us paramount information on the NS composition. In fact this physics has triggered intense research recently PhysRevC.87.015806 ; PhysRevC.98.035804 ; PhysRevC.98.065804 ; PhysRevC.95.015801 .

When one of the neutron stars in a binary system gets close to its companion just before merging, a mass quadrupole develops as a response to the tidal field induced by the companion. This is known as tidal polarizability PhysRevD.80.084035 ; PhysRevD.80.084018 and can be used to constrain neutron star macroscopic properties PhysRevD.77.021502 , which in turn, are obtained from appropriate equations of state (EOS).

In a binary system the induced quadrupole moment $Q_{ij}$ in one neutron star due to the external tidal field ${\cal E}_{ij}$ created by a companion compact object can be written as PhysRevD.77.021502 ,

[TABLE]

where, $\lambda$ is the tidal deformability parameter, which can be expressed in terms of dimensionless $l=2$ quadrupole tidal Love number $k_{2}$ as

[TABLE]

To obtain $k_{2}$, we have to simultaneously solve the TOV equations and find the value of $y$ from the following differential equation

[TABLE]

with its coefficients given by

[TABLE]

and

[TABLE]

where $\varepsilon$ and $p$ are the energy density and pressure profiles inside the star. Then we can compute the Love number $k_{2}$, which is given by

[TABLE]

where $y_{R}=y(r=R)$, $C=M/R$ is the compactness of the star and $R$ is its radius.

The tidal deformability $\Lambda$ (i.e., the dimensionless version of $\lambda$) is connected with the compactness parameter $C$ through

[TABLE]

In the next section, $\Lambda_{1}$ and $\Lambda_{2}$ refer to the the values of each one of the neutron star in the binary system.

III.1 The inner and outer crust effects on the tidal polarizability

It is obviously expected that the crust thickness and constitution affect the second Love number and consequently, the tidal polarizability. The neutron star crust is usually divided into two different parts: the outer crust and the inner crust. In jorge2018 , the authors investigated the impact of the crust by considering simple expressions. For the outer crust, the region where all neutrons are bound to finite nuclei, a crystal lattice calculation that depends on the masses of different nuclei was performed. For the inner crust, where the pasta phase is expected to exist, the authors used a polytropic EOS that interpolates between the homogeneous core and the outer crust. The conclusion of their work is that, for a fixed compactness, the second Love number is sensitive to the inner crust, but as the tidal polarizability scales as the fifth power of the compactness parameter, the overall impact is minor. In the present work, we use the BPS EOS bps1971 for the outer crust and the pasta phase for the inner crust and test their effects on the deformability of the NS. We can see the differences between the total equations of state (outer crust + inner crust + liquid core) presented here and the one used in jorge2018 in Fig.1 One can see that the outer crusts are coincident and the liquid core of all EOS are very similar. Hence, most of the differences reside on the inner crust and around the crust-core transition region. Notice that the pasta phase is present only in the QMC and QMC$\omega\rho$ models.

In the next section, we obtain the results for the second Love number and the tidal deformability and discuss them.

IV Neutron star cooling

We now turn our attention to the thermal evolution of neutron stars whose composition is described by the models discussed in this manuscript. The cooling of neutron stars is governed by the emission of neutrinos from their core, and of photons from the surface. All thermal properties of the neutron star from the neutrino emission to the heat transport, depend on their microscopic composition, as different compositions lead to different processes, or alter the rate at which certain processes take place. Evidently the specific heat and thermal conductivity also depend on the microscopic composition of the star. Furthermore global and macroscopic properties of the star, such as size, spin, crystalline structure (of the crust) also play a major role in the thermal evolution of the star. Cooling studies are thus a superb way of bridging the gap between the micro and macroscopic realms of neutron stars.

In recent years there has been great advances in the observation of thermal properties of compact objects (Beloin2016a, ; SafiHarb2008, ; Zavlin1999, ; Pavlov2002, ; Mereghetti1996, ; Zavlin2007, ; Pavlov2001, ; Gotthelf2002, ; McGowan2004, ; Klochkov2015, ; McGowan2003, ; McGowan2006, ; Possenti1996, ; Halpern1997, ; Pons2002, ; Burwitz2003, ; Kaplan2003, ; Zavlin, ; Ho2015, ).

The equations that describe thermal energy balance and transport inside a spherically symmetric general relativistic star are given in 2006NuPhA.777..497P ; 1999Weber..book ; 1996NuPhA.605..531S , with ($G=c=1$)

[TABLE]

In eqs. (57)-(58) $l$ is the luminosity, $\phi$ is the metric function, $m$ is the mass as a function of the radial distance $r$ and $T$ is the temperature. The equations above also depend on microscopic properties such as the neutrino emissivity $\epsilon_{\nu}$, the specific heat $c_{v}$ and the thermal conductivity $\kappa$. A comprehensive review of the details of such quantities can be found in references (Yakovlev2000, ).

One also needs the appropriate boundary conditions for a complete solution of eqs. (57)-(58) and these are given by the heat flow at the center of the star, which obviously needs to vanish, as well as surface conditions which connect the stellar surface luminosity to the heat flowing from the mantle to the star atmosphere. The atmosphere conditions may depend on several circumstances such as, surface magnetic field and/or heavier elements due to accreted matter at the neutron star formation (Gudmundsson1982, ; Gudmundsson1983, ; Page2006, ). In this work we consider the traditional scenario, as described in the aforementioned references, and only briefly consider the possibility of accreted matter in the atmosphere.

In the next section, we numerically solve eqs. (57)-(58) for a set of stars of each model, spanning low to high mass neutron stars, and then confront our results to observed data and draw our conclusions.

V Results and Discussion

In the present section we study the effects of the different versions of the QMC and MQMC models presented earlier on quantities related to gravitational wave observables and on the neutron star cooling.

V.1 GW170817 constraint results

We next show the results of the tidal deformability $\Lambda$ and the Love number $k_{2}$ obtained from the QMC/QMC$\omega\rho$ models, with and without pasta phases, and the three versions of the MQMC.

In Fig. 2(top) we display the Love number $k_{2}$ as a function of the compactness of the neutron star. In Fig. 2(bottom) we can see the variation of $\Lambda$ with the mass of the star. As already pointed out in jorge2018 , the curves do not collapse into one single curve because of their different dependence on $y$. When the same models are used to compute the tidal polarizability, their differences are enhanced and only two models, namely, QMC$\omega\rho$ and MQMC3 give results that are consistent with LIGO results for the canonical star. The pasta phase plays a very modest role and its influence is practically unnoticed, what can be seen if one compares the curves obtained from QMC with and without pasta and QMC$\omega\rho$ with and without pasta.

In Fig. 3 we show the tidal deformabilities of each neutron star in the binary system. $\Lambda_{1}$ is associated to the neutron star with mass $m_{1}$, which corresponds to the integration of every EoS in the range $1.37\leqslant m/M_{\odot}\leqslant 1.60$ obtained from GW170817. The mass $m_{2}$ of the companion star is determined by solving $\mathcal{M}_{c}=1.188M_{\odot}=\frac{(m_{1}m_{2})^{3/5}}{(m_{1}+m_{2})^{1/5}}$ PRL.119.161101 . We notice that only the QMC$\omega\rho$ model (with and without pasta phases), and the MQMC3 parametrization show values in between the confidence lines.

V.2 Cooling results

We now present the results of our thermal evolution studies, obtained by the numerical solution of eqs. (57) and (58). We begin by showing the cooling of a set of stars covering a wide range of masses for the QMC model. Those are shown in Fig. 4.

As shown in Fig. 4 all stars exhibit fast cooling, characterized by a sharp drop in their surface temperature at the age of $\sim 100$ years. Such fast cooling is a manifestation of a prominent presence of the direct Urca process (DU) in the stellar core, which is indeed the case for these stars. The DU process is extremely efficient in exhausting the star thermal energy, leading to the manifested fast cooling Prakash1992 . This, in turn indicates that such thermal evolution is mostly incompatible with observed data (unless the DU process is suppressed), as we will discuss in the following.

We now show in Fig. 5 the thermal evolution of QMC neutron stars with the pasta phase describing the inner crust.

We can see that the pasta phase has little effect on the overall cooling behavior of the star. This is to be expected as the pasta phase occupies only a small region. A more careful analysis lead us to the conclusion that the pasta phase had the minor effect of aiding in the core-crust thermalization, leading to a core-crust thermal relaxation on average $\sim 3.5$ years faster. This seems reasonable, as the pasta phase smooths the core-crust transition, thus, one can expect that will also smooth out heat propagation in between these regions. We note, however, that these are rough estimates and a more detailed calculation is warranted, one in which the thermal and conductivity properties of the pasta phase are explored in more details.

We now repeat the calculations above, but for the QMC$\omega\rho$ model, without pasta (Fig. 6) and with pasta (Fig. 7).

We can see that, as happens in the previous case (QMC), the pasta phase has little effect on the overall cooling - serving only to delay the thermalization by a few years. Differently than the previous model, however, in this case all stars exhibit slow cooling, meaning that there is no sharp drop in temperature when the core and crust become thermally coupled. The reason for that is that in this particular model the proton fraction at the core of the star, even at larger densities of heavier stars, is low enough to prevent the DU process from taking place, thus leading to a substantially slower thermal evolution.

Now we investigate the cooling of neutron stars described by the MQMC model in order to see if the different prescription used in the EoS calculations will play any role in the cooling. We show in Figs. 8 - 10 the cooling of neutron stars of different masses calculated for the different MQMC models studied.

We can see by the results shown in Figs. 8 - 10 that the overall cooling behavior associated with each model is, qualitatively, the same. For the MQMC models 1 and 2 we see that all stars exhibit fast cooling, with MQMC 3 displaying slow cooling (except for high mass stars).

It is worth mentioning that we see a clear correlation between the slope of the nuclear symmetry energy ($L_{0}$) and a fast/slow cooling behavior, with stars with higher values of $L_{0}$ having faster cooling (QMC, MQMC1 and MQMC2) whereas lower values of $L_{0}$ leading to slower cooling scenario (QMC$\omega\rho$ and MQMC3).

So far in our thermal investigation we have completely ignore effects of pairing, as we were interested in probing possible differences in the thermal behavior of the different models studied, and the inclusion of pairing could potentially murk the results. As it turns out all models behave similarly, with major differences connected with the different values of the symmetry energy slope rather than more specific minutia about the models. We now need to make sure that our model is in agreement with observed data, or at the very least is comparable to other results in the recent literature Negreiros2015a ; Negreiros2015b ; Beloin2018 ; Negreiros2018a ; Raduta2018 - for that we need to include pairing. Pairing in neutron stars may potentially occur in three different manners, neutron-neutron singlets (${}^{1}S_{0}$), neutron-neutron triplet (${}^{3}P_{2}$) and proton-proton singlet (${}^{1}S_{0}$). Neutron-neutron singlets are known to form mainly in the neutron-free region of the crust, whereas its triplet counterpart should form mainly in the lower densities of the core. As for the proton-proton pairs there is still great uncertainty as to how deep into the star it may take place. For a review of pairing in neutron stars we refer the reader to Yakovlev2000 . In this work we adopt a fairly standard pairing scenario, similar to the CCDK model of reference Page2004 . Below we show in (Figs. 11) the critical temperature as a function of neutron and proton Fermi momenta that were used to calculate the pairing effects in the neutron star cooling processes (suppression of the DU process, appearance of Pair-Breaking-Formation (PBF) process near $T_{c}$, and the modification of the specific heat of paired particles).

As discussed above, the presence of pasta does not change (qualitatively) the cooling behavior of the star, thus henceforth we will consider only the QMC with the pasta phase, as we believe that this is a more appropriate description for the core-crust transition. We now revisit the results shown in Figs. 5 - 10 taking into account the superfluidity model described above.

We begin by showing the thermal evolution of the QMC model (with pasta), which is shown in Fig. 12.

Fig. 13 shows the surface temperature evolution for the stars under the QMC$\omega\rho$ model with nucleonic pairing.

The results illustrated in Figs. 12 and 13 show us that, as expected, the presence of nucleon pairing leads to slower cooling (as neutrino emission processes involving paired nucleons are exponentially suppressed, once the superfluid phase is achieved). For the QMC$\omega\rho$ this is less relevant, as in this case the DU process was already suppressed due to the nature of this equation of state itself, which is connected to lower proton fractions. For this model, however, we see that a "second knee" appears near the 1000 year age. This second temperature drop is associated with the onset of neutron supefluidity on the outer layers of the star, which gives rise to a burst of neutrino emissions (associated with the Pair Breaking-formation process), leading to a sudden temperature drop. This phenomena has been used in previous studies to explain the apparent temperature drop of the neutron star in Cas A Page2011a ; Shternin2011 , although Cas A itself has been put into question recently Posselt2018 .

The results of the QMC (pasta) model (fig. 12) are more interesting for the purposes of this research. We see that with the presence of pairing, the cooling exhibits a "broader spectrum", with higher masses associated with faster cooling, smaller masses with lower, and intermediate masses in between. Such behavior makes this model ideal to confront against observational data, as it is more likely to fall into the range of observed neutron star temperatures.

We have also considered the thermal evolution of the MQMC stars with superfluidity. In here we show only the cooling of the stars of the MQMC1 and MQMC2, which are shown in Figs. 14 and  15. The cooling of superfluid stars of the MQMC3 model is omitted as they exhibit a similar behavior of that shown in Fig. 12.

These results show us that both MQMC1 and MQMC2 models lead to satisfactory thermal evolution behavior, leading to a large band covered by different stars masses that could potentially describe the observed data. As it is usual in thermal evolution studies, one cannot (usually) pin point the composition of the star by such investigations but rather rule out models that are less likely to describe observed data. We now compare the band of possible thermal evolutions (spanned by low-to-high mass stars) in the two models we have deemed best (in the cooling context) – these are shown in Figs 16 and 17.

We also confront the results of Figs 16 and 17 with most prominent thermal observed data (Beloin2018, ; SafiHarb2008, ; Zavlin1999, ; Pavlov2002, ; Mereghetti1996, ; Zavlin2007, ; Pavlov2001, ; Gotthelf2002, ; McGowan2004, ; Klochkov2015, ; McGowan2003, ; McGowan2006, ; Possenti1996, ; Halpern1997, ; Pons2002, ; Burwitz2003, ; Kaplan2003, ; Zavlin, ; Ho2015, ). – as shown in Fig. 18, from where we see that the major difference between the two models is that the MQMC2 spans a wider band, covering more of the low temperature region. We do not however have any data available on such region, possibly due to instrument sensitivity limitations. What we can see, however, is that both models have difficulties in matching the higher temperatures observed. This could be due to the presence of some unaccounted heating process or even to the uncertainties on the estimation of the star age. This can be remedied somewhat by taking into account a more sophisticated model for the neutron star atmosphere. For the present model we have only considered standard neutron star atmosphere Gudmundsson1983 . At the neutron star formation, some accretion might take place Potekhin1997 , depositing some materials onto the surface of the star. In order to account for that we consider the QMC model (with pasta and superfluidity) and allow for an accreted atmosphere of $\Delta M=10^{-10}M$. The results are shown in Fig. 19 and they show that the accreted atmosphere allow the surface to reach higher temperatures, making the agreement with observed data a little better. There is price paid, however, that is the speed up of cooling at later ages - this can be explained due to higher temperatures achieved by the star surface, that also enhance black-body photon emission (which is the dominant cooling mechanism in the late stages of evolution). We also see that a modest accreted atmosphere is not enough to explain the higher temperature stars – which is in agreement with other studies (see Negreiros2010 for instance). One must also keep in mind that there is great uncertainty regarding the age of the star Page2004 .

VI Final Remarks

In the present work we have used five different versions of the QMC model to compute the Love number and tidal polarizability and confront them with GW170817 constraints and also to investigate the cooling process they describe. Two of the models are based on the original bag potential structure and three versions consider a harmonic oscillator potential to confine the quarks. The bag-like models also incorporate the pasta phase used to describe the inner crust of neutron stars.

We have compared our EOS with the FSUGarnet jorge2018 and checked that the outer crusts are coincident and the liquid core of all EOS are very similar, most of the differences residing on the inner crust and around the crust-core transition region, where the pasta phase was included. Our results point to the fact that the pasta phase plays only a minor role in the calculation of the dimensionless tidal polarizability and just two of the models (QMC$\omega\rho$ and MQMC3) give results that lie within the expected range of the canonical star $\Lambda_{1.4}$ and, at the same time, appear inside the region delimited by the confidence lines. At this point we see no clear correlation between the model symmetry energy or its slope and the fact that it satisfies (or not) the constraints investigated.

We have then used the very same models to study the cooling process and have verified that, in this case, there is a clear correlation between the slope of the symmetry energy and the velocity of the cooling process, i.e., models with higher (lower) values of the slope produce fast (slow) cooling. The pasta phase, once again, plays a minor role in the cooling process by speeding up the process in about 3.5 years if superconductivity is not considered. The models that seems to better describe the observational data are the QMC (with pasta) and MQMC2, but none can reach the region of very high temperatures.

In fact, our studies show that we cannot pin down one unique model that can, at the same time, describe astrophysical quantities and a perfect description of the possible cooling processes and more observational data is necessary.

Acknowledgements.

This work is a part of the project INCT-FNA Proc. No. 464898/2014-5, was partially supported by CNPq (Brazil) under grants 301155.2017-8 (D.P.M.), 310242/2017-7 and 406958/2018-1 (O.L.), 308486/2015-3 (T.F.) and 433369/2018-3 (M.D.), by Capes-PNPD program (C.V.F.), and by Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP) under the thematic projects 2013/26258-4 (O.L., T.F.) and 2017/05660-0 (T.F.). R.N. also acknowledges that this project was partly funded by FAPERJ, under grant E-26/203.299/2017.