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

(2.1a)$$\begin{gather} (\partial_t + \mathcal{L}_{\boldsymbol{u}}) \nu ={-} {\rm d} \left( p - \tfrac{1}{2} | \boldsymbol{u}|^{2} \right), \end{gather}$$

(2.1b)$$\begin{gather}{\mathrm{div}} \boldsymbol{u} = 0. \end{gather}$$

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

(2.2)\begin{equation} g = \textrm{d} \zeta \otimes \textrm{d} \zeta + \left(4 + 8 \cos \theta \zeta + (3+2 \cos(2 \theta)) \zeta^{2}\right) \textrm{d} \theta\otimes \textrm{d} \theta. \end{equation}

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

(2.3)\begin{equation} |g| = 4 + 8 \cos \theta \zeta + \left(3+2 \cos(2 \theta)\right) \zeta^{2}. \end{equation}

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

(2.4)\begin{equation} \boldsymbol{u} = u^{\zeta} \partial_\zeta + u^{\theta} \partial_\theta, \end{equation}

the corresponding momentum 1-form is found from $\nu =g(\boldsymbol {u},\cdot )$ as

(2.5)\begin{equation} \nu = \nu_\zeta \,\textrm{d} \zeta + \nu_\theta \,\textrm{d} \theta = u^{\zeta} \,\textrm{d} \zeta + |g| u^{\theta} \,\textrm{d} \theta. \end{equation}

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

(2.6)\begin{equation} \mu ={\pm} |g|^{1/2} \,\textrm{d} \zeta \wedge \textrm{d} \theta. \end{equation}

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

(2.7)\begin{equation} \boldsymbol{u} \lrcorner \mu ={\pm} |g|^{1/2} (u^{\zeta} \,\textrm{d} \theta - u^{\theta} \,\textrm{d} \zeta), \end{equation}

and apply ${\rm d}$ to obtain the incompressibility condition in the explicit form

(2.8)\begin{equation} {\mathrm{div}} \boldsymbol{u} = {|g|^{{-}1/2}} \left(\partial_\zeta ( |g|^{1/2} u^{\zeta}) + \partial_\theta (|g|^{1/2} u^{\theta}) \right) = 0. \end{equation}

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

(2.9)\begin{equation} (\partial_t + \mathcal{L}_{\boldsymbol{u}}) \,\textrm{d} \nu = 0, \end{equation}

showing that the vorticity 2-form $\textrm {d} \nu$ is transported by the flow. A vorticity density $\omega$ can be defined by

(2.10)\begin{equation} \omega \mu = \textrm{d} \nu. \end{equation}

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

(2.11)\begin{equation} \omega ={\pm} |g|^{{-}1/2} \left(\partial_\zeta (|g| u^{\theta}) - \partial_\theta u^{\zeta} \right). \end{equation}

Incompressibility ensures that $(\partial _t + \mathcal {L}_{\boldsymbol {u}}) \mu = 0$, hence $\omega$ is also transported by the flow:

(2.12)\begin{equation} (\partial_t + \mathcal{L}_{\boldsymbol{u}}) \omega =0. \end{equation}

The velocity $\boldsymbol {u}$ can be related to the vorticity by means of a stream function $\psi$, with

(2.13a,b)\begin{equation} u^{\zeta} ={-} |g|^{{-}1/2} \partial_\theta \psi \quad \text{and}\quad u^{\theta} = |g|^{{-}1/2} \partial_\zeta \psi, \end{equation}

so as to satisfy the incompressibility condition (2.8). Equivalently, this can be written as

(2.14)\begin{equation} \textrm{d} \psi ={-} \boldsymbol{u} \lrcorner \mu . \end{equation}

With (2.13a,b), (2.12) becomes

(2.15)\begin{equation} \partial_t \omega + |g|^{{-}1/2} \partial(\psi,\omega)=0, \end{equation}

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

(2.16)\begin{equation} \omega = {\rm \Delta} \psi = {|g|^{{-}1/2}} \left(\partial_\zeta \left(|g|^{1/2} \partial_\zeta \psi\right) + \partial_\theta \left(|g|^{{-}1/2} \partial_\theta \psi \right) \right), \end{equation}

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

(2.17a,b)\begin{equation} \omega(\zeta,\theta+{\rm \pi},t) ={-}\omega(-\zeta,\theta,t) \quad \text{and}\quad \psi(\zeta,\theta+{\rm \pi},t) ={-}\psi(-\zeta,\theta,t). \end{equation}

In particular,

(2.18)\begin{equation} \psi(\zeta,{\rm \pi},t)={-}\psi(-\zeta,0,t), \end{equation}

which provides a first boundary condition for (2.16). The no-normal-flow condition

(2.19)\begin{equation} u^{\zeta}(\zeta={\pm} a,\theta,t)=0 \end{equation}

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

(2.20)\begin{equation} \psi (\zeta ={\pm} a, \theta,t)={\pm} C(t). \end{equation}

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

(2.21)\begin{align} \varGamma &= \int_{\partial \mathcal{M}} \nu = \int_0^{\rm \pi} \left(\nu_\theta(a,\theta,t)+\nu_\theta({-}a,\theta,t)\right) \textrm{d} \theta \nonumber\\ &= \int_0^{\rm \pi} \left({|g(a,\theta)|^{1/2}} \partial_\zeta \psi (a,\theta,t)+{|g({-}a,\theta)|^{1/2}} \partial_\zeta \psi ({-}a,\theta,t)\right) \textrm{d} \theta. \end{align}

The equation

(2.22)\begin{equation} \frac{\mathrm{d} \varGamma}{\mathrm{d} t}=0 \end{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

(2.23)\begin{equation} \mathcal{C}_f = \int_{\mathcal{M}} f(\omega) \mu = \int_{{-}a}^{a} \int_0^{\rm \pi} f(\omega) |g|^{1/2} \,\textrm{d} \zeta \,\textrm{d} \theta, \end{equation}

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

(2.24)\begin{equation} \mathcal{E} ={-} \tfrac{1}{2} \int_M \psi \omega \mu ={-} \tfrac{1}{2} \int_{{-}a}^{a} \int_0^{\rm \pi} \psi \omega |g|^{1/2} \,\textrm{d} \zeta\,\textrm{d} \theta \end{equation}

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

(3.1)\begin{equation} \partial_t \omega + |g|^{{-}1/2} \partial(\psi,\omega) = \varepsilon {\rm \Delta} \omega, \end{equation}

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

(3.2)\begin{equation} {\rm \Delta} \psi_{har} = 0 \quad \textrm{with}\ \psi_{har}(\zeta={\pm} a,\theta)={\pm} 1, \end{equation}

then, at each time step, solve the Poisson equation (2.16) with homogeneous boundary condition

(3.3)\begin{equation} {\rm \Delta} \psi_{hom} = \omega \quad \textrm{with}\ \psi_{hom}(\zeta={\pm} a,\theta,t)=0. \end{equation}

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

(3.4)\begin{equation} \psi = \psi_{hom} + C(t) \psi_{har}, \end{equation}

with $C(t)$ given by

(3.5)\begin{equation} C(t) = \left( \varGamma(t) - \varGamma_{hom}(t) \right)/\varGamma_{har} \end{equation}

to ensure that the circulation associated with $\psi$ is $\varGamma (t)$, and (2.20) holds.

Figure 2. Visualisation of the non-zero entries in the matrix discretising the Laplacian operator on the Möbius strip. The unusual antidiagonal bands in the top-right and bottom-left corners result from the twist of the Möbius strip.

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.

Figure 3. Propagation of a vortex along the boundary of the Möbius strip: the vortex travels anticlockwise along the entire boundary, with a vorticity vector field $\omega \boldsymbol {N}$ (shown by the arrow in the centre of the vortex) reversing direction as the along-strip coordinate $\theta$ increases by ${\rm \pi}$ (see also supplementary movie 1 available at https://doi.org/10.1017/jfm.2021.581).

We next examine a version of shear instability by simulating the evolution of the initial vorticity density

(3.6)\begin{equation} \omega=\zeta (a^{2}-\zeta^{2}). \end{equation}

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.

Figure 4. Shear instability on a Möbius strip: the (pseudo-scalar) vorticity density $\omega$ is shown at six successive time intervals, starting with a perturbation of the field (3.6) (see also movie 2).

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

(3.7)\begin{equation} \omega_0(\zeta,\theta) = f(\zeta,\theta)-f(2a-\zeta,\theta)+f(2a+\zeta,\theta+{\rm \pi})-f(-\zeta,\theta+{\rm \pi}). \end{equation}

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}$.

Figure 5. Decaying turbulence on a Möbius strip: the (pseudo-scalar) vorticity density $\omega$ is shown at six successive equal time intervals (see also movie 3).

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.