Interplay between $\Delta$ Particles and Hyperons in Neutron Stars

TL;DR

This paper investigates how including $\Delta$ isobars in neutron star models affects their structure, composition, and observable properties, providing new constraints on their interactions and implications for gravitational wave observations.

Contribution

It introduces a relativistic mean field model incorporating $\Delta$ isobars that aligns with nuclear and astrophysical data, offering novel insights into neutron star composition and constraints on $\Delta$ interactions.

Findings

$\Delta$ isobars significantly influence neutron star composition and structure.

Stronger $\Delta$-meson interactions are required to support 2$M_{\odot}$ neutron stars.

Including $\Delta$ isobars leads to smaller radii and lower tidal deformability, consistent with GW170817 data.

Abstract

We analyze the effects of including $\Delta(1232)$ isobars in an equation of state (EoS) for cold, $\beta$-stable neutron star matter, employing relativistic nuclear mean field theory. The selected EoS reproduces the properties of nuclear matter and finite nuclei and, in the astrophysical context, allows for the presence of hyperons in neutron stars having masses larger than 2$M_{\odot}$. We find that the composition and structure of neutron stars is critically influenced by the addition of the $\Delta$ isobars, which allows us to constrain their interaction with the meson fields taking into account astrophysical information. Imposing that the EoS is stable and ensures the existence of 2$M_{\odot}$ neutron stars, as well as requiring agreement with data of $\Delta$ excitation in nuclei, we find that, in the absence of other mechanisms stiffening the EoS at high densities, the interaction of the $\Delta$ isobars with the sigma and omega meson fields must be at least 10\% stronger than that of the nucleons. Moreover, the neutron star moment of inertia turns out to be sensitive to the presence of $\Delta$ isobars, whereas the inclusion of $\Delta$ isobars in the EoS allows for smaller stellar radii and for a lower value of the tidal deformability consistent with the analysis of the GW170817 merger event.

Tables

Table 1. Table 1: Maximum mass, radius at maximum mass, central density, radius at M = 1.4 ​ M ⊙ 𝑀 1.4 subscript 𝑀 direct-product M=1.4M_{\odot} , and onset densities for the Δ − superscript Δ \Delta^{-} baryon and the Λ Λ \Lambda hyperon. The results are shown for the FSU2H EoS with hyperons but without Δ Δ \Delta ’s, and for the EoS with Δ Δ \Delta ’s for the cases x σ ​ Δ = x ω ​ Δ = 0.8 subscript 𝑥 𝜎 Δ subscript 𝑥 𝜔 Δ 0.8 x_{\sigma\Delta}=x_{\omega\Delta}=0.8 with x ρ ​ Δ = 1 subscript 𝑥 𝜌 Δ 1 x_{\rho\Delta}=1 , x σ ​ Δ = x ω ​ Δ = x ρ ​ Δ = 1 subscript 𝑥 𝜎 Δ subscript 𝑥 𝜔 Δ subscript 𝑥 𝜌 Δ 1 x_{\sigma\Delta}=x_{\omega\Delta}=x_{\rho\Delta}=1 , and x σ ​ Δ = x ω ​ Δ = 1.15 subscript 𝑥 𝜎 Δ subscript 𝑥 𝜔 Δ 1.15 x_{\sigma\Delta}=x_{\omega\Delta}=1.15 with x ρ ​ Δ = 1 , 2 , 0 subscript 𝑥 𝜌 Δ 1 2 0 x_{\rho\Delta}=1,2,0 .

Model	$M_{\text{max}}/M_{\odot }$	$R\left(M_{\text{max}}\right)$ (km)	$n_c$ (fm-3)	$R(1.4M_{\odot }$) (km)	Onset $n_{{\Delta }^{-}}\left({\mathrm{fm}}^{-3}\right)$	Onset $n_{\Lambda }\left({\mathrm{fm}}^{-3}\right)$
FSU2H	2.01	12.13	0.87	13.30	-	0.33
$x_{\sigma \Delta }=x_{\omega \Delta }=0.8,x_{\rho \Delta }=1$	1.75	12.40	0.68	13.03	0.28	0.41
$x_{\sigma \Delta }=x_{\omega \Delta }=x_{\rho \Delta }=1$	1.87	12.23	0.73	13.00	0.27	0.40
$x_{\sigma \Delta }=x_{\omega \Delta }=1.15,x_{\rho \Delta }=1$	2.00	11.87	0.85	12.97	0.26	0.39
$x_{\sigma \Delta }=x_{\omega \Delta }=1.15,x_{\rho \Delta }=2$	2.00	12.00	0.85	13.27	0.30	0.35
$x_{\sigma \Delta }=x_{\omega \Delta }=1.15,x_{\rho \Delta }=0$	2.00	11.76	0.85	12.62	0.21	0.41

Equations

$$\mathcal{L}$$

$${\mathcal{L}}_b$$

$${\mathcal{L}}_{\Delta }$$

$${\mathcal{L}}_l$$

$${\mathcal{L}}_m$$

$$ϵ$$

$$P=\sum _i{\mu }_i{n}_i-ϵ,$$

$${\mu }_b$$

$${\mu }_{\Delta }$$

$${\mu }_l$$

$$0$$

$${\mu }_i$$

$$n=\sum _{i=b,\Delta }{n}_i,$$

$$\frac{dP(r)}{dr}$$

$$\frac{dm(r)}{dr}$$

$$\lambda =\frac{2}{3}\frac{R^5}{G}{k}_2.$$

$$\tilde{\lambda }=\frac{1}{26}\left(\frac{M_1+12M_2}{M_1}{\lambda }_1+\frac{M_2+12M_1}{M_2}{\lambda }_2\right),$$

$$\tilde{\Lambda }=\frac{32}{G^4(M_1+M_2{)}^5}\tilde{\lambda }.$$

$$x_{\sigma \Delta }=\frac{g_{\sigma \Delta }}{g_{\sigma N}}{x}_{\omega \Delta }=\frac{g_{\omega \Delta }}{g_{\omega N}}{x}_{\rho \Delta }=\frac{g_{\rho \Delta }}{g_{\rho N}}.$$

$$0\lesssim x_{\sigma \Delta }-x_{\omega \Delta }\lesssim 0.2.$$

$$U_{\Delta }=-g_{\sigma \Delta }\bar{\sigma }+g_{\omega \Delta }\bar{\omega }+g_{\rho \Delta }{I}_{3\Delta }\bar{\rho },$$

Full text

dense matter – equation of state – stars: neutron – gravitational waves

1 Introduction

Neutron stars (NSs) are the most compact objects known without an event horizon (Shapiro & Teukolsky, 1983). Their cores contain strongly interacting matter that is several times denser than the matter at the center of nuclei. The properties of NSs primarily depend on the equation of state (EoS) of isospin-rich strongly interacting cold dense matter. Even though the EoS around saturation density can be determined from several experimental sources, mainly from the properties of nuclei and from the analysis of heavy-ion collisions, the large densities at the center of NSs have not been reached experimentally, and the characteristics of their innermost core remain largely unknown.

In order to model the EoS for predicting the structure and composition of matter in these extreme conditions, a variety of phenomenological theories, which take into account effective particle interactions, phase transitions and general relativity, have been developed, see e.g. Shapiro & Teukolsky (1983); Glendenning (1997); P. Haensel A.Y. Potekhin (2007). By providing unique constraints on the tidal deformability of NSs, the recent GW170817 detection by the LIGO-Virgo collaboration of gravitational waves from the merger of a binary NS system (Abbott et al., 2017, 2018, 2019) has stimulated further interest to explore the sensitivity of the NS properties on the EoS at large baryon density and isospin asymmetry. Since the tidal effects are strongly dependent on the stellar compactness, a measurement of the tidal deformability offers insights into the underlying EoS and the possibility to discriminate among EoSs that predict similar maximum NS masses but different NS radii.

The composition of a NS is driven by the so-called $\beta$-equilibrium condition, which establishes equilibrium among weak interaction processes, imposing charge neutrality and baryon number conservation. A NS consists mostly of neutrons, protons, and electrons. However, due to the large values of the nucleon chemical potential at the large inner densities, the conversion of nucleons to hyperons is energetically favorable. The appearance of hyperonic degrees of freedom relieves the Fermi pressure exerted by baryons and makes the EoS softer, which generally leads to a reduction of the maximum possible mass of NSs below $2M_{\odot}$, in disagreement with some observations (Demorest et al., 2010; Antoniadis et al., 2013), especially with the very recent value of $2.17^{+0.11}_{-0.10}$ solar masses reported for the millisecond pulsar J0740+6620 in Cromartie et al. (2019). Nevertheless, it is still possible to fine tune the models to fulfill the conditions imposed by recent observations and many theoretical studies regarding the presence of hyperons in the core of neutron stars have been conducted to achieve this goal, see Chatterjee & Vidaña (2016) and references therein. In this respect, relativistic mean field (RMF) models are especially powerful because they offer an effective and practical way to deal with the intricacies of the nuclear many-body problem at the high densities present in NSs (Glendenning, 1997; P. Haensel A.Y. Potekhin, 2007). These models rely on a few parameters—essentially meson-baryon coupling constants—that reproduce well-known properties of ordinary nuclei and assume plausible extrapolations to higher densities within a natural covariant framework.

Comparatively, less work has been done to study the presence of $\Delta$(1232) isobars in NSs, as early studies indicated that these particles would only appear at much higher densities than the typical ones in NS cores (Glendenning, 1985). A strong correlation between the appearance of $\Delta$ isobars and the $\rho$-meson coupling constant or, relatedly, the value of the symmetry energy in strongly interacting matter was pointed out. This type of study has been revisited in recent years (Drago et al., 2014b; Cai et al., 2015; Zhu et al., 2016; Li et al., 2018) in the light of the recent tight constraints obtained on the symmetry energy at saturation density and its slope, permitting a better extrapolation of the models to describe the dense and highly isospin-asymmetric conditions of NS matter. Indeed, the recent theoretical studies have shown that $\Delta$ particles may appear between 2 to 3 times nuclear matter saturation density, in competition with the onset of hyperons, and still have the potential to support $2M_{\odot}$ stars and to produce canonical $1.4M_{\odot}$ stars with relatively small radii below $\sim$13 km, as suggested by the latest astrophysical analyses, see e.g. Fortin et al. (2015); Ozel et al. (2016); Özel, Feryal and Freire, Paulo (2016); Nättilä et al. (2016); Lattimer & Prakash (2016); Abbott et al. (2018). The presence of $\Delta$ baryons may also lead to reduced tidal deformabilities (Li & Sedrakian, 2019) and bring them closer to the values measured from the gravitational waves emitted in the binary NS merger GW170817 (Abbott et al., 2017, 2018, 2019). The situation is, however, far from being settled. For instance, a recent study based on the quark-meson-coupling model, in which the $\Delta$-nucleon mass splitting emerges from a spin- and density-dependent one-gluon-exchange interaction, predicts a strongly repulsive $\Delta$ potential in dense matter, thus preventing the $\Delta$ isobars from appearing at all in NSs (Motta et al., 2019).

The aim of the present work is to analyze the interplay between $\Delta$ particles and hyperons in NSs and the implications for the properties of the star compared with the astrophysical observations, using the EoS of the recent fine-tuned RMF model FSU2H of Tolos et al. (2017a, b). This model allows for the presence of hyperons in NSs, while reproducing the properties of infinite nuclear matter and finite nuclei, the EoS constraints deduced from heavy-ion collisions, the symmetry energy and its density slope, and is compatible with the $2M_{\odot}$ observations as well as with relatively small radii for canonical NSs of $1.4M_{\odot}$ (also see Negreiros et al. (2018) for an application of the model to cooling simulations). We will show that the presence of $\Delta$ isobars, which competes with that of hyperons, does not compromise the good features of the model for reasonable values of the $\Delta$ coupling parameters. Indeed, we find that for compact stars that harbor strangeness in their dense interiors it is important to take into account the $\Delta$ degrees of freedom in the stellar composition in order to satisfy at the same time the observational constraints of massive NSs, small stellar radii, and the tidal deformability values obtained by LIGO-Virgo from the GW170817 event. We also find that the consideration of these astrophysical bounds favors a $\sim$10%–30% deeper potential for the $\Delta$ isobars than for the nucleon in nuclear matter at normal density, an aspect of the $\Delta$ interaction that thus far is not fully resolved from the laboratory experiments (Ericson & Weise, 1988; Horikawa et al., 1980; Nakamura et al., 2010; Lehr et al., 2000; O’Connell & Sealock, 1990). Finally, we have noticed that the moment of inertia of NSs is sensitive to the presence of the $\Delta$ baryons. The latter contribute to improve the agreement with recent EoS-independent constraints (Landry & Kumar, 2018) on the moment of inertia of the primary component of the double pulsar PSR J0737-3039.

The paper is organized as follows. In Section 2 we describe the main features of the model of the equation of state, whereas in Section 3 we summarize the calculation of the mass, radius and tidal deformability of NSs. In Section 4 we discuss our results and our conclusions are given in Section 5.

2 Equation of state

The RMF theory (Serot & Walecka, 1986) provides a Lorentz-covariant description of the microphysics of the EoS of dense matter (Glendenning, 1997). In this work we study the properties of NSs by computing the EoS with the hyperonic RMF model FSU2H that we introduced in Tolos et al. (2017a, b) supplemented by the $\Delta$(1232) degrees of freedom. Altogether, the model contains the octet of the spin-$1/2$ light baryons ($n$, $p$, $\Lambda$, $\Sigma^{-,0,+}$, $\Xi^{-,0}$), as in the FSU2H case, plus the quartet of the spin-$3/2$ $\Delta$ baryons ($\Delta^{-,0,+,++}$). In the RMF theory the baryons are represented in terms of Dirac spinors and the interactions among the baryons are described through the exchange of virtual mesons. The mesons of our model comprise a scalar-isoscalar $\sigma$ meson that provides the baryon-baryon attraction, a vector-isoscalar $\omega$ meson that provides the repulsion among baryons, a vector-isovector $\rho$ meson that accounts for the isospin dependence of the interaction, and a vector-isoscalar $\phi$ meson that couples to the hyperons. Including electrons and muons for the leptons present in cold catalyzed neutron star matter, the Lagrangian density can be written as

[TABLE]

where (we use units $\hbar=c=1$ throughout this work):

[TABLE]

The subscripts $b$ and $\Delta$ run over the baryon octet and the $\Delta$ quartet, respectively, whereas the subscript $l$ runs over the electrons and muons. The field $\Psi_{b}$ is the baryon Dirac field, $\Psi_{\Delta}$ is the Rarita-Schwinger field for the $\Delta$ isobars, and $\psi_{l}$ is the lepton Dirac field. The strong interaction couplings between the different baryons and mesons are denoted by $g$, whereas the $\vec{I}$ vector stands for the isospin operator. The $\mathcal{L}_{m}$ term contains the contribution from the free meson fields. The quantities $\Omega_{\mu\nu}$, $\vec{R}_{\mu\nu}$ and $P_{\mu\nu}$ denote the strength tensors of the $\omega$, $\rho$ and $\phi$ meson fields, respectively. The $\kappa$ and $\lambda$ couplings associated with the non-linear self-interactions of the $\sigma$ field (Boguta & Bodmer, 1977) are responsible for reducing the value of the nuclear matter incompressibility to within the experimental range. The main effect of the $\zeta$ self-coupling of the $\omega$ meson is to soften the EoS in the high-density sector (Mueller & Serot, 1996). Finally, the $\Lambda_{\omega}$ coupling of the quartic isoscalar-isovector interaction between the $\omega$ and $\rho$ mesons allows one to adjust the density dependence of the nuclear symmetry energy (Horowitz & Piekarewicz, 2001; Chen & Piekarewicz, 2014, 2015).

From the above Lagrangian one derives in the standard way the energy density of the system. Working in the mean-field approximation, where the meson fields are replaced by their ground-state expectation values, namely, $\bar{\sigma}=\langle\sigma\rangle$, $\bar{\omega}=\langle\omega^{0}\rangle$, $\bar{\rho}=\langle\rho^{0}_{3}\rangle$, and $\bar{\phi}=\langle\phi^{0}\rangle$ (under time-reversal symmetry and charge conservation, only the time components of the vector fields and the third component in isospin space of the $\vec{\rho}$ field contribute), the result is (Tolos et al., 2017a, b)

[TABLE]

The contributions $\sum_{b}\epsilon_{b}$, $\sum_{\Delta}\epsilon_{\Delta}$, and $\sum_{l}\epsilon_{l}$ from the baryons of the octet, the $\Delta$ isobars, and the leptons are obtained as analytical expressions of the Fermi momenta of the particles and their effective masses ($m_{b}^{*}=m_{b}-g_{\sigma b}\bar{\sigma}$, $m_{\Delta}^{*}=m_{\Delta}-g_{\sigma\Delta}\bar{\sigma}$, $m_{l}^{*}=m_{l}$), cf. Tolos et al. (2017a, b).

The pressure of the system can be obtained as

[TABLE]

where $n_{i}$ are the densities of the different particles and $\mu_{i}$ are their chemical potentials. The latter take the expressions

[TABLE]

where $k_{Fb}$, $k_{F\Delta}$ and $k_{Fl}$ are the Fermi momenta of the baryons, $\Delta$ isobars and leptons, respectively, and $I_{3}$ denotes the third component of the isospin operator, with the convention that for protons we have $I_{3p}=+1/2$. The conditions of charge neutrality and $\beta$-equilibrium of the matter in the NS core impose the following requirements on the particle densities $n_{i}$ and chemical potentials $\mu_{i}$:

[TABLE]

where $q_{i}$ and $b_{i}$ are the electrical charge and baryon number of particle $i$. For a given baryon density $n$ such that

[TABLE]

the equations of motion for the meson fields and the different particle species are solved self-consistently, under the restrictions presented in Eqs. (8)–(10), in order to obtain the chemical potential and the corresponding density of each of the species. One is thereby able to determine the composition of the stellar matter as a function of the baryonic density and the relation between the energy density and pressure (EoS) of the system.

For the nucleon and hyperon interactions we adopt the RMF parametrization FSU2H (Tolos et al., 2017a, b), with the values of the coupling constants and masses given in Table 1 of Tolos et al. (2017b). This model predicts realistic properties for ordinary nuclear matter, with a saturation density $n_{0}=0.1505$ fm*-3*, saturation energy per particle $E/A=-16.28$ MeV, incompressibility $K=238$ MeV, symmetry energy $E_{\textrm{sym}}=30.5$ MeV, and density slope of the symmetry energy $L=3n_{0}(\partial E_{\textrm{sym}}(n)/\partial n)|_{n_{0}}=44.5$ MeV. It also provides a good description of the ground-state properties of atomic nuclei (Tolos et al., 2017a, b). In the hyperonic sector, we have determined (Tolos et al., 2017a, b) the values of the couplings between the hyperons and the vector mesons using SU(3)-flavor symmetry, the vector dominance model and ideal mixing for the $\omega$ and $\phi$ mesons (Schaffner & Mishustin, 1996; Banik et al., 2014; Miyatsu et al., 2013; Weissenborn et al., 2012; Colucci & Sedrakian, 2013). In the case of the coupling of the $\Lambda$-hyperon to the $\phi$-meson, we have allowed for a reduction of 20% from its SU(3) value in order to reproduce the $\Lambda\Lambda$ bond energies (Ahn et al., 2013). Regarding the couplings of the hyperons to the scalar $\sigma$-meson, we have adjusted them for describing laboratory data on hypernuclei, in particular the hyperon optical potential. The values that we consider in the present work for the couplings $g_{\sigma\Delta}$, $g_{\omega\Delta}$, and $g_{\rho\Delta}$ between the $\Delta$ resonance and the mesons are discussed later in the section of results. As a final condition on the EoS with the $\Delta$ degrees of freedom, we will demand the resulting EoS to be stable (i.e., $dP/d\epsilon>0$ for all considered densities) and stiff enough as to ensure the existence of NSs with masses larger than 2$M_{\odot}$, in agreement with recent observations. These two restrictions will have a crucial influence on limiting the possible values for the $\Delta$ couplings.

The calculation of the structure and properties of a NS, such as the mass and the radius, requires the knowledge of the relation between pressure and energy density not only in the uniform matter of the liquid core but also in the solid crust of the star. The EoS for the inner crust with the FSU2H interaction has been recently computed in Providência et al. (2019). This allows us to describe with FSU2H the EoS for NSs in a unified way at the level of the core and the inner crust of the star. For the external layers of the crust, i.e., the outer crust, we employ the widely used EoS of Baym-Pethick-Sutherland (BPS) (Baym et al., 1971), which is strongly constrained by nuclear physics data.

3 Stellar structure

The structure of NSs is affected by general relativistic effects. Imposing hydrodynamic equilibrium of the star within General Relativity yields the set of differential equations known as Tolman-Oppenheimer-Volkoff (TOV) equations:

[TABLE]

where $G$ is the gravitational constant, $r$ is the distance to the center of the star and $m(r)$ is the mass enclosed within a sphere of radius $r$. Once the EoS of stellar matter is supplied, the TOV equations can be integrated from the origin with initial conditions $m(0)=0$ and an arbitrary value for the central energy density $\epsilon(0)=\epsilon_{c}$ until the pressure vanishes. The point $R$ where $P(R)=0$ defines the radius and total mass, $M=m(R)$, of the star.

Pressure plays a fundamental role in the determination of the structure of relativistic stars, as it is ultimately responsible for the existence of a limiting mass in these objects (Glendenning, 1997). The specific values of the mass and radius of NSs are linked to the details of the EoS, a fact that offers a way to test the validity of the different models (and to constrain their parameters), by requiring that they can reproduce the observed heavy NS masses of around 2$M_{\odot}$ (Demorest et al., 2010; Antoniadis et al., 2013) and predict radii below $\sim$13 km for canonical NSs of $1.4M_{\odot}$, as suggested by recent astrophysical extractions of NS radii (Fortin et al., 2015; Ozel et al., 2016; Nättilä et al., 2016; Özel, Feryal and Freire, Paulo, 2016; Lattimer & Prakash, 2016; Abbott et al., 2018). In particular, the maximum NS mass predicted by an EoS strongly depends on the stiffness of the core, which is known to become softer with the appearance of hyperons (Tolos et al., 2017a) and $\Delta$ isobars (Drago et al., 2014b).

In a coalescing binary system of two NSs, the component stars experience a quadrupole deformation due to the gravitational tidal field caused by their compact companion. The tidal deformation effect on a NS is characterized by the so-called tidal deformability $\lambda$, which is defined as the ratio between the induced quadrupole moment and the external tidal field (Flanagan & Hinderer, 2008). The tidal deformability $\lambda$ can be expressed in terms of the dimensionless second tidal Love number $k_{2}$ and the radius of the star as follows (Flanagan & Hinderer, 2008; Hinderer, 2008):

[TABLE]

The Love number $k_{2}$ is calculated by solving an additional differential equation self-consistently with the integration of the TOV equations (Hinderer, 2008; Hinderer et al., 2010; Postnikov et al., 2010). The global tidal effect of the two NSs in an inspiraling binary is described by the mass-weighted tidal deformability (Flanagan & Hinderer, 2008; Hinderer et al., 2010; Damour & Nagar, 2009):

[TABLE]

where $\lambda_{1}$ and $\lambda_{2}$ are the tidal deformabilities of the two component stars and $M_{1}$ and $M_{2}$ are their masses (if $M_{1}=M_{2}$, then $\tilde{\lambda}=\lambda_{1}=\lambda_{2}$). The same quantity can be expressed in dimensionless form as

[TABLE]

During the early stages of an inspiral, the phase of the gravitational wave signal is determined to leading order by $\tilde{\Lambda}$ (Flanagan & Hinderer, 2008; Favata, 2014). As the tidal deformabilities depend on the EoS—recall Eq. (13), measurements of $\tilde{\Lambda}$ from the gravitational waves emitted in a binary NS merger, such as the GW170817 event (Abbott et al., 2017, 2019), are a potential probe of the underlying EoS of the matter of compact stars.

4 Results and discussion

The EoS presented in Sec. 2 depends on the strong interaction couplings of the baryons to the meson fields (denoted by $g_{mb}$ with $m=\sigma,\omega,\rho,\phi$) and on the parameters of the self-interactions of these fields ($\kappa$, $\lambda$, $\zeta$ and $\Lambda_{\omega}$). We employ the values of the FSU2H model (Tolos et al., 2017a, b) (see Table 1 of Tolos et al. (2017b)), with the proviso that the incorporation of the $\Delta$ degrees of freedom introduces three new free parameters, i.e., $g_{\sigma\Delta}$, $g_{\omega\Delta}$ and $g_{\rho\Delta}$. As is customary, in the discussion that follows we will treat the values of these parameters in terms of the ratios to the meson-nucleon couplings:

[TABLE]

Even though the couplings of the $\Delta$ isobars with the meson fields are poorly constrained due to the limited existence of experimental observations, some information is available that may be used to restrict their values within some range. While most of the phenomenological analyses advocate for an attractive $\Delta$-nucleus potential, no consensus has been reached on its actual size. It has been claimed to be only a few tenths of MeV at the nuclear surface from pion-nucleus scattering reactions (Ericson & Weise, 1988; Horikawa et al., 1980; Nakamura et al., 2010), or practically the same as that of the nucleon from photo-absorption reactions (Lehr et al., 2000), or up to a 30% larger from inclusive electron-nucleus scattering data (O’Connell & Sealock, 1990). Moreover, as mentioned e.g. in Drago et al. (2014b), phenomenological analyses of electron-nucleus reactions in the region of the $\Delta$ excitation within a relativistic quantum hadrodynamics mean field model provided the following constraint (Wehrberger et al., 1989):

[TABLE]

Taking this phenomenology into consideration, we will explore how the composition and structure of NSs is affected by small modifications of the coupling of the $\Delta$ baryon to the $\sigma$ and $\omega$ mesons around values similar to those of the nucleon, i.e., taking $x_{\sigma\Delta}$ and $x_{\omega\Delta}$ within the range $[0.8-1.2]$. No information is available for the coupling of the $\rho$ meson to the $\Delta$ and, for the time being, we will adopt the value $x_{\rho\Delta}=1$, although the star properties for large variations of this unconstrained parameter will also be explored.

In Fig. 1, we study the dependence of the onset density (the lowest density at which a particle first appears) of the $\Delta^{-}$ and the $\Lambda$ hyperon in $\beta$-stable matter as a function of $x_{\omega\Delta}$, for different values of $x_{\sigma\Delta}$ and taking $x_{\rho\Delta}=1$. Among the four $\Delta$ particles, the $\Delta^{-}$ is the first one to appear since, being negatively charged, it can replace a neutron and an electron at the top of their Fermi seas. The $\Lambda$ is the hyperon that first appears, owing to its lowest mass and to the fact that the $\Sigma^{-}$ feels a repulsive interaction (Kohno et al., 2006).

For all chosen values of $x_{\sigma\Delta}$, namely $x_{\sigma\Delta}=0.8$ (dotted lines), $x_{\sigma\Delta}=1$ (dashed lines), and $x_{\sigma\Delta}=1.15$ (solid lines), we find that the $\Delta^{-}$ appears at smaller densities than the $\Lambda$ hyperon for the smaller values of $x_{\omega\Delta}$. As the value of $x_{\omega\Delta}$ increases, the onset density of the $\Delta^{-}$ also increases and that of the $\Lambda$ decreases, leading to the existence of a certain $x_{\omega\Delta}$, for the given $x_{\sigma\Delta}$, at which both the $\Delta^{-}$ and the $\Lambda$ appear at the same density. For larger $x_{\omega\Delta}$, the $\Lambda$ hyperon appears before the $\Delta^{-}$, at 0.336 fm*-3*, which corresponds to the onset density of the $\Lambda$ in the FSU2H model without $\Delta$ isobars (Tolos et al., 2017a). We note that, if the constraint of Eq. (17) is imposed, then we limit the range of validity of the curves in Fig. 1 to cases in which the $\Delta^{-}$ always appears first, delaying the appearance of the $\Lambda$ hyperon at a higher density than in the absence of $\Delta$ isobars, in agreement with the outcome of Drago et al. (2014b). Note also that the curves are not very much extended to the left, i.e., $x_{\sigma\Delta}-x_{\omega\Delta}$ may not be reaching its upper boundary of   0.2, because we encounter situations in which $dP/d\epsilon<0$ and the EoS becomes unstable.

The results of Fig. 1 can be easily understood from the dependence of the $\Delta$ baryon potential, defined as

[TABLE]

on the $\sigma$ and $\omega$ fields [see also Eq. (6)]. If $g_{\sigma\Delta}$ increases, the potential of the $\Delta^{-}$ becomes smaller, as the contribution of the $\sigma$ field is attractive, thus favoring the appearance of the $\Delta^{-}$ at smaller densities, while large values of $g_{\omega\Delta}$ increase the $\Delta^{-}$ potential because of the repulsive $\omega$ field, contributing to the opposite effect.

This is more clearly seen in Fig. 2, where the potentials of the nucleons ($U_{N}=-g_{\sigma N}\bar{\sigma}+g_{\omega N}\bar{\omega}+g_{\rho N}I_{3N}\bar{\rho}$) in pure neutron matter, shown in the upper left panel, are compared to those of the $\Delta$ particles ($\Delta^{-}$, $\Delta^{0}$, $\Delta^{+}$, and $\Delta^{++}$) in the other panels. The conditions of pure neutron matter resemble those found in $\beta$-stable NS matter. Results are shown for three combinations of the ratio parameters: $x_{\sigma\Delta}=x_{\omega\Delta}=0.8$ (upper right panel), $x_{\sigma\Delta}=x_{\omega\Delta}=1$ (lower left panel), and $x_{\sigma\Delta}=x_{\omega\Delta}=1.15$ (lower right panel). We have taken $x_{\rho\Delta}=1$ in all cases. As clearly seen in the upper left panel of Fig. 2, the nucleon potentials are negative for densities $n\lesssim 0.4$ fm*-3*, due to the dominance of the attractive $\sigma$ field in this region, and, conversely, they are positive at higher densities when the repulsive $\omega$ field overcomes the attractive $\sigma$ field. Note that the potentials of the different charged baryons in each isospin family differ due to the isovector term $g_{\rho}\bar{\rho}$ associated to the $\rho$ meson, which is negative in pure neutron matter. Therefore, the proton potential is smaller than that of the neutron. Similarly, when $x_{\sigma\Delta}=x_{\omega\Delta}=x_{\rho\Delta}=1$ (lower left panel of Fig. 2), out of the four charged $\Delta$ states, the $\Delta^{0}$ and $\Delta^{+}$ potentials coincide with those of the neutron and the proton, respectively, while the $\Delta^{-}$ state is the one that feels a higher potential and the $\Delta^{++}$ feels the lowest one. From Fig. 2 we observe that if $x_{\sigma\Delta}=x_{\omega\Delta}=0.8$ (upper right panel of Fig. 2), the $\Delta^{0}$ and $\Delta^{+}$ potentials are in size more moderate than the corresponding neutron and proton ones, being less attractive up to densities $n\sim 0.4$ fm*-3* and less repulsive beyond that point. In particular, the smaller repulsion at higher densities is the origin for the instabilities that are found in $\beta$-stable NS matter, as we will see. These instabilities are delayed in density if the $\Delta$ couplings are the same as the nucleon ones but not enough to sustain NSs with 2$M_{\odot}$. These massive NSs can be obtained when we take $x_{\sigma\Delta}=x_{\omega\Delta}=1.15$ (lower right panel of Fig. 2), which produces more attractive $\Delta$ potentials at lower densities and more repulsive ones at higher densities.

In Fig. 3 we compare the particle fractions in $\beta$-stable matter as functions of the baryon density for the FSU2H EoS without including $\Delta$ isobar degrees of freedom (first panel) and for $x_{\sigma\Delta}=x_{\omega\Delta}=0.8$ (second panel), 1 (third panel) and 1.15 (fourth panel), keeping in the latter three cases $x_{\rho\Delta}=1$. One clearly sees the delay in the appearance of the $\Lambda$ hyperon caused by the presence of the $\Delta^{-}$. As the couplings $x_{\sigma\Delta}=x_{\omega\Delta}$ increase, the $\Delta^{-}$ and $\Delta^{0}$ occur at a lower density since their corresponding potential is lower due to the fact that, in this density region, the size of the attractive $\sigma$ field is still larger than that of the repulsive $\omega$ field. Conversely, the appearance of the $\Delta^{+}$ and $\Delta^{++}$ takes place at larger density as $x_{\sigma\Delta}=x_{\omega\Delta}$ is increased, since these two $\Delta$ components enter at densities for which the repulsive $\omega$ field dominates over the attractive $\sigma$ field.

In general, our results are consistent with the findings discussed in previous works, although with some noticeable differences. In the case, for instance, of $x_{\sigma\Delta}=x_{\omega\Delta}=x_{\rho\Delta}=1$ we find an onset density for the $\Delta^{-}$ of 0.27 fm*-3*, somewhat larger than the value 0.22 fm*-3* reported in Li et al. (2018) and sensibly smaller than the value 0.4 fm*-3* quoted in Drago et al. (2014b). The origin of these differences must be analyzed from the condition $\mu_{n}+\mu_{e}=\mu_{\Delta^{-}}$, which, at the onset density, can be written as $g_{\rho\Delta}\bar{\rho}+\sqrt{k_{Fn}^{2}+m_{n}^{*^{2}}}+\mu_{e}=m_{\Delta^{-}}^{*}$. As discussed in Drago et al. (2014b), the fulfillment of this condition depends on a delicate balance between the value of the nuclear parameter $L$ (density slope of the symmetry energy), which controls the rise of the difference $\mu_{n}-\mu_{p}$ and, hence, of $\mu_{e}$, and the value of the term $g_{\rho\Delta}\bar{\rho}$, which is negative. A large $L$ value favors a lower onset density of the $\Delta^{-}$ as $\mu_{e}$ increases, but the negative field $\bar{\rho}$ will also be larger, with the opposing effect. One must also consider the particular value of the coupling strength $g_{\rho N}$. In general, models that have a larger $L$ value (implying a stiffer nuclear symmetry energy), also have a more moderate $g_{\rho N}$ value, see Tolos et al. (2017a), so that they may end producing $\Delta^{-}$ particles at relatively low densities, as is the case in Li et al. (2018). Models with lower $L$ values (softer nuclear symmetry energy), which usually come with larger $g_{\rho N}$ strengths, give rise to larger $\Delta^{-}$ onset densities, as found here and in Drago et al. (2014b). It is thus fair to say that the model dependence situates the appearance of the $\Delta^{-}$ in a density range $[0.2-0.4]$ fm*-3*, also in line to what is found in the hyperon-free model of Cai et al. (2015). It is also interesting to notice that the $\Delta^{-}$ partially replaces the role of the $\Sigma^{-}$ and $\Xi^{-}$ in compensating the positive charge of the protons, which, in the densities explored, leads to the absence of these hyperons in the composition of the star when the $\Delta$ isobars are considered. The shifting of the hyperon onset to higher densities when $\Delta$ isobars are present is a well-known fact (Drago et al., 2014a), but there is also a clear model dependence in this delay, as the models discussed in Drago et al. (2014b); Li et al. (2018) permit the appearance of the $\Xi^{-}$ below $n=0.8$ fm*-3*.

In Fig. 4 we show the dependence of the pressure versus the baryon density in the case of $\beta$-stable NS matter for the four models discussed before. With respect to the model without $\Delta$ isobars, we clearly see that the case $x_{\sigma\Delta}=x_{\omega\Delta}=0.8$ becomes unstable beyond four times nuclear matter density and therefore is unable to sustain NS masses of 2$M_{\odot}$, as we will see. The instability is somewhat delayed for the case $x_{\sigma\Delta}=x_{\omega\Delta}=1$, whereas the case $x_{\sigma\Delta}=x_{\omega\Delta}=1.15$ is well behaved for all densities, giving rise to a somewhat softer EoS at intermediate densities and a stiffer EoS at higher densities than the model without $\Delta$ isobars.

We now explore the effect of the parameter $x_{\rho\Delta}$, which has been so far kept to $x_{\rho\Delta}=1$. Since this parameter is not constrained by experiments, a few studies (Cai et al., 2015; Li et al., 2018) have varied its value within a sizable range. In Fig. 5 we show three EoSs that correspond to fixing $x_{\sigma\Delta}=x_{\omega\Delta}=1.15$ and taking $x_{\rho\Delta}=0$, 1, and 2. Since the $\rho$ field is negative, it makes the $\Delta^{-}$ less attractive, and hence larger values of $x_{\rho\Delta}$ postpone the appearance of the $\Delta^{-}$ state to larger densities. We observe that the change of $x_{\rho\Delta}$ between 0 and 2 affects the EoS mostly in the density range $[0.2-0.4]$ fm*-3*. In this density range, the EoS is softer if $x_{\rho\Delta}$ is lower and, consequently, intermediate size stars acquire a smaller radius as they are easier to compress. At high densities above $\sim$0.4 fm*-3*, the EoS, in contrast, is somewhat stiffer for a lower $x_{\rho\Delta}$ value. The compensating effect between the change of the pressure at densities below and above $\sim$0.4 fm*-3* leaves the maximum mass of heavy stars almost intact with the modification of $x_{\rho\Delta}$.

These effects are better seen in Fig. 6, where we compare the mass-radius relation obtained for the FSU2H EoS (without inclusion of $\Delta$’s) with those obtained taking the strength ratios $x_{\sigma\Delta}=x_{\omega\Delta}=1.15$, for which the 2$M_{\odot}$ limit is fulfilled, and $x_{\rho\Delta}=0$, 1 and 2. We observe that, as the value of $x_{\rho\Delta}$ is reduced, so is the radius of the star, an effect that is magnified for stars with masses around $1.5M_{\odot}$ reaching a reduction of 0.7 km, whereas the maximum mass remains at 2$M_{\odot}$ for the three $x_{\rho\Delta}$ values.

In Table 1 we summarize some NS properties obtained with the different models explored in the present paper. In the case of using the same meson couplings as those of the nucleons, the presence of $\Delta$ isobars decreases the maximum mass, as well as the radius of a canonical star of $M=1.4M_{\odot}$, with respect to the EoS without $\Delta$ isobars. Increasing the $x_{\sigma\Delta}$ and $x_{\omega\Delta}$ parameters so as to fulfill the maximum 2$M_{\odot}$ condition, gives rise to more compact stars, especially if a small value for $x_{\rho\Delta}$ is also considered, hence producing canonical NSs with realistic radii (Tolos et al., 2017a).

In Fig. 7 we summarize the constraints which the strength ratios $x_{\sigma\Delta}$ and $x_{\omega\Delta}$ must fulfill to ensure the existence of NSs compatible with the stability criterion and with the observational results. Stable EoSs are generated for parameters lying in the shaded blue region. The green region covers the configurations giving maximum NS masses higher than $2M_{\odot}$. The region enclosed between the red lines indicates the combinations of parameters that fulfill the experimental constraints for the difference between $x_{\sigma\Delta}$ and $x_{\omega\Delta}$ (Wehrberger et al., 1989), see Eq. (17). The dark brown region corresponds to the $x_{\sigma\Delta}$ and $x_{\omega\Delta}$ values for which all the constrains are satisfied. We notice that this region implies that the interaction between $\Delta$ isobars and the $\sigma$ and $\omega$ meson fields must be 10%–30% larger than in the case of nucleons. For $x_{\sigma\Delta}\simeq x_{\omega\Delta}$ this implies a larger (more negative) $\Delta$ potential than the nucleonic potential in symmetric nuclear matter at normal density. Thus, whereas to date the laboratory experiments are not conclusive on the size of the $\Delta$ potential with respect to the nucleonic one (Ericson & Weise, 1988; Horikawa et al., 1980; Nakamura et al., 2010; Lehr et al., 2000; O’Connell & Sealock, 1990), the astrophysical information appears to clearly advocate for a more attractive potential for the $\Delta$ than for the nucleon.

In the next figure, Fig. 8, we show our results for the mass-weighted tidal deformability $\tilde{\lambda}$ of a NS binary system [see Eq. (14)] as a function of the chirp mass $\mathcal{M}=\left(M_{1}M_{2}\right)^{3/5}/\left(M_{1}+M_{2}\right)^{1/5}$, where $M_{1}$ and $M_{2}$ are the masses of the two stars of the binary. The calculations are performed with our baseline EoS, i.e., the FSU2H EoS without $\Delta$ isobars, and with the FSU2H EoS including the $\Delta$ baryon with couplings $x_{\sigma\Delta}=x_{\omega\Delta}=1.15$ and $x_{\rho\Delta}=0$. We plot the predictions for the cases $M_{2}=0.7M_{1}$ and $M_{2}=M_{1}$, accounting for the mass asymmetry range compatible with the component masses $M_{1}\in(1.36,1.60)M_{\odot}$ and $M_{2}\in(1.16,1.36)M_{\odot}$ reported by the LIGO-Virgo collaboration for the GW170817 NS merger (Abbott et al., 2019). For comparison we also display in Fig. 8 the constraint on $\tilde{\lambda}$ from the GW170817 event at the measured chirp mass of $\mathcal{M}=1.186^{+0.001}_{-0.001}M_{\odot}$ (Abbott et al., 2019). The latest analysis of the GW170817 data by the LIGO-Virgo collaboration published in Abbott et al. (2019), for the case of low-spin priors (which are more consistent with the spin distribution of the population of galactic NS binaries), has obtained for this NS merger a dimensionless tidal deformability $\tilde{\Lambda}=300^{+420}_{-230}$ at 90% confidence level (the new result refines the upper bound $\tilde{\Lambda}\lesssim 800$ reported in the initial data analysis (Abbott et al., 2017)). Assuming $M_{2}=0.7M_{1}$ and the measured $\mathcal{M}$ value, Eq. (15) that relates $\tilde{\Lambda}$ with $\tilde{\lambda}$ leads to the constraint on $\tilde{\lambda}$ that we show by the vertical arrow in Fig. 8 (for $M_{2}=M_{1}$ the lower and upper values of this constraint on $\tilde{\lambda}$ are shifted downwards by 10%).

We can see in Fig. 8 that the result for $\tilde{\lambda}$ of the EoS with hyperons but without $\Delta$ isobars lies in the upper boundary of the GW170817 event. The inclusion of the $\Delta$ baryons has the effect of reducing the mass-weighted tidal deformability for all values of interest of the chirp mass. In particular, at the measured chirp mass $\mathcal{M}=1.186M_{\odot}$ the calculation of the tidal deformability with the EoS with $\Delta$ isobars is compatible with the GW170817 constraint on $\tilde{\lambda}$. The improvement is due to the fact that the incorporation of $\Delta$ isobars softens the EoS at $\sim$1.25–2.5 times the saturation density $n_{0}$. This softening is responsible for smaller radii (the Love number $k_{2}$ is also somewhat smaller with $\Delta$ isobars) and, hence, it produces smaller tidal deformability values. Our conclusion is in agreement with a recent study on the presence of $\Delta$ isobars using another type of density functionals that also include hyperons (Li & Sedrakian, 2019). Thus, the GW170817 event is consistent with a merger of a binary NS system having hyperons and $\Delta$ isobars in the stellar core.

We conclude by displaying in Fig. 9 the NS moment of inertia calculated from the EoSs without $\Delta$’s and with $\Delta$’s (with couplings $x_{\sigma\Delta}=x_{\omega\Delta}=1.15$ and $x_{\rho\Delta}=0$). The results have been obtained by solving Hartle’s slow-rotation differential equation for the moment of inertia in general relativity coupled to the TOV equations (Hartle, 1967). Astronomical observations of binary pulsars may provide information on the moment of inertia of NSs, ultimately offering possible limits on the underlying EoS (Lattimer & Schutz, 2005). The only double-pulsar system known to date is PSR J0737-3039. The mass of its primary component PSR J0737-3039A, or pulsar A, is of 1.338$M_{\odot}$. It is expected that a precise measurement regarding the moment of inertia of this slowly-rotating pulsar will be obtained in the near future from radio observations of PSR J0737-3039 (Burgay et al., 2003; Lyne et al., 2004). In a recent work, Landry & Kumar (2018) use approximately universal relations among NS observables (Yagi & Yunes, 2013a, b), known as the binary-Love and $I$-Love relations (recall that the tidal deformability is related to the Love number), to determine a range of $1.15_{-0.24}^{+0.38}\times 10^{45}$ g cm2 for $I$ of pulsar A from the 90% credible limits on the tidal deformability of a 1.4$M_{\odot}$ NS reported by LIGO-Virgo from the GW170817 event (Abbott et al., 2019). Landry & Kumar (2018) also facilitate a wider range $I\leq 1.67\times 10^{45}$ g cm2 for $I$ of pulsar A from the less restrictive upper limit on the tidal deformability obtained by the LIGO-Virgo collaboration in their initial data analysis of GW170817 in Abbott et al. (2017). We have plotted in Fig. 9 these two constraints reported in Landry & Kumar (2018) for $I$ of pulsar A (we show both of the constraints because there may be a certain dependence on the analysis and on the assumed boundaries for $\Lambda$ of a 1.4$M_{\odot}$ star). We can see in Fig. 9 that the moment of inertia is an observable highly sensitive to the existence of $\Delta$ isobars in the composition of the core of the star. The presence of these particles, in particular for smaller values of the coupling of the $\Delta$ to the $\rho$-meson, improves the agreement with the constraints on the moment of inertia of pulsar A, whose measurement is anticipated within the next few years.

5 Conclusions

We have studied the implications of the appearance of $\Delta$ isobars in the interior of NSs. Within the newly developed RMF model FSU2H (Tolos et al., 2017a, b) which incorporates hyperons, we have employed phenomenological analyses on electron-nucleus reactions in the region of the $\Delta$ excitation in order to constrain the couplings between the $\Delta$ isobars and the attractive $\sigma$ meson and the repulsive $\omega$ meson to the range $0\lesssim x_{\sigma\Delta}-x_{\omega\Delta}\lesssim 0.2$. In this manner, we have analyzed the composition, mass, radius, and tidal deformabilities of NSs as we incorporate the $\Delta$ degrees of freedom.

The presence of the $\Delta^{-}$ isobar in the composition of the star systematically postpones the emergence of the $\Lambda$ hyperon to higher densities. Only if one relaxes the constraint on the coupling of the $\Delta$ isobars to the $\sigma$ and $\omega$ mesons, the $\Lambda$ hyperon can appear for lower densities than the $\Delta^{-}$. In fact, the $\Delta^{-}$, $\Lambda$ and $\Delta^{0}$ baryons appear consecutively as the density of the stellar core increases, whereas the $\Delta^{+}$ and $\Delta^{++}$ baryons enter the composition at even higher densities. This is due to the different $\beta$-stability conditions to be fulfilled by the $\Delta$ charge states, together with the different charge-dependent single-particle potentials. It is to be noted that the $\Delta^{-}$ partially plays the role of the $\Sigma^{-}$ and $\Xi^{-}$ hyperons in compensating the positive charge of the proton, leading to the absence of these hyperons in the composition of the NS in our EoS. Moreover, the simultaneous increase of the coupling ratios $x_{\sigma\Delta}$ and $x_{\omega\Delta}$ gives rise to lower values of the potential for the $\Delta^{-}$ and the $\Delta^{0}$ states around densities of 2–$3\,n_{0}$, allowing for their earlier appearance in the composition, in contrast to a delay of the entrance of the $\Delta^{+}$ and $\Delta^{++}$ states, as their single-particle potentials become more repulsive for their onset densities above $3\,n_{0}$.

In order to fulfill the observations of massive NSs of $2M_{\odot}$ while having a stable solution for the EoS for all densities studied, we have chosen $x_{\sigma\Delta}=x_{\omega\Delta}=1.15$. Compared to the original EoS without $\Delta$ degrees of freedom, with these values we have obtained a softer EoS at intermediate densities that becomes stiffer at higher densities. As a consequence, NSs with smaller radii are obtained while still reproducing the $2M_{\odot}$ observations. Indeed, smaller values of the unknown $x_{\rho\Delta}$ coupling between the $\Delta$ isobars and the isovector $\rho$ meson lead to smaller stellar radii compared to those obtained in the original EoS. This is due to the extra softening of the EoS for a smaller $x_{\rho\Delta}$ coupling, as the $\Delta^{-}$ state becomes more attractive and appears at lower densities. These small values for radii favour smaller tidal deformabilities, more consistent with the value derived from the recent LIGO-Virgo gravitational-wave detection GW170817 accounting for the merger of two NSs (Abbott et al., 2017, 2018, 2019). We also find that the moment of inertia of slowly-rotating NSs displays a significant dependence on the presence of $\Delta$ isobars in the star’s interior, which leads to a reduction of the calculated values for the moment of inertia.

We have extended our study to determine the parameter space of $x_{\sigma\Delta}$ and $x_{\omega\Delta}$ values compatible with a stable EoS, the $2M_{\odot}$ observations and the phenomenological analyses on electron-nucleus reactions. We have found that the interaction between the $\Delta$ baryon and the $\sigma$ and $\omega$ fields must be 10–30% larger than in the case of nucleons, which, in turn, implies that the astrophysical information suggests that the $\Delta$ potential is more attractive than the nucleon potential at normal nuclear matter density.

We end noting that the triple requirement imposed by large maximum NS masses, small stellar radii and the low value of the tidal deformability derived from the GW170817 binary NS merger is a very demanding challenge for the nuclear models of the EoS, especially when the models include hyperons. In this work we have seen that the consideration of $\Delta$ degrees of freedom significantly contributes to help reconcile the theoretical EoS with all of the observational constraints.

L.T. acknowledges support from the FPA2016-81114-P Grant from Ministerio de Ciencia, Innovacion y Universidades, Heisenberg Programme of the Deutsche Forschungsgemeinschaft under the Project Nr. 383452331 and PHAROS COST Action CA16214. A.R., M.C. and C.G. acknowledge support from Grant No. FIS2017-87534-P from MINECO and the project MDM-2014-0369 of ICCUB (Unidad de Excelencia María de Maeztu) from MINECO. C.G. also acknowledges Grant BES-2015-074210 from MINECO.