  An accurate description of low-density nuclear matter is crucial for
explaining the physics of neutron star crusts. In the density range between
approximately 0.01 fm$^{-3}$ and 0.1 fm$^{-3}$, matter transitions from
neutron-rich nuclei to various higher-density pasta shapes, before ultimately
reaching a uniform liquid. In this work, we introduce a variational Monte Carlo
method based on a neural Pfaffian-Jastrow quantum state, which allows us to
model the transition from the liquid phase to neutron-rich nuclei
microscopically. At low densities, nuclear clusters dynamically emerge from the
microscopic interactions among protons and neutrons, which we model based on
pionless effective field theory. Our variational Monte Carlo approach
represents a significant improvement over the state-of-the-art auxiliary-field
diffusion Monte Carlo method, which is severely hindered by the fermion-sign
problem in this low-density regime and cannot capture the onset of clusters. In
addition to computing the energy per particle of symmetric nuclear matter and
pure neutron matter, we analyze an intermediate isospin-asymmetry configuration
to elucidate the formation of nuclear clusters. We also provide evidence that
the presence of such nuclear clusters influences the amount of protons in the
crust compared to protons in beta-equilibrated, neutrino-transparent matter.
