Vorticity‐stream function approaches are inappropriate to solve the surface Navier‐Stokes equation on a torus

The growing interest in interfacial flows in engineering and biological applications demands to numerically solve the incompressible surface (Navier)‐Stokes equation. A common approach, based on a vorticity‐stream function formulation reformulates the equation into a system of scalar‐valued surface partial differential equations, for which various numerical methods have been developed. However, on surfaces, which are not simply connected, this approach is not appropriate. We here explain the underlying situation and provide details and examples.


Introduction
In a recent letter [1] Pearce et al. investigate the turbulent dynamics of a two-dimensional active nematic liquid crystal on the surface of a torus. The underlying model combines an incompressible surface Navier-Stokes equation with friction and active forcing with a surface Landau-de Gennes model for nematic liquid crystals. To solve the surface Navier-Stokes equation a vorticity-stream function approach is considered. The approach is based on a decomposition of the velocity field in divergenceand curl-free parts. While appropriate on simply connected surfaces, this is not sufficient on non-simply connected surfaces, such as the considered torus. As a consequence of the topology also non-trivial harmonic parts, velocity fields which are divergence-and curl-free, exist. They are not represented by the vorticity-stream function approach, see [2] for an example.
where v denotes the velocity, D Dt the covariant material time derivative, ∆ B the Bochner Laplacian, grad S the surface gradient, div S the surface divergence, K the Gaussian curvature, p the surface pressure, Q a Q-tensor, ρ the material density, η the dynamic viscosity, ζ the friction coefficient and α the strength of the external forcing. Eq. (1) is supplemented with the incompressibility condition div S v = 0.
We will consider three cases: Firstly, no external contributions are considered, i. e. ζ = 0 and α = 0. With the exception of some special initial conditions for which the solution converges to zero, any solution of the surface Navier-Stokes equation on a torus converges to a stationary Killing vector field, which contains non-trivial harmonic parts, see [2,3]. Secondly, we investigate damped flow, i. e. ζ = 0 and α = 0. Clearly, the velocity converges to zero. But in contrast to the argumentation in [1] this means that the harmonic and the anti-harmonic part vanish identically over time with a strong influence of the harmonic contributions. This is the same setup as in the Supplementary Information of [1], where it is argued that harmonic contributions are negligible. Thirdly, the external contributions are considered, i. e. ζ = 0 and α = 0. In this case we demonstrate on an example that the harmonic parts do not vanish over time. There is no reason to assume that the harmonic solution components are negligible. In fact neglecting them changes the velocity qualitatively in an nonphysical manner.

Notation
Let Ker div S be the space of divergence free vector fields on the surface S. This space can be decomposed w. r. t. the L 2 inner product into parts with and without curl-components by using the Hodge theorem, i. e.
where the last considers the harmonic parts. Here, we explicitly note that this decomposition does not hold in the sense of the local inner product, but in the sense of the L 2 inner product. In the following we restrict to the surface of the torus with radii ratio ξ := a/b with major and minor radius a and b, respectively, and denote the two basis vector fields by ∂ θ x and ∂ φ x as well as the two basis harmonic vector fields by Θ and Φ, as considered in [1]. It is thereby noted that these harmonic vector fields are locally orthogonal to each other w. r. t. the local inner product. Thus, an orthogonal projection of a vector field into the space of harmonic vector fields Π H : T 1 S → H(S) can be considered using basic orthogonalization techniques, i. e.
where T 1 S denotes the space of vector fields on S. Using this, the anti-harmonic projection Π ⊥ H :

Flow with no external contributions
Let ζ = 0 and α = 0. The solution of eq. (1) converges to a Killing vector field for t → ∞ (except for some special initial data with symmetric behavior which implies a vanishing solution), see [2]. These vector fields exist on rotationally symmetric surfaces and are characterized by a vanishing deformation tensor. In the present case of a torus, the Killing vector fielddenoted by v K -can be determined by v K = α K ∂ φ x with a constant α K and describes a rigidly rotating torus. It can be easily verified that Π H v K = 0 and thus, non-trivial Killing vector fields contain non-trivial harmonic parts, which cannot vanish asymptotically in general. Cancelling out the harmonic part for a rigidly rotating torus would magically generate internal friction and the velocity dissipates to zero without any physical reason.

Damped flow
If we consider external friction, i. e. ζ = 0, and no additional forces, i. e. α = 0, the harmonic part of the solution vanishes exponentially. But the anti-harmonic part do the same and we do not see any reasons to only highlight this behavior for the harmonic part. In particular, it can be verified that the damped Killing solution v K (t) = v K (0) exp (−ζ/ρt) solves eq. (1) with the Killing vector field from above as initial condition. Thus, the harmonic part Π H v K and the anti-harmonic part Π ⊥ H v K vanish exponentially with the same order (in time). Figure 1 shows this behavior, where the kinetic energy of the full damped Killing vector field as well as its harmonic and anti-harmonic part is shown. The harmonic contribution is even larger than the anti-harmonic part. Moreover, the ratio between these parts can be determined by which only depends on the radii ratio ξ of the torus and is monotonically increasing for ξ > 1. Furthermore, ξ > ((7 + 2 · 19 1/2 )/6) 1/2 ≈ 1.619 yields r H > 1 and thus larger harmonic parts.

Damped and forced flow
Considering the full equation (1) with ζ = 0 and α = 0. The force div S Q feeds harmonic parts into the solution, since div S Q has harmonic contributions itself in general. To highlight this behavior, we construct an appropriate Q-tensor. Generally, a Q-tensor can be obtained by Q = grad S V, where V ∈ H(S). This holds as tr grad S V = div S V = 0 and grad S V, ε = − rot S V = 0 with the Levi-Civita tensor ε. Thus, grad S V is a Q-tensor and we obtain div S Q = KV. For the latter identity we used the Weizenböck identity ∆ B v = ∆ dR v + Kv with the Laplace-deRham operator ∆ dR , which can be obtained by [2]. Consider the special choice Q = −b 2 (ξ 2 − 1) grad S Φ. Thus, we obtain div S Q = −b 2 (ξ 2 − 1)KΦ and the resulting harmonic force is given by Π H div S Q = Φ. This means we constantly feed the system with the harmonic field αΦ such that the harmonic part of the solution cannot vanish. In this situation we expect a balance of internal and external friction with the applied force, such that the solution converges to a stationary nontrivial vector field for t → ∞. We further expect that the reached steady state solution does not have negligible harmonic contributions. To demonstrate this we discretize eq. (1) in time by using an implicit Euler scheme and in space by using a symmetric ansatz and finite differences in poloidal direction. Figure 2 shows the kinetic energy of the full solution as well as the harmonic and anti-harmonic components. Indeed the harmonic contributions are not negligible and are even larger than the anti-harmonic parts in this simple example. The harmonic-anti-harmonic ratio r H | t=60 ≈ 1.374 in the steady state regime also reflects this behavior.

Wrong assumption
In the Supplementary Information of [1] the authors restrict to the case ζ = 0, α = 0 and argue that the harmonic solution components decay to zero due to the friction term in the surface Navier-Stokes equation. As discussed above, this is also true for the anti-harmonic part. However, the argumentation in [1] is based on eq. (S43) of [1]. In the following we will show that this equation is not correct. Let V ∈ H(S) be a harmonic vector field, i. e. a linear combination of the basis harmonic vector fields Θ and Φ. Eq. (S43) of [1] reads in the present notation With div S (V ⊗ V) = ∇ V V = 1 2 grad S V, V and div S (grad S V+grad S V T ) = −∆ dR V+2KV = 2KV, this reduces to Taking the curl of eq. (2) results in which cannot be true in general. On the considered torus rot S K points along the toroidal direction, whereas V is harmonic and points in toroidal and poloidal direction, generally. As a consequence, eq. (S43) of [1] does not hold.

Conclusion
The letter [1] was accepted for publication knowing the arguments presented above. Furthermore Phys. Rev. Lett. refused to publish these arguments as a Comment on [1], as they do not show that the results in [1] are physically wrong. We indeed cannot judge on the physical impact, but believe that computational science should only be done with mathematically appropriate approaches and hope that this note will help for future research. While these consequences of the Hodge decomposition should be known in the mathematics community, also other examples can be found, where the vorticity-stream function approach is carelessly used, see e.g. [4][5][6].
To solve the active nematic liquid crystal model considered in [1] on general curved topographies requires an approach which also handles the harmonic parts of the velocity field and directly acts on the velocity and pressure variables. Numerical approaches can be found in [2] using discrete exterior calculus (DEC) and in [3,[7][8][9] using surface finite elements for each component of an extended velocity field in the embedding space and a penalization or Lagrange parameter for the normal component. A third approach considers a local Monge parametrizations [10]. However, non of these approaches has yet been applied to the complex problem of an active liquid crystal, considered in [1]. www.gamm-proceedings.com