Schematic diagram of our Hamiltonian Monte Carlo method. The procedure starts with a white noise initial condition field s, then uses a nested, double-mesh particle-mesh (PM) integrator to obtain the final particle phase-space coordinates (x_final and v_final). The likelihood is implemented in Python using the JAX package (Bradbury et al. 2018), which automatically implements the gradient. We compute the gradient of the likelihood and back-propagate it through the double-mesh integrator steps to the initial condition field. Using this gradient, we then do a Hamiltonian Monte Carlo (HMC) update on the initial condition field and iterate the whole procedure many times to obtain our sample of initial condition fields.