Exploring the free-energy landscape of a rotating superfluid

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.


I. INTRODUCTION
Helium remains liquid at temperatures lower than any other substance, where it displays quantum mechanical behaviour on a macroscopic scale.Liquid 4 He undergoes a so-called λ-transition to a 'superfluid' phase as the temperature drops below 2.2K, at which point it exhibits a range of highly counter-intuitive properties including a complete absence of internal friction 1 .The irrotational nature of a superfluid results in the formation of vortex lines in a rotating container, each with quantised circulation in units of h/m (where h is Planck's constant and m is the mass of a 4 He atom), rather than a solid body rotation, such that the rotation rate of the vortex crystal matches Ω, the imposed angular velocity 2,3 .The vortex crystal is an exact equilibrium solution of the governing equations in a rotating frame, with the state observed at long time being a global minimiser of the free energy.Initial theoretical predictions of this phenomenon were verified experimentally 4 , with more recent experiments in other configurations (isomorphic to the two dimensional Euler equations) allowing for an increasingly detailed picture 5 where a relaxation onto the minimising state is preceded by the transient appearance of other (presumably unstable) vortex crystals 6 .
Campbell and Ziff 7 performed an extensive numerical study to search for the free-energy minimising state of a collection of identical point vortices for various numbers of vortices, 1 ≤ N ≤ 30 (and some special cases with larger N ).In their approach, Campbell and Ziff converged vortex crystals (relative equilibria of the equations of motion) by initialising point vortices on a series a) andrew.cleary@ed.ac.uk b) jacob.page@ed.ac.uk of concentric circles and performing gradient descent on the free energy to iteratively update the vortex positions until a minimum was reached.This numerical approach should be contrasted with other analytical approaches to find vortex crystals, which often rely on assuming certain symmetries or geometric constraints [8][9][10] and hence have often been restricted to fairly modest values of N .Recent reviews 11,12 have highlighted the need for a systematic study of the low-energy equilibria that emerge as the number of vortices increase.We tackle this problem here by combining numerical methods typically employed in the dynamical systems approach to turbulence 13 , and newer optimisation techniques that have emerged in the field of deep learning 14 .
The dynamical systems picture of turbulence envisions fluid motion as a trajectory in a high dimensional state space, pin-balling between simple invariant solutions 15,16 .These ideas gained significant traction following the discovery of a myriad of unstable travelling waves on the laminar-turbulent boundary in a pipe 17,18 and with the discovery of an unstable periodic orbit embedded in the turbulent attractor of a minimal Couette flow 19 .Finding these exact solutions relies on both (1) a sensible method for generating initial guesses for simple invariant sets and (2) a Newton-Raphson approach to converge them 13,20 .Point (1) here has been a significant constraint on applying these ideas at high Reynolds numbers, while point (2)  restricts the search to (relative) equilibria and periodic orbits -the method is not appropriate for computing connecting orbits for example.
Recent work in dynamical systems for fluids has benefited from the rapid advances in data driven methods and machine learning 21,22 .
This has included new methods for generating plausible guesses for simple invariant solutions 23 , including using deep learning architectures 24 .However, it is perhaps the develop-ment of automatic differentiation (AD) which offers the most compelling opportunities when combined with more traditional numerical simulation techniques.Time integration routines for systems of differential equations are essentially a sequence of function compositions involving elementary operations on an input (state) vector, and AD exploits this fact to compute exact derivatives via repeated application of the chain rule 14,25 .AD is therefore a powerful tool in optimisation problems, and differentiable numerical solvers can be straightforwardly combined with neural network architectures.The JAX library 25 in particular has gained significant popularity, leading to exciting advances in molecular dynamics 26 and numerical methods for computational fluid dynamics 27,28 .Gradient based optimisation has already shown great promise as a tool for generating robust guesses for unstable periodic orbits 29 .
In this paper we develop a new AD-based solver for point vortex evolution in a range of geometries, which we combine with a traditional Newton-GMRES-Hookstep algorithm to comprehensively tackle the problems outlined in the recent reviews 11,12 .We are able to assemble an extremely large (tens of thousands) collection of vortex crystals near the free-energy minima for a range of number of vortices 7 ≤ N ≤ 50, and discover new continuous families of equilibria which were previously unknown.Moreover, the flexibility of AD makes it straightforward to compute both dynamical (constant free-energy orbits) and non-dynamical pathways between nearby states to allow for detailed exploration of the freeenergy landscape.
The rest of this work is structured as follows.In §II, we formalise the two dimensional point vortex model in a rotating disk, and describe our approach to search for low energy vortex crystals.In §III, we assemble large numbers of vortex crystals for numbers of vortices 7 ≤ N ≤ 50 and discuss the new continuous families of double-ringed configurations which emerge.In §IV, we outline and apply an approach to compute energy-minimising pathways between locally stable configurations.In §V, we compute homoclinic orbits for saddles close to the energy minimising states and examine their transient dynamics.We conclude in §VI with a summary of our results.

A. Physical problem
We consider two-dimensional, irrotational flow in a disk of radius R, in which we place N point vortices of equal circulation, Γ.The disk rotates at constant angular velocity, Ω.Each point vortex moves at a velocity due to the induced velocity from (1) all of the others and (2) the image vortices.If vortex α is located at a dimensional complex position z * = z * α = x * α + iy * α , then its image sits at z * = R 2 /z * α , where the overbar represents the complex conjugate.Lengths are non-dimensionalised with the disk radius, R, and time with a reference timescale R 2 /Γ, so that the (dimensionless) evolution equation for each vortex in a frame rotating with the disk is then: where the prime on the summation indicates that β = α is excluded and ω := R 2 Ω/Γ is the non-dimensional disk rotation rate.Equation ( 1) is equivariant under rotations about the origin of the disk, R θ : (z 1 , . . ., z N ) → (exp(iθ)z 1 , . . .exp(iθ)z N ).There is also symmetry under permutation of the vortices, which we must carefully account for when labelling unique configurations or searching for connecting orbits.
If f t (z) is the time forward map of equation ( 1), then equilibria are solutions for which f t (z * ) = z * ∀t.These are relative equilibria (REQ) in the lab frame, rotating with the disk at angular velocity ω.These REQ, or 'vortex crystals', are stationary points of the free energy, which plays the role of a Hamiltonian in the rotating frame.We have scaled the free energy (2) to match the definition in Campbell and Ziff 7 (note that their 'rotation rate' variable is related to ours via ω CZ ≡ 2πω).The equations of motion follow from ẋα = (1/4π)∂ yα F, ẏα = −(1/4π)∂ xα F. The experimentally observed state in a superfluid is expected to be the global minimiser of F 3,7 .At this point, searching for free-energy minimising solutions for a given N also requires the selection of a disk angular velocity, ω, and hence documenting the solutions for a range of ω becomes a daunting task.Campbell and Ziff 7 have shown that the problem is substantially simplified if the image vortices are neglected, which can be justified for 'moderate' rotation rates ω (e.g. at some point beyond the critical ω c at which the equilibrium of interest can be maintained in the disk).The reasoning is that in order to match an increasing disk angular velocity, the vortices must increasingly bunch together near the origin so that the streamlines of the induced velocity are approximately circular at the boundaries and the images play only a small perturbative role enforcing the no-penetration boundary condition.In the absence of images, an equilibrium at one value of ω can then be trivially rescaled to another value via the simple relation Without images, the free energy is simply, where the first term is the free-space Hamiltonian, While ω can now be set arbitrarily, it is helpful to use an ω-independent label to identify and order the relative equilibria.A suitable observable is proposed in Campbell and Ziff 7 , who subtract the free energy of a continuum model, F c , evaluated at the same ω to find where the constant b = 0.74875.Equation ∆f is independent of ω at equilibrium as ω|z i | 2 = constant.We focus on the image-less system for the remainder of this paper, apart from §III B where boundary effects are reintroduced to examine their impact on a certain class of free-space vortex crystals.

B. Relative equilibria guess generation
With the removal of images, it is convenient to work in the stationary lab frame, where the vortex velocities can now be simply expressed in terms of the free-space Hamiltonian H 0 : ẋα = 1 4π We solve equation (5) using our fully differentiable point vortex solver jax-pv.Time integration is performed with the symplectic second order Runge-Kutta scheme at a fixed time step of δt = 10 −3 .The numerical solver is built on the JAX library 25 , which allows for efficient computation of the gradient of the time-forward map f t (x), where x := (x 1 , y 1 , . . ., x N , y N ), with respect to the initial positions of the vortices.The JAX framework allows us to also leverage efficiency benefits such as just-in-time compilation and auto-vectorisation.The ability to minimise loss functions which implicitly depend on the time-forward map is central to this work, and we use it to construct a wide range of REQ candidates which are then converged with a Newton-GMRES-Hookstep solver.
Our fully differentiable formulation allows for explicit 'targeting' of REQ with specific physical properties, and we initially seek candidates by minimising a balance of two loss functions.The first contribution acts to minimise the free energy, while the second demands that the inter-vortex distances αβ (t) = x α (t) − x β (t) are unchanged after a suitably long simulation: We combine the loss functions in the following manner: where the contribution from the inter-vortex distances is evaluated at two points along the orbit to encourage fixed inter-vortex distances throughout the computation (constraints at additional times were also tried with minimal impact on performance) and the integration time is fixed at T = 10.Prior to optimisation each guess is initialised from a uniform distribution x α , y α ∼ U [0, √ N ].We re-centre the configuration such that the centre of vorticity coincides with the origin, and we rescale the configuration throughout optimisation such that its average rotation rate is fixed at ω (t) = π/2T .The continual updating of ω ensures that the pattern traces out a quarter rotation regardless of how the inter-vortex distances are adjusted in gradient descent.Converged crystals have ω|z i | 2 = constant so can be rescaled arbitrarily and compared to known results using the the ω-independent observable ∆f (4).
During optimisation the gradient ∇ x(0) L is passed to an Adam optimiser 30,31 with initial learning rate 10 −2 .We then take guesses for which L T I ≤ 10 −3 within a maximum N opt = 2500 optimiser steps and attempt to converge them with a Newton-GMRES-Hookstep solver for relative equilibria 13 .Over the course of the optimisation procedure we anneal κ from 1 to 0 in increments of 1/N opt .In general, the 'good' guess criteria L T I ≤ 10 −3 is achieved while κ is still close to 1 (typically after approximately 150-300 optimiser steps),where the loss function (8) is dominated by the free energy.However, the gradual and slow introduction of the L I terms help to reduce the number of optimiser steps required, especially at higher values of N , and yields a much richer family of REQ than the L F term alone, which tends to produce the global minimiser.This can be thought of as an explorationexploitation trade-off, as in the Reinforcement Learning literature 32 .
At the Newton-GMRES stage, we consider guesses to be converged when where θ = −ωT .After the Newton convergence of a large set of vortex crystals, the final step in the process is to extract all unique REQ.This must be done carefully to avoid classifying symmetric copies (e.g.under rotation or permutation) as unique crystals.When comparing two configurations, both are first scaled to the same value of H 0 .This is done by scaling one of the configurations x → γx, with where ∆H 0 is the difference in Hamiltonian of both configurations.Then, all N (N − 1)/2 inter-vortex distances ..,N,j=i+1,...,N are computed and sorted for each of the two configurations, where N is the number of vortices in each configuration and k = 1, 2 is a label for the first and second configuration, respectively.The configurations are classed as symmetric copies if This uniqueness condition is different from the 'ring labelling' condition from Campbell and Ziff 7 , as we do not restrict ourselves to REQ based on concentric rings.

C. Targeted search
As N increases we anticipate the number of REQ will grow dramatically, hence we augment the random search method described above with another approach to search for 'missing' low energy states starting from a library of converged solutions.The appropriate loss function here is which seeks to converge onto REQ lower than some target free energy F < F * via the sigmoidal contribution σ(•), which (approximately) will 'fire' if F > F * .We set the hyperparameter ν = 100, the large value favouring lower energy states.We apply a recursive approach to optimise for (12), beginning from a converged solution, x * , at F = F * , where we select F * to coincide with the free energy of a previously converged state.Typically we do this for all states with a corresponding ∆f below some threshold value ∆f max .We take the converged solution and perturb each vortex position according to (x α , y α ) → (x α + εx α , y α + εy α ), where ε = 0.1 and x α , y α ∼ U [0, 1].We do this 20 times to generate 20 new candidate equilibria, before then performing gradient descent on (12) with an initial δ = 10 −1 .For this REQ data set augmentation we find that the simpler AdaGrad optimiser 33 is particularly effective.Our initial learning rate here is η = 0.1.Since we begin this optimisation process with a small perturbation to a converged equilibrium, the initial value of the inter-vortex loss (our usual convergence metric) is typically L T I 10 −5 .To ensure that we generate new candidates and not simply repeats of our starting solution, we require that the optimiser performs at least 25 steps.During this process the loss component L T I increases while the vortex positions adjust to reduce the large sigmoidal contribution to L S .We pass solutions for which L T I ≤ 10 −4 to a Newton solver.Finally, new unique solutions with ∆f < ∆f max are identified via (11).

III. VORTEX CRYSTALS IN THE IMAGE-LESS SYSTEM A. Crystal distributions
We begin with an initial test case at N = 7, where it has been rigorously established that there are exactly 12 REQ 34 .In this case we seek to converge all known solutions and hence do not wish to restrict ourselves to low free energies, so we use a version of loss function (8) with κ = 0 throughout optimisation.The 'random search' approach yields 11 of the 12 known solutions, with only the high-F co-linear state requiring special attention -we find this solution by initialising our vortices on a straight line -presumably because its basin of attraction in the optimisation does not favour the random vortex initialisations.The REQ solutions and corresponding free-energy labels are reported in figure 1.To generate the first 11 solutions we started with 1000 random initialisations, from which we converged 188 states in total.The most frequently computed REQ was the energy minimising crystal at ∆f = 0.10749 (54 times), while the REQ at ∆f = 4.96133 was only computed once.Setting κ = 0 is an important modification when attempting to compute REQs at high energies: if κ is annealed as previously described, then we are unable to find the crystal at ∆f = 3.49227, even after 8000 random initialisations.This is likely due to the very similar structure at the slightly lower value of ∆f = 3.49226, which is preferred by the L F term in (8).
We now consider large numbers of guesses for N = 10, 11, . . ., 20, as well as N = 30 and N = 50.A summary of our initial random search using loss function ( 8) is reported in table I.We searched extensively at N ∈ {10, 12, 16, 20, 30, 50}, while other vortex numbers were considered to explore the possible appearance of degenerate families of vortex crystals that we discuss later in §III B. Notably, we are able to generate very large numbers of viable guesses via automatic differentiation (success rates at this stage are almost all 0.98).As is often found in turbulence 20,35 , the success rate of the Newton algorithm is a bottleneck, particularly at higher N , although even at N = 50 we are still able to converge around a quarter of our AD guesses.
The number of unique convergences identified (as classified by equation ( 11)) is also described in table I.These numbers should be interpreted in the context of the raw number of samples considered for each N .For example, the large numbers of guesses considered at N = 10 or N = 12 and the very low 'success rate' in the unique solutions column indicates that we may have found most (potentially all) REQ in those cases, while the high success rates at N = 30 and N = 50 indicates that we have found only a small subset of all available crystals.
We examine the converged solutions further in figure 2, where we plot both the number of Newton convergences and the number of unique convergences below ∆f = 1 against the number of samples over the course of the computation for N ∈ {10, 20, 30}.At both N = 10 and N = 20 the number of unique convergences is found to level off while the number of Newton convergences remains roughly constant.The plateauing of the number of unique solutions is observed for roughly an order of magnitude more samples when N = 20 compared to the N = 10 case.At higher N = 30 we continue to observe roughly constant linear growth in the number of unique solutions even at O(10 5 ) guesses.
We also explore the distribution of states in as a function of ∆f for N ∈ {10, 20, 30} in the lower panels of figure 2. Clearly the density of states increases enormously with increasing N , and by N = 30 there is some indication that there is an exponential increase in the number of REQ with increasing free energy.Note that the explicit search for low free-energy solutions in our loss (8)  is likely the reason for the sudden drop in numbers of convergences beyond ∆f ≈ 1.
Due to the indications in table I and the top-right panel of figure 2 that we are far from identifying all REQ at N = 30, we supplement our random search there with the targeted sweep described in section IIC.We set ∆f max = 0.5, below which we initially had 19 converged solutions (see inset in figure 2).Each of these 19 solutions is used to define a value of F * in loss function (12), and we generate 20 perturbed states for each, or 680 new guesses in total.Of these 680 guesses, 190 unique vortex crystals were converged, of which 17 were new solutions for which ∆f < ∆f max .Additional solutions found this way are indicated in red on the histogram in figure 2; the lowenergy sweep was necessary to find the lowest free-energy state, which differs from that reported previously 7 .Note that the Newton solver also produces a large number of states for which the converged ∆f > ∆f max -see the red detail in figure 2.
We report the four energy minimising states for N ∈ {10, 12, 14, 16, 18, 20, 30} in figure 3.This figure includes some states identified by Campbell and Ziff 7 , but with many additional solutions.In figure 3  A handful of solutions shown in figure 3 are stable centres (including, as expected, the energy minimising states), though the majority are saddles.The stability of all states is determined by computing the eigendecomposition of the Jacobian of the vortex velocities (the right hand side of equation ( 1) without images), in a frame rotating at angular velocity ω with the crystal.Most of the saddle solutions in the vicinity of the free-energy minimum have low dimensional unstable manifolds.Stability properties of the 10 lowest free-energy solutions for all N considered are reported in appendix A, along with additional details of the states themselves.All states possess a neutral eigenvector (eigenvalue λ 0 = 0) corresponding to a rotation through the origin.
For some values of N we observe a large number of 'unique' states (as identified by the symmetry-invariant observable (11)) with identical values of ∆f .These distinct states all possess an additional neutral eigenvector, and appear to be smoothly connected to one another -i.e.given any two crystals with equal ∆f which are 'close', we can always interpolate to find another intermediate equilibrium between them.Vortex crystals exhibiting this behaviour are shown in red in figure 3, where we are actually reporting only one of a large number of converged solutions at that particular ∆f .To our knowledge these continuous families of vortex crystals have not been documented before, and we now explore them in more detail.

B. Degenerate double ring solutions
We examine both 'neutral' eigenvectors (λ (1,2) 0 = 0) for the energy minimising state at N = 16 in figure 4, which is the lowest N at which we found a continuous family of solutions.In the top left panel of that figure we see that one of these eigenvectors corresponds to the rotational symmetry in the system (red arrows) as expected.The second eigenvector, which is identified with blue arrows, TABLE I: Summary of computations for various numbers of vortices.Note that the success rate shown in brackets in the columns is defined relative to the raw number in the previous column.The success rate in the final column should be interpreted differently to those preceding it -a low success rate here for a large number of samples suggests that most, or potentially all, unique states may have been obtained.is markedly different: in the co-rotating frame it appears that the inner ring rotates while the outer ring wobbles in a non-axisymmetric manner.The dependence of the free energy on the logarithm of the relative distances between vortices does in principle allow for more complex behaviour than uniform rotation/translation: Consider a continuous family of solutions parameterised by some arc-length s (formally, a group orbit).The condition that Clearly, this condition holds for pure rotation as d s ij = 0 ∀i, j, but non-trivial combinations of vortex displacements are also allowed.As a further confirmation of a continuous family of crystals, we also verified that the second directional derivative of the free energy along the relevant neutral eigendirection was zero.In the free-energy landscape, we can think of the family of solutions as a smooth, continuous valley along which the crystal structure varies in a non-trivial manner (in contrast to the rotational symmetry).Computationally, we are able to take any member of a continuous family of vortex crystals and move smoothly through the valley by perturbing it along the relevant neutral eigendirection, x * → x * + εx 0 .We set ε = 5×10 −4 and construct a sequence of connected states {Px * i } by re-converging the perturbed solution with a Newton solver before repeating the process, performing a pullback [36][37][38] P along the rotational symmetry direction (see appendix B 1) to account for any drift in the rotational-symmetry valley introduced by the Newton solver.The connected sequence of states is parameterised by an arclength, s n := s n−1 + Px * n − Px * n−1 2 , and we stop the calculation once we reach a symmetrically equivalent copy of the starting solution.
The family of solutions obtained via arclength continuation at N = 16 is reported in the remaining panels of figure 4. In the top-right panel of the figure we measure the sum of the differences in inter-vortex dis-tances, relative to the initial state, identifying a 'period' of ∆s ≈ 0.5.Note that we start this continuation from a state with mirror symmetry (see red markers in lower panels of figure 4).The vortex configurations that make up one arclength-period of the continuous family are visualised in the lower panels of figure 4, where we choose a reference frame in which the angle of one of the innerring vortices is fixed.We observe that the outer ring vortices move smoothly relative to the inner ring until a permuted copy of the starting solution is reached.The permutations are shown in the final panel of the figure: the vortex labels of the inner ring are rotated through one 'click' while the outer ring vortices are permuted by two 'clicks'.Note that the vortices on the inner ring are slightly perturbed as we move through the family.If the arclength process is continued through many cycles, the outer ring of vortices is observed to move through a smoothed approximation to a pentagon -matching the number of vortices on the inner ring -which is shown with the thick grey line in the bottom left panel of figure 4.
We searched our library of solutions for other examples of continuous families of crystals, and all configurations exhibiting this behaviour are summarised in figure 5.The The valley has a periodicity of approximately ∆s ≈ 0.5 and we start from a mirror-symmetric state (see red markers in bottom panels).(Bottom left) The movement of the vortices along one cycle (∆s ≈ 0.5) through the family is shown in colour.For ease of visualisation we fix the location of one of the inner ring vortices as we move through the family.Note the small motion of the inner ring highlighted in the inset.The grey outline in the background shows a continuation of the family over many cycles, and traces out a pentagon shape to match the vortices on the inner ring.(Bottom right) Arrows show the permutation of vortices that has occurred after one cycle.
states are always 'double ring' configurations, and there may or may not be a central vortex.In all cases, following the continuous valley of solutions results in the outer ring moving relative to the inner on a smooth approximation to a geometric shape specified by the number of vortices on the inner ring.For example, if there are six vortices on the inner ring, the outer vortices approximately move around a hexagon.The continuous families of solutions identified for certain double ring configurations exist for the unbounded problem (no images).It is natural to query what happens when the images are included, particularly at lower rotation rates, where the proximity of the image vortices will have a non-negligible impact in the problem.To examine the importance of the images we apply the random search algorithm from §II B to target strict equilibria in the non-inertial frame.As the vortices are now constrained to a unit radius disk, the guesses are initialised from a uniform distribution x α , y α ∼ U [−0.7, 0.7].As ω is increased, the vortices in a free-energy minimising equilibrium are brought closer to the axis of rotation so that the induced velocities are sufficiently high to match the imposed rotation rate.To compensate for these larger induced velocities, and to be consistent with our scaling in the unbounded domain, we set the integration time FIG.6: Examples of the discrete, unique energy minimising equilibria at ω CZ = 20 (left) and ω CZ = 40 (right).The contraction of the vortices towards the axis of rotation for the higher rate of rotation is evident when comparing the two plots.All 9 unique equilibria at ω CZ = 20 are shown in various colours.These 9 equilibria are arranged in 5 'bands'; any equilibria in the same band are very close to overlapping.All 1128 unique equilibria computed at ω CZ = 40 are shown in grey, illustrating the approach to the unbounded limit as ω CZ is increased.
T = 2π/ω.We also fixed the annealing parameter κ = 1 as we were only interested in computing the energy minimiser in the rotating disk domain.Otherwise, the random search algorithm and optimisation parameters are unchanged from §II B.
In the full problem with images, we find that the continuous families are replaced by a discrete of states with equal free energy.For example at N = 16 with ω CZ = 20 (the critical rotation rate is 19 < ω c CZ < 20 7 ) we converged 64 energy minimising solutions from 1000 samples, of which 9 were unique solutions.Increasing slightly to ω CZ = 22, we converged 36 energy minimising equilibria from the same number of samples, of which 12 were unique solutions.The very low rate of unique convergences suggests that we may have found all unique equilibria at these rotation rates.At a much higher value ω = 40, the images have a much smaller impact in the problem and we find a very high rate of unique convergences, generating 2020 energy minimising equilibria from 4250 samples, of which 1128 are unique solutions.The states found at ω CZ = 20 and ω CZ = 40 are summarised in figure 6.At the lower rotation rate the outer ring vortices for the various configurations appear to be approximately uniformly spaced and arranged in a series of five 'bands' while some locations (relative to the inner ring) are completely inaccessible.At the higher rotation rate we obtain a very close approximation to the continuous family reported in figure 4.

IV. ENERGY-MINIMISING PATHWAYS BETWEEN RELATIVE EQUILIBRIA
In a number of cases we have found that the freeenergy minimising solution could actually be any one of a continuous family of states -e.g. this is true for all N considered in figure 5 apart from the 17 vortex crystal.Experimental realisations of these configurations [4][5][6] typically show a sequence of unstable states en route to the minimising solution, or even transient behaviour at long time 6 .It is natural to question whether there is a 'preferred' state within the continuous family which is more likely to be observed in an experiment.More generally, it would be useful to map out and visualise the free-energy landscape around local minima, as global minima with a particularly shallow free-energy barrier may not be robust to small perturbations to the system, while a local minimum with high barrier energy might be observed in practice.To explore the freeenergy landscape around a particular vortex crystal, we adapt a form of the doubly nudged elastic band (DNEB) method 39 , to find non-dynamical, free-energy minimising pathways between REQ.This method has been shown to work successfully in multiple settings, such as structural self assembly 40 and rearrangements of the Lennard-Jones cluster 39 .The DNEB methodology should be contrasted with the simple approach adopted by Campbell and Ziff 7 , who searched for minimising pathways with the constraint that the pathway from one state to the other coincides with a monotonic outward radial motion of one vortex from an inner ring to an outer ring, which we will see is often not the case.
To approximate the energy-minimising pathway between two REQs, x * a and x * b , the DNEB method seeks an energy-minimising chain of N I intermediate states {x i } i=1,...,N I .For the point vortex system in the unbounded domain, we define a 'potential energy' as the sum of the free energies along the chain: Minimising V alone would likely result in states 'sliding down' to one of the two minima 39 , and to counter this effect we add a fictional elastic band potential which inserts frictional high-dimensional springs between adjacent states: where x * a ≡ x 0 and x * b ≡ x N I +1 .For standard Hookean springs, we would write L G (x i , x i−1 ) = x i − x i−1 2 ; in the vortex pathway problem it is necessary to modify this term to account for both the rotational and permutation symmetries in the system, which we describe below.Energy minimising pathways are then determined by minimising an overall loss function To account for the permutation symmetry we define a permutation-invariant observable by centering a twodimensional Gaussian field over each vortex.For an Nvortex configuration x * , this observable is defined as (17) The hyperparameter σ is set as σ = 10 −2 A/N , where A = L 2 is the area of the computational grid on which G is evaluated; typically we set L = 1.25 max |x * | and evaluate G on a square grid with N x = N y = 64 grid points.An example of observable (17) is reported in figure 7 for the N = 16 system.The vortex configuration is scaled so that it rotates at ω = π/2T with T = 10 to match our convention from §II B. The Gaussian observable is then used to define the spring term appearing in equation ( 16): where the state x i−1 has been rotated through some angle θ.The angle used in between each pair of states on the chain is determined simply via Putting this all together the overall loss function for finding minimum energy pathways between states ( 16) reads The two states x a and x b on either side of the chain are kept fixed throughout the optimisation process.We solve the outer optimisation problem (20) using an Ada-Grad optimiser with a fairly aggressive initial learning rate η = 0.5, and classify our chain as converged when the relative difference between sequential iterations of the optimiser was less than 10 −6 .At each outer optimisation step we must also solve the inner optimisation problem (19) for the crystal angles, which is done in O(10) steps with an Adam optimiser (learning rate η = 10 −2 .Throughout we set the spring constant to k = 10 −2larger values caused the intermediate states close to the REQs to clump together-and the number of interpolating states to N I = 200.Further details of the optimisation, particularly how gradients are used in the DNEB method, are provided in appendix B 2. To initialise the chain we linearly interpolate between the two vortex crystals x * a and x * b , where we rotate x * b such that it maximally overlaps with x * a .This is done as a pre-processing step rather than as part of the chain optimisation, similar to the approach in Goodrich et al. 40 . We now apply the methodology outlined above to find minimising pathways into REQ 16 c 1 in the N = 16 case, looking for pathways from the two next lowest energy centres (see figure 3).These pathways are reported in figure 8.This case was also considered in Campbell and Ziff 7 but their pathways featured discontinuities in the free energy due to the constraint on how the optimisation was performed, i.e. radial motion of a vortex between rings.No such discontinuities are found here, and the smooth routes from the nearby local minima 16 c 2 and 16 c 4 both pass through unstable REQ (these states on the chain were converged in a couple of Newton steps).The barrier energy from 16 c 2 is considerably higher than that from 16 c 4 (note the correspondence of the pattern 16 c 2 with that seen in self assembly experiments of Grzybowski, Stone, and Whitesides 6 ).
The state at 16 c 1 is the simplest example of a degenerate continuous family of solutions discussed in §III B. While one example from the family is considered in figure 8, we can use our methodology to determine whether there is a preferred state when computing energy minimising chains from nearby local minima.To do this, we initialise 50 states uniformly along one period of the valley of solutions (see discussion in §III B) and compute pathways from the local minimum 16 4 to these 50 states.The 50 'final' states from the solution valley are shown in grey in figure 9. We then assess whether the pathway has entered the solution valley before continuing along the constant free-energy surface to the final state.We determine an 'entry point' by selecting the last interpolation state on the pathway which is within a relative difference of 10 −5 of ∆f ( 16 c 1 ), using the Newton solver to then converge this state to a REQ.The 50 entry points determined this way are shown in colour for all pathways considered in figure 9, the tight clustering of points clearly indicating a preferred state.If no preferred entry state existed, we would expect the coloured entry states to be spread evenly along the grey final states.Instead, the coloured entry states are clustered around one preferred entry state.The entry states become even more tightly clustered around the preferred entry state if the tolerance on the relative difference of ∆f (16 c 1 ) is increased to 10 −4 .The energy pathways optimisation does not involve any dynamics and is easily deployed at much higher N .For example, we consider a pathway between centres in the N = 30 system in figure 10, which as expected passes through an unstable saddle before a very slow descent into 30 c 1 .We can adapt the methodology here to explore the dynamics connected to the weakly unstable saddles that deform the free-energy landscape around the local minimiser; modifying the Gaussian-based loss (16) to find homoclinic orbits, exploiting the ability of the AD implementation to efficiently differentiate through trajectories.

V. HOMOCLINIC CONNECTIONS
Unlike convergence of many simple invariant solutions (equilibria, REQ, periodic orbits) there is no generally adopted methodology for computing heteroclinic or homoclinic orbits.So-called shooting methods have found some connections in plane Couette flow 41 , while a recently-proposed variational method has been used to compute connections for the one-dimensional Kuramoto-Sivashinsky equation 42 .Here, we outline an optimisation method exploiting a fully differentiable solver.

A. Candidate orbit generation
Generally speaking, we find only a single REQ at any specific energy level (see table II in appendix A), hence we restrict our search to homoclinic orbits.To generate candidate trajectories we perform a simple recurrent flow analysis 20 to identify returns to a starting equilibrium x * .Given the large number of vortices we anticipate that connections are likely to be between permuted copies of x * , hence we use the permutation invariant observable (equation 17) and modify the associated loss so that is our measure of similarity used to flag possible connecting orbits.We seed the recurrent flow analysis calculations with initial conditions x 0 = x * + εx , where x ∈ E u (x * ) (the unstable subspace of the REQ being considered) and we write where {v i } are unstable eigenvectors of x * , before normalising |x | = 1.We set the constant ε = 10 −5 .A coarse initial search with the constants c i ∈ {0, ±1} is performed, and we use the same discretisation grid for G as described in §IV.We searched over a time horizon t ∈ [0, 1000], computing R(x 0 , t) every time unit.Cases where R ≤ 0.1 are saved as candidates for connecting orbits.
FIG. 9: Common entry point (colours) into continuous family of solutions (grey), defined via a tolerance of 10 −5 on the relative difference between ∆f of the entry point and ∆f (16 c 1 ).The grey crosses identify 50 equally-spaced equilibria from with the continuous family of solutions that defines REQ 16 c 1 .An energy-minimising pathway was computed from REQ 16 4 to each.The non-dynamical pathway always approaches the configuration identified in red (which is a mirror symmetric state) before moving along the valley to the respective final state in grey.The clustering of the entry points for each pathway (coloured from blue to red0 is so tight that only the orange-and red-coloured states are visible.1 .A saddle point from our library of solutions is again found along the pathway (see details of these solutions in appendix A).

B. Convergence
Suitable guesses are converged using gradient-based optimisation.
First, a modified Jonker-Volgenant variant 43 of the Hungarian algorithm 44 is used to compute vortex permutations at the end of the connecting orbit.The Hungarian algorithm is a combinatorial optimisation algorithm that solves the linear assignment problem in O(N 3 ).The linear assignment problem in this setting is P = arg min P trace(PC), (23)   where P is a permutation matrix which permutes the rows of a matrix C whose elements C ij are the Euclidean distances between the i th vortex of the initial REQ state x * and the j th vortex in the final state on the connecting orbit candidate f t (x 0 ).Then, for the suitable candidates, f t (x 0 ) ≈ Px * .Once we have identified a candidate connection and determined the required vortex permutations, the second stage is to converge the connecting orbit.As the permutation symmetry is now respected, we minimise the Euclidean loss function: The idea here is that we must remain in the unstable subspace of x * by fixing the perturbation size ε, but the specific direction is unknown.We minimise (24) at fixed integration time t (discussed further below) using an Adam optimiser with an initial learning rate of η c = 10 −3 to determine the vector of constants c, and an Adam optimiser with an initial learning rate of η θ = 10 −2 to determine the optimal θ.The initial integration time t is set by the criteria R ≤ 0.1 from the candidate orbit generation step.We then carefully increase t throughout the convergence process.Each time t is increased, the connecting orbit is converged to a smaller tolerance level.In practice, we increase the integration time by 1% 10 times, while logarithmically decreasing the tolerance from 10 −3 to 10 −5 .We allow a maximum number of 1000 optimiser steps every time we converge the connection to a smaller tolerance.When L t E < 10 −5 , we consider the connection converged, which we confirm by verifying that the dynamics are linear in the vicinity of x * .
As a proof of concept, we successfully computed the only connection in the integrable 3-vortex system 45 , which is a homoclinic connection between the collinear state via a permutation of vortices.We then found a homoclinic connection in the non-integrable 4-vortex system, also between two permuted collinear states (with a two-dimensional unstable subspace).This connection is reported in figure 11, where we show both the initial evolution of L G that flagged the guess along with the 'converged' L G .The vortex evolution is shown in a co-rotating frame, the two (initially) inner-most vortices moving outwards along straight lines as the two outer vortices loop back to replace them.Many of unstable saddles reported in §II B (see also appendix A) have low dimensional unstable manifolds, and often feature along free-energy minimising pathways between minima as observed in §IV.The method outlined above for homoclinic orbits gives us a further tool with which to explore the range of dynamical events we may expect to see in a dissipative relaxation onto the free-energy minimising state.We report one such example in figure 12, where we find a homoclinic orbit for a weakly unstable saddle (REQ 30 4 ) in the N = 30 system, with a 2-dimensional unstable subspace.

VI. CONCLUSION
In this paper we have introduced a new computational approach to efficiently identify large numbers of relative equilibria in point vortex dynamics, focusing on free-FIG.12: A homoclinic connection for REQ 30 4 is achieved via a permutation of vortices, and is shown here in a rotational symmetry reduced frame.The connection evolves (grey outline) from REQ 30 4 (red crosses) to its permuted copy (blue crosses).
energy minimising states in a rotating cylinder.Using a combination of modern optimisation techniques with a fully differentiable solver, we were able to assemble tens of thousands of low-energy vortex crystals for a wide range of vortex numbers N , with some evidence that we may have found all low-energy states for lower values N 20.
As part of the search for energy minimisers we discovered a new set of double-ringed continuous families of vortex crystals in unbounded domains, in which the outer ring can be smoothly deformed relative to the inner.These crystals are often the energy minimiser when they occur.The introduction of a solid wall breaks up the continuum of solutions into a set of discrete states with identical free energies, though the number of such states increases dramatically as the disk rotation rate is increased and the role of the image vortices diminishes.
Finally, we developed computational methods to explore the free-energy landscape in the vicinity of the global minimiser, finding both non-dynamical, minimumenergy pathways and dynamical homoclinic orbits for unstable saddles.When the energy minimiser is a continuous family of states, the non-dynamical pathways from nearby local minima indicate a common entry point into the family, and are suggestive of a preferred state experimentally.Future work could exploit the differentiability of the solver here to directly connect simple invariant solutions for point vortices to more complex data e.g. to assess the role of vortex crystals in the decay of twodimensional turbulence 46 , or adapt the connecting orbit search to more complex dissipative systems like the Navier-Stokes equations.sends every point on our manifold of possible point vortex configurations M to a representative point on a slice of the manifold M ⊂ M. The N vortex system has a 2N dimensional manifold and a continuous symmetry group G (rotation).This slice is then a 2N − 1 dimensional submanifold which intersects all the group orbits of G of x transversally and at most once.The slice is defined as the hyperplane normal to the tangent of G at the slice fixing point x : The group tangent at x is defined as τ = Tx , where T is the generator of infinitesimal rotations for our system, such that G(θ) = exp (θT) is an element of G. T is a 2N × 2N matrix as: .
The pullback operator P is defined by a group element G(θ), and it remains to find an appropriate expression for θ(x).The slice condition in equation (B2) can be simplified to (G(θ)x) T Tx = 0, (B3) as T T = −T is antisymmetric so that x T Tx = 0. We now use equation (B3) to find an expression for θ(x).

DNEB method
Gradients of the loss in (20) are 'nudged' to counteract two common issues: 'sliding-down' and 'corner-cutting'.'Sliding-down' occurs when the intermediate points move towards the minima.The nudging approach counteracts this by only considering the portion of the gradient of V orthogonal to the direction of the pathway.'Cornercutting' occurs when the path is highly curved, and the spring potential causes the intermediate states to cut out this highly curved bend.The nudging approach counteracts this by only considering the portion of the gradient of Ṽ parallel to the direction of the pathway.The discussion here follows that given in Trygubenko and Wales 39 and the supplementary information for Goodrich et al 40 .
Formalising this, let τi be the unit tangent vector to the chain at x i : τi =      (j − i)(x j − x i ) x j − x i If exactly one of the neighbours j satisfies F(x j ) > F(x i ) x i+1 − x i−1 x i+1 − x i−1 Otherwise (i.e. is a local extremum) (B4) Let ∇ i be the gradient with respect to the state vector of x i .
The nudged elastic band (NEB) method uses to optimise L P .The doubly-nudged elastic band (DNEB) method only projects out some of the g⊥ i term so that is instead used to optimise L P .

FIG. 1 :
FIG.1:All relative equilibria for the point vortex system with N = 7, in order of increasing ∆f .
there are new states at N ∈ {10, 14, 18, 20, 30} (identified in detail in the figure caption).Most notably, all four solutions at N = 30 in the figure are new, including a new global minimum of the free energy.

FIG. 2 :
FIG. 2: (Top) Unique convergences below ∆f = 1 (red) and the number of raw convergences including repeats (blue) as a function of the number of samples.(Bottom) Histograms of the distribution of states as a function of the free energy with bin width ∆f = 0.08.Results are shown for numbers of vortices N = 10 (left), N = 20 (centre) and N = 30 (right).The red contributions to the distribution at N = 30 indicates extra solutions found by applying the targeted search for low energy REQs using equation (12).

FIG. 3 :
FIG.3:The four energy minimising states for N = {10, 12, 14, 16, 18, 20, 30} returned from the random search method detailed in §II B (with additional states found at N = 30 via the targeted search described in §II C).The solutions plotted in red denote members of continuous solution families (see §III B).The states are labelled by their number of vortices, with a subscript denoting their ordering by ∆f , while superscript 'c' is used to identify centres.Solutions 10 3 , 14 2 , 18 3 , 20 3 , 20 4 and all crystals at N = 30 shown here are new (not listed in Campbell and Ziff 7 ), including a new global minimiser 30 c 1 .

FIG. 4 :
FIG. 4: Continuous family of crystals at N = 16 (the free-energy minimiser).(Top left) The two neutral eigenvectors overlayed on one state from the solution family.Red arrows correspond to the continuous rotational symmetry while the blue arrows correspond to the new non-trivial direction.(Top right) Sum of the difference of intervortex distances between a crystal in the family and a starting state as a function of arclength, s.The valley has a periodicity of approximately ∆s ≈ 0.5 and we start from a mirror-symmetric state (see red markers in bottom panels).(Bottom left) The movement of the vortices along one cycle (∆s ≈ 0.5) through the family is shown in colour.For ease of visualisation we fix the location of one of the inner ring vortices as we move through the family.Note the small motion of the inner ring highlighted in the inset.The grey outline in the background shows a continuation of the family over many cycles, and traces out a pentagon shape to match the vortices on the inner ring.(Bottom right) Arrows show the permutation of vortices that has occurred after one cycle.

1 FIG. 5 :
FIG.5: Continuous energy minimising valleys for N ∈ {16, 17, 18, 20, 21}.These were the only states exhibiting this behaviour in our solution library.Here we show in red a mirror symmetric state in the family, along with the states reached via arclength continuation in grey.Note the tendency of the outer rings to approximately follow a geometric shape set by the number of vortices on the inner ring.In all examples we use a reference frame in which the location of one of the inner ring vortices is fixed.

4 FIG. 8 :
FIG.8: Energy minimising pathways from local minima in the N = 16 system into the global minimum.Two pathways are shown here on the left and right of the figure respectively, where the the centres correspond to solutions reported in figure3.Both pathways move through a saddle point which is also a state identified in the random search (see details of these solutions in appendix A).

FIG. 10 :
FIG. 10: Energy minimising pathway at N = 30 pathway from the next lowest energy centre REQ 30 c 5 to the minimiser 30 c1 .A saddle point from our library of solutions is again found along the pathway (see details of these solutions in appendix A).

1 FIG. 11 :
FIG.11:(Top)The evolution of min θ L G over approximately 23 time units for the 4 vortex system, for both the candidate connection (blue) and the converged homoclinic connection (red).(Bottom) Snapshots at various stages of the converged homoclinic connection in a co-rotating frame.The result of the connection is that the inner vortices are swapped with the outer vortices.
vector x is a 2N dimensional vector (x 1 , y 1 , x 2 , y 2 , . . ., x N , y N ).Expanding the exponential of the above generator, one gets the rotation matrix which rotates each vortex (x i , y i ) counter-clockwise by θ: