Density-dependent van der Waals model under the GW170817 constraint

TL;DR

This paper introduces a density-dependent van der Waals model that accurately describes high-density nuclear matter and neutron star properties, aligning with recent astrophysical and gravitational wave observations.

Contribution

It develops a generalized van der Waals model with density-dependent attraction, successfully applying it to neutron star structure and merger event constraints.

Findings

Neutron star mass range: 1.97-2.07 solar masses.

Model predictions agree with GW170817 tidal deformability bounds.

Mass-radius relations match recent X-ray data.

Abstract

We propose a density-dependent function for the attractive interaction in the original van der Waals model to correctly describe the flow constraint at the high-density regime of the symmetric nuclear matter. After a generalization to asymmetric nuclear matter, it was also possible to study the stellar matter regime from this new model. The mass-radius relation for neutron stars under $\beta$-equilibrium is found to agree with recent X-ray observations. The neutron star masses supported against gravity, obtained from some parametrizations of the model, are in the range of $(1.97-2.07)M_{\odot}$, compatible with observational data from the PSR J0348+0432 pulsar. Furthermore, we verify the reliability of the model in predicting tidal deformabilities of the binary system related to the GW170817 neutron star merger event and find a full agreement with the new bounds obtained by the LIGO/Virgo collaboration.

Equations

$$\mathcal{E}(\rho )=(1-b\rho ){\mathcal{E}}_{\mathrm{id}}^{\ast }({\rho }^{\ast })-a{\rho }^2,$$

$${\mathcal{E}}_{\mathrm{id}}^{\ast }({\rho }^{\ast })=\frac{\gamma }{2{\pi }^2}{\int }_0^{k_F^{\ast }}dk{k}^2(k^2+M^2{)}^{1/2},$$

$$Z_{\mathrm{CS}}(\eta )=\frac{1+\eta +{\eta }^2-{\eta }^3}{(1-\eta {)}^3}=1+\sum _{j=0}^{\infty }(j^2+3j){\eta }^j,$$

$$\mathcal{E}(\rho )=f(\eta ){\mathcal{E}}_{\mathrm{id}}^{\ast }-a{\rho }^2,$$

$$\mathcal{E}(\rho )=[1-b(\rho )\rho ]{\mathcal{E}}_{\mathrm{id}}^{\ast }({\rho }^{\ast })-a(\rho ){\rho }^2,$$

$$b(\rho )=\frac{1-f(\eta )}{\rho }=\frac{1}{\rho }-\frac{1}{\rho }\exp \left[\frac{-(4-\frac{3b\rho }{4})\frac{3b\rho }{4}}{{\left(1-\frac{3b\rho }{4}\right)}^2}\right].$$

$$a(\rho )$$

$$P(\rho )$$

$$P_{\mathrm{id}}^{\ast }=\frac{\gamma }{6{\pi }^2}{\int }_0^{k_F^{\ast }}\frac{dk{k}^4}{(k^2+M^2{)}^{1/2}}.$$

$$K$$

$${\Sigma }^{\prime }(\rho )$$

$$\mu (\rho )$$

$$a(\rho )$$

$$\mathcal{E}(\rho ,y)$$

$${\rho }_p^{\ast }=\frac{{\rho }_p}{1-b(\rho )\rho },{\rho }_n^{\ast }=\frac{{\rho }_n}{1-b(\rho )\rho }.$$

$$P(\rho ,y)=P_{\mathrm{id}}^{\ast }-a(\rho ){\rho }^2+\rho \Sigma (\rho ,y)+d(2y-1{)}^2{\rho }^2,$$

$${\mu }_{p,n}(\rho ,y)=\frac{\partial \mathcal{E}}{\partial {\rho }_{p,n}}$$

$$-2a(\rho )\rho \pm 2d(2y-1)\rho ,$$

$${\Sigma }^{\prime }(\rho ,y)$$

$$+(b^{\prime \prime }\rho +b^{\prime })P_{\mathrm{id}}^{\ast }({\rho }_p^{\ast },{\rho }_n^{\ast })-a^{\prime \prime }{\rho }^2-2a^{\prime }\rho ,$$

$$\displaystyle\mathcal{S}(\rho)\simeq\mathcal{S}_{2}(\rho)=\frac{1}{8}\frac{\partial^{2}(\mathcal{E}/\rho)}{\partial y^{2}}\Big{|}_{y=\frac{1}{2}}=\mathcal{S}_{kin}^{*}(\rho)+d\rho,$$

$$L(\rho )=3\rho \frac{\partial \mathcal{S}}{\partial \rho }=\xi (\rho )L_{kin}^{\ast }(\rho )+3d\rho ,$$

$$L_{kin}^{\ast }(\rho )=\frac{k_F^{\ast 2}}{3E_F^{\ast }}\left(1-\frac{k_F^{\ast 2}}{2E_F^{\ast 2}}\right)=2{\mathcal{S}}_{kin}^{\ast }\left(1-\frac{3{\mathcal{S}}_{kin}^{\ast }}{E_F^{\ast }}\right)$$

$$L(\rho )=3\mathcal{S}(\rho )+{\mathcal{S}}_{kin}^{\ast }(\rho )\left\{2\xi \left[1-\frac{3{\mathcal{S}}_{kin}^{\ast }(\rho )}{E_F^{\ast }(\rho )}\right]-3\right\}.$$

$${\mathcal{E}}_T(\rho ,{\rho }_e,y)$$

$$P_T(\rho ,{\rho }_e,y)$$

$$k_2$$

$$Q(r)$$

Full text

van der Waals model, neutron star matter, GW170817 calculations

1 Introduction

The applicability of hadronic equation(s) of state (EoS) goes from the description of superheavy nuclei to the structure of neutron stars. Therefore, a complete understanding of nuclear physics, as well as the astrophysics involved in these environments, is needed. The hadronic EoS is one of the most important ingredients to correctly predict, for instance, the mass of a neutron star. This object is one of the densest in the visible universe, having its density around $5-6$ times the nuclear saturation density (Lattimer & Prakash, 2004). It also provides a unique natural laboratory to test the hadronic EoS profile at the high-density regime (Lattimer & Prakash, 2004). Previously, some of the few constraints on the EoS, coming from the astrophysical context, were the values of the neutron star mass and the canonical star radius. Nowadays, the recent observation of the gravitational waves detected from a binary neutron star merger event, named as GW170817 (Abbot et al., 2018) provides an opportunity to constrain various properties of hadronic EoS, and to search for more realistic ones, which will give a clear picture of the physics of a neutron star. The internal structure of such an object is a controversial subject since there are indications of many interesting phenomena like kaon production (Glendenning & Schaffner-Bielich, 1998, 1999; Gupta & Arumugam, 2012; Glendenning, 1985), hyperons emergence (Glendenning & Schaffner-Bielich, 1998; Glendenning & Moszkowski, 1991; Ambartsumyan et al., 1960), and the transition to a deconfinement quark phase (Collins & Perry, 1975). Some of them are directly affected by the value of the bulk parameters, at the saturation density, given by the hadronic models, such as the symmetry energy ($J$), the incompressibility ($K_{0}$), the skewness parameter ($Q_{0}$), and related quantities. Many works establish constraints on these quantities, see, for instance, Ref. (Dutra et al., 2014). As an example, the incompressibility, which defines the stiffness of the EoS, has a value in the range of $K_{0}=(240\pm 10$) MeV, as found in Refs. (Colo et al., 2004; Todd-Rutel & Piekarewicz, 2005; Agrawal et al., 2005), predicted from the isoscalar giant monopole resonance of the heavy nuclei, or even the range of $250~{}\mbox{MeV}\leqslant K_{0}\leqslant 315~{}\mbox{MeV}$, more recently obtained in Ref. (Stone et al., 2014) from a reanalysis of updated data on isoscalar giant monopole resonance energies of Sn and Cd isotopes. Some ranges for $K_{0}$ are found through a leptodermous expansion of the finite nucleus incompressibility, with $K_{0}$ as one of the terms. However, many works point out to the drawbacks of such a procedure, see, for instance, Refs. (Shlomo & Youngblood, 1993; Chen et al., 2009; Pearson et al., 2010). The current consensus regarding the value of $K_{0}$ is $220~{}\mbox{MeV}\leqslant K_{0}\leqslant 260~{}\mbox{MeV}$, as one can see in a very recent review on this subject in Ref. (Garg & Colo, 2018), for instance. There are also a lot of uncertainties in the density dependence of some of the bulk parameters, such as the symmetry energy, especially at the high-density regime (Steiner et al., 2005; Li et al., 2008; Lynch et al., 2009; Li et al., 2014; Horowitz et al., 2014; Baldo & Burgio, 2016; Li, 2017; Li et al., 2019).

In this manuscript, we discuss the applicability of our proposed density-dependent van der Waals (DD-vdW) model, namely, an improved version of the previous van der Waals (vdW) model applied to nuclear matter (Vovchenko et al., 2015a, b; Vovchenko, 2017), that takes into account density-dependent repulsive and attractive interactions. Though the vdW model takes the description of the nuclear interactions in a simple way, its predictive capacity is as promising as other realistic relativistic (Walecka, 1974; Boguta & Bodmer, 1977; Serot & Walecka, 1979) and nonrelativistic models (Skyrme, 1959; Vautherin & Brink, 1972), since it can reproduce some basic properties such as the saturation at a specific density, and the liquid-gas phase transition in symmetric nuclear matter. However, it presents some limitations like the causality violation at the low-density regime, which makes it inappropriate to describe stellar matter. In the DD-vdW model, this problem is circumvented. The neutron star calculations are performed, including the comparison of the tidal deformabilities of a binary neutron star system with the corresponding data related to the GW170817 event, recently reported in Ref. (Abbot et al., 2018).

This manuscript is organized as follows: in Sec. 2 we give an outline of the theoretical formalism that establishes the basis for a density-dependent version of the vdW model. In Sec. 3, we discuss its causality limitation and show how its new version, proposed here, describes nuclear matter at higher densities, including the compatibility with the flow constraint predicted in Ref. (Danielewicz et al., 2002). A generalization to asymmetric systems is also developed, and special attention is given to stellar matter and tidal deformability calculations related to the GW170817 event. We finish the manuscript with a summary and concluding remarks in Sec. 4.

2 Density dependent vdW model

By following the calculations performed at zero temperature regime in Refs. (Vovchenko et al., 2015a, b), the energy density of an infinite symmetric nuclear matter (SNM) system in the grand canonical ensemble, can be obtained from the vdW model, as

[TABLE]

where $\mathcal{E}^{*}_{\rm id}$ is the kinetic energy of a relativistic ideal Fermi gas of nucleons of mass $M=938$ MeV, given by

[TABLE]

with $k_{F}^{*}=(6\pi^{2}\rho^{*}/\gamma)^{\frac{1}{3}}$, and $\rho^{*}=\rho/(1-b\rho)$. The degeneracy factor, $\gamma=4$ for SNM and $b=16\pi r^{3}/3$, being $r$ the nucleon hard-sphere radius. The parameters $a$ and $b$ are related to the strength of attractive and repulsive interactions, respectively, and can be obtained by forcing the model to present a bound state at a particular density. In the case of infinite nuclear matter, the binding energy of such state is given by $\mathcal{E}(\rho_{0})/\rho_{0}-M=-B_{0}$ at the saturation density $\rho_{0}$. Usually, in the mean-field models, the binding energy per particle ($B_{0}\approx 16.0$ MeV) at the saturation density $\rho_{0}\approx 0.16$ fm*-3* are well established observables, and constrain the value of the parameters to $a=328.93$ MeV fm3 and $b=3.41$ fm3.

From the structure of the excluded volume mechanism, it is clear that the vdW model has a density range of $\rho<\rho_{\mbox{\tiny max}}$, with $\rho_{\mbox{\tiny max}}=b^{-1}$. For $b=3.41$ fm3, one has $\rho_{\mbox{\tiny max}}=1.83\rho_{0}$. In order to avoid such limitations, the repulsive term was modified in Ref. (Vovchenko, 2017), by applying the Carnahan-Starling (CS) (Carnahan & Starling, 1969) method of excluded volume, in which the pressure of hard-core particles of radius $r$ is given by $P=\rho TZ_{\rm CS}(\eta)$, with

[TABLE]

and $\eta=b\rho/4$. By following this method, the first eight virial expansion coefficients are obtained, unlike the traditional vdW excluded volume procedure, where only two of them are recovered since for this case $Z(\eta)=(1-4\eta)^{-1}$. By taking the CS procedure to the vdW model, the nuclear matter energy density in the grand canonical ensemble is written as (Vovchenko, 2017)

[TABLE]

with $f(\eta)=e^{-(4-3\eta)\eta/(1-\eta)^{2}}$ and $\rho^{*}=\rho/f(\eta)$. We name this model as vdW-CS one. For this model, it is found $a=347.02$ MeV fm3 and $b=4.43$ fm3 with maximum density given by $\rho_{\mbox{\tiny max}}=4b^{-1}=5.64\rho_{0}$ (Vovchenko, 2017).

In an effective way, the vdW-CS model can be seen as a density-dependent model by rewriting Eq. (4) as

[TABLE]

with $\rho^{*}=\rho/[1-b(\rho)\rho]$ and

[TABLE]

Thus, the repulsive interaction becomes a density-dependent function. We generalize this idea to the attractive part, regulated by the $a$ parameter, by making it depending on $\rho$ as well, i.e., we assume $a\rightarrow a(\rho)$ in Eq. (5). Therefore, the original vdW model is given by the particular case in which $b(\rho)=b$ and $a(\rho)=a$. For the vdW-CS model, $b(\rho)$ is given by Eq. (6) and $a(\rho)=a$.

From the perspective of a density-dependent model, it is also possible to use Eq. (5) to describe all real gases models analyzed in Ref. (Vovchenko, 2017), namely, Redlich-Kwong-Soave (RKS), Peng-Robinson (PR), and Clausius-2 (C2), by identifying

[TABLE]

and

[TABLE]

We denote here the last model by Clausius-2 (C2) because it is a two parameters ($a$, $b$) version of the original Clausius model. A three parameters ($a$, $b$, $c$) version of such model is studied in Ref. (Vovchenko et al., 2018).

From Eq. (5), the pressure of the system can be written as

[TABLE]

with

[TABLE]

Here, $a^{\prime}$ and $b^{\prime}$ are the first density derivatives of $a(\rho)$ and $b(\rho)$, respectively. Furthermore, the density dependence of the interactions gives rise to a rearrangement term in the pressure given by $\Sigma(\rho)=b^{\prime}\rho P^{*}_{\rm id}-a^{\prime}\rho^{2}$. The incompressibility of SNM is calculated as the derivative of pressure with respect to $\rho$ as

[TABLE]

where $K^{*}_{\rm id}=3k_{F}^{*2}/\sqrt{k_{F}^{*2}+M^{2}}$, and

[TABLE]

Finally, the chemical potential is obtained from

[TABLE]

with $\mu^{*}_{\rm id}=E_{F}^{*}=(k_{F}^{*2}+M^{2})^{1/2}$. We remind the readers that the rearrangement term $\Sigma(\rho)$, originated from the density dependence of the attractive and repulsive interactions, is an essential quantity in order to preserve the thermodynamic consistency of the model. It is straightforward to see that equations of state presented here with $\Sigma(\rho)$ included, give rise to the relation $P+\mathcal{E}=\mu\rho$, indicating the consistency.

3 Results and Discussions

Since the structure and EoS of the density-dependent vdW model were presented, we are now able to analyze in more details the real gases of Ref. (Vovchenko, 2017), namely, those given by Eqs. (7)-(9). At present, a diverse set of nonrelativistic and relativistic mean-field (RMF) models describe quite well nuclear and neutron star matter (Dutra et al., 2014, 2012; Lourenço et al, 2007; Lourenço et al., 2016). Hence, it is also important to see how the density-dependent vdW models describe the same environment.

3.1 Causality analysis from flow constraint

Our first analysis is concerning the flow constraint proposed in Ref. (Danielewicz et al., 2002), in which limits on the pressure-density relationship of SNM and its curvature can be obtained from the experimental data on the motion of ejected matter in the energetic nucleus-nucleus collisions. Measurements of the particle flow in the collisions of ${}^{197}\rm Au$ nucleus at incident kinetic energy per nucleon varying from about $0.15$ to $10$ GeV is used in Ref. (Danielewicz et al., 2002), that extrapolated for pressure at the range of $2.0\leqslant\rho/\rho_{0}\leqslant 4.6$, at the zero temperature. In principle, the CS excluded volume procedure enables the models to reach values for $\rho_{\mbox{\tiny max}}$ greater than those found by using the traditional excluded volume method. However, the causality is violated at densities lower than $\rho_{\mbox{\tiny max}}$.

We remind the reader that the conventional excluded volume technique for nucleons are taken into account for repulsion at short distances by treating them as rigid spheres in a nonrelativistic context. In a relativistic framework, as in the case of nuclear matter at high densities, the Lorentz contraction of such hard sphere nucleons should be implemented to avoid causality violation for any density regions (Bugaev, 2008). Thus, the effect of Lorentz contraction in nucleons can be seen as a decreasing function of their radius. Effectively, this is the case of the CS method, since it treats the excluded volume parameter as a density-dependent quantity, see Eq. (6). The method identifies the nucleon as a sphere in which the radius is depending on the density, at least as a first approach. However, such a density-dependent is not enough to completely avoid the causality violation. The sound velocity is still greater than $1$ (units of $c=1$) in CS approach, but at higher densities in comparison with the case in which conventional excluded volume procedure (fixed radius) is taken.

In Fig. 1, we show the squared sound velocity, $v_{s}^{2}=\partial P/\partial\mathcal{E}=K/9\mu$, as a function of the ratio $\rho/\rho_{0}$ for the real gases models presented in the last section.

It is clear that the repulsive interaction plays an important role in the causal limit, since it induces its violation at higher densities. The physical reason of such a result is that the CS method weakens the repulsive interaction as a function of density, producing results closer to the ones of an ideal gas of massive point nucleons. The more $b(\rho)$ decreases, the more the nucleons behave like point particles, since $b(\rho)\rightarrow 0$ indicates structureless objects.

Even by applying the CS excluded volume method in the real gases models, the causality is still broken in the range of $2.0\leqslant\rho/\rho_{0}\leqslant 4.6$ of the flow constraint. In Fig. 1 we see that $v^{2}_{s}=1$ at $\rho/\rho_{0}=2.5$ for the C2-CS model, for instance. In order to circumvent this limitation, we implement a modification in the attractive interaction of the model. A new proposed form for the parameter $a(\rho)$, is given by,

[TABLE]

It is inspired in the C2-CS model, in which causality is violated at higher densities in comparison with the remaining models. Its repulsive interaction remains the same, i.e., the CS excluded volume method is considered. The model with the new proposal for $a(\rho)$ is named here as the DD-vdW model, with the couplings $a(\rho)$ and $b(\rho)$ given by Eqs. (15) and (6), respectively. Notice that for the particular cases of $n=0$, and $n=1$, the vdW-CS and C2-CS models are reproduced, respectively. The squared sound velocity for some values of $n$ is displayed in Fig. 2.

Notice that the effect of the $n$ power in $a(\rho)$, Eq. (15), is to weaken the strength of the attractive interaction. The more $a(\rho)$ decreases, the more the DD-vdW model approaches to the free Fermi gas of massive particles. The combined effects of the density-dependent parameters $a(\rho)$ and $b(\rho)$ favor the model to move the causality violation to higher densities, as one can see in Fig. 2.

3.2 DD-vdW model in symmetric nuclear matter

The $n$ power in Eq. (15) directly affects the incompressibility value at the saturation density, $K_{0}\equiv K(\rho_{0})$. Therefore, we can use this parameter to control the $K(\rho_{0})$ quantity. We use this procedure to submit the DD-vdW model to the flow constraint of Ref. (Danielewicz et al., 2002). The resulting parametrizations of this model with different $K_{0}$ values are shown in Fig. 3.

From this figure, we see that the DD-vdW model satisfies the flow constraint for parametrizations presenting $242.4~{}\mbox{MeV}\leqslant K_{0}\leqslant 251.5~{}\mbox{MeV}$. All the curves of Fig. 3 are compatible with the causal limit, i.e., for each curve exhibited one has $v_{s}^{2}<1$. The parametrization with $K_{0}=251.5$ MeV, for instance, is causal for $\rho/\rho_{0}\leqslant 4.9$. Furthermore, all these parametrizations are also in agreement with the restriction of $220~{}\mbox{MeV}\leqslant K_{0}\leqslant 260~{}\mbox{MeV}$ (Garg & Colo, 2018). It is worth to mention that in the recent work of Ref. (Vovchenko, 2017), the results pointed out that none of the other real gases models produces $K_{0}$ inside the aforementioned range for $K_{0}$. Finally, it is also clear from Fig. 3 that the weakening effect of the interactions leads the DD-vdW model to the direction of the free Fermi gas, with full agreement with the flow constraint. Another approach performed in the van der Waals model that also makes it to satisfy the flow constraint, and also the maximum mass observational data for neutron stars, includes induced surface tension in its formulation. For details, see (Sagun et al., 2018; Bugaev et al., 2019) and references therein.

3.3 Asymmetric matter formulation

In order to proceed for a complete analysis of the infinite nucleonic matter, it is necessary to take the isospin asymmetry system into account, i.e., the EoS for $y\equiv\frac{\rho_{p}}{\rho}\neq\frac{1}{2}$. Here $\rho_{p}$ is the proton density. For the original vdW model, a generalization for different hard sphere particles was performed in Ref. (Vovchenko et al., 2017), where EoS were developed for neutron-proton and nucleons-$\alpha$ particles systems. Here, in order to avoid the emergence of more free parameters, we consider the same density-dependent couplings $a(\rho)$ and $b(\rho)$ of the symmetric matter, Eqs. (15) and (6), for protons and neutrons. Therefore, the individual components are distinguished only by their respective kinetic energies. However, as already discussed in Ref. (Vovchenko et al., 2017) in the case of $a(\rho)=a$ and $b(\rho)=b$, the symmetry energy at the saturation density calculated from this approach, $J\equiv\mathcal{S}(\rho_{0})$, presents an underestimate value. In order to avoid such a limitation, we propose a new term in the energy density of the DD-vdW model, proportional to the squared difference between protons and neutrons densities, $\rho_{3}=\rho_{p}-\rho_{n}=(2y-1)\rho$, as widely used in some RMF models (Dutra et al., 2014). This new term can be seen as a simulation of the $\rho$ meson exchange in its simple form, i.e., a minimal coupling between this meson and the finite particle nucleon (excluded volume included). Thus, the asymmetric nuclear matter energy density is given by,

[TABLE]

where $\mathcal{E}^{*}_{\rm id}(\rho^{*}_{p},\rho^{*}_{n})=\mathcal{E}^{*p}_{\rm id}(\rho^{*}_{p})+\mathcal{E}^{*n}_{\rm id}(\rho^{*}_{n})$, for $\mathcal{E}^{*i}_{\rm id}(\rho^{*}_{i})$ following the same form as in Eq. (2) with $\gamma=2$, $k_{F}^{*}\rightarrow k_{F}^{*i}$ and $\rho^{*}\rightarrow\rho^{*}_{i}$ ($i=p,n$). The different densities are related to each other by

[TABLE]

The couplings are written as in the previous case, see Eqs. (6) and (15), respectively, for $b(\rho)$ and $a(\rho)$.

From Eq. (16), we can derive the remaining quantities for the asymmetric nuclear matter. The expressions are,

[TABLE]

and

[TABLE]

for the pressure, and chemical potentials for protons (upper sign) and neutrons (lower sign), respectively. Furthermore, $P^{*}_{\rm id}(\rho^{*}_{p},\rho^{*}_{n})=P^{*p}_{\rm id}(\rho^{*}_{p})+P^{*n}_{\rm id}(\rho^{*}_{n})$ is written as in Eq. (11) with $\gamma=2$, $k_{F}^{*}\rightarrow k_{F}^{*i}$ and $\rho^{*}\rightarrow\rho^{*}_{i}$. The “ideal” chemical potentials are $\mu^{*i}_{\rm id}=E^{*i}_{F}$. The rearrangement term and its derivative with respect to the density are given by $\Sigma(\rho,y)=b^{\prime}\rho P^{*}_{\rm id}(\rho^{*}_{p},\rho^{*}_{n})-a^{\prime}\rho^{2}$, and

[TABLE]

with $K^{*i}_{\rm id}(\rho^{*}_{i})$ defined as in the symmetric nuclear matter case.

By defining the DD-vdW model generalized to asymmetric matter in this way, there are only four free parameters to be adjusted, namely, $a$, $b$, $n$, and $d$. The first three ones are already determined from the symmetric case (reproducing the values of $\rho_{0}$, $B_{0}$, and $K_{0}$). The remaining free parameter, $d$, is adjusted in order to correctly reproduce a bulk quantity of asymmetric nuclear matter, namely, the symmetry energy. This quantity measures the change in binding of the nucleon system as the proton to neutron ratio is changed at a fixed value of the density, $\mathcal{S}(\rho)=E(\rho,0)-E(\rho,1/2)$, where $E(\rho,y)$ is the energy per particle. A detailed analysis of the quantity is quite important for understanding many aspects of different isospin asymmetric systems, from astrophysics to finite nuclei. Furthermore, the symmetry energy slope at saturation density provides the dominant contribution to the pressure in neutron stars, as well as affects the neutron skin thicknesses of heavy nuclei (Horowitz & Piekarewicz, 2001; Bhuyan, 2015; Bhuyan et al., 2018). For a recent review regarding the importance of $\mathcal{S}(\rho)$, see Ref. (Baldo & Burgio, 2016). By considering $E(\rho,y)\simeq E(\rho,1/2)+\mathcal{S}_{2}(\rho)(1-2y)^{2}+\mathcal{O}[(1-2y)^{4}]$, one can take $\mathcal{S}(\rho)\simeq\mathcal{S}_{2}(\rho)$ as a good approximation in order to compute the symmetry energy. The $\mathcal{S}(\rho)$ of the DD-vdW model can be written as,

[TABLE]

with $\mathcal{S}_{kin}^{*}(\rho)=k_{F}^{*2}/(6E_{F}^{*})$ and $k_{F}^{*}=(3\pi^{2}\rho^{*}/2)^{\frac{1}{3}}$.

One can notice that the excluded volume procedure directly affects the kinetic part of $\mathcal{S}(\rho)$ as well as all in other thermodynamical quantities. The determination of the $d$ parameter is straightforward since we use the analytical expression of $\mathcal{S}(\rho)$ for this aim by imposing our model to present consistent values for the symmetry energy at saturation density, $J$. We find $d$ by constraining $J$ to the range of $25~{}\mbox{MeV}\leqslant J\leqslant 35~{}\mbox{MeV}$. This range was established in order to encompass data obtained from various terrestrial nuclear experiments and astrophysical observations, such as isospin diffusion, neutron skins, pygmy dipole resonances, modes of decay near the drip-line, transverse flow, mass-radius relations, and torsional crust oscillations of neutron stars. One can find a collection of these data in Ref. (Li & Han, 2013).

The slope parameter, i.e., the symmetry energy slope as a function of density, is obtained from Eq. (21) as

[TABLE]

where

[TABLE]

and $\xi(\rho)=[{1+b^{\prime}\rho^{2}}]/[{1-b(\rho)\rho}]$. The advantage of the specific form of the last term included in the Eq. (16) is the possibility of an analytical relationship between $\mathcal{S}(\rho)$ and $L(\rho)$ for all densities. Notice that because $d=(\mathcal{S}-\mathcal{S}_{kin}^{*})/\rho$, it is possible to write

[TABLE]

At the saturation density, this equation is used to calculate $L(\rho_{0})\equiv L_{0}$. The quantity $L_{0}$ is of great interest for constraining the EoS of asymmetric nuclear matter in many hadronic models (Baldo & Burgio, 2016; Santos et al., 2015, 2014; Holt & Lim, 2018). For the DD-vdW model, we find $L_{0}$ in the range of $63.4~{}\mbox{MeV}\leqslant L_{0}\leqslant 96.5~{}\mbox{MeV}$, by taking into account the constraint of $25~{}\mbox{MeV}\leqslant J\leqslant 35~{}\mbox{MeV}$ and the range of $242.4~{}\mbox{MeV}\leqslant K_{0}\leqslant 251.5~{}\mbox{MeV}$ obtained from the flow constraint analysis. The obtained values for $L_{0}$ are in agreement with the constraint of $25~{}\mbox{MeV}\leqslant L_{0}\leqslant 115~{}\mbox{MeV}$ used in Refs. (Dutra et al., 2014; Li & Han, 2013).

3.4 Stellar matter

A neutron star (NS) is a very compact object composed not only by neutrons, but also by protons and leptons. Different reactions such as the $\beta$ decay, namely, $n\rightarrow p+e^{-}+\bar{\nu}_{e}$ and its inverse process $p+e^{-}\rightarrow n+\nu_{e}$, take place in the interior of a NS. For densities in which the electron chemical potential exceeds the muon mass value, the reactions $e^{-}\rightarrow\mu^{-}+\nu_{e}+\bar{\nu}_{\mu}$, $p+\mu^{-}\rightarrow n+\nu_{\mu}$ and $n\rightarrow p+\mu^{-}+\bar{\nu}_{\mu}$ may be energetically allowed. In this case, muons can also emerge. Here, we consider that all neutrinos escape from the star. By taking these assumptions into account, one can write total energy density and pressure of the stellar system, respectively, as

[TABLE]

and

[TABLE]

where, $\mathcal{E}(\rho,y)$ and $P(\rho,y)$ are given in the Eqs. (16) and (18), respectively. The chemical equilibrium and the charge neutrality conditions are given by $\mu_{n}(\rho,y)-\mu_{p}(\rho,y)=\mu_{e}(\rho_{e})$ and $\rho_{p}(\rho,y)-\rho_{e}=\rho_{\mu}(\rho_{e})$, with $\mu_{p}$ and $\mu_{n}$ defined in Eq. (19). Furthermore, one has $\mu_{e}=(3\pi^{2}\rho_{e})^{1/3}$, $\rho_{p}=y\rho$, $\rho_{\mu}=[(\mu_{\mu}^{2}-m_{\mu}^{2})^{3/2}]/(3\pi^{2})$, and $\mu_{\mu}=\mu_{e}$, for $m_{\mu}=105.7$ MeV and massless electrons. Thus, for each input density $\rho$, the quantities $\rho_{e}$ and $y$ are calculated by solving the restrictions for the chemical potentials and densities simultaneously. The output is used to compute $\mathcal{E}_{T}(\rho,\rho_{e},y)$ and $P_{T}(\rho,\rho_{e},y)$ as a function of the density.

In Fig. 4, we show the EoS of neutron star matter under $\beta$-equilibrium for DD-vdW parametrizations consistent with the flow constraint. All the $P_{T}\times\mathcal{E}_{T}$ curves in the figure are consistent with the causal limit, i.e., the parametrizations are restricted to density ranges in which $v_{s}^{2}<1$. In these curves, one can observe the effects of the bulk parameters. For example, it is clear that $K_{0}$ plays a major role in the EoS of the neutron star matter, exactly as in the flow constraint, see Fig. 3. A more careful inspection shows that the results marginally depend on the symmetry energy. It is also verified that the curves are also compatible with the results found by Steiner (Steiner et al., 2010) and Nättiliä (Nättiliä et al., 2016).

The mass-radius diagrams of spherically symmetric neutron stars are found from the solutions of the Tolman-Oppenheimer-Volkoff (TOV) equations (Tolman, 1939; Oppenheimer & Volkoff, 1939). To solve such equations, we take the $\beta$-equilibrated energy density and pressure under chemical equilibrium and charge neutrality given by the DD-vdW model, along with the Baym-Pethick-Sutherland (BPS) equation of state (Glendenning, 2000; Baym et al., 1971) in the low density regime, namely, at $0.1581\times 10^{-10}$ fm${}^{-3}<\rho<0.008907$ fm*-3*, in order to specifically describe the NS crust. Certainly, a more profound study regarding a detailed description of the crust (nonuniform and clustered matter) of the neutron star (Oertel et al., 2017; Ducoin et al., 2011; Atta & Basu, 2014; Rueda et al., 2014; Carreau et al., 2019), even considering the pasta structure in the inner crust, is important and addressed to future works. Here we used a simple equation of state, the BPS one, connected to the DD-vdW model as a first approach to describe the stellar matter. In Fig. 5, we show the mass-radius diagrams obtained from the DD-vdW model.

From the figure, one can see that the DD-vdW parametrizations are consistent with the findings related to the PSR J0348+0432 pulsar, namely, $M_{\mbox{\tiny NS}}=(2.01\pm 0.04)M_{\odot}$ (Antoniadis et al., 2013), and also with the calculations performed by Steiner (Steiner et al., 2010) and Nättiliä (Nättiliä et al., 2016). Another observational range for the NS mass is given in Ref. (Demorest et al., 2010) and is related to the PSR J1614-2230 pulsar, in which $M_{\mbox{\tiny NS}}=(1.97\pm 0.04)M_{\odot}$. These values were recently modified (shifted down) to $M_{\mbox{\tiny NS}}=(1.928\pm 0.017)M_{\odot}$ in Ref. (Fonseca et al., 2016). The range of $K_{0}$ for the parametrizations compatible with the observational constraint of $M_{\mbox{\tiny NS}}\sim 2M_{\odot}$ is $247.8~{}\mbox{MeV}\leqslant K_{0}\leqslant 251.5~{}\mbox{MeV}$, with the symmetry energy at the saturation density limited to $25~{}\mbox{MeV}\leqslant J\leqslant 35~{}\mbox{MeV}$. Once again, we restricted the curves to a density range compatible with the causal limit. For example, for the parametrization in which $J=25$ MeV and $K_{0}=247.8$ MeV, the causal limit allows the density range up to $\rho_{/}\rho_{0}=5.66$. On the other hand, the maximum density reaches a value of $\rho_{/}\rho_{0}=5.00$ for the parametrization in which $J=35$ MeV and $K_{0}=251.5$ MeV. In the figure, we also display the predictions of the LIGO/Virgo Collaboration for the radii related to the masses of the binary neutron star system of the GW170817 event. We see that the DD-vdW parametrizations are also consistent with this constraint.

3.5 Tidal deformability calculations

In the stellar matter context, the tidal deformability (TD) is the measure of the deformation in an NS due to an external field. In case of a binary NS system, the TD in one star is due to the gravitational field created by its companion. In a brief mathematical language, one can say that $Q_{ij}=-\lambda{\varepsilon}_{ij}$ is the relationship between the TD ($\lambda$), the quadrupole moment $Q_{ij}$ developed in the NS, and the external tidal field ${\varepsilon}_{ij}$ (Hinderer, 2008; Flanagan & Hinderer, 2008; Hinderer et al., 2010; Damour et al., 2012). The dimensionless TD ($\Lambda$) is written in terms of the Love number $k_{2}$ as $\Lambda=2k_{2}/(3C^{5})$, where $C=M_{\mbox{\tiny NS}}/R$ is the compactness of the NS and $R$ is its radius. In terms of $y_{R}\equiv y(R)$, the Love number is given by

[TABLE]

where $y(r)$ is found through the solution of $r(dy/dr)+y^{2}+yF(r)+r^{2}Q(r)=0$ simultaneously with the TOV equations, with $F(r)=\{1-4\pi r^{2}[\epsilon(r)-p(r)]\}/f(r)$,

[TABLE]

and $f(r)=1-2m(r)/r$. The inputs for $\epsilon(r)$ and $p(r)$ are given by Eqs. (25)-(26), and the resulting neutron star mass for each radius $R$ is $M_{\mbox{\tiny NS}}=m(R)$.

In Fig. 6 we display the dimensionless TD of the single NS as a function of its mass, for those DD-vdW parametrizations predicting massive NS as found in Fig. 5. From the figure, one can notice that $\Lambda$ decreases nonlinearly as the NS mass increases. Furthermore, the TD of a canonical neutron star, $\Lambda_{1.4}$, can also provide a constraint on the EoS. The recent LIGO/Virgo detection of gravitational waves suggests values for $\Lambda_{1.4}$ inside the range of $70\leqslant\Lambda_{1.4}\leqslant 580$ (Abbot et al., 2018), which help us to test DD-vdW model. Our findings show that all parametrizations presented in Fig. 5 predict $\Lambda_{1.4}$ inside this range. Quantitatively, the maximum value of this quantity for the DD-vdW parametrizations analyzed is $\Lambda_{1.4}=527$ (for $K_{0}=251.5$ MeV and $J=35$ MeV), completely inside the limits from the LIGO/Virgo collaboration. We can also verify from Fig. 6 that $\Lambda_{1.4}$ increases as both $K_{0}$ or $J$ increase.

In Fig. 7, we show the dimensionless tidal deformabilities $\Lambda_{1}$ and $\Lambda_{2}$ for a binary NS system having primary mass $m_{1}$ and secondary mass $m_{2}$, respectively. At low orbital and gravitational frequency, the time evolution of the latter is determined by a defined combination given by ${\cal M}=(m_{1}m_{2})^{3/5}/(m_{1}+m_{2})^{1/5}$, with ${\cal M}$ being the chirp mass. In this work we fixed the chirp mass at ${\cal M}=1.188M_{\odot}$ and run $m_{1}$ in the range of $1.36\leqslant m_{1}/M_{\odot}\leqslant 1.60$ (corresponding to $1.17\leqslant m_{2}/M_{\odot}\leqslant 1.36$) according to the data from Ref. (Abbot et al., 2017, 2018). In the figure, we also present two lines corresponding to the 90% and 50% confidence limits obtained from the LIGO/Virgo collaboration coming from the analysis of the GW170817 event (Abbot et al., 2018). It is interesting to note that the tidal deformabilities for all analyzed DD-vdW parametrizations are inside the 90% credible region. This interesting result incentive us to perform further studies concerning our new developed DD-vdW model also for other hadronic environments.

4 Summary and Concluding Remarks

We developed a density-dependent van der Waals (DD-vdW) model by adopting the Carnahan-Starling method of excluded volume over the original vdW model. We also proposed a specific density dependence for the attractive part of the interaction by adding the $n$ power shown in Eq. (15). Such a new parameter was included as a generalization of the Clausius model proposed in Ref. (Vovchenko, 2017) to weaken the attraction and make our proposed model get closer to the Fermi free gas. We have shown in Figs. 2 and 3 that such a weakening is important in both cases, namely, to push the break of the causal limit to higher densities, and to impose agreement with the high-density behavior of the thermodynamical pressure established in Ref. (Danielewicz et al., 2002) (flow constraint). This constraint is significant and widely used to built or even to select relativistic hadronic models (Dutra et al., 2014; Kumar et al., 2018). Furthermore, the incompressibility at the saturation density, $K_{0}$, is also controlled by fixing $n$. In our model we can ensure values for $K_{0}$ in the range of $220\,\mbox{MeV}\leqslant K_{0}\leqslant 260~{}\mbox{MeV}$, according to the current consensus for the value of this quantity (Garg & Colo, 2018).

We also performed a generalization in the model in order to describe asymmetric systems. It was done by introducing a new term in the energy density with one more free parameter, adjusted to reproduce the symmetry energy at the saturation density. This new term can be seen as simulating the $\rho$ meson exchange between the finite nucleons (excluded volume included). Even with such a simple assumption, the model was shown to be compatible with some asymmetric nuclear matter constraints. Indeed, such a coupling can be improved, even with the inclusion of more free parameters, but we consider it quite suitable for a first approach.

Equations of state for neutron star matter under charge neutrality and $\beta$-equilibrium condition were also calculated with the DD-vdW model. It was shown that the mass-radius diagrams are in good agreement with the x-ray observations performed by Steiner (Steiner et al., 2010) and Nättilä (Nättiliä et al., 2016). The model is also compatible with the predictions for the binary neutron star system concerning the mass-radius diagrams obtained by the LIGO/Virgo collaboration (Abbot et al., 2018). We also verified that the model predicts massive stars in the range of $(1.97-2.07)M_{\odot}$, in agreement with observational data from Ref. (Antoniadis et al., 2013). Furthermore, we calculated the dimensionless tidal deformability of a single neutron star as a function of its mass and found a maximum value of $\Lambda_{1.4}=527$ for the canonical star. This value and other ones also obtained from the model are entirely consistent with the range of $70\leqslant\Lambda_{1.4}\leqslant 580$ recently obtained by the LIGO/Virgo collaboration (Abbot et al., 2018). Moreover, we have calculated the tidal deformabilities $\Lambda_{1}$ and $\Lambda_{2}$ related to a binary system and verified that the results obtained are fully compatible with the observational boundaries established from the analysis of the neutron star merger event GW170817 (Abbot et al., 2018).

In summary, the results obtained from the calculations performed with our new proposed DD-vdW model are in agreement with the available predictions for symmetric, asymmetric and stellar matter, including data from the recent GW170817 event.

Acknowledgments

This work is a part of the project INCT-FNA Proc. No. 464898/2014-5, partially supported by Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) under grants 310242/2017-7 and 406958/2018-1 (O. L.) and 433369/2018-3 (M.D.), by FAPESP under thematic projects 2014/26195-5 (M. B.), 2013/26258-4 (O. L.), and 2017/05660-0 (M. B., M. D., O. L.), and by the National key R&D Program of China under grant 2018YFA0404402. O. L. and M. D also thank J. Piekarewicz for fruitful discussions and useful suggestions. O. L. , M. D. and B. M. S. dedicate this paper to Antonio Delfino Jr., from Federal Fluminense University (UFF), Brazil, who passed away in July 2019.