Vortex dynamics on a M\"obius strip

We consider the dynamics of a two-dimensional incompressible perfect fluid on a M\"obius strip embedded in $\mathbb{R}^3$. The vorticity-streamfunction formulation of the Euler equations is derived from an exterior-calculus form of the momentum equation. The non-orientability of the M\"obius strip and the distinction between forms and pseudo-forms this introduces lead to unusual properties: a boundary condition is provided by the conservation of circulation along the single boundary of the strip, and there is no integral conservation for the vorticity or for any odd function thereof. A finite-difference numerical implementation is used to illustrate the M\"obius-strip realisation of familiar phenomena: translation of vortices along boundaries, shear instability, and decaying turbulence.


Introduction
In this paper, we examine the dynamics of a two-dimensional incompressible inviscid fluid confined to a Möbius strip. The main motivation is curiosity about the impact of the nonorientability of the Möbius strip on the phenomenology of vortex dynamics. A secondary, more frivolous motivation is the wish to produce interesting fluid-dynamical pictures. The system is physically realisable: soap films, for example, can be used as proxys for two-dimensional fluids (e.g. Couder & Basdevant, 1986;Couder et al., 1989); with a suitable wire frame they take the topology of a Möbius strip (e.g. Courant, 1940;Goldstein et al., 2010). The experiments in this paper are however exclusively numerical.
The Möbius strip is a topological object but fluid dynamics requires a geometry, which we need to choose. The simplest geometry is that of a flat strip, imposing antisymmetry conditions across a 'seam' to achieve the necessary twist. However, the flat Möbius strip cannot be embedded in R 3 , a necessary condition for physical realisability. We therefore have chosen another widely used model, the ruled surface constructed by rotating the midpoint of a segment along a horizontal circle while rotating the segment in the vertical at half the rate. This surface, which we denote by M, is parameterised as with −a ≤ ζ ≤ a and 0 ≤ θ ≤ π (or rather θ interpreted as taken modulo π). The property x(ζ, θ + π) = x(−ζ, θ) reflects the twist of the Möbius-strip topology. The parameterisation and shape of M are illustrated in Fig. 1. The mean curvature of M does not vanish, so it is not a minimal surface as formed by a soap film. Its Gaussian curvature does not vanish either, so it is not developable and cannot be made of an inextensible material. Closed-form expressions Figure 1: Möbius strip in R 3 and its parameter space for a = 0.25. The Möbius strip is coloured according to the value of |g| 1/2 , with |g| the metric determinant in (4), which ranges from about 1.5 (dark blue) to 2.5 (light yellow) and whose variations are associated with curvature effects.
are available for Möbius strips in R 3 that are either minimal (Odehnal, 2016) or developable (see Schwarz, 1990) and could replace (1) to gain realism. The equilibrium shape taken by an inextensible material twisted as a Möbius strip can also be obtained, but this requires solving an energy-minimisation problem numerically (Starostin & van der Heijden, 2007). The non-orientability of M is manifested in the fact that the basis (∂ ζ x, ∂ θ x) of its tangent space has opposite orientations for θ = 0 and θ = π. Equivalently, if a normal vector N is added to ∂ ζ x and ∂ θ x to make up a right-handed basis of R 3 , it satisfies N (θ = π) = −N (θ = 0). Non-orientability is a global property, not discernible at the level of a local patch of the strip; incompressible fluid dynamics is also global, through the pressure field. Nonetheless, for narrow strips such as the one in Fig. 1, and with the scale of fluid structures set by the strip's width, we can expect the non-orientability to have a limited impact on the dynamics, which differs from that in a flat cylinder (channel) primarily because of the strip's curvature. The main effect of non-orientability arises from the fact that the sense of rotation of vortices is reversed when they are transported once along the entire length of the strip, that is, when θ increases by π. We will illustrate this by considering the dynamics of a single vortex propagating along the strip's boundary. A mathematical consequence is that vorticity cannot be meaningfully integrated over the Möbius strip; as a result, the familiar connection between circulation and vorticity integral is lost, and the conservation of integrals of functions of the vorticity (Casimirs) is restricted to even functions.
The paper is structured as follows. In §2, we derive the vorticity-streamfunction formulation of the Euler equations for an incompressible perfect fluid on M in terms of the coordinates (ζ, θ). We start with a convenient formulation of the Euler equations that highlights the dual role of the vector and 1-form fields associated with velocity, and treats vorticity as a 2-form. We carefully distinguish between forms (such as the vorticity 2-form) and pseudo-forms (such as the area 2-form, the streamfunction and a pseudo-scalar version of the vorticity which we term 'vorticity density'). The distinction is crucial for a non-orientable surface such as M. We derive the boundary conditions for the vorticity-streamfunction formulation; these include a condition involving the circulation along the (single) boundary of M. In §3 we describe a numerical implementation of the equations which we employ to examine the propagation of a vortex along the strip's boundary, shear instability and decaying turbulence. Section 4 summarises the main features of vortex dynamics on a non-orientable surface. Two appendices provide mathematical details.

Euler equations on the Möbius strip
A coordinate expression for the Euler equations governing the dynamics of an incompressible perfect fluid on an arbitrary surface was derived by Zermelo (1902) (see Zermelo (2013) for a commented English translation). Modern formulations with a focus on point vortices and closed surfaces can be found in, e.g., Kimura (1999), Boatto & Koiller (2015), Dritschel & Boatto (2015), Ragazzo & de Barros Viglioni (2017) or Gustafsson (2019). While Zermelo's expressions can be readily particularised to the Möbius strip, it is instructive to start from a general, coordinate-free version of the Euler equations. A convenient form is Here u is the velocity (vector) field, ν is the momentum 1-form, obtained from u and the metric g as ν = g(u, ·), L u is the Lie derivative and d the exterior derivative (see e.g. Schutz, 1980;Arnold & Khesin, 1998;Besse & Frisch, 2017;Gilbert & Vanneste, 2018). The metric on the Möbius strip M is found by pulling back the Euclidean metric dx ⊗ dx + dy ⊗ dy + dz ⊗ dz of R 3 , that is, by using (1) to express dx, dy and dz in terms of dζ and dθ to find This is a diagonal tensor, indicating that the coordinates (ζ, θ) are orthogonal (i.e. the basis vectors ∂ ζ and ∂ θ are orthogonal) with determinant |g| = 4 + 8 cos θ ζ + (3 + 2 cos(2θ)) ζ 2 .
It follows from (3) that the vectors ∂ ζ and ∂ θ have norms |∂ ζ | = 1 and |∂ θ | = |g| 1/2 . Writing the velocity field in components as the corresponding momentum 1-form is found from ν = g(u, ·) as Introducing (6) into (2a) with L u = u ζ ∂ ζ + u θ ∂ θ and using the commutation L u d = dL u readily gives the momentum equation in terms of u ζ and u θ . This is worked out in Appendix A. The divergence in (2b) is defined by the relation div u µ = d(u µ) = L u µ, where µ is the volume form associated with the metric g (area form since M is two-dimensional) and denotes the interior product or contraction. The area form is given by This is a pseudo-2-form, which changes sign when the orientation of the basis (∂ ζ x, ∂ θ x) changes, as indicated by the notation ±. Since M is non-orientable, a change in the orientation of the basis is unavoidable; it occurs across θ = 0 mod π in the parameterisation (1). It can be verified that (7) is the projection of the Euclidean volume dx ∧ dy ∧ dz on the normal N = (N 1 , N 2 , N 3 ) to the Möbius strip, that is, and apply d to obtain the incompressibility condition in the explicit form The momentum and incompressibility equations (35) and (9) can be conveniently replaced by the vorticity equation. This is best derived directly from (2a): applying d and using that d and L u commute and that d 2 = 0 gives showing that the vorticity 2-form dν is transported by the flow. A vorticity density ω can be defined by ωµ = dν.
This is a pseudo-scalar, which, like µ, changes sign with the orientation of the basis on M.
In particular, ψ(ζ, π, t) = −ψ(−ζ, 0, t), which provides a first boundary condition for (17). The no-normal-flow condition and (14) imply that ψ depends only on t along the boundary of the strip. Because of the pseudo-scalar nature of ψ, this constancy translates into This provides the second boundary condition for the Poisson equation (17). To determine C(t), we need to revert to the momentum formulation (2). It implies conservation of the circulation along the (single) boundary, given by The equation dΓ dt = 0 completes the formulation. The description of our numerical implementation in §3 makes it clear how (23) determines C(t) for use in (21). The existence of infinitely many conserved integrals is a key feature of two-dimensional inviscid incompressible fluids. Here the non-orientablity of the Möbius strip leads to important differences compared with the familiar, orientable case. This is because, while pseudo-2-forms can be integrated over non-orientable surfaces, proper 2-forms cannot -the result of such an integration would not be coordinate independent (e.g. Frankel, 2004). Thus while the area of the Möbius strip M µ is well defined, the integral M dν of the vorticity 2-form is not (and cannot be related to the circulation Γ). As a result, conserved Casimirs functions can only be defined for arbitrary even functions f . The evenness of f is required to ensure that the change of sign of ω which accompanies a change of basis orientation has no effect on f (ω), so that f (ω)µ is a pseudo-2-form. The energy is also well defined and conserved.

Numerical simulations 3.1 Numerical implementation
We solve the vorticity-streamfunction equations (16)-(17) numerically. To control the generation of fine scales that inevitably accompanies non-trivial flows, we modify (16), adding a small-scale dissipation to obtain with ε > 0 small. This requires an additional boundary condition, taken as ω = 0 on the strip boundary ∂M, corresponding to a no-stress condition. With dissipation, the circulation Γ along ∂M changes. We derive the form of dΓ/dt in Appendix B where it is given as (41). Note that the dissipation model chosen is not standard Newtonian friction: the complete Navier-Stokes vorticity equation involves geometric terms omitted from the right-hand side of (26) (see Appendix B, Gilbert et al., 2014;Gilbert & Vanneste, 2021). In the simulations of §3.2, ε is taken small enough the results are representative of the limit ε → 0 except for scales close to the grid size.
To solve (17) and (26) numerically, we discretise ω and ψ on a regular grid in (ζ, θ) coordinates. The Laplacian ∆ is discretised using a straightforward 5-point stencil, accounting for the condition (19) at the seam θ = 0 and assuming homogeneous conditions for ζ = ±a corresponding to the strip's boundary ∂M. The resulting matrix is sparse and has an interesting structure displayed in figure 2, with the standard band structure modified by anti-diagonal blocks in the top-right and bottom-left corners. These result from the twist of condition (19). The boundary condition (21) is implemented as follows. We first compute a numerical approximation to the (time-independent) harmonic function ψ har solving ∆ψ har = 0 with ψ har (ζ = ±a, θ) = ±1, then, at each time step, solve the Poisson equation (17) with homogeneous boundary condition We compute the circulations Γ har and Γ hom (t) associated with each solution by direct evaluation of (22), and we obtain the total circulation Γ(t) from the differential equation (41). The solution of the Poisson equation satisfying the boundary condition (21) is then calculated as with C(t) given by to ensure that the circulation associated with ψ is Γ(t) and (21) holds. We use a splitting scheme to advance (26) in time, with a three-step Adams-Bashforth scheme for the nonlinear advection relying on the Arakawa (1966) discretisation of the Jacobian, Figure 3: Propagation of a vortex along the boundary of the Möbius strip: the vortex travels counterclockwise along the entire boundary, with a vorticity vector field ωN (shown by the arrow in the centre of the vortex) reversing direction as the along-strip coordinate θ increases by π (see also movie 1). and a backward Euler scheme for the viscous diffusion. For the simulations reported below, the half-width of the Möbius strip is set to a = 0.25; the spatial discretisation uses a 200 × 1000 grid in (ζ, θ); the time step is 0.2 and the dissipation parameter is = 4 × 10 −7 . With this small value, energy is approximately conserved over the duration of the simulations presented next.

Simulations
We start by simulating the dynamics of a single vortex patch initialised near the boundary of the Möbius strip. The vortex is expected to travel along the boundary, forming a dipole with its notional image as is familiar in the planar case. The Möbius strip case is intriguing in that the strip has a single boundary, and the direction of rotation of the fluid can be reversed by transporting the vortex once along the strip. The motion of the vortex is in fact straightforward: the vortex travels continuously along the entire boundary. When the along-strip coordinate θ of the vortex has increased by π, so that the vortex returns to the same segment θ = const, the sign of ω is reversed but the vortex initially near the boundary at ζ = a is then near ζ = −a and the propagation direction is unchanged. This is illustrated in the snapshots shown in figure 3 and (better) in the accompanying movie 1. These visualise the vorticity 2-form dν which, unlike ω, is a true geometric object, independent of the choice of coordinate basis. The field dν can be represented as the (true) vector field ωN , normal to M and defined intrinsically as dual to dν via the R 3 volume-form dx ∧ dy ∧ dz. The visualisation in figure 3 and movie 1 shows ωN where |ω| is maximum, that is, at the centre of the vortex, as well as the scalar field |ω| on the strip. The propagation of a vortex along the edge of the Möbius strip makes the single-sidedness of the strip manifest and can be viewed as a fluid-dynamical equivalent of Escher's famous print Möbius Strip II showing ants crawling along the strip.
We next examine a version of shear instability by simulating the evolution of the initial vorticity density ω = ζ(a 2 − ζ 2 ).
This is not a stationary solution of the Euler equations because the associated streamfunction Figure 4: Shear instability on a Möbius strip: the (pseudo-scalar) vorticity density ω is shown at six successive time intervals, starting with a perturbation of the field (31) (see also movie 2).
depends on θ through the metric terms that appear in the Poisson equation (17). Its evolution is nonetheless very slow compared to that which occurs when a small θ-dependent perturbation is added -we use a perturbation of the form ω ∝ (a 2 − ζ 2 ) sin(8θ). The rapid evolution in the presence of a perturbation can be interpreted as resulting from shear instability. Note the relative complexity of (31) compared with the quadratic vorticity profile that is often employed to demonstrate shear instability in planar flows (a quadratic vorticity has a gradient vanishing at a single point, the inflection point of the velocity profile). This is imposed by the topological constraint that 2-forms, here dν = ωdx ∧ dy, must vanish somewhere on non-orientable surfaces. Figure 4 shows a series of snapshots of the pseudo-scalar vorticity density ω at equal time intervals. For distributed vorticity distributions, this field is easier to visualise than the vector field ωN shown in figure 3, but it has the drawback of an abrupt sign change across the seam θ = 0, visible near the top right in the representation of the Möbius strip in figure 4. We emphasise that this discontinuity is artificial and that the physical fields ν and dν (or equivalently ωN ) and the direction of rotation of the fluid are continuous. The snapshots of figure 4 and the accompanying movie 2 show that shear instability on the Möbius strip evolves very much as in planar channel, with the amplification of the initial wavy perturbation of the constant-vorticity lines leading to the formation of rows of counter-rotating vortices. Our final simulation, with results displayed in figure 5 and movie 3, illustrates the dynamics of decaying turbulence on the Möbius strip. We build a small-scale, random initial vorticity density ω 0 (ζ, θ), enforcing its pseudo-scalar property by first constructing a random function f (ζ, θ) that is 4a-periodic in ζ and 2π-periodic in θ using a standard Fourier series method. The initial vorticity density is then taken as The periodicity of f ensures that ω 0 vanishes at the boundary and satisfies condition (18) at the seam θ = π, as required for a pseudo-scalar. The choice of peak wavenumbers 8π/a and 80π for f leads to the reasonably isotropic, small-scale ω 0 shown in the first panel of figure 5.
In addition to an initial vorticity density, we also impose an non-zero (weak) initial circulation which leads to an overall drift of the vortices best seen in movie 3. The evolution shown in figure 5 and movie 3 is characteristic of an inverse energy cascade, with vortex mergers that are local and hence unaffected by the topology of the strip. Over time scales much larger than those explored in our simulation, the scale of the vorticity field may become such that the evolution is more strongly influenced by the topology of M.

Summary
This paper examines how the Euler equations governing the dynamics of two-dimensional inviscid fluids can be formulated on a non-orientable surface, using a Möbius strip embedded in R 3 as a concrete example. The velocity vector field u, the corresponding momentum 1-form ν and the pressure field p are all geometrically intrinsic objects, which are independent of the choice of coordinate basis regardless of whether the surface is orientable or not. Pseudo-fields, which depend on coordinates through the orientation of the basis, only appear when the Euler equations are simplified to their vorticity-streamfunction formulation. While the vorticity dν is a proper 2-form, the more convenient ω is a pseudo-scalar because its definition via dν = ωµ involves the area pseudo-2-form µ. Similarly, the streamfunction ψ is defined via the relation dψ = −u µ and is a pseudo-scalar. With our choice of coordinates (ζ, θ), the difference between pseudo-scalars and scalars is simply that the former satisfy the twist condition (19) across θ = 0 mod π instead of continuity. In the vorticity-streamfunction formulation of the Euler equations, the condition that ψ depends on time only along each connected component of the boundary must be supplemented by additional conditions tracing back to the original momentum formulation. On the Möbius strip, the boundary is simply connected and the single function of time C(t) that needs to be determined is found straighforwardly from conservation of the circulation along the single boundary. The main dynamical effect of non-orientability is that the sign of ω can be changed by transport: a vortex patch rotating clockwise transported once along the strip returns rotating counterclockwise, as demonstrated in the numerical experiment of figure  3. This is why odd functions of the vorticity density are not conserved in an integral sense. Loosely speaking, the number of integral conservation laws on the Möbius strip is half what it is on an orientable surface. It would be interesting to examine the consequences of this for the long-time dynamics of turbulence, e.g. as described by statistical mechanics. It would also be interesting to study purely non-dissipative dynamics on a Möbius strip by considering point vortices or contour dynamics.