Filling the Mass Gap: How Kilonova Observations Can Unveil the Nature of the Compact Object Merging with the Neutron Star

During the O1 and O2 observing runs, the LIGO Scientific Collaboration and Virgo Collaboration (LVC) detected gravitational wave (GW) signals from 10 coalescing stellar-mass black hole binaries (BHBH) and a neutron star binary system (NSNS), the latter accompanied by a multiwavelength electromagnetic (EM) counterpart (Abbott et al. 2017b, 2017c; The LIGO Scientific Collaboration & the Virgo Collaboration 2018a). At the time of writing, as the third observing run (O3) is in progress, the LVC reported the detection of two probable black hole–neutron star (BHNS) binary merger candidates (S190814bv—The LIGO Scientific Collaboration & the Virgo Collaboration 2019a, and S190910d—The LIGO Scientific Collaboration & the Virgo Collaboration 2019b), plus candidates with a lower probability of being actual astrophysical events.⁵ Before the beginning of O3, the estimated BHNS detection rate for this run was in the range of 0.04–12 yr⁻¹ (Dominik et al. 2015). At the time of writing, there are no indications of observed EM counterparts associated with these candidates (Coughlin et al. 2019; for S190814bv see, e.g., Klotz 2019; Soares-Santos 2019; Srivastav 2019, for S190910d see, e.g., Crisp 2019; Pereyra 2019).

From a theoretical point of view, BHNS mergers can be accompanied by an EM counterpart as in the NSNS case. This occurs when the NS is (at least partially) tidally disrupted before crossing the BH event horizon (Shibata & Taniguchi 2011). Tidal disruption is favored in binaries with low mass ratio q = M₁/M₂ and large NS tidal deformability Λ_NS (corresponding to a small NS mass and/or to a "stiff" equation of state). A high black hole spin,⁶ which brings the last stable circular orbit of the binary closer to the BH horizon, also greatly enhances the tidal disruption (Shibata & Taniguchi 2011; Kawaguchi et al. 2016; Foucart et al. 2018; Barbieri et al. 2019a, 2019b). The unbound NS material ("ejecta") is thought to produce kilonova emission (Lattimer & Schramm 1974; Li & Paczyński 1998; Metzger 2017). Moreover, Shapiro (2017), Paschalidis (2017), and Ruiz et al. (2018) showed that after a BHNS merger a relativistic jet can be launched, powering a short gamma-ray burst (sGRB; Eichler et al. 1989; Narayan et al. 1992) and GRB afterglow emission (Sari et al. 1998; D'Avanzo et al. 2018; Ghirlanda et al. 2018; Salafia et al. 2019).

The BH mass distribution observed so far in coalescing binaries is broad (The LIGO Scientific Collaboration et al. 2018b), extending up to ${50}_{-10.2}^{+16.6}\,{M}_{\odot }$, with the lightest BH carrying a mass ${7.6}_{-2.1}^{+1.3}\,{M}_{\odot }$, close to the mean BH mass in observed Galactic X-ray binaries of ∼7.8 ± 1.2 M_⊙ (Özel et al. 2010). Double NS systems observed so far carry masses in the interval 1.165 M_⊙–1.590 M_⊙ (Zhang et al. 2019), and the NS with the highest and best estimated mass is the radio pulsar J0740+6620 with ${M}_{\mathrm{NS}}={2.14}_{-0.09}^{+0.10}\,{M}_{\odot }$ in a low-mass binary (Cromartie et al. 2019). Thus, observations appear to indicate a discontinuity between the observed mass distributions of NSs and stellar BHs, called "mass gap," located approximately between ≈3 M_⊙ (the maximum NS mass inferred from causality arguments) and ∼5 M_⊙ (Lattimer & Prakash 2001). However, Thompson et al. (2019) recently reported the discovery of a BH with mass ${3.3}_{-0.7}^{+2.8}\,{M}_{\odot }$ in a noninteracting binary system with a red giant.

The mass spectrum of compact objects depends sensitively on the mass of the carbon–oxygen core at the end of stellar evolution, on the compactness of the collapsing core at bounce and on the supernova (SN) explosion engine. Belczynski et al. (2012) and Fryer et al. (2012) showed that, in the presence of a significant amount of fallback, explosions happening over a large interval of post-bounce times lead to a continuous range in remnant masses. By contrast explosions happening predominantly within a few hundreds of milliseconds after bounce, characterized by negligible amounts of fallback material, produce more easily the mass gap. Interestingly, at the time of writing, the LVC reported event candidates with binaries having at least one component in the mass gap (The LIGO Scientific Collaboration & the Virgo Collaboration 2019c, 2019d).

It is known that the binary chirp mass M_c, a combination of the masses of the two components, is one of the best measured parameters encoded in the GW signal. It is the prime parameter used to identify in low-latency searches the nature of the binary—whether the system hosts two NSs, two stellar BHs, or a BH and an NS. Interestingly, we note that if the NS and BH mass spectra join to form a continuum, i.e., there is no "mass gap" between BH and NS mass distributions (as Thompson et al. 2019 seem to indicate), there exists a range of values of the chirp mass M_c where the nature of the binary cannot be identified uniquely based on the chirp mass only (see also Mandel et al. 2015). In particular, hereafter we call "ambiguous" the chirp masses whose values are compatible with either an NSNS or a light BHNS system (see Figure 1).

Zoom In Zoom Out Reset image size

In this Letter, we aim at answering the following question: can EM observations of coalescing binaries in this "ambiguous" chirp mass interval help to disentangle their nature and narrow down the uncertainties on the existence of a mass gap? To this purpose, we study the properties of the kilonova emission of NSNS and BHNS systems which fall in this "ambiguous" chirp mass interval, using the semianalytical model presented in Barbieri et al. (2019a, 2019b).

The binary chirp mass is defined as

$$\begin{eqnarray}&&{M}_{{\rm{c}}}=\displaystyle \frac{{\left({M}_{1}{M}_{2}\right)}^{3/5}}{{\left({M}_{1}+{M}_{2}\right)}^{1/5}},\end{eqnarray} \tag{ 1 }$$

where M₁ and M₂ are the masses of the two component stars (we take M₁ ≥ M₂). LVC public alerts follow a classification scheme to communicate probabilistic estimates of the nature of the merging system to the community. The scheme classifies as "BNS" any system with both masses M₁ and M₂ smaller than 3 M_⊙; as "BBH" any system with both M₁ and M₂ larger than 5 M_⊙; as "NSBH" any system with M₁ > 5 M_⊙ and M₂ < 3 M_⊙, and as "MassGap" any system with at least one component carrying a mass between 3 and 5 M_⊙. An additional "Terrestrial" category is defined to represent triggers that are not of astrophysical origin (i.e., false alarms).

In this work we follow a slightly different classification. We assume the SFHo equation of state (EoS), for which the maximum mass of a nonrotating NS is ${M}_{\mathrm{NS}}^{\max }=2.058\,{M}_{\odot }$ (Steiner et al. 2013). We also fix the minimum NS mass to ${M}_{\mathrm{NS}}^{\min }=1\,{M}_{\odot }$ (∼10% lower than the value found in Suwa et al. 2018). We thus classify as "NSNS" those systems with both M₁ and M₂ between 1 and 2.058 M_⊙ (yellow region in Figure 1); "BHNS" those with M₁ > 5 M_⊙ and M₂ < 2.058 M_⊙ (green region); "BHBH" those with both masses above 5 M_⊙ as in the LVC classification (purple region). Considering that compact objects populating the mass gap have masses larger than ${M}_{\mathrm{NS}}^{\max },$ we assume these to be stellar-origin BHs. In Figure 1 we divide the "MassGap" region in three subregions: "BH+gap" for those systems with a BH above 5 M_⊙ and a BH in the gap; "gap+gap" for those with two BHs in the gap; "gap+NS" for those with a BH in the gap and an NS.

Two limiting values of the chirp mass can be identified: ${M}_{{\rm{c}},\min }^{\mathrm{gapNS}}=1.233\,{M}_{\odot }$ is the chirp mass corresponding to a gap+NS binary with ${M}_{\mathrm{NS}}={M}_{\mathrm{NS}}^{\min }$ and ${M}_{\mathrm{BH}}={M}_{\mathrm{NS}}^{\max }$ (red line). All GW events with chirp mass below ${M}_{{\rm{c}},\min }^{\mathrm{gapNS}}$ are NSNS mergers. Similarly, ${M}_{{\rm{c}},\max }^{\mathrm{NSNS}}=1.792\,{M}_{\odot }$ is the chirp mass corresponding to an NSNS binary with both NSs having the maximum allowed mass (blue line). Events with chirp mass above ${M}_{{\rm{c}},\max }^{\mathrm{NSNS}}$ cannot be produced by an NSNS merger. Events with chirp mass between ${M}_{{\rm{c}},\min }^{\mathrm{gapNS}}$ and ${M}_{{\rm{c}},\max }^{\mathrm{NSNS}}$ can be either NSNS or gap+NS mergers (green–orange lines), i.e., they are "ambiguous."

During the final phase of an NSNS merger, tidal forces lead to a partial disruption of the stars, producing an outflow of neutron-rich material. When the crusts of the two NS impact each other, compression, shock heating, and potentially neutrino ablation cause an additional outflow (Bauswein et al. 2013; Hotokezaka et al. 2013; Radice et al. 2016; Dietrich et al. 2017; Beloborodov et al. 2018). The released NS material can be divided into two components: the dynamical ejecta, gravitationally unbound, that leave the merger region, and a bound component, which forms an accretion disk around the merger remnant. On longer timescales, other outflows originate from the disk: faster ejecta produced by magnetic pressure and neutrino-matter interaction during the initial neutrino-cooling-dominated accretion phase (we call these "wind ejecta"), and slower but more massive ejecta produced by viscous processes in the disk, especially during the advection-dominated phase (Dessart et al. 2009; Metzger & Fernández 2014; Perego et al. 2014; Just et al. 2015; Siegel & Metzger 2017—we call these "secular ejecta"). Substantial differences in the ejecta properties arise depending on the postmerger scenario (see, i.e., Kawaguchi et al. 2019).

In order to calculate dynamical ejecta and disk mass from an NSNS merger we adopt the fitting formulae reported in Radice et al. (2018), which are calibrated on a large suite of high-resolution GRHD simulations.⁷ Both quantities depend on the NS masses and tidal deformabilities. We also adopt their formula for the dynamical ejecta mass-weighted average asymptotic velocity v_dyn.

The NS tidal disruption can occur also in BHNS mergers. If the NS is disrupted outside the innermost stable circular orbit, then the released material remains outside the BH in the form of a crescent (e.g., Kawaguchi et al. 2016), otherwise the NS plunges directly onto the BH. We adopt the fitting formula from Foucart et al. (2018) to calculate the total mass remaining outside the BH, M_out. This quantity depends on the BH mass and spin, and on the NS mass, tidal deformability and baryonic mass ${M}_{\mathrm{NS}}^{{\rm{b}}}$. We adopt the formulae in Kawaguchi et al. (2016) to calculate the dynamical ejecta mass and average velocity v_dyn in this case. M_dyn depends on the BH mass and spin, the NS mass, baryonic mass, and compactness C_NS, and on the angle ι_tilt between the binary total angular momentum and the BH spin. We assume ι_tilt = 0. v_dyn depends only on the mass ratio q = M_BH/M_NS. We proceed as in Barbieri et al. (2019b) to calculate C_NS from Λ_NS and ${M}_{\mathrm{NS}}^{{\rm{b}}}$ from M_NS and C_NS. We then calculate the disk mass as ${M}_{\mathrm{disk}}=\max [{M}_{\mathrm{out}}-{M}_{\mathrm{dyn}},0]$. As in Barbieri et al. (2019b) we assume that M_dyn cannot exceed 30% M_out, considering recent BHNS simulations presented in Foucart et al. (2019).

In what follows, we conservatively assume the BH to be nonspinning (χ_BH = 0), corresponding to the worst condition for ejecta production.⁸

In Figure 2 we show the dynamical ejecta and disk masses on the (M₁, M₂) plane along lines of constant M_c. We limit the y-axis to ${M}_{\mathrm{NS}}^{\max }$ as we focus on systems that contain at least one NS. It is apparent that NSNS configurations compatible with "ambiguous" chirp masses do not produce dynamical ejecta (upper panel of Figure 2). In this parameter region, the fits from Radice et al. (2018) predict the absence of this kind of ejecta, due to the prompt collapse of the remnant to a BH. Conversely, BHNS configurations can more easily produce dynamical ejecta. These systems have small mass ratios q < 5 and low-mass (large Λ_NS) NSs, which is the optimal condition to produce massive ejecta in BHNS mergers (as shown in Barbieri et al. 2019a, 2019b). The same arguments hold for disk masses (bottom panel in Figure 2). Note that the value of M_disk predicted by the fitting formula for NSNS systems in the considered range is set by the lower limit indicated in Radice et al. (2018), which is M_disk = 10⁻³. For BHNS configurations, instead, disks with masses up to ∼7 × 10⁻² M_⊙ are produced.

Zoom In Zoom Out Reset image size

It is important to note that when the differences of the ejecta mass for the BHNS and NSNS cases are substantial, they are larger than the systematic errors. The uncertainties on the fitting formulae for NSNS are ${\rm{\Delta }}{M}_{\mathrm{dyn}}^{\mathrm{NSNS}}=0.5\,{M}_{\mathrm{dyn}}+5\times {10}^{-5}\,{M}_{\odot }$ and ${\rm{\Delta }}{M}_{\mathrm{disk}}^{\mathrm{NSNS}}=0.5\,{M}_{\mathrm{disk}}+5\times {10}^{-4}\,{M}_{\odot }$ (Radice et al. 2018). The uncertainties on the fitting formulae for BHNS are ${\rm{\Delta }}{M}_{\mathrm{dyn}}^{\mathrm{BHNS}}=0.2\,{M}_{\mathrm{dyn}}$ (Kawaguchi et al. 2016) and ΔM_out = 0. 1 M_out (Foucart et al. 2018). Therefore, being ${M}_{\mathrm{disk}}^{\mathrm{BHNS}}={M}_{\mathrm{out}}-{M}_{\mathrm{dyn}}^{\mathrm{BHNS}}$, we assume its uncertainty to be ${\rm{\Delta }}{M}_{\mathrm{disk}}^{\mathrm{BHNS}}=\sqrt{{({\rm{\Delta }}{M}_{\mathrm{out}})}^{2}+{({\rm{\Delta }}{M}_{\mathrm{dyn}}^{\mathrm{BHNS}})}^{2}}$. We define ${\sigma }_{\mathrm{dyn}}=\sqrt{{({\rm{\Delta }}{M}_{\mathrm{dyn}}^{\mathrm{NSNS}})}^{2}+{({\rm{\Delta }}{M}_{\mathrm{dyn}}^{\mathrm{BHNS}})}^{2}}$ and ${\sigma }_{\mathrm{disk}}=\sqrt{{({\rm{\Delta }}{M}_{\mathrm{disk}}^{\mathrm{NSNS}})}^{2}+{({\rm{\Delta }}{M}_{\mathrm{disk}}^{\mathrm{BHNS}})}^{2}}$. We indicate, as the pink shadowed area in Figure 2, the regions where the differences in the mass of dynamical ejecta and disk for the BHNS and NSNS cases are greater than or equal to σ_dyn and σ_disk, respectively. In these regions the ejecta mass differences are larger than the systematic errors.

Figure 3 summarizes the differences between two representative NSNS and BHNS systems with "ambiguous" chirp masses (cases II and III in the figure), and also a "GW170817-like" NSNS case, for comparison. For the latter we consider an NSNS system with masses 1.46 M_⊙ and 1.27 M_⊙.

Zoom In Zoom Out Reset image size

Merger (I) produces relatively low-mass dynamical ejecta at all latitudes, with a preferentially equatorial angular distribution ∝ sin²θ, where θ is the polar angle (Perego et al. 2017). The accretion disk is massive and ∼20% of its mass is unbound in the form of secular ejecta, with a similar angular distribution as for the dynamical ones, while ∼5% of its mass goes into the wind ejecta, mostly confined in the polar region (θ < π/3 rad—Perego et al. 2017). After the merger, an intermediate state with a hyper-massive NS could exist before collapsing to a BH (gray central object represented in Figure 3(I)). The strong neutrino winds produced in this state interact with the ejecta, increasing the electron fraction Y_e or, equivalently, lowering the opacity.

We consider a system with 2 M_⊙ and 1.6 M_⊙ stars (II) as our representative NSNS merger in the "ambiguous" chirp mass range. As explained above, in this case we expect no dynamical ejecta and a low-mass accretion disk, resulting in low-mass wind and secular ejecta. The merger remnant collapses promptly to a BH. The absence of an intermediate hyper-massive NS state implies little neutrino wind, giving a low Y_e in the ejecta (Kawaguchi et al. 2019).

Finally, as BHNS merger in the "ambiguous" chirp mass range we consider a system with M_BH = 3 M_⊙ and M_NS = 1.1 M_⊙ (III). In BHNS mergers, the dynamical ejecta have a crescent-like shape, extending into half of the equatorial plane and limited to the region with θ < 0.3 rad (Kawaguchi et al. 2016). In the considered system, the dynamical ejecta and accretion disk are massive. Due to the lack of a neutrino wind, the fraction of accretion disk flowing into wind ejecta is lower than in the NSNS case (we assume ∼1%). The disk fraction that goes into secular ejecta is the same as in the NSNS case. As a consequence, the secular ejecta are massive, while the wind ejecta have low mass. The ejecta Y_e is lower than the NSNS case.

Therefore, with the ejecta properties being substantially different for the NSNS and BHNS cases in the "ambiguous" chirp mass range, we expect the kilonova light curves to present important differences as well.

The neutron-rich material ejected in NSNS and BHNS mergers is an ideal site for r-process nucleosynthesis, which produces the heaviest elements in the universe (Lattimer & Schramm 1974; Eichler et al. 1989; Korobkin et al. 2012; Wanajo et al. 2014). The synthesized nuclei are unstable and they decay radioactively, powering the kilonova emission (Li & Paczyński 1998; Metzger et al. 2010; Kasen et al. 2013).

We compute the kilonova light curves using the composite semianalytical model presented in Barbieri et al. (2019a, 2019b) (in part based on Grossman et al. 2014; Martin et al. 2015; Perego et al. 2017). For the NSNS cases we assume the model parameters (ejecta geometry, opacity, and the fractions of M_disk that go into wind and secular ejecta) as in Perego et al. (2017). The model has been tested on GW170817: using the parameters inferred for this event (Abbott et al. 2017a; Perego et al. 2017), we obtain multiwavelength light curves consistent with the observed ones (Villar et al. 2017, paper in preparation). For BHNS systems we assume the same model parameters as in Barbieri et al. (2019a, 2019b) (based on Just et al. 2015; Kawaguchi et al. 2016; Fernández et al. 2017).

The kilonova light curves are highly degenerate with respect to binary parameters. Thus, it is impossible to infer the system properties from the kilonova light curve alone. This degeneracy can be (at least partially) broken using information from GW analysis. In particular, the measurement of the binary chirp mass reduces the number of parameters by one. Leaving, i.e., M₁ as a free parameter, M₂ is constrained by the measured M_c.

In Figure 4 we show the envelope of the kilonova light curves expected from NSNS and BHNS mergers, for four selected values of the chirp mass. We consider emission in the g (509 nm) and K (2143 nm) bands and the figure shows the absolute magnitude as a function of time. For all the "ambiguous" chirp masses the fitting formulae from Radice et al. (2018) in NSNS mergers predict no dynamical ejecta and a minimum allowed disk mass M_disk = 10⁻³ M_⊙. Thus we have a single light curve for NSNS mergers, and we can expect that these events would not produce kilonovae brighter than shown in Figure 4.

Zoom In Zoom Out Reset image size

For BHNS mergers there exists a range of light curves for each M_c, arising from the different combinations of the component masses, producing different ejecta properties. For M_c = 1.45 M_⊙ (panels (1a)–(1b)) all kilonovae from BHNS mergers are much brighter at every time than that from NSNS mergers. Therefore a single observation in one of these bands would allow us in principle to distinguish the nature of the merging system.

At higher values of the chirp mass, there is only a small overlap between the BHNS and NSNS cases, at the bottom of the BHNS envelope. Therefore, except for observed light curves at low absolute magnitudes, it should be always possible to distinguish the nature of the merging system by the sole kilonova brightness. We note that the disentangling of the nature of the binary is optimal when M_c = 1.45 M_⊙ (panels (1a)–(1b)). In this case, as shown in Figure 2, the mass interval of the ejecta from BHNS mergers is the narrowest, and this in turn leads to the narrowest spread in the kilonova light curves.

The prediction of BHNS kilonovae as bright or brighter than NSNS ones is presented in Kawaguchi et al. (2019). They find that BHNS kilonovae are brighter in the near-infrared K-band, due to the smaller electron fraction Y_e in the ejecta owing to the lack of strong neutrino irradiation from the central remnant. In the i band, Kawaguchi et al. (2019) find that NSNS configurations ending with the formation of a supermassive NS leads to brighter kilonovae than the BHNS case. This is due to the strong neutrino emission produced in this case, that increases Y_e in the ejecta. However, in their study they compare sundry BHNS and NSNS configurations not selected on the bases of the chirp mass. By contrast, in our work, we compare BHNS and NSNS mergers at fixed chirp mass. This requirement restricts the NSNS binary configurations to cases producing no dynamical ejecta and very low mass disks. Therefore, whatever the value of Y_e in the ejecta from the NSNS merger is, the mass propelled in the merger is so low that almost all the BHNS kilonovae are brighter, at all wavelengths.

For other comparisons between NSNS and BHNS merger outcomes and studies on distinguishing the nature of merging compact binaries see Hinderer et al. (2019) and Coughlin & Dietrich (2019), who considered an unconventional BH companion with mass of ∼1.4 M_⊙, thus below the maximum NS mass.

As a visual comparison we also show the kilonovae for BHNS binaries having an NS with a "representative" mass of 1.4 M_⊙ (aqua/magenta lines). For M_c = 1.45 M_⊙ such a binary does not exist, while for M_c = 1.75 M_⊙ it is fated to a direct plunge, thus there is no kilonova.

We remark that the kilonova light curves from BHNS are inferred assuming nonspinning BHs (χ_BH = 0). As explained in Barbieri et al. (2019a, 2019b), increasing the BH spin (fixing all the other parameters) leads to more massive ejecta and, consequently, more luminous kilonovae. Therefore, if the BHs have a nonzero spin, our argument would be even stronger. As an example, for χ_BH ≳ 0.5 all light curves from BHNS kilonovae would be brighter than those from NSNS binaries in each band and at any time, in this critical range of "ambiguous" chirp masses.

The detection of a BHNS coalescence could be the next ground-breaking discovery in multimessenger astronomy. At the time of writing, there are promising GW candidates detected by the LVC during the observation run O3. The associated detection of an electromagnetic signal from these new GW sources might contribute to our understanding of the physical processes that power the multiwavelength EM emission (Gompertz et al. 2018; Rossi et al. 2019, and references therein).

From the GW signal, one of the best constrained parameters in low latency is the binary chirp mass, a combination of the masses of the two components. This parameter is currently used to classify the binary, whether the system hosts two NSs, two stellar BHs, or a BH and an NS. In the present Letter, we point out that in absence of a "mass gap" between the NS and BH mass distributions (as Thompson et al. 2019 seem to indicate), there exists a range of M_c (as shown in Figure 1) for which it is not possible to distinguish the nature of the binary on the basis of the chirp mass measurement alone.⁹ For the SFHo EoS adopted in this analysis, we find that the values of the chirp mass between 1.233 M_⊙ and 1.792 M_⊙ are compatible either with NSNS and BHNS systems.

In this Letter we show that the observation of the kilonova emission from these systems can break the degeneracy in the "ambiguous" chirp mass range, and constrain the nature of the merging system. We find that kilonova emission shows substantial differences in the luminosity and temporal evolution in NSNS and BHNS systems (see Figure 4). In particular, the BHNS case is far more luminous than the NSNS case, even when the BH is nonspinning. This happens because in this "ambiguous" M_c range the NSNS configurations represent the worst cases for ejecta production, while the BHNS configurations represent the best ones. It is important to note that, when the differences of the ejecta mass for the BHNS and NSNS cases are substantial, they are larger than the systematic errors in the modeling. Therefore, observing the kilonova associated with such an event is of fundamental importance to break the degeneracy on the nature of the merging system. Furthermore, if the observed kilonova is compatible with a BHNS merger, this would provide evidence in support of the existence of low-mass BHs, filling the "mass gap."

This work illustrates the potential of multimessenger observations of compact binary mergers, and the importance of an efficient exchange of information between the GW and EM communities (Biscoveanu et al. 2019; Margalit & Metzger 2019).

We thank F. Zappa and S. Bernuzzi for sharing EoS tables. The authors acknowledge support from INFN, under the Virgo-Prometeo initiative. O.S. and G.G. acknowledges the Italian Ministry for University and Research (MIUR) for funding through project grant 1.05.06.13.