## Abstract

We construct multisymplectic formulations of fluid dynamics using the inverse of the Lagrangian path map. This inverse map, the ‘back-to-labels’ map, gives the initial Lagrangian label of the fluid particle that currently occupies each Eulerian position. Explicitly enforcing the condition that the fluid particles carry their labels with the flow in Hamilton's principle leads to our multisymplectic formulation. We use the multisymplectic one-form to obtain conservation laws for energy, momentum and an infinite set of conservation laws arising from the particle relabelling symmetry and leading to Kelvin's circulation theorem. We discuss how multisymplectic numerical integrators naturally arise in this approach.

## 1. Introduction

A system of partial differential equations (PDEs) is said to be *multisymplectic* if it is of the formwhere each of the two-formsis closed. Here *z* is an ordered set of dependent variables, total differentiation with respect to each independent variable *q*^{α} is denoted by the subscript *α* after a comma, and the Einstein summation convention is used.

The closed two-form *κ*^{α} is associated with the independent variable *q*^{α}; it is analogous to the symplectic two-form for a Hamiltonian ordinary differential equation. Hence, there is a symplectic structure associated with each independent variable. In the first of a series of papers, Bridges (1997) pioneered the development of multisymplectic systems, showing that the rich geometric structure that is endowed by the symplectic two-forms can be used to understand the interaction and stability of nonlinear waves. For many important PDEs, the multisymplectic formulation has revealed hidden features that are important in stability analysis. In order to preserve at least some of these features in numerical simulations, Bridges & Reich (2001) introduced multisymplectic integrators, which generalize the symplectic methods that have been widely used in numerical Hamiltonian dynamics. Hydon (2005) showed that multisymplectic systems of PDEs may be derived from Hamilton's principle whenever the Lagrangian is affine in the first-order derivatives and contains no higher-order derivatives. This can usually be achieved by introducing auxiliary variables to eliminate the derivatives.

The aim of this paper is to provide a unified approach to produce multisymplectic formulations of fluid dynamics based on the inverse map. Our approach covers all fluid dynamical equations that are written in Euler–Poincaré form (Holm *et al.* 1998), i.e. all equations which arise due to the advection of fluid material. First, we use the inverse map to form a canonical Euler–Lagrange equation (following the Clebsch representation given in Holm & Kupershmidt (1983), which supports all possible initial conditions for vorticity, in particular, flows in three dimensions with non-zero helicity). This approach is different to that of Marsden *et al.* (2001) where the forward map is used. Then the Lagrangian is made affine in the space and time derivatives by using constraints that introduce additional variables. Following Hydon (2005), we obtain a one-form quasi-conservation law which, when it is pulled back to the space and time coordinates, gives conservation laws for momentum and energy. We also obtain a two-form conservation law that represents conservation of symplecticity; when this is pulled back to the spatial coordinates, it leads to a conservation law for vorticity. The multisymplectic version of Noether's Theorem yields an infinite space of conservation laws from the particle relabelling symmetry for fluid dynamics; these conservation laws imply Kelvin's circulation theorem. The conserved momentum that is canonically conjugate to the back-to-labels map plays a key role in the derivation of the conservation laws. The corresponding velocity is the convective velocity, whose geometric properties are discussed in Holm *et al.* (1986). The connection between particle relabelling symmetry and circulation preservation was first addressed in a multisymplectic context in Bridges *et al.* (2005).

In this paper, we show how the above constructions are made in general, illustrating this with examples. We also discuss how multisymplectic integrators can be constructed using these methods. Sections 2 and 3 review the relations among multisymplectic structures, the Clebsch representation and the momentum map associated with particle relabelling. Section 4 shows how to construct a multisymplectic formulation of the Euler–Poincaré equation for the diffeomorphism group (EPDiff), and derives the corresponding conservation laws, including the infinite set of conservation laws that yield Kelvin's circulation theorem. The section then extends this formulation to the Euler–Poincaré equation with advected quantities. This is illustrated by the incompressible Euler equation, showing how the circulation theorem arises in the multisymplectic formulation. Section 5 sketches numerical issues in the multisymplectic framework. Finally, §6 summarizes and outlines the directions for future research.

## 2. Review of multisymplectic structures

This section reviews the formulation of multisymplectic systems and their conservation laws, following Hydon (2005).

A system of PDEs is multisymplectic provided that it can be represented as a variational problem with a Lagrangian that is affine in the first derivatives of the dependent variables:(2.1)

For convenience, we restrict attention to autonomous systems where *H* depends on the dependent variables only; this can always be accomplished by introducing new dependent variables as in the Hamiltonian case.

The Euler–Lagrange equations are then(2.2)where the functions(2.3)are coefficients of the multisymplectic structure matrix. We define the (closed) symplectic two-forms(2.4)and obtain the structural conservation law (Bridges 1997) given by(2.5)Hydon (2005) showed that the Poincaré Lemma leads to a one-form quasi-conservation law(2.6)whose exterior derivative is (2.5).

Every one-parameter Lie group of point symmetries of the multisymplectic system (2.2) is generated by a differential operator of the form(2.7)Noether's Theorem implies that if *X* generates variational symmetries, that is, if for some functions *B*^{α}, then the interior product of *X* with the one-form quasi-conservation law yields the conservation law(2.8)This is the multisymplectic form of Noether's theorem.

Every multisymplectic system (2.1) is invariant under translations in the independent variables ** q**. For each of these symmetries, Noether's theorem yields a conservation lawSuch conservation laws can equally well be obtained by pulling back the quasi-conservation law (2.6) to the base space of independent variables. Commonly, the independent variables are spatial position

**and time**

*x**t*. Pulling back (2.6) to these base coordinates yields the energy conservation law from the

*t*component, and the momentum conservation law from the remaining components. We shall see the form of these conservation laws for fluid dynamics in §§4 and 5.

## 3. The inverse map and Clebsch representation

### (a) Lagrangian fluid dynamics and the inverse map

Lagrangian fluid dynamics provides evolution equations for particles moving with a fluid flow. This is typically done by writing down a flow map from some reference configuration to the fluid domain *Ω* at each instance in time. As the fluid particles cannot cavitate, superimpose or jump, this map must be a diffeomorphism.

For an *n*-dimensional fluid flow, the flow map given by specifies the spatial position at time *t* of the fluid particle that has *label* . The *inverse map* gives the label of the particle that occupies position ** x** at time

*t*as the function . The Eulerian velocity field gives the velocity of the fluid particle that occupies position

**at time**

*x**t*as follows:Each label component satisfies the advection law(3.1)Here

_{,t}and

_{,i}denote differentiation with respect to

*t*and

*x*

_{i}, respectively. We use Cartesian coordinates and the Euclidean inner product,1 so we shall not generally distinguish between ‘up’ and ‘down’ indices; summation from 1 to

*n*is implied whenever an index is repeated.

### (b) Clebsch representation using the inverse map

A canonical variational principle for fluid dynamics may be formulated by following the standard Clebsch procedure using the inverse map (Seliger & Whitham 1968; Holm & Kupershmidt 1983). The Clebsch procedure begins with a functional of the Eulerian fluid velocity , which is known as the *reduced Lagrangian* in the context of Euler–Poincaré reduction (Holm *et al.* 1998). One then enforces stationarity of the action under the constraint that equation (3.1) is satisfied by using a vector of *n* Lagrange multipliers, which is denoted as *π*. These Lagrange multipliers are the conjugate momenta to in the course of the Legendre transformation to the Hamiltonian formulation. One may choose to be solely the kinetic energy, which depends only on . More generally, will also depend on thermodynamic Eulerian variables such as density, whose evolution may also be accommodated by introducing constraints. These constraints are often called the ‘Lin constraints’ (Serrin 1959). This idea was also used in reformulating London's variational principle for superfluids (Lin 1963).

*(Clebsch variational principle using inverse map).* The Clebsch variational principle using the inverse map iswhere are Lagrange multipliers which enforce the constraint that particle labels are advected by the flow. In bounded domains, the boundary conditions are together with any boundary conditions associated with the definition of the functional .

Taking the indicated variations leads to the following equations:(3.2)(3.3)(3.4)where(3.5)and the variational derivative is defined by

*(Clebsch representation)* In the language of fluid mechanics, the expression (3.2) for the spatial momentum in terms of canonically conjugate variables is an example of a ‘Clebsch representation,’ which expresses the solution of the EPDiff equations (see below) in terms of canonical variables that evolve by standard canonical Hamilton equations. This has been known in the case of fluid mechanics for more than 100 years. For modern discussions of the Clebsch representation for ideal fluids, see Holm & Kupershmidt (1983) and Marsden & Weinstein (1983). In the language of geometric mechanics, the Clebsch representation is a momentum map.

*(Vorticity structure)* In general, one may obtain equivalent equations by enforcing the advection of any physical quantity, for example, such as circulation density (Nazarenko 1997). However, the number of constraint equations must equal the number of spatial dimensions in order to support all possible momentum distributions (and hence vorticity distributions). For example, if one only enforced the advection of one component of the fluid labels then the momentum would take the form which does not cover all possible initial conditions. For more details, see Holm & Kupershmidt (1983).

*(Local and global existence)* The standard Clebsch representation of fluid mechanics (which only has one scalar constraint) is not globally defined in certain circumstances, for instance, if the helicity is non-zero. Indeed, it is not locally defined at points where the vorticity vanishes (see Graham & Henyey 2000). However, by using sufficiently many constraints, these problems can be overcome. This is the approach that we take.

### (c) Particle relabelling

As the physics of fluids should be independent of the labelling of particles, one may relabel the particles (by a diffeomorphism of the flow domain) without changing the dynamics. This is called the *particle relabelling symmetry*. Noether's theorem applied to this symmetry leads to the Kelvin circulation theorem; see Holm *et al.* (1998) for a modern description.

### (d) Clebsch momentum map

A *momentum map* is a map from the cotangent bundle of the configuration manifold *Q* to the dual ^{*} of the Lie algebra of a Lie group *G* that acts on *Q*. The momentum map is defined by the formula,(3.6)where and . In this formula, is the infinitesimal generator of the action of *G* on *Q* associated with the Lie algebra element *ξ*, and is the natural pairing of an element of with an element of .

For the case of equations (3.2–3.4), elements of *Q* are given the labels ** l** while elements of are given by conjugate pairs of labels with their conjugate momenta.

*The Clebsch relation* *(3.2)* *defines a momentum map for the right action* *of the diffeomorphisms of the domain Ω on the back-to-labels map l*. 2

The spatial momentum in equation (3.2) may be rewritten as a map from the cotangent bundle to the dual of the vector fields (*Ω*) given by . In other words, maps the space of labels and their conjugate momenta to the space of one-form densities on *Ω*. The map may be associated with the *right action* of smooth invertible maps (diffeomorphisms) *η* of the back-to-labels maps ** l** by composition of functions, . The infinitesimal generator of this right action is obtained from its definition as(3.7)in which the vector field is tangent to the curve of diffeomorphisms

*η*

_{s}at the identity

*s*=0. Thus, pairing the map with the vector field yieldswhere is the

*L*

^{2}pairing of an element of (a one-form density) with an element of (a vector field). This is the defining relation for (3.6) to be a momentum map. ▪

Being the cotangent lift of the action of Diff(*Ω*), the momentum map is equivariant and Poisson. In other words, substituting the canonical Poisson bracket into the momentum map definition yields the Lie–Poisson bracket on the space of ** m**; see Holm & Kupershmidt (1983) and Marsden & Weinstein (1983) for more explanation, discussion and applications. The momentum map property of the Clebsch representation guarantees that the canonically conjugate variables may be eliminated in favour of the spatial momentum

**. Before its momentum map property was understood, the use of the Clebsch representation to eliminate the canonical variables in favour of Eulerian fluid variables was a tantalizing mystery (Seliger & Whitham 1968).**

*m*Note that the right action of Diff(*Ω*) on the inverse map is not a symmetry. In fact, as we shall see, the right action of Diff(*Ω*) on the inverse map generates the fluid motion itself.

### (e) Advected quantities

To construct more general fluid equations, we shall include advected quantities *a* whose flow rules are defined by(3.8)where is the Lie derivative. Such advected variables typically arise in the potential energy or the thermodynamic internal energy of an ideal fluid. For example, advected scalars *s* (e.g. salinity) satisfy(3.9)and advected densities *ρ*d*V* satisfy(3.10)A more extensive list of different types of advected quantity is given in Holm *et al.* (1998).

We write the reduced Lagrangian *ℓ* as a functional of the Eulerian fluid variables ** u** and

*a*, and add further constraints to the action

*S*to account for their advection relations,(3.11)The Euler–Lagrange equations, which follow from the stationarity condition , are(3.12)(3.13)(3.14)where the diamond operator is defined as the dual of the Lie derivative operation with respect to the

*L*

^{2}pairing. Explicitly, under integration by parts(3.15)

The map to the spatial momentum in equation (3.12)(3.16)is again a momentum map, this time for the semidirect product action of the diffeomorphisms on . Again the momentum map property allows the canonical variables to be eliminated in favour of the Eulerian quantities. As a result, eliminating the variables ** l**,

**and**

*π**ϕ*leads to the Euler–Poincaré equation with advected quantities

*a*.

### (f) Elimination theorem

Eliminating the canonically conjugate variables, and *ϕ* produces an equation of motion for , which is constructed in the proof of the following theorem:

(Elimination theorem with advected quantities) *The labels* *, their conjugate momenta* *and the conjugate momentum (ϕ) to the advected quantities (a) may be eliminated from equations* *(3.12–3.14)* *to obtain the weak form of the Euler–Poincaré equation with advected quantities*

Take the time-derivative of the inner product of with a test function ** w**, ▪

These Euler–Poincaré equations with advected quantities cover all conservative fluid equations which describe the advection of material. For a large collection of examples, see Holm *et al.* (1998).

### (g) Example: *EPDiff(*H^{1})

To give a concrete example, consider EPDiff with being the norm for ** u**. This is the

*n*-dimensional Camassa–Holm (CH) equation (Camassa & Holm 1993; Holm

*et al.*1998; Holm & Marsden 2004), which has applications in computational anatomy (Miller

*et al.*2002; Holm

*et al.*2004). This system has the reduced Lagrangianwhich defines the norm for

*λ*>0. The EPDiff equation amounts toWhen

*n*=1, these reduce to the CH equation,

### (h) Example: incompressible Euler equations

As an example, consider the reduced Lagrangian for the incompressible Euler equationsHere is the ratio of the local fluid density to the average density over *Ω*; this is governed by the continuity equation (3.10). The pressure *p* is a Lagrange multiplier that fixes the incompressibility constraint *ρ*=1. The variational derivatives in this case areConsequently, the Euler–Poincaré equations becomeand rearrangement gives the Euler fluid equations,

## 4. Inverse map multisymplectic formulation for the Euler–Poincaré equations of fluid dynamics

As we now have a canonical variational principle for fluid dynamics *via* the inverse map, one may obtain its multisymplectic formulation by extending the phase space so that the Lagrangian is affine in the space and time derivatives. In this section, we show how to do this for EPDiff(*H*^{1}) as discussed in §3, and then show how to generalize this to equations with advected quantities, using the incompressible Euler equations as an example.

### (a) Affine Lagrangian for *EPDiff(*H^{1})

After introducing the inverse map constraint, the Lagrangian becomesAny high-order derivatives and nonlinear functions of first-order derivatives must now be removed from the Lagrangian to make affine. We introduce a tensor variablethis relationship may be enforced by using Lagrange multipliers. However, it turns out that the multipliers can be eliminated and the Lagrangian becomes(4.1)which is now affine in the space and time derivatives of ** u**,

*W*,

**and**

*l***.**

*π*### (b) Multisymplectic structure

The Euler–Lagrange equations for the affine Lagrangian (4.1) areThese equations possess the following multisymplectic structure as in equation (2.2):where , and

### (c) One-form quasi-conservation law

For our multisymplectic formulation of EPDiff(*H*^{1}), the independent variables areand the dependent variables are(4.2)where *i, j* and *k* range from 1 to *n*. Comparing (4.1) with (2.1) gives the following non-zero components :Therefore the one-form quasi-conservation law amounts to(4.3)The exterior derivative of this expression yields the structural conservation law(4.4)

### (d) Conservation of energy

For EPDiff(*H*^{1}), the *t*-component of the pullback of the one-form conservation law (4.3) givesIn terms of ** u** and its derivatives, this amounts towhereThis is the energy conservation law for EPDiff(

*H*

^{1}).

### (e) Conservation of momentum

Similarly the conservation law that is associated with translations in the *x*_{i}-direction iswhich amounts to the momentum conservation law

### (f) Conservation of vorticity

Next, consider the coefficient of each in the pullback of the structural (two-form) conservation law (4.4). This iswhich amounts toOne can regard this as a vorticity conservation law for EPDiff(*H*^{1}); it is a differential consequence of the momentum conservation law.

### (g) Particle relabelling symmetry

As we discussed in §3, fluid equations in general, and EPDiff in particular, are invariant under relabelling of particles. In the context of the inverse map variables, relabelling is accomplished by the action of the diffeomorphism group defined byThe corresponding infinitesimal action of the vector fields is thenand the cotangent lift of this action isTo obtain the symmetry generator (2.7), we extend the above action to first derivatives as follows:The relabelling symmetries are variational, becauseNoether's theorem then gives the conservation lawA conservation law exists for each element ** ξ** of , so particle relabelling generates an infinite space of conservation laws.

### (h) Circulation theorem

To see how the particle relabelling conservation laws relate to conservation of circulation, note that if *ρ* is any density that satisfiesthenIf we pick a loop *C*(*t*) which is advected with the flow, thenFor a vector field ** ξ** which is tangent to the loop at time 0, and satisfies on the loop, thenfor all times

*t*, and one finds(4.5)The momentum formula (3.2) giveswhich is the circulation theorem for EPDiff.

### (i) Advected quantities

To extend this method to more general equations with advected quantities is very simple: take the Lagrangian obtained from equation (3.11) and add variables to represent higher-order derivatives. For the sake of brevity, we shall compute one example, the incompressible Euler equations, and briefly discuss the implications for the circulation theorem.

### (j) Multisymplectic form of incompressible Euler equations

We start with the reduced Lagrangianwhere *p* is the pressure and *ρ* is the relative density, and add dynamical constraints to form the LagrangianThis Lagrangian is already affine in the first-order derivatives, so the Euler–Lagrange equations are automatically multisymplectic in these variables, and take the formwhere the quantityis negative of the Hamiltonian density.

### (k) Circulation theorem for advected quantities

The conservation law for particle relabelling follows exactly as in §4 and we obtain equation (4.5) as before. The difference is that now the momentum formula (momentum map) isand so one obtainsFor the incompressible Euler equations, *a* is the relative density *ρ*, sowhich leads to the circulation theorem

*(‘Hidden’ conservation law)* Yakhot & Zakharov (1997) discuss a ‘hidden’ conservation law, which amounts to the preservation of the integral of over label space. Similarly the label-space integral of any differentiable function over label space is preserved. This corresponds to the local conservation lawwhich is a consequence of (3.13) and the incompressibility condition.

## 5. A note on multisymplectic integrators

In this section, we discuss briefly how to produce multisymplectic numerical integrators, using the inverse map formulation given in this paper. We note in particular that the multisymplectic method may have discrete forms of some of the particle relabelling symmetries and hence we will obtain a method that has discrete conservation laws.

### (a) Variational integrators

A multisymplectic integrator for a PDE is a numerical method which preserves a discrete conservation law for the two-form given in equation (2.4) (Bridges & Reich 2001). As described in (Hydon 2005), a discrete variational principle with a Lagrangian that is affine in first-order differences automatically leads to a set of difference equations which conserve the multisymplectic two-form. This now makes it very simple to construct multisymplectic integrators for fluid dynamics using the inverse map formulation: one simply replaces the spatial and time integrals in the action with numerical quadratures, replaces the first-order derivatives by differences, and takes variations following the standard variational integrator approach (Lew *et al.* 2004). While the method will preserve the discrete conservation law for the two-form *κ*, the one-form quasi-conservation law will not be preserved in general, and hence the other conservation laws will not be exactly preserved.

### (b) Discrete relabelling symmetry

As ** π** and

**are still continuous in the discretized equations, the multisymplectic integrator will have discrete particle relabelling symmetries analogous to those given in §4, provided that discretization maintains the symmetry. For example, take the following centred differencing scheme for the advection constraint of the Lagrangian**

*l**L*(choosing one space dimension for simplicity):where are the approximations to and where

*n*+1/2 and

*i*+1/2 indices indicate averaging over time levels

*n*and

*n*+1, and averaging over space-levels

*i*and

*i*+1, respectively. If we choose the relabelling generated by , then the corresponding cotangent-lifted symmetry isFollowing the variational integrator programme described in Lew

*et al.*(2004), the discrete form of Noether's theorem will give rise to discrete conservation laws for the multisymplectic method. It is an open question as to whether it is possible to obtain differencing schemes that are invariant under more general transformations of the label variables.

## 6. Summary and outlook

### (a) Summary

This paper describes a multisymplectic formulation of Euler–Poincaré equations (which are, in essence, fluid dynamical equations with a particle relabelling symmetry). We have used the inverse map to obtain a canonical variational principle following Holm & Kupershmidt (1983). As noted in Hydon (2005), a multisymplectic formulation can be obtained by choosing variables such that the Lagrangian is at most linear in the first-order derivatives and contains no higher-order derivatives. We have shown how to construct the multisymplectic formulation for the Euler–Poincaré equations for diffeomorphisms, using the example of the EPDiff(*H*^{1}) equations, and how to extend the method to the Euler–Poincaré equations with advected quantities. These equations encompass many fluid systems, including incompressible Euler, shallow-water, Euler-alpha, Green-Naghdi, perfect complex fluids, inviscid magnetohydrodynamics, etc.

The techniques of Hydon (2005) have led to conservation laws for these systems, including the usual multisymplectic conservation laws for energy and momentum plus an infinite set of conservation laws, which arise from the particle relabelling symmetries of fluid dynamics. We have highlighted the connection between these latter conservation laws and Kelvin's circulation theorem, and showed that multisymplectic integrators based on this formulation may have discrete conservation laws associated with some these symmetries.

### (b) Outlook

In §5 of this paper, we have discussed the possibility of developing multisymplectic integrators for fluids using this framework. These ideas may aid the future development of integrators that have conservation laws for vorticity and circulation, which are desirable for numerical weather prediction and other applications. It is undoubtedly simple to construct such integrators, but the issue of accuracy with time arises whenever the flow is strongly mixing and numerical errors make the label field ** l** very noisy. At present, this restricts the applicability of these methods to problems such as nonlinear internal wave propagation, which do not entangle fluid labels. One potential solution would be to try to construct an integrator where the conjugate variables can be eliminated; this would only be possible if there were a finite-dimensional Lie algebra that could be used to approximate the space of vector fields.

In a different direction, we believe that multisymplectic integrators would be especially apt for applications of EPDiff to template matching in computational anatomy (Holm *et al.* 2004). The matching problem is an initial–final value problem. In such problems, space and time may be treated on an equal footing, just as in the multisymplectic formulation.

## Acknowledgements

The work of D.D.H. was partially supported by the Royal Society of London Wolfson Award and the US Department of Energy Office of Science ASCR. The authors would like to thank the referees for their suggested improvements to the paper.

## Footnotes

↵This is only done for clarity and the equations are easily extended to the case when the domain

*Ω*is a curved manifold.↵As we discuss later, this right action contrasts with fluid particle relabelling, which arises by the left action of the diffeomorphisms on the inverse map.

- Received February 27, 2007.
- Accepted June 14, 2007.

- © 2007 The Royal Society