Curved vortex surfaces in four-dimensional superfluids: II. Equal-frequency double rotations

As is well-known, two-dimensional and three-dimensional superfluids under rotation can support topological excitations such as quantized point vortices and line vortices respectively. Recently, we have studied how, in a hypothetical four-dimensional (4D) superfluid, such excitations can be generalised to vortex planes and surfaces. In this paper, we continue our analysis of skewed and curved vortex surfaces based on the 4D Gross-Pitaevskii equation, and show that certain types of such states can be stabilised by equal-frequency double rotations for suitable parameters. This work extends the rich phenomenology of vortex surfaces in 4D, and raises interesting questions about vortex reconnections and the competition between various vortex structures which have no direct analogue in lower dimensions.


I. INTRODUCTION
Quantum vortices are topological excitations that play an important role in the physics of superfluids [1][2][3][4][5][6][7][8][9].Such vortices are characterised by the quantized circulation of the superfluid around a local density depletion, which is called the "vortex core".As is well-studied, in a twodimensional (2D) superfluid, the vortex core corresponds effectively to a 0-dimensional point, while for a threedimensional (3D) superfluid, the core extends into a onedimensional line or ring.As vortices are excitations, they are associated with an energy cost, but can be stabilised by either rotating the superfluid [3,4], or equivalently by engineering artificial magnetic fields [10][11][12][13].
Recently, we began in Ref. [14] to investigate the possible phenomenology of vortex structures in a fourdimensional (4D) superfluid, by studying a 4D generalisation of the Gross-Pitaevskii equation (GPE) including rotation [1].Interestingly, the extension to 4D considerably enriches the possible vortex structures as there are fundamental differences between rotations (or equivalently, magnetic fields) in different numbers of spatial dimensions.As we shall review further below, in 2D and 3D, all rotations are "simple rotations" that can be characterised by a single rotation plane and rotation frequency, while in 4D, generic rotations are "double rotations", meaning that two completely orthogonal planes of rotation, and hence two rotation frequencies, can be identified.In Ref. [14], we explored how the simplest case of a double rotation with equal frequencies can stabilise a new type of vortex structure in which the vortex core consists of two rigid orthogonal planes intersecting at a point, with no direct analogue in lower dimensions.
In this paper, we shall go further to explore what happens as the two rotation frequencies in a 4D system are made unequal.As we shall show, this can lead to 4D vortex structures with cores composed of skewed nonorthogonal surfaces which curve to avoid the expected intersection point.We shall present both analytical and numerical calculations based on the 4D generalised GPE under rotation, and we shall develop and numerically test a theory to explain the skewed vortex planes in terms of a simplified competition between the rotational energy and the hydrodynamic vortex-vortex interaction terms.For unequal-frequency double rotations, we find skewed vortex surfaces that can be lower in energy than a pair of rigid orthogonal vortex planes [14] for our system sizes and parameters.This lays the groundwork for a followon work in Ref. [15], which will apply a similar analysis to the case of equal-frequency double rotation.
Looking further ahead, we note that we are studying a minimal 4D mathematical model, which is motivated as a natural extension of the standard GPE description of 2D and 3D superfluids.In the future, it will also be very interesting to explore if similar structures can be found in more experimentally-realistic models, connecting with recent theoretical and experimental advances in probing higher-dimensional physics [16][17][18][19][20][21][22][23][24][25][26][27][28][29][30], e.g. based on techniques such as "synthetic dimensions" [18,, using which a 4D atomic quantum Hall system has recently been experimentally realized [56].More generally, the rich phenomenology of curved vortex surfaces that we have begun to explore raises the possibility of finding other exotic topological excitations, such as closed vortex surfaces.Some of our results also suggest that vortices can lose some of their individual character in 4D, as the curved surfaces that we have found do not easily decompose into two separate but intersecting vortex states, unlike in Ref. [14].This opens interesting questions, for example, about what would happen at even higher rotation frequencies, where we may expect the number of vortices to become large.
In this paper, we shall begin in Section II by reviewing the basic physics of quantum vortices in 2D and 3D superfluids.Then in Section III, we shall discuss in more detail the different possible types of rotations in four spa-arXiv:2309.08015v2[cond-mat.quant-gas]19 Sep 2024 tial dimensions, including simple, double and isoclinic rotations.As part of this section, we shall discuss the rotation planes of isoclinic rotations in detail, as this will be useful for later sections of the paper.In Section IV, we shall briefly review the 4D GPE and the physics of completely orthogonal 4D vortex planes, which were studied in 4D doubly-rotating superfluids with equal rotation frequencies in Ref. [14].We shall also introduce the numerical methods to be used throughout this work.In Section V, we shall then derive the hydrodynamic energy of a pair of non-orthogonal skewed 4D vortex planes that intersect at a point.As we then develop in Section VI, our hydrodynamic theory predicts that, in a doubly-rotating 4D superfluid with two unequal rotation frequencies, a pair of rigid vortex planes will become skewed so as to start aligning with the higher frequency and with each other.Our numerical results verify this predicted skewness at large distances, but also show that the vortex surfaces curve near the origin to avoid the intersection point.Finally, in Section VII, we shall summarize our results and discuss possible future extensions.

II. REVIEW OF SUPERFLUID VORTICES
In this section, we shall briefly review the basic properties of quantized vortices in 2D and 3D superfluids.We shall begin by introducing the GPE and reviewing how this leads to the structure of a 2D superfluid vortex, before discussing the energy of such a state within a finite system.We shall then briefly discuss systems with multiple 2D vortices, before reviewing some aspects of vortex physics in 3D superfluids.This will lay the groundwork for the discussion of 4D vortices in Section IV onwards.

II.1. Structure of a 2D Superfluid Vortex
We consider a system of weakly-interacting bosons in the absence of an external potential as can be described by the time independent Gross-Pitaevskii equation [ where ψ is the complex order parameter, m is the mass of the particle, g is the interaction strength, and µ is the chemical potential.Note that a harmonic trapping potential is also often included in the GPE, as this is present in many ultracold gas experiments [1][2][3][4][5][6][7]57]; however, for simplicity we will neglect all such effects and restrict our discussion to infinite systems and finite systems with hard-wall boundary conditions, as specified below.
From the order parameter, it is possible to directly define the superfluid density, phase, and velocity fields as ρ = |ψ| 2 , S = Arg ψ and v = ℏ m ∇S, respectively [1].Importantly, the form of the velocity field implies that the circulation of the superfluid around any closed loop C has to be quantized with the circulation being given by where the integral on the right hand side gives the change in the superfluid phase around the loop C.This would be zero if S is single-valued, but since S is a phase and its value is only defined modulo 2π, so the superfluid circulation is quantized generally as [1][2][3][4]9] where k is an integer called the winding number.Note that continuous deformations of the loop C within the superfluid will not change the integer winding number, since v varies continuously so long as ψ is non-zero.This means that k ∈ Z is a topological invariant, which will be equal to zero whenever C can be continously contracted to a point.Hence, a simply-connected superfluid (i.e. one in which all loops are contractible) cannot circulate.If a simply-connected superfluid is rotated, it therefore cannot behave as a rigid body but will instead form internal holes, called vortices, where the density goes to zero and around which the phase winds by a quantized amount [1][2][3][4][5][6][7]57].In 2D, the structure of a rotationallysymmetric vortex is described by where (r, θ) are 2D polar coordinates centered on the vortex core.The winding number k corresponds to the topological charge of the vortex, and is conventionally taken to be positive for vortices and negative for antivortices.From the above ansatz, the vortex velocity field is then [1] where θ is the unit vector pointing along the θ direction.As required by Eq. ( 3), this velocity field falls off as 1/r, and reverses direction when the sign of the winding number, k, is flipped.The angular momentum of the vortex in Eq. ( 4) is also quantized as with respect to the axis through the center of the vortex core, with ẑ being the unit vector pointing out of the 2D plane, along the z direction.More generally, in a system with axial symmetry, the angular momentum is quantized only for an on-axis vortex [2].
To complete this description of the vortex state [Eq.( 4)], the density-profile can be obtained numerically by solving the GPE [Eq.(1)].When so doing, it is common to define the uniform background density as n and then rescale ψ → √ nψ and r → ξr, where ξ is the healing length, which satisfies ℏ 2 /mξ 2 = gn = µ, and which physically is the distance over which ρ typically varies.(Note that often a factor of 1/2 is included in the definition of ξ [1].)Under these rescalings, Eq. ( 1) becomes dimensionless as which is then solved numerically [1,2] to obtain a dimensionless real-valued density profile function f k (r) = |ψ|.
While the obtained f k (r) has no closed form, it can be shown that it asymptotically vanishes towards the vortex core as f k (r) = O(r |k| ) when r → 0, and asymptotically returns to the profile of a homogeneous superfluid as f k (r) = 1 − O(r −2 ) when r → ∞ [4], with a crossover between these two behaviours around the healing length.

II.2. Energy of a 2D Superfluid Vortex
Within a hydrodynamic description [1,2,9], the energy of a 2D vortex in the absence of an external potential is made up of a kinetic contribution and an inter-particle interaction contribution where the integrals are carried out over the area of the 2D system, which we consider to be circular with radius R and hard-wall boundary conditions.Both of the above contributions can be estimated analytically by assuming that the density is zero within the vortex core, which we take to be circular with a radius of our healing length ξ, and constant otherwise across the system.(Note that other definitions for the size of the vortex core can be used [1,9]; however, the following argument is only approximate and such changes do not significantly affect the relative scaling and importance of the different energy terms.)For such a simple, so-called "hollow-core", model for a vortex, the total energy difference, ∆E, between the homogeneous and the vortex state is straightforward to calculate; firstly, in the hydrodynamic kinetic energy (Eq.( 8)) introduced above, the first ("quantum pressure") term vanishes so that the extra kinetic energy due to a vortex is given simply by [1] where n is the constant density within the system outside of the vortex core.Similarly, the interaction energy (i.e. the energy needed to make a hole in the superfluid) can be calculated as The latter can clearly be neglected for a large system with R ≫ ξ, meaning that the total hydrodynamic energetic cost of a vortex can be approximated as [9] A more accurate estimate of the energy cost relative to the uniform state can be found using the dimensionless numerical density-profile function, f k (r), in the vortexstate ansatz [Eq.( 4)].Using the grand canonical energy at fixed chemical potential µ takes care of corrections to the background density n from the core depletion.Then, using that, µ = ℏ 2 /mξ 2 , and n = N/πR 2 , in the uniform system, where N is the total number of bosons, we may write this numerical vortex energy as which importantly is the same functional form as the simple hydrodynamic estimate [Eq.(12)], up to the numerical prefactor within the logarithm.(Note that if the healing length had instead been defined including a factor of 1/2 as ℏ 2 /2mξ 2 = gn = µ, then this numerical prefactor becomes 1.46 [1].)Vortices can be energetically stabilised by rotation (or equivalently an artificial magnetic field) [3,4].In a rotating reference frame, the GPE [Eq.( 1)] becomes where L = −iℏr × ∇ is the (3D) angular momentum operator, and Ω is the rotation frequency vector [1].In 2D, we can assume that Ω = Ωẑ, and hence the energy reduction from rotation is given by ∆E rot = Ω⟨ Lz ⟩.As discussed above, vortices carry a finite amount of angular momentum and so are favoured by rotation.
To leading order, we can assume the superfluid has a constant density and neglect the depletion of the core, so that the energy reduction from rotation can be calculated as [1] (This approximation cannot be applied to the calculation of the hydrodynamic energy of a single vortex [Eq.(10)] as the 1/r 2 dependence of the integrand gives a singular contribution from the area around r = 0.) As can be seen, this term reduces the energy of a state containing a vortex for which the circulation is aligned with the rotation, and raises the energy of a state (with opposite k) that is anti-aligned with the rotation.For a vortex to be energetically stabilised, the reduction in energy must be greater than (or equal to) the cost of making a vortex within the same approximations (e.g.Eq. ( 12)).This leads to an estimate of the critical frequency of [1] i.e. this is the minimal rotation frequency needed to stablise a vortex with winding number k.Using the energy for the numerical vortex profile (Eq.( 13)) leads to a more accurate calculation for this frequency as Note also that this critical frequency will depend on any external potentials that are present, and so will be different, e.g. with a harmonic trap [2].However, in this paper, we will focus on untrapped systems with hard-wall boundary conditions, as mentioned above.

II.3. Multiple Vortices in a 2D Superfluid
As can be seen from Eq. ( 17), the critical rotation frequency is proportional to the winding number k, meaning that higher frequencies are required to stabilise vortices with higher winding numbers.However, by comparing the hydrodynamic energy with the rotation energy, it can be seen that, even at higher frequencies, it will always be energetically unfavourable (in the absence of additional external potentials) to produce a multiply charged vortex (i.e. with |k| > 1) as compared to multiple singly-charged vortices (with |k| = 1) [9].
The above argument also suggests that a pair of similarly-charged 2D vortices will interact repulsively, as it is energetically unfavourable to bring them together and merge them into a single vortex with a higher winding number.Indeed, it can be shown that, in a sufficiently large system, the interaction energy between a pair of well-separated vortices, with charges k 1 and k 2 respectively, can be approximated as [2] where ∆r is the distance between the two vortex cores.
As can be seen, this is attractive for oppositely-charged vortices (i.e. a vortex and anti-vortex pair) but repulsive for vortices with the same sign.In an infinite system, a pair of vortices can therefore continually lower their energy by moving apart, while a vortex and anti-vortex pair can lower their energy by coming together and annihilating.Note that Eq. ( 18) is derived under the approximation that the density is constant everywhere in the system, i.e. ignoring the density depletion at the vortex core.Consequently, the calculated vortex-vortex interaction energy [Eq.(18)] is only valid for separations ∆r ≫ ξ, and attempting to take the limit ∆r → 0 gives a logarithmic divergence.In reality, when a pair of vortices with winding numbers k 1,2 come together, they combine into a vortex with winding number k 1 + k 2 .One can still obtain this correct result from Eq. ( 10) if we consider the vortices to be combined once their separation is similiar to the healing length ∆r ∼ ξ.This is consistent with the constant density approximation, as the latter amounts to ignoring variations on the scale of ξ or below (except in the presence of a trap).
As the rotation frequency increases therefore above the critical frequency [Eq.(17)], it will be energetically favourable to have more and more singly-charged vortices in the system.The effectively repulsive interactions between these vortices then mean that, at high enough rotation frequencies, the lowest energy state in the rotating frame exhibits a uniform array of vortices, known as an Abrikosov lattice [2,58].

II.4. Vortices in 3D Superfluids
The above discussion can be straightforwardly generalised to describe vortices in a 3D superfluid [1,2,9].In 3D, a vortex core can be approximated as an extended 1D line, which must either begin and end on the surface of the system, or else form a closed loop within the superfluid.The former structures are often referred to as "vortex lines" or "vortex filaments", while the latter are typically called "vortex rings" [9,[59][60][61][62][63][64].As our paper is concerned with the lowest-energy vortex structures to be stabilised by rotation, we shall hereafter focus on vortex lines, although it would also be very interesting to study the analogue of vortex rings in higher dimensions.
In the simplest case, a cylindrically-symmetric vortex line in 3D can be described [9] e.g. by in cylindrical polar coordinates (r, θ, z), where we have assumed that the rotation axis lies along the z direction and that the rotation is sufficiently strong so as to align and straighten the vortex core.In the absence of an additional potential, the vortex structure is then invariant along the z direction and the dimensionless density profile is given by the radial function found numerically from the 2D GPE.Consequently, a 3D vortex line has the same velocity field as a 2D vortex [Eq.( 5)], as well as the same critical frequency (in a cylindrical system) [9].The latter point can be easily appreciated by noting that, in this case, the 2D calculation for the hydrodynamic energy follows through identically up to an overall multiplicative factor in both Eq. ( 12) and Eq. ( 15), to represent the height of the system [9].Similarly, when the rotation frequency becomes much higher than the critical frequency, many vortices enter the 3D system, and should eventually form a vortex lattice analogous to that in 2D, except with the vortex cores extended as straight lines along the rotation axis [8,9].It is also worth noting that, unlike in 2D, the shape and orientation of a 3D vortex line can depend, for example, on both the choice of rotation axis as well on the geometry and boundary conditions of the system [60].For example, in 3D there can be a competition between aligning the vortex core with the rotation axis in order to capitalise on energy reduction from rotation, and minimising the length of the vortex core in the superfluid so as to minimise the interparticle interaction energy.
Another new phenomenon that emerges in 3D is the reconnection of vortex lines [65]; when two vortex lines are made to intersect in 3D, they will generically reconnect and move apart so as to remove the intersection point.Note that there are some special cases of metastable stationary states in 3D with intersecting vortex lines [66].As we shall review later in Section IV, a key difference between 3D and 4D superfluids is that, in the latter case there can be an intersection point between two vortex planes in a stationary state which is energetically stabilised by double rotation [14].However, as we shall go on to explore in Section V onwards, we can also find stationary states with curved vortex surfaces, in which the vortex core curves spatially in order to avoid the intersection point.Analogies between these surfaces and reconnections in 4D will be further explored in Ref. [15].

III. ROTATIONS IN 4D
In order to further lay the groundwork for our discussion of 4D vortex structures in Section IV, we shall now review the different types of rotations that are possible with four spatial dimensions, comparing these with 2D and 3D systems.We shall begin by introducing the concepts of simple, double and isoclinic rotations, before discussing the possible rotation planes of 4D isoclinic rotations in more mathematical detail.As we shall see, this will be relevant when considering the effects of rotation in a generalised 4D GPE in later parts of this paper.

III.1. Simple, double, and isoclinic rotations
In two dimensions, rotations are completely specified by their centre and rotation angle.The centre is the one fixed point of the rotation, while all other points are angularly displaced about the centre by the rotation angle.Represented as a matrix, any rotation of 2D space will be given as where α ∈ (−π, π] is the angle of rotation and we are defining the origin as the centre of rotation, as we will throughout this paper.Similarly, rotations in three dimensions are commonly described in terms of their axis and angle of rotation.The axis is both the line of points fixed by the rotation and the centre about which the rotation occurs.One can equally define rotations in 3D by their plane of rotation, which is orthogonal to the axis of rotation.All rotations in 3D are just 2D rotations of their plane of rotation, with the third direction left unchanged.This is obvious from the matrix representation of a 3D rotation, which can always be brought into the following form via a suitable choice of basis.The rotation plane is left invariant by the rotation but not pointwise invariant, unlike the axis.This means that points on the rotation plane remain on it after the rotation, but are rotated about the rotation axis.Just as we can extend 2D rotations into a third direction to define 3D rotations, we may generate rotations of 4D space by extending 3D rotations into a fourth direction.In the 3D case this gave us every possible rotation, up to a change of basis.However, in 4D we can only generate a proper subset of rotations by extending our 3D definitions in this way.Members of this subset are commonly termed "simple" rotations, since they reduce to the familiar three and two dimensional cases.Simple rotations have a single rotation plane just as in the 3D case, but are centred around a plane of fixed points as opposed to an axis.This fixed plane is completely orthogonal to the rotation plane, by which we mean that every vector in one plane is orthogonal to every vector in the other.In a matrix representation, any simple rotation of 4D space can take the following form in a suitable basis Note that in 4D, there are six Cartesian coordinate planes, meaning that the rotation group SO(4) of fourdimensional space has six generators, physically describing angular momentum.For this reason, the representation of these generators (and hence of angular momentum) as spatial vectors does not work in 4D, as it does in 3D.Moreover, generic elements of SO(4) are so-called "double" rotations.These new types of rotations occur simultaneously through two completely orthogonal planes of rotation (e.g. the x-y and z-w planes), each with their own rotation angle.Represented as a matrix, any double rotation can be brought into the form by a suitable change of basis.This matrix form makes it clear that a double rotation can be thought of as two simultaneous simple rotations: in this case a rotation of angle α ∈ (−π, π] in the x-y plane, and one of angle β ∈ (−π, π] in the (z-w) plane.This means that any point on the x-y or z-w plane will remain on it but be rotated around the origin by an angle α or β, respectively.Points not on either rotation plane are rotated by an angle whose magnitude is strictly between |α| and |β| [67], assuming that |α| < |β|.Consequently, the origin is the only fixed point, as long as neither rotation angle is zero.If either angle vanishes, we recover simple rotations as a special case of double rotations.Besides simple rotations there is another very important special class of double rotations, called "isoclinic" rotations, which will play an important role in the rest of this paper.These are the double rotations where both rotation angles are equal up to a sign, such as M (α, α) and M (α, −α).They come in two types known as right handed and left handed based on the relative senses of rotation in the two planes.For example, M (α, α) is a left isoclinic rotation of the x-y and z-w planes, while M (α, −α) is a right isoclinic rotation of these planes.All left isoclinic rotations commute with all right isoclinic ones, and any rotation of 4D space can be decomposed into a product of a left isoclinic rotation and a right isoclinic rotation [67].However, this is not unique, as M = M L M R can also be written as M = (−M L )(−M R ), where M L and M R denote left and right isoclinic rotations respectively.

III.2. Rotation Planes of an Isoclinic Rotation
In later sections of this paper, we will find it useful to take advantage of various mathematical properties of isoclinic rotations in our analysis of vortices in 4D superfluids.For that reason, we shall now discuss these special types of rotations in greater detail, focusing in particular on how to identify the rotation planes of left and right isoclinic rotations respectively.
Recall that a general double rotation will rotate a vector through an angle with magnitude between |α| and |β|.However, for an isoclinic rotation α = ±β, so every vector is displaced by the same given rotation angle, meaning that there is an infinite number of rotation planes.Each of these rotation planes can be described as the span of an arbitrary vector v and its image under the rotation (i.e.either M L v or M R v), which means that every point in R 4 lies on one of these rotation planes [67].However, this does not imply that every possible 2D plane is a rotation plane (except for very special cases, as mentioned below), nor does it mean that these rotation planes are unique: any completely orthogonal pair of them can be used as a basis to define the particular isoclinic rotation.For example, from Eq. ( 23) we can see that two of the rotation planes of M (α, α) are given e.g. by the x-y and z-w rotation planes, although these are not the only rotation planes as we shall see below.This is in contrast to generic double rotations (i.e.M (α, β), with α ̸ = β ̸ = 0, π), which have only two unique rotation planes as discussed above.
Our aim is now to mathematically identify the rotation planes of a given left isoclinic rotation, which we shall denote as M L ≡ M (α, α), i.e. we chose our basis such that this particular rotation has the form given in Eq. ( 23) with α = β.As we shall see, an easy way to find the corresponding rotation planes is then to use the complex representation C 2 to represent R 4 , such that the Cartesian position vector (x, y, z, w) T is represented by (x + iy, z + iw) T .Note that the natural inner product in C 2 , given in Cartesian coordinates by (25) contains the inner product in R 4 as its real part.This means that any unitary matrix acting on C 2 will be equivalent to some orthogonal matrix acting on R 4 .However, the converse is not necessarily true as unitary matrices preserve both the real and imaginary parts of the complex inner product, while orthogonal transformations need only preserve the real inner product.Nevertheless, we can say that if an orthogonal transformation of R 4 (e.g.such as a 4D rotation) is represented by a matrix in the complex representation then that complex matrix is automatically unitary.To see this, note that the norm on C 2 agrees with the norm on R 4 , that is As any orthogonal transformation will preserve the R 4 norm, so its C 2 representation will preserve the corresponding complex norm, meaning that if that representation is a complex matrix it must therefore be a unitary matrix.
Returning to the particular case of double rotations, we see that in the complex representation, Eq. ( 23) becomes: which is indeed unitary.It is also clear that the desired left isoclinic rotation M L can simply be represented in C 2 as e iα times the identity.We will now use this to show how to construct and parametrise the rotation planes of this left isoclinic rotation using this complex representation, before also discussing the case of right isoclinic rotations.Note that, in the following, rather than Cartesian coordinates, we shall primarily use double polar coordinates (r 1 , θ 1 , r 2 , θ 2 ), which are defined by r 1 e iθ1 ≡ x + iy, and r 2 e iθ2 ≡ z+iw, such that the complex position vector becomes (r 1 e iθ1 , r 2 e iθ2 ) T in the complex representation.
In general, a 2D plane in R 4 can be defined as the set of solutions to a pair of simultaneous linear equations (e.g.x = 0 together with y = 0 defines the z-w plane passing through the origin).In contrast, in the complex representation we can define a plane using a single equation which is linear in r 1,2 e iθ1,2 and their complex conjugates.In other words, given four complex numbers (a 1 , a 2 , b 1 , b 2 ), the equation a 1 r 1 e iθ1 + a 2 r 2 e iθ2 + b 1 r 1 e −iθ1 + b 2 r 2 e −iθ2 = 0 (28) defines a plane passing through the origin, and any such plane can be defined (not uniquely) in this way.(For example, the above z-w plane can now be defined simply either as r 1 e iθ1 = 0 or equivalently as r 1 e −iθ1 = 0.) To get back to the real representation we then just take the real and imaginary parts of the complex equation.Note that we included the complex conjugates r 1,2 e −iθ1,2 in Eq. ( 28) so that the complex equation can have the same number of parameters as the two real equations.
As the rotation planes of M L are invariant under M L , to find these rotation planes we must find the equations of the form (28) that are also invariant in this way.The action of the left isoclinic rotation, M L , is given by θ 1,2 → θ 1,2 +α as introduced above, so that the image of Eq. ( 28) under M L is given by For this to reproduce Eq. ( 28), we require that, in general, either a 1,2 = 0 or b 1,2 = 0. Note that for the special angles of α = 0, π, we recover Eq. ( 28) irrespective of the values of (a 1 , a 2 , b 1 , b 2 ), meaning that every single plane is a rotation plane of M L for these cases.However, these special cases are trivial as they physically correspond to, respectively, no rotation or to flipping the direction of all axes simultaneously.Focusing therefore on the general case, we identify two possibilities: either However, we can map the latter equation onto the former by taking the complex conjugate of both sides and identifying b * 1 = a 1 , and b * 2 = a 2 .Therefore, both cases are the same and so the rotation planes of M L are given by the solutions to the equation a 1 r 1 e iθ1 + a 2 r 2 e iθ2 = 0, for arbitrary complex numbers a 1,2 .
In later sections of this paper, we will want to sometimes work in a coordinate system defined in relation to an arbitrary completely orthogonal pair of these rotation planes (which we will denote as P 1 and P 2 ), in the same way that the coordinates (r 1 , θ 1 , r 2 , θ 2 ) are defined in relation to the x-y and z-w planes.We shall therefore now go through how such a coordinate system can be defined.To begin, let P j be given by the solutions to a j1 r 1 e iθ1 + a j2 r 2 e iθ2 = 0, (32) for j = 1, 2, and let our coordinate system defined with respect to these planes be given by (r that the coefficients a jk are not all independent: once the plane P 1 is chosen, P 2 is already fixed as the orthogonal complement of P 1 .For now we will not consider this constraint, but will effectively derive it later by comparing the primed and unprimed coordinate systems.We now proceed to define the primed coordinates.Recalling that the x-y and z-w planes are defined by r 2 e iθ2 = 0 and r 1 e iθ1 = 0 respectively, we see that in our new coordinate system the planes P 1,2 should be given by r ′ 1,2 e iθ ′ 1,2 = 0 respectively.Given Eq. ( 32), a simple way to achieve this is to define our coordinates as follows Let us now determine the way in which the coefficients a jk are constrained.This can be done by noting that the unprimed coordinates are an orthonormal system; for this to also be true of the primed coordinates, we must have that the total distance from the origin is preserved, that is This condition is equivalent to requiring that the coefficients a jk furnish a unitary with φ 1,2,3 ∈ [0, 2π), and η ∈ [0, π/2).However, the U (1) factor e iφ3 is redundant since this represents a left isoclinic rotation of the planes, just like the original rotation M L but with a different angle.Since such a rotation leaves all the rotation planes invariant we may discard it, and we are left with the following The above matrix is the general expression for a member of U (2)/U (1) = SU (2), i.e a special unitary 2 × 2 matrix.We can interpret this family of matrices as the group of right isoclinic rotations [68].To see this note that Eq. ( 35) is an expression in the complex representation for rotations that commute with the given left isoclinic rotation, M L .(To see this, note that these are complex linear transformations and so they commute with i, while M L in the complex representation is simply multiplication by e iα .)However, it is also well-known that all left isoclinic rotations commute with all right isoclinic rotations [67], whereas two isoclinic rotations of the same sense (or two generic double rotations) will only commute if they share the same rotation planes, as introduced above.In going from Eq. (35) to Eq. (36) we have factored out those left isoclinic rotations which commute with M L , as they take the same form as M L in the chosen basis and so have the same rotation planes.The matrix in Eq. ( 36) therefore is a representation of the right isoclinic rotations, as these are the remaining rotation matrices that commute with M L .We can also further simplify Eq. ( 36) by factoring out the subgroup of right isoclinic rotations which take the form M R ≡ M (α, −α) in our chosen basis.Letting φ 2 = φ 1 + φ, with φ ∈ [0, 2π), we obtain cos η −e −iφ sin η e iφ sin η cos η e iφ1 0 0 e −iφ1 , (37) where the second matrix can be recognised as M R [c.f.Eq. ( 27)].Since the diagonal factor M R just corresponds to initial rotations within the x-y, and z-w planes, it is redundant in describing the coordinate transformation [Eq.(36)] from these planes to the arbitrary rotation planes, P j , of M L .We can therefore discard it such that our final expression for the general transformation is where η ∈ [0, π/2] and φ ∈ [0, 2π).This is now the general form for a coordinate transformation from a fixed pair of rotation planes (e.g. the x-y, and z-w planes) to all other rotation planes of M L , with all redundant parameters removed.Interestingly, it is clear from Eq. ( 38) that φ is undefined when η = 0, because the off-diagonal elements vanish.Moreover, a careful analysis shows that this also occurs at the other endpoint, η = π/2, as here the diagonal elements vanish and so we can eliminate φ as follows 0 −e −iφ sin η e iφ sin η 0 = 0 − sin η sin η 0 This means that when η = 0 or π/2 every value of φ ∈ [0, 2π) gives the same completely orthogonal pair of rotation planes.In other words, φ and η parameterise a 2-sphere (S 2 ), with the north and south pole given by η = 0 and π/2, respectively.Effectively, in going from Eq. (36) to Eq. ( 38), we have just taken the quotient SU (2)/U (1) = S 3 /S 1 = S 2 , which is the celebrated Hopf fibration [68].We can therefore say that the space of all rotation planes of any given left isoclinic rotation M L is topologically equivalent to a 2-sphere (S 2 ).For example, the orbits of the two points x = ±1, y = z = w = 0 under this set of transformations are cos η −e −iφ sin η e iφ sin η cos η i.e. x = ± cos η, y = 0, z = ± cos φ sin η, w = ± sin φ sin η, which are the north and south hemispheres of the 2-sphere given by x 2 + z 2 + w 2 = 1, y = 0. Now that we have finished this derivation, we will conclude this section with a few observations.Firstly, there is a much quicker way of deriving Eq. (35), that the transformations from a fixed pair of rotation planes of M L to all other rotation planes are represented in C 2 by the unitary group, based on the following argument.Let the rotation which takes one rotation plane of M L into another be given by U .The fact that M L acts on all its rotation planes in the same way means that M L should be invariant under a change of basis by the rotation U .In other words, we have the equation U −1 M L U = M L , which means that M L and U commute.In the complex (C 2 ) representation M L is simply e iα multiplied by the 2 × 2 identity matrix [c.f.Eq (27)], and so in this representation, U must commute with i.The most general way for U to satisfy this is if U is simply a complex matrix.However, we know that U is also a rotation, and we derived earlier that if a 4D rotation (or more generally an orthogonal transformation) in the C 2 representation is given by a matrix, then that matrix is unitary.
Secondly, note that all of the arguments of this section can also be applied to a given right isoclinic rotation M R ≡ M (α, −α), provided we use a different complex representation of R 4 where the position vector is given by (x + iy, z − iw) T = (r 1 e iθ1 , r 2 e −iθ2 ) T .Transformations between rotation planes of M R -which form the left isoclinic subgroup of SO(4) -then take the same form as Eq. ( 36), but with e iθ2 and e iθ ′ 2 replaced by their complex conjugates.Then, the left isoclinic rotations of the form M L ≡ M (α, α) can be factored out, just as rotations of the form of M R could be factored out of the Eq (36).Thus we can obtain an equivalent expression to Eq. ( 38) for the general coordinate transformation from a fixed pair of rotation planes of M R to all other rotation planes, with all redundant parameters removed.It is explicitly given by Again this equation is identical to Eq. ( 38) except that the second element of each of the position vectors is replaced by its complex conjugate.
Finally, consider a fluid undergoing constant rigid rotation associated with a left isoclinic rotation in time, i.e. taking M L (t) ≡ M (α(t), α(t)) with α(t) = ωt, where ω is a constant frequency.Such a system can be described using any completely orthogonal pair of the rotation planes of M L (t); to see this, we define one such pair by r j e iθj = 0 and another by r ′ j e iθ ′ j = 0, with j = 1, 2, where the primed and unprimed coordinates are related by Eq. (38).Then we take the R 4 gradient of Eq. ( 38) as follows and then take the complex inner product of this equation with Eq. ( 38) to obtain Expanding the inner product and evaluating the gradient we end up with a complex equation with a real part given by and an imaginary part given by where the hats above symbols indicates the unit vectors in those directions.The first equation shows that the R 4 position vector takes the same form in both bases, as expected.The second equation is less trivial, and can be physically interpreted as equating two velocity fields; once Eq. ( 44) is multiplied by the frequency ω, the RHS is a velocity field describing rigid left isoclinic rotation through the unprimed planes, while the LHS describes the same thing through the primed planes.That these two are equal shows that either pair can be used to describe such a rigidly rotating fluid, and therefore the fluid exhibits symmetry with respect to all right isoclinic rotations.However, as mentioned above, a superfluid does not behave as a rigid body under rotation, but instead forms quantized vortices [1].Indeed, for all the 4D superfluid vortex states that we shall study in the remainder of this paper, the 4D velocity field is significantly different to that of rigid rotation (Eq.( 44)).Instead, for these states, the SU (2) symmetry, associated with the set of equivalent rotation planes, is naturally broken.This leads to degeneracies between states that are oriented with respect to the different rotation planes of an isoclinic double rotation.We shall see this, first of all, in the next section where we review the case of orthogonal 4D vortex planes, which we previously studied in Ref. [14].

IV. ORTHOGONAL VORTEX PLANES IN 4D SUPERFLUIDS
So far, we have reviewed the well-known physics of vortices in 2D and 3D, and introduced the different types of rotation that become possible in 4D systems.We shall now combine these ideas in order to discuss some of the vortex structures that can emerge in a 4D superfluid under double rotation.In this section, we shall focus, in particular, on the case of orthogonal vortex planes, which we earlier described in Ref. [14].After briefly reviewing the main findings of this previous work, we shall proceed to re-derive the hydrodynamic energy of two completely orthogonal vortex planes within a hyperspherical system, and to introduce our numerical methods, illustrating these by presenting numerical results that complement those already published.The intention of this section is to establish a basis of comparison for when we extend our discussion to non-orthogonal 4D vortex planes in Section V.

IV.1. Structure of 4D Orthogonal Vortex Planes
As in Section II, we want to consider a superfluid described by the GPE without external potentials, but now with atoms free to move in four spatial dimensions.In the absence of rotation, a generalised 4D GPE can be written in the same form as in lower dimensions [Eq.( 1)], namely [14]: except now with ∇ 2 corresponding to the 4D Laplacian.This serves as a minimal model in which to explore 4D vortex physics, and is a plausible mathematical description of low-temperature interacting bosons in a hypothetical universe with four spatial dimensions [14,[69][70][71].In the future, it will also be interesting to consider a more tailored model moving towards a realistic 4D experiment, based, for example, on adding one or more synthetic dimensions to an ultracold bosonic gas [17, 31-41, 55, 56].However, the form of such a model will depend strongly on the details of the particular experimental implementation chosen, and will likely include other effects, such as lattices, unusual interaction terms and asymmetries between real and synthetic dimensions.These more experimental models therefore go beyond our current work, but raise interesting opportunities for future research as discussed briefly in our conclusions in Section VII.While the above 4D GPE is identical to that in lower dimensions, in the rotating frame [c.f.Eq. ( 14)], we have to be more careful as the angular momentum operator can no longer be treated as a vector as discussed in Section III.Instead, in 4D, the angular momentum operator is a 4x4 antisymmetric tensor, with components Lγδ that correspond to the angular momentum in the γ-δ plane (with γ, δ ∈ {x, y, z, w}).The general form of the rotating-frame GPE then takes the form where Ω γδ is the rotation frequency associated with the γ-δ plane, and the sum runs over the six different Cartesian planes in 4D.
As can be seen from Eq. ( 46), the simplest situation is when there is only one plane with a non-zero rotation frequency, e.g.Ω xy ≡ Ω ̸ = 0.This corresponds to the case of simple rotation, which can be understood as a usual three-dimensional rotation extended into 4D [c.f.Eq. ( 22)].As we previously showed in Ref. [14], this sort of rotation can stabilise a single "vortex plane", where the dimensionless order parameter can be described by: where (r 1 , θ 1 ) are plane polar coordinates in the plane of rotation (e.g.x-y), and f k (r) is independent of the coordinates not involved in rotation (e.g.z, w) such that the radial function is the same as that found numerically from the 2D GPE [c.f.Section II].In this case, the vortex core is a single plane (defined e.g. by x = 0, y = 0), as was also verified numerically in Ref. [14].Physically, this can be understood as the natural extension of point vortices from 2D and line vortices from 3D into 4D, as the extra dimension plays no role.
In contrast, double rotations are an intrinsically 4D (or higher) phenomenon and so can lead to much richer vortex physics, as will be our focus in the remainder of this article.Specifically, we will focus on the 4D GPE in a doubly rotating frame where Ω j and Lj are respectively the rotation frequency and angular momentum operator in plane j.For example, we could choose plane 1 as the x-y plane (i.e.Ω 1 ≡ Ω xy , L1 ≡ Lxy = −iℏ(x∂ y − y∂ x )), and plane 2 as the z-w plane (i.e.Ω 2 ≡ Ω zw , L2 ≡ Lzw = −iℏ(z∂ w −w∂ z )).Note that such a set-up is related to certain 4D quantum Hall models in which a nontrivial second Chern number is generated by applying magnetic fields in two completely orthogonal planes [17][18][19][20]72].
To proceed, for simplicity we will henceforward adopt double polar coordinates (r 1 , θ 1 , r 2 , θ 2 ), defined by (x, y, z, w) = (r 1 cos θ 1 , r 1 sin θ 1 , r 2 cos θ 2 , r 2 sin θ 2 ), such that Lj = −iℏ∂ θj .As L1 and L2 describe a double rotation, they commute with each other, meaning that we are able to look for simultaneous eigenstates of both angular momentum operators.As we showed in Ref. [14], for suitable equal-frequency (Ω 1 = Ω 2 ) rotations, a reasonable ansatz for the (dimensionless) ground-state is where k 1 and k 2 are the integer winding numbers in the respective rotation planes and f k1,k2 (r 1 , r 1 ) describes the 4D superfluid density profile, which is assumed to just be a function of the radii of the two planes.This ansatz describes a pair of completely orthogonal vortex planes which intersect at the origin, with the superfluid circulating simultaneously and independently in the two rotation planes with a velocity field given by corresponding to a superposition of 2D vortex-velocity fields in each rotation plane [c.f.Eq. ( 5)].Such a vortex structure is therefore topologically characterised by the Z × Z topological winding numbers [14].Note that this ansatz preferentially picks out the x-y and z-w planes; however, equal frequency rotations are isoclinic and hence have an infinite number of rotation planes, as we discussed in Sec.III.2.This means that suitable ansatzes could be defined with respect to any of these planes, and our choice is arbitrary.The function f k1,k2 (r 1 , r 1 ) in Eq. ( 49) can be found numerically from solving the 4D GPE [Eq.(48)].As we previously showed [14], this function appears to be close to a product ansatz f k1,k2 (r 1 , r 1 ) ≈ f k1 (r 1 )f k2 (r 2 ) where f ki (r i ) is the 2D density profile of a vortex with winding number k i in plane i; however, this separable approximation fails significantly near the intersection of the vortex planes near the origin due to the intrinsic nonlinearity of the GPE equation.Before presenting a numerical example of such a vortex structure, we shall first discuss the associated hydrodynamic energy and critical frequency, in an extension of the standard textbook discussion for 2D vortices that was presented in Section II.

IV.2. Hydrodynamic Energy of Completely Orthogonal Vortex Planes
As it will be helpful in the following sections, we shall here derive the hydrodynamic energy for a pair of completely orthogonal vortex planes in a 4D superfluid.As discussed in Ref. [14], we have previously studied the energy of the orthogonal-vortex structure [Eq.(49)] in a "duocylinder" geometry defined by hard-wall boundaries at r 1 = R 1 and r 2 = R 2 , where R j are the radii in the j = 1, 2 planes.In this geometry, the energy could be approximated by a decomposition as where E kj (R j ) is the 2D energy [Eq.( 13)] associated with having a vortex with winding number k j in a 2D disc of radius R j and hard-wall boundary conditions.In this paper, we shall focus on a 4D hypersphere (or "4D ball") geometry, which is defined by having hardwall boundaries at r 2 1 + r 2 2 = R 2 , where R is the hyperspherical radius.This geometry is theoretically interesting in the following sections as it preserves the symmetry of isoclinic rotations [c.f.Section III] unlike the duocylinder geometry, which has boundary conditions that preferentially pick out two planes as being special.
Nevertheless, as we shall now show, we can also approximate the energy of completely orthogonal vortex planes in such a hypersphere as a decomposition into a sum of the energies of each individual plane, in an analogous manner to Eq. ( 51).To see this, we consider the hydrodynamics of a simplified "hollow core" vortex model, similar to that used in 2D as reviewed in Section II.Specifically, we will consider a pair of completely orthogonal vortex planes, which intersect at the origin, and we will approximate the density profile as zero within one healing length of each vortex core and at the system boundary, and equal to a constant N/V everywhere else, where N is the particle number and V = π 2 R 4 /2 is the volume of the 4D ball.Following Section II, we can again neglect contributions to the hydrodynamic energy from density variations [c.f.Eq. ( 10)] and from interparticle interactions [c.f.Eq. ( 11)], leaving only where in the second equality, we have used that the velocity field is well-described by a sum of the individual velocities for each vortex plane [c.f.Eq. ( 50)].Furthermore, as v j lies in plane j, it follows that v 1 • v 2 = 0, i.e. that the hydrodynamic vortex-vortex interaction term vanishes, leaving us with where we have used that 2 ) from Eq. ( 50).
The integration region is symmetric with respect to swapping r 1 and r 2 , so each term in the integrand gives the same result, just with a different coefficient k 2 j .This means we really only need to consider one of these terms, and with the order of integration we have above, it is easiest to compute the 1/r 2 1 term.We will also rescale r j → Rr j and evaluate the θ integrals to give which to leading order in ξ/R gives us This corresponds to a sum like that in Eq. ( 51), but for the hydrodynamic energy in this simplified constantdensity model in a 4D hypersphere.Note that this equation is very similar to the point vortex energy in 2D [Eq.( 12)], except we have k 2 1 + k 2 2 instead of k 2 , and there is a geometric factor of 2 coming from the difference between the area of a disk and the 4D volume of a hypersphere.Again, just as in the 2D case, we can obtain a more accurate energy for these orthogonal intersecting vortices by using the dimensionless numerical densityprofile function f k1,k2 (r 1 , r 2 ) from our ansatz [Eq.(49)].Using the grand canonical energy relative to the uniform state with a chemical potential given by µ = ℏ 2 /mξ 2 , we obtain numerically via a fitting procedure to the form of Eq. ( 57), with the logarithmic prefactor as the fit parameter.This result, without the initial factor of 2, was obtained the exact same way in our previous paper [14] for a duocylinder geometry given by r 1 ≤ R 1 , and r 2 ≤ R 2 .In that instance, however, the analytical calculation yielded precisely this form since the boundary conditions in the two planes were decoupled, and the energy integral therefore decomposed into a sum in the two planes.In the hyperspherical geomtery, there are less than leading order terms that we have ignored that do not have the same form.Therefore we believe that the result for a spherical geometry [Eq.(58)] is more approximate.The energy reduction from equal-frequency double rotation (Ω 1 = Ω 2 = Ω) can be calculated as where we have used that each vortex plane independently contributes angular momentum equal to N ℏk j [c.f.Eq.( 15)].This gives a critical frequency of to stabilise an orthogonal pair of k = 1 vortex planes.Naturally, we will present most of our results with frequencies in units of Ω c .Note that Ref. [14] worked in units of Ω 2D c , the 2D critical frequency [Eq.( 17)] for a k = 1 point vortex in a disk of radius R. To convert between these unit conventions, we may use the fact Ω c = 2Ω 2D c .

IV.3. Numerical Methods and Results
In this section we will briefly describe the numerical methods which are then used to support our analytical results throughout the rest of this paper.We shall then illustrate these methods with an example of a structure with orthogonal vortex planes, so as to complement the results previously presented in Ref. [14] and to provide a basis for comparison with later sections of this paper.
As in our earlier work [14], the imaginary time evolution method (ITEM) is used to find solutions of the 4D GPE with double rotation [Eq.(46)].We use second order finite differences in space and a first order explicit FIG. 1. Vortex core of orthogonal intersecting planes in the final state of the ITEM with parameters Ωxy = Ωzw = 1.5Ωc, and ∆x = 0.2ξ, giving R ≈ 8.3ξ.Note that the coordinates we are plotting against are rotated relative to those used in the numerics, in order to show both planes at once.discretisation in time.All calculations are performed on a Cartesian grid within a 4D hypersphere of radius roughly equal to N grid ≈ 41 gridpoints, with a hardwall boundary condition imposed on the boundary points (defined as the points with fewer than 8 nearest neighbours).This then corresponds to a total number of gridpoints roughly equal to 1.4 × 10 7 .The spatial step size for most calculations is set to ∆x = 0.5ξ, which ensures a large system of radius R ≈ 21ξ to reduce the importance of boundary effects.
We calculate the predicted critical frequency Ω c from Eq. ( 60) with R set to N grid ∆x − ξ, subtracting one healing length in order to approximately account for the boundary region.As we shall see later, our numerical results suggest that, for ∆x = 0.5ξ, a more accurate value for the critical frequency of 0.9Ω c .This is likely due to a combination of finite size effects and the fact that Eq. ( 58) is an approximate result based on fitting to the function form of Eq. (57).
Our initial states are constructed in terms of a density profile and phase profile, with a degree of noise (up to 20% of the background value) then added to the real and imaginary parts of ψ.The initial density profile is chosen to be homogeneous except at the boundary where it smoothly goes to zero, while the initial phase factor is determined by the vortex configuration we expect to see at low energy for the chosen parameters.The ITEM is deemed to have converged once the relative variations in the particle number, N (calculated as the sum of |ψ| 2 ∆x 4 ), and chemical potential (the sum of the LHS of Eq. ( 46) multiplied by ψ * ∆x 4 /N ), from one iteration to the next reach below 10 −10 .Once the ITEM reaches this threshold accuracy level we output the state and calculate the energy (using a finite difference version of Eqs ( 8) and ( 9)).We also determine the coordinates of all points making up the vortex core and separately output these, where we deem a point to be in the core if it is more than one healing length from the boundary and if the modulus of the order parameter at that location is less than the spatial resolution ∆x/ξ.This latter criterion is motivated by the fact that the order parameter goes to zero linearly as one approaches a singly charged vortex core [1].In order to then plot the vortex core, we supplement a 3D scatter plot (showing the x, y, and z coordinates) with colour (representing the w coordinate).
To illustrate this numerical method, in Fig 1 we show the vortex core structure obtained from ITEM under equal-frequency double rotation with parameters Ω xy = Ω zw = 1.5Ω c , and ∆x = 0.2ξ, giving R ≈ 8.3ξ.Here, the initial state (before adding noise) was chosen to have a density profile that was homogeneous within the system away from the boundaries, and a phase profile given by arctan 2(y, x) + arctan 2(w, z), corresponding to the phase-winding expected from Eq. ( 49) with k 1 = k 2 = 1.Note that both the Cartesian grid and the initial phase profile numerically break the hyperspherical symmetry, and the symmetry associated with isoclinic rotation.However, we expect the grid effects to be small, and the phase profile simply picks out one of several degenerate states for equal frequency rotation.As in Ref. [14], the resulting stationary state is found to contain a vortex core structure consisting of a pair of intersecting completely orthogonal planes (here, corresponding to the x-y and z-w planes).Note that in plotting Fig 1, we have rotated our coordinates according to in order to better depict both planes at the same time.This visualization of the vortex-core structure complements the results previously presented in Ref. [14] for the phase and density profile of the stationary state, and will serve as a useful basis for comparison to results obtained in later sections of this paper under other parameters and initial conditions.

V. NON-ORTHOGONAL VORTEX PLANES
Having introduced vortex structures composed of completely orthogonal vortex planes in the previous section, we shall now consider the possibility of a pair of nonorthogonal vortex planes in a 4D superfluid.As we shall show later in this paper, such non-orthogonal vortex planes are a natural candidate for the low energy configuration of a 4D superfluid doubly rotating at unequal frequencies.In preparation, we shall therefore derive in this section the total hydrodynamic energy of a pair of non-orthogonal vortex planes.
As in Section II and Section IV, we shall neglect the contributions from density variations and from interparticle interactions, such that the kinetic energy can be approximated as where we have again assumed that the velocity field can be decomposed as a sum of the velocity fields associated with each of the (now non-orthogonal) vortex planes separately.From this it can be seen that in general, we can split the total hydrodynamic energy up into a sum of the energies of each individual vortex plane [c.f.Section IV.2] together with the hydrodynamic vortex-vortex interaction, which is given by As we discussed above and in Ref. [14], the velocity fields v 1,2 induced by two orthogonal planes were themselves everywhere orthogonal, v 1 •v 2 = 0, meaning that this hydrodynamic vortex-vortex interaction between the planes vanished.However, as we shall now show, this is not true for non-orthogonal vortex planes, meaning that the vortex-vortex interaction term is non-zero in general.
In order to find E vv , we will start from the assumption that a pair of non-orthogonal vortex planes can be described by the dimensionless ansatz where the primed and unprimed coordinates are related by a double rotation given by a matrix M defined below, such that r ′ = M r, and k 1,2 are the winding numbers of the two vortices while σ j = sign(k j ).The function g ensures that the density returns to the homogeneous value when both r 1 and r ′ 2 are large compared to the healing length, and is given by g , where f k1,k2 is the dimensionless profile associated with the ansatz for orthogonal vortex planes in Eq. (49).In particular, g is always positive, and must satisfy the asymptotic relations where f k is the dimensionless 2D vortex profile described in Sec II.1.These asymptotics physically are the requirement that far from one of the vortex planes, the density profile is determined purely by the remaining one.
Concretely, the state in Eq. ( 64) contains vortex planes along x = y = 0 and z ′ = w ′ = 0 respectively, which intersect at the origin.In order to keep this ansatz general while minimizing the number of parameters, we will refer to appendix A, where we derive the general form of M required to describe the tilting of a plane in R 4 .The result is that, without loss of generality, we may choose the following form with α 1,2 ∈ [0, π/2), such that To use this result in describing our skewed vortex planes we must assume that the vortices exist within a spherically symmetric 4D superfluid of radius R, such that Note that having a pair of orthogonal vortex planes, as discussed in the previous section, corresponds to taking α 1 = α 2 = 0, such that M becomes an identity matrix and r ′ = r.Given the spherical geometry we will assume, in analogy with Section IV, that the velocity fields induced by each vortex have the following simple forms Let's first consider the special case where the matrix M is a simple rotation, meaning one of the angles α 1,2 is equal to zero.In this case it is easier to use the Cartesian representation of Eq. ( 69), which is Without loss of generality we can choose α 2 = 0 such that we have w ′ = w and ŵ′ = ŵ, which is of course orthogonal to both x and ŷ.Therefore, the dot product between the velocity fields is given by Using Eq. ( 68) we have that ẑ′ = sin α 1 x + cos α 1 ẑ, and therefore the interaction energy is given by Keeping in mind that only x and z appear in z ′ , we can see that the above integrand is an odd function of both y and w.This integral therefore vanishes, since our chosen geometry is symmetric with respect to both of these coordinates.In order to get a non-zero interaction potential we would need the superfluid to occupy a region that is asymmetric in both the y and w directions.This may be an interesting avenue for future work but is beyond the scope of this paper.
We shall now derive the form of the hydrodynamic interaction energy [Eq.(63)] in the special cases where M is an isoclinic rotation.There are two main reasons for this choice of rotation: firstly, an isoclinic tilt allows us to derive an analytic form for the interaction using the integral transform into non-orthogonal double polar coordinates derived in Appendix B; and secondly, in Sec VI we will use the results we derive here to investigate possible lowenergy vortex configurations in a superfluid doubly rotating at unequal frequencies, and we obtain predictions that agree closely with numerics when the frequencies are not too high.Isoclinic rotation here corresponds to the condition α 2 = να 1 , with ν = ±1 denoting whether M is left (−) or right (+) isoclinic.Let us therefore define η ≡ α 1 for simplicity and proceed.Using Eq. ( 68), we see that the primed coordinates now take the form where we have applied the shorthand c = cos η, s = sin η.
The equations of the plane z ′ = w ′ = 0 then become This is tilted away from the plane z = w = 0 by an angle η, so the angular separation of the planes z ′ = w ′ = 0, and x = y = 0 is π/2 − η.
Alternatively, we can derive the angle between the planes as follows.First, note that the plane z ′ = w ′ = 0 is spanned by the unit vectors x′ and ŷ′ , and that these are related to the unprimed basis vectors by This lets us define arbitrary unit vectors in the x ′ -y ′ plane and z-w plane as follows where ϕ ′ and ϕ are arbitrary angles between 0 and 2π.We can then find the angle between these vectors in the usual way using the dot product When ϕ ′ −ϕ = π/2 we obtain the minimum angle between the two planes, which is given by arccos(sin η) = π/2 − η.In the rest of the paper, we will sometimes refer to the angle η as the "skewness" between two planes, since it measures how far from orthogonal those two planes are; when η = 0, the two planes are orthogonal, and when η reaches its maximum value of π/2 the two planes coincide.
Continuing, we substitute Eq. ( 73) into Eq.( 64) and find that our ansatz is given by where we have suppressed the arguments of g for brevity.
Note that if we have ν = σ 1 σ 2 , then x+σ 1 iy = x+νσ 2 iy, and the planes are skewed in such a way that they are beginning to perfectly align, while ν = −σ 1 σ 2 corresponds to pure anti-aligning.We therefore expect that ν = sign(k 1 k 2 ) will give rise to a repulsive interaction between the planes, while ν = − sign(k 1 k 2 ) will lead to a attractive interaction [c.f.Section II.3].
In order to compute the hydrodynamic vortex-vortex interaction energy, we will first rewrite Eq. ( 73) in double polar coordinates as As the hydrodynamic interaction energy density depends on ρv 1 • v ′ 2 [c.f.Eq. ( 63)], we must find an expression for the dot product θ1 • θ′ 2 under the assumption of velocity fields of the form in Eq. ( 69).To do this, we will start by taking the vector gradient of Eq. ( 82) as where we have used the primed coordinate system on the LHS and the unprimed coordinate system on the RHS.Then taking the dot product of both sides with θ1 gives Dividing through by ie iθ ′ 2 and then taking the real part of both sides gives , and the vortex-vortex interaction is given by where B d (R) denotes the d-dimensional ball of radius R centred at the origin, which is our chosen geometry.This integral will be as difficult to compute in the primed coordinate system as the unprimed one; however, we can greatly simplify the integrand by using the non orthogonal coordinate system defined by (r 1 , θ 1 , r ′ 2 , θ ′ 2 ), at the cost of complicating the integration limits.In Appendix B we derive the integral transformation into this non-orthogonal coordinate system for the region B 4 (1).To use this result we must make the substitution r → Rr so that the integral is over B 4 (1) rather than B 4 (R), and we make the substitution θ ′ 2 = θ + νθ 1 to simplify the cosine.Altogether, we then have FIG. 2. A geometric/algebraic description of the two terms in Eq. ( 88).(Left) The first term corresponds to the region where sr1 < cr2, meaning that the complex coordinate r ′ 2 e iθ ′ 2 encircles the origin.This means that r ′ 2 always reaches a minimum value of 0, and θ ′ 2 spans a full period, i.e. r ′ 2 ∈ [0, r+], The second term corresponds to sr1 ≤ cr2, which means that r ′ 2 e iθ ′ 2 no longer winds around the origin, so r ′ 2 has some minimum value r− which can be positive, and θ ′ 2 no longer spans a full period.This in turn means that θ ′ 2 takes values in the interval [θ1 − θ * , θ1 + θ * ], where θ * ≤ π/2 with equality occuring when sr1 = cr2 and the blue dotted circle passes through the origin.The limits on r1 are derived in Appendix B, and come from combining the inequalities on sr1 and cr2 with the spherical geometry where the prefactor A = νk 1 k 2 R 2 ℏ 2 /m, and the limits θ * and r ± are given by While this integral transform is derived in great detail in Appendix B, a quick pictoral explanation for the form taken by Eq. ( 88) is given in Fig 2 .In particular, the fact that there are two distinct terms in this equation with different integration limits is directly related to whether r ′ 2 e iθ ′ 2 encircles the origin (first term) or not (second term).Now that we are using the natural coordinates for this problem, we can proceed to evaluate the integrals.Similar to Section II.1, we will ignore the vortex core, and approximate the density as constant ρ = N/V , where N is the particle number and V = π 2 R 4 /2 is the 4D volume of the system.This approximation works provided the angle η is not too close to π/2 as will be discussed later.The vortex-vortex interaction energy is then given as where the new prefactor A ′ = 4νk 1 k 2 N ℏ 2 /πmR 2 .After further algebraic steps detailed in Appendix B, we then obtain the final result We now see that our expectation regarding the sign of E vv was correct: the overall sign is given by sign(νk 1 k 2 ) such that the interaction is positive (i.e.repulsive) when the planes are skewed in an aligning sense (ν = sign(k 1 k 2 )), and negative (i.e.attractive) when they are anti-aligning (ν = − sign(k 1 k 2 )).
Combining this with the result from Sec IV.2, we have that the total hydrodynamic energy of the nonorthogonal vortex plane state is Note, however, that we just derived the interaction term under the constant density approximation, while the individual hydrodynamic energies of the vortices was calculated using a hollow core model.This hollow core was needed to remove the unphysical 1/r 2 singularity around each vortex that gives a divergent contribution to the energy, which in a mathematical sense is why vortices have cores.In contrast, the interaction energy density only goes as 1/r 1 r ′ 2 which is not singular once integrated.Recall, as per the discussion in Sec II.3, that the same is true of point vortices in 2D: their interaction can be approximated with a constant density, but this fails to give a finite answer for the individual energies.As previously discussed, the correct answer can still be obtained if we take the vortices to be combined once their separation is of the order of ξ or below.
The question arises whether we can recover the expression for vortex combination in this 4D case.Here we have an angle between the vortex planes, given by π/2 − η, and we see that the interaction energy E vv diverges in the limit η → π/2 under the constant density approximation.This also occurred when using this approximation for point vortices in 2D as their separation distance ∆r approached zero, so we see that the angular separation π/2 − η is playing a similar role here in 4D as ∆r did in 2D.In contrast to 2D, however, we do not have a unique value for the separation distance between the vortices, so coming up with a criterion for when they have combined seems difficult.
We will identify a natural separation as follows: each plane makes a circle of intersection with the boundary of the hypersphere, and we argue that the maximum reasonable value for the distance between the planes should be given by the minimum distance between these circles.This distance is the length of the most direct straight line between the two planes at the boundary, and -as we derive in Appendix C -is given by √ 2R (1 − sin η) 1/2 .This result can be obtained by naively applying the cosine rule in an analogy to lines in 2D.Setting this distance less than or equal to ξ, we have which rearranges to We want this inequality in terms of cos η, since this is what appears in the interaction.Therefore we square both sides, which is safe as they are each non-negative, giving which rearranges to To leading order in ξ/R, we can therefore say that the vortices are combined once cos η = ξ/R or less.Substituting this into Eq.( 94) (using that ν 2 = 1) we obtain which is the correct result for a vortex plane with winding number k 1 + νk 2 .This concludes our discussion of the energetic properties of non-orthogonal vortex planes, the results of which we will shall now use to construct a model of superfluids doubly rotating at unequal frequencies.

VI. UNEQUAL FREQUENCY DOUBLE-ROTATION
In this section, we will consider the behaviour of a 4D superfluid undergoing constant double rotation with unequal frequencies, given by Ω xy = Ω zw + ∆Ω in the lab (x, y, z, w) frame.We will take ∆Ω > 0 without loss of generality, and we will also assume the superfluid occupies a hyperspherical (4D ball) region for simplicity.Note that this breaks the isoclinic SU(2) symmetry [c.f.Sec.III], which is associated with the manifold of different rotation planes when the frequencies are equal.We will begin in Section VI.1 by considering what happens to a single vortex plane in this set-up, before discussing the case of two vortex planes in Section VI.2.We shall then present our numerical results in Section VI.3.

VI.1. A single vortex plane under unequal frequency double rotation
In this section, we shall develop our intuition by considering the simple case of a single vortex plane in a system with unequal frequency double rotation.We will assume that this plane remains rigid but allow it to arbitrarily tilt in order to optimise its energy.The energy of the superfluid in the rotating frame is reduced by the amount relative to an inertial frame.Since Ω xy > Ω zw it is natural to presume that the lowest energy occurs when the vortex plane lies along the z-w plane -thereby inducing rotation in the x-y plane -such that ⟨ Lxy ⟩ = N ℏ and ⟨ Lzw ⟩ = 0.The converse case of a vortex plane spanning the x-y plane would certainly be higher in energy, as this state would have ⟨ Lxy ⟩ = 0 and ⟨ Lzw ⟩ = N ℏ; however, it is not obvious that the superfluid energy decreases monotonically as the vortex plane is tilted from the x-y plane to the z-w plane.That is, the lowest energy overall could occur when the vortex plane is oriented somewhere between these two limits.Such a state would have positive values of both ⟨ Lxy ⟩ and ⟨ Lzw ⟩, and would be given as follows where the primed coordinates are to be defined shortly, and the function g is given by g(r 2 ) = const × f 1 (r/ξ)/r, with f 1 the dimensionless density profile of a point vortex in 2D (see Sec II.1).The primed coordinates are defined such that the vortex plane is given by x ′ = y ′ = 0, soas derived in appendix A -we may assume without loss of generality that the primed coordinates are given by A more useful form of the order parameter for this state is given in polar coordinates, r ′ 1 e iθ ′ 1 ≡ x ′ + iy ′ , as follows In order to use this, we must express the transformation into the primed coordinate system in polar coordinates as well.Taking the combination x ′ + iy ′ gives and switching to polar coordinates and rearranging gives The terms on the second line are proportional to e −iθj , and therefore generate angular momentum counter to the external rotation (since Ω xy L xy = −iℏΩ xy ∂ θ1 and similarly for zw).To maximize the energy reduction from rotation we should eliminate these terms.We therefore set α 1 = α 2 ≡ η to obtain the following An expression of this same exact form [Eq. ( 30)] was derived in Sec III.2 in the context of finding the rotation planes of any left isoclinic rotation in the x-y and z-w planes, i.e. a rotation M L generated by L+ = Lxy + Lzw .Here, however, r ′ 1 e iθ ′ 1 = 0 encodes the vortex plane, so we see that the vortex plane always lies in a rotation plane of M L , regardless of the value of η.This vortex generates angular momentum in the plane orthogonal to itself, which is also a rotation plane of M L .This is all essentially summed up by the (easily verifiable) fact that r ′ 1 e iθ ′ 1 is an eigenfunction of the sum of the angular momenta, L+ , despite not being an eigenfunction of either component.The superfluid containing this tilted vortex therefore has the same value of ⟨ L+ ⟩ = N ℏ for every value of η.We can exploit this knowledge by rewriting the rotational energy [Eq.( 100)] as Since the first term is constant with respect to η we therefore maximize E rot by simply maximizing the second term, proportional to ⟨ Lxy ⟩.This clearly occurs at η = 0, where r ′ 1 e iθ ′ 1 = r 1 e iθ1 , and so the initial intuition was correct: a single perfectly rigid vortex plane will always want to fully align with the higher frequency.

VI.2. Two vortex planes under unequal frequency double rotation
Armed with this knowledge we now seek to find the optimal configuration of two rigid vortex planes in a doubly rotating superfluid with unequal rotation frequencies.Each vortex will want to align its angular momentum with the x-y plane as much as possible to gain from the larger rotation frequency in this plane.However, as we showed in Sec V, vortex planes will interact with each other hydrodynamically once they are not orthogonal.This interaction will limit how close together (in orientation) each vortex can be and the competition between this effect and the rotational energy will determine the optimal orientation of each plane respectively.This is a simplified model of the situation, and it is worth briefly discussing the approximations we are making.
We are going to again assume a constant density profile given by ρ = N/V , thereby ignoring the vortex core.As previously discussed in Section II.1, one needs to account for the core to avoid a divergent hydrodynamic energy cost of a vortex -see, for example, Eq. ( 10) with ξ taken to zero.However, in this case we are only interested in how the energy varies with the orientation of each vortex plane.This means we can ignore any terms which do not vary as the planes tilt.If we denote the velocity field induced by each plane by v j , with j = 1, 2 respectively, then the hydrodynamic energy can be expanded as in Section V as The first term is the individual hydrodynamic cost of each vortex, which diverges if we ignore the core by assuming a constant density.However, this term does not vary with orientation due to the spherical symmetry of the boundary.On the other hand, the second term, which is the hydrodynamic interaction between the planes, depends on their relative orientation but does not diverge in a constant density approximation, as explained in Sec V.
We therefore ignore the constant first term, keeping only the second term which can safely be approximated using a constant density.This constant density approximation also allows us again to ignore the energy contributions from quantum pressure (the first term in Eq. ( 8)), and the bosonic interaction (Eq.( 9)) as is also done in Section V.
As we are assuming a constant density n = N/V , we may take the order parameter for this configuration to be where the acute (ŕ) and grave (r) coordinate systems are to be defined shortly.Note that this assumes the vortex planes remain flat and intersecting at the origin.
In our numerical results, we do find some curvature and an avoided crossing near the origin (see Section VI.3), but these seem to have only a small effect on the energies.We will investigate the phenomena of curved vortex cores with avoided crossings in more in [15].Recall from Sec VI.1 that the vortex planes will want to stay on one of the rotation planes of left isoclinic rotations, M L , generated by L+ .We will therefore use our result [Eq.(38)] from Sec III.2 for the general form of such rotation planes relative to a fixed basis (in this case the lab basis, x, y, z, w).Using this result for each of the acute and grave coordinates, we have where η 1,2 ∈ [0, π/2] and φ 1,2 ∈ [0, 2π), with φ j undefined when η j = 0 or π/2.The location of each vortex plane is then given by ŕ1 e i θ1 = 0 and r2 e i θ2 = 0, respectively.The parameters η 1,2 denote the angle that each plane makes with the x-y (resp.z-w) plane, while φ 1,2 denote the direction of this tilt.The vortex at ŕ1 e i θ1 = 0 is tilted by an angle η 1 away from the x-y plane, while the vortex at r2 e i θ2 = 0 is tilted by π/2 − η 2 off of the same plane.Since the two vortex planes are indistinguishable, we can define the acute and grave coordinates such that the former vortex is closer to the x-y plane than the latter, which translates to the following constraint on the angles Note that η 1 = η 2 = 0 corresponds to the configuration that we have previously studied [14] -a completely orthogonal pair of vortex planes spanning the rotation planes of the superfluid: the x-y and z-w planes, respectively.A change of basis in either of these planes redefines the φ j variables as follows e −iα 0 0 e −iβ cos η j e iφj sin η j −e −iφj sin η j cos η j e iα 0 0 e iβ = cos η j e i(φj −α+β) sin η j −e −i(φj −α+β) sin η j cos η j (114) We will choose a basis in which −φ 2 = φ 1 ≡ φ, leaving us with three free parameters, (η 1 , η 2 , φ), describing the orientation of the two vortex planes relative to the lab frame.In terms of these parameters the two planes are defined by the zeroes of the following complex coordinates ŕ1 e i θ1 = cos η 1 r 1 e iθ1 + e iφ sin η 1 r 2 e iθ2 , r2 e i θ2 = cos η 2 r 2 e iθ2 − e iφ sin η 2 r 1 e iθ1 , restated here for clarity.We will now find and then minimise the sum of the rotational and hydrodynamic energies of the superfluid with respect to these three variables.

VI.2.1. Rotational Energy
Firstly, we will calculate the rotational energy, which is the expectation value of −Ω xy Lxy − Ω zw Lzw in the state ψ [c.f.Eq. ( 100)].Since this is a first order differential operator we may use the product rule on Eq. ( 110), e.g. for each angular momentum component as ψ * Ω xy Lxy ψ = n e −i θ1 Ω xy Lxy e i θ1 + e −i θ2 Ω xy Lxy e i θ2 , (118) such that the rotational energy density is simply the sum of contributions from each vortex independently.From Section VI.1, recall that the favourable possible orientations of the two planes are limited to those which are planes of rotation of the isoclinic rotation generated by L+ .This means that -as in the single vortex case we just considered -we can rewrite the rotational energy density as where the extra factor of 2 in the first term on the RHS arises because we now have two vortex planes instead of one [c.f.Eq. ( 108)].This means that to proceed we simply have to evaluate the angular momentum of each vortex in the x-y plane.In other words, we must compute the following integral Lxy e i θ1 e i θ1 d 4 r; (120) this is carried out in Appendix D, with the final result that The calculation for the e i θ2 term follows identical logic and so we simply state the result, which is Putting these together, remembering that n = N/V and V = π 2 R 4 /2, we have for the rotational energy of two rigid, intersecting vortex planes under unequal frequency double rotation, with η 1 (resp.η 2 ) denoting the angle that the first (resp.second) plane is tilted compared to the x-y (resp.z-w) rotation plane.

VI.2.2. Vortex-Vortex Interaction Energy
Secondly, we consider the hydrodynamic vortex plane interaction previously derived in Sec V, which, as we state again here, is calculated from As this depends on the dot product between the velocity fields of each individual vortex, this energy is entirely dependent on the skewness, η, of the two planes, that measures how far from being mutually orthogonal the vortex planes are (see Sec V for details).The most direct way to find η is to take the acute and grave coordinates and define a rotation transforming between them.Remembering that each of the vortex planes is a rotation plane of a left isoclinic rotation M L , we can use our result from Sec III.2 that the transformation between them has the following form [c.f.Eq. ( 36)] r1 e i θ1 r2 e i θ2 = e iϕ1 cos η e iϕ2 sin η −e −iϕ2 sin η e −iϕ1 cos η where η ∈ [0, π/2] and ϕ 1,2 ∈ [0, 2π).Substituting Eq. ( 112) into Eq.( 111) gives a relation of this form; specifically, examining the top-left entry of the matrix allows us to relate η to η 1 , η 2 , and φ as follows Taking particular values of φ in the above equation will give us intuition for how this parameter corresponds to the direction of the tilt, as mentioned previously, and hence allow us to deduce the form of the vortex-vortex interaction.In particular, we will examine the cases in which φ is a multiple of π/2, as this renders the RHS of Eq. ( 126) real and non-negative, allowing us to simplify this equation by choosing ϕ 1 = 0. Specifically, when φ = 0 or π we have while φ = π/2 or 3π/2 gives Note that cos(η 1 + η 2 ) is non-negative due to the constraint on η 1,2 [Eq.( 113)].
To understand these special cases let's look at the equations of the vortex planes directly.Substituting φ = 0 into Eqs (115) and (116) gives such that the two planes are defined by while the same procedure for φ = π/2 gives ŕ1 e i θ1 = 0 : In the φ = 0 case, Eqs ( 131) and ( 132) describe a pair of lines in the 2D (x, z) subspace, and identical lines found by taking the first lines and sending (x, z) → (y, w).Similarly, when φ = π/2 we have a pair of lines in the (x, w)  131) and ( 132), which describe a pair of lines in x, z space and an identical pair of lines in y, w space, plotted simultaneously in the left panel.When φ = π/2 we have instead Eqs ( 133) and ( 134), describing a pair of lines in x, w space and their reflection about the horizontal axis in y, z space, which we plot on the same graph here by transforming z, w → w, −z.Equivalent pictures for φ = π and 3π/2 respectively can be found by reflection about the vertical axis.
subspace, and a pair in the (y, z) space which are related to the first pair by (x, w) → (y, −z).For φ = π, 3π/2 note that φ → φ + π transforms the equations for the planes simply by z → −z and w → −w.
Since the equations separate into lines in 2D subspaces this way, we can visualise the planes by simply plotting these lines, as shown in Fig 3 .From the left panel of this figure it is clear that when φ = 0 (or π) the vortex planes are tilted away from the coordinate planes in the same direction.Similarly, the right panel shows that when φ = π/2 or 3π/2 the vortices are tilted in opposite directions, and hence towards each other.Other values of φ interpolate between these two scenarios such that the two vortex planes are not tilted along a common direction.This visual understanding also agrees with the two expressions for η, given in Eqs ( 127) and (128).Now we can use physical intuition to deduce -without any further calculation -the most energetically favourable value for φ for any fixed values of the parameters η 1,2 .Recalling that the rotational energy was independent of φ, we need only consider the interaction potential between the two vortices, given by which is positive and therefore repulsive.This result was derived in Sec V in the case that ϕ 1,2 = 0, however these angles do not affect the interaction since they can be absorbed into the definition of θ1,2 and θ1,2 .Since this interaction is repulsive we can clearly see that it is maximised when φ = π/2 or 3π/2 as the vortex planes are tilted directly toward one another.Equally, in the other case where φ = 0 or π the vortices are tilted in the same direction and the interaction energy cost is minimised.Therefore, we can set η = |η 1 −η 2 | and proceed with find-ing the minimum energy as a function of the remaining parameters η 1,2 .For concreteness we will also set φ = 0, but note that φ = π provides an equivalent solution with the same energy.

VI.2.3. Finding the minimum
Our final step is to add the rotational and the vortexvortex interactions energies together and to minimise the resulting sum.Firstly we will define quantities that will make the calculation simpler.Let E ⊥ rot = N ℏ(Ω xy +Ω zw ) denote the reduction in energy due to rotation of the state with orthogonal vortex planes along the x-y and z-w planes [c.f.Eq. ( 59)].We then define a dimensionless energy density relative to E ⊥ rot , given by and a dimensionless frequency difference ω = R 2 ℏ∆Ω/2ξ 2 µ.
Note that in units of the critical frequency Ω c [Eq. ( 60)], this dimensionless frequency is given by ω = ln(2.07R/ξ)∆Ω/Ωc .We then must find the minimum of the following 137) Note that we have not needed to include the absolute value on the RHS of Eq. ( 127), since cos|x| = cos x.Additionally, the logarithmic divergences as η 1 − η 2 → π/2 are unphysical as there the vortex planes coincide and the constant density approximation that we took in Sec V fails.
Setting the derivatives of this energy to zero gives us the following simultaneous equations Firstly, examining the sign of the terms in each of these equations (recalling that η 1,2 < π/2), we must have that Physically, this is because if η 1 is greater than η 2 then the force from the repulsive interaction acts in the same direction as the force from the rotational energy.Therefore we can eliminate the absolute value in Eq. ( 127), such that Secondly, Eqs (138) and (139) together imply sin 2η 1 = sin 2η 2 .There are two ways to satisfy this, i.e. by taking The former case is precisely the condition that the two planes are orthogonal, which eliminates the interaction term.Substituting Eq. (142) into Eqs (138) and (139) then leads to the result η 1 = η 2 = 0, the state we have previously studied [14].This state has an energy of ε = 0 by definition, since ε was defined relative to this state.Moreover, we also note that if we substitute Eq. (142) into Eq.( 137), we see that any state with η 1 = η 2 has energy ε = 0. Interestingly, these orthogonal states are still all degenerate despite the isoclinic symmetry being broken when ω ̸ = 0.The stationary point at (η 1 , η 2 ) = (0, 0) is therefore a saddle point, since it has this line of constant energy passing through it.The latter case is much more interesting as it arises from competition between the interaction and rotational energies.Note that the relation between the tilt angles [Eq.( 143)] ensures that the two planes are symmetrically tilted with respect to the rotation planes of the superfluid; what we mean by this is that each vortex makes the same angle with the x-y plane, and also with the z-w plane.This can be seen by considering Fig 3 and noting that the vortices are each tilted away from the x-y plane by an angle of π/2 − η 1 , and η 2 respectively.When η 1 = π/2 − η 2 these two angles are equal, and the same is of course true with the angles the vortices make with the z-w plane.Using Eqs (141) and ( 143) we can write FIG. 5. Numerical vortex core in the final state of the ITEM-evolved 4D GPE under double rotation [Eq.( 46)].The spatial step size was ∆x = 0.5ξ, corresponding to R ≈ 20.6ξ.Rotation frequencies were Ωzw = 0.85Ωc, and Ωxy ≈ 1.43Ωc, chosen such that the predicted skewness was precisely η = 40 • .The initial phase profile was given by that of the ansatz described in Sec VI.2, with this predicted value of η and added noise.The first two panels show two different views of the core in (x, y, z) space, with points coloured according to their w value.The overall structure resembles the skew planes of the ansatz but the second panel clearly shows how the core curves away from these planes to form an avoided crossing [15].The final panel shows the points projected down into (x, z) space -again with w shown as colour -as well as the lines z = ±x tan(π/4 + η/2), which the theory predicts that the core points should lie along from this perspective.The agreement between these analytical lines and numerical points is very good.
both the angles η 1,2 in terms of the skewness as follows which makes it clear that η 1 ≤ π/4 and η 2 ≥ π/4.At this point it is worth restating our ansatz, since it now only depends on η.Recall that the order parameter is defined as [Eq.( 110)] ψ = n 1/2 e i θ1 e i θ2 , with these angles defined by Eqs ( 115) and (116).Substituting φ = 0 and the above equations for η 1,2 , Eqs ( 115) and (116) become where we have used that cos(π/4±η/2) = sin(π/4∓η/2).Note that we can now clearly see that two planes are arranged symmetrically, in the sense that after a rotation of angle π in the x-y plane (θ 1 → θ 1 + π) the two equations Eq. ( 146) and (147) swap and hence the two vortex planes swap (note that this is also true for a π rotation in the z-w plane, up to a shift in the angles θ1 and θ2 ).This symmetry can also be seen in the the equations for the vortex planes, Eqs (131) and (132) which are now given by where + and − refer to the planes given by r2 e i θ2 = 0 and ŕ1 e i θ1 = 0, respectively.From these equations we can actually see that this configuration is invariant under a π rotation in any one of the six coordinate planes.
We also now see from these equations that both vortices are closer in angle to the z-w plane than they are to the x-y plane.An interesting consequence of this is that when η = 0 the orthogonal state we get does not consist of vortices spanning the x-y and z-w planes.Instead the vortices occupy a pair of diagonal (in terms of the lab frame) planes, given by z = ±x and w = ±y.This doesn't matter when the frequency difference is zero, as then the rotation is isoclinic and these diagonal planes are also rotation planes [c.f.Section III], but for any other value of ∆Ω the only rotation planes are the x-y and z-w planes so it is perhaps surprising that none of these states ever occupy them.
Substituting Eqs (144) and (145) into Eq.( 138) gives the following relation between ω and the optimal skewness which rearranges to the following quadratic equation for sin η This has only one solution for sin η in the interval [−1, 1], given by Physically this means that the optimal skewness vanishes in the limit that ω (i.e. the frequency difference in rotation) goes to zero, corresponding to the situation where the two vortex planes become completely orthogonal, as expected.In the opposite limit that ω becomes very large, this formula instead predicts that sin η → 1 and hence η → π/2, meaning that the angle between the two planes goes to zero.This corresponds physically to the two planes both aligning with the z−w plane so as to maximise the energetic reduction due to the higher rotation frequency Ω xy > Ω zw .However, this limit should also be treated with caution, as at high enough frequencies, we expect that it will become energetically favourable to introduce more vortices and/or more complicated vortex structures, as briefly discussed in Appendix E. We also expect there will be other contributions to the energy, which we have neglected here; for example, our assumption of a constant density profile will break down when the two vortex planes become very close together.Now we can find the energy of this optimally skewed state; using Eqs (141), (143), and (145) the energy [Eq.( 137)] becomes Rearranging Eq. ( 149) we can quickly find that cos 2 η = 2 sin η/ω, and after using cos(π/2 + η) = − sin η we have everything in terms of sin η.Substituting Eq. ( 151) then gives the energy density for the optimal skewed statesin terms of ω as For small and large ω we have the following asymptotics Recall that ω ∝ R 2 ∆Ω, therefore these limits can be reached by decreasing (resp.increasing) either ∆Ω or the radius R.
This energy [Eq.( 154)] is also negative for all ω > 0, which means our simplified model has predicted that this tilted vortex plane state is lower energy than the orthogonal state for any frequency difference ∆Ω > 0. Fig 4 shows the energy landscape as a function of both tilt angles for dimensionless frequency difference ω = 2.As expected, we see a line of constant energy along η 1 = η 2 , a minimum energy along the line η 1 = π/2 − η 2 , and a range of tilt angles for which the energy is negative.

VI.3. Numerical Results
We will now compare the above analytical predictions for tilted vortex planes to numerical results obtained using the methods described in Sec IV.3.We choose an initial phase profile identical to that of our non-orthogonal vortex ansatz [Eq.(110)], where the acute and grave coordinates are given by Eqs ( 146) and (147), respectively, with a chosen value for the skewness, η.We then use Eq. ( 149) to calculate the frequency difference, ∆Ω, that will energetically favour the chosen value of η if the model is accurate, and then we run the ITEM with this value of ∆Ω on our initial state with added noise.
Once the ITEM is converged we compare the geometry of the vortex core in the numerical final state with that of our predictions.Fig 5 shows the numerical vortex core for ∆x = 0.5ξ, which corresponds to a system radius of R ≈ 20.6ξ, and with frequencies Ω zw = 0.85Ω c , and Ω xy ≈ 1.43Ω c , corresponding to a predicted skewness of η = 40 • .The first two panels show two different rotations of the core in (x, y, z) space, with the points coloured according to their w value (see the colourbar on the far right).Already we can see that, at large distances from the origin, the vortex cores look like the predicted tilted planes, symmetrically arranged with respect to the rotation planes of the superfluid.The third panel in Fig 5 shows a side-on view where the vortex core appears approximately as a pair of lines, just as in Fig 3 .On top of these data points we have plotted the lines z = ± tan(π/4 + η/2)x [c.f.Eq. ( 148)] which are the predicted lines on which the numerical data should lie.As can be seen in the figure, there is excellent agreement between these numerical final states and our analytic predictions.However, note that near the origin we seemost prominently in the second panel of the figure -an avoided crossing structure [15], as also discussed further below.
In addition to the above qualitative comparisons of the numerical and predicted vortex cores, we have made a quantitative analysis of the accuracy of our predicted energy [Eq.( 154)].To do this we performed the ITEM for a range of different frequencies Ω zw and ∆Ω, again using our ansatz to determine the initial phase, and then calculating the energy of each final state.The procedure for these calculations was as follows; we fixed a value of Ω zw , then ran the ITEM with ∆Ω = Ω c on our prescribed initial state.From the final state we calculated the energy, and then to speed up calculations we used this final state as the initial state for the next ITEM run with ∆Ω = 0.9Ω c .This process was repeated down to the isoclinic point, ∆Ω = 0. Finally, many of these loops were run at once with different values of Ω zw , so that we could explore an area in frequency space rather than just a line.The results for these energies are shown in Fig 6, with each value of Ω zw corresponding to a different colour.On top of these points we have plotted lines given by performing a single fit of this data over the area in frequency space to a redimensionalized version of Eq. ( 154) given by FIG. 6. Numerical (points) and analytical (solid lines) results for the energy of aligning non-orthogonal vortex states -as a function of Ωzw and ∆Ω -described in this section, as well as numerical results for the energy of the orthogonal vortex state (dotted lines).The resolution was set to ∆x = 0.5ξ, giving a system radius of R ≈ 20.6ξ.All of the analytical lines were generated by a single fit of the numerical points to Eq. ( 157), with E0 and R as fitting parameters.The fit produced E0 ≈ 0.71µN and R ≈ 19.5ξ, and agrees excellently with the numerical data.The numerical zero-vortex ground state in this system has an energy of E hom ≈ 0.67µN , which suggests that the numerical critical frequency in this system is very close to 0.9Ωc, where Ωc is our predicted value [Eq.( 60)].
where E 0 is simply the energy of the system when the vortex planes are orthogonal and there is no external rotation.Since ξ and µ are both external parameters in the numerics, the only parameters in the fit were R and E 0 , the energy at Ω zw = Ω xy = 0. Furthermore, R is not truly a free parameter as it is determined up to boundary effects by the radius of the simulated region, which is roughly 20.6ξ.The fit produced a value of R ≈ 19.5ξ, which is consistent if we estimate the size of the boundary region to be roughly equal to ξ.Additionally, the agreement between the numerical points and lines from the fit is excellent.We also attempted to track the energy of the theoretically predicted skew plane branch at a frequency of Ω zw = 1.3Ω c , by performing two of these iterative ITEM runs from ∆Ω = Ω c down to ∆Ω = 0. Interestingly, as we decreased ∆Ω for these runs the states we obtained increasingly diverged from the skew plane states, even down to the isoclinic point at ∆Ω = 0.For plots of these states at the isoclinic point, see Appendix E. The dotted line shows the numerical energy of the orthogonal state [14] as a function of the frequencies.This was found by running the ITEM once, with Ω zw = Ω xy = 1.2Ω c , and then calculating the energy of this fixed final state for different values of Ω zw and ∆Ω.As shown in Fig 6, the energies of the orthogonal state form dotted straight lines that meet the fit lines tangentially at ∆Ω = 0, which is exactly as predicted since the skewness of the tilted state is approaching zero in this limit.Finally, note that the numerical zero-vortex ground state in this system has an energy of E hom ≈ 0.67µN , which is roughly equal to the energy of the skew and orthogonal vortex states when Ω zw = 0.9Ω c , and ∆Ω = 0.This means that the value of the numerical critical frequency in this system is very close to 0.9Ω c , where Ω c is our predicted value [Eq.( 60)].
To get an idea of the size of the avoided crossing as a function of the frequencies, we have taken each state represented in Fig 6 and calculated the minimum distance between its vortex core and the origin, which we denote r min .This is then plotted in Fig 7, which shows that, in general, the avoided crossing decreases in size with both Ω zw and ∆Ω.Note that these lines are not perfectly smooth and also change their ordering as ∆Ω changes, suggesting that there are multiple metastable branches being sequentially followed by our numerical states.Nevertheless, these must be very close together in energy, since the fit in Fig 6 is very good, and the long range core structure has good agreement with the predicted state.
Minimum distance rmin of the vortex core from the origin -as a function of Ωzw and ∆Ω -for the nonorthogonal numerical states whose energy is shown in Fig 6.
This quantity captures the size of the avoided crossing region seen in Figs 5 and similar final states.As shown above, this region generally shrinks as either frequency increases.These data are not particularly smooth, which may be due to some sampling error from the discretisation of the Cartesian grid.
Lastly, we have further tested our analytical results by using a different initial phase profile in the numerics.We still use the tilted plane ansatz [Eq.( 110)], but with different values of η 1 and η 2 than those predicted.Instead of the predicted values, η 1 + η 2 = π/2 [Eq.( 143)], we chose η 1 = δ, η 2 = η + δ, where η is the skewness of the theoretically predicted configuration, and δ is a small angle added to ensure all symmetries are broken.The corresponding planes still have skewness given by η, but are now asymmetric with respect to the planes of the external rotation.This initial state (with added noise) converges to the same final state as Fig 5 under the same value of all parameters (∆x = 0.5ξ, corresponding to R ≈ 20.6ξ Ω zw = 0.85Ω c , Ω xy ≈ 1.43Ω c , corresponding to η = 40 • ), showing that this predicted state is likely the ground state in this regime.However, at a higher value of Ω zw = 1.25Ω c , we find very different final states of the ITEM depending on which initial phase profile is used.Specifically, using the predicted phase profile (with added noise) we find the same final states as before, with just as good agreement to the analytics.Using the asymmetric phase profile (with added noise) described above, we find very different vortex core structures, with slightly higher energies than the theoretical states (see Appendix E).

VII. CONCLUSIONS
In the first part of this section, we will focus on summarizing the main conclusions of our paper and highlighting the open questions that directly follow on from our study.In the second part of this section, we will then briefly discuss the more general future outlook for research into topological excitations in 4D superfluids, looking beyond the physics of the minimal model studied in this paper.

VII.1. Summary
In this paper, we have demonstrated that stationary states of the 4D GPE under unequal-frequency double rotation can host complicated vortex core structures consisting of skew planes and curved surfaces.This work generalises Ref. [14], which focused on completely orthogonal and rigid vortex planes, and lays the groundwork for Ref. [15], which will explore whether a different configuration of tilted vortex planes can be energetically favoured under equal-frequency double rotation.
In more detail, we firstly showed in Sec.V that nonorthogonal vortex planes interact hydrodynamically, in contrast to the case of orthogonal planes, and we analytically derived the form of this interaction in the special cases of planes related by an isoclinic or simple rotation.Understanding this interaction potential in the general double rotation case would be an interesting topic for future work, allowing our analytics to extend to the full configuration space of non-orthogonal vortex planes.
Secondly, we approached the problem of a 4D hyperspherical superfluid doubly rotating at unequal frequencies, under the assumption that the vortices remain rigid planes.In Sec.VI.1 we showed that a single vortex plane in this system will always want to fully align with the higher of the two rotation frequencies.Then, using this result, we tackled the case of an intersecting vortex pair in Sec.VI.2, proposing a non-orthogonal such pair as an ansatz for the ground state.This was based on the observation that for both planes to benefit from the higher frequency they had to be skewed in a purely aligning sense, thereby inducing repulsive vortex-vortex interaction from our result in Sec.V. We built an analytic model based simply on the balance between these two energies, which predicted that a skew configuration of vortex planes could indeed have lower energy than an orthogonal one.With this model we were able to find which of these configurations was optimal, and calculated the predicted tilt-angles and energy.Comparing these results with numerics, we found excellent agreement, despite the fact that we did not account for the avoided crossing of these states that was seen numerically, but which will be discussed further in Ref. [15].At high frequencies we also found that more exotic states with highly curved vortex surfaces could appear (see Appendix E), suggesting the ground state of a doubly rotating superfluid is in general very complex.
The results presented in this paper show that the physics of vortex surfaces in 4D can be incredibly rich, even in the absence of dynamics.The fact that curved and tilted vortex surfaces can be stable and exist at low energies in such a minimal model is a dramatic departure from the physics of lower dimensions under rotation, suggesting a vast configuration space to explore and investigate in the future.

VII.2. General Outlook
In this paper, we have focused on a minimal theoretical model for a 4D superfluid based on the 4D GPE under rotation.This is motivated as the simplest extension of textbook 2D and 3D superfluids [c.f.Section II] into higher spatial dimensions [14].However, in the future, it will be both interesting and relevant to go beyond this simple model to study more realistic systems with the aim of making an experimental proposal, and to explore the even richer vortex physics that will likely emerge.
However, as we discussed in our previous paper [14], the 4D GPE that we have considered [Eq.(48)] is a toy model lacking elements which are necessary for experimental relevancy to current synthetic dimension approaches.For example, we have considered a purely hypothetical four-dimensional space that is isotropic and continuous, and we have chosen a hyperspherical hardwall boundary to preserve the rotational symmetry.However, the motion, inter-particle interactions, and boundary conditions along any synthetic dimensions can differ from that in real space.In practice, it is likely that an experiment may contain both real and synthetic dimensions, which would break SO(4) rotational symmetry.This will affect the behaviour of the tilted and curved vortex planes that we have studied, adding in additional physics that will compete with the rotational and hydrodynamic energies that we have considered.Additionally, most synthetic dimension implementations are discrete with hard-wall boundary conditions, and hence are best described by tight binding models on a lattice.It is then important to consider how many synthetic lattice sites are spanned by the typical length scales of the problem.If the answer is many, then a continuum approximation can be appropriate in the mean-field regime.If not, then a tight binding model must be used and rich physics can be expected to arise from competition between these length scales and the synthetic lattice spacing.
Moreover, synthetic dimensions can also have features that are rarely seen in typical tight binding models, and which in themselves warrant further research.These can include nonuniform hoppings, limited system sizes, nonequilibrium effects from external driving and long-range interactions along the synthetic dimension [32,33].All of these are details that should be considered to make this work more experimentally relevant but they also depend strongly on the experiment in question.For this reason, and for simplicity, we have studied a minimal extension of 3D superfluid physics into 4D, in order to begin investigating what is possible in higher dimensions.
An obvious direction of future work is therefore to connect these results to experiment, by studying more complex models that take experimental details into account.We hope that such research can build upon our work by using similar techniques and ansatzes, and that more physical models will yield even richer behaviour.One simple modification that could still have interesting effects is to keep the continuous, isotropic 4D GPE model but to change the geometry to one which breaks the rotational symmetry and better reflects the boundary conditions in a synthetic dimension.A possible choice would be to pick out one or two directions as "synthetic" and give them independent hard wall boundaries (i.e.w ∈ [−L, L], while retaining a rotationally symmetric geometry in the remaining coordinates. There are many other interesting avenues for extending our research, aside from making the model more relevant to experiment.Our numerical stationary states with curved vortex surfaces raises the interesting possibility that other stationary states under rotation could contain closed vortex surfaces that do not meet the boundary of the system.These would be the four dimensional generalisation of vortex loops (including links and knots) in 3D [122][123][124][125][126][127].Additionally, there is a far richer classification of closed surfaces [128] than of closed loops, suggesting there could be more possible closed vortex configurations in 4D.
It would also be interesting to study vortex surface configurations for even higher rotation frequencies.The presence of intersection, curvature, and avoided crossings in our vortex core results suggests that vortices can lose their individual character in 4D.It is therefore not entirely clear, even in some low energy stationary states, whether we can meaningfully assign an integer to the number of vortices in the system.In lower dimensions the number of vortices becomes very large in the rapidly rotating limit, where the vortices form an Abrikosov lattice [8].Investigating the limit of high frequency in one or both planes of rotation in 4D is therefore an interesting and open problem, due to the more malleable nature of the vortex core(s).
This work can also be extended to consider more interesting order parameters in 4D.Certain phases of spinor condensates in 3D are known to host non-Abelian vortices [129,130], which have more interesting behaviour when they intersect and reconnect.Given that intersection and reconnection are also relevant for the behaviour of vortex planes, it is natural to ask what phenomena would arise for non-Abelian vortices in 4D.Finally, this work also represents a small step towards the strongly interacting fractional quantum Hall effect in 4D [131,132], thanks to the analogy between a rotating superfluid and a quantum Hall system [3].
We want to derive the simplest rotation to describe a plane tilting in 4D without loss of generality.Consider the plane P defined in 4D Cartesian coordinates as the set of solutions to x = y = 0, and another plane P ′ as the image of P under a double rotation.We will represent P ′ as the set of solutions to x ′ = y ′ = 0, where the primed coordinates are related to the original coordinates by double rotation with matrix M , that is, r ′ = M r.It will be useful to write this in a block form such that where A, B, C, D are the 2 × 2 blocks of M .Rotations in 4D generally have six free parameters, but we can reduce this down to two for the matrix M by exploiting the symmetry of P under certain rotations, and by using our freedom to choose a basis.Firstly, using the following shorthand for a 2D rotation matrix R(ϕ) = cos ϕ − sin ϕ sin ϕ cos ϕ , note that, for arbitrary ϕ 1,2 , we can redefine M to be without changing P ′ .The reason for this is that the initial rotation we have added is a double rotation in the x-y and z-w planes [c.f.Section III], which leaves the plane P invariant, such that the combined transformation results in the same transformed plane P ′ .Secondly, we will use another double rotation in the x-y and z-w planes to change basis, as follows Denoting this matrix as R(ϕ 3 , ϕ 4 ), we have that M → R(ϕ 3 , ϕ 4 )M R(−ϕ 3 , −ϕ 4 ) under this transformation.Combining this with the redefinition from Eq. (A2) we can write, for arbitrary ϕ j , j = 1, 2, 3, 4 A4) without any loss of generality.Note that we have made the shifts ϕ 1 → ϕ 1 + ϕ 3 , and ϕ 2 → ϕ 2 + ϕ 4 for simplicity.We can use these four free parameters to transform the upper left (A) and lower right (D) blocks into diagonal 2 × 2 matrices.To see this, start by expanding the product in Eq. (A4) Denoting the elements of A in the standard fashion and employing the shorthand s j = sin ϕ j , c j = cos ϕ j , the off-diagonal elements of A ′ = R(ϕ 3 )AR(ϕ 1 ) are given by Setting these both to zero and taking the sum and difference of the two gives the following simultaneous equations where a 1,2 and d 1,2 now denote the only non-zero elements of the upper left and lower right blocks after these blocks have been made diagonal.To proceed further, we will first focus on the upper right B block.Normalization of the first two rows of M can be ensured, without loss of generality, by the following form such that orthogonality of the first two rows now implies This has sin α 1 = 0 or sin α 2 = 0 as special cases, which we ignore for now since these each lead to a simple rotation of the plane P [c.f.Section III].What we will derive instead is the general case for double rotation by requiring β 2 = β 1 + π/2, and this general case will actually include the simple rotation as a special case.Proceeding, we have Orthogonality of the last two columns gives Again, we have a special case, given by α 2 = α 1 , which will give an isoclinic rotation of the plane P [c.f.Section III].We will ignore this solution for now, and again find that it can be found as a particular case of the remaining solution.We therefore require either cos β 1 = 0 or sin β 1 = 0.This leads to the following two forms respectively, up to an unimportant common sign in the upper right block which can be absorbed in to the definition of α 1 and α 2 .Furthermore, these two forms are related to each other by a change of basis and redefinition of parameters.We therefore choose the second form without loss of generality.Orthonormality of the columns and the last two rows now allows us to determine the remaining unknowns, such that we finally have as is used in the main text.

Appendix B: Integration in skew double polar coordinates
In this appendix we derive an integral transformation from 4D Cartesian coordinates (x, y, z, w) into a non-orthogonal coordinate system given by (x, y, z ′ , w ′ ), where the primed coordinates form another Cartesian framed related to the unprimed one by a double rotation and apply this to the vortex-vortex interaction energy.
If we are only interested in preserving the relationship between the two planes defined by x = y = 0 and z ′ = w ′ = 0 respectively, then without loss of generality we can choose this double rotation to have the form as derived in Appendix A. As discussed in the main text, here we will only deal with the special case where this double rotation is isoclinic, such that α 2 = να 1 , with ν = ±1.From here on we will employ the shorthand c = cos α 1 , s = sin α 1 .We will derive this integration over non-orthogonal coordinates for the case of a 4D ball of unit radius, since this geometry preserves the symmetry between the primed and unprimed coordinates.The primed coordinates can be introduced into the integral using Dirac deltas as follows where B d (R) is the ball of radius R centred at the origin in d dimensions, and here R 2 = 1 2 − x 2 − y 2 .Our goal now is to eliminate z and w by evaluating the corresponding integrals.This in turn will define the limits of integration for their primed counterparts.However, this is more easily accomplished in double polar coordinates, whereby where we have defined z 0 = (z ′ −sx)/c, w 0 = (w ′ −νsy)/c, r 2 0 = z 2 0 + w 2 0 , and θ 0 = arctan(z 0 , w 0 ).We now integrate out r 2 and θ 2 as follows where we have used that Θ(R − r 0 ) = Θ(R 2 − r 2 0 ) since both R and r 0 are non-negative.In double polar coordinates we have that R 2 = 1 − r 2 1 , and Substituting this into the Theta function on the RHS of Eq. (B5) gives where we have used that Θ(c 2 •) = Θ(•).Altogether this gives where we have also used the spherical symmetry to restrict the upper limit of r ′ 2 to 1, by comparison to that of r 1 .(This is unnecessary, since the step function will ultimately control the limits of whichever radius is integrated over first, but it makes the equivalence between the primed and unprimed coordinates fully clear.) From now we will assume that the primed coordinates will be integrated over first, so let us make the substitution θ ′ 2 = θ + νθ 1 , treating θ 1 as a constant within the θ ′ 2 integral, in order to simplify the cosine.The limits of the θ integral will be (−νθ 1 − π, −νθ 1 + π), but this is arbitrary since we are integrating over a full circle, so we can just as easily write θ ∈ (−π, π).In order to figure out exactly how the step function translates into integration limits, consider the inequality it enforces This form is ideal for integrating over θ first, but it will actually be easier to integrate over r ′ 2 first.For this reason we will rewrite Eq. (B9) by completing the square for r This inequality has no solutions for r ′ 2 where the RHS is negative, so we immediately obtain as a constraint for θ.Note that this constraint is trivially satisfied whenever c 1 − r 2 1 1/2 > sr 1 , which occurs when r 1 < c.Given this condition for θ, we can then satisfy the inequality (B10) when r ′ 2 ∈ (r − , r + ), where The last step is to enforce the constraint r ± ≥ 0, since r ′ 2 cannot be negative.Rearranging each inequality gives where χ = c 1 − r 2 1 1/2 /sr 1 .Note that the quantity on the RHS of both of these inequalities is always nonnegative, so r − ≥ 0 requires θ ∈ [−π/2, π/2], while r + ≥ 0 is automatically satisfied in this same region.With this consideration of the sign of the LHS in mind, we can square both inequalities and rearrange to find Combining all of this with the inequality (B11), gives us two separate integration regions.We have and where θ * = arcsin χ.Finally, we can write the full result as As stated in the main text, the vortex-vortex interaction energy is then given as where we have introduced J 1 and J 2 as shorthand to denote the two integrals.We will now deal with each of these integrals separately; focusing on the first term, we have where we have carried out the integral over r ′ 2 .The second integral on the RHS of this equation can be shown to vanish as follows which works for any arbitrary function F .This leaves us with We now turn to the second term of Eq. (B21), which depends on This has the form of the vanishing term in J 1 , except that the θ limits now do not cover a full period.In fact, the limits do not cover even half a period since θ * ≤ π/2 (consider Fig 2 with the blue dotted circle passing through the origin), with the consequence being that this term now contributes.To compute it we will apply the substitution sr 1 sin θ = c 1 − r 2 1 1/2 sin u, to give Combining these results then gives the final result as stated in the main text.
Appendix C: The minimum distance between two circles of common centre and radius in 4D Consider a pair of circles (C 1 , C 2 ) in R 4 with the same radius and centre, but occupying different planes, and let these two planes be related by an isoclinic rotation.In this appendix we will derive an expression for the minimum distance between two such circles.Without loss of generality, we may encode the two circles in the following vector equations: The vector between an arbitrary point on C 1 and an arbitrary point on C 2 is given by d = r C1 − r C2 .All we have to do is compute the length of this vector and minimize it with respect to θ 1 and θ ′ 2 .Evaluating the modulus squared of d, we have The minimum value of d is therefore √ 2R(1−s) 1/2 , which occurs when θ 1 = νθ ′ 2 + π.

Appendix D: Evaluation of rotational energy integrals
As stated in the main text, in order to calculate the rotational energy of two vortex planes under unequal frequency double rotation, we must compute the integral (Eq.( 120)): as well as the corresponding integral for e i θ2 .Here, we will only show the direct calculation of the first integral, since the second follows identical logic.To begin, we will consider acting with Lxy ≡ −iℏ∂ θ1 on Eq. (115) as Lxy ŕ1 e i θ1 ŕ1 e i θ1 = ℏ cos η 1 r 1 e iθ1 cos η 1 r 1 e iθ1 + e iφ sin η 1 r 2 e iθ2 , (D2) where we have also divided through by Eq. (115).Then by using the product rule, we can see that the desired integrand in Eq. ( 120) can be expressed as Lxy e i θ1 e i θ1 = ℏr 1 e i(θ1−θ2−φ) r 1 e i(θ1−θ2−φ) + tan η where the limits reflect that the 4D hypersphere is bounded by r 2 1 + r 2 2 = R 2 .Since both r j are nonnegative we can safely rewrite the step function as Θ(r 2 1 − tan 2 η 1 r 2 2 ), which allows us to make the substitutions u j = r 2 j /R 2 , such that the integral then becomes (D8) It is now much easier to compare the step function to the integration limits, as we have that u 1 > tan 2 η 1 u 2 and 0 < u 1 < 1 − u 2 , whilst 0 < u 2 < 1. Clearly tan 2 η 1 u 2 is non-negative, so we can make this value the new lower limit for u 1 provided it is not greater than the upper limit of 1−u 2 .This will be true for a certain range of u 2 values which satisfy the inequality tan 2 η 1 u 2 ≤ 1 − u 2 , (D9)  That is, we fixed Ωzw, and varied ∆Ω starting from Ωc and stepping down to the isoclinic point in units of 0.1Ωc.We used the final state of each run as the initial state of the next.These two states started out very close to the predicted skew planes at high frequency difference, however as we decreased ∆Ω these planes began to increasingly curve, up to the maximum amount shown in this figure at ∆Ω = 0.These states both have energy given by E ≈ 0.65µN , which is approximately equal to the orthogonal state energy at these frequencies, so we would need higher precision to tell if these are degenerate with the orthogonal state or if they are metastable excited states.
tively.Evaluating this, we obtain In this appendix, we will present extra numerical results to supplement those in the main text.Some of these data are from simulations not mentioned in the text, and others are additional data from simulations described in the text, to aid in explanation and visualisation.
Fig 8 shows the final state vortex core for a run of the ITEM with frequencies Ω zw = 1.25Ω c , and Ω xy ≈ 1.35Ω c , which gives a predicted skewness angle of η = 10 • .Just as in Fig 5, the agreement between theory and numerics is still very good apart from the avoided crossing region near the origin.In particular, the third panel shows a side on view in which the vortex cores lies roughly along a pair of lines, on top of which we have plotted the theoretically predicted lines we expect.As can be seen, the data is still very close to the predicted values, although interestingly there is a small degree of asymmetry visible in the third panel of the data points are not symmetric about the vertical axis -that can't be seen in Fig 5 .However, without further investigation we cannot tell whether this is due to numerical inaccuracy or some genuine physical phenomena., with a spatial resolution of ∆x = 0.2ξ, compared to 0.5ξ for the previous figures.The agreement between these data and our predictions is as good as before, as can be seen in the third panel, and the data does have mirror symmetry about the vertical axis.However, there are also some interesting qualitative differences to the lower resolution results.For example, the avoided crossing between the planes has a different orientation to that in the previous two figures, as can be seen by comparing the first two panels.Finally, the size of the reconnected region is smaller, which is most likely due to the smaller system size, but could be due to the higher resolution.The former x-y plane is now so curved that it has forced the other surface to become tilted and displaced, and there is an avoided crossing between them as visible in the first panel.Note that the second panel shows that parts of this surface appear to be parallel to each other.Here, the energy is approximately the same as the predicted skew plane energy at E ≈ 0.64µN .served at the isoclinic point for identical parameters.The frequencies were Ω zw = Ω xy = 1.3Ω c , and the spatial step size was ∆x = 0.5ξ, which corresponds to a radius of R ≈ 20.6ξ.Both of these final states were the last iteration in a sequence of ITEM runs, starting from ∆Ω = Ω c down to ∆Ω = 0 in steps of 0, 1Ω c , with Ω zw = 1.3Ω c fixed.The final state of each run was used as the initial state of the next, so that we could follow the evolution of a particular energy branch.We were attempting to track the predicted skew plane states for a fixed Ω zw = 1.3Ω c and find their energies [c.f.Fig 6], and the two final states at ∆Ω = Ω c did in fact closely correspond to these predicted planes.However, both of these states began to deviate from the theoretical states as we decreased ∆Ω, becoming more and more curved all the way to the isoclinic point.These isoclinic curved states have approximately the same energy as each other, and as the orthogonal state at the same frequencies (E ≈ 0.65µN ).Since any energy differences are below the precision of our numerics we would need more accuracy to investigate this.If these are indeed low lying excitations, this may be due to the degeneracy associated with isoclinic symmetry, in which case we may expect many more such states.
Finally, Fig. 11 shows the vortex cores obtained using the asymmetric phase profile (with added noise), for ∆Ω values corresponding to η = 10 • , 20 • and 40 • , respectively.As mentioned in the main text, these are slightly higher in energy than the theoretical states studied in the main text.Already, in the 10 • case we can see that our assumption of flat vortex planes is broken as there is some long-range curvature of the core.This is then exacerbated as the frequency ∆Ω increases, with the 20 • state clearly showing that there is almost no overall tilt of the former planes, but instead the plane at z = w = 0 has begun to buckle in an approximately threefold symmetric pattern, curving towards the other surface in different directions for different values of the angle θ 1 .Finally, in the 40 • figure, this curvature has become so extreme that the former plane at x = y = 0 appears to have become tilted and displaced as well as curved, leading to a sizeable avoided crossing where these two surfaces come together, as can be seen in the first panel.Interestingly, the second panel appears to show three parts of the vortex core surface that are parallel to each other.This suggests that this bizarre state may be limiting towards a state with multiple vortex planes parallel to the z-w plane but separated in the x-y plane.This is the expected lowest energy state for the case of high frequency simple rotation in the x-y plane (i.e. with Ω zw = 0) so it seems reasonable that it should also be the ground state when Ω xy ≫ Ω zw .However, we are not quite reaching this limit in the 40 • case, as there we have Ω xy ≈ 1.46Ω zw .We therefore tentatively describe this strange set of states as an instance where the frequencies are large enough that the system wants more than two vortex planes but not enough for three.The planes can curve in order to become larger, thereby fitting a larger area of vortex surface in the system.Whether this is a correct description or not, it is clear that the behaviour of vortex surfaces in 4D is incredibly rich, and there is much more to be explored.

FIG. 3 .
FIG.3.A possible configuration of tilted vortex planes relative to the x-y and z-w planes, visualised as lines in 2D, for special values of φ.When φ = 0 the planes are given by Eqs (131) and (132), which describe a pair of lines in x, z space and an identical pair of lines in y, w space, plotted simultaneously in the left panel.When φ = π/2 we have instead Eqs (133) and (134), describing a pair of lines in x, w space and their reflection about the horizontal axis in y, z space, which we plot on the same graph here by transforming z, w → w, −z.Equivalent pictures for φ = π and 3π/2 respectively can be found by reflection about the vertical axis.

FIG. 4 .
FIG.4.Difference in dimensionless energy density [Eq.(137)] between the skew vortex plane ansatz and the orthogonal configuration as a function of the tilt angles η1,2.The dimensionless frequency difference has been set to ω = 2, for this value the minimum occurs at a skewness of η ≈ 38 • as calculated from Eq. (151).Contour lines are included as guides to the eye to highlight the degeneracy along the line η1 = η2, and the minimum in the upper left.We have omitted values of ε > 2 to avoid the logarithmic divergences at (0, π/2) and (π/2, 0).

FIG. 8 .
FIG. 8.As in Fig 5, except with rotation frequencies of Ωzw = 1.25Ωc and Ωxy ≈ 1.35Ωc, corresponding to a predicted skewness of η = 10 • .The overall structure is essentially the same as in Fig 5, and the agreement with theory is still very good as shown in panel three.However, this panel also shows that the data has a slight mirror asymmetry in the vertical axis that is not visibly present in Fig 5 or accounted for in the theory.

FIG. 9 .r 1 −π dθ 1 Lxy e i θ1 e i θ1 = ℏ |ζ|=1 −ir 1 dζ r 1 ζ + tan η 1 r 2 = 4 re i θ1 e i θ1 = 4ℏπ 2 R 0 r 2 dr 2 √R 2 −r 2 2 0r 1
FIG. 9.As inFig 8,  except with a higher resolution of ∆x = 0.2ξ, and with Ωxy ≈ 2.25Ωc, corresponding to a predicted skewness of η = 45 • .The agreement with the theoretical planes is still excellent, as seen in the third panel, and there is no visible asymmetry about the vertical axis.Interestingly, the avoided crossing at the origin has a different orientation to that seen in Figs5 and 8, and appears to be smaller.

4 r
D10) in which case the above integral including the step function is equivalent toB 4 (R) d Lxy e i θ1 e i θ1 = ℏπ 2 R 4where c and t are shorthand for cos η 1 and tan η 1 , respec-

FIG. 10 .
FIG.10.(Top and bottom) Two different unusual vortex core structures observed in final states of ITEM runs with the same parameters.The frequencies were set at the isoclinic point, Ωzw = Ωxy = 1.3Ωc, and the spatial resolution was equal to ∆x = 0.5ξ, corresponding to system radius of R ≈ 20.6ξ.These simulations were the final iteration in a series of ITEM runs similar to those that generated Fig 6.That is, we fixed Ωzw, and varied ∆Ω starting from Ωc and stepping down to the isoclinic point in units of 0.1Ωc.We used the final state of each run as the initial state of the next.These two states started out very close to the predicted skew planes at high frequency difference, however as we decreased ∆Ω these planes began to increasingly curve, up to the maximum amount shown in this figure at ∆Ω = 0.These states both have energy given by E ≈ 0.65µN , which is approximately equal to the orthogonal state energy at these frequencies, so we would need higher precision to tell if these are degenerate with the orthogonal state or if they are metastable excited states.

Fig 9
Fig 9 shows the final state for η = 45• , with a spatial resolution of ∆x = 0.2ξ, compared to 0.5ξ for the previous figures.The agreement between these data and our predictions is as good as before, as can be seen in the third panel, and the data does have mirror symmetry about the vertical axis.However, there are also some interesting qualitative differences to the lower resolution results.For example, the avoided crossing between the planes has a different orientation to that in the previous two figures, as can be seen by comparing the first two panels.Finally, the size of the reconnected region is smaller, which is most likely due to the smaller system size, but could be due to the higher resolution.