1. 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 a curiosity about the impact of the non-orientability 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 proxies for two-dimensional fluids (e.g. Couder & Basdevant Reference Couder and Basdevant1986; Couder, Chomaz & Rabaud Reference Couder, Chomaz and Rabaud1989); with a suitable wire frame they take the topology of a Möbius strip (e.g. Courant Reference Courant1940; Goldstein et al. Reference Goldstein, Moffatt, Pesci and Ricca2010). 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 $\mathbb {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 $\mathcal {M}$, is parameterised as
with $-a \le \zeta \le a$ and $0 \le \theta \le {\rm \pi}$ (or, rather, $\theta$ interpreted as taken modulo ${\rm \pi}$). The property $\boldsymbol {x}(\zeta ,\theta +{\rm \pi} )=\boldsymbol {x}(-\zeta ,\theta )$ reflects the twist of the Möbius-strip topology. The parameterisation and shape of $\mathcal {M}$ are illustrated in figure 1. The mean curvature of $\mathcal {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 are available for Möbius strips in $\mathbb {R}^{3}$ that are either minimal (Odehnal Reference Odehnal2016) or developable (see Schwarz Reference Schwarz1990) and could replace (1.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 Reference Starostin and van der Heijden2007).
The non-orientability of $\mathcal {M}$ is manifested in the fact that the basis $(\partial _{\zeta } \boldsymbol {x},\partial _\theta \boldsymbol {x})$ of its tangent space has opposite orientations for $\theta =0$ and $\theta ={\rm \pi}$. Equivalently, if a normal vector $\boldsymbol {N}$ is added to $\partial _{\zeta } \boldsymbol {x}$ and $\partial _\theta \boldsymbol {x}$ to make up a right-handed basis of $\mathbb {R}^{3}$, it satisfies $\boldsymbol {N}(\theta ={\rm \pi} )=-\boldsymbol {N}(\theta =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 figure 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 $\theta$ increases by ${\rm \pi}$. 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 the 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–stream function formulation of the Euler equations for an incompressible perfect fluid on $\mathcal {M}$ in terms of the coordinates $(\zeta ,\theta )$. 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 stream function and a pseudo-scalar version of the vorticity which we term ‘vorticity density’). The distinction is crucial for a non-orientable surface such as $\mathcal {M}$. We derive the boundary conditions for the vorticity–stream function formulation; these include a condition involving the circulation along the (single) boundary of $\mathcal {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.
2. 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 (Reference Zermelo1902) (see Zermelo (Reference Zermelo2013) for a commented English translation). Modern formulations with a focus on point vortices and closed surfaces can be found in, for example, Kimura (Reference Kimura1999), Boatto & Koiller (Reference Boatto and Koiller2015), Dritschel & Boatto (Reference Dritschel and Boatto2015), Ragazzo & de Barros Viglioni (Reference Ragazzo and de Barros Viglioni2017) or Gustafsson (Reference Gustafsson2019). 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 $\boldsymbol {u}$ is the velocity (vector) field, $\nu$ is the momentum 1-form, obtained from $\boldsymbol {u}$ and the metric $g$ as $\nu = g(\boldsymbol {u},\cdot )$, $\mathcal {L}_{\boldsymbol {u}}$ is the Lie derivative and $\textrm {d}$ the exterior derivative (see e.g. Schutz Reference Schutz1980; Arnold & Khesin Reference Arnold and Khesin1998; Besse & Frisch Reference Besse and Frisch2017; Gilbert & Vanneste Reference Gilbert and Vanneste2018).
The metric on the Möbius strip $\mathcal {M}$ is found by pulling back the Euclidean metric ${\textrm {d}x} \otimes {\textrm {d}x} + {\textrm {d}y} \otimes {\textrm {d}y} + \textrm {d} z \otimes \textrm {d} z$ of $\mathbb {R}^{3}$, that is, by using (1.1) to express ${\textrm {d}x}$, ${\textrm {d}y}$ and $\textrm {d} z$ in terms of $\textrm {d} \zeta$ and $\textrm {d} \theta$ to find
This is a diagonal tensor, indicating that the coordinates $(\zeta ,\theta )$ are orthogonal (i.e. the basis vectors $\partial _\zeta$ and $\partial _\theta$ are orthogonal), with determinant
It follows from (2.2) that the vectors $\partial _\zeta$ and $\partial _\theta$ have norms $|\partial _\zeta | =1$ and $|\partial _\theta |=|g|^{1/2}$. Writing the velocity field in components as
the corresponding momentum 1-form is found from $\nu =g(\boldsymbol {u},\cdot )$ as
Introducing (2.5) into (2.1a) with $\mathcal {L}_{\boldsymbol {u}} = u^{\zeta } \partial _\zeta + u^{\theta } \partial _\theta$ and using the commutation $\mathcal {L}_{\boldsymbol {u}} \,\textrm {d} = \textrm {d} \mathcal {L}_{\boldsymbol {u}}$ readily gives the momentum equation in terms of $u^{\zeta }$ and $u^{\theta }$. This is worked out in Appendix A.
The divergence in (2.1b) is defined by the relation ${\mathrm {div}} \boldsymbol {u} \,\mu = \textrm {d} (\boldsymbol {u} \lrcorner \mu ) = \mathcal {L}_{\boldsymbol {u}} \mu$, where $\mu$ is the volume form associated with the metric $g$ (here area form since $\mathcal {M}$ is two-dimensional) and $\lrcorner$ 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 $(\partial_\zeta,\partial_\theta)$ changes, as indicated by the notation $\pm$. Since $\mathcal {M}$ is non-orientable, a change in the orientation of the basis is unavoidable; it occurs across $\theta = 0 \mod {\rm \pi}$ in the parameterisation (1.1). It can be verified that (2.6) is the projection of the Euclidean volume ${\textrm {d}x} \wedge {\textrm {d}y} \wedge \textrm {d} z$ on the normal $\boldsymbol {N}=(N^{1},N^{2},N^{3})$ to the Möbius strip, that is, $\mu = {\boldsymbol {N}} \lrcorner \, {\textrm {d}x} \wedge {\textrm {d}y} \wedge \textrm {d} z = N^{1} {\textrm {d}y} \wedge \textrm {d} z + N^{2} \,\textrm {d} z \wedge \mathrm {d} x + N^{3} \,\mathrm {d} x \wedge {\textrm {d}y}$. With (2.6) we compute
and apply ${\rm d}$ to obtain the incompressibility condition in the explicit form
The momentum and incompressibility equations (A3) and (2.8) can be conveniently replaced by the vorticity equation. This is best derived directly from (2.1a): applying ${\rm d}$ and using that ${\rm d}$ and $\mathcal {L}_{\boldsymbol {u}}$ commute and that ${\rm d}^{2}=0$ gives
showing that the vorticity 2-form $\textrm {d} \nu$ is transported by the flow. A vorticity density $\omega$ can be defined by
This is a pseudo-scalar, which, like $\mu$, changes sign with the orientation of the basis on $\mathcal {M}$. Applying ${\rm d}$ to (2.4) and using (2.6) gives the explicit expression
Incompressibility ensures that $(\partial _t + \mathcal {L}_{\boldsymbol {u}}) \mu = 0$, hence $\omega$ is also transported by the flow:
The velocity $\boldsymbol {u}$ can be related to the vorticity by means of a stream function $\psi$, with
so as to satisfy the incompressibility condition (2.8). Equivalently, this can be written as
With (2.13a,b), (2.12) becomes
where $\partial (\psi ,\omega ) = \partial _\zeta \psi \partial _\theta \omega - \partial _\theta \psi \partial _\zeta \omega$ is the Jacobian operator. Substituting (2.13a,b) into (2.11) then gives the vorticity–stream function relation
where we have identified the Laplacian ${\rm \Delta}$. A completely geometric derivation of (2.13a,b) and (2.16), which avoids coordinates but requires the introduction of the Hodge * operator, is given in Appendix B.
Equations (2.15) and (2.16), together with the form (2.3) of the metric determinant $|g|$ constitute the vorticity–stream function formulation of the Euler equations on the Möbius strip. In the familiar way, the vorticity density $\omega$ is the dynamical variable from which the stream function $\psi$ is derived by inverting the Poisson equation (2.16). This requires boundary conditions. Because $\omega$ and $\psi$ are pseudo-scalars, they satisfy
In particular,
which provides a first boundary condition for (2.16). The no-normal-flow condition
and (2.13a,b) imply that $\psi$ depends only on $t$ along the boundary of the strip. Because of the pseudo-scalar nature of $\psi$, this constancy translates into
This provides the second boundary condition for the Poisson equation (2.16). To determine $C(t)$, we need to revert to the momentum formulation (2.1). It implies conservation of the circulation along the (single) boundary, given by
The equation
completes the formulation. The description of our numerical implementation in § 3 makes it clear how (2.22) determines $C(t)$ for use in (2.20).
The existence of infinitely many conserved integrals is a key feature of two-dimensional inviscid incompressible fluids. Here the non-orientability 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 Reference Frankel2004). Thus, while the area of the Möbius strip $\int _\mathcal {M} \mu$ is well defined, the integral $\int _\mathcal {M} {\rm d} \nu$ of the vorticity 2-form is not (and cannot be related to the circulation $\varGamma$). 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 $\omega$ which accompanies a change of basis orientation has no effect on $f(\omega )$, so that $f(\omega ) \mu$ is a pseudo-2-form. The energy
is also well defined and conserved.
3. Numerical simulations
3.1. Numerical implementation
We solve the vorticity–stream function equations (2.15) and (2.16) numerically. To control the generation of fine scales that inevitably accompanies non-trivial flows, we modify (2.15), adding a small-scale dissipation to obtain
with $\varepsilon >0$ small. This requires an additional boundary condition, taken as $\omega =0$ on the strip boundary $\partial \mathcal {M}$, corresponding to a no-stress condition. With dissipation, the circulation $\varGamma$ along $\partial \mathcal {M}$ changes. We derive the form of $\mathrm {d} \varGamma /\mathrm {d} t$ in Appendix B where it is given as (B6). 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 (3.1) (see Appendix B, Gilbert, Riedinger & Thuburn (Reference Gilbert, Riedinger and Thuburn2014) and Gilbert & Vanneste (Reference Gilbert and Vanneste2021)). In the simulations of § 3.2, $\varepsilon$ is taken small enough that the results are representative of the limit $\varepsilon \to 0$ except for scales close to the grid size.
To solve (2.16) and (3.1) numerically, we discretise $\omega$ and $\psi$ on a regular grid in $(\zeta ,\theta )$ coordinates. The Laplacian ${\rm \Delta}$ is discretised using a straightforward five-point stencil, accounting for the condition (2.18) at the seam $\theta =0$ and assuming homogeneous conditions for $\zeta =\pm a$ corresponding to the strip's boundary $\partial \mathcal {M}$. The resulting matrix is sparse and has an interesting structure displayed in figure 2, with the standard band structure modified by antidiagonal blocks in the top-right and bottom-left corners. These result from the twist of condition (2.18). The boundary condition (2.20) is implemented as follows. We first compute a numerical approximation to the (time-independent) harmonic function $\psi _{har}$ solving
then, at each time step, solve the Poisson equation (2.16) with homogeneous boundary condition
We compute the circulations $\varGamma _{har}$ and $\varGamma _{hom}(t)$ associated with each solution by direct evaluation of (2.21), and we obtain the total circulation $\varGamma (t)$ from the differential equation (B6). The solution of the Poisson equation satisfying the boundary condition (2.20) is then calculated as
with $C(t)$ given by
to ensure that the circulation associated with $\psi$ is $\varGamma (t)$, and (2.20) holds.
We use a splitting scheme to advance (3.1) in time, with a three-step Adams–Bashforth scheme for the nonlinear advection relying on the Arakawa (Reference Arakawa1966) discretisation of the Jacobian, 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 \times 1000$ grid in $(\zeta ,\theta )$; the time step is $0.2$ and the dissipation parameter is $\epsilon = 4 \times 10^{-7}$. With this small value, energy is approximately conserved over the duration of the simulations presented next.
3.2. 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 $\theta$ of the vortex has increased by ${\rm \pi}$, so that the vortex returns to the same segment $\theta ={const.}$, the sign of $\omega$ is reversed but the vortex initially near the boundary at $\zeta =a$ is then near $\zeta =-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 ${\rm d} \nu$ which, unlike $\omega$, is a true geometric object, independent of the choice of coordinate basis. The field ${\rm d} \nu$ can be represented as the (true) vector field $\omega \boldsymbol {N}$, normal to $\mathcal {M}$ and defined intrinsically as dual to ${\rm d}\nu$ via the $\mathbb {R}^{3}$ volume-form ${\textrm {d}x} \wedge {\textrm {d}y} \wedge \textrm {d} z$. The visualisation in figure 3 and movie 1 shows $\omega \boldsymbol {N}$ where $|\omega |$ is maximum, that is, at the centre of the vortex, as well as the scalar field $|\omega |$ 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
This is not a stationary solution of the Euler equations because the associated stream function depends on $\theta$ through the metric terms that appear in the Poisson equation (2.16). Its evolution is nonetheless very slow compared with that which occurs when a small $\theta$-dependent perturbation is added – we use a perturbation of the form $\omega ' \propto (a^{2}-\zeta ^{2}) \sin (8 \theta )$. The rapid evolution in the presence of a perturbation can be interpreted as resulting from shear instability. Note the relative complexity of (3.6) 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 $\textrm {d} \nu = \omega \,{\textrm {d}x} \wedge {\textrm {d}y}$, must vanish somewhere on non-orientable surfaces. Figure 4 shows a series of snapshots of the pseudo-scalar vorticity density $\omega$ at equal time intervals. For distributed vorticity distributions, this field is easier to visualise than the vector field $\omega \boldsymbol {N}$ shown in figure 3, but it has the drawback of an abrupt sign change across the seam $\theta = 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 $\nu$ and $\textrm {d} \nu$ (or equivalently $\omega \boldsymbol {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 $\omega _0(\zeta ,\theta )$, enforcing its pseudo-scalar property by first constructing a random function $f(\zeta ,\theta )$ that is $4a$-periodic in $\zeta$ and $2{\rm \pi}$-periodic in $\theta$ using a standard Fourier series method. The initial vorticity density is then taken as
The periodicity of $f$ ensures that $\omega _0$ vanishes at the boundary and satisfies condition (2.17a,b) at the seam $\theta ={\rm \pi}$, as required for a pseudo-scalar. The choice of peak wavenumbers $8{\rm \pi} /a$ and $80{\rm \pi}$ for $f$ leads to the reasonably isotropic, small-scale $\omega _0$ shown in the first panel of figure 5. In addition to an initial vorticity density, we also impose a 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 $\mathcal {M}$.
4. 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 $\mathbb {R}^{3}$ as a concrete example. The velocity vector field $\boldsymbol {u}$, the corresponding momentum 1-form $\nu$ 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–stream function formulation. While the vorticity $\textrm {d} \nu$ is a proper 2-form, the more convenient $\omega$ is a pseudo-scalar because its definition via $\textrm {d} \nu = \omega \mu$ involves the area pseudo-2-form $\mu$. Similarly, the stream function $\psi$ is defined via the relation $\textrm {d} \psi = - \boldsymbol {u} \lrcorner \mu$ and is a pseudo-scalar. With our choice of coordinates $(\zeta ,\theta )$, the difference between pseudo-scalars and scalars is simply that the former satisfy the twist condition (2.18) across $\theta = 0 \mod {\rm \pi}$ instead of continuity. In the vorticity–stream function formulation of the Euler equations, the condition that $\psi$ 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 straightforwardly from conservation of the circulation along the single boundary. The main dynamical effect of non-orientability is that the sign of $\omega$ can be changed by transport: a vortex patch rotating clockwise transported once along the strip returns rotating anticlockwise, 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 of what it is on an orientable surface. It would be interesting to examine the consequences of this for the long-time dynamics of turbulence, for example 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.
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2021.581.
Acknowledgements
J.V. thanks S.G. Llewellyn Smith for pointing out useful references.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Momentum equation
We derive the momentum equation in coordinates $(\zeta ,\theta )$. The derivation illustrates the convenience of working directly with forms rather than their components, as advocated by Frankel (Reference Frankel2004). Starting from the coordinate-free momentum equation (2.1a) we compute
and note that
to find the two components
where $\mathcal {L}_{\boldsymbol {u}} = u^{\zeta } \partial _\zeta + u^{\theta } \partial _\theta$ since it acts on scalars.
Appendix B. Geometric derivation
The vorticity–stream function formulation can be derived without resorting to coordinates. Following Kimura (Reference Kimura1999) and Ragazzo & de Barros Viglioni (Reference Ragazzo and de Barros Viglioni2017), this is best carried out using the Hodge * operator which, in two dimensions, is determined by the following:
where $\alpha$ is a 1-form and $\alpha ^{\sharp }$ its dual vector, with $\alpha =g(\alpha ^{\sharp },\cdot )$. In coordinates, we have
Noting that $** \alpha = - \alpha$ (e.g. Frankel Reference Frankel2004), we can rewrite (2.14) and (2.16) as
which identifies the Laplacian acting on (pseudo-)scalars as ${\rm \Delta} = * \textrm {d} * \textrm {d}$.
This form of the Laplacian is also useful to compute the change in the circulation along the boundary that is introduced by dissipation. For simplicity, we take the our dissipative model to be
where ${\rm \Delta} = \textrm {d}*\textrm {d}* + *\textrm {d}*\textrm {d}$ is the Laplace–de Rham operator. We emphasise that this is not the standard (Navier–Stokes) viscous dissipation which, instead, involves a different Laplacian defined via the divergence of the stress tensor (Gilbert et al. Reference Gilbert, Riedinger and Thuburn2014; Gilbert & Vanneste Reference Gilbert and Vanneste2021). Integrating over the (closed) boundary $\partial \mathcal {M}$ and noting that the outer exact differential in ${\rm \Delta}$ integrates to zero gives
leading to the coordinate expression
on using (B2) and that $\omega =0$ on $\partial \mathcal {M}$.