  The equation of state (EoS) of neutron matter plays a decisive role in our
understanding of the properties of neutron stars as well as the generation of
gravitational waves in neutron star mergers. At sufficient densities, it is
known that the appearance of hyperons generally softens the EoS, thus leading
to a reduction in the maximum mass of neutron stars well below the observed
values of about 2 solar masses. Even though repulsive three-body forces are
known to solve this so-called "hyperon puzzle", so far performing \textit{ab
initio} calculations with a substantial number of hyperons has remained
elusive. In this work, we address this challenge by employing simulations based
on Nuclear Lattice Effective Field Theory with up to 232 neutrons (pure neutron
matter) and up to 116 $\Lambda$ hyperons (hyper-neutron matter) in a finite
volume. We introduce a novel auxiliary field quantum Monte Carlo algorithm,
allowing us to simulate for both pure neutron matter and hyper-neutron matter
systems up to 5 times the density of nuclear matter using a single auxiliary
field without any sign oscillations. Also, for the first time in {\em ab
initio} calculations, we not only include $N\Lambda$ two-body and $NN\Lambda$
three-body forces, but also $\Lambda\Lambda$ and $N \Lambda\Lambda$
interactions. Consequently, we determine essential astrophysical quantities
such as the mass-radius relation, the speed of sound and the tidal
deformability of neutron stars. Our findings also confirm the existence of the
$I$-Love-$Q$ relation, which gives access to the moment of inertia of the
neutron star.
