Hyperons in Neutron Stars

: I review the issues related to the appearance of hyperons in neutron star matter, focusing in particular on the problem of the maximum mass supported by hyperonic equations of state. I discuss the general mechanism that leads to the formation of hyperons in the core of neutron stars and I review the main techniques and many-body methods used to construct an appropriate equation of state to describe the strongly interacting system of hadrons hosted in the core of neutron stars. I outline the consequences on the structure and internal composition of neutron stars and also discuss the possible signatures of the presence of hyperons in astrophysical dynamical systems like supernova explosions and binary neutron star mergers. Finally, I brieﬂy report about the possible important role played by hyperons in the transport properties of neutron star matter and on the consequences of neutron star cooling and gravitational wave instabilities induced by the presence of hyperons.


Introduction
The birth of hyperon physics dates back to 1947 when Rochester and Butler observed the first hyperon by studying cosmic rays in a gas chamber [1]. Some years later, in 1953, Danysz and Pniewski discovered the first hypernucleus [2] by the analysis of an emulsion experiment. Today, hypernuclear physics represents a very active research field whose main goal is to explain the features of systems that include hyperons starting from the underlying interactions that involve nucleons and hyperons. Unfortunately, several technical issues affect the experimental research about hyperons. The short lifetimes of hyperons and the low beam fluxes produced in scattering experiments do not allow constraints in the nucleon-hyperon (NY) interaction as much as the nucleon-nucleon (NN) one. NY scattering has been studied mainly in the p − Λ channel [3,4]; few data points exist in the p − Σ − channel [5,6]. No direct data scattering information is currently available in the hyperon-hyperon (YY) channel. Besides the scattering experiments, additional informations on hypernuclear physics is provided by the binding energies of hypernuclei, which are bound nuclear systems where one or more nucleons are replaced by hyperons. Today, more than 40 Λ-hypernuclei have been identified, while there is no clear evidence of the existence of Σ − hypernuclei. Ξ − hypernuclei have been observed in some events [7][8][9][10]. The study of hypernuclei represents a crucial step towards the understanding of hypernuclear interactions.
Besides the experimental and theoretical interest in the systems mentioned above, the study of hyperon physics is strongly correlated with the study of physics neutron stars. The latter poses a tough challenge that tests our ability to understand the behaviour of matter under extreme conditions of density and temperature. The huge density variation in the range ∼(1 − 10 15 ) g/cm 3 expected in neutron stars requires the modelling of systems in very different physical conditions, like heavy neutron-rich nuclei arranged to form a lattice structure, as in the outer crust of the star, or a neutral charge system of interacting hadrons and leptons forming a quantum fluid, as in the stellar core [11]. The description of such a variety of systems, of considerable interest for nuclear physics as well as for astrophysics, needs a challenging theoretical effort and an accurate knowledge of the interactions between the particles present inside the star. The bulk properties of neutron stars (e.g., gravitational mass, radius, moment of inertia, mass-shed frequency, tidal deformability, etc.) chiefly depend on the equations of state (EOS) describing the microscopic properties of stellar matter. The EOS of the dense matter is also one of the main inputs for the study of various astrophysical phenomena related to neutron stars like core-collapse supernovae (CCSN) [12] and binary neutron star mergers (BNSMs) [13][14][15]. The recent detection of gravitational waves (GWs) by the LIGO-Virgo collaboration [16][17][18][19], from several systems of merging neutron stars, has strongly increased interest in dense matter physics.
The behaviour of hyperons in a vacuum is very different from what occurs in the highdensity matter, as in the core of neutron stars. Simple thermodynamical considerations suggest that hyperons are produced in neutron star matter when the baryonic density becomes larger than a certain threshold, which depends on the features of the interaction between particles. If the baryonic density, namely the sum of the neutron and proton numerical densities, exceeds the threshold value, the system starts to convert some nucleons into hyperons. In addition, the huge density of the system and the action of the Pauliblocking mechanism prevent hyperon decay. Conversely, in a vacuum, hyperons are unstable and decay on a typical timescale set by the weak interaction. Typical processes that lead to the formation of Λ and Σ − hyperons are: p + e − → Λ + ν e and n + n → p + Σ − , where n and p stand for neutron and proton, while e − is the electron with corresponding electron neutrino ν e . In the following, I focus on cold neutron star matter; this condition is approximately realized after about 30 s from the star birth when the temperature in the core drops below 1 MeV [11]. In addition, on these timescales, the neutrino mean-free path becomes much larger than the star radius, and thus, in weak equilibrium reactions, the neutrinos chemical potential can be assumed to be zero. The hyperon formation processes mentioned above are energetically favourable in neutron star matter because they allow a decrease in both the total energy and pressure of the system. A well-known problem related to this scenario is that the softening of the EOS, induced by the formation of hyperons, leads to a quite low neutron star maximum mass, which may be inconsistent with observations. This issue is known in the literature as the "hyperon puzzle" of neutron stars. One of the main tasks of this work is to review the various mechanisms proposed to solve this issue.
According to various theoretical models, that will be discussed in this review, the hyperons expected to appear in neutron star matter are those contained in the baryonic octect: n, p, Λ, Σ − , Σ 0 , Σ + , Ξ 0 , Ξ − and, in addition, the ∆ quartet (∆ − , ∆ 0 , ∆ + , ∆ ++ ) of the decuplet. The specific population of hyperons according to a given model depends on the features of the interactions between nucleons and hyperons. Note that the present review is focused on the formation of hyperons in neutron star matter and it is, therefore, assumed that cold neutron stars are in general hyperonic stars. However, it should be stressed that other scenarios, which involve a quark deconfinement phase transition, cannot be currently excluded [20][21][22][23][24][25]. Another important point that should be clear is that such alternative scenarios are not incompatible with the presence of hyperons in neutron stars [26].
The paper is organised as follows, I first discuss the many-body methods and corresponding interactions used to derive the hyperonic matter EOS employed for this task. Specifically, I first discuss relativistic mean field approaches and then move to the nonrelativistic microscopic ones. In both cases, I show the main features about the resulting neutron star structure and composition when hyperons are included in the system. In the last part of this review, I report the possible signatures of the appearance of hyperons in astrophysical systems like BNSMs and SNe, as well as some possible important roles played by the presence of hyperons in hot young neutron stars.

Many-Body Methods
The first study of hyperons in neutron stars dates back to the work of Ambartsumyan and Saakyan in 1960 [27]. Later on, this research was considered by several other authors, who employed various many-body methods to carry out their investigations. These different techniques can be divided into two main categories: phenomenological and microscopic approaches. The first ones are mainly based on relativistic mean-field (RMF) theory, while the second is built in a non-relativistic framework, where hadrons interact with each other through a potential. Among the phenomenological approaches, some authors have also tried to implement non-relativistic Skyrme interactions [28] that have been widely used in the pure nucleonic sector to describe both finite nuclei and infinite nuclear matter. In the following, I review the main features of these methods and discuss their predictions about the structure of neutron stars.

Relativistic Mean-Field Models
RMF models are based on a relativistic Lagrangian density inspired by quantum field theory. Hadron fields are assumed to interact with mesons by a minimal Yukawa coupling. In this approach, all baryon and meson fields are treated as classical fields, and the various operators entering the Lagrangian density are replaced by their expectation values calculated in the nuclear medium. A typical Walecka-type Lagrangian density can be written as: where σ * is the baryon effective mass. Ψ B and ψ l are the baryon and lepton Dirac fields, respectively, and σ, ω, and ρ represent the scalar, vector, and isovector meson fields, which describe the nuclear interaction. The two additional hidden strangeness mesons σ * and φ do not couple with nucleons using SU(6) symmetry (g σ * N = 0, g φN = 0). The scalar σ * meson corresponds to the physical f 0 (975) meson, while the vector φ one corresponds to the φ(1020) (meson mass is reported in parenthesis). The φ meson is interpreted at the quark level as an ss state (being s (s), the strange quark (antiquark)), while there is no unique interpretation for the σ * meson. The σ * and φ mesons were originally introduced in Ref. [29], to solve the problem of the attractive ΛΛ bond (see discussion below for additional details and definitions). Lepton masses and bare baryon masses are denoted by m l and m B , respectively. The coupling constants of i-meson (where i = σ, ω, ρ) with a baryon B are denoted by g i,B , where the index B runs over the eight baryons of the baryonic octect: n, p, Λ, Σ − , Σ 0 , Σ + , Ξ − , and Ξ 0 . The couplings g i,B set the strength of the interaction between two baryons and a meson. The sum on l is over electrons (e − ) and muons (µ − ). The ∆-quartet can also be included in this formalism as discussed later. For simplicity, in (1), I set g i,N ≡ g i for nucleons.
Λ ω is a coupling constant of the mixed ρω-term [30] introduced to adjust the behaviour of symmetry energy near saturation density (n 0 = 0.16 fm −3 ). The coupling of Λ N allows us to finally consider a mixing between the σ-and ρ-meson. The constants k and ξ are the weights of the σ and ω self-interaction terms and τ B is the isospin operator, namely the vector of the three Pauli matrices. The mesonic field tensors are given by their usual antisymmetric expressions: In this notation, the b-field is associated to the ρ-meson.
Note that, so far, various classes of RMF models have been developed according to the terms retained in the Lagrangian density (1). In some models, the couplings g σ , g ω , and g ρ have a density dependence, and where terms proportional to k, λ, ξ, Λ ω , and Λ N are set to zero, while in other models g σ , g ω , and g ρ are constant, and all terms of the above Lagrangian are kept. In the following, I refer to the first class of models as densitydependent (DD) mean-field models and to the seconds as non-linear (NL) mean-field models. In the following, I discuss these two kinds of RMF models, considering them as representative cases and then briefly report other possible choices.
In the case of DD models, the density dependence of the couplings is usually taken in the form: with x = n/n 0 and: for the isoscalar coupling and: for the isovector couplings. The values of the parameters a m , b m , c m , d m , and the density dependence in the above functions were originally proposed in Ref. [31] in order to reproduce the density behaviour of the σ-, ω-, and ρ-couplings found in Dirac-Brueckner-Hartree-Fock calculations [32,33]. More in general, the meson-nucleon coupling constants can be fixed in order to reproduce nuclear experimental observables like the saturation properties of symmetric nuclear matter [20], finite nuclei properties (binding energies, radii, etc.), constraints from heavy ion collision experiments, or even neutron star properties [31,[34][35][36][37][38][39][40].
Concerning the meson-hyperon and the strange meson-hyperon coupling constants g ωH , g ρH , g σH , g σ * H , and g φH , one of the first models that attempted to determine these meson-hyperon couplings was proposed by Glendenning [41], using Λ hypernuclear data to fix the σ − Λ coupling and assuming simple rescaling relations between the couplings between the other mesons and the hyperons. I note that in this model, the presence of the σ * , and φ mesons was not considered. Later on, other ways to fix the meson-hyperon couplings were proposed; one popular choice assumes that the ω-and ρ-hyperon couplings are determined by SU(6) symmetry relations: The coupling constants g σH of the hyperons with the scalar meson σ are then adjusted to the potential depths U H felt by a hyperon H in symmetric nuclear matter at saturation density which is given by the relation: In the above expression ω 0 and σ 0 are the values of the fields ω and σ in symmetric nuclear matter at saturation density. According to Equations (5)-(7), all hyperon couplings g σH , g ωH , g ρH H=Λ,Σ,Ξ are finally determined once the coupling constants g σ , g ω , g ρ of the nucleon sector are provided, together with the values of U H . Concerning the hyperon-single particle potential depths in symmetric nuclear matter at saturation density, typical values used are in the following ranges: U Λ = (−28; −32) MeV, U Σ = (−30; 30) MeV, and U Ξ = (−18; 18) MeV. Note that apart from the Λ-hyperon, whose value of U Λ can be obtained from the extrapolation of the separation energy of Λ-hypernuclei in the limit of infinite mass number [42], the other hyperons of the experimental scenario is quite uncertain. For the Σ-hyperon, no information from hypernuclei is currently available; some experiments on Σ − -atoms and inclusive (π − , K + ) [43][44][45][46] spectra on medium to heavy targets seem to indicate a repulsive Σ potential. For the Ξ-hyperon there are some experimental data from the hypernuclei: 15 Ξ − C and 12 Ξ − Be. The latter was produced in (K − , K + ) reaction on a 12 C target [7], while the first, known as the Kiso event [8], corresponds to an intermediate state in the reaction Ξ − + 14 N → 15 Ξ − C → 10 Λ Be + 5 Λ He. Both sets of data indicate an attractive Ξ single particle potential compatible with a 1p state, and consequently bound Ξ-hypernuclei. More recently, in an emulsion-counter hybrid experiment performed at J-PARC [9], they observed a Ξ − absorption event on 14 N, which decayed into twin single-Λ hypernuclei. Furthermore, in this case, the Ξ − energy level was compatible with a 1p state. Finally, in Ref. [10], they reported evidence for a 15 Ξ − C hypernucleus with the Ξ − hyperon in a 1s state. The experimental scenario is even worse in the hyperon-hyperon (YY) sector where the only informations come from double-Λhypernuclei and no direct data scattering is available. Such uncertainties clearly reflect the setting of YY-meson couplings. In order to set the σ * -meson couplings, a popular choice adopted is to fix it by considering the ΛΛ coupling according to the bond energy in double ΛΛ hypernuclei: . ∆B ΛΛ can thus be determined once the double-and single-Λ binding energies are experimentally measured. Old experiments of Λ emulsion suggested a value of ∆B ΛΛ ∼ 5 MeV [47,48]. One should, however, note that the identification of some double-Λ hypernuclei in these experiments was ambiguous and, therefore, some caution should be used when employing the data of these experiments to constrain the ΛΛ interaction. More recently, at KEK [49], the 6 ΛΛ He double-Λ hypernucleus was identified without ambiguity and a value of ∆B ΛΛ ∼ 1.1 ± 0.2 +0. 18 −0.11 MeV was derived. This last value has been updated to ∆B ΛΛ ∼ 0.67 ± 0.17 MeV, considering the variation of the Ξ-mass [50].
The g σ * Λ is related to ∆B ΛΛ by: where U Λ Λ is the Λ potential depth in Λ matter at saturation density, and like in Equation (7), the label 0 indicates the mean-field value of the meson fields. Note that both Equations (7) and (8) are correct in the case of RMF models without a density dependence in the coupling constants. For the general case of DD RMF models, one should account for extra terms involving the derivative of the couplings (see for instance Ref. [51]). Once the value of g σ * Λ has been obtained by the above procedure, the use of SU(6) fixes the couplings of the σ * meson with the other hyperons: For the φ-meson couplings, using again SU(6) symmetry relations one has: This setting of YY-meson couplings has been implemented, for instance, in Ref. [36] (assuming B ΛΛ = 5 MeV), using as base model for the nucleonic sector of both the TM1 RMF parametrisation of the NL Walecka model and a modification (TM1-2 model) proposed again in Ref. [36], where the nucleon-meson couplings were adjusted to better fit the flow constraint in heavy-ion-collision experiments by Danielewicz [52]. I remark that, besides the representative example reported above, several other possibilities for the determination of the NY-meson and the strange meson-hyperon couplings have been proposed. In Ref. [51], the NY-meson couplings where calibrated on the separation energies of Λ-hypernuclei. For a fixed value of the ω-Λ coupling taken both from SU(6) symmetry and by considering an extreme value for which the symmetry is broken, the authors determined the σ-Λ coupling by minimising the spread between the theoretical prediction given by the RMF calculations on several hypernuclei, and the corresponding experimental value. Concerning the value of the coupling constants of the Λ-hyperon to the hiddenstrangeness mesons σ * and φ, they were fixed using the measured ΛΛ bond energy of 6 He ΛΛ by KEK. From this analysis, it turned out that the usual value used in the literature of U Λ Λ = −5 MeV was not consistent with hypernuclear data. I note that the resulting value of U Λ varied in the range of −30 to −35 MeV, according to the specific parametrisation adopted for the nuclear matter Lagrangian density. In Ref. [53], the σ-Λ coupling was fixed using hypernuclear data from a single Λ-hypernuclei, while an upper bound on the σ-Σ coupling was obtained from the requirement that the lower bound on the maximum mass of a neutron star be around 2M . It is important to note that in some works, different versions of the Lagrangian density (1) have been proposed. In Ref. [54], a nonlinear chiral SU(3) model was developed to study the generation of baryon mass. Later, the model was extended to study finite nuclei [55] and neutron stars [56]. In addition, the model was also extended to finite temperature with the inclusion of the ∆ quartet [57]. In this model, baryons interact through the usual one-meson-exchange with the inclusion of an effective δ-meson plus a dilaton field χ, which was included to mimic the effect of a gluon condensate, and a ζ-meson that corresponds to the strange quark condensate. The inclusion of the δ-meson turned out to be important to achieve a good description of asymmetric nuclear matter and, more generally, leads to a larger repulsion in dense neutron-rich matter and to a definite splitting of proton and neutron effective masses [58][59][60]. It is interesting to note that in hyperonic matter, the effect is the opposite, namely, there is a softening of the equation of state due to the effect of the δ meson [61].
Another RMF model employed in the literature, is the quark-meson coupling (QMC) [62,63]. In this model, baryons are described in the framework of the MIT bag model and the properties of nuclear matter can be self-consistently calculated by the coupling of scalar (σ) and vector (ω) fields to the quarks within the nucleons, rather than to the nucleons themselves. As a result of the scalar coupling to the quarks, the internal structure of the baryon is modified with respect to the free case.
Another ingredient that has been considered in the Lagrangian density (1) is the possibility to have a kaon and/or other meson condensation in neutron star matter. The possibility of pion condesation was proposed independently by Migdal [64] and Sawyer [65], pointing out that pions could condense in dense matter due to attractive p-wave interactions with nucleons. Later, several authors contributed to this research [66][67][68][69] with contradictory conclusions. As stressed in Ref. [69], the pion condensation in neutron star matter should be revisited using modern experimental data in the pion-nucleon interaction, as well as most advanced many-body techniques. Concerning the kaon condensation, if the K − meson develops sufficient attraction in dense matter, it could be energetically more favorable, above some critical density, to neutralise the positive charge with antikaons rather than with electrons. Such possibility has been proposed by Kaplan and Nelson [70] and investigated extensively by several authors [71][72][73][74][75][76][77][78][79][80][81]. All these works agree that the kaon condensation leads to a softening of the hadronic EOS of neutron star matter. In addition, Ref. [82] argued that the presence of hyperons may shift the onset of kaon condensation to very large densities, not expected to be reached in neutron stars. I remark that there are still uncertainties in theKN (K = K − , K 0 ) interaction, especially about the so-called sigma term (Σ KN ), which describes the kaon-nucleon interaction as s-wave. The sigma term is poorly known and, unfortunately, it turns out to be one of the main ingredients to determine the onset of kaon condensation in neutron stars. This issue is one of the main challenges for theoretical investigations. New studies on the kaon-atom system carried out on at Siddharta experiment [83] will probably improve the status of this research.
The ρ condensation in neutron stars has been considered in Refs. [84,85] assuming Brown-Rho-type hadronic mass modifications at high density [86]. The Brown-Rho conjecture is based on the assumption that hadrons (excluding the pseudo-Goldstone bosons like pions) experience a mass reduction in nuclear matter, which is proportional to the in-medium quark condensate. The idea behind the ρ − -meson condensation is quite similar to the one that concerns the kaon condensation-at some density, the role played by electrons should be replaced by ρ − -mesons for generating a charge-neutral system. In Ref. [87], they discussed the effects of a magnetic field on the ρ-condensation; the authors concluded that the magnetic field has a relevant effect in triggering such processes.
Note that some studies point out the important role that may be played by the ∆baryon in neutron star matter, especially considering the shift they could produce of the hyperon onset. The ∆-baryon is part of the baryon decuplet and is in a isospin quartet state; old calculations [88] suggested that ∆-baryons do not appear in cold β-stable neutron star matter. Recently, the issue of ∆ formation in cold matter has been reconsidered in several works [89][90][91][92][93][94][95][96][97], where they have critically discussed the influence and uncertainties of the meson-baryon couplings in this specific channel. In Ref. [89], according to a NL RMF model, it was found that the presence of the ∆-baryon shifts the onset of the other hyperons to a larger density and accordingly gives rise to a very strong softening of the EOS of the system. The authors concluded that this mechanism may trigger a quark deconfinement phase transition. The authors of Ref. [95] essentially reached the same conclusion using a different hadronic model. Finally, in Ref. [96], using the quark meson coupling model (QMC), it was found that the appearance of ∆ baryons is suppressed by strongly repulsive many-body forces that modify the internal structure of baryons at large densities. Unfortunately, just one microscopic calculation [98], based on the microscopic Brueckner-Hartree-Fock (BHF) approach currently exists on the single particle potential of ∆-baryon in symmetric nuclear matter. In addition, no β-stable calculation in a microscopic framework has ever been carried out. These kind of calculations, based on a different many-body approaches, are important to independently test the assumptions made in the RMF theory and to benchmark such results.
In the Hartree approximation, the equations of motion of the various fields can be easily obtained starting from a typical Lagrangian density, like (1), and the Euler-Lagrange equations. In this approximation, all of the fields retain only a density dependence. The total set of equations is therefore reduced to a set of nonlinear ordinary equations, which can be solved numerically with a multidimensional Newton-Raphson method. I note, however, that some works have been successfully considered, such as the inclusion of the Fock-term in the equations of motion [99][100][101]. In this case, the solution of the resulting system of equations is much more complex and requires the solution of several coupled integral equations. In cold β-stable hyperonic matter for a given baryonic density, the full set of equations is closed by weak equilibrium relations between the chemical potentials of hadrons and leptons, and by the charge neutrality condition. A final remark that to make is that almost all RMF approaches employed in the literature assume a phenomenological one-meson exchange model in the Lagrangian density. This hypothesis regarding the form of the interaction is mainly dictated by the request of simplicity rather than a real physical motivation. Recently, in Ref. [102], they developed an extension of RMF theory beyond the one-meson exchange approximation, including contributions from two-loop diagrams arising from the exchange of isoscalar and isovector mesons between nucleons. In the future, it would be very interesting if similar calculations could be extended to the hyperonic sector.

Neutron Star Structure of RMF Models
In this section, I consider representative RMF models described in the previous section, to show results about the mass-radius relation obtained, solving the Toleman-Oppheneimer-Volkov (TOV) equations and the corresponding β-stable matter composition. In Figure 1, I report several mass-radius curves together with the observed pulsars PSR J0348+0432 with mass 2.01 ± 0.04 M [103] and PSR J0740+6620 [104], which are represented by yellow and purple bands, respectively, indicating the uncertainty on the measurement. Recently, the mass of PSR J0740+6620 has been updated to 2.08 ± 0.07 M in Ref. [105]. Together with the PSR J1614-2230 with mass 1.97 ± 0.04 M [106], these pulsars represent the largest neutron star masses ever observed. It is interesting to note that the error bars on these measurement is very small. This has been possible by the use of the Shapiro delay effect [107]. I note that other indirect constraints on the maximum mass of neutron stars came very recently from the analysis of the BNSM event GW170817 [108][109][110]. Previous works agree that the maximum mass of neutron stars cannot be larger than ∼2.2 M in perfect agreement with direct observations. A very tight constrain on the EOS of neutron star matter would be provided by the combined mass-radius measurements for a series of neutron stars. Such very difficult tasks have been considered by several authors [111][112][113][114]. In Ref. [113], using spectroscopic techniques, a constraint of R > 11.1 km was placed on the radius on the millisecond pulsar PSR J0437-4715 at the 3σ of the confidence level. In Ref. [112], from the analysis of thermal emissions during thermonuclear X-ray bursts from neutron stars in low-mass X-ray binaries has been deduced at a range of 11.5 < R/km < 13 for neutron stars with mass between (1.3-1.8) M . In Ref. [111], the authors reported a study of spectroscopic radius measurements of twelve neutron stars obtained during thermonuclear bursts, or in quiescence. It was concluded that a 1.5 M neutron star has a radius in the range 10.1 < R/km < 11.1. Very recently, new constraints on the mass radius relation have been provided by the Neutron Star Interior Composition Explorer (NICER) experiment. A constraint on the radius measurement on PSR J0740+6620 has been recently proposed in Ref. [115], where it was put as a limit of R < 13.7 −1.5 +2.6 km, and in Ref. [116], where it was concluded that R < 12.39 −0.98 +1.3 km. The above results derived from two independent analyses of the X-ray data taken by NICER and the X-ray Multi-Mirror (XMM-Newton) observatory. Previously, NICER reported estimates of the mass and radius of the isolated 205.53 Hz millisecond pulsar PSR J0030+0451 obtained using a Bayesian inference approach to analyze its energy-dependent thermal X-ray waveform [117,118]. The best estimates of mass and radius turned out to be: M = 1.44 +0. 15 −0.14 M R = 13.02 +1.24 −1.06 km at 68% of the confidence level. Future more refined measurements of this kind may offer the opportunity to further constrain the EOS of high density matter. Mass-radius relations for several EOSs models discussed in the text. The observed pulsars PSR J0348+0432 [103] and PSR J0740+6620 [104], which are represented by yellow and purple bands, respectively, indicate the uncertainty on the measurement. See text for details.
The black continuous line in Figure 1 is the mass-radius relation derived in Ref. [101] in a relativistic Hartree-Fock (RHF) calculation using the QMC model in the SU(3)-flavour symmetry. Coupling constants and form factors were taken from the Nijmegen-Soft-Core 2008 interaction. The authors showed that the only hyperon that appears in β-stable neutron star matter, according to the chosen model setting, is the Ξ − one. This is shown in the right panel of Figure 2, where the particle composition for this model is reported.
The maximum mass predicted is consistent with a 2M neutron star. A similar conclusion was obtained also by Huber [119], considering a Walecka model in HF approximation and using hyperon-meson couplings compatible with hypernuclei calculations performed in Hartree approximations. The authors of Ref. [101], pointed out that the same result cannot be replicated if the SU(6)-quark symmetry is considered. The mechanism that allows an EOS hard enough to support a 2M neutron star is due to several competing effects. The Fock-term, in general, tends to soften the EOS, except in the channel that takes into account the tensor coupling. The last one acts in a such way as to shift the hyperon onset to a larger density. In addition, medium modification of the baryon structure provided by the QMC model produces a repulsive effect which gets the full stiffer. The red dotted line (FSU2R) in Figure 1 has been obtained in Ref. [120] using the IUFSU RMF parametrisation including hyperons. The model does not take into account the effect of the σ * -meson while strong repulsion is provided by the effective YY-interaction given by the φ-meson couplings. Note that one of the first systematic studies on the effect of the variation of hyperon single particle potential depths in β-stable neutron star matter was carried out in Ref. [121], employing a standard σωρ Walecka model with the inclusion of the φ-meson repulsion in the hyperonic sector. In this work, the authors concluded that the main effect on maximum neutron star mass is not provided by the inclusion of YY repulsion caused by the φ-meson. The blue long-dashed and green-dashed lines in Figure 1, represent the gravitational mass-radius relation found in Ref. [122] using the DD2 [40] and DDME2 [39] parametrisations with hyperonic degrees of freedom. At nucleonic level DD2 and DDME2, the EOSs differ in the fit performed to fix the parameters of the model. In the DD2 case, parameters were fixed in order to reproduce the properties of selected double magic atomic nuclei and additionally requiring that the neutron skin thickness in 208 Pb be well reproduced. In the case of the DDME2 model, the properties of twelve spherical nuclei were used to adjust the interaction Lagrangian and, in addition, employed some data on the excitation energies of isoscalar monopole and isovector dipole giant resonances in spherical nuclei. Concerning the hyperon single particle potentials in symmetric nuclear matter, they were fixed as follows: U Ξ = −17.5 MeV for the DD2 parametrisation and U Ξ = −18.78 MeV for the DDME2 parametrisation; while for both models, it was assumed that U Σ = +30 MeV and U Λ = −32 MeV. The Ξ-single particle potential depth was fixed as suggested by the observation of Ξ-hypernuclei [7,44], while the U Λ value was obtained, fixing the single-Λ hypernuclei separation energies in the s and d shells as detailed in Ref. [51]. These two EOSs have the nice feature of being consistent with the existing hypernuclear data and to reproduce the constraints of a two solar mass neutron star. The corresponding particle composition for the two models is shown in Figure 3. The nucleonic part of the EOS does not provide significant variations to the concentration of hyperons. Note that, in this case, there is a richer aboundance of hyperons compared to the RHF calculation by Katayama and collaborators. This scenario clearly reflects the strong uncertainties on the hypernuclear interactions, especially at a large baryonic density. The maroon and purple dot-dashed lines in Figure 1 represent the gravitational mass-radius relation found in Ref. [57] using the SU(3) chiral model without the inclusion of ∆-isobars (blue line) and with the inclusion of ∆-isobars (green line). The inclusion of the ∆-degrees of freedom makes the EOS softer, which results in a maximum mass below the 2M limit. The particle composition for the model without ∆'s is shown in the left panel of Figure 2. Among the considered EOSs, those based on the Chiral SU(3) model give rise to radii larger than those based on the density-dependent models DD2Y and DDMEY, the RHF by Katayama, and the FSU2H model. All models considered seem more consistent with the radii measurements reported by NICER in Ref. [116]. The formation of hyperons leads to an EOS softer than that with a a pure nucleonic content. Consequently, hyperonic neutron stars have smaller radii compared to the nucleonic ones. This discussion shows that more precise mass-radius measurements should be provided in order to get more stringent constraints on the EOS of matter at large densities and, consequently, on the neutron star structure. Finally, we report, for comparison, the mass-radius relation for the old GM3Y model proposed by Glendenning [41] in order to reconcile the canonical value of a 1.4M neutron star with single-Λ hypernuclei separation energy. This last EOS is clearly inconsistent with recent astrophysical observations; however, it is worth analysing it in order to understand the missing physical mechanism responsible for the stiffening of the EOS at a large baryonic density. Note that, in all other EOSs considered, one of the key ingredients to get repulsion at a large density is the YY repulsion produced by the exchange of the φ-vector meson. This mechanism does not alter the onset of hyperons in β-stable hadronic matter because the φmeson does not couple with nucleons. The φ-meson effect was not taken into account in the GM3Y model, which indeed resulted in a soft EOS. I want to stress that some authors also suggest other mechanisms to provide stiffening to the hyperonic matter EOS. In Ref. [123], they proposed a study of hypernuclear matter within a relativistic density functional theory with density-dependent couplings. The hyperon-scalar meson couplings were obtained by allowing for mixing and breaking of SU(6) symmetry keeping the nucleonic coupling constants fixed. The authors showed that in a restricted parameter space of the coupling constants, massive neutron stars with mass ∼2.25M could be obtained. Maslov [124] proposed an RMF model where both baryon masses and coupling constants carry a σ-meson dependence. The couplings of hyperons with vector-mesons were fixed by SU(6)-symmetry. The parameters of the model allow for reproduced constraints from both heavy-ion-collisions and astrophysics. The authors pointed out that the hyperon puzzle can be solved in this framework if a dependence of the φ-meson mass on the σ field is introduced. This last assumption gives rise to an increased stiffening of the EOS. Finally, Bednarek [125] proposed inclusion in the Lagrangian density quartic terms involving the φ-, ρ-, and ω-meson fields, plus mixed terms of the form φ 2 ω 2 , φ 2 ρ 2 , and ρ 2 ω 2 , getting an EOS able to support a 2M neutron star. A similar approach was also used by Taurines [126] and Gomes [127], who reached the same conclusion.
The RMF formalism recently discussed for cold hadronic matter has been extended at finite temperatures in several works [128][129][130][131][132]. Finite temperature EOSs are fundamental in order to study the properties of proto-neutron stars and investigate the evolution of dynamical systems like BNSM and SNe. I briefly discuss these issues in the last part of this review.
I conclude this discussion on the EOSs based on RMF models by making some general comments on the EOS of hyperonic matter based on this approach. All the EOSs proposed are able to explain the existence of hyperonic massive neutron stars, and include some specific mechanisms to provide repulsion at large baryonic density. RMF models are sufficiently flexible to take care of these essential ingredients in a rather easy way. In the next section, I will take up microscopic approaches, which are much more complex to implement from a numerical point of view. However, I want to stress that similar repulsive mechanisms required by RMF models to get the EOS of hyperonic matter stiff, are also present in microscopic frameworks and have their origin in the many-body forces between nucleons and hyperons. A tedious technical complication concerning microscopic approaches is that many-body forces are quite difficult to implement. Nonetheless, they represent an essential ingredient that naturally emerges due to the very same nature of baryons. The latter are not point-like particles, and their internal structure gets modified when baryons are near each other. The modification of the internal structure of baryons, from a classical point of view, produces some polarisation effects, which are similar to those that originate in a three-body system earth-sun-satellite. Due to the action of the gravitational force and the finite size of the earth and sun, matter distribution inside the earth loses its spherical symmetry, giving rise to tidal forces. A satellite orbiting around the earth is consequently affected by these tidal forces; the latter would not be present if the size of the earth and sun could be considered point-like. This is an example of a three-body force realised in classical physics.

Microscopic Approaches
Microscopic approaches are based on a non-relativistic framework and the interaction between particles is described by a potential. One of the main difference between phenomenological and microscopic approaches is that in the second ones, the interaction is fixed forever from the very beginning of each calculation. Among the microscopic manybody methods, we can distinguish between diagrammatic methods like: Brueckner-Hartree-Fock [133], Many-body perturbation theory [134], Self-Consistent-Green-Function [135,136], and exact approaches like: Diffusion Monte Carlo [137], Green function Monte Carlo (GFMC) [138], and Auxiliary-Field-Monte Carlo (AFDMC) [139]. According to the diagrammatic methods, the thermodynamical properties of the considered system are calculated perturbatively as the sum of specific contributions that can be represented as diagrams and can be evaluated according to well defined rules. Exact approaches are, in principle, able to solve the many-body problem in an exact way.
Currently, to the best of my knowledge, only the BHF and AFDMC have been extended to include hyperonic degrees of freedom. I now briefly review the main aspects of both methods. Additional details can be found in the quoted literature.
The main idea behind the BHF approach is that the repulsive hard core, present in the baryon-baryon interaction, does not allow a conventional perturbation theory to calculate the energy per particle of the system. The BHF approach designs a method to get the matrix elements of the interaction of such a value that the perturbation theory converges quite quickly. The bare baryon-baryon interaction is replaced by the so called G-matrix, which describe the interaction between two baryons in the presence of a surrounding medium. The G-matrices are obtained by solving the coupled-channel Bethe-Goldstone equation, written schematically as: where the first (last) two subindices indicate the initial (final) two-baryon states compatible with a given value S of the strangeness, namely nucleon-nucleon (NN) for S = 0 and nucleon-hyperon (NY) for S = −1; V B 1 B 2 ,B 3 B 4 is the bare baryon-baryon interaction; Q B i B j is the Pauli operator, that prevents the intermediate baryons B i and B j from being scattered to states below their respective Fermi momenta; ω is the starting energy, which corresponds to the sum of the nonrelativistic single-particle energies of the interacting baryons. The single-particle energy of a baryon B i is given by: where M B i denotes the rest mass of the baryon, and the single-particle potential U B i represents the average field "felt" by the baryon owing to its interaction with the other baryons of the medium. In the BHF approximation, U B i is calculated through the "on-shell energy" G-matrix, and is given by: where n B j (| k|) is the occupation number of the species B j , and the index A indicates that the matrix elements are properly antisymmetrised when baryons B i and B j belong to the same isomultiplet. I note that in order to determine the single-particle potentials when solving the Bethe-Goldstone equation, two possibilities are usually adopted. The first one is the so-called continuous prescription, where the single particle potential is assumed to be continuous above the Fermi momentum. The second possibility is to adopt the so-called gap choice: in this case the single particle potential is set to zero above the Fermi momentum. It has been shown by the authors of Ref. [140] that the contribution to the energy per particle from the three-hole line diagrams, which include correlations between three particle, is minimised employing the continuous prescription. The inclusion of threebaryon forces in the framework of the BHF approach follows a strategy implemented in the case of pure nuclear matter calculations where a density-dependent NN interaction is obtained from the original NNN force by averaging over the coordinates of one of the nucleons [141,142]. In the case of a generic BBB interaction, such an average is performed with respect to one of the baryons. The exact treatment of three-baryon forces in the BHF approach would require the solution of several coupled Bethe-Faddev equations in the medium. This task is currently too difficult to deal with. The Bethe-Goldstone equation is usually solved numerically using the inversion matrix technique [143] after a partial waves expansion. I note, however, that Ref. [144] proposed a method to solve Equation (10) without partial wave expansion.
I now turn to a discussion on the DMC method and one of its extensions, the AFDMC, that is suitable for applications in nuclear systems. The general idea of the DMC method is to solve, in a stochastic way, the time-dependent many-body Schrödinger equation. A formal solution of this equation in the imaginary time (it = τ) is given by: where G(R, R , τ) is the Green's function of the operator H + ∂ ∂τ that can be expanded as: {φ n } is a complete set of the eigenvectors of H, while the global coordinate R specifies a walker, namely the 3N spatial coordinates of the particles. Now, if we consider a generic trial wave function Ψ T as linear combination of {φ n }, the imaginary time evolution is given by: where E 0 is the exact ground state. If Ψ(R, τ) is not orthogonal to Ψ(R, 0), for large τ, it converges to the ground state Ψ 0 of H. The evolution can then be done by solving the integral of Equation (13) in the limit of infinite imaginary time. For a non-interacting system, it can be shown that the Green's function is a simple Gaussian with variance proportional to τ. Such a Gaussian represents a Brownian diffusion of a set of particles with a dynamic governed by random collisions. In the case of a generic interaction, several complications arise because, in general, the kinetic and interaction operators do not commute. However, by using the Trotter formula of order (∆τ) 3 , valid in the limit ∆τ → 0, the Green's function can be broken up into the product of a diffusion term containing the interaction, and a free Green function. The resulting integral can be finally worked out in a Monte Carlo way by propagating the particle coordinates by sampling a path according to the diffusion term in the integral. Several aspects about the sampling of the states are important for the correct implementation of DMC. The interested reader can find in the quoted literature, further details about these issues. The AFDMC is an extension of the diffusion DMC method to deal with interactions with spin-isospin dependence. The main issue against the direct implementation of the DMC to nuclear systems is that the nuclear Hamiltonian contains terms that are quadratic in the spin-isospin operator. This prevents solving the total wave function as a product of single-particle spin-isospin states, which is essential for the implementation of the method. The idea of AFDMC is to rewrite the Green's function in such a way as to transform the quadratic dependence on spin and isospin operators into a linear dependence by means of the Hubbard-Stratonovich transformation. The direct implementation of the GFMC is not feasible for an infinite system due to the unfavourable scaling behaviour of the computational time. In contrast, in the AFDMC approach, the propagator is rewritten by applying the Hubbard-Stratonovich transformation using auxiliary fields, which changes the scaling behaviour favourably at the cost of additional integrations over auxiliary fields.
Considering the trial wave function, a simple choice for Ψ(R, 0) is given by a Slater Determinant of one-body space-spin orbitals multiplied by a purely central Jastrow correlation operator: In quantum Monte Carlo simulations for light nuclei, a more realistic correlation operator, that includes spin-dependent operators, is generally used. Usually, the operators employed are those present in the two-and three-body potential in the Hamiltonian. This requires, however, a complete sum over the spin-isospin degrees of freedom. To avoid an exponential growth with N in the calculation time, a sample of the spin-isospin degrees of freedom is usually performed, rather than a complete sum. In order to do this, the imaginary time propagator exp(−(H − E 0 )τ) is applied to the walkers and a sample on the final configuration is carried out. It is, therefore, necessary to provide an efficient method for sampling the states produced when they are operated on by the imaginary time propagator.
the walker remains the same in this limit. The fermion sign problem, which affects fermionic systems on the lattice, in the considered case, becomes a phase problem as the overlap of the walkers with the complex trial function. In order to deal with this, the path of the walkers is constrained to regions where the real part of the overlap with the trial function is positive [145]. For spin independent potentials, this reduces to the so called fixed-node approximation [146]. For more details and a recent review on the applications of AFDMC to nuclear systems, see Refs. [147,148]. The AFDMC has been extended to finite and infinite systems, which includes hyperons in Ref. [149]. Despite its complexity, the AFDMC method has been successfully used to describe finite nuclear systems [150,151], as well as infinite nuclear matter [152,153]. A benchmark study between the BHF and the AFDMC approaches for pure neutron matter has been reported in Ref. [154].

The NY and YY Interactions
The NY and YY interactions have been developed mainly by the Nijmegen and Jülich groups. The first potentials have been first developed by the Nijmegen group and are known as hard core potentials [155,156]. Later, in 1989, the two groups constructed new soft core potentials, that included both the NΛ and NΣ channels, using the available scattering data in the NY sector. The two potentials are based on the meson-exchange model. Specifically, the Nijmegen Soft Core 1989 (NSC89) [157] potential takes into account the one meson exchange process considering the following exchanged mesons: π, η, η , ρ, ω, φ, , and S * Regge trajectories. In addition, the J = 0 contributions from the tensor f , f , A 2 , and pomeron exchange are also considered. The meson-baryon coupling constants are calculated via SU(3) symmetry, using the coupling constants of the NN analysis as input. The interactions described by this potential are limited to the NΛ and NΣ channels. The potential is regularised using a Gaussian cutoff in order to assure a regular behaviour at short distances.
The potential provided by the Jülich group in 1989 (Ju89) [158] differs from the NSC89 by also considering the two-meson exchange processes. The mesons included in the potential are: π, ρ, K, K * . All the possible combinations of two meson exchange diagrams built using the previous mesons are considered. In addition, the excitations of ∆ and Σ * baryons in the NY-meson vertices of the various Feynman diagrams are also included. Coupling constants and cutoffs of this potential involving N and ∆ baryons, are taken from the Bonn NN potential [159], while the couplings for hyperons are related to the nucleonic ones using the SU(6) symmetry. The free parameters of the Ju89 potential are the values of cutoff at vertices involving strange particles. Such freedom is fixed using the available data scattering. Later, in 2005, a new version of this NY potential was proposed [160]. In this potential, the various contributions to the NY scattering amplitude are restricted to those of the baryonic octect.
The first evolution of the NSC89 potential is represented by the Nijmegen Soft Core 1997 (NSC97) [161], which also includes the YY interaction in channels with strangeness S = −2, −3, −4 of the baryonic octect. A new version of the NSC97 potential was the Extended Soft Core 2008 (ESC08) [162][163][164]. In this potential, the two-meson exchange process is implemented. The ESC models describe the NN, NY, and YY interactions in a unified way using broken flavour SU(3) symmetry. The NSC97 and ESC08 potentials are characterised by different choices of the magnetic vector F/(F + D) ratio, α m v , which produces different scattering lengths in the NΛ and NΣ channels. This freedom allows a description in a similar way to the NY scattering data. The latest versions of the baryonbaryon potentials by the Nijmegen group are the extended-soft-core baryon-baryon models (ESC16) [165,166]. The main new features of these potentials are the inclusion of the axialvector meson potentials, a zero in the scalar-and axial-vector meson form factors, and a treatment of the strange scalar κ-meson within the scheme of the Gell-Mann-Okubo mass relations resulting from considering the κ-meson as a broad one like the ρ and mesons. In contrast to ESC08, in the ESC16, potential is not considered a medium strong flavorsymmetry breaking of the coupling constants. In addition, the multiple-gluon exchanges are added contributions, due to the odd number of gluon exchanges. Finally, a novel contribution is represented by the inclusion of structural effects, due to the quark core of the baryons.
A very promising approach to describe the NY and YY interaction is provided by Chiral Effective Field Theory (ChEFT). This theory derives a generic BB interaction and, more generally, many-body baryon interactions according to a well defined scheme. ChEFT for low-energy quantum chromodynamics (QCD) opened a new and systematic way to describe the interaction between baryons. In 1960, following the ideas developed by Schwinger [167], Gell-Mann and Levy [168] proposed a linear realisation of chiral symmetry. This model was elaborated to solve an issue related to the the pion-nucleon scattering length, which came out two orders of magnitude too large when the pseudo-scalar pionnucleon interaction was employed. In the linear sigma model, the problem was fixed by including the contributions of a fictitious σ-meson. Such a solution was not, however, very satisfying, due to the reliance on cancellations of large terms and the fictitious character of the σ-meson. Later, inspired by a suggestion of Schwinger [169], Weinberg worked out a general theory of non-linear realisations of chiral SU(2) × SU(2) [170,171], which was soon generalised to arbitrary groups in the papers of Callan, Coleman, Wess, and Zumino [172,173]. The considerable advantage of using ChEFT lies in the fact that twobody, as well as many-body, baryonic potentials can be calculated order by order, according to a well defined scheme based on a low-energy effective QCD Lagrangian, which retains the main symmetries of QCD and the approximate chiral symmetry. This ChEFT is based on power counting in the ratio q/Λ χ , where q denotes a low-energy scale identified with the magnitude of the three-momenta of the external baryon; whereas, Λ χ ∼ 1 GeV denotes the chiral symmetry breaking scale. In the pure nucleonic sector, single-pion exchange and multi-pion exchange give the long-and intermediate-range part of the nuclear interaction; whereas, the short-range components are included via contact terms. Within this approach, the details of the QCD dynamics are contained in the parameters, the so called low-energy constants (LECs), which are fixed by low-energy experimental data. Note that the relativistic treatment of baryons in ChEFT leads to some problems because the time-derivative of a relativistic baryon field generates a factor E ∼ M (M being the baryon mass), which is not small compared to the chiral-symmetry breaking scale (indeed, M/Λ χ ∼ 1). In addition, the baryon mass does not vanish in the chiral limit. A solution to this problem has been proposed in Refs. [174,175], where the so called Heavy-Baryon chiral EFT has been developed. The basic idea is to treat the baryons as heavy static sources in such a way that the momentum transfer between baryons by pion exchange is small compared to the baryon mass. The expansion is, thus, carried out in terms of these small momenta over the baryon mass. For further details see the quoted references. For a review on these issues see, for instance, Ref. [176]. It is important to remark that although there is a vast number of NN and NNN interactions derived so far in ChPT, NY, and YY interactions, they have been constructed only by the Jülich-Bonn-Munich group within this framework [177][178][179][180]; it was developed first by the NY interaction at the leading order (LO) [177], and next-to-leading order (NLO) [178]. Then, the YY interaction was constructed at NLO [180]. Recently, the NY potential at NLO, reported in Ref. [178], has been revised by using a different strategy to fit the LECs appearing at this order of ChEFT [179]. Specifically, in the potential reported in Ref. [179], some of the S-wave LECs were inferred from the NN sector via the underlying SU(3) flavour symmetry; in this way, only a reduced number of LECs were needed to be determined from the empirical information in the NY sector.
A totally different (and in a way, more fundamental) method to derive the baryonbaryon interaction has been recently provided by lattice QCD calculations performed by the HALQCD and NPLQCD collaborations. The two groups adopt different approaches to the problem. The HALQCD collaboration [181] tries to extract the various baryon-baryon interactions from the Nambu-Bethe-Salpeter wave function calculated on the lattice in terms of the derivative expansion, which is shown to reproduce the scattering phase shifts correctly below the inelastic threshold. Using this definition, it is possible to formulate a method to extract the potential in lattice QCD. Results obtained in the three-flavour lattice [182][183][184] for the NN potential show that they reproduce the qualitative features of the phenomenological NN potentials, namely, attractive wells at long and medium distances, central repulsive cores at short distances, and strong tensor forces with a negative sign. At long distances, it was observed that the ranges of the tail structures in the central and tensor forces become longer at lighter quark masses. At short distances, the repulsive cores in the central forces are found to be enhanced at lighter quark masses. This could be explained by the short-range repulsion due to the one-gluon-exchange in the quark model, whose strength is proportional to the inverse of the constituent quark mass. Such studies extended to the case of hyperon forces with the same lattice setup reveal that the nature of the repulsive core is well-described by the quark Pauli blocking effect together with the one-gluon-exchange effect [182,183,[185][186][187][188].
In some recent works [189,190], the HALQCD collaboration reported the study of ΛΛ, NΞ, and NΩ systems near the physical mass of the pion (m π ∼ 146 MeV) and the kaon (m K ∼ 525 MeV). For further details, and a recent review about the results and methods of the HALQCD collaboration, see Ref. [191].
The strategy of the NPLQCD collaboration is to numerically evaluate path integrals representing Euclidean correlation functions using Monte Carlo sampling methods. In recent work [192], the NPLQCD collaboration performed a study of BB scattering, including baryons of the octect, up to a channel with strangeness S = −4. In Ref. [192], they adopted values of pion and kaon masses of 450 and 596 MeV, respectively. At the unphysical pion mass, the authors found bound BB systems, but in the spin-triplet ΣN and ΞΞ channels. In addition, using results from some previous studies at larger pion masses m π ∼ 806 MeV [193], the authors performed a naive extrapolation of the results of the binding energies of the two-baryon systems at the physical point. In this way, a comparison with experimental values and phenomenological predictions was possible. The results for ground-state energies of two-nucleon systems were found to be compatible with the experimental values. Furthermore, stronger evidence for the existence of bound states in the ΞΞ ( 1 S 0 ) and NΞ ( 3 S 1 ) channels was observed. Finally, the determination of scattering parameters at low energies of two-baryon systems allow constraints on the LO and NLO interactions of a pionless EFT, for both the SU(3) flavor-symmetric and symmetry-breaking interactions.
In conclusion, lattice approaches are very promising in the attempt to derive the interaction between baryons from a more fundamental point of view. In addition, they provide a powerful tool to calculate nuclear and hypernuclear observables using a unique approach comparing those usually employed for studies in nuclear and hypernuclear physics.

Results from Microscopic Calculations
I now discuss the main works where the microscopic BHF and AFDMC approaches have been used to describe infinite hyperonic matter. A study of the onset of the Λ and Σ hyperons in β-stable neutron star matter was first proposed in Ref. [194] in the framework of a BHF approach. The NN force employed was the AV14 [195] potential, supplemented by a density-dependent-interaction derived by the UIX three-nucleon force [196]; concerning the NY interaction, where the NSC89 potential was adopted. The conclusion of this study was that hyperons are expected to appear at 2-3 times normal saturation density in cold β-stable neutron star matter. In Ref. [197], they carried out a study of hypernuclear matter based on the NN, NY, and YY interactions provided by the NSC97 interaction. The first EOS calculation of β-stable hyperonic matter based on the BHF approach was performed in Refs. [198,199]. These works concluded that the appearance of hyperons strongly softens the EOS of neutron star matter giving rise to maximum masses of the order of 1.3M . Some subsequent studies [200,201] used some different updated versions of the NY and YY interactions without changing the previous conclusion. In Ref. [202], Djapo and collaborators reported several BHF calculations based on a large variety of NY interactions provided by the Nijmegen (NSC89, NSC97) and Jülich (Ju89, Ju04) groups supplemented by a phenomenological three-nucleon force, and using the V low−k technique to reduce the high-momentum component of the baryon-baryon interaction. Furthermore, in this case very low maximum neutron star masses were found. This scenario was well summarized in Ref. [203], where several hyperonic EOSs based on different combinations of NY, NN, and NNN interactions were compared. The authors of Ref. [203], concluded that some new mechanism should be considered in order to solve the neutron star hyperon puzzle.
This quite intricate scenario, has been improved in the last years when, following an old idea by Takatsuka and collaborators [204,205], some authors tried to incorporate the effects of many-body forces between nucleons and hyperons. Note that such baryonic many-body forces can play an important role is a quite natural one; keeping in mind the studies carried out in pure nucleonic systems (both finite and infinite). It is well known that it is not possible to correctly reproduce binding energies of finite nuclei or good nuclear matter at saturation properties unless to introduce a three-nucleon force [206,207]. From a microscopic point of view, such many-body forces originate from the fact that nucleons have a quite complicated quark structure, which strongly changes when nucleons are near each other. Recently, a consistent theory of three-baryon forces in SU(3)-flavour symmetry, involving the octet baryons has been successfully developed in the framework of ChEFT [208]. Note that the ChEFT of many-body nucleonic forces was discussed earlier in Ref. [209]. Later works derived explicit expressions of the many-body nuclear force at various orders of ChEFT [210][211][212].
A preliminary study on the impact of hyperonic three-body forces on neutron star matter was proposed in Ref. [213], where they developed a model based on a microscopic BHF approach of hyperonic matter employing the Argonne V18 [214] NN force and the Nijmegen NY soft-core NSC89 force, supplemented by some phenomenological densitydependent contact terms, to set a numerical lower and upper limit to the effect of hyperonic three-body forces on the maximum mass of neutron stars. Assuming that the strength of these forces is either smaller or as large as the pure nucleonic ones, the results of this work showed that, although hyperonic three-body-forces can reconcile the maximum mass predicted by microscopic approaches with the canonical value of 1.4-1.5M , they are unable to provide the repulsion needed to make the predicted maximum masses compatible with the recent observations of massive neutron stars.
In Ref. [215], it was reported that a Monte Carlo calculation of neutron matter with nonvanishing concentrations of Λ-hyperons including neutron-neutron (nn), neutronneutron-neutron (nnn), neutron-Λ (nΛ), and neutron-neutron-Λ (nnΛ) interactions. Specifically, in order to describe the nn and nnn interactions, the authors employed the AV6' and UIX potentials. Concerning the nΛ interaction, the authors used a phenomenological potential fitted to the available NΛ scattering data, while they considered a purely central repulsive nnΛ force of the same form of the one present in the UIX interaction. The authors of Ref. [215] proposed two parametrisations (hereafter I and II) for the nnΛ force. Such parametrisations were developed in Ref. [149], with the aim of reproducing the separation energies of several single-Λ hypernuclei. The resulting EOSs obtained using the settings described above are reported in the right panel of Figure 4. The green line refers to the pure neutron matter system, while the red and blue dashed curves correspond to the EOSs obtained, employing the nΛ potential alone, and the full nΛ+nnΛ (I) interaction, respectively. The onset of the Λ hyperon is shown in the left panel of the same figure for both calculations. The blue lines in the left panel of Figure 4 refer to the case where the full nΛ+nnΛ (I) interaction was adopted. Using parametrisation II (black dots in the right panel of Figure 4), the authors found that the effect of the hyperonic three-body force is so large that Λ-hyperons do not appear in neutron matter up to a density corresponding to the central one, of a two-solar-mass neutron star. According to such calculations, this rather simplified version of the hyperon puzzle is solved by the fact that hyperons are not formed in the core of neutron stars due to the very strong repulsive effect of three-baryon forces. According to parametrisation I (blue line in Figure 5), hyperons appear around 2n 0 , giving rise to a very soft EOS and to a quite low maximum mass of the order ∼1.34M , not consistent with observations. A similar result was also obtained using of the sole nΛ potential, where the resulting maximum neutron star mass turned out to be below 1M (red line in Figure 5).   In Ref. [215] was reported a Monte Carlo calculation of neutron matter with n nishing concentrations of Λ-hyperons including neutron-neutron (nn), neutron-ne neutron (nnn), neutron-Λ (nΛ) and neutron-neutron-Λ (nnΛ) interactions. Speci in order to describe the nn and the nnn interactions the authors employed the AV the UIX potentials. Concerning the nΛ interaction, the authors used a phenomenol potential fitted to the available NΛ scattering data, while they considered a purely c repulsive nnΛ force of the same form of the one present in the UIX interaction. The au of Ref. [215] proposed two parametrizations (hereafter I and II) for the nnΛ force. parametrizations were developed in Ref. [149] with the aim of reproducing the sepa energies of several single-Λ hypernuclei. The resulting EOSs obtained using the s described above are reported in the right panel of Fig. 4. The green line refers to th neutron matter system while the red and blue dashed curves correspond to the Ref. [215]). The observed pulsars PSR J0348+0432 [103] and PSR J0740+6620 citecro19, which are represented by yellow and purple bands, respectively, indicate the uncertainty on the measurement. See text for details. Some years later in Ref. [216], they performed a new calculation in the framework of non-relativistic BHF approaches using realistic NN and NNN interactions derived in chiral effective field theory (χEFT), supplemented by NΛ and NNΛ interactions. Specifically, the two-body NN interaction was used for the local chiral potential presented in Ref. [217] at the next-to-next-to-next-to-leading order (N3LO) of ChEFT, which includes the ∆(1232) isobar in the intermediate states of the NN scattering. Regarding the NNN force, the authors made use of the potential derived in Ref. [210] at the next-to-next-to-leading-order (N2LO) in the local version reported in Ref. [218,219]. Note that this NNN force takes into account the possibility of the ∆-excitation at the Nπ vertex. The low energy constants of the NNN interaction were fixed as discussed in Ref. [219], where it was shown that a good description of nuclear matter can be achieved using that setting. These interactions have been employed in Ref. [220] to calculate the β-stable EOS of nuclear matter and the structure of neutron stars. It was found that a neutron star maximum mass of 2.07M , in agreement with the largest measured neutron star masses. The resulting EOS has also been recently used in Ref. [221] to simulate the merging of two equal mass neutron stars. This EOS has been extended to finite temperatures in Ref. [222]; the finite temperature version of this EOS has been employed in Ref. [223] for the study of several asymmetric mass simulations of BNSs.
Concerning the NY interaction, the authors of Ref. [216], employed two versions (a and e) of the NΛ NSC97 meson-exchange interaction. Finally, for the NNΛ force, the authors adopted the NNΛ three-baryon force recently derived by the Jülich-Bonn-Munich group in the framework of χEFT [208]. In Ref. [216], they provided two parametrizations of this NNΛ interaction: NNΛ 1 and NNΛ 2 ; set to reproduce the depth of the Λ single particle potential in symmetric nuclear matter at saturation density, varying the value of the only low energy constant present at this order in the so called decuplet saturation approximation. According to Ref. [42], the authors adopted: U Λ = −28, −30 MeV.
The composition of β-stable neutron star matter according to this model is shown in the left panel of Figure 6 for the models NSC97a and NSC97a+ NNΛ 1 . Qualitatively similar results were obtained by employing the other NY and NNY interactions. The continuous lines show the results when only NΛ, in addition to NN and NNN forces, are taken into account, whereas the dashed ones also include the contribution of the NNΛ force. The effect of the latter is twofold. First, it shifts the onset of the Λ-hyperon to slightly larger baryonic densities. The second effect, may be the most important one, is that the NNΛ force strongly reduces the abundance of Λ particles at large baryonic densities leading to a stiffening of the EOS comparing to the case in which the NNΛ force is not included. This can be appreciated in the right panel of Figure 6, where the total pressure P is reported as a function of the total energy density ε. Consequently, the mass of the neutron star, and in particular its maximum value, increases. This is shown in Figure 7, where the mass-radius relation is reported for the models NSC97a and NSC97e, with and without the inclusion of the NNΛ force. The black line corresponds to the case of pure nucleonic matter included as reference. It is remarkable that the maximum masses obtained, including the NNΛ force are compatible with the largest measured masses of ∼2M [103,104,106]. An additional interesting feature found in Ref. [216] is that the overbiding of the Λ separation energies in heavy hypernuclei found using the sole NΛ interaction in microscopic calculations [224,225], is strongly reduced, adding the NNΛ 1 or the NNΛ 2 interactions. These results are consistent with the findings of Refs. [149,215].  [216]. See text for details.  [216]. The observed pulsars PSR J0348+0432 [103] and PSR J0740+6620 [104], which are represented by yellow and purple bands, respectively, indicate the uncertainty on the measurement. See text for details.
The only calculations employing chiral NY interactions in infinite hyperonic matter have been reported in Refs. [226,227], using the BHF approach, and considering the case of nonvanishing Λ particle concentrations in symmetric nuclear matter and pure neutron matter. These works agree with the conclusion of Refs. [215,216] that the NNΛ interaction becomes very repulsive at large baryonic densities, and that this effect provides an important stiffening on the EOS of hyperonic matter at least in the restricted scenario considered. The single particle potentials of Λ and Σ −,0,+ hyperons have also been evaluated in Refs. [228,229] using the chiral NLO NY interaction derived in Ref. [178], and including the effect of a NNY interaction fixed using the decuplet saturation. In Ref. [229], they pointed out the importance of including NNΛ-NNΣ transition channels in order to get a good description of Λ-single particle potential in symmetric nuclear matter. Comparing to NSC97 and Ju04 NY potentials, the interactions derived at NLO in ChEFT have a unique feature to predict a repulsive value of the potential depth in symmetric nuclear matter for the Σ − hyperon. This result is more consistent with the present experimental scenario. Note, however, that the most recent versions of the NY potential of the Nijmegen group, namely the ESC08 and ESC16, show a similar feature.
In the framework of microscopic approaches, it should be mentioned that the results of Ref. [230], again based on the BHF approach, and employing the NSC08 NN and NY interactions supplemented by a phenomenological gaussian repulsive three-baryon force arising from a model of multi-pomeron-exchange (MPP). Such force was included universally in all BB channels keeping the same form and strength. The strengths of MPPs were adjusted in order to reproduce the experimental data in 17 O + 17 O scattering [230]. At a microscopic level, this mechanism mimics a multi-gluon interaction. The authors of Ref. [230] found that the inclusion of this repulsion provides the required stiffness to the EOS necessary to get a maximum neutron star mass compatible with two solar masses. The same conclusion was also obtained in Ref. [231], where it has been reported that a variational calculation of the hyperonic matter EOS employing the same phenomenological three-body force adopted in Ref. [230]. Concerning the nucleonic sector, the authors used the AV18+UIX interaction. Note that in the calculations reported in Refs. [230,231], besides the Λ-hyperon, the Σ − hyperon was also included.
According to these calculations, hyperons appear at a density around 2n 0 when the universal baryonic three-body force is taken into account; however, the total hyperon fraction remains low enough to keep the EOS stiff. This result is consistent with the findings of Ref. [216]. According to Ref. [215], the hyperon puzzle can be solved in a microscopic framework only if hyperons never appear in neutron star matter. Other calculations, together with a deeper understanding of BB and hyperonic three-body-forces are clearly needed to better outline this scenario.
I will end this part of the review focused on the neutron star structure predicted by microscopic approaches, discussing the attempt to include hyperonic degrees of freedom in the framework of the Dirac-Brueckner-Hartree-Fock approach. The Dirac-Brueckner-Hartree-Fock approach [232,233] is, in a way, an extension of the non-relativistic Brueckner-Hartree-Fock one. In the relativistic framework, the self-energies and interactions of the various particles are calculated according to the solutions of the Dirac equation in the HF approximation. In the last ones, the bare nuclear masses are replaced by effective masses, which can be calculated in a self-consistent way, therefore acquiring a strong in-medium effect. One of the main advantages of this approach is that the effect of a certain class of three-body forces is effectively included at the two-body level. This particular class of three-body forces takes into account the nucleon-anti-nucleon pair formation in the medium. This mechanism is able to provide good saturation properties of nuclear matter and a stiff nucleonic EOS able to support a two solar mass neutron star [234] at the nucleonic level. The first attempt to include hyperonic degrees of freedom in the DBHF framework has been reported by Sammarruca in Ref. [235], considering the effect of fixed hyperonic concentrations on the energy per particle of the system, and employing the Ju04 NY interaction. To the best of my knowledge, the only calculation of β-stable hadronic matter, in this framework, has been carried out by the authors of Ref. [236] who, using an effective NY interaction, showed that the repulsive mechanisms that account for the stiffness of the EOS at the nucleonic level, act in a similar way in the hyperonic sector. This repulsive mechanism, therefore, seems able to produce a large mass neutron star in agreement with observations.

Possible Signatures of the Presence of Hyperons in Astrophysical Dynamical Systems
In this section, I discuss the possible signatures of the presence of hyperons in dynamical systems like BNSMs and SNe. The typical thermodynamical conditions predicted by numerical simulations in full general relativity for these systems are very different to those of cold neutron star matter. Ranges of temperature and density found in BNSMs are about: 0 < T/MeV < 100 and 10 < /(g/cm 3 ) < 3 × 10 15 [237], while for SNe: 0 < T/MeV < 50 and 10 < /(g/cm 3 ) < 4 × 10 14 [12]. Note that the precise range for the above quantities depend on the total mass and mass ratio in the case of BNSMs and on the progenitor star mass in the case of SNe. These physical conditions would suggest that due to thermal effects, a relevant fraction of hyperons can be present, even around the subnuclear densities compared to the zero temperature case. However, this consideration is not always true and indeed, if neutrinos are trapped in neutron star matter, the opposite scenario is realised. In neutrino-trapped matter, the electron chemical potential is diminished by the neutrino. In β-stable nuclear matter, this process leads to a softening of the EOS and, consequently, the onset of hyperons shift to a larger density. I note, however, that once hyperons are formed, the presence of a neutrino sea leads conversely to a stiffening of the EOS. Neutrino diffusion strongly depends on density, temperature, and composition of stellar matter. Such thermodynamical conditions can rapidly change in the dynamical evolution of astrophysical systems. In order to establish if neutrinos are trapped or not, the neutrino diffusion timescale should be confronted with the typical dynamical timescale of the considered system [237][238][239]. Accurate studies of neutrino propagation in dynamical systems are very demanding from a computational point of view and represent a very active and still progressing research field. The impact of hyperons in BNSMs was studied for the first time in Ref. [240], using a version of the TM1 EOS that includes the Λ hyperon. Therefore, it was pointed out that the formation of hyperons strongly affect the dynamics of the hypermassive neutron star formed after the merging process. In particular, the softening of the EOS, induced by the formation of hyperons, leads to a faster collapse of BH. A very important point stressed in Ref. [240] is that the presence of hyperons is imprinted in GWs; therefore, GW observations can, in principle, provide a chance to explore the composition of neutron star matter. In Ref. [241], a similar analysis, based on a different EOS, was carried out. The EOS used in this case was the BHB [132], an extension of the nucleonic DD2 EOS that takes into account the presence of Λ-hyperon. The authors show that physical effects introduced by the presence of hyperons change the qualitative dynamics of the remnant evolution, but they are not identifiable as a signature in the GW frequency, with the exception of possible black-hole formation effects. Concerning the EOS softening, this is instead, incoded in the GW luminosity and phase, and is detectable.
Considering CCSN simulations, the role of hyperons was analysed in Ref. [242] employing the TM1 EOS with Λ hyperons, in a 1D hydrodinamical simulation of spherical core collapse and explosion without neutrino transport. Considering several models of progenitor stars, the authors found that hyperons appear around the core bounce and gain explosion energy by (3)(4))% in each model. In Ref. [243], a similar analysis showed that hyperons appear off the center of the system and later prevail at the center when the central density becomes high enough. Compared with a pure nucleonic EOS, the authors found small differences in the luminosity and in the average energies, while the neutrino emission stopped, much earlier in the case of the hyperonic EOS. A better treatment of the neutrino transport was proposed in Ref. [244], where the inclusion of pion condensation was also taken into account. The results were in agreement with the findings of the previous works. Similar conclusions were also obtained in Ref. [245].
Another scenario in which hyperons may play an important role, is in the neutron star cooling after the star birth. Neutron stars are formed in Type-II, Ib, or Ic SNe [12]. As stated before, thermal effects and neutrino trapping shift the onset of hyperons to a larger baryonic density, leading to a softening of the EOS. A neutron star, just after birth, is very hot in its core and the neutrino mean free path is much smaller than the star radius. In these conditions, neutrinos are trapped and form a neutrino sea. The baryonic mass of the neutron star with a good approximation remains constant along the neutron star evolution and, therefore, it can be considered as an independent variable. As discussed in Refs. [11,[246][247][248], the presence of both hyperons and neutrinos, and the subsequent deleptonisation process caused by neutrino escape, leads to a scenario of metastable protoneutron stars. If only nucleons and neutrinos are present in the system, there is a window of baryonic masses that can be stabilised after deleponisation while, in the case of nucleons and hyperons embedded in neutrino-trapped matter, the deleptonisation produces a window of baryonic masses unstable against gravitational collapse leading to a faster creation of BHs [11]. These different evolution scenarios clearly affect the neutrino emission signal from the protoneutron stars. The presence of hyperons may also potentially influence the cooling of newborn neutron stars. Neutron stars cool according to processes like: n → p + l +ν l p + l → n + ν l (17) The first two processes are known as direct (or fast) Urca processes and take place only above a certain proton fraction around 11%. The remaining two processes are known as modified (or slow) Urca processes, and can occur according to arbitrary matter composition. In the above equations, N stands for a generic nucleon (neutron or proton). In addition to these processes, other mechanisms can contribute; the main ones are the bremsstrahlung: and the Cooper pair formation: Cooper pair formation is strongly sensitive to temperature, and can occur only for temperatures below a critical threshold. In addition, the process is strongly sensitive to the nucleon paring gaps. I note that several calculation of nucleon pairing gaps have been proposed in the literature, considering different many-body approaches and nuclear interactions [249][250][251][252][253]. However, if hyperons are present in neutron star matter, other direct processes like B → B + l +ν l , or B → N + l +ν l , can contribute. This additional fast cooling mechanism may give rise to temperature cooling curves not compatible with observations. It should be remarked, however, that these processes are also strongly influenced by hyperon paring gaps and, unfortunately, not few calculations have been proposed in the literature on this issue [254][255][256]. For further details on the cooling of neutron stars, and for a proper discussion on transport problems related to it, see Refs. [257,258].
A final aspect regarding the physics of neutron stars, in which hyperons are expected to be relevant, concerns the so-called r-mode stability. Rotating neutron stars are known to be stable against collapse to BH if the rotational frequency is smaller than a maximum frequency set by the so called Keplerian frequency [259]. At this frequency, the centrifugal force would overcome gravity at the star's equator and, as a consequence, the star starts to eject matter from the poles. However, some mechanical instabilities affect the structure of a neutron star; thereby, not allowing it to reach the Keplerian frequency. One of these instabilities, known as r-mode instability, has received a lot of attention due to the fact that hot, rapidly rotating systems (like young neutron stars) subject to this instability are GW-emitters [260]. The r-mode is a toroidal oscillation mode, whose restoring force is the Coriolis force. The detection of GWs produced by this mechanism may provide important informations about the internal structure of neutron stars. The emission of GWs acts in such a way as to amplify the r-mode, while viscous processes tend to damp it. In order to establish which process dominates, one should compare the typical timescale of GW emissions with the one set by viscosity. If the latter is considerably slower than the first, the system has enough time to produce GWs; these, due to conservation of angular momentum, will carry away the imprint of the rotational frequency of the source. For hot systems, the main contribution to viscosity is given by the bulk viscosity, which is generated by gradients of pressure and density in stellar matter. These perturbations tend to lead matter out of β-equilibrium. As a consequence, a certain amount of energy is dissipated by the system to restore the weak equilibrium. In nucleonic matter, the main channels contributing to the bulk viscosity are the direct and modified Urca processes. However, the presence of hyperons opens the way to several new channels, like weak non-leptonic processes N + N → B + N, direct and modified hyperonic Urca, or even strong processes like B + B → N + B. First, studies on the effect of hyperons on the bulk viscosity can be found in Refs. [261][262][263]. In Ref. [264], the presence of antikaon condensation in regard to the bulk viscosity was analysed. The authors of Ref. [264] concluded that the presence of kaon condensation leads to a suppression of the bulk viscosity. A study of bulk viscosity due to nonequilibrium weak processes in superfluid hyperonic matter was proposed in Ref. [265], while the effects of magnetic fields were investigated in Ref. [266]. Such effects are important only for huge values of magnetic fields (of the order of 10 17 Gauss), which are not expected to exist inside neutron stars. In conclusion, a general remark reported by all works is that the presence of hyperons induces a smaller r-mode instability due to the growth of the bulk viscosity. Now, a final remark about the transport coefficients of hyperonic matter. In this sector, few calculations have been completed [267]. These calculations are relevant for the very important application in astrophysical systems, like BNSM and SNe. A microscopic study of thermal conductivity and shear viscosity of hyperonic matter have been very recently proposed in Ref. [267], in the framework of BHF approaches using versions a and e of the NSC97 NY and YY interactions. The effect of hyperonic three-body forces was not included. The authors found that neutrons dominate the total thermal conductivity over the whole range of densities explored and that, due to the onset of the Σ − hyperon, and to the subsequent deleptonisation of the neutron star core, they dominate the shear viscosity at high density. In the case of the pure nucleonic matter, the lepton contribution turns out to always be dominant.

Conclusions
Here, I have reviewed the role played by hyperons in neutron star matter. I have analysed the main issues related to the formation of hyperons in neutron star matter, considering the results provided by both phenomenological and microscopic approaches that have been (so far) extended to include hyperonic degrees of freedom. All models agree that the softening of the EOS of neutron star matter, caused by the formation of hyperons in β-stable neutron star matter, require specific repulsive mechanisms to support the large M max . In RMF models, these mechanisms can be realised in various ways, like taking into account the effect of the φ-meson repulsion in the YY sector, or introducing non linear self-interaction terms in the Lagrangian density, concerning the φ-meson and/or mixed terms involving the other vector mesons. In microscopic approaches, to get a stiff hyperonic EOS, it seems that a crucial role is played by the many-body forces between nucleons and hyperons. In some works, it was shown that the inclusion of these many-body forces is essential to provide repulsion to the hyperonic matter EOS [149,216,230,231]. In addition, it seems that the description of medium and heavy Λ-hypernuclei strongly improve when NNΛ forces are added to the bare NΛ interaction [149,216]. Unfortunately, due to the very few available scattering data, strong uncertainties persist in the knowledge of the NY and YY interactions. Clearly, these uncertainties are reflected and definitely amplified when the same interactions are used in applications to hyperonic matter. Hopefully, this scenario will improve in the near future. I have also discussed the possible signatures of the presence of hyperons in dynamical systems like BNSMs and SNe. Multimessanger physics represents, in this line, a very powerful tool to investigate the properties of matter under extreme conditions. Future detections of gravitational waves and electromagnetic counterparts from compact systems will definitely improve our knowledge on the EOS and composition of neutron star matter and will contribute to provide a deeper understanding on the physics of hyperons. In addition, the third generation of GW detectors, like the Einstein telescope will probably give the chance to observe the post merger signal of BNSMs. Measurements of this kind would provide precious information about the composition of the core of neutron stars.
I have finally, very briefly, discussed the consequences of the formation of hyperons in the cooling of neutron stars and in the growth of gravitational instabilities in rapidly rotating hot neutron stars. An important point in the present review, is that I did not considered the possibility of a quark-phase transition in neutron star matter. This scenario has been investigated by several authors, and represents a possible scenario, especially in hot deleptonised matter. Such conditions are those typical of BNSMs, and it has been shown in some works that the presence of hyperons may favour the phase transition to quark matter [25,26,89,268,269]. I must note, however, that this very interesting scenario should still be systematically tested in BNSMs simulations.  Acknowledgments: It's a pleasure to acknowledge Ignazio Bombaci for several discussions on the role of hyperons in neutron stars.

Conflicts of Interest:
The author declares no conflict of interest.