occurence of numerical instabilities becomes prohibitive. One might argue that, in any case, performing calculations at $\rho\gtrsim 2\rho_0$ is not legitimate given the framework of chiral effective field theory with its limited range of applicability. On the other hand, the non-perturbative ladder summations of the Brueckner calculation lead beyond chiral perturbation theory, and the input interactions in the hyperon-nucleon sector are tested for total momenta $P$ up to and even above 800 MeV, whereas the neutron Fermi momentum even at central neutron star densities around $5\,\rho_0$ does not exceed 600 MeV, the typical cutoff scale in the calculations. We can therefore assume that at least qualitative extrapolations towards such high densities can be performed.

For extrapolations to densities beyond the technically (numerically) accessible $\rho_c \simeq 3.5\,\rho_0$, we argue as follows. Around $\rho_c$, one finds clear evidence in both symmetric nuclear matter and neutron matter that the density dependence of $U_\Lambda$ has a prominent quadratic behaviour, $U_\Lambda(\rho) \propto \rho^2$. This is expected as the three-body terms begin to dominate at high densities and grow as $\rho^2$. The two-body terms at high density increase linearly with $\rho$, and this expected subleading behaviour is also observed in the numerical results already at densities 