  The equilibrium state of a superfluid in a rotating cylindrical vessel is a
vortex crystal -- an array of vortex lines which is stationary in the rotating
frame. Experimental realisations of this behaviour typically show a sequence of
transient states before the free-energy minimising configuration is reached.
Motivated by these observations, we construct a new method for a systematic
exploration of the free-energy landscape via gradient-based optimisation of a
scalar loss function. Our approach is inspired by the pioneering numerical work
of Campbell & Ziff (Phys. Rev. B 20, 1979), and makes use of automatic
differentiation (AD) which crucially allows us to include entire solution
trajectories in the loss. We first use the method to converge thousands of
low-free-energy relative equilibria for vortex numbers in the range $10 \leq N
\leq 30$, which reveals an extremely dense set of mostly saddle-like solutions.
As part of this search, we discover new continuous families of relative
equilibria (in the unbounded domain) which are often global minimisers of the
free energy. These continuous families all consist of crystals arranged in a
double-ring configuration, and we assess which state from the family is most
likely to be observed experimentally by computing energy-minimising pathways
from nearby local minima -- identifying a common entry point into the family.
Finally, we develop an approach to compute homoclinic orbits and use it to
examine the dynamics in the vicinity of the minimising state by converging
connections for low-energy saddles.
