  In recent years, neural quantum states have emerged as a powerful variational
approach, achieving state-of-the-art accuracy when representing the
ground-state wave function of a great variety of quantum many-body systems,
including spin lattices, interacting fermions or continuous-variable systems.
However, accurate neural representations of the ground state of interacting
bosons on a lattice have remained elusive. We introduce a neural backflow
Jastrow Ansatz, in which occupation factors are dressed with translationally
equivariant many-body features generated by a deep neural network. We show that
this neural quantum state is able to faithfully represent the ground state of
the 2D Bose-Hubbard Hamiltonian across all values of the interaction strength.
We scale our simulations to lattices of dimension up to $20{\times}20$ while
achieving the best variational energies reported for this model. This enables
us to investigate the scaling of the entanglement entropy across the
superfluid-to-Mott quantum phase transition, a quantity hard to extract with
non-variational approaches.
