Superfluid Vortices in Four Spatial Dimensions

Quantum vortices in superfluids have been an important research area for many decades. Naturally, research on this topic has focused on two and three-dimensional superfluids, in which vortex cores form points and lines, respectively. Very recently, however, there has been growing interest in the quantum simulation of systems with four spatial dimensions; this raises the question of how vortices would behave in a higher-dimensional superfluid. In this paper, we begin to establish the phenomenology of vortices in 4D superfluids under rotation, where the vortex core can form a plane. In 4D, the most generic type of rotation is a"double rotation"with two angles (or frequencies). We show, by solving the Gross-Pitaesvkii equation, that the simplest case of equal-frequency double rotation can stabilise a pair of vortex planes intersecting at a point. This opens up a wide number of future research topics, including unequal-frequency double rotations; the stability and reconnection dynamics of intersecting vortex surfaces; and the possibility of closed vortex surfaces.

Quantum vortices are fundamental topological excitations of superfluids, which have been widely studied for many years [1][2][3][4][5][6][7]. Unlike a lot of many-body phenomena, vortices can be understood at the mean-field level through the Gross Pitaevskii equation (GPE) [1]. A superfluid vortex consists of a local density depletion within the "vortex core", around which the superfluid circulates. In 2D and 3D superfluids, this vortex core forms a point and a line respectively, as sketched in Fig 1. Vortices have an associated energy cost, but can be stabilised by rotation of the superfluid [2,3], or equivalently by artificial magnetic fields [8][9][10].
The potential of synthetic dimensions for reaching 4D with (for example) ultracold bosonic atoms [23,33] motivates the question of how superfluid vortices behave in higher dimensions. In this paper, we take an initial step in this direction by exploring the 4D GPE under rotation, with local atom-atom interactions. This is chosen as a * bdm375@student.bham.ac.uk Quantum vortices in superfluids have been an important research area for many decades. Naturally, research on this topic has focused on two and three-dimensional superfluids, in which vortex cores form points and lines, respectively. Very recently, however, there has been growing interest in the quantum simulation of systems with four spatial dimensions; this raises the question of how vortices would behave in a higher-dimensional superfluid. In this paper, we begin to establish the phenomenology of vortices in 4D superfluids under rotation, where the vortex core can form a plane. In 4D, the most generic type of rotation is a "double rotation" with two angles (or frequencies). We show, by solving the Gross-Pitaesvkii equation, that the simplest case of equal-frequency double rotation can stabilise a pair of vortex planes intersecting at a point. This opens up a wide number of future research topics, including unequal-frequency double rotations; the stability and reconnection dynamics of intersecting vortex surfaces; and the possibility of closed vortex surfaces.
Quantum vortices are fundamental topological excitations of superfluids, which have been widely studied for many years [1][2][3][4][5][6][7]. Unlike a lot of many-body phenomena, vortices can be understood at the mean-field level through the Gross Pitaevskii equation (GPE) [1]. A superfluid vortex consists of a local density depletion within the "vortex core", around which the superfluid circulates. In 2D and 3D superfluids, this vortex core forms a point and a line respectively, as sketched in Fig 1. Vortices have an associated energy cost, but can be stabilised by rotation of the superfluid [2,3], or equivalently by artificial magnetic fields [8][9][10].
The potential of synthetic dimensions for reaching 4D with (for example) ultracold bosonic atoms [23,33] motivates the question of how superfluid vortices behave in higher dimensions. In this Letter, we explore this by studying the 4D GPE equation under rotation, with local * bdm375@student.bham.ac.uk

2D
3D 4D 1. (Colour online) Sketch of minimal vortex structures, stabilised for different system dimensionalities (columns) and types of rotation (rows). Here, "simple" and "double" indicate rotations with one or two planes of rotation respectively, as discussed in the text. In 2D and 3D, only simple rotations exist, stabilising vortex cores as a point and line, respectively, about which the superfluid rotates (black arrow). In 4D space (shown as 3D cross-sections coloured according to w value), both types of rotation exist, leading to a richer vortex phenomenology. In 4D, equal-frequency double rotations can lead to a new type of vortex configuration consisting of two vortex planes intersecting at a point, while simple rotations stabilise a single vortex plane. In these sketches, a vortex plane appears either as a line persisting for all w (lines of varying colour), or as a plane for a particular w value (purple disc), depending on the rotation plane. Note that in the 4D column we have omitted the arrow indicating superfluid motion.
atom-atom interactions. The 4D GPE equation can be justified physically as a description of low-temperature interacting bosons in 4D [48][49][50][51]. It is also a minimal model, which allows us explore the simplest examples of 4D vortices without complications that depend on how the synthetic dimension is implemented [25][26][27][28][29][30][31][32]35]. The inclusion of more realistic experimental details, such as lattices, anisotropies, and long-range interactions in the FIG. 1. Sketch of minimal vortex structures, stabilised for different system dimensionalities (columns) and types of rotation (rows). Here, "simple" and "double" indicate rotations with one or two planes of rotation respectively, as discussed in the text. In 2D and 3D, only simple rotations exist, stabilising vortex cores as a point and line, respectively, about which the superfluid rotates (black arrow). In 4D space (shown as 3D cross-sections coloured according to w value), both types of rotation exist, leading to a richer vortex phenomenology. In 4D, equal-frequency double rotations can lead to a new type of vortex configuration consisting of two vortex planes intersecting at a point, while simple rotations stabilise a single vortex plane. In these sketches, a vortex plane appears either as a line persisting for all w (lines of varying colour), or as a plane for a particular w value (purple disc), depending on the rotation plane. Note that in the 4D column we have omitted the arrow indicating superfluid motion.
minimal model, which naturally extends a standard textbook problem to 4D in order to establish basic aspects of 4D vortex physics. More realistic models for experiments will depend on the specific synthetic-dimension implementation chosen, and are likely to include other effects, such as lattices and unusual interactions with respect to the synthetic dimension, that will further enrich the possible vortex states, but will go beyond the current work. We also note that while our main motivation for studying the 4D GPE is as an initial stepping-stone towards pos-sible synthetic-dimension experiments, this model is also plausible as a description of low-temperature interacting bosons in a hypothetical 4D universe (see Appendix A, and [48][49][50]), and so is of mathematical interest for generalising classic results about superfluid vortices to higher dimensions.
To investigate vortices in 4D, we must first appreciate that rotations (or equally, magnetic fields) in higher dimensions can have a fundamentally different form; all rotations in two and three dimensions are so-called "simple rotations", while in 4D, generic rotations are "double rotations" [51]. This difference will be discussed in more detail later, but can be understood in brief by noting that in 2D/3D every rotation has a single rotation plane and angle, while in 4D there can be two independent rotation planes, e.g. the xy and zw planes, each with their own angle of rotation.
In this paper, we show that equal-frequency double rotation of a 4D superfluid can stabilise a vortex structure formed by two vortex planes intersecting at a point, while a simple rotation stabilises a single vortex plane, as sketched in Figure 1. We obtain our results, firstly by using a phase ansatz to numerically solve an effective 2D radial equation, and secondly by numerically solving the full 4D GPE under rotation. This generalisation of superfluid vortices to higher dimensions opens up many avenues of future research, such as questions concerning the unequal-frequency case; reconnections of vortex planes; possible curvature of vortex surfaces; and more realistic setups capturing experimental details.

I. REVIEW OF VORTICES IN 2D AND 3D SUPERFLUIDS
We begin by reviewing the basic properties of 2D and 3D vortices, in order to lay the groundwork for our discussion of 4D superfluids. We consider systems of weaklyinteracting bosons described by a complex order parameter, ψ, which obeys the time-independent GPE with no external potential [1] where m is the particle mass, g is the interaction strength, and µ is the chemical potential. A hydrodynamic description can be obtained from this equation by substituting ψ = √ ρe iS , where ρ is the superfluid density, and S is the phase [1]. The velocity field, v = m ∇S, is irrotational wherever S is well behaved. A consequence of this property is that a superfluid supports quantized vortices. This can be seen by noting that the superfluid circulation around a closed loop C is quantised as where [∆S] C is the phase winding [3]. Since ψ is singlevalued, we must have [∆S] C = 2πk, where k is the in-teger winding number (or vortex charge) [1]. Smoothly deforming the loop cannot change k as long as vortices are avoided. This can only be true if v diverges like 1/r as the distance r from a vortex core goes to zero. Since particles cannot have infinite velocity ρ must vanish in this same limit. The region of density depletion is known as the vortex core; in 2D, this is localised around a point, and in 3D around a line, as shown in Fig 1. More generally, vortices must be localised in two directions.
As is well known, the density profile around the vortex core can be calculated directly by applying the GPE to a homogenous superfluid with a single vortex [1]. By defining the uniform background density n, the healing length ξ can be introduced, which satisfies 2 /mξ 2 = gn = µ [52], and which physically is the distance over which ρ typically varies. Hereafter, we rescale r → ξr, and ψ → √ nψ such that Eq (1) becomes dimensionless as A rotationally symmetric vortex state in 2D has the form ψ = f k (r)e ikθ , where (r, θ) are polar coordinates centred on the vortex core, f k (r) is real, and k is the winding number. Substituting this into Eq 3 gives [1] where ∆ r = ∂ 2 /∂r 2 + (1/r)∂/∂r. This equation has no closed-form solution, but does admit the asymptotic forms f k (r) = O(r |k| ) as r → 0, and f k (r) = 1−O(r −2 ) as r → ∞ [3]. The crossover between these two behaviours occurs at around the healing length. Note that a straight vortex line in an otherwise homogeneous and isotropic 3D superfluid has this same profile, with (r, θ) defined in the plane perpendicular to the vortex line [1]. Using this density profile the energy cost of a vortex relative to the ground state can be evaluated. For a singly charged vortex (k = 1) the energy can be written as where N is number of bosons, and R is the radius of the superfluid in the plane orthogonal to the vortex core. Eq (5) is valid in any number of dimensions. Vortices can be energetically stabilised by rotation (or equivalently an artificial magnetic field), whereby Eq (1) is modified in 3D by adding the term −Ω·Lψ to the left hand side, with L = −i r × ∇ the angular momentum operator, and Ω the frequency vector [1]. This term reduces the energy of a state containing a vortex aligned with the rotation, making it more energetically favourable.

II. SIMPLE AND DOUBLE ROTATIONS
Given the intrinsic link between rotation and vortices, we will now discuss the different types of rotations possi-ble in 4D, as compared to lower dimensions, in preparation for our discussion of vortices in 4D superfluids below.
In three dimensions or fewer, every rotation is "simple"; this means that the rotation is specified by a rotation angle α ∈ (−π, π], and a plane of rotation which is unique up to translation. Under rotation, the points on the plane of rotation remain on the plane, but are displaced through the angle α. Generalising to D dimensional space, simple rotations have D − 2 eigenvectors with eigenvalue one, all of which are orthogonal to every vector in the rotation plane. For example, a rotation about the z axis in 3D has the xy plane (defined by z = 0) as its rotation plane, and fixes any point along the z axis. We may write this as a matrix in the standard basis as We can think of this as a rotation of 2D space (spanned by x and y) extended into a third (z) direction. Similarly, simple rotations in 4D can be thought of as rotations of 3D space extended into a fourth direction. Labelling the fourth axis as w, our previous example becomes a rotation about the zw plane (defined by x = y = 0), given in matrix form by and I is the 2D identity. Note that there are six Cartesian coordinate planes in 4D, so the rotation group SO(4) has six generators, and the representation of these generators (which physically describe angular momentum) as spatial vectors no longer works in 4D as it does in 3D. The set of fixed points of a simple rotation in 4D are a plane, not a line, and this fixed plane is completely orthogonal to the plane of rotation, by which we mean that every vector in one plane is orthogonal to every vector in the other. In contrast to 2D and 3D, in four dimensions, we can also have "double rotations", which generically have only one fixed point, and two completely orthogonal planes of rotation each with a corresponding rotation angle [51]. To visualise this, consider a double rotation in the xy and zw planes represented by the matrix [53] for angles α, β ∈ (−π, π]. For those familiar with certain 4D quantum Hall models, this is analogous to generating a second Chern number by applying magnetic fields in two completely orthogonal planes [12,13,23,24,54]. Double rotations are in fact the generic case of rotations in 4D, as if either α or β = 0, the rotation reduces to the special case of simple rotation discussed above [51]. From here on we will refer to the two planes of rotation as planes 1 and 2 respectively and focus only on so-called "isoclinic" double rotations for which α = β. Before continuing, it is worth noting that isoclinic rotations have an additional symmetry. To see this, we remember that, as introduced above, generic double rotations have one fixed point and two planes of rotation, with corresponding angles α, β ∈ (−π, π]. Vectors in R 4 which do not lie in these rotation planes are displaced through an angle between α and β [51]. However, if α = β, then this means that every vector is displaced by the same angle. As a consequence, for a given isoclinic rotation there is a continuum of pairs of completely orthogonal planes that can each be though of as the two planes of rotation. In other words, isoclinic rotations therefore no longer have two unique planes of rotation, although they still have a single fixed point. However, numerically we break this degeneracy since the phase winding of our initial state picks out the xy and zw planes in particular. We can also anticipate that a more experimental model would likely break this symmetry too, e.g. through the inclusion of lattices or through inherent differences between real and "synthetic" spatial dimensions.

III. VORTEX PLANES IN 4D
Now that we have discussed some of the geometry of rotations in 4D we are ready to study the associated vortex physics. As above, we consider a superfluid described by the GPE in the absence of external potentials, but now with atoms free to move in four spatial dimensions.
The simplest case to consider is that of a 4D superfluid under a constant simple rotation. As shown in Eq (7), a simple rotation can be viewed as a 3D rotation extended into a fourth dimension, hence stabilising a vortex plane, as sketched in Fig 1. The corresponding order parameter profile is ψ = f k (r 1 )e ikθ1 , where (r 1 , θ 1 ) are plane polar coordinates in the plane of rotation, and f k (r) is the solution of Eq (4). As this is independent of the other two coordinates, the vortex core becomes a plane; this is di-rectly analogous to the extension of a point vortex in 2D into a line in 3D. We have verified this result numerically, as shown in Appendix B.1. This can be understood as the natural extension of vortices into 4D, as the extra dimension plays no role, and the vortex plane is homotopically characterised by a Z topological winding number, as in 2D and 3D. For a more detailed discussion of homotopy classification of vortex planes in 4D, see Appendix C.
In contrast we expect that double rotation, being an intrinsically 4D (or higher) phenomenon, will lead to more interesting vortex configurations. To address this problem, we look for the ground states of the 4D GPE in a doubly rotating frame where Ω j and L j are the rotation frequency and angular momentum operator in plane j. In Cartesian coordinates (x, y, z, w), L 1 = −i (x∂ y − y∂ x ), and L 2 = −i (z∂ w − w∂ z ). For simplicity we will 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 ), The simple rotation case discussed before corresponds to Ω 2 = 0, where the vortex core spans plane 2. In this paper we focus on equalfrequency doubly rotating superfluids, that is Ω ≡ Ω 1 = Ω 2 . The fact that L 1 and L 2 generate a double rotation means that they commute. We may look for a solution which is a simultaneous eigenstate of both angular momentum operators; therefore we propose an ansatz for the ground state under rotation of the form where f (r 1 , r 2 ) is real and the k j are integer phase winding numbers in each rotation plane. This phase profile corresponds to the superfluid circulating in both planes simultaneously, about both vortex cores. We have suppressed the dependence of f on each k j for brevity, and in all numerical results both winding numbers are one. This state exhibits a phase singularity when either r j = 0, so we require f (0, r 2 ) = f (r 1 , 0) = 0 from the same reasoning as in 2D and 3D. In other words, this describes a pair of completely orthogonal vortex planes that intersect at a single point as illustrated in Fig 1, and which are characterised by Z × Z topological winding numbers (see Appendix C). Intersection of two planes at a point is only possible in 4D or higher and, in fact, is the generic case in 4D. This is in contrast with 3D, where the intersection of lines is a special case, and so vortex lines intersect and reconnect at specific times [55][56][57][58].
To examine our ansatz, we now proceed to numerically solve for the density profile, under this phase constraint. Substituting the ansatz [Eq (10)] into the GPE [Eq (1)] in 4D, and de-dimensionalising in the same way as in the 2D case, we obtain the following equation for f (r 1 , r 2 ) (11) where ∆ rj = ∂ 2 /∂r 2 j + (1/r j )∂/∂r j . Since each vortex produces only a local density depletion, we expect that f (r 1 , r 2 ) ∼ f k2 (r 2 ) as r 1 → ∞ and equally for (1 ↔ 2), where f k (r) is the point vortex solution of Eq (4). Note that this limiting "boundary condition" can be satisfied by a product, f k1 (r 1 )f k2 (r 2 ), of 2D density profiles in each plane. However, this form fails to solve the full equation due to the non-linear f 3 term. This product form therefore gives a natural approximation to compare to, and we expect it to fail significantly only in the vicinity of the origin, where both f kj (r j ) differ appreciably from unity.
To verify this, and find the full density profile, we have solved Eq (11) by imaginary time evolution within a discretised grid in (r 1 , r 2 ) space with hard-wall boundary conditions at a radius R = 100ξ in each plane (r j = R), and at the origin in each plane (r j = 0). The latter condition is required due to the centrifugal term diverging at the vortex cores; consequently the precise location of the vortex cores was an assumption in these calculations. We used a forward Euler time-discretization and second order finite differences in space. We chose a large value of R compared to ξ so that we could examine the vortex cores within a homogeneous region. (Future studies could include the effect of additional trapping potentials, such as harmonic traps along some or all directions.) We were able to achieve a resolution of 0.05ξ, and the calculations were converged until the relative change in chemical potential and particle number over one timestep converged below 10 −14 .
The results for k 1 = k 2 = 1 are shown in Fig 2(a), where we observe the expected local density depletion around the vortex cores when either r 1 = 0 or r 2 = 0. We also compare our numerical solution with the product approximation, f 1 (r 1 )f 2 (r 2 ), in Fig 2(b); we observe that the product approximation is very accurate except within a distance of roughly ξ from the intersection point, as expected. Immediately around the intersection, the product approximation fails, overestimating the density by a factor of about 4/3.
Just as in the 2D case we can use our calculation of the density profile to find the energy of this vortex configuration relative to the state with no vortices. Defining independent radii R j in each plane, such that r j ≤ R j , we find numerically (see Appendix B.4) that the energy is approximately given as where E k (R) is the single-vortex energy given in Eq (5). This can be understood from the superfluid kinetic energy ρv 2 d 4 r, which is the main contribution to the energy of a vortex. The velocity field is given by v = v 1 +v 2 where v j = kj rjθ j is the velocity induced by vortex j. As v j lies in plane j, we see that v 1 · v 2 = 0 and so the hydrodynamic vortex-vortex interaction term, ρv 1 ·v 2 d 4 r, vanishes. The total kinetic energy integral therefore splits into a sum of the individual kinetic energies. Note that this argument relies on the assumptions that the two vortex cores have no curvature and are completely orthogonal to each other. In order to confirm the existence and stability of the intersecting vortex plane state we have performed imaginary time evolution with the 4D GPE under both simple and double rotation [Eq (9)] directly on a 4D Cartesian grid within a 4D ball of radius R = 8.25ξ with a hard-wall boundary. A hyper-sphere rather than a hyper-cube was chosen as the majority of the 4D volume of a hyper-cube is taken up by regions "in the corners", that is, outside of the hyper-sphere that just fits inside. This allowed us to relax our above constraint on the phase profile, at the cost of smaller numerical system sizes. Again, we used the forward Euler method for time-discretization and second order finite differences in space. We were able to obtain resolutions of up to 0.2ξ, and by repeating simulations at different resolutions, we checked that our main conclusions were qualitatively insensitive to the coarse-graining of the numerics. At the system sizes and resolutions we have been able to reach, the homogeneous region extends over a few healing lengths. The calculations were converged to an accuracy threshold of 10 −12 .
A benefit of performing calculations with all four coordinates is that we were able to test our ansatz by allowing the phase to evolve, and by removing the boundary condition at r j = 0 mentioned previously. More precisely, we used an initial state with homogeneous density away from the edge of the ball, and a phase profile given by arctan2(y, x) + arctan2(w, z), for the doubly rotating case, and arctan2(y, x) for the singly rotating case. We tested the robustness of our results to noise (up to 20% of the background value) added to the real and imaginary parts of the initial ψ. Note that we measure the applied frequency in units of the critical frequency of a single vortex in a homogeneous 2D disk of the same radius as our 4D ball; this is given (in our units) by [59] Ω 2D crit = µ log(2.07R/ξ) For the results shown in Fig 3 both the frequencies of rotation used were roughly 2.5Ω 2D crit . Further work could investigate the effect of double rotation with unequal frequencies.
For a suitable range of frequencies Ω we find good agreement between the stationary state obtained from the full 4D numerics and our ansatz for two intersecting vortex planes, as shown in Fig 3. Panel (a) shows that the phase profile of the state after relaxation perfectly agrees with that of the ansatz. Panels (b) and (c) show the density and phase profiles, respectively, for the 2D cut in which y = w = 0. As can be seen the density drops to zero along the lines x = 0 and z = 0, corresponding to the intersections of each vortex core with the plane of the cut, as expected. Further two dimensional cuts of this state are given in Appendix B.2.

IV. DISCUSSION AND CONCLUSIONS
In this paper, we have shown that the simple rotation of an idealised 4D superfluid can stabilise a vortex plane, while equal-frequency double rotations can lead to two vortex planes intersecting at a point which do not interact hydrodynamically. This significantly extends the phenomenology of superfluid vortices, demonstrating that new effects can emerge in higher spatial dimensions even within mean-field theory.
It is important to note that we have studied an idealised model, which allows us to explore vortex physics in 4D without experimental details that depend on how the synthetic dimension is implemented [25][26][27][28][29][30][31][32]35]. The main differences between our work and possible experiments are, firstly, that the majority of practical implementations would lead to (tight-binding) lattice models, whereas we have considered four continuous dimensions as a theoretical first step. Adding a lattice should introduce rich additional effects particularly when the lattice spacing is comparable to or greater than other length scales. However when this spacing is very small, it should be possible to approximate a lattice model with a continuum model in the mean-field regime as we have considered here. Furthermore, synthetic-dimension schemes can include unusual effects, which are very dependent on the specific experimental implementation. In terms of the tight binding description previously mentioned, these complications can include position-dependent hopping strengths, limited numbers of sites, and long-range interactions [25][26][27][28][29][30][31]35]. For the sake of generality as well as simplicity we have therefore chosen an idealised model, which can then be adapted in different ways for promising experimental scenarios in further work.
We also note that Eq (1) has SO(4) (4D rotational) symmetry, which would be broken in any experiment due to inequivalence of the synthetic and real spatial dimensions. Numerically, we break this symmetry with the phase anzatz, which was assumed in the radial case, and imposed on the initial state in the Cartesian case. However, we do still assume an SO(2) symmetry in each of the xy and zw planes to obtain the effectively 2D radial equation [Eq (11)], and simplify the corresponding numerics. In the Cartesian case, we also chose a boundary condition (a hard-wall at some radius from the origin) that preserved these in-plane symmetries. In synthetic dimension experiments, on the other hand, the most common boundary condition is an open boundary condition which is independent of the other dimensions [25][26][27][28][29][30][31][32]35]. Hence a more experimentally relevant geometry would involve one or more dimensions which have their own independent hard-wall boundary conditions -for example a "spherinder" boundary specified by {r ∈ R 4 |x 2 + y 2 + z 2 = R 2 } ∪ {r ∈ R 4 |w = ±L} for some R and L. Investigating the effect of breaking one or both of these in-plane rotational symmetries geometrically is an interesting and natural next step for future work.
As well as a first step towards understanding future experimental models, this work also opens up many interesting theoretical research directions. Natural next steps include the study of 4D superfluids doubly rotating at unequal frequencies, and 4D generalisations of previously studied questions from 2D and 3D [2,3]. Firstly, closed vortex surfaces in 4D would naturally generalise the vortex loops that arise in 3D [1], but with potentially an even richer classification when non-orientability and surfaces of higher genus are included [60]. Secondly, vortex lines in 3D are known to dynamically reconnect upon intersection [55][56][57][58], whereas here we have shown that completely orthogonal intersecting vortex planes in 4D form a stationary state stabilised by rotation. It is an open question whether vortex planes reconnect if they are not completely orthogonal, and this question could have relevance to the general case of unequal-frequency double rotation. For example, intuitively, we would expect an adiabatic change from Ω 2 = Ω 1 to Ω 2 > Ω 1 would cause the vortex in plane 2 (inducing rotation in plane 1) to tilt towards plane 1 to benefit from the now larger rotational energy discount in plane 2. Finally, in the longer term this work opens up questions related to the inclusion of strong interactions and the 4D fractional quantum Hall effect, as well as the study of models with more interesting order parameter spaces [61,62], potentially hosting non-Abelian vortices.
The 4D GPE is a natural and mathematically simple generalisation of the 3D GPE, allowing for easy comparison to superfluid vortex physics in lower dimensions. In this section, we also point out that the 4D GPE can be motivated as the proper description of interacting bosons in a hypothetical 4D universe, and so is an interesting theoretical model in its own right. As is well known, the use of the GPE to describe a system of interacting bosons relies on taking the Hartree-Fock approximation and replacing the interaction potential by a contact (Dirac delta) potential. The latter trick is in turn justified by looking at the low energy limit of the solutions for two-particle scattering. In this limit, the solutions are spherically symmetric (s-wave) and correspond to solutions for a contact interaction with the same scattering length as the original potential. While this argument is usually applied only in three dimensions and below, it has also been generalised to arbitrary dimensions [48][49][50], showing that the dimensionality only affects the contact interaction strength, and the form of the short-range singularities that must be removed from the scattering equation. The interaction strength can be considered arbitrary due to scale invariance of the GPE in the absence of an external potential, and the singularities have no effect on the GPE. Hence, it can be concluded that the GPE should be a valid mean-field description of interacting bosons at low energy in 4D.

B.1 Simple Rotations
As described in Section III, we expect that a simple rotation should be able to stabilise a single vortex plane, extending the concept of 2D point vortices and 3D line vortices straightforwardly to four-dimensional systems. Assuming the rotation is in plane 1 (as defined in Section II), this would correspond to a condensate wavefunction of the form: with f (0, r 2 ) = 0, and such that this wave-function approximately takes the form ψ ∝ (x+iy) near the vortex core.
We have verified this minimal vortex structure numerically by performing imaginary time evolution on the full 4D GPE under simple rotation in the plane orthogonal to the expected vortex core (i.e. [Eq (9)] with Ω 1 = 0, Ω 2 = 0). The corresponding density and phase profiles for the numerical stationary state are shown for selected 2D cuts in Fig A1. Here, the initial state was chosen as detailed in Section II and the rotation frequency was chosen as 2Ω 2D crit . These numerical calculations were performed within a discretized 4D hypersphere of radius 8.25ξ, and with resolution 0.5ξ.
As can be seen in Fig A1, the observed density and phase profiles are in good agreement with the single vortex plane [Eq A1]. In particular, the density is depleted for the plane defined by z = 0 and w = 0, as is expected for a single vortex plane that approximately takes the form ψ ∝ (x + iy) near the vortex core. Depending on the 2D cut, this vortex plane either appears as a point [see (e) and (g)], as a line [see (a) and (c)] or as a plane [not shown]. Furthermore, around the vortex plane, the superfluid rotates, as can be seen from the winding of the phase in panels (f) and (h) and from the phase jumps in (b) and (d).

B.2 Double Rotations
As we have shown, the double rotation of a 4D superfluid can stabilise a new type of vortex configuration consisting of two vortex planes intersecting at a point. In Fig A2, we plot the density and phase profiles for additional 2D cuts of the numerical stationary state presented in Fig 3. As can be seen, these profiles have a much richer structure as compared to the case of a single vortex plane shown in Fig A1, as the phase winds simultaneously around both vortex cores with two independent winding numbers. This is also in contrast to 3D systems where two vortex lines may intersect and reconnect over time, but a pair of intersecting vortices is not stabilised by rotation as a stationary state of the system.

B.3 Cuts of the radial profile
As discussed in Section III and shown in Fig 2 (b), we have numerically verified for the solution of the radial equation [Eq (11)] that far from the intersection point of the vortex planes the corresponding density profile is well approximated by a product state of the 2D vortex profiles. To visualise this in an alternative way, we have plotted in Fig A3 (a) cuts of Fig 2 (a) for specific values of r 2 , and then rescaled these by f 1 (r 2 ) in Fig A3 (b). As shown the rescaled curves approach f 1 (r 1 ) for large values of r 2 , verifying the approximation as expected.

B.4 Energy calculation for two intersecting vortex planes in a 4D superfluid
Here, we numerically verify [Eq (12)], which predicts that the energy cost of two intersecting and completely orthogonal vortex planes in a 4D superfluid can be decomposed as a sum of the individual kinetic energies associated with each vortex plane in isolation.
Firstly, we used the numerical solution of the 4D radial density profile presented in Fig 2 to calculate the energy of the intersecting vortex planes as a function of system size in each plane. We then produced a fit of this energy to the functional form of [Eq (12)], with the coefficient of R j /ξ inside the logarithm as the fitting parameter. From this we obtained 2.06 which is very close to the known coefficient of 2.07 (in our units) within the logarithmic form of the vortex energy in 2D and 3D [59]. This shows that the energy of our numerical solution for the radial equation is consistent with being a sum of two individual vortex energies.
Secondly, we performed further simulations on a Cartesian 4D grid, with the same parameters as Fig 3, except for the convergence accuracy which was chosen to be 10 −10 to speed up calculations. We repeated these calculations for different values of Ω ≡ Ω 1 = Ω 2 , ranging between two and three times Ω 2D crit , in order to numerically verify the expected dependence of the energy on the rotation frequency. Here we used three different initial states: one with no phase winding, one with "simple" winding in one plane, and one with "double" winding in two planes. The resulting values for E and µ as a function of Ω are shown in Fig A4, given in units of µ 0 (the chemical potential of a homogeneous state with no vortices or hard-walls but the same number of particles). We obtain straight lines for each of these data series, showing that each state has well defined angular momentum.
For the case with no phase winding, we find that E/µ 0 N = 0.931 and µ/µ 0 = 1.622 are constants which do not depend on frequency, as expected; this data series is therefore plotted with a straight line joining the dots as a guide to the eye. For the double winding case, we have performed a linear fit, obtaining E/µ 0 N = 1.119 − 0.083Ω/Ω 2D crit and µ/µ 0 = 1.822 − 0.083Ω/Ω 2D crit . The gradient, −0.083, is equal to −2Ω 2D crit /µ 0 , meaning that this is the expected gradient of −2 corresponding to particles having one unit of angular momentum in each plane of rotation. For the simple winding case, we fix the gradient to be half that of the double winding line, since this state has angular momentum in only one of the two planes, and perform a linear fit with only the y intercept as a free parameter. We then obtain E/µ 0 N = 1.023 and µ/µ 0 = 1.722 when Ω = 0. This gives an energy cost of 0.188 = 1.119 − 0.931 for the intersecting vortex planes and 0.092 = 1.023 − 0.931 for the single plane, as compared to the state with no vortices. We expect from [Eq (12)] that these energy costs should be related by a simple factor of two for this geometry, and indeed we find numerically that 0.188 − 2 × 0.092 0. The parameters and discretization are detailed in Sec. III. This discretization is reflected in the pixelation, particularly at the boundaries of the plots. The observed density and phase profiles are in good agreement with our numerical ansatz [Eq (10)], which approximately takes the form ψ ∝ (x + iy)(z + iw) near the vortex cores.

Appendix C HOMOTOPY THEORY FOR 4D VORTEX PLANES
Topological excitations, such as vortices, are characterised by topological invariants through homotopy theory. In this approach, the set of allowed topological charges for a given topological defect is given by the set of homotopy classes of maps from a region enclosing the defect to the order parameter manifold. Furthermore, the associated group structure of this set determines the rules for combining two such defects into one.
In 4D, a plane is enclosed by a circle, just like a line in 3D, or a point in 2D, such that the corresponding homotopy group (for a complex order parameter) is  Fig 2 (a), given by fixed values of r2. (b) As in (a) but rescaled by the 2D vortex density profile f1(r2); note the convergence to f1(r1) for large values of r2, showing that the order parameter can be well approximated as a product of 2D vortex profiles away from the intersection point. Close to the intersection point, this approximation breaks down, as can be seen from the deviation between these rescaled cuts in this region.   (9)] with different initial phase profiles. The lines correspond to fits and guides to the eye, respectively, as detailed in the text. The gradient and intercept of these lines give the angular momentum and energy at zero frequency, respectively, of each state, which agree with expected behaviour. π 1 (S 1 ) = Z. This group is the same as for vortices in lower dimensions, and tells us that each vortex has an integer winding number, and that when two vortices combine their winding numbers combine additively. For the case of two intersecting vortex planes the enclosing region is a 2D torus, such as the product of a circle in the xy plane and another circle in the zw plane. The corresponding homotopy group is therefore given by the set of homotopy classes of maps from S 1 × S 1 to S 1 , which is isomorphic to Z × Z [63]. This simply means that each vortex plane has its own winding number, and the two are independent, as expected for two vortices.
Note that this topological classification is the same as for a pair of linked vortex lines in 3D, which can also be enclosed by a torus. The configuration of 4D intersecting planes therefore offers a simple way to realise the homotopy classification of linked vortex lines within the ground state of a simple 4D GPE model. In the future, it would be interesting to generalise this model to more complicated order parameters, such as those realised in the various phases of spinor BECs [61], as then the homotopy group would gain a richer structure, as has been studied in the context of linked line defects in liquid crystals [62].