maximally repulsive (while avoiding instabilities of the numerical procedure): $(H_1,H_2) = (-2.2,0)$ and $(-2.47,1.2)$ in units of the squared inverse pion decay constant, $f^{-2} =4.8\,$fm$^2$. These are the parameters of choice for extrapolations to higher densities, as we are seeking mechanisms and conditions for avoiding a softening of the equation-of-state that would be in conflict with neutron star observations.

The computations of $U_\Lambda$ based on self-consistent solutions of the Bethe-Goldstone equation (\[eq:BGE\]) together with eq.(\[eq:U\]) involve a momentum cutoff not only in the input hyperon-nucleon potentials but also in the momentum-dependent single-particle potentials $U(k)$ that appear in energy denominators and are integrated over intermediate momenta $k$. We apply a cutoff by a factor $\exp[-(k/\lambda_{spp})^6]$ and choose $\lambda_{spp}$ between 600 and 650 MeV, close to the cutoff in the NLO hyperon-nucleon potential, while not altering the low-energy behavior. The ranges of these cutoffs, together with the two possible choices for the pair of constants ($H_1,H_2$), are reflected in the uncertainty estimates shown in the figures.

The computations of $U_\Lambda$ for both symmetric nuclear matter and pure neutron matter can be carried out safely (i.e. with well converging and numerically stable results) up to densities $\rho\simeq \rho_c=3.5\,\rho_0$. At even higher densities the 