category: literaturenote citekey: bressanmetriplecticrelaxationequilibria2025 title: Metriplectic relaxation to equilibria authors: "Bressan, C.; Kraus, M.; Maj, O.; Morrison, P. J." year: 2025 date: 2025-00-00 2025 url: "https://arxiv.org/abs/2506.09787v1" zotero_key: HFDIK49Q zotero_storage: N5SF22AN collections: transpordi_planeerimine folder: 001_artiklid firstAuthor: "Bressan, C."

status: converted

Metriplectic relaxation to equilibria

C. Bressan¹ , M. Kraus^a , O. Maj^a , P. J. Morrison^b

^a Max Planck Institute for Plasma Physics, Boltzmannstrasse 2, 85747, Garching, Germany. ^b Department of Physics and Institute for Fusion Studies, University of Texas at Austin, Austin TX 78712-1060, USA.

Abstract

Metriplectic dynamical systems consist of a special combination of a Hamiltonian and a (generalized) entropy-gradient flow, such that the Hamiltonian is conserved and entropy is dissipated/produced (depending on a sign convention). It is natural to expect that, in the long-time limit, the orbit of a metriplectic system should converge to an extremum of entropy restricted to a constant-Hamiltonian surface. In this paper, we discuss sufficient conditions for this to occur, and construct a class of metriplectic systems that satisfy one such condition. The proposed metriplectic structure is inspired by the Landau operator for Coulomb collisions in plasmas, which is included as special case. We apply these results to the construction of relaxation methods for the solution of equilibrium problems in fluid dynamics and plasma physics.

Keywords: Dissipative Dynamical Systems, Metriplectic Systems, Hamiltonian Systems, Lyapunov theorem, Landau Collision Operator, Magnetohydrodynamics, Fluid Dynamics

| 1 | | Introduction | 2 | | |---|---------------------------------------------------------------------------|-----------------------------------------------------------------------------|----|--| | | 1.1 | Overview of equilibrium calculations
| 3 | | | | 1.2 | Overview of metriplectic relaxation | 4 | | | 2 | Metriplectic dynamics and variational principles for equilibria of fluids | | | | | | | and plasmas | 7 | | | | 2.1 | Metriplectic dynamics
| 7 | | | | 2.2 | Examples of equilibrium problems
| 10 | | | | | 2.2.1
Reduced Euler equations
| 10 | | | | | 2.2.2
Axisymmetric MHD equilibria
| 13 | | | | | 2.2.3
Beltrami fields
| 16 | | | 3 | | On the relaxation of metriplectic systems | 18 | | | | 3.1 | Finite-dimensional systems: Lyapunov stability
| 18 | | | | 3.2 | Finite-dimensional systems: the Polyak– Lojasiewicz inequality
| 25 | | | | 3.3 | Infinite-dimensional systems: tentative generalization of the PL inequality | 35 | |

Preprint submitted to Communications in Nonlinear Science and Numerical Simulation June 12, 2025

| 4 | | Two examples: metric double brackets and projectors | 37 | | |---|--------------------------------|-------------------------------------------------------------|----|--| | | 4.1 | Metric double brackets | 40 | | | | 4.2 | Projector-based metric bracket
| 46 | | | 5 | Collision-like metric brackets | | | | | | 5.1 | General construction of collision-like brackets
| 50 | | | | 5.2 | Collision-like brackets based on div–grad operators
| 58 | | | | 5.3 | Collision-like brackets based on curl–curl operators | 60 | | | | 5.4 | Application to the reduced Euler equations
| 60 | | | | 5.5 | Application to Grad-Shafranov equilibria
| 64 | | | 6 | | Diffusion-like metric brackets | 67 | | | | 6.1 | General construction of the brackets
| 67 | | | | 6.2 | Diffusion-like brackets based on div–grad operators | 69 | | | | 6.3 | Diffusion-like brackets based on curl–curl operators | 70 | | | | 6.4 | Application to nonlinear Beltrami fields
| 71 | | | 7 | | Conclusions | 76 | | | | Appendix | A
Bilinear forms and Leibniz identity | 77 | | | | Appendix | B
On continuous bilinear forms on Hilbert spaces | 79 | | | | Appendix | C
A Lagrangian variational principle for Beltrami fields | 80 | |

1. Introduction

There are two main purposes of this paper: investigate sufficient conditions for metriplectic relaxation (reviewed in Section 2.1) to occur and use metriplectic relaxation to find equilibria of fluid dynamics and plasma physics systems. Metriplectic dynamical systems, as introduced in [1, 2, 3], are designed to formally converge to an extremum of entropy while being restricted to a constant-Hamiltonian surface. Here we more rigorously examine the conditions for such relaxation and then construct and investigate metriplectic systems that achieve the relaxation for finding equilibria of a collection of fluid and plasma physical systems. In the remainder of this section, in Section 1.1 we first give an overview of the challenges and previous methods for calculating equilibria, followed in Section 1.2 by our and others' previous metriplectic relaxation methods.

Email addresses: michael.kraus@ipp.mpg.de (M. Kraus), omar.maj@ipp.mpg.de (O. Maj), morrison@physics.utexas.edu (P. J. Morrison)

¹Formerly at the Max Planck Institute for Plasma Physics, Boltzmannstrasse 2, 85747, Garching, Germany.

1.1. Overview of equilibrium calculations

The calculation of equilibria of physical systems often leads to ill-posed nonlinear problems, where the ill-posedness is due to the nonuniqueness of the solution. Additional constraints are needed to define uniquely the equilibrium of interest, depending on the application at hand. In some situations, prescribing enough constraints to determine a unique equilibrium may not be straightforward. This lack of uniqueness for equilibrium problems is precisely discussed below in Section 2.2 for examples taken from fluid dynamics and magnetohydrodynamics (MHD): equilibria of the Euler equations in vorticity form reduced to two dimensions [4, p.488], axisymmetric MHD equilibria [5], linear and nonlinear Beltrami fields.

In some cases, after providing additional physical constraints, the equilibrium conditions can be reformulated as a well-posed mathematical problem. This is the case, for instance, for the Euler equations, axisymmetric MHD equilibria, and linear Beltrami fields.

In more complicated situations, such as for nonlinear Beltrami fields and full MHD equations in three dimensions, the problem of computing an equilibrium point has no good solution. A mathematical perspective on the intricacies of this problem can be found in the introduction of the paper by Bruno and Laurence [6]. A large fraction of MHD equilibrium calculations in three dimensions are based on a reformulation of the problem in which one assumes that the magnetic field is tangent to a family of nested toroidal surfaces [7, 8]. This strategy has been extremely successful for the design of stellarators [9], since it reduces significantly the complexity of the problem. While this approach gives an acceptable representation of the magnetic field in the core of a stellarator with modest computational cost, it cannot account for more complicated magnetic field configurations, such as those with magnetic islands and chaotic field lines, due to the built-in foliation of the domain by toroidal surfaces. Yet magnetic-field islands and chaotic regions are relevant in practice and calculations of MHD equilibria beyond the paradigm of nested toroidal surfaces are needed. An iterative procedure for the calculation of general MHD equilibria has been outlined by Grad [10], and a similar iterative procedure is implemented in the PIES code [11]. Such iteration schemes are purely heuristics: there is no theoretical control on the convergence.

Another approach to the computation of equilibria is based on artificial relaxation. Relaxation methods solve the Cauchy problem for a fictitious dissipative evolution law that contains a tailored dissipation mechanism. If the dissipation mechanism is well designed, the solution of the Cauchy problem, with a given (well-prepared) initial condition, exists globally in time, has a limit for t → +∞, and the limit is an equilibrium of the considered physical system. The dynamical evolution itself might not be physical, but the solution should converge to a physical equilibrium as fast as possible. Some care may be taken to preserve important properties of the solution. For the specific case of MHD, for instance, one can evolve the magnetic field B according to Faraday's equation,

$$\partial_t B + c \operatorname{curl} E = 0, \quad B(0, x) = B_0(x),$$

but with a properly chosen effective electric field E, which does not need to have physical meaning. (Note, Gaussian units are used throughout this paper, with c being the speed of light in free space.) If the initial condition B⁰ satisfies div B⁰ = 0, then div B = 0 for all t ≥ 0. Probably the most intuitive relaxation method can be obtained by choosing E ∝ −U × B, where the advecting velocity field U solves the viscous MHD momentum balance equation [12]. The idea of this method is physically intuitive: magnetic energy is converted into kinetic energy and dissipated by viscosity. Since there is no resistivity, magnetic helicity is preserved and this provides a lower bound for the dissipation of magnetic energy [13, 14]. The evolution of the system is not physically consistent (the viscosity term is usually very simple and resistivity is zero), but the relaxed state is guaranteed to be an ideal MHD equilibrium.

From the application point of view, relaxation methods constructed in this way are not fully satisfactory because of the following reasons: (1) Faraday's equation reduces to Lie-dragging of the magnetic field and therefore preserves the magnetic flux and the field-line topology of the initial condition (frozen-in law [15, 5]); hence not all equilibrium points are accessible from a given initial condition [12] and this might be an advantage or a disadvantage depending on the intended application. (2) The method does not offer any mechanism to control important properties of the equilibrium such as the pressure profile, i.e., the constant value of the pressure on the surfaces tangent to the magnetic field (flux surfaces) or the shear of the magnetic field. (3) Another drawback is that the relaxation mechanisms based solely on viscosity do not lead to the shortest path from the initial condition down to an equilibrium point.

An example of a relaxation method based on viscosity is implemented in the HINT code [16, 17]. In HINT, pressure is relaxed with an ad hoc algorithm, in a separate step, during the magnetic field relaxation. If resistivity is accounted for in the relaxation of the magnetic field, then the topology of the magnetic field lines can change, but with finite resistivity helicity is no longer preserved and there is no lower bound for the magnetic energy.

Another relaxation method that seeks a faster way to relax the magnetic energy is based on the variational principle for the equilibrium conditions (reviewed in Appendix C for the case of Beltrami fields). This method has been proposed by Chodura and Schl¨uter [18], cf. also Moffatt [14, sec. 8.2], and it can be specialized to the case of nonlinear Beltrami fields [19]. The idea is to Lie drag both the magnetic field and the pressure with an advecting velocity field U chosen to guarantee the maximum decay rate of the magnetic energy. These ideas are strictly related to the modern theory of optimal transport of differential forms [20].

The idea of Lie dragging both the magnetic field and the pressure has been exploited in the SIESTA code as well [21], where each displacement of both magnetic field and pressure is generated by an infinitesimal "Lie dragging step".

1.2. Overview of metriplectic relaxation

Metriplectic dynamics and its concomitant variational principles for equilibria will be thoroughly reviewed in Section 2. In this subsection we give a overview of the paper while describing the usage in this work, where we explore using artificially constructed metriplectic dynamical systems in order to construct relaxation methods for the calculation of equilibrium points. Metriplectic dynamics was introduced by Morrison [1, 2, 3] as a generalization of noncanonical Hamiltonian dynamics with the aim of including dissipative phenomena. (See [22, 23, 24, 25, 26] for recent developments.) The equation of evolution is constructed in terms of two algebraic structures: a Poisson bracket [4], which is antisymmetric and defines the Hamiltonian part of the equation, and a metric bracket, which is symmetric and accounts for dissipation. In addition to the brackets, a Hamiltonian function H and an entropy function S are given, satisfying appropriate compatibility conditions. As a direct consequence of the construction, the Hamiltonian H is conserved and the entropy S is dissipated (more precisely, it is nonincreasing). The fact that entropy is a monotonic function of time quantifies the dissipation of the system; whether entropy grows or decays is just a matter of a sign convention and we choose to dissipate entropy, which in this work is regarded as a Lyapunov function [27] more than a physical entropy. Many physically relevant mathematical models have been found to possess a metriplectic structure. For instance the Vlasov-Maxwell-Landau system [3], various fluid mechanical systems [2, 28, 24, 26], visco-resistive MHD [29], and dissipative extended MHD [22]. This structure can also be exploited in order to design numerical schemes that preserve the key features of the physical model [30, 31]. Here instead we shall use metric brackets as equilibrium solvers.

Not to be confused with metriplectic relaxation is a relaxation method that uses the Hamiltonian structure. This method, which is based on squaring the Poisson bracket, was introduced in [32, 33] for two-dimensional vortical motion of neutral fluids; it was later generalized so as to make it more effective and work in a broader context in Flierl and Morrison [34] and it was applied to reduced MHD problems by Furukawa and coworkers [35, 36, 37, 38] (see [39] for recent review). This approach, which has been named "double bracket" dynamics or simulated annealing, has different properties as compared to metriplectic dynamics: with double brackets, the Hamiltonian is dissipated while the system evolves on the a specific hypersurface (a symplectic leaf) determined by the constants of motion built into the Poisson brackets (Casimir invariants). Another early approach is that of Brockett who used a version of the double bracket for matrices constructed out of the commutator [40] (see also the work of Bloch and coauthors [41, 42]).

Since metriplectic dynamics dissipates entropy but preserves the Hamiltonian, one might expect that a global solution, if it exists, will approach a minimum of the entropy function restricted to the level set of the Hamiltonian that contains the initial condition. Specifically, if we denote by u a point in the phase space V of the system, we might expect that the solution of a generic metriplectic system with Hamiltonian H(u), entropy S(u), and initial condition u(0) = u⁰ ∈ V would converge to a solution of the variational principle

$$\min{\mathcal{S}(u) \colon u \in V, \ \mathcal{H}(u) = \mathcal{H}(u_0)}. \tag{1}$$

Variational principles of the form (1) are physically relevant since, the resulting state of a physical relaxation mechanism can, in certain cases, be well approximated by the solution of such a minimum entropy principle. Typically one envisages a situation where the physical evolution of the system is nearly ideal, that is, dissipation mechanisms are very small, but such small nonideal effects are sufficient to dissipate one of the ideal constant of motion, while causing negligible variations in the other constants of motion. This is referred to as selective decay and plays an important role in the relaxation of fluids and plasmas to a self-organized state [43, 44]. In MHD for instance, linear Beltrami fields are found as a result of processes that dissipate magnetic energy, while (approximately) preserving magnetic helicity as argued by Woltjer [45]. The precise physical relaxation mechanism has been discussed in detail by Taylor and it is referred to as Woltjer-Taylor relaxation [46, 47, 48]. In Section 2.2 we shall see that equilibrium problems can often be reformulated as a variational principle of the form (1): among the many possible solutions of the equilibrium conditions, the variational principle (1) selects only those that minimize entropy on a constant energy hypersurface, and thus reduces significantly the issue of nonuniqueness of the equilibrium.

When the solution of a metriplectic system has a limit for t → +∞ and the limit is a solution of (1), we say that the system has completely relaxed. Unfortunately, complete relaxation does not always happen. It depends on the null space of the metric brackets, which is defined precisely in Section 2.1. We shall demonstrate complete relaxation (or lack thereof) by means of numerical experiments. We shall propose and test a particular class of brackets modeled upon Morrison's brackets for the Landau collision operator [1, 3] and show by means of numerical experiments that these new brackets completely relax an initial condition. These new brackets are referred to as collision-like metric bracket.

Collision-like brackets have the disadvantage of generating integro-differential evolution equations, that are usually computationally expensive (although efficient methods exist [49]). In an attempt to reduce the computational cost of the relaxation methods, we have introduced a simplified version of the collision brackets that are local and thus lead to pure partial differential equations that have the structure of diffusion equations. We refer to these simplified brackets as diffusion-like. We are able to recover known relaxation methods, such as the metriplectic bracket based on Nambu dynamics of [42], which was also obtained and proposed in the context of vortex dynamics in [50, 51], and the method of Chodura and Schl¨uter [18], as special cases of diffusion-like brackets. All of the brackets used in this paper follow naturally from the inclusive 4-bracket construction given in [23].

The remainder of the paper is structured as follows. As noted above, in Section 2 we cover review material: we recall the precise definition of metriplectic systems in Section 2.1 and describe the equilibrium problems that we consider as test cases together with corresponding variational principles in Section 2.2. Section 3 presents some mathematical results on relaxation and the relaxation rate for metriplectic systems, including results that we are unable to find in the literature. In Sections 3.1 and 3.2 we prove extensions of the Lyapunov stability theorem and the Polyak– Lojasiewicz condition for the rate of relaxation, for finite-dimensional systems with the inclusion of constraints, respectively, while in Section 3.3 we make some comments on the extension of these results to infinitedimensional systems. In Section 4, we address simple examples of metric brackets and study the issue of complete relaxation, both analytically and numerically. Section 4.1 describes brackets, which we refer to as metric double brackets, that fail to completely relax, while Section 4.2 describes projection based brackets that do completely relax. Collision-like brackets are introduced in Section 5 with theory presented in Sections 5.1, 5.2, and 5.3 and numerical experiments in Sections 5.4 and 5.5 for the Euler equations in vorticity form and for the Grad-Shafranov MHD equilibria, respectively. For these applications complete relaxation of the solution is critical. Section 6 is dedicated to diffusion-like brackets, with their general construction and various forms given in Sections 6.1, 6.2, and 6.3. We shall see that their properties make them most suitable for applications such as the calculation of nonlinear Beltrami fields, which are considered in Section 6.4, for which the property of complete relaxation is not needed. Finally, we conclude in Section 7.

2. Metriplectic dynamics and variational principles for equilibria of fluids and plasmas

In this section, we review the necessary background material. First we briefly recall the definition and basic properties of metriplectic dynamics. Then we formulate and discuss examples of equilibrium problems, and their variational formulation. These problems will be used as a test bed for the metriplectic relaxation method proposed in this work.

2.1. Metriplectic dynamics

Metriplectic dynamics, being a special kind of continuous-time dynamical system, is determined by a phase space and a law governing the evolution in time of a point in the phase space. In this work, we are mainly concerned with infinite-dimensional metriplectic systems with a phase space given by a Banach space V of functions over a domain Ω ⊂ R ^d with values in R ^N , where d, N ∈ N. (We use "domain" as a shorthand for "open, connected set".) We shall however briefly discuss finite-dimensional examples for sake of clarity and simplicity. In the finite-dimensional case the phase space is chosen to be an open connected subset Z ⊆ R ⁿ, n ∈ N, with coordinates z = (z i ) n ⁱ=1.

For any F ∈ C 1 (V ), we denote by DF(u) its Fr´echet derivative [52] at the point u ∈ V . We shall also need the functional derivative of F. If W is another Banach space with a nondegenerate pairing ⟨·, ·⟩^V ^×^W : V × W → R, we can define the functional derivative δF/δu of F in W as the unique element of W, if it exists, such that [53]

$$D\mathcal{F}(u)v = \left\langle v, \frac{\delta \mathcal{F}(u)}{\delta u} \right\rangle_{V \times W}, \quad \forall v \in V.$$ (2)

Unless otherwise stated, in this work we assume V ⊆ L 2 (Ω, µ; R ^N ), W = L 2 (Ω, µ; R ^N ), and the pairing is the L ² product with respect to a given measure dµ = m(x)dx, where m is a smooth and integrable function and dx the Lebesgue measure (volume element) on Ω,

$$(u,v){L^2} = \int{\Omega} u(x) \cdot v(x) d\mu(x), \quad u,v \in L^2(\Omega,\mu;\mathbb{R}^N).$$

(The nontrivial measure is needed in order to accommodate cases such as Grad-Shafranov equilibria discussed below.)

A metriplectic dynamical system is specified by giving two functions H, S ∈ C^∞(V ), namely the Hamiltonian and the entropy, respectively, together with compatible Poisson and metric brackets on V . We recall the definitions [54, 1, 4, 53, 3].

A Poisson bracket on V is a bilinear antisymmetric map

$${\cdot,\cdot}: C^{\infty}(V) \times C^{\infty}(V) \to C^{\infty}(V),$$

such that, for any F, G, and H in C^∞(V ),

$${\mathcal{F}, \mathcal{GH}} = {\mathcal{F}, \mathcal{G}}\mathcal{H} + \mathcal{G}{\mathcal{F}, \mathcal{H}},$$ (3a)

$${\mathcal{F}, {\mathcal{G}, \mathcal{H}}} + {\mathcal{G}, {\mathcal{H}, \mathcal{F}}} + {\mathcal{H}, {\mathcal{F}, \mathcal{G}}} = 0. \tag{3b}$$

Equations (3a) and (3b) are referred to as the Leibniz identity and the Jacobi identity, respectively. Then a Poisson bracket defines a Lie algebra structure on C^∞(V ) which in addition is a derivation in each argument.

A metric bracket on V is a bilinear symmetric map

$$(\cdot,\cdot)\colon C^{\infty}(V)\times C^{\infty}(V)\to C^{\infty}(V),$$

such that, for any F in C∞(V ),

$$(\mathcal{F}, \mathcal{F}) \ge 0.$$

By definition, the Poisson bracket must satisfy Leibniz and Jacobi identities. Leibniz identity, in particular, implies at least formally that the bracket can be written in term of the functional derivatives of its arguments, cf. Appendix A for a precise definition. Usually, the symmetric bracket does not need to satisfy any condition other then bilinearity, symmetry, and positive semidefiniteness. However, if one requires the symmetric bracket to satisfy Leibniz identity,

$$(\mathcal{F}, \mathcal{GH}) = (\mathcal{F}, \mathcal{G})\mathcal{H} + \mathcal{G}(\mathcal{F}, \mathcal{H}),$$ (4)

then both the Poisson and the symmetric brackets have a similar representation, i.e.,

$${\mathcal{F}, \mathcal{G}} = \sum{k,l=1}^{N} \int{\Omega} \int{\Omega} \frac{\delta \mathcal{F}(u)}{\delta u{k}}(x) \mathcal{J}{kl}(u; x, x') \frac{\delta \mathcal{G}(u)}{\delta u{l}}(x') d\mu(x') d\mu(x), \tag{5a}$$

and analogously

$$(\mathcal{F}, \mathcal{G}) = \sum{i,j} \int{\Omega} \int_{\Omega} \frac{\delta \mathcal{F}(u)}{\delta ui}(x) \mathcal{X}{ij}(u; x, x') \frac{\delta \mathcal{G}(u)}{\delta u_j}(x') \, d\mu(x') \, d\mu(x), \tag{5b}$$

where the functional derivatives are computed with respect to the L ² product with a given measure µ on Ω. The kernels J (u) and K (u) define an anti-symmetric and a symmetric, positive semidefinite operator, J(u) and K(u), respectively. In finite dimensions, Eqs. (5) take the form

$$\begin{split} {\mathcal{F},\mathcal{G}} &= J^{kl}(z) \frac{\partial \mathcal{F}(z)}{\partial z^k} \frac{\partial \mathcal{G}(z)}{\partial z^l}, \ (\mathcal{F},\mathcal{G}) &= K^{ij}(z) \frac{\partial \mathcal{F}(z)}{\partial z^i} \frac{\partial \mathcal{G}(z)}{\partial z^j}, \end{split}$$

where the sum over repeated indices is implied, J(z) is an antisymmetric contravariant tensor, and K(z) is a symmetric positive semidefinite contravariant tensor over the domain Z ⊆ R ^N . In particular J is referred to as the Poisson tensor.

The evolution equation for a metriplectic system u(t) ∈ V is formulated as an evolution equation for arbitrary functions of u(t), that is,

$$\frac{d\mathcal{F}}{dt} = {\mathcal{F}, \mathcal{H}} - (\mathcal{F}, \mathcal{S}), \quad \text{for all } \mathcal{F} \in C^{\infty}(V),$$ (6a)

where H, S ∈ C^∞(V ) are the Hamiltonian and entropy functions, respectively, whereas {·, ·} and (·, ·) are the Poisson and metric bracket on V , respectively, satisfying the compatibility conditions

$${\mathcal{F}, \mathcal{S}} = 0, \quad (\mathcal{F}, \mathcal{H}) = 0, \quad \text{for all } \mathcal{F} \in C^{\infty}(V).$$ (6b)

If both brackets satisfy the Leibniz identity, the strong form of the evolution equation then reads

$$\partial_t u = J(u) \frac{\delta \mathcal{H}(u)}{\delta u} - K(u) \frac{\delta \mathcal{S}(u)}{\delta u}.$$

In general, both J(u) and K(u) have a nontrivial null spaces. Per definition, the null space of a bracket is identified with that the corresponding operators J(u) and K(u). In finite dimensions, the null space of J(z) is spanned by the functional gradient of Casimir invariants and equilibria from the variational principle align with those of the equations of motion [4]. We note, however, that there are subtleties at points where the rank of J(z) changes [55] and equilibria at such points can possesses nearby behavior that is not Hamiltonian [56, 57]. The null space of K(z) contains at least the gradient of the Hamiltonian, due to the compatibility condition (6b); similar remarks on rank changing could apply. In general, the vector field J(u)δH/δu can be viewed as a generalization of Hamiltonian flow ω(u) −1 δH/δu where the inverse of the symplectic operator ω(u) is replaced by a possibly noninvertible Poisson operator J(u). On the other hand, the vector field −K(u)δH/δu can be viewed as generalization of a gradient flow −G(u) −1 δH/δu where the inverse of the metric operator G(u) is replaced by a symmetric, positive semidefinite operator. Metriplectic dynamics combines the (generalized) symplectic and gradient flows. Typically, the symplectic part describes the ideal dynamics, while the gradient flow accounts for a nonideal relaxation mechanism. For this reason, accepting a slight abuse of terminology, we refer to symmetric brackets (·, ·) with the Leibniz property as metric brackets. We shall always tacitly assume that H and S are not functionally dependent, i.e., δS(u)/δu is not everywhere parallel to δH(u)/δu, otherwise the metric bracket part vanishes identically.

We are mostly interested in the dynamical systems generated by metric brackets, i.e., we shall drop the Poisson bracket part,

$$\frac{d\mathcal{F}}{dt} = -(\mathcal{F}, \mathcal{S}), \quad \text{for all } \mathcal{F} \in C^{\infty}(V), \tag{7a}$$

where S is the entropy function and the metric bracket satisfies the compatibility condition

$$(\mathcal{F}, \mathcal{H}) = 0$$ , for all $\mathcal{F} \in C^{\infty}(V)$ , (7b)

where H is the Hamiltonian.

A solution of either (6) or (7) satisfies

$$\frac{d\mathcal{H}(u)}{dt} = 0, \qquad \frac{d\mathcal{S}(u)}{dt} \le 0, \tag{8}$$

that is, both system (6) and (7) dissipate entropy on the surface of constant energy (Hamiltonian).

Because of (8) one may expect that solutions of the variational principle (1) are necessarily equilibria of a metriplectic system. Indeed this is the case and it follows from the method of Lagrange multipliers [53], which gives a necessary condition for (1): if u is a solution of (1), then there is a constant λ ∈ R (the Lagrange multiplier) such that

$$DS(u) - \lambda D\mathcal{H}(u) = 0, \quad \mathcal{H}(u) = \mathcal{H}_0.$$ (9)

Alternatively, one can write the Lagrange condition using the functional derivative, if they exist,

$$\frac{\delta \mathcal{S}(u)}{\delta u} - \lambda \frac{\delta \mathcal{H}(u)}{\delta u} = 0, \quad \mathcal{H}(u) = \mathcal{H}_0.$$

Equations (9) constitute a system of two equations for the pair (u, λ) ∈ V × R. Let us assume that the set of solutions

$$\mathfrak{C}_{\mathcal{H}_0} := { u \in V : \exists \lambda \in \mathbb{R} \text{ such that } (u, \lambda) \text{ solves } (9) },$$

is nonempty (CH⁰ ̸= ∅), the restriction S to CH⁰ is bounded from below, and the minimum is attained, that is, there are points u^e ∈ CH⁰ where S(ue) = min{S(u): u ∈ CH⁰ }. Then, the solutions of (1) correspond to those points u^e of CH⁰ where S attains its minimum. The set CH⁰ is the set of constrained critical points of S, that is, of critical points of S restricted to the energy surface H(u) = H0.

If the symmetric bracket satisfies the Leibniz identity, the Lagrange condition together with either the compatibility condition (6a) or (7b) imply that any point in C^H⁰ , i.e., any constrained critical point of S, is an equilibrium point of the metriplectic system, and thus, in particular, any solution u^e of (1) is necessarily an equilibrium point [42, sec. 4.1]. The converse, however, is not true, since all constrained critical points are equilibria, not just the minima. In addition, there can be equilibrium points of either (6) or (7) that are not constrained critical points of S. One example is given in Section 4.1 below.

The fact that the set of equilibrium points of a metriplectic system can, in general, be (much) larger than the set C^H⁰ of constrained critical points of S can be an obstruction to convergence of an orbit u(t) to a solution of (1), and this is important in some (but not all) applications. In the next section, we review a few physically relevant equilibrium problems and discuss their relation to variational principles of the form (1). The main result of the paper is the construction of appropriate metric brackets that relax a given initial condition to a solution of such equilibrium problems.

2.2. Examples of equilibrium problems

In this section, we review the examples of equilibrium problems that we shall use as test cases for metriplectic relaxation. All considered test problems are mathematically ill-posed, because they admit multiple solutions. In some cases, the ill-posedness can be mitigated by adding additional physics constraints. We shall also discuss the variational principles for the considered equilibrium problems.

2.2.1. Reduced Euler equations

We begin with the Euler equations reduced to two dimensions [4, p.488 and references therein], which is the simplest of a hierarchy of models including the reduced MHD model [58]. Let x = (x1, x2) be Cartesian coordinates in a bounded domain Ω ⊂ R 2 with a sufficiently regular boundary ∂Ω. In Ω, we consider an incompressible flow U = (U1, U2) = ∂2ϕ, −∂1ϕ , given in terms of a stream function ϕ(x) with ∂ⁱ = ∂/∂xⁱ . Then div U = 0 and from definition of the scalar vorticity ω := ∂1U² − ∂2U¹ one obtains the Poisson equation

$$-\Delta \phi = \omega \text{ in } \Omega, \quad \phi = 0 \text{ on } \partial \Omega,$$ (10)

where ∆ is the Laplace operator in R 2 . Vice versa, given ω, we can solve (10) for ϕ and reconstruct the flow U. Hence, the incompressible Euler equations in two-dimensions

$$\partial_t U + U \cdot \nabla U = -\nabla p$$ , div $U = 0$ ,

with p being the pressure field, amount to an evolution equation for the scalar vorticity ω,

$$\begin{cases} \partial_t \omega + [\omega,\phi] = 0, & \text{in } \Omega, \ -\Delta \phi - \omega = 0, & \text{in } \Omega, \ \phi = 0, & \text{on } \partial \Omega, \end{cases}$$

where [ω, ϕ] := ∂1ω ∂2ϕ − ∂2ω ∂1ϕ is the canonical Poisson bracket in R . This model is referred to as the reduced Euler equations.

The phase space V of the reduced Euler equations is the space of vorticity fields, i.e., u = ω, and ϕ = −∆ −1 ^Ω,⁰ ω is regarded as a function of ω, given by the inverse of the Laplacian on Ω with homogeneous Dirichlet boundary conditions.

The equilibrium problem for the reduced Euler equations then reads

$$\begin{cases} [\omega, \phi] = 0, & \text{in } \Omega, \ -\Delta \phi - \omega = 0, & \text{in } \Omega, \ \phi = 0, & \text{on } \partial \Omega. \end{cases}$$ (11)

Problem (11) admits many solutions and therefore is mathematically ill-posed. One can in fact construct a large class of solutions upon noticing that [ω, ϕ] = 0 implies that ω is constant on the isolines (contours) of the potential ϕ. The contours may have many connected components, each one with a possibly different topology (i.e., homeomorphic to a different model space) and the constant value of ω on different connected components may be different. Given a function f ∈ C 1 (R), we may set ω = λf(ϕ) with a normalization factor λ ∈ R to be determined; in this way we assign the same value of ω to all connected components of the same contour of ϕ. This is a special case which is considered here for sake of simplicity. Then we consider the problem: find (ϕ, λ), with λ ̸= 0, such that

$$\begin{cases} -\Delta \phi = \lambda f(\phi), & \text{in } \Omega, \ \phi = 0, & \text{on } \partial \Omega. \end{cases}$$ (12)

Any solution (ϕ, λ) of (12) yields a solution ω = λf(ϕ) of the equilibrium problem (11). The case λ = 0 leads to the trivial equilibrium ω = 0 and it is not considered.

Problem (12) is a "eigenvalue problem" for a semilinear elliptic equation. If f ′ (y) ≤ 0 (respectively, f ′ (y) ≥ 0) for all y, the equation has a unique solution for any λ ≥ 0 (respectively, λ ≤ 0) [59]. In the other cases, the solution may not exist for all λ; if a solution exists, uniqueness is not guaranteed, e.g., for degenerate eigenvalues. For instance, when f(y) = y, problem (12) reduces to the standard eigenvalue problem for the Laplace operator with Dirichlet boundary conditions; then, we have discrete, possibly degenerate, positive eigenvalues λⁿ > 0, n ∈ N, each with a corresponding finite set of eigenfunctions ϕn,k depending on the multiplicity of λn. For λ ≤ 0, the trivial solution ϕ = 0 is the unique solution.

This, in particular, shows that problem (11) is ill-posed, since there is a rich set of solutions for each choice of f, and many choices of f are possible. In order to mitigate the nonuniqueness problem, in practice the function f, which will be referred to as the equilibrium profile, is prescribed, and among the solutions of (12), the one with the lowest λ is considered. The reformulation of the equilibrium problem (11) into the nonlinear eigenvalue problem (12) with fixed f is good enough in practice. There are efficient iterative algorithms [60] for the solution of this type of eigenvalue problems with the lowest λ.

An alternative reformulation of the equilibrium problem (11) is possible, based on a variational principle of the form (1). For H⁰ > 0 and s ∈ C 2 (R) satisfying s ′′(y) ̸= 0, y ∈ R, let us consider the problem

$$\min{S(\omega): \mathcal{H}(\omega) = \mathcal{H}_0}.$$ (13)

where

$$S(\omega) = \int{\Omega} s(\omega)dx$$ , $\mathcal{H}(\omega) = \frac{1}{2} \int{\Omega} |\nabla \phi|^2 dx$ , (14)

are the entropy and Hamiltonian functions, respectively, with ϕ regarded as a function of ω defined via (10). The condition on the entropy profile s(y) implies that s ′ (y) is a strictly monotonic function, either decreasing or increasing. In general, s(y) will be chosen ad hoc and does not necessarily have a physical meaning. The Hamiltonian H is the kinetic energy of the incompressible fluid, since |∇ϕ| ² = |U| , where U is the flow velocity.

Problem (9) in this case reads: Find (ω, λ), such that

$$\begin{cases} s'(\omega) - \lambda \phi = 0, & \mathcal{H}(\omega) = \mathcal{H}_0, \ \text{with } \phi \text{ solution of (10).} \end{cases}$$ (15)

The set C^H⁰ of constrained critical points, cf. Section 2.1, amounts to

$$\mathfrak{C}_{\mathcal{H}_0} = {\omega \colon (\omega, \lambda) \text{ solves } (15) \text{ for some } \lambda \in \mathbb{R}}.$$

All elements of the set C^H⁰ are solutions of the original equilibrium problem (11). In fact, for λ ̸= 0, we have [ω, ϕ] = λ −1 [ω, s′ (ω)] = 0. The case λ = 0 is somewhat special, since s ′ (ω) = 0 with s ′′(y) ̸= 0 implies that ω(x) = ωc(x) = y^c = constant, where y^c is the unique zero of s ′ (y); the corresponding potential ϕ^c is given by the solution of problem (10) with constant right-hand side. There is therefore only one solution for λ = 0 and this carries the energy ¹ ∥∇ϕc∥ 2 ^L² = H^c and it is always an equilibrium, since ω^c is constant. If H⁰ = Hc, this solution belongs to C^H⁰ , otherwise λ = 0 is not a possible value for the Lagrange multiplier.

Under the hypothesis s ′′(y) ̸= 0, s ′ (y) is monotonic, thus an invertible function of y ∈ R. Problem (15) is related to the eigenvalue problem (12) with equilibrium profile given by f(y) = (s ′ ) −1 (y). Precisely, if (ω, λ) is a solution of (15) with λ ̸= 0, we can define ˜ω = λω and ϕ˜ = λϕ, and obtain from (15)

$$-\Delta \tilde{\phi} = \lambda \omega = \lambda (s')^{-1} (\lambda \phi) = \lambda f(\tilde{\phi}),$$

which shows that ϕ˜ solves (12) with the eigenvalue being the same as the Lagrange multiplier λ.

Among all these equilibria, the minimization in (13) selects those with minimum entropy, thus mitigating the nonuniqueness problem, as shown in the following special case. As in Yoshida and Mahajan [44], we choose s(y) = y ²/2, hence s ′ (y) = y and solutions of (15) must necessarily solve the eigenvalue problem for the Laplace operator on Ω with homogeneous Dirichlet boundary conditions. From standard theory [61], we know that there is an orthonormal basis {ϕj,k} in L 2 (Ω) of eigenfunctions, −∆ϕj,k = λjϕj,k, labeled by j ∈ N⁰ with k = 1, . . . , d^j counting the multiplicity and d^j being the dimension of the eigenspace corresponding to the eigenvalue λ^j > 0. Then, the set CH⁰ comprises all and only the vorticity fields of the form

$$\omega_j = \lambda_j \phi_j, \quad \phij = \sum{k=1}^{dj} a{j,k} \phi_{j,k},$$

for any j ∈ N⁰ and a^j = (aj,k) ∈ R ^d^j satisfying the energy constraint

$$\lambda_j a_j^2 = 2\mathcal{H}_0.$$

The energy constraint in particular fixes the length of the vector a^j but not its direction. The function S restricted to C^H⁰ is given by

$$S(\omega_j) = \mathcal{H}_0 \lambda_j$$ ,

which is bounded from below, since the spectrum of −∆ is bounded from below and H⁰ is a constant. The solutions of (13) correspond to the eigenfunctions with minimum eigenvalue (the ground states). Usually the eigenspace corresponding to the lowest eigenvalue is one-dimensional hence we have two solutions that differ only by the sign of the vorticity. In this example, the variational principle (13) picks the solution of (12) with the lowest λ.

In this paper, we use metriplectic dynamics in order to solve the variational principle (13). For this particular application it is essential that the orbit of the chosen metriplectic dynamical system relaxes completely to a constrained entropy minimum.

2.2.2. Axisymmetric MHD equilibria

A similar equilibrium problem arises from axisymmetric ideal magnetohydrodynamic (MHD) equilibria of electrically conducting fluids. For MHD, the general equilibrium condition with zero flow amounts to [10, 62, 63, 18, 6].

$$J \times B = c \nabla p$$ , $\operatorname{curl} B = 4\pi J/c$ , $\operatorname{div} B = 0$ , (16)

where B and J are vector fields and p is scalar field. Physically B and J are the magnetic field and the electric current density, respectively, while p is the fluid pressure. As noted above, Gaussian units are used with c being the speed of light in free space. The first equation, in particular, expresses the force balance between the Lorentz force J × B/c and the pressure gradient ∇p. The force balance implies the necessary conditions

$$B \cdot \nabla p = 0$$ , $J \cdot \nabla p = 0$ , (17)

that is, p is constant on the field lines of both B and J.

For axisymmetric solutions, i.e., solutions that have rotational symmetry around an axis, we introduce cylindrical coordinates (r, φ, z) around the symmetry axis, and from div B = 0 it follows that [5],

$$B = \chi \nabla \varphi + \nabla \psi \times \nabla \varphi,$$

$$4\pi J/c = \operatorname{curl} B = -\Delta^* \psi \nabla \varphi + \nabla \chi \times \nabla \varphi,$$

where ∆^∗ = r[∂r(r ⁻1∂r)] + ∂ z is a linear elliptic second-order differential operator in (r, z) coordinates, the Grad-Shafranov operator, whereas χ(r, z), p(r, z), and ψ(r, z) are real-valued scalar functions. The operator ∇ is the full three-dimensional gradient. Then axisymmetric equilibria with zero flow must satisfy the conditions

$$-\Delta^* \psi \nabla \psi - \chi \nabla \chi - 4\pi r^2 \nabla p = [\psi, \chi] r \nabla \varphi,$$

$$[\psi, p] = 0, \quad [\chi, p] = 0,$$

where the brackets [·, ·] are the canonical Poisson brackets in the (r, z)-plane, e.g., [χ, ψ] = r∇φ·(∇ψ×∇χ) = ∂rχ∂zψ−∂rψ∂zχ. The first equilibrium condition expresses the force balance, while the latter two follow from the necessary conditions (17), respectively. With homogeneous Dirichlet boundary conditions, we can formulate the problem

$$\begin{cases} u\nabla\psi - \chi\nabla\chi - 4\pi r^2\nabla p = 0, & \text{in } \Omega, \ [\psi, p] = 0, & [\psi, \chi] = 0, & \text{in } \Omega, \ -\Delta^*\psi - u = 0, & \text{in } \Omega, \ \psi = 0, & \text{on } \partial\Omega, \end{cases}$$ (18)

where Ω ⊂ R⁺ × R is a domain in the (r, z) plane satisfying r > 0 in the closure Ω (so that Ω is bounded away from the singularity of cylindrical coordinates r = 0). The auxiliary variable u is related to the φ component of the current density, since J^φ = J · ∇φ/|∇φ| = −c∆^∗ψ/(4πr) = cu/(4πr). Without specifying other constraints this problem is ill-posed in the same way as problem (11): the vanishing of the two Poisson brackets in (18) implies that, if ∇ψ, ∇p, and ∇χ are all nonzero, both χ and p are constant on the level sets of ψ. However, the functional relation between χ, p and ψ is undetermined. Fortunately, providing such information is straightforward. If we prescribe χ = √ λF(ψ) and p = λG(ψ) for given functions F, G ∈ C (R) and λ > 0, we obtain u = λ[(F ²/2)′ (ψ) + 4πr2G′ (ψ)] = λf(r, ψ) and the equilibrium problem reduces to

$$\begin{cases} -\Delta^* \psi = \lambda f(r, \psi), & \text{in } \Omega, \ \psi = 0, & \text{on } \partial\Omega, \end{cases}$$ (19)

which is referred to as the Grad-Shafranov equation, and is the analog of (12). (We remark that in realistic applications proper care must be taken to assign physically meaningful values of χ and p to different connected components of the ψ-contours.) Equation (19) is an "eigenvalue problem" for a semilinear elliptic equation and thus the same remarks about the well-posedness of Eq. (12) are valid here. Solving this eigenvalue problem is the standard way of computing axisymmetric MHD equilibria in tokamaks [64, 60, 65]. As for the Euler equations, the Grad-Shafranov problem can be considered solved, and it is used in this work as a benchmark problem.

As was the case of reduced Euler equilibria, solutions of the axisymmetric MHD equilibrium conditions (18) can be characterized by a variational principle of the form (1). For a state variable we choose u(r, z) = (4π/c)rJφ(r, z) defined over the bounded domain Ω ⊂ R⁺ × R, with r > 0 on Ω, and we consider the measure dµ = r ⁻1drdz on Ω. We assume that the profiles F, G ∈ C (R) are given so that f(r, y) := (F ²/2)′ (y)+ 4πr2G′ (y) satisfies f ′ y (r, y) ̸= 0. Therefore the map y 7→ f(r, y) is monotonic and has an inverse s ′ y (r, y) for any fixed r. After integration in y we find a function s ∈ C 2 (R⁺ × R), with s ′′ y (r, y) ̸= 0 and such that s ′ (r, ·) ⁻¹ = f(r, ·). Then we consider the problem

$$\min{\mathcal{S}(u)\colon \mathcal{H}(u) = \mathcal{H}_0},\tag{20}$$

with entropy and Hamiltonian

$$S(u) = \int{\Omega} s(r, u)d\mu$$ , $\mathcal{H}(u) = \frac{1}{2} \int{\Omega} |\nabla_{r,z}\psi|^2 d\mu$ , (21)

where ψ is regarded as a function of u given by the solution of the linear elliptic problem

$$-\Delta^* \psi = u, \quad \psi|_{\partial\Omega} = 0.$$ (22)

Here H amounts to the magnetic energy stored in the poloidal component ∇ψ × ∇φ of the magnetic field.

The functional derivatives are defined with respect to the L -product with the measure µ on Ω, so that

$$\frac{\delta \mathcal{S}(u)}{\delta u} = s'_y(r, u), \quad \frac{\delta \mathcal{H}(u)}{\delta u} = \psi.$$

Equation (9) gives

$$\begin{cases} s'_y(r, u) - \lambda \psi = 0, & \mathcal{H}(u) = \mathcal{H}_0, \ \text{with } \psi \text{ solution of (22)}. \end{cases}$$ (23)

The set C^H⁰ of constrained critical points is then

$$\mathfrak{C}_{\mathcal{H}_0} = {u \colon (u, \lambda) \text{ solves (23) for some } \lambda > 0},\$$

and among the elements of C^H⁰ , those with minimum entropy are solutions of the entropy principle (1).

Under the assumption on s(r, y), each element u ∈ C^H⁰ with λ ̸= 0 corresponds to an axisymmetric equilibrium. In order to see this, we have to find the fields p and χ corresponding to u. Given the profiles F and G from which s has been derived, let us define √

$$\chi = \sqrt{\lambda}F(\lambda\psi), \quad p = \lambda G(\lambda\psi),$$

together with the scaled variables ˜u = λu and ψ˜ = λψ. One can check that ˜u, ψ˜ together with χ and p given above solve conditions (18). We also have that ψ˜ is a solution of the Grad-Shafranov equation.

We remark that in general, the variational principle (20) can be formulated for generic profiles s(r, y) and as long as s ′′ y (r, y) ̸= 0, we can still find f(r, ·) = s ′ y (r, ·) ⁻¹ and a correspondence between (23) and (19). However, for general profiles, f(r, y) cannot be written in terms of F(y) and G(y), since f may not be a quadratic function of r; therefore, one cannot always find χ and p and thus an equilibrium in the sense of (18).

As in the case of the reduced Euler equations, for the application of metriplectic dynamics to the solution of (20) it is essential that the orbits of the metriplectic system completely relax to a constrained entropy minimum.

2.2.3. Beltrami fields

Let Ω ⊂ R ³ be a bounded simply connected domain, with sufficiently regular boundary ∂Ω, and let n: ∂Ω → R ³ be the outward unit normal to Ω. A vector field B : Ω → R 3 is called a Beltrami field if it satisfies the conditions

$$\begin{cases} (\operatorname{curl} B) \times B = 0, & \operatorname{div} B = 0, & \operatorname{in} \Omega, \ n \cdot B = 0, & \operatorname{on} \partial \Omega. \end{cases}$$ (24)

Solutions of (24) satisfy the ideal MHD equilibrium conditions (16) with constant pressure, hence Beltrami fields are force-free MHD equilibria [19].

A special class of Beltrami fields is given by the solutions of the eigenvalue problem for the curl operator: find (B, λ), λ ∈ R, such that

$$\begin{cases} \operatorname{curl} B = \lambda B, & \text{in } \Omega, \ n \cdot B = 0 & \text{on } \partial \Omega. \end{cases}$$ (25)

The eigenfunctions of this problem will be referred to as linear Beltrami fields and they are necessarily divergence-free. Since there is a countable family of eigenvalues of the curl operator, each corresponding to a finite-dimensional space of eigenfunctions [66, 67], both problem (24) and (25) are mathematically ill-posed due to nonuniqueness of the solution (although for (25), this is standard, since it is an eigenvalue problem).

Linear Beltrami fields can be characterized by a variational principle of the form (1), which has been proposed by Woltjer [45] and later applied to self-organized states in fusion plasmas by Taylor [46]. This variational principle is also central in multi-region relaxed MHD [68, 69]. Generalizations of Woltjer's variational principle have been proposed by Dixon and co-workers, including the free-boundary case [70].

For the formulation of the variational principle, let H⁰ ∈ R, and let V be the space of L ² vector fields u = B on Ω such div B = 0 and n · B = 0 on ∂Ω. Then we consider the problem

$$\min{\mathcal{S}(B)\colon \mathcal{H}(B) = \mathcal{H}_0},\tag{26}$$

with the entropy and Hamiltonian functions given by

$$S(B) = \frac{1}{2} \int_{\Omega} |B|^2 dx, \quad \mathcal{H}(B) = \frac{1}{2} Hm(B) := \frac{1}{2} \int{\Omega} A \cdot B dx, \tag{27}$$

where A is the vector potential for B defined as the unique solution of the problem

$$\begin{cases} \operatorname{curl} A = B, & \text{in } \Omega, \ \operatorname{div} A = 0, & \text{in } \Omega, \ n \times A = 0, & \text{on } \partial \Omega. \end{cases}$$ (28)

The necessary condition for solutions of (26) reads

$$\begin{cases} B - \lambda A = 0, & \mathcal{H}(B) = \mathcal{H}_0, \ \text{with } A \text{ given by (28)}, \end{cases}$$ (29)

where λ ∈ R is the Lagrange multiplier. The set of constrained critical points of the entropy function is

$$\mathfrak{C}_{\mathcal{H}_0} = {B \colon (B, \lambda) \text{ solves } (29) \text{ for some } \lambda \in \mathbb{R}}.$$

From (28) and (29), one has that each constrained critical point B ∈ CH⁰ is a linear Beltrami field since curl B = λ curl A = λB, and the same holds true for the corresponding potential, curl A = λA.

The formal analysis of the variational problem proceeds as in the example of the reduced Euler equations with a linear profile s ′ (y). Condition (29) implies that B ∈ V satisfies the eigenvalue problem (25). From the theory of this eigenvalue problem [66], we know that there is an orthonormal basis of eigenfunctions corresponding to the eigenvalue λ = λ^j ∈ R \ {0} labeled by j ∈ N⁰ and with k = 1, . . . , d^j counting the multiplicity. As in the case of the Euler equations, the set CH⁰ comprises all and only the vector fields B^j = λjA^j with A^j = P k aj,kAj,k, a^j = (aj,k)^k ∈ R dj , and that satisfy the constraint H(B) = H0, which is 2H0λ^j = |a^j | 2 . Different from the case of the reduced Euler equations, both H⁰ and the eigenvalues of the curl operator can be negative. Since |a^j | ² > 0, if H⁰ < 0 (resp. H⁰ > 0), only the eigenfunctions B^j with λ^j < 0 (resp. λ^j > 0) belong to C^H⁰ . The entropy function restricted to C^H⁰ amounts to

$$S(B_j) = \mathcal{H}_0 \lambda_j > 0,$$

and it is always positive even if H⁰ and λ^j can be negative. It is possible to show that the minimum is attained [71]. Hence the solutions of (26) are the ground states for the curl operator.

The variational principle (26) asks for divergence-free fields B with minimum energy (1/2)∥B∥ 2 ^L² = S(B) subject to the constraint that magnetic helicity Hm(B) = 2H(B) is held constant. Since magnetic helicity is a global constraint, i.e. H(B) takes values in R, λ is a constant in R and thus the variational principle selects linear Beltrami fields. Magnetic helicity H^m is related to the topology of the field lines of B [72]. Specifically magnetic helicity is the average asymptotic linking number of the field lines [73, 74]. Nonlinear Beltrami fields (24) on the other hand can be obtained from a variational principle with a much stronger constraint, i.e., Beltrami fields are energy minima constrained to configurations of the field B that are continuous deformation of a prescribed field B0. This constraint preserves the topology of the field lines, and thus magnetic helicity as well [75]. An overview of this variational principle is given in Appendix C for sake of completeness.

In this work, we discuss metriplectic dynamical systems on the space V of divergencefree fields. We address in particular convergence of an orbit to nonlinear Beltrami fields. This is an example for which complete relaxation of the orbit is not a desired property, since complete relaxation would leads to linear Beltrami fields. We shall address instead brackets that have a class of invariants much richer than just the Hamiltonian.

From a computational point of view, the direct numerical solution for Beltrami fields is possible by a variety of techniques [76, 69, 77]. Here, we shall not attempt to compare the performance of these methods with that of the metriplectic relaxation. Our aim is rather to study the convergence of metriplectic systems on a physically relevant problem. We note however that relaxation methods for force-free equilibria are common in several applications [19] and our study eventually aims at improving the rate of convergence of such methods.

3. On the relaxation of metriplectic systems

In this section, we present some remarks and mathematical results pertaining to equilibrium points of metriplectic systems, their stability, and sufficient conditions for convergence of nearby orbits.

Metriplectic dynamical systems dissipate entropy at constant energy, cf. (8). It is therefore central to ask whether a solution u(t) of a metriplectic system, defined for t ∈ [0, +∞), has a limit for t → +∞, and whether the limit, when it exists, is a minimum of entropy on the surface of constant Hamiltonian H(u) = H⁰ = H(u0), u⁰ = u(0) being the initial state. This is the variational principle (1). For some applications, e.g., the equilibrium problems introduced in Sections 2.2.1 and 2.2.2, complete relaxation of the orbit is essential. Recall, this means that the solution of the metriplectic system has a limit as t → +∞ and the limit is a solution of (1). In general for applications it is also useful to know the rate of convergence to the limit.

In this section we address implications of general metriplectic structure, i.e., properties (8), for the complete relaxation of the orbit. A first implication is that, since S is monotonically nonincreasing along an orbit, it is a candidate for a Lyapunov function [27, 78, 79]. This argument is standard for nondegenerate gradient flows, but adaptation is needed for metriplectic systems, which have degeneracy. In Section 3.1 we first recall some standard arguments for Lyapunov stability, and then discuss their adaptation to metriplectic systems. Then, in Section 3.2, we consider another classical tool in the theory of nondegenerate gradient flows, the Polyak– Lojasiewicz condition. We develop the details in the finite-dimensional setting, but make some comments on extension to infinite dimensions in Section 3.3.

3.1. Finite-dimensional systems: Lyapunov stability

Let us start by recalling the Lyapunov stability theorem for finite-dimensional dynamical systems. Consider a generic vector field X : Z → R ⁿ on a domain Z ⊆ R n, n ∈ N. We only assume that X is locally Lipschitz continuous (locally Lipschitz for short), that is, every point z⁰ ∈ Z has a neighborhood Uz⁰ where

$$|X(z) - X(z')| \le L_{z0}|z - z'|, \quad z, z' \in \mathcal{U}{z_0},$$

for a constant Lz⁰ > 0, possibly depending on z0. This assures existence and uniqueness for the ordinary differential equation system dz/dt = X(z).

Let z^∗ ∈ Z be an equilibrium point, i.e., X(z∗) = 0. A continuous function L: O → R defined on an open subset O ⊆ Z containing z∗, and differentiable in O \ {z∗} is a Lyapunov function for X, if

$$X(z) \cdot \nabla \mathcal{L}(z) \le 0, \quad z \in \mathcal{O} \setminus {z_*};$$ (L1)

$$\mathcal{L}(z*) = 0 \text{ and } \mathcal{L}(z) > 0, \quad z \neq z*.$$ (L2)

Such an L is a strict Lyapunov function for X, if it is a Lyapunov function and satisfies

$$X(z) \cdot \nabla \mathcal{L}(z) < 0, \quad z \in \mathcal{O} \setminus {z_*}.$$ (L3)

Condition (L2) in particular implies that z^∗ is an isolated minimum of L in O.

The the Lyapunov stability theorem [27, 53, 80, 81] states that, if a Lyapunov function exists, then z^∗ ∈ O is a stable equilibrium point. By definition this means that for any ε > 0 there is δ > 0, such that |z⁰ − z∗| < δ implies ^z(t) ⁻ ^z^∗ < ε for all ^t ^≥ 0, with z(t) being an integral curve of X with initial condition z(0) = z0. In addition, if the Lyapunov function is strict, z^∗ is an asymptotically stable equilibrium point, that is, z^∗ is a stable equilibrium point, in the sense defined above, and δ > 0 can be chosen so that limt→+^∞ z(t) = z∗, for all initial conditions z⁰ satisfying |z⁰ − z∗| < δ.

Although the Lyapunov stability theorem is well known [27, 53, 80, 81], we recall the proof for sake of completeness, since the other results in Section 3 rely on the same ideas. The various arguments available in the literature differ essentially only in the final step in the proof of asymptotic stability. Here we follow Moretti [81], which we find particularly clear. Recall that Br(z) = {z ′ ∈ R ⁿ : |z ′ − z| < r} denotes the open ball of radius r > 0 in R n.

Theorem 1 (Lyapunov stability). Let X : Z → R ⁿ be a locally Lipschitz vector field, z^∗ ∈ Z an equilibrium point of X, O ⊆ Z an open subset containing z∗, and L: O → R continuous in O and differentiable in O \ {z∗}.

(i) If L is a Lyapunov function for X, then z^∗ is a stable equilibrium point.
(ii) If L is a strict Lyapunov function for X, then z^∗ is an asymptotically stable equilibrium point.

Proof. (i) Stability. Since O is open, we can choose ε > 0 so small that the ball Bε(z∗) is contained in the neighborhood O. On the boundary ∂Bε, L(z) > 0 because of (L2). Let L^ε := minz∈∂Bε(z∗) L(z). The minimum exists since ∂Bε(z∗) is compact. Let us now choose δ > 0 so small that L(z) < L^ε for z ∈ Bδ(z∗). This is possible since L is continuous and L(z∗) = 0. For any z⁰ ∈ Bδ(z∗), let z(t), t ∈ [−τ0, +τ0] be an integral curve of the field X with initial condition z(0) = z0. Assumption (L1) implies that,

$$\mathcal{L}(z(t)) \leq \mathcal{L}(z0) < \mathcal{L}{\varepsilon}$$ , for all $t \in [0, \tau_0]$ .

(ii) Asymptotic stability. If L is a strict Lyapunov function, then in particular, it is a Lyapunov function and thus z^∗ is a stable equilibrium point. We can choose ε⁰ > 0 and a δ⁰ > 0 such that any integral curve z(t) with z(0) ∈ B^δ⁰ (z∗) stays in B^ε⁰ (z∗).

We want to show that for any z⁰ ∈ B^δ⁰ (z∗) and for any ε ∈ (0, ε0) there is a time T^ε > 0 such that the integral curve z(t) with initial condition z(0) = z⁰ satisfies |z(t)−z∗| < ε for all t > Tε.

Uniqueness of the orbit passing through a given point implies that, if z⁰ ̸= z∗, then |z(t) − z∗| > 0 (since z(t) = z^∗ is a solution) and thus dL z(t) /dt < 0 for all t > 0. Therefore function t 7→ L z(t) is strictly monotonically decreasing, and it is bounded form below, hence it has a limit

$$\mathcal{L}(z(t)) \to \ell = \inf_{\substack{t \ge 0 \ 19}} \mathcal{L}(z(t)) \ge 0,$$

The limit must be ℓ = 0. If not, L z(t) ≥ ℓ > 0 for all t > 0, and continuity of L implies that there a radius r ∈ (0, ε0) such that |z(t) − z∗| > r. This leads to a contradiction, since, if z(t) stays in the compact region R = {z : r ≤ |z − z∗| ≤ ε} for all t ≥ 0, then with −M = maxz∈^R X(z) · ∇L(z) < 0,

$$\mathcal{L}\big(z(t)\big) = \mathcal{L}(z_0) + \int_0^t X\big(z(s)\big) \cdot \nabla \mathcal{L}\big(z(s)\big) ds \le \mathcal{L}(z_0) - Mt.$$

For t > L(z0)/M > 0, we have L z(t) < 0, which is impossible. Hence, the limit must be ℓ = 0, that is, for every λ > 0 there is T^λ > 0 such that L z(t) < λ for t > Tλ.

We claim that for any ε > 0 we can find λ = λ^ε > 0 such that

$$A{\lambda} := { z \in B{\varepsilon0}(z*) \colon \mathcal{L}(z) < \lambda },$$

is contained in the ball Bε(z∗), i.e. A^λ ⊂ Bε(z∗). If this is the case, corresponding to λε, there is a time T^ε > 0 such that, for all t > Tε, z(t) ∈ Aλ^ε ⊂ Bε(z∗), which is the thesis. Therefore it remains to prove the claim. Let us assume that the claim is false, i.e. there is a value ε^∗ > 0 such that, for all λ, there is at least one point ˜z^λ that satisfies the conditions L(˜zλ) < λ and ε^∗ ≤ |z˜^λ − z∗| ≤ ε0. Upon choosing λ = 1/n for n ∈ N we obtain a sequence zⁿ = ˜z1/n, which belongs to a compact set. Hence there is a converging subsequence zn^k → z˜∗, and ε^∗ ≤ |z˜^∗ − z∗| ≤ ε⁰ so that necessarily L(˜z∗) > 0. On the other hand, we have L(zn^k ) < 1/n^k → 0 and, by continuity of L, L(˜z∗) = 0, which is a contradiction.

We want to apply Theorem 1 to metriplectic vector fields of the following form:

$$X(z) = J(z)\nabla \mathcal{H}(z) - K(z)\nabla \mathcal{S}(z). \tag{30}$$

For comparison, we also address the case of a standard nondegenerate gradient flows

$$X(z) = -\nabla S(z),$$ (31)

with entropy S ∈ C 2 (Z).

First, recall that in the case of nondegenerate gradient flows (31), with S ∈ C 2 (Z), the function L(z) = S(z) − S(z∗) is a strict Lyapunov function in a neighborhood of any isolated, local minimum z^∗ of S [80, Proposition 15.0.2]. Therefore, the Lyapunov stability theorem implies that any integral curve of the gradient flow with initial condition near an isolated entropy minimum converges to that entropy minimum for t → +∞ (asymptotic stability). This result is a local version of the property we called complete relaxation.

Now, consider the case of a metriplectic vector field (30). If z^∗ ∈ Z is an equilibrium point, the function L(z) = S(z) − S(z∗) satisfies condition (L1), because of the general properties of metriplectic systems, cf. (8). If the entropy function S has an isolated minimum at z^∗ ∈ Z, L(z) is a Lyapunov function of the system in a neighborhood O of z∗. Therefore z^∗ is a stable equilibrium. However, an orbit z(t) can converge to z^∗ only if the initial condition z⁰ = z(0) and the local entropy minimum z^∗ belong to the same energy isosurface, i.e., H(z0) = H(z∗). This is a necessary condition that follows from the continuity of H: if z(t) → z∗, passing to the limit in H(z0) = H(z(t)) yields H(z0) = H(z∗). This observation leads to the following conclusion.

Proposition 2. Let X : Z → R ⁿ be a locally Lipschitz vector field on a domain Z ⊆ R n, z^∗ ∈ Z an equilibrium point of X, and H, S ∈ C (Z) such that

1. X(z) · ∇H(z) = 0 and X(z) · ∇S(z) ≤ 0, z ∈ Z;
2. ∇H(z∗) ̸= 0.

Then, if L = S − S(z∗) satisfies (L2) in a neighborhood O ⊆ Z of z∗, there is at least one point z ′ ∈ O \ {z∗} such that X(z ′ ) · ∇S(z ′ ) = 0.

Proof. By contradiction, let us assume that X · ∇S < 0 in O \ {z∗}. Then hypothesis 1 implies that L(z) = S(z) − S(z∗) satisfies the conditions for a strict Lyapunov function for X. For the Lyapunov stability theorem, there exists δ ′ > 0 such that for any z⁰ with |z⁰ − z∗| < δ′ the orbit z(t) of the dynamical system dz/dt = X(z) with initial condition z(0) = z⁰ exists for all t ≥ 0 and z(t) → z^∗ as t → +∞. Hypothesis 1 also implies that H is a constant of motion and it is continuous, therefore H(z0) = H(z∗).

From hypothesis 2 and the continuity of the derivative ∇H we can find a ball of radius δ ′′ > 0 around z^∗ where ∇H ̸= 0.

We choose δ < min{δ ′ , δ′′}. In the ball of radius δ centered at z^∗ there is at least one point z⁰ such that H(z0) ̸= H(z∗). If not, then H is constant in the ball, and this is not possible since ∇H ̸= 0. On the other hand, since |z⁰ − z∗| < δ < δ′ we must have H(z0) = H(z∗), which is a contradiction.

For metriplectic systems, hypothesis 1 is verified, cf. (8). Hypothesis 2 holds away from critical points of H, which are usually isolated. Therefore, this proposition implies that the entropy S of a metriplectic system in most cases (specifically under hypothesis 2) cannot be used as a strict Lyapunov function, since it is not strictly decaying everywhere in O \ {z∗}. Without a strict Lyapunov function, asymptotic stability of local entropy minima does not follows directly from the usual Lyapunov theorem. This is in stark contrast with the case of pure gradient flows (31) discussed above.

Asymptotic stability with Lyapunov functions that are not strictly dissipated have been addressed by De Salle and Lefschetz [82], who showed that any integral curve, defined for t ≥ 0, in a bounded strict sublevel set of the Lyapunov function must approach the largest invariant set contained in the region where X(z) · ∇L(z) = 0. More specific results for systems with a conserved energy and a dissipated entropy were considered by Beretta [83] with applications to quantum thermodynamics. Here we apply similar ideas to metriplectic systems. Suppose a locally Lipschitz vector field X : Z → R ⁿ on a domain Z ⊆ R ⁿ has k constants of motion I 1 , . . . , I ^k ∈ C^∞(Z), for 1 ≤ k < n, which means X(z) · ∇Iα(z) = 0, for α ∈ {1, . . . , k}. The functions I ^α are independent at z ∈ Z if

$$\nabla \mathcal{I}^1(z), \dots, \nabla \mathcal{I}^k(z)$$ are linearly independent in $\mathbb{R}^k$ . (32)

It is convenient to define the function I := I , . . . , I k ∈ C^∞(Z, R k ). Then (32) is equivalent to rank ∇I(z) = k.

If the constants of motion are independent at an equilibrium point z^∗ ∈ Z, where X(z∗) = 0, then the local submersion theorem [53, Theorem 2.5.13] allows us to find a neighborhood U of z^∗ where the level sets

$$\mathcal{U}_{\eta} := { z \in \mathcal{U} \colon \mathcal{I}(z) = \eta },\$$

are closed submanifolds of U, for η is a neighborhood of η^∗ := I(z∗). Since I ^α are constants of motion, X is tangent to Uη, so that the dynamical system can be reduced locally to each submanifold Uη. This leads to the following straightforward result, which is essentially a finite-dimensional version of Theorem 2 in [83].

Proposition 3. Let X : Z → R ⁿ be a locally Lipschitz vector field on a domain Z ⊆ R n, I = (I 1 , . . . , I k ) ∈ C∞(Z, R k ) be constants of motion with 1 ≤ k < n, z^∗ ∈ Z be an equilibrium point of X, rank ∇I(z∗) = k, η^∗ := I(z∗), and let U^η^∗ = {z ∈ U : I(z) = η∗}, where U is the neighborhood of z^∗ given by the local submersion theorem. If there is L ∈ C (U) satisfying

$$X(z) \cdot \nabla \mathcal{L}(z) \le 0,$$ $z \in \mathcal{U},$ (L1')

$$\mathcal{L}(z*) = 0 \text{ and } \mathcal{L}(z) > 0,$$ $z \in \mathcal{U}{\eta*} \setminus {z*},$ (L2')

then for any sufficiently small ε > 0, there is a δ > 0 such that the integral curve z(t) of X with initial condition z(0) = z⁰ ∈ Bδ(z∗) ∩ U^η^∗ is defined for all t ≥ 0, and z(t) ∈ Bε(z∗) ∩ Uη^∗ . If in addition

$$X(z) \cdot \nabla \mathcal{L}(z) < 0, \quad z \in \mathcal{U}{\eta} \setminus {z_},$$ (L3')

then z(t) → z^∗ for t → +∞.

Proof. Since rank ∇I(z∗) = k, the local submersion theorem [53] allows us to find open subsets N ⊆ I(Z) and Z ′ ⊆ R ⁿ−^k ∼= ker ∇I(z∗), together with local coordinates (η, ζ) ∈ N × Z′ , defined in an open, connected subset U ⊆ Z containing the point z∗, such that the inverse coordinate map φ: N × Z′ → U is a C^∞-diffeomorphism and satisfies

$$I(\varphi(\eta, \zeta)) = \eta.$$

Therefore, the sets U^η = {z ∈ U : I(z) = η}, with η ∈ N , are closed submanifolds parameterized by ζ 7→ φ(η, ζ). Let us denote by (η∗, ζ∗) ∈ N ×Z′ the point corresponding to z^∗ = φ(η∗, ζ∗).

In this local coordinate system the integral curves of the vector field X in U solve

$$\frac{d\eta}{dt} = X \cdot \nabla \eta = 0, \quad \frac{d\zeta}{dt} = X_{\eta}(\zeta) := (X \cdot \nabla \zeta) \circ \varphi(\eta, \zeta),$$

since X · ∇η = X · ∇I = 0. The field X^η defines a tangent vector on U^η for any η ∈ N . Since ∇ζ is C^∞, one can check that X^η is a locally Lipschitz continuous function of ζ. Then, with L^η := L ◦ φ(η, ·), one has

$$(X \cdot \nabla \mathcal{L}) \circ \varphi = X{\eta}^{a} \frac{\partial \mathcal{L}{\eta}}{\partial \zeta^{a}},$$

where the sum over a ∈ {1, . . . , n − k} is implied. Hence hypotheses (L1′ )-(L3′ ) are equivalent to Lyapunov conditions (L1)-(L3) for the function L^η^∗ and for the dynamical system dζ/dt = X^η^∗ on U^η^∗ , and the claim follows from Theorem 1.

Proposition 3 establishes stability and asymptotic stability for orbits with initial conditions in U^η^∗ , which being a lower dimensional set has zero Lebesgue measure in the phase space. This issue was addressed in [83] by assuming z^∗ is part of a continuous family of equilibria.

At last, we address metriplectic vector fields. We know that there is at least one constant of motion, the Hamiltonian H. More generally, let X be a metriplectic vector field with k constants of motion I 1 , . . . , I k , 1 ≤ k < n, and for definiteness I ^k = H. Differently from the case of Proposition 3, we assume that (32) holds in a whole subdomain U. Under these conditions, the submersion theorem [53, Theorem 3.5.4] establishes that the sets

$$\mathcal{U}_{\eta} := { z \in \mathcal{U} : \mathcal{I}(z) = \eta }, \tag{33}$$

is a closed submanifold of U for any η ∈ I(U), and the restriction of the entropy to such submanifolds, i.e., S|^U^η , is smooth.

We make the following hypothesis on the entropy in U:

$$\mathcal{U}$$ is bounded, $\overline{\mathcal{U}} \subset \mathcal{Z}$ , $\exists z_m \in \mathcal{U}$ , where $\mathcal{S}(z_m) < \inf{{\mathcal{S}(z) : z \in \partial \mathcal{U}}}$ . (34)

This ensures that there is a minimum of the entropy in the interior of the subdomain U.

Lemma 4. Let X : Z → R d be a vector field on a domain Z ⊆ R d , and let U ⊂ Z and S ∈ C (Z) satisfy (34). If in the subdomain U, X is locally Lipschitz and X·∇S ≤ 0, then there is an open, non-empty subset O ⊆ U such that, for any z⁰ ∈ O the integral curve z(t) of X with initial condition z(0) = z⁰ can be prolonged to the interval t ∈ [0, +∞) and z(t) ∈ O for all t ≥ 0.

Proof. Per hypotheses S ∈ C 1 (Z), therefore S ∈ C 1 (U), and we have S^b := inf∂^U S > S(zm) > −∞. The set

$$\mathcal{O} := {z \in \mathcal{U} : \mathcal{S}(z) < \mathcal{S}_b}$$

is non-empty, since z^m ∈ O, and open, since S(z) is continuous. Let z : (τ1, τ2) → U, τ¹ < 0 < τ2, be the maximal solution of the initial value problem

$$dz/dt = X(z), \quad z(0) = z_0 \in \mathcal{O}.$$

The solution exists given that X is locally Lipschitz continuous in U. Since X · ∇S ≤ 0, the function t 7→ S z(t) is C ¹ and non-increasing, hence S z(t) ≤ S(z0) < Sb, and thus z(t) ∈ O for all t ∈ [0, τ2).

We show that τ² = +∞. With this aim we rely on a standard argument from the theory of ordinary differential equations, which we report in full for sake of completeness. If τ² is finite, let {tn}n∈^N be a sequence, tⁿ ∈ (τ1, τ2), and tⁿ → τ² as n → +∞. Then

$$|z(t_n) - z(tm)| \le \max{z \in \overline{U}} |X(z)| |t_n - t_m|,$$

for all m, n ∈ N. Since {tn} is a convergent sequence, {z(tn)} ⊂ O is a Cauchy sequence and must have a limit ¯z ∈ O as n → ∞. Passing to the limit in the inequality S z(tn) ≤ S(z0) < S^b yields S(¯z) ≤ S(z0) < Sb, hence ¯z ∈ O. We can use the limit point ¯z as an initial condition, we can extend the solution z(t) beyond τ² contradicting the fact that z(t) is the maximal solution. Therefore τ² = +∞ as claimed.

The set O is given by

$$\mathcal{O} = \left{ z \in \mathcal{U} \colon \ \mathcal{S}(z) < \inf_{\partial \mathcal{U}} \mathcal{S} \right}$$

This lemma shows that, if S and U satisfy (34) and S is non-increasing along the integral curves of a vector field X, then we can find a positively invariant subset O. This lemma can be applied directly to both standard gradient flows and metriplectic systems.

Proposition 5. Let X = J∇H−K∇S be a metriplectic vector field on a domain Z ⊆ R n, with 1 ≤ k < n constants of motion I = (I 1 , . . . , I k ) ∈ C∞(Z, R k ), rank ∇I = k in a subdomain U satisfying (34), and let O be the open set given by Lemma 4. If in addition it holds that

$$\forall \eta \in \mathcal{I}(\mathcal{O}), \text{ on } \mathcal{U}{\eta} \coloneqq {z \in \mathcal{U} \colon \mathcal{I}(z) = \eta}, \, \mathcal{S}|{\mathcal{U}{\eta}} \text{ has a unique critical point } z{\eta},$$

which is a strict local minimum,

and

$$\ker K = \operatorname{span} { \nabla \mathcal{I}^1, \dots, \nabla \mathcal{I}^k } \text{ in } \overline{\mathcal{U}}.$$ (L3")

then for any integral curve z(t) of X with z(0) = z⁰ ∈ O, limt→+^∞ z(t) = z^η ∈ O, where η = I(z0).

Proof. First we show that z^η ∈ Uη∩O. Given that η ∈ I(O), the set Uη∩O is non-empty. If z^η ̸∈ Uη∩O, there cannot be any extremum of S|^U^η in Uη∩O since z^η is the only critical point in Uη. This implies that both the minimum and the maximum of S restricted to Uη∩O are attained on the boundary Uη∩∂O, where we have S = S^b = inf{S(z): z ∈ ∂U}. We deduce that S|^U^η is constant and equal to S^b in U^η ∩ O, but per definition, points in O satisfy S < Sb. Therefore, the only critical point z^η of S|^U^η must be in U^η ∩ O.

The function S|^U^η is smooth and its critical points necessarily satisfy the Lagrange condition

$$\nabla S(z) = \sum{\alpha=1}^{k} \lambda{\alpha} \nabla \mathcal{I}^{\alpha}(z), \quad \mathcal{I}(z) = \eta, \tag{35}$$

for λ^α ∈ R [53, Theorem 3.5.27]. Hypothesis (L2′′), states, in particular, that this can only happen at the point zη.

Condition (L2′′) also requires z^η to be a strict local minimum of S|^U^η . Then, z^η must be an equilibrium point of X. In order to see this, let us consider the integral line zˆη(t) ∈ O of the vector field X with initial condition ˆzη(0) = zη. Given that z^η ∈ O, ˆzη(t) is defined for t ≥ 0 (Lemma 4). Since I is a constant of motion, ˆzη(t) ∈ Uη, and continuity implies that for any ε > 0, there is δ > 0 such that |t| ≤ δ implies |zˆη(t) − zη| < ε. If zˆ^η is not identically equal to z^η when t ∈ [0, δ), then at some t ′ ∈ (0, δ), ˆzη(t ′ ) ̸= z^η and S zˆη(t ′ ) > S(zη) since z^η is an entropy local minimum. On the other hand, dissipation of entropy requires S zˆη(t ′ ) ≤ S zˆη(0) = S(zη). Therefore, it must be ˆzη(t) = z^η identically for t ∈ [0, δ), which implies X(zη) = 0 and z^η is an equilibrium.

Given z⁰ ∈ O, the integral curve of X with initial condition z(0) = z⁰ exists for all t ≥ 0 (Lemma 4). Let us define the function

$$\mathcal{L}(z) = \mathcal{S}(z) - \mathcal{S}(z_{\eta}),$$

where η = I(z0). Then ∇L(z) = ∇S(z) and

$$X(z) \cdot \nabla \mathcal{L}(z) = -\nabla \mathcal{S}(z) \cdot K(z) \mathcal{S}(z) \le 0,$$ (36)

so that L satisfies condition (L1′ ) of Proposition 3 in the open set O defined above.

Since z^η is a local minimum of S|^U^η there is a neighborhood of z^η on U^η where S(z) > S(zη) for z ̸= zη. But by definition of the submanifold topology, this neighborhood has the form U^η ∩ V with V an open subset of U. Then, condition (L2′ ) holds true with U replaced by V.

At last we show that condition (L3′ ) is also satisfied. If not, one could find a point z˜^η ∈ U^η such that ˜z^η ̸= z^η and ∇S(˜zη) ∈ ker K(˜zη), but then assumption (L3′′) implies that ˜z^η satisfies the Lagrange conditions (35) and this contradicts the uniqueness of the critical point zη. Therefore, inequality (36) must be strict, which is condition (L3′ ).

In summary, the point z^∗ = zη, η = I(z0), and the function L(z) = S(z)−S(zη) on the neighborhood V satisfy all the hypotheses of Proposition 3, including (L3′ ). Therefore, there is a δ > 0 such that any integral curve of X with initial condition in U^η ∩ Bδ(zη) converges to z^η as t → +∞.

We claim that, for any integral curve, z(t) ∈ O, t ≥ 0 and z(0) = z⁰ ∈ U^η ∩ O, there must be a time t^δ such that z(tδ) ∈ U^η ∩ Bδ(zη). If this is the case, then necessarily z(t) → zη, which is the thesis. Therefore it remains to prove this last claim.

The sets U^η and U^η \Bδ(zη) are compact in R ⁿ. In view of conditions (L2′′) and (L3′′), on the set U^η \Bδ(zη) we have X·∇S < 0. Let −M := sup{X(z)·∇S(z): z ∈ U^η \Bδ(zη)}. Then M > 0 since the supremum is attained, and if the integral curve z(t) is such that z(t) ∈ U^η \ Bδ(zη) for all t ≥ 0,

$$S(z(t)) = S(z_0) + \int_0^t (X \cdot \nabla S)(z(s)) ds \le S(z_0) - Mt.$$

This leads to a contradiction as in the proof of the classical Lyapunov stability theorem. Hence there must be a time t^δ at which z(tδ) ∈ U^η ∩ Bδ(zη).

Condition (L2′′) ensures that (L2′ ) holds for the function L(z) = S(z) − S(zη) on Uη. The uniqueness of the critical point of S together with (L3′′) ensures that entropy is strictly dissipated by requiring that the metric bracket generated by the tensor K is "specifically degenerate", by which we mean it preserves the invariants I ^α only. We use the terminology "minimally degenerate" to mean that the only degeneracy is that associated with H.

We remark that the limit point z^η ∈ U^η is the unique solution of the optimization problem

$$\min{\mathcal{S}(z)\colon z\in\mathcal{O},\ \mathcal{I}(z)=\mathcal{I}(z_0)},\$$

hence Proposition 5 establishes a local version of the desired complete-relaxation property for all orbit starting in the set O.

3.2. Finite-dimensional systems: the Polyak– Lojasiewicz inequality

Nondegenerate gradient flows and gradient descent methods have been studied under the assumptions that the entropy function S satisfies the classical Polyak– Lojasiewicz (PL) condition [84, 85] and ∇S is Lipschitz continuous, the latter being a natural hypothesis.

A function S ∈ C 1 (Z) satisfies the PL condition in a non-empty subset Q ⊆ Z if S^∗ := inf{S(z): z ∈ Q} > −∞ and there exists a constant κ > 0 such that

$$\frac{1}{\kappa} |\nabla \mathcal{S}(z)|^2 \ge \mathcal{S}(z) - \mathcal{S}_*, \quad z \in Q.$$ (PL)

If S satisfies the PL inequality, the associated gradient flow has the following properties:

• If z(t) ∈ Q, t ∈ [0, +∞), is an integral curve of the gradient flow (31), then S z(t) converges exponentially to S^∗ as t → +∞. This follows from

$$\frac{d}{dt} \left[ \mathcal{S} \big( z(t) \big) - \mathcal{S}* \right] = - \left| \nabla \mathcal{S} \big( z(t) \big) \right|^2 \le -\kappa \left[ \mathcal{S} \big( z(t) \big) - \mathcal{S}* \right]$$

which implies

$$\left[\mathcal{S}(z(t)) - \mathcal{S}_\right] \le \left[\mathcal{S}(z0) - \mathcal{S}\right]e^{-\kappa t},$$

where z⁰ ∈ Q is the initial condition at t = 0. The exponential rate of convergence is given by the constant κ in the inequality. Unfortunately convergence of the entropy values alone does not directly imply convergence of the orbit z(t) to a limit for t → +∞ (cf. the counterexample of Palis and de Melo [86]).

• All equilibrium points of the gradient flow (31) are global minima of S (not just critical points). In fact, if z^e ∈ Q is an equilibrium point, ∇S(ze) = 0 and the PL condition implies 0 = |∇S(ze)| ² ≥ κ S(ze) − S^∗ ≥ 0, hence S(ze) = S∗.

Any strongly convex function satisfies condition (PL), cf. the short proof reported below, and thus are included as a special case. In general, (PL) is weaker than strong convexity. In fact, it has been shown that the classical PL condition is weaker than several other conditions introduced in order to address the convergence of gradient descent methods [87]. In addition, no assumption is made on the global minimum of S, which, in particular, does not need to be an isolated point. For instance, the function S(z) = (|z| ²−1)² satisfies the PL condition with κ = 16r 2 0 in the domain S = {z : |z| > r0} for any r⁰ ∈ (0, 1/2), and attains its minimum on the sphere |z| = 1; hence there is no isolated minimum. On the other hand, this function does not satisfies the PL condition on the whole space R ⁿ, because of the critical point at z = 0, which is a local maximum.

Polyak [84] under the additional (natural) hypothesis that the vector field X = −∇S is Lipschitz continuous (not just locally Lipschitz), established convergence of the gradient flow trajectories to the global entropy minimum. Specifically, Polyak's result in the notation used here amounts to the following.

Theorem 6 (Polyak 1963, Theorem 9 [84]). Let z⁰ ∈ Z, ρ, κ, L > 0, and S ∈ C 1 (Z) be such that (PL) is satisfied and ∇S is Lipschitz continuous in the closed ball Bρ(z0) ⊂ Z with Lipschitz constant L, and γ = √ 8Lφ0/(ρκ) ≤ 1, where φ⁰ = S(z0)−S∗. Then there is z^∗ ∈ Bγρ(z0) such that the integral curve z(t) of the gradient flow (31) with z(0) = z⁰ is defined for all t ≥ 0 and ^z(t) ⁻ ^z^∗ ^≤ γρe⁻κt/²

We shall now give the details concerning the Polyak– Lojasiewicz condition. Although well-known, we recall for convenience of the reader the short proof of the fact that a strongly convex entropy satisfies inequality (PL).

Proof: Strongly convex imply PL. A strongly convex function S ∈ C (Z), Z ⊆ R ⁿ with parameter α > 0 satisfies

$$S(z') \ge S(z) + (z'-z) \cdot \nabla S(z) + \frac{\alpha}{2} |z'-z|^2,$$

for any z, z′ ∈ Z. The right-hand side is bounded from below by its infimum over z ′ which is attained at z ′ = z − α ⁻¹∇S(z) and thus

$$S(z') \ge S(z) - \frac{1}{2\alpha} |\nabla S(z)|^2.$$

This inequality holds for any z ′ , and thus implies that S^∗ = inf{S(z): z ∈ Z} > −∞. Taking the infimum over z ′ yields

$$S_* \ge S(z) - \frac{1}{2\alpha} |\nabla S(z)|^2$$ .

This can be rearranged to give the PL inequality with constant κ = 2α.

In Polyak's original formulation of the theorem, the size ρ of the ball is determined by the condition γ ≤ 1 in terms of the entropy at the initial position z0. For instance, S(z) = z ²/2 satisfies the hypothesis with κ = 2 and L = 1; given any z⁰ ∈ R ⁿ, γ ≤ 1 is equivalent to ρ ≥ |z0|, hence Bρ(z0) is large enough to contain the unique global minimum z^∗ = 0. We give a slightly different statement of Polyak's result.

Proposition 7. Let S ∈ C 1 (Z), X = −∇S Lipschitz continuous in a a subdomain U ⊂ Z satisfying (34), and let O be the open, non-empty subset given by Lemma 4. If S satisfies (PL) in U, then for any integral curve z(t) of X with initial condition z(0) = z⁰ ∈ O there is z^∗ ∈ O and a constant θ > 0 depending on z⁰ such that S(z∗) = S∗,

$$\left|z(t)-z*\right| \leq \theta e^{-\kappa t/2}\,, \quad \left[\mathcal{S}!\left(z(t)\right)-\mathcal{S}\right] \leq \left[\mathcal{S}!\left(z0\right)-\mathcal{S}\right] e^{-\kappa t}\,, \quad \textit{for } t \geq 0\,.$$

Next we give a proof of Proposition 7, which is based on the original argument due to Polyak [84]. The only difference consists in the use of Lemma 4 to establish the existence of a global solution. We give the proof in details since similar ideas are then needed below for the metriplectic case.

Proof of Proposition 7. Lemma 4 ensures that for any z⁰ ∈ O there is an integral curve z(t) ∈ O of the gradient flow X = −∇S passing through z(0) = z⁰ and this is defined for all t ≥ 0. As a consequence of (PL), φ(t) = S z(t) − S^∗ decays to zero exponentially, i.e., φ(t) ≤ φ0e ⁻κt, with φ⁰ = φ(0).

If there is a finite time 0 ≤ t <¯ +∞ at which φ(t¯) = 0, then z(t¯) is a minimum for the entropy S and thus ∇S z(t¯) = 0 so that z(t¯) is an equilibrium point. We deduce z(t) = z(t¯) = z⁰ for all t ≥ 0 and the thesis holds true with z^∗ = z0.

As for the non-trivial case φ(t) > 0 for all t ≥ 0, we distinguish two key steps. Step 1. For any t1, t² ≥ 0, t¹ < t2, we have

$$S(z(t_1)) - S(z(t2)) = \int{t_1}^{t_2} |\nabla S(z(t))|^2 dt.$$

Following Polyak, we estimate the right-hand side from below. First the triangular inequality gives

The second term on the right-hand side is bounded by

$$\left|\nabla S(z(t)) - \nabla S(z(t_1))\right| \le L|z(t) - z(t_1)|,$$

where L > 0 is the Lipschitz constant of ∇S. In addition,

$$\begin{split} \frac{d}{dt} |z(t) - z(t_1)| &\leq \left| \frac{d}{dt} |z(t) - z(t_1)| \right| \leq \left| \nabla \mathcal{S} \big( z(t) \big) \right| \ &\leq \left| \nabla \mathcal{S} \big( z(t_1) \big) \right| + \left| \nabla \mathcal{S} \big( z(t) \big) - \nabla \mathcal{S} \big( z(t_1) \big) \right| \ &\leq \left| \nabla \mathcal{S} \big( z(t_1) \big) \right| + L |z(t) - z(t_1)|. \end{split}$$

The Gr¨onwall's inequality then gives

$$\left|z(t) - z(t_1)\right| \le \frac{1}{L} \left|\nabla \mathcal{S}(z(t_1))\right| \left[e^{L(t-t_1)} - 1\right],$$

and thus

$$\left|\nabla S(z(t)) - \nabla S(z(t_1))\right| \le \left|\nabla S(z(t_1))\right| \left[e^{L(t-t_1)} - 1\right].$$

If L(t² − t1) ≤ log 2, the term on the right-hand side is ≤ ∇S z(t1) and thus

$$\left|\nabla S(z(t))\right| \ge \left|\nabla S(z(t_1))\right| \left[2 - e^{L(t-t_1)}\right].$$

With t² = t¹ + (log 2)/L and t¹ arbitrary, we obtain

$$S(z(t_1)) - S(z(t_2)) \ge \frac{1}{\alpha L} |\nabla S(z(t))|^2$$

where

$$\frac{1}{\alpha} \coloneqq \int_0^{\log 2} \left[2 - e^s\right]^2 ds$$

is a positive numerical constant. At last, we use the fact that S z(t) is non-increasing and t¹ < t² so that

$$S(z(t_1)) - S(z(t_2)) = [S(z(t1)) - S] - [S(z(t2)) - S]$$

$$\leq [S(z(t1)) - S*],$$

and deduce that, for any t¹ ≥ 0,

$$\left|\nabla \mathcal{S}(z(t_1))\right|^2 \le \alpha L\left[\mathcal{S}(z(t1)) - \mathcal{S}*\right] = \alpha L\varphi(t_1). \tag{37}$$

We have established that, under the hypotheses, ∇S along the orbit is controlled by the entropy decay.

Step 2. Given arbitrary points in time t1, t² ≥ 0, t¹ < t2,

$$\left|z(t_1) - z(t2)\right| \le \int{t_1}^{t_2} \left|\nabla S(z(t))\right| dt,$$

and (37) together with φ(t) ≤ φ0e ⁻κt yields

$$|z(t_1) - z(t_2)| \le \frac{\sqrt{4\alpha L\varphi_0}}{\kappa} \left[ e^{-\kappa t_1/2} - e^{-\kappa t_2/2} \right].$$ (38)

This inequality can be used to show that for any sequence {tn}n∈^N with tⁿ → +∞ as n → +∞, z(tn) is a Cauchy sequence on O. Therefore there is a point z^∗ ∈ O such that z(t) → z^∗ as t → +∞. Continuity of S implies S(z∗) = S^∗ as claimed. At last, z^∗ ∈ O, for, if not, then z^∗ ∈ ∂O and S(z∗) > S(z0) ≥ S(z∗), which is a contradiction (the first inequality is strict). Passing to the limit t² → +∞ in inequality (38) yields

$$\left|z(t1) - z*\right| \le \theta e^{-\kappa t_1/2},$$

with $$\theta = \sqrt{4\alpha L\varphi_0}/\kappa$$ and any $t_1 \ge 0$ .

We now generalize this result to the case of finite-dimensional metriplectic systems. Under the same conditions as in Proposition 5, a simple generalization of (PL) for metriplectic vector fields reads: for any η ∈ I(U), inf{S(z): z ∈ Uη} =: S^η > −∞ and there is a constant κ^η > 0, depending on η, such that

$$\frac{1}{\kappa{\eta}} \nabla \mathcal{S}(z) \cdot K(z) \nabla \mathcal{S}(z) \ge \mathcal{S}(z) - \mathcal{S}{\eta}, \quad z \in \mathcal{U}_{\eta}. \tag{PL'}$$

Differently from (PL), which is a condition on S, inequality (PL′ ) involves both the entropy function and the bracket (S, S) = ∇S · K∇S, which gives the entropy decay rate. Similarly to (PL), if inequality (PL′ ) holds, then the metriplectic vector field has the following properties:

• If z(t) ∈ Uη, t ∈ [0, +∞), is an integral curve of X with initial condition z⁰ = z(0), then S z(t) → S^η as t → +∞, with exponential convergence. In fact, (PL′ ) implies

$$\left[\mathcal{S}(z(t)) - \mathcal{S}_{\eta}\right] \leq \left[\mathcal{S}(z0) - \mathcal{S}{\eta}\right] e^{-\kappa t},$$

which follows as in the case of standard gradient flows.

• If z^e ∈ U^η is an equilibrium point, then it is necessarily a global minimum of S|^U^η that is, S(ze) = Sη. This follows from the fact that at an equilibrium point necessarily ∇S(ze) · K(ze)∇S(ze) = 0, and (PL′ ) implies S(ze) = Sη.

We can now state the analog of Proposition 7 for the case of metriplectic vector fields. However, we consider only the dissipative part, that is, a metriplectic vector field of the form (7), without the symplectic part, and make an additional assumption that is sufficient to ensure that the orthogonal projection onto ker K(z) is smooth in z. Specifically, we assume that

there is a constant $$r > 0$$ , such that, for any $z \in \mathbb{Z}$ ,
zero is the only eigenvalue of $K(z)$ in the interval $[-r, r]$ . (39)

Hypothesis (39) implies that, for any z, only one eigenvalue of K(z), the zero eigenvalue, belongs in the disk of radius r > 0 centered at zero in the complex plane.

Proposition 8. Let X = −K∇S be a metriplectic vector field of the form (7) on a domain Z ⊆ R ⁿ, with K satisfying (39) and let I = (I 1 , . . . , I k ) ∈ C∞(Z, R k ), 1 ≤ k < n, be such that K∇I^α = 0. Assume that rank ∇I = k in a subdomain U satisfying (34) and let O be the open set given by Lemma 4. If (PL′ ) holds in U, then for any integral curve z(t) of X with initial condition z(0) = z⁰ ∈ O there is a point z^η ∈ O and a constant θ^η depending on z0, such that η = I(zη) = I(z0), S(zη) = S^η and

$$\left|z(t)-z{\eta}\right| \leq \theta{\eta}e^{-\kappa{\eta}t/2}\,, \quad \left[\mathcal{S}\big(z(t)\big)-\mathcal{S}{\eta}\right] \leq \left[\mathcal{S}\big(z{0}\big)-\mathcal{S}{\eta}\right]e^{-\kappa_{\eta}t}\,, \quad for \ t \geq 0\,.$$

The proof of proposition 7 can be adapted to this case. The key point is replacing ∇S with (I − π0)∇S, where π⁰ is the orthogonal projector onto the ker K.

Proof of Proposition 8. Per hypothesis, both K and S are smooth on Z and thus X is Lipschitz continuous on any bounded subset and in particular on U. Then Lemma 4 gives an open subset of O ⊆ U such that for any point z⁰ ∈ O there is an integral curve z(t) ∈ O of X through the point z⁰ = z(0), defined for all t ≥ 0. Along the orbit I α z(t) = I ^α(z0) = η ^α, α ∈ {1, . . . , k}, and thus z(t) ∈ U^η ∩ O.

At a point z^η ∈ U^η where S(zη) = S^η the function S|^U^η attains its minimum and thus S must satisfy (35). It follows from the assumption K∇I^α = 0 that X(zη) = 0. Therefore, if there is t¯≥ 0 such that S z(t¯) = Sη, then X z(t¯) = 0 and z(t) = z(t¯) = z⁰ for all t ≥ 0. In this case the statement of the proposition is trivially true.

Let us now consider the non-trivial case S z(t) > Sη, t ≥ 0. We shall follow the same two steps as in Polyak original proof, with the necessary changes to account for the degeneracy of the metriplectic flow. This requires a preliminary step in which we establish the needed properties of the matrix K.

Step 0. For any point z ∈ Z, let Ki(z) > 0 be the i-th non-zero eigenvalues of K(z) and πi(z) the orthogonal projector on ker K(z)−Ki(z)I), i.e. on the eigenspace of K(z) corresponding to the eigenvalue Ki(z). Then

$$K(z) = \sum_{i} K_i(z)\pi_i(z).$$

Since K(z) is symmetric any eigenvalue (including zero) is semisimple [88, Appendix 3.I]. We assumed that the closed disk of radius r > 0 centered in 0 ∈ C contains only the zero eigenvalue of K(z) for all z ∈ Z, hence Ki(z) > r. It follows that [88, Theorem 3.I.1] the orthogonal projector π0(z) onto ker K(z) is a smooth function of z ∈ Z. In addition we have that, for any vector Z ∈ R n,

$$Z \cdot K(z)Z = \sum{i} K{i}(z)Z \cdot \pi{i}(z)Z$$ $$\geq r \sum{i} Z \cdot \pi{i}(z)Z = rZ \cdot (I - \pi{0})Z.$$

If Z = ∇S(z), we deduce

$$\nabla S(z) \cdot K(z) \nabla S(z) \ge r \left| \left( I - \pi_0(z) \right) \nabla S(z) \right|^2. \tag{40}$$

We also have X = −K∇S = −K(I − π0)∇S, and |X(z)| ≤ ∥K(z)∥^F I − π0(z) ∇S(z) where ∥K(z)∥^F is the Frobenius norm, which is a continuous function of z. For any z in the compact set U, ∥K(z)∥^F ≤ R = max{∥K(z ′ )∥^F : z ′ ∈ U}, so that

$$|X(z)| \le R |(I - \pi_0(z)) \nabla S(z)|. \tag{41}$$

Since π⁰ is smooth, Y (z) = I − π0(z) ∇S(z) is smooth. This is the component of ∇S orthogonal to ker K.

From (41) we can also deduce that Y z(t) ̸= 0 for, if not, then (PL′ ) implies S z(t) = Sη, which is the trivial case.

We shall show that Polyak's argument can be repeated with Y instead of ∇S.

Step 1. Upon using (40), given 0 ≤ t¹ < t2,

$$S(z(t_1)) - S(z(t2)) \ge r \int{t_1}^{t_2} |Y(z(t))|^2 dt.$$

On the other end, (41) yields

$$\frac{d}{dt} |z(t) - z(t_1)| \le |X(z(t))| \le R|Y(z(t))|$$

Since Y ∈ C^∞(Z, R ⁿ), it is in particular Lipschitz continuous on O. Let L > 0 be the Lipschitz constant of Y . The same argument of step 1 in the proof of Proposition 7 can be repeated leading to

$$|Y(z(t))|^2 \le \frac{\alpha RL}{r} [S(z(t)) - S],$$ (42)

where α is the same constant defined in the proof of Proposition 7.

Step 2. We have already shown that

$$\varphi{\eta}(t) \le \varphi{0,\eta} e^{-\kappa_{\eta} t},$$

where φη(t) = S z(t) − Sη, and φ0,η = φη(0). For any 0 ≤ t¹ < t2, inequality (42) gives

$$\begin{split} \left|z(t_1) - z(t2)\right| &\leq R \int{t_1}^{t2} \left|Y\big(z(t)\big)\right| dt \ &\leq R \sqrt{\frac{\alpha R L}{r}} \int{t_1}^{t2} \sqrt{\varphi{\eta}(t)} dt \ &\leq \frac{1}{\kappa{\eta}} \sqrt{\frac{4\alpha R^3 L \varphi{0,\eta}}{r}} \left[e^{-\kappa_{\eta} t1/2} - e^{-\kappa{\eta} t_1/2}\right]. \end{split}$$

It follows that z(t) has a limit z^η ∈ O for t → +∞ and the limit satisfies S(zη) = limt→+^∞ S z(t) = Sη, hence z^η ∈ O. In addition, I(zη) = limt→+^∞ I z(t) = I(z0), hence z^η ∈ U^η ∩ O. At last, passing to the limit t² → +∞ yields

$$|z(t1) - z{\eta}| \le \theta{\eta} e^{-\kappa{\eta} t_1/2},$$

with constant

$$\theta{\eta} = \frac{1}{\kappa{\eta}} \sqrt{\frac{4\alpha R^3 L \varphi_{0,\eta}}{r}}.$$

This is the claimed inequality.

This exponential convergence result becomes less useful when the constant θ^η is large, which can happen when either r or κ^η are small. In the former case, an eigenvalue of K becomes small at least in some region of the domain. In the latter case, the metric bracket is small even where entropy is far from the minimum.

This result follows from a minimal modification of Polyak's argument for standard gradient flows. One should note that here ∇I^α are assumed to be in the kernel of K, and this assumption is consistent with the requirement (7b) for H.

Unlike the results of Section 3.1, convergence results based on Polyak inequalities do not require the uniqueness of the minimum entropy state. On the other hand, the precise point on the set of minima at which each orbit converges depends on the specific orbit. We illustrate inequality (PL′ ) with a few examples.

Example 1. Let the phase space be Z = R ⁿ, n ∈ N and n ≥ 2, with coordinates z = (z i )i , s = (si)ⁱ and h = (hi) ∈ R ⁿ be a covariant vectors, K = (Kij ) be a symmetric, positive-semidefinite, contravariant tensor on R ⁿ such that Kijh^j = 0, and let σ = (σij )ij be a symmetric positive definite, covariant tensor. We define the dissipative part of the metriplectic vector field X = −K∇S, with Hamiltonian H(z) = hiz i , and entropy S(z) = siz ⁱ + σij z i z j . We further assume that the null space of K coincides with the line spanned by h, i.e., Kω = 0 implies ω = λh for some λ ∈ R. Then the metric bracket defined by K is "minimally degenerate", in the sense defined above.

We claim that this metriplectic system satisfies condition (PL′ ), and thus Proposition 8 applies.

In order to show this, let us first observe that the change variables z 7→ z˜ = σ 1/2 z transform the system into an analogous one with σij replaced by δij and with h and s replaced by σ ⁻1/2h and σ −1/2 s, respectively. Hence it is enough to discuss the case σij = δij . We can also assume |h| ² = 1, because the normalization of h only changes the value of the Hamiltonian but not its isosurfaces. Then, for any η ∈ R, U^η = {z ∈ R ⁿ : H(z) = η} is the plane given by hiz ⁱ = η. A point z ∈ U^η can be written as z = ηh + z⊥, with z^⊥ = z − (h · z)h. Given η ∈ R, we can use Lagrange multipliers to compute the constrained entropy minima Sη: we search for (λ, z) such that (with ∂ⁱ = ∂/∂zⁱ ) ∂iS(z) = λ∂iH(z) with H(z) = η, which is equivalent to

$$\begin{cases} z = \lambda h - s, \ h \cdot z = \eta. \end{cases}$$

The solution (λη, zη) is readily found,

$$\lambda{\eta} = \eta + h \cdot s, \quad z{\eta} = \lambda{\eta} h - s = \eta h - s{\perp},$$

where s^⊥ = s − (h · s)h. Therefore, if z ∈ Uη,

$$S(z) = s \cdot (\eta h + z{\perp}) + \frac{1}{2} |\eta h + z{\perp}|^2 = S(z{\eta}) + \frac{1}{2} |z{\perp} + s_{\perp}|^2.$$

On the other hand, since Kh = 0, for any z ∈ Uη,

$$\nabla \mathcal{S}(z) \cdot K \nabla \mathcal{S}(z) = (z+s) \cdot K(z+s) = (z{\perp} + s{\perp}) \cdot K(z{\perp} + s{\perp}) \ge K1 |z{\perp} + s_{\perp}|^2,$$

where K¹ > 0 is the smallest eigenvalue of K restricted to (ker K) ^⊥, the orthogonal of its kernel. We deduce

$$\nabla S(z) \cdot K \nabla S(z) \ge 2K1 [S(z) - S{\eta}], \quad z \in \mathcal{U}_{\eta},$$

where S^η = S(zη) is the constrained minimum of the entropy. This is condition (PL′ ) with κ^η = 2K1.

We remark that, at least in this case, the condition on the kernel of K being "minimal" is crucial for the modified PL condition. In fact, if there is a vector h ′ ∈ R ⁿ, orthogonal to h, and such that Kh′ = 0, then (z^⊥ +s⊥)·K(z^⊥ +s⊥) = 0 for any non-zero z^⊥ +s^⊥ ∝ h

This example is simple enough that an analytical solution of the integral curves of X can be obtained. In fact, the equation for the new variable y = s + z amounts to the linear system dy/dt = −Ky with initial condition y⁰ = s + z0. If H(z0) = h · z⁰ = η, we must have h · y⁰ = η + h · s. Upon representing y on the basis of the unit eigenvectors {ei} n−1 ⁱ=0 of K, with e⁰ = h being the eigenvector that corresponds to the zero eigenvalue, we obtain

$$y(t) = (h \cdot y0)h + \sum{i \ge 1} e^{-t\lambda_i(K)} c_i e_i,$$

where λi(K) > 0 are the positive eigenvalues of K and cⁱ = eⁱ · y0. We deduce

$$|z(t) - z_{\eta}| \le |z0 - z{\eta}|e^{-tK_1}$$

with K¹ = mini≥1{λi(K)}.

Hence, in this case we have exponential convergence to the equilibrium point, with convergence rate being half of the constant in (PL′ ).

Example 2. With the same metric bracket and Hamiltonian as in Example 1, let us consider the entropy function

$$S(z) = \frac{|z|^2}{1+|z|^2}, \quad z \in \mathcal{Z} = \mathbb{R}^n.$$

As before, this entropy is rotationally symmetric, with a global minimum at z = 0, but it is not a convex function.

Since H is the same as in Example 1, z ∈ U^η if and only if z = ηh + z^⊥ and we compute

$$\nabla \mathcal{S}(z) \cdot K \nabla \mathcal{S}(z) = 4 \frac{z{\perp} \cdot K z{\perp}}{(1+|z|^2)^4} \ge 4K1 \frac{|z{\perp}|^2}{(1+|z|^2)^4}.$$

where K¹ is defined in Example 1. We can use Lagrange multipliers in order to compute minima of the entropy constrained to U^η with the result that there is a unique minimum at z^η = ηh and

$$S{\eta} = S(z{\eta}) = \frac{\eta^2}{1 + \eta^2}.$$

Then, we compute

$$S(z) - S{\eta} = \frac{|z{\perp}|^2}{(1 + \eta^2)(1 + |z|^2)},$$

from which we deduce

$$\frac{|z{\perp}|^2}{(1+|z|^2)^4} = \frac{1+\eta^2}{(1+|z|^2)^3} [\mathcal{S}(z) - \mathcal{S}{\eta}], \quad z \in \mathcal{U}_{\eta},$$

hence, for any R > 0, on the ball |z| < R, we have

$$\nabla S(z) \cdot K \nabla S(z) \ge \kappa{\eta} [S(z) - S{\eta}], \quad z \in \mathcal{U}{\eta} \cap B{R}(0),$$

with constant κ^η = 4K1(1 + η )/(1 + R² ) . Therefore the modified PL condition is satisfied on balls of arbitrary large radius R, even though the entropy is not convex, but the constant as a function of the radius R is not uniformly bounded away from zero.

Example 3. As a last example, we consider a strongly nonlinear case with a bracket built from an orthogonal projection onto the hyper-plane perpendicular to the gradient of the Hamiltonian. This particular metric structure will play a key role in the following, even in the infinite-dimensional cases in fluid and plasma dynamics.

Given s ∈ R ⁿ, on the open half-space Z = {z ∈ R ⁿ : z · s < 0}, let us consider the field X(z) = −K(z)∇S(z), with

$$K(z) \coloneqq |z|^2 I - z \otimes z$$ , $\mathcal{H}(z) \coloneqq \frac{1}{2} |z|^2$ , and $\mathcal{S}(z) \coloneqq s \cdot z$ .

Since K(z) is proportional to the projector onto the subspace normal to ∇H(z), we have K(z)∇H(z) = 0 and K(z) is a symmetric positive semidefinite tensor; hence, X is metriplectic with a trivial symplectic part. We also stress that K(z) is minimally degenerate since K(z) is strictly positive definite on the subspace normal to ∇H(z).

The constant-energy surfaces are spheres, for any η > 0,

$$z \in \mathcal{U}_{\eta} \iff z = \sqrt{2\eta}\zeta, \quad \zeta \in S^{n-1}, \quad \zeta \cdot s < 0,$$

where points ζ on the (n − 1)-dimensional sphere S ⁿ−¹ are identified with unit vectors in R ⁿ. The entropy restricted to U^η amounts to S|^U^η (ζ) = ^√ 2ηs · ζ and the minimum S^η = − p 2ηs² is attained at ζ = −s/|s|. The same result is of course obtained by means of Lagrange multipliers that lead to the system

$$\begin{cases} s = \lambda z, \ \frac{1}{2}|z|^2 = \eta, \end{cases} \text{ and } z \cdot s < 0.$$

Then we compute, for z = 2ηζ ∈ Uη,

$$\begin{split} \nabla \mathcal{S}(z) \cdot K(z) \nabla \mathcal{S}(z) &= |z|^2 |s|^2 - (z \cdot s)^2 = \left( \sqrt{2\eta} |s| - \sqrt{2\eta} s \cdot \zeta \right) \left( \sqrt{2\eta} |s| + \sqrt{2\eta} s \cdot \zeta \right) \ &\geq \sqrt{2\eta} |s| \left[ \mathcal{S}(z) - \mathcal{S}_{\eta} \right], \end{split}$$

which is inequality (PL′ ).

It should be noted that, if one drops the condition z · s < 0, so that Z = R ⁿ, then the metric system cannot satisfy (PL′ ) since S|^U^η has two critical points: a minimum z ^η with z η · s < 0 and a maximum z + ^η with z + η · s > 0. The metric bracket (S, S) = ∇S · K∇S vanishes at both points, but S(z + η ) − S^η > 0.

3.3. Infinite-dimensional systems: tentative generalization of the PL inequality

A version of the Lyapunov stability theorem, valid for the case of infinite-dimensional systems, is available under suitable coercivity assumptions on the Lyapunov function [53]. Such assumptions are needed to compensate for the lack of compactness. For instance, the closed unit ball is not compact in a Banach space. Compactness of closed and bounded sets in R ⁿ is used repeatedly in the classical proofs in finite dimensions (cf. Sections 3 and 3.2). Analogously, the (PL) condition can be extended to infinitedimensional systems. A more difficult point is the existence of a global-in-time solution to the equation defining the dynamical system, under reasonable hypothesis [79]. In the infinite-dimensional setting, this means proving the existence of a global solution for highly nonlinear partial differential equations, which is often difficult and requires special treatment for each individual case. Nonetheless, under the assumption that a global-in-time solution exists, one can think of extending Propositions 5 and 8 to infinite dimensions, but we leave the details for future work. Here we merely state the infinitedimensional version of (PL′ ).

Consider a metriplectic system on a Banach space V as introduced in Section 2.1. We assume that this system has a finite (for simplicity) family of constants of motion I ∈ C^∞(V, R k ), that satisfies the hypotheses of the submersion theorem [53, Theorem 3.5.4] in an open set U ⊆ V . In particular, the operator DI(u) must be surjective with split kernel for any u ∈ U. Then U^η = {u ∈ V : I(u) = η} are closed submanifolds of U for any η ∈ I(U), as in the finite-dimensional case. Since there is at least one invariant, namely the Hamiltonian, we have k ≥ 1. Then the equivalent of condition (PL′ ) reads: S^η := inf{S : z ∈ Uη} > −∞ and there exists a constant κ^η > 0 depending on η, such that

$$\frac{1}{\kappa{\eta}}(\mathcal{S}, \mathcal{S}) \ge \mathcal{S} - \mathcal{S}{\eta}, \quad \text{on } \mathcal{U}_{\eta}.$$ (PL")

If inequality (PL′′) is fulfilled, the exponential convergence of the entropy follows as in the finite-dimensional case. Also, (S, S)(u) = 0 on U^η only if u is a global minimum of S restricted to Uη, i.e. S(u) = Sη.

A necessary condition for (PL′′) can be stated for the special class of metric brackets that satisfy

$$(\mathcal{F}, \mathcal{F})(u) = 0 \iff D\mathcal{F}(u) = \sum{\alpha} \lambda{\alpha} D\mathcal{I}^{\alpha}(u),$$ (43)

for some constant λ^α ∈ R. This generalizes condition (L3′′) for specifically degenerate brackets. Then, condition (PL′′) is satisfied only if critical points of S|^U^η are global minima. In fact, if u ∈ U^η is a critical point of S|^U^η , the theory of Lagrange multiplier [53, Theorem 3.5.27] gives DS(u) = P α λαDI ^α(u), hence (S, S)(u) = 0. But if S(u) > Sη, inequality (PL′′) is violated.

Beyond these preliminary considerations, the mathematical analysis of (PL′′) exceeds the scope of this paper. We conclude with an example of a infinite-dimensional metric bracket that is specifically degenerate and satisfies (PL′′). We shall give physically relevant examples in Section 4.

Example 4. In this example (cf. [23]) we proceed formally. On the space V of smooth functions from T → R, where the torus T := R/2πZ is identified with the interval Ω = [0, 2π] with periodic boundary conditions, we consider the metric bracket given by

$$\left(\mathcal{F},\mathcal{G}\right) = \int_{0}^{2\pi} \left(\frac{\delta \mathcal{F}(u)}{\delta u}\right)' \left(\frac{\delta \mathcal{G}(u)}{\delta u}\right)' dx,\tag{44}$$

where u ′ (x) = du(x)/dx denotes the derivative of u : T → R, and the functional derivatives are computed with respect to the standard L ² product (cf. Section 2.1). We assume that the functions F(u) and G(u) are regular enough for their functional derivative to exist and be sufficiently smooth. The Hamiltonian and entropy functions are given by

$$\mathcal{H}(u) = \int_0^{2\pi} u(x) dx \quad \text{and} \quad \mathcal{S}(u) = \frac{1}{2} \int0^{2\pi} |u(x)|^2 dx = \frac{1}{2} ||u||{L^2(\Omega)}^2.$$

Condition (7b) is satisfied since δH(u)/δu = 1. After integration by parts, the strong form of (7a) amounts to the heat equation

$$\begin{cases} \partial_t u = \partial_x^2 u, & (t, x) \in [0, +\infty) \times [0, 2\pi], \ u(t, 0) = u(t, 2\pi), & t \in [0, +\infty), \ u(0, x) = u_0(x), & x \in [0, 2\pi]. \end{cases}$$ (45)

First we show that (44) is a minimally degenerate bracket, i.e., the null space is spanned by δH(u)/δu. In fact, (F, F) = 0 implies (δF(u)/δu) ′ = 0 and thus δF(u)/δu = λ where λ ∈ R is constant. Hence, δF(u)/δu = λδH(u)/δu, which proves property (43), with the Hamiltonian being the only invariant.

The manifolds of constant Hamiltonian,

$$\mathcal{U}_{\eta} = { u : \mathcal{H}(u) = \eta \in \mathbb{R} },\$$

consist of functions with the same average over [0, 2π]. They are affine spaces, rather than generic manifolds. The critical points of entropy restricted to the constant-Hamiltonian spaces are determined by

$$u(x) = \lambda, \quad \int_0^{2\pi} u(x)dx = \eta.$$

Therefore, for any η ∈ R there is only one critical point, that is, the constant function

$$u_{\eta}(x) = \eta/(2\pi).$$

The entropy of u^η is S^η = S(uη) = η ²/(4π). The Fourier series representation,

$$u(x) = \sum_{n \in \mathbb{Z}} u_n e^{inx}, \quad u_n \in \mathbb{C},$$

yields that H(u) = η if and only if u⁰ = η/(2π), hence

$$S(u) - S{\eta} = \pi \sum{n \neq 0} |un|^2 \ge 0, \quad u \in \mathcal{U}{\eta},$$

with equality only if u = uη. This shows that u^η is a global minimum of S restricted to Uη. Upon using again the Fourier series representation, one finds

$$(S,S)(u) = ||u'||{L^2(\Omega)}^2 = 2\pi \sum{n\neq 0} n^2 |un|^2 \ge 2\pi \sum{n\neq 0} |u_n|^2$$

and

$$(S, S)(u) \ge 2[S(u) - S{\eta}], \quad u \in \mathcal{U}{\eta},$$

which is condition (PL′′) with κ^η = 2.

In fact, the solution of (45) can be readily written in terms of a Fourier series as

$$u(t,x) = \sum{n \in \mathbb{Z}} e^{inx - n^2 t} u{0,n},$$

where u0,n ∈ C are the Fourier coefficients of the initial condition. With η = H(u0) = 2πu0,0, we deduce

$$||u(t) - u{\eta}||{L^{2}(\Omega)} \le e^{-t} ||u(0) - u{\eta}||{L^{2}(\Omega)},$$

which shows exponential relaxation toward the entropy minimum on Uη, with the energy η being determined by the initial condition.

In conclusion, this metric system satisfies the generalized Polyak– Lojasiewicz condition (PL′′) and all orbits completely relax exponentially to a solution of (1) with exponential convergence rate given by κη/2, where κ^η is the constant in (PL′′). We observe that the convergence rate is the same as the one predicted in Proposition 8 for finite-dimensional systems.

4. Two examples: metric double brackets and projectors

In this section, we discuss two special cases of metriplectic systems of the form (7), which we call metric double bracket and projector bracket systems. We shall see that for the metric double brackets treated in Section 4.1, the dissipation mechanism does not completely relax the state of the system (in the sense of Section 3), while for the projector brackets of Section 4.2 complete relaxation is achieved. We shall discuss and compare the properties of these two metric brackets on the basis of the insights gained in Section 3.

We select two benchmark equilibrium problems, and we attempt to construct a relaxation method to solve them by using the two considered metric brackets. The benchmark problems are: the reduced Euler equations (cf. Section 2.2.1), but with Dirichlet boundary conditions replaced by periodic boundary conditions, and an analytically solvable model derived from the reduced Euler equations. In both cases periodic boundary conditions give rise to an additional invariant other than the Hamiltonian, and this allows us to examine cases in which the concept of "specifically degenerate" bracket is not equivalent to "minimally degenerate". We shall return to the original problem with Dirichlet boundary conditions later.

Therefore in both cases, the domain is T 2 := (R/2πZ) ² with coordinates x = (x1, x2), and the phase space V is the space of smooth functions v : T ² → R. As usual, T 2 is identified with the square Ω = [0, 2π] ² with periodic boundary conditions. On such a periodic domain, the scalar vorticity ω = −∆ϕ (cf. Section 2.2.1 for the definitions) must have zero average, i.e.,

$$\omega{\Omega} \coloneqq \frac{1}{4\pi^2} \int{\Omega} \omega(x) dx = 0.$$

We use systematically the subscript Ω to denote the average over the domain Ω. We choose the whole space V as the phase space and define the vorticity by

$$\omega = u - u_{\Omega},$$

but we impose the additional constraint

$$\mathcal{M}(u) = \int{\Omega} u(x)dx = 4\pi^2 u{\Omega} = \mathcal{M}_0 \in \mathbb{R},$$

as well as energy conservation H(u) = H⁰ ∈ R. Given u ∈ V , and thus ω, the Poisson equation determines the stream function ϕ modulo a constant, which we set to zero; hence (10) is replaced by

$$-\Delta \phi = u - u{\Omega}, \quad \phi{\Omega} = 0.$$ (46)

(Equivalently, we could have chosen the phase space to be the subspace of functions satisfying u^Ω = 0 and u = ω.)

Both the considered benchmark problems can be formulated as variational problems: given a Hamiltonian function H, and a regular value η = (M0, H0) for the two invariants I = (I 1 , I ) := (M, H), find

$$\min{S(u): I(u) = \eta},$$ (47a)

with entropy

$$S(u) = \frac{1}{2} \int{\Omega} \omega^2 dx = \frac{1}{2} ||u - u{\Omega}||_{L^2(\Omega)}^2.$$ (47b)

The two test cases differ by the choice of the Hamiltonian. In summary,

• Analytical test case: Given a function h: T ² → R, let

$$\mathcal{H}(u) = \int{\Omega} h \,\omega \, dx = (h - h{\Omega}, u)_{L^{2}(\Omega)}. \tag{48}$$

The solutions of (47) with (48) are equilibria of the linear advection equation

$$\partial_t u + [h, u] = 0,$$

with [f, g] = ∂1f ∂2g − ∂1g∂2f.

• Reduced Euler test case: We again use S as in (47b), but now

$$\mathcal{H}(u) = \frac{1}{2} \int{\Omega} |\nabla \phi|^2 dx = (\phi, u){L^2(\Omega)},$$ (49)

i.e., we assume (14) with s(y) = y ²/2, and we assume ϕ depends on u via (46). Solutions of (47) with (49) are equilibria of the reduced Euler equations on the flat torus T 2

In both cases, problem (47) can be solved analytically. In order to compute the solutions, we first find the set of critical points of entropy restricted to U^η = {u: I(u) = η}, i.e.,

$$\mathfrak{C}{\eta} := {u \colon D\mathcal{S}(u) = \sum{\alpha} \lambda_{\alpha} D\mathcal{I}^{\alpha}(u), \, \mathcal{I}(u) = \eta}.$$ (50)

Then we find the minimum of S on Cη. We summarize here the results, which will be used to assess the metriplectic relaxation methods.

Solution for the analytical test case: For the case of Eq. (48), the set C^η is given by

$$u - u_{\Omega} = \lambda_1 + \lambda2 (h - h{\Omega}), \quad \mathcal{I}(u) = \eta.$$ (51)

Then λ¹ = 0 and u^Ω = M0/(4π ). Upon multiplying by h − h^Ω and integrating over Ω, we deduce

$$\mathcal{H}_0 = \lambda2 ||h - h{\Omega}||_{L^2(\Omega)}^2$$ ,

from which we can compute λ2. Hence, the set C^η contains the following single point:

$$u_{\eta} = \frac{\mathcal{M}_0}{4\pi^2} + \frac{\mathcal{H}0}{|h - h{\Omega}|{L^2(\Omega)}^2} (h - h{\Omega})$$ (52)

and the value of entropy on this unique critical point is

$$S_{\eta} = \min{S(u) : \mathcal{I}(u) = \eta} = \frac{\mathcal{H}0^2/2}{|h - h{\Omega}|_{L^2(\Omega)}^2}.$$ (53)

One can check that this is the minimum of the entropy on Uη.

Solution for the reduced Euler equations: For the case of (49), elements of C^η satisfy

$$\begin{cases} u - u_{\Omega} = \lambda_1 + \lambda_2 \phi, & \mathcal{I}(u) = \eta, \ \text{with } \phi \text{ solution of (46).} \end{cases}$$ (54)

Since ϕ^Ω = 0, we have λ¹ = 0 and u^Ω = M0/(4π 2 ), as before. Then (ϕ, λ2) must be a solution of the eigenvalue problem

$$-\Delta \phi = \lambda2 \phi, \quad \phi{\Omega} = 0,$$

which is readily solved in terms of Fourier series. We find that the set C^η consists of vorticity fields u−u^Ω = ω = λ2ϕ, with ϕ being an eigenfunction of −∆ corresponding to the eigenvalue λ² > 0, and with norm ∥ϕ∥ ^L2(Ω) = 2H0/λ2. Then the entropy evaluated on the constrained critical points amounts to λ2H0. It follows that the entropy minimum on U^η corresponds to the lowest non-trivial eigenvalue, which is λ² = 1. The corresponding stream function must be of the form

$$\phi(x) = a_1 \cos(x_1 + \theta_1) + a_2 \cos(x_2 + \theta_2),$$

with arbitrary phase shifts θ1, θ2, and with coefficient a1, a² determined by the condition ∥ϕ∥ 2 ^L2(Ω) = 2H0. Since λ² = 1, this amounts to a 2 ¹ + a ² = H0/π² . Thus,

$$\omega(x) = \phi(x) = \frac{\sqrt{\mathcal{H}_0}}{\pi} \left[ \cos \theta_0 \cos(x_1 + \theta_1) + \sin \theta_0 \cos(x_2 + \theta_2) \right], \tag{55}$$

with arbitrary phases θ0, θ1, and θ² in [0, 2π).

From the analytical solution, (55), we deduce that the entropy minimum constrained to I(u) = η is not attained at an isolated point, but on a family of points parameterized by three phases. The constrained minimum value of the entropy is given by

$$S_{\eta} = \min{S(u): \mathcal{I}(u) = \eta} = \mathcal{H}_0,$$ (56)

since λ² = 1.

4.1. Metric double brackets

The first example is given by the metric double bracket [42, 50, 51] (not to be confused with the double bracket of Flierl and Morrison [34] discussed in the introduction). We recall the general definition first, but quickly restrict the discussion to the examples.

In general, metric double brackets originate from a Lie algebra. If g is a Lie algebra with Lie brackets [·, ·], let the vector space V be its dual, V = g ∗ . The functional derivative of a function f ∈ C∞(g ∗ ) is computed with respect to the duality pairing between g ∗ and g, that is, δf(u)/δu ∈ g is the unique element of g such that Df(u)v = ⟨δf(u)/δu, v⟩ for all v ∈ V = g ∗ . Under these conditions, it is well-known [53] that the Lie bracket in g induces two Poisson brackets in C^∞(g ∗ ), namely, {f, g}^± = ±⟨u, [ δf(u) δu , δg(u) δu ]⟩, so that both (g ∗ , {·, ·}±) are Poisson manifolds. However, if in addition g is equipped with a positive definite bilinear form γ : g × g → R, on the space of smooth functions C^∞(g ∗ ), we can also define the symmetric bracket

$$(f,g) = \gamma \Big( \Big[ \frac{\delta f(u)}{\delta u}, \frac{\delta h(u)}{\delta u} \Big], \Big[ \frac{\delta g(u)}{\delta u}, \frac{\delta h(u)}{\delta u} \Big] \Big),$$

for any fixed function h ∈ C^∞(g ∗ ). One can readily check that this is a metric bracket preserving the Hamiltonian function h.

Formally at least, this construction can be extended to infinite-dimensional systems. Let V and W be Banach spaces, with a nondegenerate duality pairing ⟨·, ·⟩^V ^×^W : V × W → R, and let W be equipped with (i) a symmetric positive definite bilinear form γ : W ×W → R and (ii) a bilinear antisymmetric operation [·, ·] : W ×W → W. (For the purpose of defining the metric bracket we do not need to require [·, ·] to be a Lie bracket, that is, we can relax the Jacobi identity.) Then, given a fixed Hamiltonian H ∈ C^∞(V ), we can construct the bilinear form (·, ·) : C^∞(V ) × C^∞(V ) → C^∞(V ) given by [50, Eq. (2.9)]

$$(\mathcal{F}, \mathcal{G}) := \gamma \Big( \left[ \frac{\delta \mathcal{F}}{\delta u}, \frac{\delta \mathcal{H}}{\delta u} \right], \left[ \frac{\delta \mathcal{G}}{\delta u}, \frac{\delta \mathcal{H}}{\delta u} \right] \Big), \tag{57}$$

where the functional derivative are evaluated with respect to the duality pairing between V and W, and thus are elements of W. We remark that when W = V ′ is the (topological) dual of V , that is, the space of continuous linear functionals on V , δF(u)/δu exists and it is equal to DF(u) for any F ∈ C (V ) and for all u ∈ V . In general, however, δF(u)/δu does not always exists for all F.

As an example, let V = W be the space of smooth functions from T ^d → R, identified with functions over Ω = [0, 2π] ^d with periodic boundary conditions. If [u, v]^J = ∇u·J∇v is a Poisson brackets on R d (not necessarily canonical), H(u) is a given Hamiltonian function, and γ is given by the standard product in L 2 (Ω), then (57) reduces to

$$(\mathcal{F}, \mathcal{G}) := \int{\Omega} \left[ \frac{\delta \mathcal{F}}{\delta u}, \frac{\delta \mathcal{H}}{\delta u} \right]{J} \left[ \frac{\delta \mathcal{G}}{\delta u}, \frac{\delta \mathcal{H}}{\delta u} \right]_{J} dx. \tag{58}$$

In the following, let d = 2, x = (x1, x2) ∈ Ω = [0, 2π] ² ⊂ R ² with periodic boundary conditions, and for J, we choose the canonical Poisson tensor in R 2 , so that [·, ·]^J = [·, ·] is the canonical bracket defined after Eq. (11).

We construct a relaxation method based on this bracket for the solution of the two test problems introduced at the beginning of this section. In both cases we consider a field u ∈ V evolving from an initial condition u⁰ according to Eqs. (7) and with bracket given by (58).

We start from problem (47) with (48) for the linear advection equation. With those choices of bracket, Hamiltonian and entropy, Eq. (7a) amounts to

$$\int{\Omega} \frac{\delta \mathcal{F}}{\delta u} \frac{\partial u}{\partial t} dx = -\int{\Omega} \left[ \frac{\delta \mathcal{F}}{\delta u}, h \right] \left[ u, h \right] dx,$$

for all F. This can be viewed as the weak form of the evolution equation, with δF/δu being the test function. After integration by parts, one obtains the evolution equation for u in strong form, namely,

$$\partial_t u = [h, [h, u]].$$

We observe that [u, h] = div(Xhu) = X^h · ∇u, where X^h = ^t (∂2h, −∂1h) is the Hamiltonian vector field generated by h with canonical Poisson bracket in R 2 . Hence

$$\partial_t u = \operatorname{div}(X_h \otimes X_h \nabla u),$$ (59)

which shows that this particular combination of metric bracket and entropy describes anisotropic diffusion, parallel to the field lines of Xh, or equivalently along the contours of the function h.

The Cauchy problem associated to (59) with periodic boundary conditions can be solved analytically. The limit for t → +∞ of the solution exists and is equal to the average of the initial condition on the contours of h. Therefore, in general the limit of the solution retains some information of the initial condition, while the completely relaxed solutions (52) only depends on the energy H⁰ = H(u0) of the initial condition u0. This implies that for generic initial conditions, the limit of the solution of a metric dynamical system is in general not a solution of the variational principle (47).

Figure 1 shows the result of the numerical solution of (59) with Hamiltonian

$$h(x) = \cos^2(x_1)\sin^2(x_2).$$ (60)

The contours of h form a periodic array of islands, cf. the black contours in the upper panels of Fig. 1. In each island, h takes the same values. The initial condition, represented as a color map, is an anisotropic Gaussian centered between two islands, i.e., u0(x) = uG(x) with

$$u_G(x) = \frac{1}{N} \exp\left[-\frac{(x1 - x{0,1})^2}{w_1^2} - \frac{(x2 - x{0,2})^2}{w_2^2}\right],\tag{61}$$

with center x⁰ = (π, π + 0.1), w¹ = 0.25, w² = 0.4, and N = 2πw1w2. The solution is obtained with a standard spectral method with Fourier basis, on a 256 × 256 uniform grid. The time integrator is the standard 4th order explicit Runge-Kutta method with time step ∆t = 10⁻⁴ . The parallel diffusion equation (59) tends to equalize the solution on the contours of h, but the dynamics at the boundary of the islands is very slow; the

Figure 1: Example of solution of Eq. (59) with Hamiltonian (60). Upper panels: initial condition and final state, compared to the contours of h (black circular curves). Lower panel: A visualization of the functional relation between h and u obtained by plotting the points (hij , uij ), with hij and uij being the values of h and u, at the node (i, j) of the computational grid. (For clarity, in the color maps we display the solution u only where u ≥ 10−⁴ .)

Hamiltonian vector field X^h at the boundary is zero and the solution remains constant on those contours (referred to as separatrices). The color map of the relaxed state is shown in Fig. 1 upper panel, while the lower panel represents the functional relation between h and the solution u by marking on the h-u plane a the point (hij , uij ) for each grid node xij . At the initial time, (black markers) there is no relation between the values of u and those of h, showing that the initial condition (61) is far from an equilibrium. As the solution evolves, all values of u sampled on the same contour of h tend to a common value, but the "condensation" of points on a line is slower for h small, that is, near the separatrices. Blue markers show to the average of the initial condition on each contour of h: one can see that the solution tends to the averages as predicted by the analytical solution. In the limit t → +∞ the relation between h and u is multi-valued with a countable set of branches, one for each island. In Fig. 1 one can distinguish the upper branch (larger values of u), corresponding to the island that contains the maximum of the initial condition. (In order to separate the two branches the center x⁰ of the Gaussian has been shifted up in the direction x2.) A second branch with lower values of u corresponds to the neighboring island. All the other islands do not overlap with the initial condition significantly and therefore appear as a line of points u = 0 for all h in Fig. 1 lower panel. In Fig. 1 the analytical solution (52) for the completely relaxed state is also shown (green crosses), and it is clearly different from the obtained equilibrium.

We now consider problem (47) with (49) for the Euler equations. For this choice of entropy and Hamiltonian, Eq. (7a) with (58) amounts to the anisotropic diffusion equation (59), but with h replaced by ϕ, which depends on u via (46). Hence the problem is nonlinear with a cubic nonlinearity, and in general no analytical solution is known (to the best of our knowledge).

Figure 2 shows an example of relaxation of an initially anisotropic vortex. The initial state is again given by (61), but now with x⁰ = (π, π), w¹ = 0.3, w² = 1.0, and N = 1, on [0, 2π] ² with a uniform mesh of 256 × 256 nodes. The time integrator is the standard 4th-order explicit Runge-Kutta method with time step ∆t = 10⁻³ . The solution relaxes to a symmetric vortex, which is an equilibrium of the Euler equations. During the evolution, the Hamiltonian is constant and the entropy is monotonically dissipated, consistent with (8). However, the entropy appears to converge to a value that is higher than its constrained minimum Sη, given in Eq. (56), and indicated by the thick horizontal line in Fig. 2. The fact that the final state is (a numerical approximation of) an equilibrium of the Euler equations can be deduced from the plot of the final state: the solution ω appears to be constant on the contours of the potential ϕ. A more quantitative indication is provided in Fig. 3, where the relation between ϕ and ω is represented. At the initial time t = 0, there is no functional relation between ω and ϕ: the values (ϕ, ω) on the computational grid do not belong to a curve. As the state relaxes, the scatter of points is reduced, and at the final time, one finds a clear functional relation.

For a completely relaxed state, we expect a linear relation between the potential and the vorticity, Eq. (55), and this is indicated by green crosses in Fig. 3. The obtained relationship is however very different, proving that the dynamical system reaches an equilibrium that does not satisfy the variational principle of minimum constrained entropy (1). From the numerical experiments in Fig. 3, one can see that the relaxed state is close to the average of the initial condition on the contours of the initial potential. The average of the initial condition gives the exact long-time limit of the solution in the case of the linear problem (59), and it is not expected to give an accurate prediction of

Figure 2: Metriplectic relaxation of a vortex toward an equilibrium of the reduced Euler equations, using (58) with (14) and s(y) = y ²/2. Top row: initial and finale state of the system; the color scheme represents the vorticity ω = u − uΩ; white lines represent the contours of the potential ϕ. Bottom row: relative error of the Hamiltonian and the value of entropy during the evolution. The thick horizontal line indicates the constrained entropy minimum S^η = H0, cf. Eq. (56).

Figure 3: Visualization of the relation between the potential ϕ and the vorticity ω, for the initial and final state of the calculation. The green crosses mark the linear relation for a minimum entropy state, Eq. (55). The data marked "averages" represent the average of the initial condition on the contours of the corresponding potential ϕ.

the relaxed state in general. Yet we observe that the solution converges to a state close to the average of the initial condition.

We conclude that the relaxation mechanism of (58) fails to capture the linear profile encoded in the choice of the entropy function. The relaxed state is an equilibrium of the reduced Euler equations, but corresponding to a profile that differs from the target one and that depends on the initial condition in a complicated way.

We can try to understand the behavior of this metric system in terms of the ideas put forward in Section 3.3. In the simple case of the linear advection equation, the Hamiltonian is given by (48). In this case (S, S) = 0 if and only if [u, h] = 0, and this is satisfied for u = f(h), for any sufficiently regular function f : R → R. In contrast the entropy minimum on U^η is attained at the unique point u^η given in Eq. (52). Hence inequality (PL′′) cannot be satisfied. Bracket (58) is not even specifically degenerate since for F(u) = (f(h), u)L2(Ω) one has (F, F) = 0 identically, and yet δF(u)/δu = f(h) is a linear combination of δI ^α(u)/δu only if f is affine. Therefore, bracket (58) with Hamiltonian (48) does not satisfy any of the (formally) sufficient conditions for relaxation to a minimum entropy state described in Section 3.3. In fact the analytical solution and the numerical experiment show that the relaxation method finds a point of the (rather large) set {u: (S, S) = 0}, instead of the unique entropy minimum (52).

For the case of reduced Euler equation, we have that (S, S) = 0 if u = f(ϕ), but now ϕ is related to u by the Poisson equation (46). A precise calculation of the set {u: (S, S) = 0} is far more complicated, but the numerical experiments show that this set is much larger than the set of points given in (55).

We conclude that brackets of the form (57) are in general not suited for the construction of relaxation methods for variational problems of the form (47). Yet they can be useful in another way as we shall see below for the nonlinear Beltrami field.

4.2. Projector-based metric bracket

We address now a construction of metric brackets based on L 2 -orthogonal projectors, which are patterned after that given for finite-dimensional systems in [3]. As before let Ω = [0, 2π] ^d and V be the space of functions on R d , 2π-periodic in each direction. Given a Hamiltonian function H, the L ² orthogonal projector onto the direction of δH(u)/δu is

$$\Pi{\mathcal{H}}(u)v := v - c(u, v) \frac{\delta \mathcal{H}(u)}{\delta u}, c(u, v) := \left| \frac{\delta \mathcal{H}(u)}{\delta u} \right|{L^{2}(\Omega)}^{-2} \left( \frac{\delta \mathcal{H}(u)}{\delta u}, v \right)_{L^{2}(\Omega)}.$$ (62)

With the projector, let us define

$$(\mathcal{F}, \mathcal{G}) := \left(\frac{\delta \mathcal{F}}{\delta u}, \Pi{\mathcal{H}} \frac{\delta \mathcal{G}}{\delta u}\right){L^{2}(\Omega)}.$$ (63)

We claim that the symmetric bi-linear form (63) satisfies (7b) and the Leibniz identity. Since Π^H is a projector ΠH(δH/δu) = 0, hence (F, H) = 0 for all F; in addition,

$$(\mathcal{F},\mathcal{F}) = \left(\frac{\delta \mathcal{F}}{\delta u}, \Pi{\mathcal{H}} \frac{\delta \mathcal{F}}{\delta u}\right){L^2(\Omega)} \geq 0,$$

since projectors are symmetric and nonnegative definite. The Leibniz identity is straightforward.

We utilize bracket (63) in (7) in order to obtain a relaxation method for the variational problems (47). For the case of the linear advection equation, energy (48) yields the evolution equation

$$\partialt u = -\left[u - u{\Omega} - \frac{\mathcal{H}(u)}{|h - h{\Omega}|{L^2(\Omega)}^2} (h - h_{\Omega})\right]$$

where u^Ω = M(u)/(4π 2 ). Since both M and H are constants of motion, the affine transformation u 7→ w = u−M0/(4π 2 )− H0/∥h−hΩ∥ 2 L2(Ω) (h−hΩ) with M⁰ = M(u0) and H⁰ = H(u0), u⁰ being the initial condition, transforms the equation into ∂tw = −w, which leads to the analytical solution

$$u(t,\cdot) = \left \frac{\mathcal{M}_0}{4\pi^2} + \frac{\mathcal{H}0}{|h - h\Omega|{L^2(\Omega)}^2} (h - h\Omega) \right + u_0 e^{-t}.$$

The term in square brackets is exactly the unique entropy minimum (52) on the manifold U^η with η = (M0, H0) being determined by the initial condition. Therefore for any initial condition u0, this metriplectic system relaxes completely to the accessible entropy minimum and with exponential convergence rate. This is the desired behavior. From the analytical solution one can see how the initial condition is quickly "forgotten", leaving only the fully relaxed state.

Let us now move to the test case of the reduced Euler equations. Bracket (63) with Hamiltonian (49) and Eq. (7) yields the evolution equation

$$\partialt u = -\left[u - u{\Omega} - \frac{\mathcal{H}(u)}{|\phi|_{L^2(\Omega)}^2}\phi\right],$$

$$46$$

with ϕ depending on u via the Poisson equation (46). The evolution of u is therefore governed by a nonlinear integral operator with a nonpolynomial nonlinearity. We begin by considering a numerical experiment. Figure 4 shows the initial and final state of a solution of the initial value problem for this equation, and Fig. 5 gives the representation of the functional relation between the potential ϕ and the vorticity ω = u − uΩ. The initial condition as well as the numerical method and the numerical parameters (grid size and time steps) are the same as in Fig. 2.

Figure 4 show the relative error in energy conservation and the entropy as a function of time. The thick horizontal line denotes the minimum entropy value in Eq. (56). This time the minimum entropy value is quickly reached by the system.

From the scatter plot in Fig. 5 one can see that the final state is characterized by a linear relation between ϕ and ω. This suggests that the projector-based metric bracket (63) relaxes the state of the system completely. This nice property comes at the price of a larger dissipation of entropy, cf. Fig. 4. As a consequence, the vorticity of the final state is significantly lower than in the initial condition. The scatter plot in Fig. 5 also shows the average of the initial condition on the contours of the initial potential (which is the same as in Fig. 3); in this case, the relaxed state bears little or no similarity to the initial condition.

In order to check if the relaxed vorticity is in agreement with the analytical solution (55), we have computed the best fit of the analytical solution (55) to the final state in Fig. 4, varying the three phases θ0, θ1, and θ2. The difference between the best fit and the relaxed state gives an estimate of the distance of the latter from the entropy minimum and it is shown in Fig 5, right-hand-side panel.

In the terminology of Section 3.3, we claim that bracket (63) is minimally degenerate. In fact, for any function F, we have (F, F) = 0 if and only if δF(u)/δu ∈ ker ΠH(u), or equivalently

$$\frac{\delta \mathcal{F}(u)}{\delta u} = \lambda \phi,$$

which is condition (43) with λ^α ̸= 0 only if I ^α = H.

If F = S, this condition is satisfied at any point of the set C^η defined in equation (50), i.e., for any constrained critical point of the entropy, not just at the minimum. It follows that the generalization (PL′′) of the PL condition cannot be true on the whole phase space. Nonetheless, with (58) we have

$$(S,S)(u) = 2S(u) - \frac{4\mathcal{H}0^2}{|\phi|{L^2(\Omega)}^2}.$$ (64)

We claim that for any u such that H(u) = H0, the following inequalities hold true:

$$2S(u) - \frac{4\mathcal{H}0^2}{|\phi|{L^2(\Omega)}^2} \ge 0, \tag{65a}$$

$$S \ge \mathcal{H}_0,$$ (65b)

$$|\phi|{L^2(\Omega)}^2 \le |\nabla\phi|{L^2(\Omega)}^2 = 2\mathcal{H}_0.$$ (65c)

Specifically, Eq. (65a) follows from the Cauchy-Schwarz inequality (which also ensure the positivity of the projector). Equation (65b) is a consequence of (56), while (65c) is

Figure 4: The same as in Fig. 2 but for the metric bracket (63).

a Poincar´e inequality on T , obtained from the solution of the Poisson problem (46) via Fourier series. Inequalities (65) imply that the evolution of the solution of the metriplectic system must be such that S(u), 1/∥ϕ∥ 2 L2(Ω) ∈ R 2 ⁺ remains within the cone

$$\mathcal{S}(u) \geq \mathcal{H}0 = \mathcal{S}{\eta}, \quad \frac{1}{2\mathcal{H}0} \leq \frac{1}{|\phi|{L^2(\Omega)}^2} \leq \frac{1}{2\mathcal{H}0} + \frac{\mathcal{S}(u) - \mathcal{S}{\eta}}{2\mathcal{H}_0^2}.$$

We observe that the set C^η is contained in the "upper boundary" of this cone, i.e., it satisfies

$$\frac{1}{|\phi|_{L^2(\Omega)}^2} = \frac{1}{2\mathcal{H}0} + \frac{\mathcal{S}(u) - \mathcal{S}{\eta}}{2\mathcal{H}_0^2}.$$

For any a ∈ (0, 1), and for any u ∈ U^η such that

$$\frac{1}{|\phi|_{L^2(\Omega)}^2} = \frac{1}{2\mathcal{H}0} + (1-a)\frac{S(u) - S\eta}{2\mathcal{H}_0^2},\tag{66}$$

Figure 5: Left-hand-side panel: the relation between ϕ and ω = u − u^Ω for the case of Fig. 4. Righthand-side panel: difference between the best fit of the exact solution (55) and the relaxed state.

we have

$$(S, S) \ge 2a[S(u) - S_{\eta}],$$

which is inequality (PL′′) with constant κ^η = 2a. Hence, if the solution stays within the shrunk cone (66), we expect exponential convergence with exponent related the angle of the cone (66).

Figure 6 shows the trajectory of the solution in the plane S(u), 1/∥ϕ∥ 2 L2(Ω) for the case of Fig. 4. The solution indeed stays within the cone defined by (65). The figure also demonstrate exponential convergence of both entropy and vorticity with exponential rate consistent with the generalized PL condition. We note that, cf. Fig. 4, the system traverses most of the trajectory in Fig. 6 quickly, and reaches the vertex of the cone already at t ≈ 5.

Figure 7 shows the same trajectory for a different initial condition, viz.,

$$u_0(x) = \cos(2x_2) + u_G(x),$$ (67)

where u^G is the Gaussian defined in (61) with 1/N = 1.8, x⁰ = (π, 3π/2), w¹ = 0.3, and w² = 1. The cosine terms shift the initial condition closer to the boundary. As a consequence the initial entropy relaxation rate is slower, but approaches ≈ 1 as the trajectory approaches the vertex of the cone.

5. Collision-like metric brackets

The specific structure of the metric bracket for the Landau collision operator was introduced and generalized in [1, 3]. Here we propose a further generalization that we use as a "template" for the construction of relaxation methods. This generalized bracket will be referred to as the collision-like metric bracket [89], as it originates from Morrison's bracket for the Landau collision operator. The resulting evolution equation is integrodifferential. We demonstrate the use of collision-like brackets for the solution of the variational principles in Eqs. (13) and (20) for equilibria of the reduced Euler equations and axisymmetric MHD, respectively. In both these applications the specification of

Figure 6: Left-hand-side panel: trajectory in the plane (S(u), 1/∥ϕ∥L2(Ω)) for the case of Fig. 4; the shaded area indicates the region of the plane excluded by inequalities (65). The dashed line indicates the upper boundary of the shrunk cone (66) with a = 1/2. The minimum entropy state corresponds to the vertex of the cone. Right-hand-side panel: evolution of the norm of the distance of ω(t) from the relaxed state, using the vorticity ω at the last point in time t = T as an approximation of the latter, and the "excess entropy" S u(t) − Sη. The semi-log scale shows exponential relaxation of entropy with exponential rate ≈ 1. This is consistent with the inequality (PL′′). The fact that ω has relaxation rate ≈ 1/2 is a consequence of the simple choice of the entropy function, cf. Eq. (47).

equilibrium profiles, i.e., the relation between the state variable u and the corresponding potential (either ϕ for the Euler equations or ψ for axisymmetric MHD), is essential. In the variational principle, the profile is encoded in the choice of the entropy; therefore, for the result of metriplectic relaxation to be consistent with the imposed profile, the metric bracket must completely relax the state of the system, in the sense made precise in Section 3. Preliminary, numerical results on these two problems have been reported in the proceedings of the joint Varenna-Lausanne workshop on "The Theory of Fusion Plasmas" [90] and in the Ph.D. thesis by Bressan [89].

5.1. General construction of collision-like brackets

We introduce the class of collision-like metric brackets in a fairly abstract way, but we make no attempt to give a mathematically rigorous definition, except for a few basic considerations. The construction given here is somewhat more general than the one proposed earlier [90, 89].

As in Section 2.1, the phase space V is a space of sufficiently regular functions v : Ω → R ^N over a bounded domain Ω ⊂ R ^d with N, d ∈ N. We always assume V ⊆ L (Ω, µ; R ^N ) =: W, and the functional derivative of a functional F ∈ C (V ) is computed with respect to the standard product in W, hence δF(u)/δu ∈ W, when it exists.

For the construction of the bracket, we consider a bounded domain O ⊂ R ⁿ equipped with a measure ν, and the space W˜ := L 2 (O, ν; R N˜ ) with n, N˜ ∈ N. For our purposes it is sufficient to assume that dν(z) = ˜m(z)dz with ˜m ∈ C^∞(O), dz being the Lebesgue measure on O. Next, we choose a (possibly unbounded) linear operator

$$P: W \to \tilde{W}, \quad \operatorname{dom}(P) = \Phi,$$

where the domain Φ is a subspace of W with a finer topology, so that V ⊆ Φ ⊆ W ⊆ Φ ′ with Φ′ being the dual of Φ (the space of continuous linear functionals on Φ), and with

Figure 7: The same as in Fig. 6, but for the initial condition (67). The black dots marked on the boundary of the cone correspond to the points of the set Cη.

continuous inclusions. Then, W has the structure of a rigged Hilbert space, with the finer space Φ containing the phase space V . In this way, both quadratic functionals like F(u) = R Ω u ²dµ and linear functionals like F(u) = R ^Ω w · udµ for w ∈ Φ are such that δF(u)/δu ∈ Φ. In addition, given a Hamiltonian H such that δH(u)/δu ∈ Φ, we choose

$$T: V \rightarrow \mathcal{B}(\tilde{W}),$$

with values in the space B(W˜ ) of bounded linear operators from W˜ into W˜ , such that T(u) is symmetric, positive semidefinite, and

$$\mathsf{T}(u)P\frac{\delta\mathcal{H}(u)}{\delta u} = 0. \tag{68}$$

In terms of the operator P and the function T, collision-like brackets that preserve H are defined by

$$(\mathcal{F}, \mathcal{G}) := \int_{\mathcal{O}} P \frac{\delta \mathcal{F}(u)}{\delta u} \cdot \mathsf{T}(u) P \frac{\delta \mathcal{G}(u)}{\delta u} d\nu.$$ (69)

Symmetry and positive semidefiniteness of the bracket is ensured by the fact that T(u) is symmetric and positive semidefinite, while Eq. (68) implies (F, H) = 0 for any F. Hence, Eq. (69) defines a metric bracket on the class of functions F such that δF(u)/δu ∈ Φ. In most cases, T(u) is the operator of multiplication by a function T(u; x) with values in the space of real symmetric, positive semidefinite N˜ × N˜ matrices.

Upon introducing the dual operator P ′ : W˜ → Φ ′ defined by

$$\langle w, P'\tilde{w} \rangle = \int_{\mathcal{O}} \tilde{w} \cdot Pw \, d\nu, \quad \text{for all } w \in \Phi,$$ (70)

where ⟨·, ·⟩: Φ × Φ ′ → R is the duality pairing between Φ and Φ′ , bracket (69) can be equivalently written as

$$(\mathcal{F}, \mathcal{G}) = \left\langle \frac{\delta \mathcal{F}(u)}{\delta u}, P' \left[ \mathsf{T}(u) P \frac{\delta \mathcal{G}(u)}{\delta u} \right] \right\rangle.$$

Figure 8: Construction of the operators P and P ′ in Eqs. (69) and (71).

If u ∈ C 1 ([0, T], V ) is a trajectory in V and F ∈ C 1 (V ) has a functional derivative δF(u)/δu in Φ, then t 7→ F u(t) is differentiable and

$$\frac{d}{dt}\mathcal{F}(u(t)) = \left\langle \frac{\delta \mathcal{F}(u(t))}{\delta u}, \partial_t u \right\rangle.$$

Therefore, given an entropy S such that δS(u)/δu ∈ Φ, Eq. (7a) for u(t) amounts to

$$\partial_t u = -P' \left[ \mathsf{T}(u) P \frac{\delta \mathcal{S}(u)}{\delta u} \right] \quad \text{in } \Phi'.$$ (71)

This construction is summarized in Fig. 8.

The operator P maps an R ^N -valued function over the d-dimensional domain Ω into an R N˜ -valued function over the n-dimensional domain O, and we are particularly interested in the case n > d, N˜ ≥ N. As we shall see, increasing the number of dimensions has some advantages that, however, come at the price of a higher computational cost. The function T will be referred to as the kernel of the bracket and, in general, it depends on both P and H, because of (68), but for simplicity, this dependence is not explicitly indicated in the notation.

There are of course many ways to choose T such that condition (68) is satisfied and we give two particularly relevant examples below. Among the various choices, we are interested in those that satisfy:

$$||w||_W \le CP ||Pw||{\tilde{W}}, \qquad w \in \text{dom}(P), \tag{72a}$$

$$\tilde{w} \in \ker \mathsf{T}(u) \cap \operatorname{ran} P \iff \tilde{w} = \lambda P \frac{\delta \mathcal{H}(u)}{\delta u}, \text{ for } \lambda \in \mathbb{R} \text{ constant},$$ (72b)

for a constant C^P . Here ker T(u) and ran P denote the null space of the operator T(u) and the range of the operator P, respectively, and they are both subspaces of W˜ . When P = ∇, Eq. (72a) is the Poincar´e inequality [91].

Formally at least, conditions (72) imply that the bracket defined in Eq. (69) is minimally degenerate in the sense of Eq. (43). In fact, (F, F) = 0 is equivalent to

$$P\frac{\delta \mathcal{F}(u)}{\delta u} \in \ker \mathsf{T}(u) \cap \operatorname{ran}(P),$$

$$52$$

and condition (72b) implies that there is a constant λ such that

$$P\frac{\delta \mathcal{F}(u)}{\delta u} = \lambda P\frac{\delta \mathcal{H}(u)}{\delta u}.$$

Then condition (72a) yields

$$\left| \frac{\delta \mathcal{F}(u)}{\delta u} - \lambda \frac{\delta \mathcal{H}(u)}{\delta u} \right|{W} \le C{P} \left| P \left[ \frac{\delta \mathcal{F}(u)}{\delta u} - \lambda \frac{\delta \mathcal{H}(u)}{\delta u} \right] \right|_{\tilde{W}} = 0,$$

which implies condition (43) with λ^α ̸= 0 only if I ^α = H.

It follows that metric brackets of the form (69), with the defining operator P and the kernel T satisfying (72) could be used to construct a relaxation method for variational problems of the form (1).

We now give examples of (69) for which conditions (72) are, at least formally, satisfied.

Example 5. Morrison's brackets for the Landau collision operator [1] and its generalization [3] can be obtained as special cases of (69). This motivates our choice of the name collision-like.

The configuration space V is the space of particle distribution functions u(x) = f(x, v), where x is the spatial position and v is the velocity of the particles. We assume that x ∈ Ω^x and v ∈ Ωv, with both domains Ωx, Ω^v ⊂ R ^D being bounded, D ≥ 2, and f(x, ·) compactly supported in Ωv, i.e., the velocity-space domain is large enough to contain all particle velocities. Then, Ω = Ω^x ×Ωv, d = 2D, and N = 1, since f is a scalar field. The measure on Ω is the Lebesgue measure dµ(x) = dxdv.

We choose O = Ω × Ωv, n = 3D, dν(x, v, v ′ ) = dxdvdv ′ , and for g ∈ L 2 (Ω, µ; R) with ∇vg ∈ L (Ω, µ, ; R ^D),

$$Pg(x, v, v') = \nabla_v g(x, v) - \nabla_v g(x, v')$$

hence P g(x, v, v ′ ) ∈ R ^N˜ with N˜ = D. With the Hamiltonian

$$\mathcal{H}(f) = \int_{\Omega} \left[ \frac{1}{2} m \mathbf{v}^2 + \mathbf{V}(\mathbf{x}) \right] f(\mathbf{x}, \mathbf{v}) d\mathbf{x} d\mathbf{v},$$

where m is the mass of the considered particle species and V is a potential energy, we have P(δH(f)/δf) = m(v−v ′ ). Therefore, a possible choice of the operator T(u) = T(f) is the multiplication by the matrix-valued kernel

$$T{L}(f; \mathbf{x}, \mathbf{v}, \mathbf{v}') = \frac{\nu{c}}{2} M(f(\mathbf{x}, \mathbf{v})) M(f(\mathbf{x}, \mathbf{v}')) U_{L}(\mathbf{v} - \mathbf{v}'), \tag{73}$$

where ν^c > 0 is a constant collision frequency, M : R⁺ → R⁺ is arbitrary, and

$$U_{\mathrm{L}}(\mathsf{v}-\mathsf{v}') \coloneqq \frac{1}{|\mathsf{v}-\mathsf{v}'|} \Big( I - \frac{(\mathsf{v}-\mathsf{v}') \otimes (\mathsf{v}-\mathsf{v}')}{|\mathsf{v}-\mathsf{v}'|^2} \Big).$$

Condition (68) holds, and bracket (69) reduces to

$$(\mathcal{F}, \mathcal{G}) = \frac{\nuc}{2} \int{\mathcal{O}} M(f(\mathbf{x}, \mathbf{v})) M(f(\mathbf{x}, \mathbf{v}')) \left[ \nabla{\mathbf{v}} \frac{\delta \mathcal{F}(f)}{\delta f}(\mathbf{x}, \mathbf{v}) - \nabla{\mathbf{v}} \frac{\delta \mathcal{F}(f)}{\delta f}(\mathbf{x}, \mathbf{v}') \right] \cdot U{\mathbf{L}}(\mathbf{v} - \mathbf{v}') \left[ \nabla{\mathbf{v}} \frac{\delta \mathcal{G}(f)}{\delta f}(\mathbf{x}, \mathbf{v}) - \nabla_{\mathbf{v}} \frac{\delta \mathcal{G}(f)}{\delta f}(\mathbf{x}, \mathbf{v}') \right] d\mathbf{x} d\mathbf{v} d\mathbf{v}', \quad (74)$$

which is the bracket of Eq. (44) in Morrison's paper [3]. One can notice that in Eq. (73) the orthogonal projection onto the direction of

$$m(\mathbf{v}-\mathbf{v}') = \nabla{\mathbf{v}} \frac{\delta \mathcal{H}(f)}{\delta f}(\mathbf{x},\mathbf{v}) - \nabla{\mathbf{v}} \frac{\delta \mathcal{H}(f)}{\delta f}(\mathbf{x},\mathbf{v}')$$

ensures property (68). In most examples considered in this paper, the kernel of the bracket is constructed from a similar projector. We also note that, in the construction given here, there is no need to use a singular distributional kernel, cf. the Dirac's delta in Eq. (45) of [3].

The operator P ′ defined in Eq. (70) acting on a function ˜w(z) = ˜g(x, v, v ′ ) can be formally computed after integration by parts, with the result that

$$P'\tilde{g}(\mathbf{x},\mathbf{v}) = -\operatorname{div}{\mathbf{v}}\Big[\int{\Omega_{\mathbf{v}}} \big(\tilde{g}(\mathbf{x},\mathbf{v},\mathbf{w}) - \tilde{g}(\mathbf{x},\mathbf{w},\mathbf{v})\big) d\mathbf{w},\Big],$$

and the evolution equation (71) takes the form

$$\begin{split} \partialt f &= \mathrm{div}\mathbf{v} \left[ \nuc \int{\Omega\mathbf{v}} M \big( f(\mathbf{x}, \mathbf{v}) \big) M \big( f(\mathbf{x}, \mathbf{w}) \big) \ &\quad \times \ U\mathrm{L}(\mathbf{v} - \mathbf{w}) \Big( \nabla\mathbf{v} \frac{\delta \mathcal{S}(f)}{\delta f}(\mathbf{x}, \mathbf{v}) - \nabla\mathbf{v} \frac{\delta \mathcal{S}(f)}{\delta f}(\mathbf{x}, \mathbf{w}) \Big) d\mathbf{w} \right], \end{split}$$

which reduces to the Landau operator for Coulomb collisions when S(f) = R f log f dxdv and M(f) = f.

As for conditions (72), given g(x, v) in the domain of the operator P, we have

$$\int{\Omega{\mathbf{v}}} \nabla_{\mathbf{v}} g(\mathbf{x}, \mathbf{v}) d\mathbf{v} = 0,$$

since, per hypothesis, g(x, ·) vanishes near the boundary of Ωv. Then,

$$|Pg|{\tilde{W}}^2 = 2|\Omega{\mathbf{v}}| \int{\Omega} \left|\nabla{\mathbf{v}} g(\mathbf{x}, \mathbf{v})\right|^2 d\mathbf{x} d\mathbf{v},$$

where |Ωv| = R Ω^v dv ′ . The minimum of ∥P g∥ 2 W˜ subject to the constraint ∥g∥^W = 1 is attained for

$$-\Delta_{\mathbf{v}}g = \lambda0 g, \qquad \int{\Omega} |g(\mathbf{x}, \mathbf{v})|^2 d\mathbf{x} d\mathbf{v} = 1,$$

where λ⁰ is the minimum eigenvalue of the Laplace operator ∆^v on the velocity domain Ω^v with homogeneous Dirichlet boundary conditions. Hence

$$||Pg||_{\tilde{W}}^2 = 2|\Omegav| \int{\Omega} |\nabla_v g(x, v)|^2 dx dv \ge 2|\Omega_v|\lambda_0||g||_W^2,$$

which shows (formally) that condition (72a) holds. As for condition (72b), a function g˜(x, v, v ′ ) belongs to ker T^L when

$$\tilde{g}(x, v, v') = \Lambda(x, v, v')(v - v'),$$

almost everywhere (a.e.) in O, where Λ(x, v, v ′ ) ∈ R. If, in addition, ˜g ∈ ran(P), we must have

$$\nabla{\mathbf{v}}g(\mathbf{x},\mathbf{v}) - \nabla{\mathbf{v}}g(\mathbf{x},\mathbf{v}') = \Lambda(\mathbf{x},\mathbf{v},\mathbf{v}')(\mathbf{v}-\mathbf{v}'),$$ a.e. in $\mathcal{O}$ .

If g and Λ satisfying this condition exist, then necessarily Λ(x, v, v ′ ) = Λ(x, v ′ , v) and, upon fixing a point v ′ = a ∈ Ωv,

$$\nabla{\mathbf{v}}g(\mathbf{x},\mathbf{v}) = \nabla{\mathbf{v}}g(\mathbf{x},\mathbf{a}) + \Lambda(\mathbf{x},\mathbf{v},\mathbf{a})(\mathbf{v}-\mathbf{a}),$$

and

$$\begin{split} \nabla{\mathbf{v}}g(\mathbf{x},\mathbf{v}) - \nabla{\mathbf{v}}g(\mathbf{x},\mathbf{v}') &= \Lambda(\mathbf{x},\mathbf{v},\mathbf{a})(\mathbf{v}-\mathbf{a}) - \Lambda(\mathbf{x},\mathbf{v}',\mathbf{a})(\mathbf{v}'-\mathbf{a}) \ &= \Lambda(\mathbf{x},\mathbf{v},\mathbf{v}')\big((\mathbf{v}-\mathbf{a}) - (\mathbf{v}'-\mathbf{a})). \end{split}$$

Then we must have

$$\big\Lambda(\mathsf{x},\mathsf{v},\mathsf{v}')-\Lambda(\mathsf{x},\mathsf{v},\mathsf{a})\big-\big\Lambda(\mathsf{x},\mathsf{v},\mathsf{v}')-\Lambda(\mathsf{x},\mathsf{a},\mathsf{v}')\big=0.$$

The two vectors v−a and v ′−a in R ^D are linearly dependent only if they are proportional and that can only happen for a set of points (x, v, v ′ ) of measure zero in O. Hence, we deduce that Λ must satisfy the necessary conditions

$$\begin{cases} \Lambda(\textbf{x},\textbf{v},\textbf{v}') - \Lambda(\textbf{x},\textbf{v},\textbf{a}) = 0, \ \Lambda(\textbf{x},\textbf{v},\textbf{v}') - \Lambda(\textbf{x},\textbf{a},\textbf{v}') = 0, \end{cases} \quad \text{a.e. in } \mathcal{O},$$

and this must hold for almost any choice of the arbitrary point a. This is possible only if

$$\Lambda(x, v, v') = \lambda = \text{constant}$$ , a.e. in $\mathcal{O}$ ,

which gives Eq. (72b). Therefore, at least formally, the bracket for the Landau operator is minimally degenerate in the sense introduced in Section 3.3.

Example 6. The bracket based on the L 2 -orthogonal projector discussed in Section 4.2 can be obtained as a special case of (69). In fact, we can utilize Eq. (69) in order to generalize (63) to the case of R ^N -valued fields.

With this aim, let Ω ⊂ R d , N ∈ N, O = Ω × Ω, and N˜ = 2N, that is, we double both the dimension of the domain and of the field. Then, with Φ = W = L 2 (Ω, µ; R ^N ) = Φ′ , W˜ = L 2 (O, ν; R N˜ ), and dν(x, x′ ) = dµ(x)dµ(x ′ ), we define

$$P: W \to \tilde{W}, \qquad Pw(x, x') = \begin{pmatrix} w(x) \ w(x') \end{pmatrix} \in \mathbb{R}^{2N}.$$

Given a Hamiltonian function H(u) with δH(u)/δu ∈ W, we choose T(u) to be the multiplication operator by the kernel

$$\mathsf{T}(u;x,x') = \kappa(u) \begin{pmatrix} |h(x')|^2 I_N & -h(x) \otimes h(x') \ -h(x') \otimes h(x) & |h(x)|^2 I_N \end{pmatrix},$$

where I^N is the N × N identity block, and h = δH(u)/δu is a short-hand notation for the functional derivative (h may still depend on u). One can check that condition (68) is satisfied. Inserting these choices into Eq. (69) yields

$$\begin{split} (\mathcal{F},\mathcal{G}) &= 2\kappa(u) \int{\mathcal{O}} \frac{\delta \mathcal{F}(u)}{\delta u}(x) \cdot \Big[ \Big| \frac{\delta \mathcal{H}(u)}{\delta u}(x') \Big|^2 \frac{\delta \mathcal{G}(u)}{\delta u}(x) \ &\qquad \qquad - \frac{\delta \mathcal{H}(u)}{\delta u}(x) \Big( \frac{\delta \mathcal{H}(u)}{\delta u}(x') \cdot \frac{\delta \mathcal{G}(u)}{\delta u}(x') \Big) \Big] d\mu(x) d\mu(x') \ &= 2\kappa(u) \Big| \frac{\delta \mathcal{H}(u)}{\delta u} \Big|{L^2}^2 \int{\Omega} \frac{\delta \mathcal{F}(u)}{\delta u}(x) \cdot \Pi{\mathcal{H}} \frac{\delta \mathcal{G}(u)}{\delta u}(x) d\mu(x), \end{split}$$

where Π^H is the L 2 -orthogonal projector defined in Eq. (62), but generalized to the case of R ^N -valued fields. If 2κ is set to the inverse of ∥h∥ 2 ^L² , this reduces to (63) when N = 1. Condition (72a) follows from

$$|Pw|{L^2(\mathcal{O},\nu;\mathbb{R}^{2N})}^2 = \int{\Omega\times\Omega} \left(|w(x)|^2 + |w(x')|^2\right) d\mu(x) d\mu(x') = 2\mu(\Omega) |w|_{L^2(\Omega,\mu;\mathbb{R}^N)}^2,$$

where µ(Ω) = R dµ. As for condition (72b), if h ̸= 0, ˜w ∈ kerT(u)∩ran(P) implies that there is w ∈ W such that ˜w = P w, and

$$|h(x')|^2 w(x) - h(x) (h(x') \cdot w(x')) = 0 \in \mathbb{R}^N,$$

for almost all (x, x′ ). Upon multiplying by h(x), we have

$$|h(x')|^2 (h(x) \cdot w(x)) = |h(x)|^2 (h(x') \cdot w(x')).$$

This is equivalent to saying that − |h(x ′ 2 , |h(x)| 2 and h(x) · w(x), h(x ′ ) · w(x ′ ) are orthogonal in R 2 , or that

$$\begin{pmatrix} h(x) \cdot w(x) \ h(x') \cdot w(x') \end{pmatrix} = \Lambda(x, x') \begin{pmatrix} |h(x)|^2 \ |h(x')|^2 \end{pmatrix},$$

for a function Λ, which is general depends on (x, x′ ) and this holds for almost all (x, x′ ). We deduce that Λ must be constant almost everywhere in O, i.e. Λ(x, x′ ) = λ for almost all (x, x′ ) and thus

$$\lambda |h(x)|^2 = h(x) \cdot w(x) \iff w(x) = \lambda h(x),$$

which implies condition (72b). Hence, metriplectic brackets constructed by means of an L 2 -orthogonal projection are minimally degenerate as already shown in Section 4.2 for the scalar case (N = 1).

Example 7. As in example 6, let O = Ω × Ω, n = 2d, with coordinates z = (x, x′ ) ∈ Ω × Ω, and dν(x, x′ ) = dµ(x)dµ(x ′ ). The operator P is given by

$$Pw(x, x') = P_L w(x, x') = Lw(x) - Lw(x'),$$

where L : W → L (Ω, µ; R N˜ ) is a possibly unbounded linear operator with domain dom(L) = Φ ⊆ W and taking values in the space of R N˜ -valued, squared-integral functions, with N˜ ≥ 2. In general, we allow N˜ to be different from N. Specific cases will be considered in Secs. 5.2 and 5.3 below.

As for the kernel, a way to satisfy condition (68) utilizes the matrix

$$Q_k(z) := |z|^2 I_k - z \otimes z, \qquad z \in \mathbb{R}^k,$$ (75)

where k ∈ N and I^k is the k×k identity block. Specifically, let T(u) be the multiplication operator by the matrix

$$\mathsf{T}(u;x,x') = \frac{1}{2}\kappa(u;x,x')Q_{\tilde{N}}\Big(P_L\frac{\delta\mathcal{H}(u)}{\delta u}(x,x')\Big),\tag{76}$$

where κ(u; .) is a positive scalar weight function satisfying the symmetry condition κ(u; x, x′ ) = κ(u; x ′ , x). (This choice of the kernel is trivial if N˜ = 1, which is excluded.) With the foregoing choices, Eq. (69) becomes

$$(\mathcal{F}, \mathcal{G}) := \frac{1}{2} \int{\Omega} \int{\Omega} \kappa(u; x, x') \Big[ PL \frac{\delta \mathcal{F}(u)}{\delta u}(x, x') \Big] \cdot Q{\bar{N}} \Big( P_L \frac{\delta \mathcal{H}(u)}{\delta u}(x, x') \Big) \Big[ P_L \frac{\delta \mathcal{G}(u)}{\delta u}(x, x') \Big] d\mu(x) d\mu(x'). \tag{77}$$

This is the form of collision-like brackets as originally formulated by Bressan et al. [90, 89]. The dual operator P ′ can be determined in terms of the dual of L, which is the linear operator L ′ : L (Ω, µ; R N˜ ) → Φ ′ defined for all ϖ ∈ L (Ω, µ; R N˜ ) by

$$\langle w, L'\varpi \rangle = \int_{\Omega} \varpi(x) \cdot Lw(x) d\mu(x), \text{ for all } w \in \Phi,$$ (78)

and we find

$$\begin{split} \int_{\mathcal{O}} \tilde{w}(x,x') \cdot PL w(x,x') d\nu(x,x') &= \int{\Omega} L w(x) \int{\Omega} \left[ \tilde{w}(x,x') - \tilde{w}(x',x) \right] d\mu(x') d\mu(x) \ &= \left\langle w, L' \Big[ \int{\Omega} \left[ \tilde{w}(\cdot,x') - \tilde{w}(x',\cdot) \right] d\mu(x') \Big] \right\rangle \ &= \left\langle w, P_L' \tilde{w} \right\rangle. \end{split}$$

In terms of the operator L ′ , the evolution equation (71) reads

$$\partial_t u = -L' \Big[ \mathbb{D}(u) L \frac{\delta \mathcal{S}(u)}{\delta u} - \mathbb{F} \Big( u, \frac{\delta \mathcal{S}(u)}{\delta u} \Big) \Big], \tag{79}$$

where we have defined

$$\mathbb{D}(u;x) := \int{\Omega} \kappa(u;x,x') Q{\tilde{N}} \Big( P_L \frac{\delta \mathcal{H}(u)}{\delta u}(x,x') \Big) d\mu(x'), \tag{80a}$$

$$\mathbb{F}(u, v; x) := \int{\Omega} \kappa(u; x, x') Q{\tilde{N}} \left( P_L \frac{\delta \mathcal{H}(u)}{\delta u}(x, x') \right) Lv(x') d\mu(x'). \tag{80b}$$

Eq. (79) has the same structure as the Landau operator for Coulomb collision discussed in Example 5, with ∇^v being replaced by L.

At last, we check conditions (72). Since dµ(x) = m(x)dx, with m continuous on the closed domain Ω, and dν(x, x′ ) = dµ(x)dµ(x ′ ), we have

$$|Pw|{\tilde{W}}^2 = \int{\mathcal{O}} \left|Lw(x) - Lw(x')\right|^2 d\nu(x, x') \ge m0^2 \int{\mathcal{O}} \left|Lw(x) - Lw(x')\right|^2 dx dx',$$

where m⁰ = min{m(x): x ∈ Ω}. For condition (72a) to hold, it is therefore sufficient that L satisfies

$$||w||_W \le CL ||Lw||{L^2(\Omega)}, \qquad \int_{\Omega} Lw(x)dx = 0,$$ (81)

so that

$$\int{\mathcal{O}}\left|Lw(x)-Lw(x')\right|^2dxdx'=2|\Omega|\int{\Omega}\left|Lw(x)|^2dx-2\right|\int_{\Omega}Lw(x)dx\bigg|^2\geq 2|\Omega|C_L^{-2}|w|_W^2,$$

which give Eq. (72a). All considered cases in the applications below satisfy condition (81).

As for the kernel of the bracket, we argue as in the case of Example 5. With the kernel given in Eq. (76), a function ˜w ∈ W˜ belongs to ker T(u) ∩ ran(P) only if there is a function w ∈ W and Λ : Ω × Ω → R, such that

$$Lw(x) - Lw(x') = \Lambda(x, x') [Lh(x) - Lh(x')],$$

where h = δH(u)/δu and necessarily Λ(x, x′ ) = Λ(x ′ , x). If we fix an arbitrary point x ′ = a ∈ Ω, we obtain an explicit expression for Lw(x),

$$Lw(x) = Lw(a) + \Lambda(x, a)X_a(x), \qquad X_a(x) := Lh(x) - Lh(a).$$

Therefore

$$Lw(x) - Lw(x') = \Lambda(x, a)X_a(x) - \Lambda(a, x')X_a(x')$$ $$= \Lambda(x, x')[X_a(x) - X_a(x')],$$

or equivalently

$$\left[\Lambda(x,x')-\Lambda(x,a)\right]X_a(x)-\left[\Lambda(x,x')-\Lambda(a,x')\right]X_a(x')=0.$$

Under the assumption that Xa(x) and Xa(x ′ ) are linearly independent for almost all (x, x′ ) ∈ Ω × Ω = O, we deduce Λ(x, x′ ) = λ = constant as in Example 5. Hence, if the operator L satisfies (81) and the Hamiltonian is such that Xa(x) and Xa(x ′ ) are linearly independent almost everywhere in Ω × Ω, the bracket (77) is minimally degenerate.

In general we observe that the brackets of the form discussed in Example 7 are all special cases of the metriplectic 4-bracket introduced recently by Morrison and Updike [23], since the kernel depends quadratically on δH/δu. Specifically, they can be obtained from the metriplectic 4-bracket constructed by utilizing the Kulkarni-Nomizu product as discussed in Section III.D.1 of Ref. [23].

Example 7 will be used as a template for the construction of dissipative operators for the solution of some of the variational problems introduced in Section 2.2. Before addressing the applications, we specialize bracket (77) for two relevant choices of the operator L.

5.2. Collision-like brackets based on div–grad operators

In Eq. (77), let us consider the case of s scalar field (N = 1) and choose

$$Lw = \nabla w,$$

$$58$$

with domain dom(L) = H¹ 0 (Ω), the space of functions in L 2 (Ω) with weak derivatives also in L 2 (Ω) and satisfying homogeneous Dirichlet boundary conditions on ∂Ω. This operator satisfies condition (81), and N˜ = d (the inequality, in particular, is the classical Poincar´e inequality for H¹ 0 , cf. Theorem 3 in Section 5.6.1 of Ref. [91]). The dual operator L ′ can be obtained from Eq. (78). Given ϖ ∈ L 2 (Ω, µ; R d ), sufficiently regular, we have

$$\langle w, L'\varpi\rangle = \int{\Omega}\varpi(x)\cdot\nabla w(x)d\mu(x) = -\int{\Omega}w(x)\operatorname{div}_{\mu}\varpi(x)d\mu(x),$$

and thus −L ′ϖ = div^µ ϖ = 1 m P i ∂xⁱ mϖⁱ is the divergence operator associated to the volume form dµ(x) = m(x)dx. For a generic ϖ, we write L ′ ⁼ ⁻divg^µ to denote that the divergence is defined weakly. As for the kernel, we choose

$$\kappa(u; x, x') = M(x, u(x))M(x, u(x')),$$

where M(x, y) > 0 is a positive function over Ω × R, which can be used to simplify the bracket [3] as in Example 5.

As for the entropy, we shall restrict to functions of the form

$$S(u) = \int_{\Omega} s(x, u(x)) dx, \tag{82}$$

where s(x, y) is a given profile, convex in y, and possibly depending on position x. Then

$$\frac{\delta S(u)}{\delta u}(x) = \partial_y s(x, u(x)),$$

and

$$\nabla \frac{\delta \mathcal{S}(u)}{\delta u}(x) = \partial_y^2 s\big(x, u(x)\big) \nabla u(x) + \partial_x \partial_y s\big(x, u(x)\big).$$

The arbitrary function M can be chosen so that [3]

$$M(x,y)\partial_y^2 s(x,y) = 1, (83)$$

under the condition ∂ y s(x, y) > 0, consistently with the convexity assumption for the profile. Eq. (83) introduces a dependence of the kernel of the bracket on the entropy function. With this choice, Eq. (79) becomes

$$\partialt u = \widetilde{\operatorname{div}}{\mu} \Big[ \mathbb{D}_s(u) \big( \nabla u + M(x, u) \partial_x \partial_y s(x, u) \big) - M(x, u) \mathbb{F}_s(u, \nabla u) \Big], \tag{84}$$

and

$$\begin{split} \mathbb{D}{s}(u;x) &= \int{\Omega} Q{2} \Big( \nabla \frac{\delta \mathcal{H}(u)}{\delta u}(x) - \nabla \frac{\delta \mathcal{H}(u)}{\delta u}(x') \Big) M \big( x', u(x') \big) d\mu(x') \ \mathbb{F}{s}(u;x) &= \int{\Omega} Q{2} \Big( \nabla \frac{\delta \mathcal{H}(u)}{\delta u}(x) - \nabla \frac{\delta \mathcal{H}(u)}{\delta u}(x') \Big) \ & \big[ \nabla u(x') + M \big( x', u(x') \big) \partial{x} \partial{y} s \big( x', u(x') \big) \big] d\mu(x'). \end{split}$$

Eq. (84) still depends on the choice of the measure µ, the entropy profile s, and the Hamiltonian function H. This choices will be specified below separately for the cases of the reduced Euler and Grad-Shafranov equations in two-dimensions (d = N˜ = 2).

5.3. Collision-like brackets based on curl–curl operators

With magnetic fields in mind, we consider an example of bracket (69) for divergencefree fields. Then, Ω is bounded, simply connected domain in R ³ with sufficiently regular connected boundary, d = N = 3, and dµ(x) = dx. We choose

$$Lw = \operatorname{curl} w$$ ,

with domain

$$\Phi = { w \in H(\text{curl}, \Omega) \cap H(\text{div}, \Omega) : \text{ div } w = 0 \text{ in } \Omega, \ n \times w = 0 \text{ on } \partial \Omega },$$

where H(curl, Ω) and H(div, Ω) are the spaces of vector fields w ∈ L 2 (Ω; R 3 ) such that curl w ∈ L (Ω; R 3 ) and div w ∈ L (Ω), respectively. Specifically, w ∈ Φ = dom(L) is a divergence-free field with zero tangential component n × w on the boundary, and n : ∂Ω → R is the outward unit normal on ∂Ω. It follows that condition (81) is satisfied: the inequality amounts to the Poincar´e inequality for divergence-free fields on a simply connected domain with connected boundary (cf. Corollary 3.51 in Ref.[92]), while

$$\int{\Omega} \operatorname{curl} w dx = \int{\partial \Omega} (n \times w) d\sigma = 0,$$

where dσ is the surface element on ∂Ω.

If ϖ ∈ L (Ω; R 3 ) is sufficiently regular, Eq. (78) and integration by parts gives

$$\langle w, L'\varpi \rangle = \int{\Omega} \varpi(x) \cdot \operatorname{curl} w(x) dx = \int{\Omega} w(x) \cdot \operatorname{curl} \varpi(x) dx,$$

hence L ′ϖ = curl ϖ, while for a generic ϖ, L ′ ⁼ curl is the weak curl operator. ^g

For the kernel (76), we choose κ = 1, and the corresponding evolution equation can be written as

$$\partial_t u = -\widetilde{\operatorname{curl}} \Big[ \mathbb{D}_v(u) \operatorname{curl} \frac{\delta \mathcal{S}(u)}{\delta u} - \mathbb{F}_v \Big( u, \frac{\delta \mathcal{S}(u)}{\delta u} \Big) \Big], \tag{85}$$

where

$$\mathbb{D}{v}(u;x) = \int{\Omega} Q_{3} \left( \operatorname{curl} \frac{\delta \mathcal{H}(u)}{\delta u}(x) - \operatorname{curl} \frac{\delta \mathcal{H}(u)}{\delta u}(x') \right) dx',$$

$$\mathbb{F}{v}(u,v;x) = \int{\Omega} Q_{3} \left( \operatorname{curl} \frac{\delta \mathcal{H}(u)}{\delta u}(x) - \operatorname{curl} \frac{\delta \mathcal{H}(u)}{\delta u}(x') \right) \operatorname{curl} v(x') dx'$$

A non-trivial factor κ of the form used in Section 5.2 can easily be accounted for [89].

The evolution equation (85) preserves the divergence-free constraint. In fact, for any φ ∈ H¹ (Ω) the function Fφ(u) = R u · ∇φdx = − R φ div udx is a constant of motion.

5.4. Application to the reduced Euler equations

In order to demonstrate the properties of the collision-like brackets, we construct a relaxation method for the solution of the variational principle (13) for equilibria of the reduced Euler equations, with entropy and Hamiltonian functions given in Eq. (14).

This is the same problem addressed in Section 4, with the difference that here we consider a bounded domain Ω with homogeneous Dirichlet boundary conditions as discussed in Section 2.2.1. Specifically, the domain is the unit square Ω = [0, 1] × [0, 1].

When, in Eq. (14), s(ω) = ω ²/2, we known from Section 2.2.1 that the solution of the variational principle (13) is related to the eigenfunction of the Laplacian operator on Ω corresponding to the smallest eigenvalue. On the unit square with homogeneous Dirichlet boundary conditions, the eigenfunctions and the corresponding eigenvalues are given by

$$\phi{m,n}(x) = \mathcal{N}{m,n} \sin(m\pi x_1) \sin(n\pi x2), \qquad \lambda{m,n} = \pi^2(m^2 + n^2),$$

labeled by two non-zero integers m, n ∈ N and normalized by a non-zero constant Nm,n. The smallest eigenvalue corresponds to the function ϕ1,1, and energy conservation allows us to determine the normalization constant N1,1, that is

$$2\mathcal{H}0 = \int{\Omega} \phi \omega dx = \lambda{1,1} |\phi{1,1}|{L^2(\Omega)}^2 = \lambda{1,1} \mathcal{N}_{1,1}^2 / 4,$$

from which we deduce the unique solution of (13) with this entropy, that is

$$\phi(x) = \left(2\sqrt{\mathcal{H}_0}/\pi\right)\sin(\pi x_1)\sin(\pi x2), \omega(x) = \lambda{1,1}\phi(x) = 4\pi\sqrt{\mathcal{H}_0}\sin(\pi x_1)\sin(\pi x2),$$ $\lambda{1,1} = 2\pi^2 \approx 19.7392.$ (86)

More general choices of the entropy density s(ω) lead to the nonlinear eigenvalue problems of the form (15). Then no analytical solution is known, but an estimate of the eigenvalue λ can be found from energy conservation in a similar way. Upon multiplying by ω Eq. (15) and integrating over Ω, we find

$$\int{\Omega} \omega s'(\omega) dx = \lambda \int{\Omega} \omega \phi dx = 2\lambda \mathcal{H}_0,$$

from which we can obtain λ provided that the integral on the left-hand side can be evaluated, e.g. from the numerical solution. Specifically, if s(ω) = ω log ω we have

$$\lambda = \frac{1}{2\mathcal{H}_0} (\mathcal{M}(\omega) + \mathcal{S}(\omega)), \tag{87}$$

where, in particular, M(ω) = R Ω ωdx. These analytical results can be used to assess the result of the proposed relaxation method.

As for the construction of the method itself, we consider the collision-like metriplectic system of Section 5.2. The state variable u(t) ∈ V ⊆ Φ is identified with vorticity, u(t, x) = ω(t, x), and evolved according to Eq. (84) with S and H given in Eq. (14).

The numerical scheme for the solution of Eq. (84) is based on H¹ -conforming finite elements, i.e., the discrete solution uh(t) belongs to the same space H¹ (Ω) that contains the phase space V . The scheme preserves the discrete Hamiltonian H(uh) (modulo roundoff errors) and we use the Crank-Nicolson scheme in time, which preserves the property of monotonic dissipation of entropy in the quadratic case (s(ω) = ω ²/2). More details on the derivation of the scheme can be found in Ref. [89].

Figure 9: Relaxation of an initial vortex with initial vorticity given in Eq. (88), according to the metriplectic system (84) with the state variable u being the fluid scalar vorticity ω. The initial condition is shown in the left-hand-side panel: the color map represents the vorticity field ω and the dashed (white) lines the contours of the corresponding streaming function ϕ, cf. Eq. (10). The relaxed state is displayed in the middle panel. The right-hand-side panel represents the functional relation between ω and ϕ for the initial condition (black dots), the final state (red bullets), and the expected linear relation for the exact solution (86) (green crosses).

Figure 10: Evolution of entropy (left-hand-side panel) and of the variation of the Hamiltonian relative to its initial value (right-hand-side panel), for the case in Fig. 9.

The obtained numerical method has been implemented in the finite-element library FEniCS [93, 94], in which the weak form of the operator in Eq. (84) can be directly specified by means of the domain-specific language UFL [95], and discretized by the finite-element component DOLFIN [96, 97]. Among the various tests [89], here we discuss in details three cases only.

Single vortex. In the simplest example, the domain Ω = [0, 1]² is discretized by a uniform grid of 64 × 64 nodes. The vorticity field u = ω is approximated in the space of second-order Lagrange elements that are available in FEniCS [94]. The entropy density is quadratic, i.e. s(ω) = ω ²/2, hence the analytical solution of the variational principle is given in Eq. (86). The initial condition is

$$u_0(x) = u_G(x)$$ , with $u_G$ defined in Eq. (61), (88)

and with parameters x0,¹ = x0,² = 1/2, w ¹ = 0.01, w 2 ² = 0.07, and N = 1. The numerical scheme provably preserves the Hamiltonian, which in this case is the kinetic energy of the fluid, and dissipates monotonically the entropy, independently of the magnitude of the time step.

Figure 11: The same as Fig. 9, but for the initial condition (89).

Figure 12: Evolution of entropy (left-hand-side panel) and of the variation of the Hamiltonian relative to its initial value (right-hand-side panel), for the case in Fig. 11.

Fig. 9 shows the initial condition, the final state, and the "scatter plot", which we use to identify functional relations between different fields, cf. the analysis in Section 4. Qualitatively, we see that the initial condition is quite far from an equilibrium of the Euler equations as the contours of the streaming function ϕ and those of vorticity ω are misaligned. Metriplectic relaxation with collision-like brackets yields a solution that appears to be an equilibrium, and from the scatter plot (right-hand-side panel in Fig. 9), one can see that the relaxed state is indeed an equilibrium characterized by the same linear relation of the exact solution (86). Therefore, the collision-like metric bracket appears to have completely relaxed the initial solution in the sense discussed in Section 3. From Fig. 10, one can verify energy conservation and entropy monotonic dissipation.

Perturbed equilibrium. We repeat the experiment of Fig. 9 with an initial condition close to a critical point of entropy restricted to the constant-Hamiltonian surface. Specifically, the initial condition is

$$u_0(x) = \sin(6\pi x_1)\sin(4\pi x_2) + u_G(x)$$ , with $u_G$ defined in Eq. (61), (89)

and with the same parameters as in Eq. (88) except for N = 100.

Fig. 11 shows the initial condition, the final state, and the usual scatter plot, which visualizes the relationship between ω and ϕ. The initial condition (Fig. 11, left-handside panel) is basically an eigenfunction of the Laplace operator corresponding to a large eigenvalue, the perturbation being hardly visible. This is confirmed by the scatter plot (Fig. 11 right-hand-side panel, black dots) in which the initial state is concentrated on

Figure 13: The same as Fig. 9, but with the entropy density s(ω) = ω log ω. For the reference solution (green crosses) we used Eq. (90) with λ computed from Eq. (87).

a straight line, with a small spread due to the Gaussian perturbation. Therefore the initial condition is close to an equilibrium. The relaxed state (Fig. 11, middle panel) is not exactly the same as in Fig. 9, since the initial value H⁰ of the Hamiltonian is different, but is it consistent with the exact solution (86) as shown in the scatter plot. Therefore, the initial condition has been completely relaxed to a solution of variational principle (13). In Fig. 12, one can see that the Hamiltonian function is preserved to machine accuracy, while entropy decays monotonically as expected. However, initially entropy appears to remain constant, due to the proximity of the initial condition to an equilibrium point.

Gibbs entropy density. We consider now a more complicated entropy function, namely, the Gibbs entropy, which is given by the entropy density s(ω) = ω log ω. In this case the expected relaxed state is determined by, cf. Eq. (15),

$$1 + \log \omega = \lambda \phi \iff \omega = e^{\lambda \phi - 1},$$ (90)

where the eigenvalue λ could be numerically estimated by means of Eq. (87). The initial condition is the same as the one in Eq. (88) except for the amplitude, which here is increased to 1/N = 10.

Fig. 13 shows the initial condition, the final state after the relaxation and the usual scatter plot. In this case the solution of the solution of the variational principle (green crosses) is computed using the numerically computed values of M(ω) and S(ω) in Eq. (87). Again we see evidence of complete relaxation of the system toward the solution of the variational principle. Fig. 14 shows the expected behavior of entropy and Hamiltonian functions. It is worth noting that since entropy is not quadratic in ω, the numerical scheme does not preserve the property of monotonic entropy dissipation, hence sufficiently small time steps must be used to ensure the the evolution of the system is approximated with sufficient accuracy.

5.5. Application to Grad-Shafranov equilibria

As a second example, we construct a relaxation method for axisymmetric MHD equilibria, cf. Section 2.2.2. Essentially, this amounts to a different iterative method for the solution of the Grad-Shafranov equation. The metriplectic relaxation method ensures

Figure 14: Evolution of entropy (left-hand-side panel) and of the variation of the Hamiltonian relative to its initial value (right-hand-side panel), for the case in Fig. 13.

conservation of the poloidal magnetic energy and monotonic dissipation of an "ad hoc" entropy, but at a higher computational cost.

Specifically, we construct a relaxation method for the variational principle (20). On a bounded domain Ω, satisfying Ω ⊂ R⁺ × R with coordinates x = (x1, x2) = (r, z), cf. Section 2.2.2, the state variable is a scalar field u(t) proportional to the toroidal component of the plasma current, u(t, r, z) = (4π/c)rJφ(t, r, z), and it is evolved in time according to Eq. (84) as in the case of reduced Euler equations.

The entropy and Hamiltonian functions are chosen as in Eq. (21), with entropy density and measure µ on Ω given by

$$s(r,y) = \frac{1}{2} \frac{y^2}{Cr^2 + D}, \quad d\mu(r,z) = \frac{1}{r} dr dz,$$ (91)

and the assumptions on the domain imply r ≥ r⁰ > 0 in Ω. Then, the first condition in Eq. (23) defines the toroidal current

$$\frac{4\pi}{c}J_{\varphi} = \lambda \Big(Cr + \frac{D}{r}\Big)\psi,$$

which is the current of the equilibrium found by Herrnegger [98] and Maschke [99], cf. also Mc Carthy [100, Section II.B]. Equivalently, since the state variable is u = (4π/c)rJφ, the condition in Eq. (23) can be written as

$$\frac{u}{Cr^2 + D} = \lambda \psi. {92}$$

In order to have a reference solution, we resort to the direct numerical solution of the Grad-Shafranov equation by using the classical iterative scheme [60, pp. 22–23, Eqs. (2.111) and (2.112)]. In the following experiments C = 0.6 and D = 0.2.

Rectangular domain. We started with a rectangular domain Ω = [1, 7] × [−9.5, +9.5], with coordinates x = (x1, x2) = (r, z), discretized by a uniform grid of 64 × 64 nodes. The initial condition u⁰ for the state variable is the same as in Eq. (88) with parameters x0,¹ = r⁰ = 4, x0,² = z⁰ = 0, w 2 ¹ = 0.5, w 2 ² = 3.2, and N = 1.

Fig. 15 shows the results of this numerical experiment. Instead or plotting the state variable u directly, the color plot represents the field u/(Cr² + D), which should be proportional to ψ if the system reaches a state consistent with Eq. (92). Qualitatively

Figure 15: Relaxation of an initial current J^φ = (c/4πr)u0, with u⁰ Gaussian, according to the evolution equation (84) applied to the Grad-Shafranov problem (20) with entropy (91) on a rectangular domain. The initial condition and the final state of the system are given in left-hand-side and middle panels, respectively, while the right-hand side panel shows the scatter plot, with the same color/symbol code as in Fig. 9. The color map represents the field u/(Cr² + D), and the white contours the flux function ψ, so that condition (92) is easily checked. Analogously the axes in the scatter plot refer to the values of the flux function ψ and the field u/(Cr² + D). Here, C = 0.6 and D = 0.2.

Figure 16: Evolution of entropy (left-hand-side panel) and of the variation of the Hamiltonian relative to its initial value (right-hand-side panel), for the case in Fig. 15.

the results are similar to those of Fig. 9 for the reduced Euler equations: the initial condition evolves toward a state consistent with Eq. (92). Fig. 16 confirms the expected behavior of the entropy and Hamiltonian functions.

Mapped domain. At last, we give an example of the relaxation method for the Grad-Shafranov equation on a non-trivial mapped domain with a smooth boundary. The domain Ω is obtained by mapping the unit disk {z ∈ C: |z| < 1} with the map defined by

$$r = a \left[ b + \frac{1}{\varepsilon} \left( 1 - \sqrt{1 + \varepsilon(\varepsilon + 2s\cos\theta)} \right) \right],$$

$$z = c \frac{e\xi s\sin\theta}{2 - \sqrt{1 + \varepsilon(\varepsilon + 2s\cos\theta)}},$$ (93)

where z = s exp(iθ) is a point in the unit disk, and the parameters are e = 1.4, ε = 0.3, a = 4, b = 3, and c = 6.3, with ξ = 1/ p 1 − ε ²/4. This map is a slightly modified version of the one used by Zoni and G¨ucl¨u [101, Eq. (3) and references therein]. The initial condition is the same as for the rectangular domain, except for r⁰ = 12, w ¹ = 0.6, and w ² = 6.0.

Figure 17: The same as in Fig. 15, but on the domain obtained mapping the unit disk with Eq. (93).

Figure 18: Evolution of entropy (left-hand-side panel) and of the variation of the Hamiltonian relative to its initial value (right-hand-side panel), for the case in Fig. 17.

The results for the mapped domain are shown in Fig. 17. The color map represents the field u/(Cr² + D), as in the previous case. The relaxed state is again consistent with Eq. (92). Fig. 18 shows the expected monotonic entropy dissipation and the conservation of the Hamiltonian to machine precision.

6. Diffusion-like metric brackets

A feature of the general collision-like bracket is that the generalized diffusion tensor and friction flux (80) are nonlocal functions of the unknown u. Therefore, their evaluation requires an integration over the whole domain Ω. This poses an issue of computational complexity even harder than that of the standard Landau collision operator, which is local in half of the variables (cf. Example 5). In order to mitigate the computational cost of a relaxation method based on these brackets, we have studied the special case of brackets (69) that corresponds to choosing O = Ω, i.e., we do not increase the size of the domain (n = d). We will refer to this case as a diffusion-like bracket, and reserve the name "collision-like" for the case n > d. As we shall see, in general, one cannot expect complete relaxation (in the sense defined in Section 3) from the diffusion-like brackets. For instance, the metric double bracket (57) is a special case of a diffusion-like bracket.

6.1. General construction of the brackets

With the same setup of Section 5.1, which is summarized in Fig. 8, let us consider the case O = Ω, ν = µ, but we allow N˜ ̸= N, so that in general W˜ ̸= W.

Then bracket (69) reduces to

$$(\mathcal{F}, \mathcal{G}) := \int_{\Omega} P \frac{\delta \mathcal{F}(u)}{\delta u} \cdot \mathsf{T}(u) P \frac{\delta \mathcal{G}(u)}{\delta u} d\mu, \tag{94}$$

where P : W → W˜ is a linear (possibly unbounded) operator, and T(u) ∈ B(W˜ ) is symmetric, positive semidefinite, bounded linear operator that satisfies (68). The evolution equation generated by the bracket (94) is formally the same as Eq. (71), but we shall see in the examples that the dual operator P ′ does not involve any integral operator.

Conditions (72) are of course still sufficient conditions for brackets of the form (94) to be minimally degenerate, but we shall see that condition (72b), in particular, is usually not satisfied in this case.

Example 8. Brackets of the form (57), i.e., metriplectic double brackets acting on scalar fields, are special cases of diffusion-like brackets (94). In order to see this, let N = 1, thus W = L 2 (Ω, µ), Φ = C^∞(Ω), and

$$Pw = \nabla w$$ .

Hence, N˜ = d and W˜ = L (Ω, µ; R d ) is the space of L ² vector fields over Ω. Given a function J(x) of class C^∞(Ω) taking values in the space of antisymmetric d×d matrices, we can define the antisymmetric bilinear operation

$$[w_1, w_2]_J = \nabla w_1(x) \cdot J(x) \nabla w_2(x).$$

As in Section 4.1, it is not necessary that [·, ·]^J satisfies the Jacobi identity. In addition let us consider a symmetric, positive-definite bi-linear form γ : W × W → R, together with the associated linear bounded, symmetric positive definite operator Γ: W → W, that is, cf. Appendix B,

$$\gamma(w_1, w2) = \int{\Omega} w_1(x) \cdot \Gamma w_2(x) d\mu(x),$$

for all w1, w² ∈ W. In terms of J and γ, we define the kernel

$$\mathsf{T}(u) = X_h(u) \circ \Gamma \circ {}^tX_h(u),$$

where Xh(u) and ^tXh(u) are the operators of multiplication by the vector fields

$$X_h(u;x) = J(x)\nabla h(x), \quad {}^tX_h(u;x) = -\nabla h(x)\cdot J(x),$$

respectively, and h = δH(u)/δu ∈ Φ = C^∞(Ω). If J is a Poisson tensor, i.e., [·, ·]^J satisfies the Jacobi identity, Xh(u; ·) is the d-dimensional Hamiltonian vector field generated by the Hamiltonian function h = δH(u)/δu. Since Xh(u; ·) ∈ C^∞(Ω; R d ), T(u) maps W˜ into itself. The operator Γ is symmetric and positive definite, hence T(u) is symmetric and positive semidefinite, and we have ∇h ∈ kerT(u) as required by condition (68).

With the foregoing choices of P and T(u), Eq. (94) becomes

$$\begin{split} (\mathcal{F},\mathcal{G}) &= \int_{\Omega} \Big( \nabla \frac{\delta \mathcal{F}(u)}{\delta u} \cdot J \nabla \frac{\delta \mathcal{H}(u)}{\delta u} \Big) \Gamma \Big( \nabla \frac{\delta \mathcal{G}(u)}{\delta u} \cdot J \nabla \frac{\delta \mathcal{H}(u)}{\delta u} \Big) d\mu \ &= \gamma \Big( \Big[ \frac{\delta \mathcal{F}(u)}{\delta u}, \frac{\delta \mathcal{H}(u)}{\delta u} \Big]_J, \Big[ \frac{\delta \mathcal{G}(u)}{\delta u}, \frac{\delta \mathcal{H}(u)}{\delta u} \Big]_J \Big), \ 68 \end{split}$$

which is Eq. (57). As discussed in Section 4.1, this bracket is in general not minimally degenerate. In fact, while condition (72a) amounts to the Poincar´e inequality and holds true on a bounded domain Ω, condition (72b) fails since for any sufficiently regular function f : R → R, any function of the form ˜w = ∇f(h), with h = δH(u)δu, belongs to kerT(u) ∩ ran(P).

6.2. Diffusion-like brackets based on div–grad operators

We address the diffusion-like version of the bracket introduced in Section 5.2. For scalar fields (N = 1) on a bounded domain Ω ⊂ R d , d ≥ 2, let P = ∇ with domain Φ = H¹ 0 (Ω). The kernel of the bracket is constructed from the matrix Qd, cf. Eq. (75),

$$\mathsf{T}(u;x) = \kappa(u;x)Q_d\Big(\nabla \frac{\delta \mathcal{H}(u)}{\delta u}(x)\Big),$$

where κ(u; x) is a positive scalar function, and T(u) is defined as the operator of multiplication by T(u; ·). Then, bracket (94) becomes [89]

$$(\mathcal{F}, \mathcal{G}) = \int_{\Omega} \kappa(u) \nabla \frac{\delta \mathcal{F}(u)}{\delta u} \cdot Q_d \left( \nabla \frac{\delta \mathcal{H}(u)}{\delta u} \right) \nabla \frac{\delta \mathcal{G}(u)}{\delta u} d\mu,$$ (95)

and the corresponding evolution equation is

$$\partialt u = \widetilde{\operatorname{div}}{\mu} \left[ \kappa(u) Q_d \left( \nabla \frac{\delta \mathcal{H}(u)}{\delta u} \right) \nabla \frac{\delta \mathcal{S}(u)}{\delta u} \right] \quad \text{in } \Phi' = H^{-1}(\Omega).$$

This is a "local version" of Eq. (84) which justifies the name "diffusion-like" for this bracket. Condition (72b) fails in the same way as in Example 8.

It is worth noting that in two spatial dimensions, d = 2, one has

$$Q_2(\nabla h) = X_h \otimes X_h = X_h^t X_h, \qquad h = \delta \mathcal{H}(u)/\delta u,$$

where X^h = Jc∇h = (−∂2h, ∂1h) is the canonical Hamiltonian vector field generated by h(x) and ^tX^h denotes its transpose. This means that the diffusion-like bracket (95) in a two dimensional domain amounts to the metric double bracket addressed in Example 8, with Γ = I being the identity operator. For d > 2 the bracket (95) is however different from the metric double brackets in Example 8. In fact, if X^h ̸= 0, the null space of the matrix Q^d is always one dimensional for any dimension d, while the null space of Xh⊗X^h is d − 1 dimensional, hence the two matrices have the same null space only if d = 2.

Nonetheless, for d ≥ 3, one can write the matrix Qd(∇h) in terms of suitable pairing of two antisymmetric operations by using the identity

$$\frac{1}{(d-2)!} \sum_{i1,\dots,i{d-2}} \epsilon_{i1,\dots,i{d-2},i,k} \epsilon_{i1,\dots,i{d-2},j,l} = \delta{ij} \delta{kl} - \delta{il} \delta{jk},$$

with ϵⁱ1,...,i^d being the completely antisymmetric symbol. We obtain

$$[Qd(\nabla h)]{ij} = |\nabla h|^2 \delta_{ij} - \partial_i h \partialj h = \sum{kl} \left[ \delta{ij} \delta{kl} - \delta{il} \delta{jk} \right] \partial_k h \partial_l h$$

$$= \frac{1}{(d-2)!} \sum_{i1, \dots, i{d-2}} \sum_{k,l} \epsilon^{i1, \dots, i{d-2}, i, k} \epsilon^{i1, \dots, i{d-2}, j, l} \partial_k h \partial_l h,$$ (96)

where ∂jh = ∂h/∂x^j , and thus, Eq. (95) can be written equivalently as

$$(\mathcal{F},\mathcal{G}) = \frac{1}{(d-2)!} \sum{\alpha} \int{\Omega} \kappa(u) \, \mathcal{E}_d^{\alpha}(\nabla f, \nabla h) \, \mathcal{E}_d^{\alpha}(\nabla g, \nabla h) \, d\mu,$$

where α = (i1, . . . , id−2) is a multi-index, f = δF(u)/δu, g = δG(u)/δu, and

$$\mathcal{E}d^{\alpha}(\nabla \varphi, \nabla \psi) = \sum{i,k} \epsilon^{i1, \dots, i{d-2}, i, k} \partial_i \varphi \partial_k \psi, \quad \alpha = (i1, \dots, i{d-2}).$$

In dimension d = 3, E³ defines a Lie bracket in R 3 . This is the standard Lie algebra structure on R ³ given by the cross product arising in the case of rigid body rotation [3].

Yet another form of this bracket makes use of the Kulkarni-Nomizu (K-N) product and the metriplectic 4-bracket structure [23]. Using the first identity in (96), we can write

$$(\mathcal{F},\mathcal{G}) = \frac{1}{2} \sum{i,j,k,l} \int{\Omega} \kappa(u) (\delta \otimes \delta)_{ijkl} \left[ \partial_i \frac{\delta \mathcal{F}(u)}{\delta u} \right] \left[ \partial_j \frac{\delta \mathcal{H}(u)}{\delta u} \right] \left[ \partial_k \frac{\delta \mathcal{G}(u)}{\delta u} \right] \left[ \partial_l \frac{\delta \mathcal{H}(u)}{\delta u} \right] d\mu,$$

where (^δ ? ^δ)ijkl = 2(δikδjl ⁻ ^δilδjk) is the K-N product of two identity tensors.

6.3. Diffusion-like brackets based on curl–curl operators

As a last example, we address the diffusion-like version of the curl-curl brackets of Section 5.3. We consider a vector field over a bounded domain Ω ⊂ R ³ with Lebesgue measure dµ(x) = dx, hence d = N = 3. We choose

$$Pw = \operatorname{curl} w,$$

so that N˜ = N = 3 and W˜ = W, with dom(P) given by

$$\Phi = { w \in H(\operatorname{curl}, \Omega) \cap H(\operatorname{div}, \Omega) : \operatorname{div} w = 0 \text{ in } \Omega, \ n \cdot w = 0 \text{ on } \partial \Omega }.$$

This differs from the space Φ considered in Section 5.3 by the "opposite" choice of boundary conditions: the normal component is set to zero instead of the tangential component. The Poincar´e inequality for the operator curl holds for this space as well [92, Corollary 3.51], so that condition (72a) holds true. (Here, we have the choice of the boundary condition since we do not need to satisfy the second identity in Eq. (81).)

As an example, let the kernel be once again constructed from Q3, defined in Eq. (75),

$$\mathsf{T}(u;x) = \kappa(u;x)Q_3\Big(\operatorname{curl}\frac{\delta\mathcal{H}(u)}{\delta u}\Big),$$

where κ(u; x) is a positive function. Even though the operator P = curl with domain dom(P) = Φ satisfies a Poincar´e inequality, in general, the kernel fails to satisfy condition (72b): a function ˜w ∈ kerT(u) ∩ ran(P) ⊂ W˜ = L (Ω; R 3 ) must satisfy ˜w = curl w, with w ∈ Φ and ˜w(x) = Λ(x)b(x), b = curl h, h = δH/δu, and thus the pair (Λ, w) must solve

$$\begin{cases} \operatorname{curl} w = \Lambda b, & \text{in } \Omega, \ \operatorname{div} w = 0, & \text{in } \Omega, \ b \cdot \nabla \Lambda = 0, & \text{in } \Omega, \ n \cdot w = 0, & \text{on } \partial \Omega. \end{cases}$$

As a special case let b be a nonlinear Beltrami field, i.e., a solution of (24) such that curl b = µb with µ(x) not a constant, then Λ = µ and w = b is a solution that violates condition (72b).

With the foregoing choices, Eq. (94) amounts to

$$(\mathcal{F}, \mathcal{G}) = \int{\Omega} \kappa(u) \operatorname{curl} \frac{\delta \mathcal{F}(u)}{\delta u} \cdot Q{3} \left( \operatorname{curl} \frac{\delta \mathcal{H}(u)}{\delta u} \right) \operatorname{curl} \frac{\delta \mathcal{G}(u)}{\delta u} dx, \tag{97}$$

and the corresponding evolution equation becomes

$$\partial_t u = -\widetilde{\operatorname{curl}} \left[ \kappa(u) Q_3 \left( \operatorname{curl} \frac{\delta \mathcal{H}(u)}{\delta u} \right) \operatorname{curl} \frac{\delta \mathcal{S}(u)}{\delta u} \right] \quad \text{in } \Phi' = H'(\operatorname{curl}, \Omega).$$

This bracket can be written in terms of the antisymmetric bilinear operator

$$\mathcal{E}3(X,Y) = [X,Y]{\mathbb{R}^3} := X \times Y, \qquad X,Y \in \mathbb{R}^3,$$

which is the standard Lie bracket in R 3 . In fact Eq. (97) can be shown to be a special case of the following (cf. [3, 51]):

$$(\mathcal{F}, \mathcal{G}) = \int{\Omega} \left[ \operatorname{curl} \frac{\delta \mathcal{F}(u)}{\delta u}, \operatorname{curl} \frac{\delta \mathcal{H}(u)}{\delta u} \right]{\mathbb{R}^{3}} \Gamma \left[ \operatorname{curl} \frac{\delta \mathcal{G}(u)}{\delta u}, \operatorname{curl} \frac{\delta \mathcal{H}(u)}{\delta u} \right]_{\mathbb{R}^{3}} dx, \tag{98}$$

where Γ(u) ∈ B L 2 (Ω; R 3 ) is a symmetric, positive definite operator; Eq. (97) is obtained for Γ(u) = κ(u), the multiplication operator by the function κ(u; ·). Equation (98) is a metric double bracket of the form (57). Applied to magnetic fields this gives a generalization of the relaxation method of Chodura and Schl¨uter [18] with constant pressure, cf. also Moffatt [14]. An explicit example will be briefly reported in Section 6.4 below.

6.4. Application to nonlinear Beltrami fields

So far we have focused on examples of equilibrium problems for which complete relaxation of the solution is essential. We have shown that a metriplectic relaxation method for such problems should be based on metric brackets that are minimally degenerate (or specifically degenerate if more than one constraint is considered), cf. Section 3. Diffusionlike brackets do not appear to be appropriate for those problems.

For sake of completeness, we address an example of equilibrium problems that are characterized as minima of a function subject to topological constraints. This is the case of nonlinear Beltrami fields, for which the variational principle is given in Lagrangian representation, cf. Section 2.2.3 and Appendix C. Full three-dimensional MHD equilibria satisfy the same type of Lagrangian variational principle.

Because of their larger null space, metric double brackets of the form (98) allow us to obtain an evolution equation that preserves the necessary constraints.

We identify the state variable u(t) with a magnetic field u(t, x) = B(t, x) on a simply connected, bounded domain Ω ⊂ R 3 . More specifically, we assume that B(t) ∈ V ⊂ Φ, where Φ is the same space defined in Section 5.3. The evolution equation for B(t) is given by Eq. (7a) and bracket (98), with Γ = I, the identity operator, for simplicity, and with entropy and Hamiltonian given in Eq. (27). Therefore, if an orbit of this metriplectic system completely relaxes, it would converge in time to a linear Beltrami field, cf. Section 2.2.3. We shall see that this, in general, cannot happen (due to the null space of the bracket), and that the limit must be a nonlinear Beltrami field.

The resulting evolution equation amounts to

$$\partial_t B = \widetilde{\operatorname{curl}} [B \times (B \times \operatorname{curl} B)],$$

where we have accounted for the identity curl[δH(B)/δB] = curl A = B. If B is sufficiently regular, we can replace curl by curl and write ^g

$$\begin{cases} \partial_t B - \operatorname{curl} \left[ V \times B \right] = 0, & \text{in } \Omega, \ V - (\operatorname{curl} B) \times B = 0, & \text{in } \Omega, \ n \cdot B = 0, & \text{on } \partial\Omega, \end{cases}$$ (99)

which shows that the magnetic field B is advected by the flow of the effective "velocity" field V . Hence, the field lines of B are frozen into the flow, i.e., they cannot change their topological properties. This is a much stronger constraint than just preservation of magnetic helicity 2H(B). Therefore, the solution can relax to a linear Beltrami field only if the initial condition is a continuous deformation of a linear Beltrami field, and this is typically not the case. The equilibrium points correspond to V = 0, which is equivalent to Eq. (24) for nonlinear Beltrami fields.

As anticipated, Eq. (99) is exactly the relaxation method of Chodura and Schl¨uter [18] with constant pressure. The method itself is therefore not new. In solar physics this relaxation method is known as the magneto-frictional method [102, 103, 104, 105], and it has been applied to the computation of force-free magnetic field in coronal active regions [19]. The bracket formalism, however, opens the way to possible generalizations by means of different choices of the kernel. This possibility will be explored in future work. Here we present a single numerical experiment obtained by a structure-preserving numerical scheme [89], which we derived by adapting the finite-element exterior calculus (FEEC) scheme of Hu et al. [106] for incompressible MHD. The scheme provably preserves the Hamiltonian (magnetic helicity), the constraint div B = 0, and the monotonic behavior of entropy (magnetic energy), but it does not preserve the topology of the field lines exactly. Similar work has been independently published by He et al. [107] earlier this year. Previously, the magnetic relaxation problem has been dealt with by means of Lagrangian [108] and finite difference [109] methods. More recently, a different kind of Lagrangian numerical scheme has been proposed [110, 111], which is based on the discretization of the domain in narrow flux tubes, each one being relaxed by a curveshortening flow in a modified metric. This interesting scheme therefore preserves the topological properties of the field lines. Yet with the domain discretized by a finite set of lines, the reconstruction of the magnetic field at arbitrary points of the domain, needs to be addressed.

Before describing the considered test case, let us address the role of magnetic-helicity conservation. In a domain Ω where the Poincar´e inequality for the curl operator holds true, magnetic helicity Hm(B) = 2H(B) provides a lower bound for magnetic energy. In fact, one has [13]

$$|Hm(B)| \le ||A||{L^2(\Omega)} ||B||_{L^2(\Omega)} \le CP ||B||{L^2(\Omega)}^2.$$ (100)

For an initial condition B⁰ with Hm(B0) = 0, it is possible that the solution of (99) with the chosen boundary conditions (B·n = 0 on ∂Ω) relaxes to a trivial field, i.e., |B(t)| → 0 for t → +∞, even if the topology of the field lines is preserved. This is the case for the class of one-dimensional solutions of (99), which are obtained, for instance, by assuming

$$B(t,x) = \begin{pmatrix} 0 \ 0 \ b(x_1) \end{pmatrix} = \operatorname{curl} \begin{pmatrix} 0 \ a(x_1) \ 0 \end{pmatrix},$$

where a ′ (x1) = b(x1), x = (x1, x2, x3) and the field is constant in (x2, x3). We have A · B = 0 and thus Hm(B) = 0. Correspondingly, equations (99) reduce to

$$\partial_t b - (b^2 b')' = 0,$$

where a prime denotes spatial differentiation. This is a standard heat equation. The steady states are solution to b b ′ = constant, which gives b(x1) = (c¹ + c2x1) 1/3 , with c1, c² being integration constants. For instance homogeneous boundary conditions for b on an interval yields the unique solution b(x1) = 0. Magnetic relaxation in one dimension has been recently considered by Yeates [112] and compared to the corresponding full MHD relaxation, thus exposing the limitations of the magneto-frictional method.

It is therefore meaningful to consider initial conditions with non-trivial magnetic helicity H(B) = ¹ ²Hm(B) ̸= 0. We construct such an initial condition from the vector potential

$$\tilde{A}(x) := a \begin{pmatrix} (n/\sqrt{m^2 + n^2}) \sin(\pi m x_1) \cos(\pi n x_2) \ -(m/\sqrt{m^2 + n^2}) \cos(\pi m x_1) \sin(\pi n x_2) \ \sin(\pi m x_1) \sin(\pi n x_2) \end{pmatrix},$$

with a ∈ R and m, n ∈ N being parameters (we shall choose m = n = 1). For any choice of the parameters, A˜ is divergence-free and a linear Beltrami field, periodic is all variables; it is an eigenvalue of curl corresponding to the eigenvalue λm,n = π(m² + n 2 1/2 . We localize this field in the unit cube Ω = [0, 1]³ by means of the cut-off function

$$\eta(x) := \chi(x_1)\chi(x_2)\chi(x_3), \qquad \chi(y) := y^2(1-y)^2, \quad y \in [0,1].$$

We have η(x) = 0 and ∇η(x) = 0 for x ∈ ∂Ω since both χ and its derivative χ ′ vanish for y = 0 and y = 1. We construct the initial condition on the domain Ω = [0, 1]³ as

$$\begin{split} A_0 &:= \eta \tilde{A}, \ B_0 &:= \operatorname{curl} A_0 = \nabla \eta \times \tilde{A} + \eta \operatorname{curl} \tilde{A}, \end{split} \tag{101}$$

and A⁰ ∈ H0(curl, Ω), B⁰ ∈ H0(div, Ω) with div B⁰ = 0. As for magnetic helicity,

$$H_m(B_0) = 2\mathcal{H}(B0) = \int{\Omega} A_0 \cdot B0 dx = \lambda{m,n} ||A0||{L^2(\Omega)}^2 > 0.$$

We can control the initial helicity by means of the parameters a ∈ R and m, n ∈ N. The magnetic field (101) is represented in Fig. 19 (top row), by means of a Poincar´e plot using the plane x² = 1/2 as a Poincar´e section. From the plot (Fig. 19, top-left panel), one can identify a rather complex topology of the field lines, with, in particular, four large islands of period two that are rendered in three dimensions in Fig 19, top-right panel, by tracing a few selected field lines for each island.

Figure 19: Poincar´e plot in the plane x1-x² and selected field lines of the magnetic field B, for the initial condition (top row) and the final state (bottom row), after the relaxation process. The initial condition is given in Eq. (101) with m = n = 1 and a = 1, while the evolution equation is the magneto-frictional method (99). The selected field lines correspond to the four large islands visible in the Poincar´e plot around x¹ = 1/2. The Poincar´e section is defined by x² = 1/2 and it is shown in light gray in the panels on the right-hand side. The initial points of the field lines are sampled differently for the initial and finals state of the field, hence they are not exactly the evolution of one another.

Figure 20: From top to bottom, time evolution of the entropy (magnetic energy), the relative variation of the Hamiltonian (magnetic helicity), the L² norm of the divergence, and of the vector j × H, where ^j ⁼ curl ^g^B is the current density (computed with the weak curl operator), and ^H is the ^L² -orthogonal projection of B onto H0(curl, Ω). The equilibrium condition for the considered numerical scheme reads j × H = 0. Corners and jumps in the time traces corresponds to restarts with larger time steps.

The time evolution of the initial condition (101) is obtained numerically by means of the FEEC scheme, which has been implemented in FEniCS as in the case of the tests reported in Section 5. Here, we use a relatively course resolution: the domain Ω = [0, 1]³ has been discretized by a uniform grid of 32³ nodes. The time step is adapted, but limited to a maximum of 10−³ . The obtained final state is represented in Fig. 19 (bottom row), again by means of a Poincar´e plot, Fig. 19, bottom-left panel. The four large islands appear to have been preserved by the relaxation process and are rendered in Fig. 19, bottom-right panel. Other large islands appear to have been preserved, but a closer analysis shows that the field line topology is not exactly preserved [89]. Indeed the numerical scheme preserves magnetic helicity only, and this alone does not completely guarantee the exact preservation of the field line topology.

In using Poincar´e plots for the visualization of the field line topology, one should address the effect of the error of the projection onto the finite-element space used for the representation of the magnetic field. In this case, B is approximated in the space of linear Raviart-Thomas elements for computations, and the discrete approximation is further projected onto the space of linear Lagrange elements, for visualization purposes (precise definitions of these finite-element spaces can be found, e.g., in the FEniCS book [94]). Lagrange elements are nodal so that the degrees of freedom coincide with the value of the field at the grid nodes and can be directly interpolated. We have qualitatively checked the effect of all those operations on the results by comparing the Poincar´e sections of the analytical field (101) with that of its projection onto the finite element space on the considered grid.

Figure 20 shows the evolution in time of the main quantities of interest, namely, the entropy, the relative variation of the Hamiltonian with respect to its initial value, and the L 2 -norms of div ^B and of the vector ^j ^× ^H, where ^j ⁼ curl ^g^B is (proportional to) the current density (computed weakly), while H is the L 2 -orthogonal projection of B ∈ H0(div, Ω) onto the space H0(curl, Ω). For sufficiently regular fields, we have j × H = (curl B) × B, but in general B ∈ H0(div, Ω) and H ∈ H0(curl, Ω) are different and the numerical equilibrium condition is j × H = 0. We recall that in this application the entropy and the Hamiltonian coincide with the magnetic energy and the magnetic helicity, respectively. From Fig. 20, we see that the qualitative properties of the relaxation method are preserved: the entropy decreases monotonically, the Hamiltonian is constant within a relative error of 10⁻10, and div B = 0 within an absolute error of 2 × 10⁻⁶ measured by the norm in L 2 (Ω). We also see that j × H decreases, which indicates that the solution is approaching a configuration that satisfies the Beltrami condition j × H = 0.

7. Conclusions

We have considered the question of whether, given a metriplectic system with Hamiltonian function H and (dissipated) entropy S, the orbit with initial condition u0, in the long-time limit, converges to a minimum of S on the surface of constant Hamiltonian {u: H(u) = H(u0)}. This question is interesting in itself, since many physical systems are metriplectic, but our work is mainly motivated by the idea of utilizing artificial metriplectic systems as relaxation methods for the computation of equilibria of fluids and plasmas.

We have shown that, in general, the answer is negative. For finite-dimensional metriplectic systems, we have given a sufficient condition, Proposition 5, under which an orbit relaxes to a constrained entropy minimum. These results are proven by means of a natural extension of the Lyapunov stability theorem for systems with constants of motion. One key assumption in Proposition 5 is that the metric bracket should be specifically degenerate, or minimally degenerate if the Hamiltonian is the only constant of motion. In addition, we have generalized the Polyak– Lojasiewicz condition for the exponential convergence of gradient flows, to the case of metriplectic system. The finitedimensional results have been extended to the infinite-dimensional case without proof, and supported by a number of examples only. In Section 4, we have studied quite in details a specific equilibrium problem for the reduced Euler equation. We constructed two relaxation methods based upon two metriplectic systems: one is specifically degenerate and the other one is not. The results of our numerical experiments with these two relaxation mechanisms can be explained in terms of the theoretical results.

A class of minimally degenerate metriplectic systems have been constructed in section 5. Such systems have been put forward as a generalization of the Landau collision operator, which is included in the class as a special case. For this reason we propose the term "collision-like" brackets. As a by-product, we obtain that the metric bracket for the Landau collision operator is minimally degenerate. We demonstrate the use of such system as the basis for relaxation methods for various equilibrium problems for both the reduced Euler equations and axisymmetric MHD equilibria, the latter being equivalent to solving the Grad-Shafranov equation. These are well-known equilibrium problems, for which various methods of solution exist, and are used here only as a proof of concept. From purely computational point of view, the direct solution of the Grad-Shafranov equation, in particular, is usually faster than the relaxation method constructed here, but the latter provides a recipe that can be adapted to more complicated equilibrium problems.

At last, we have discussed a simplified class of brackets, in which the nonlocal nature of collision-like brackets is removed. These brackets lead to diffusion-like evolution equations and reduce to a number of known brackets in special cases. Without the nonlocality characteristic of collision-like operators, these diffusion-like brackets are not minimally degenerate, but they can still be used to construct relaxation methods for equilibrium problems with stronger constraints. In the simplest case, we recover the known magnetic relaxation method of Chodura and Schl¨uter, and we give a numerical example based upon a structure-preserving numerical scheme obtained via finite-element exterior calculus (FEEC).

Appendix A. Bilinear forms and Leibniz identity

In this appendix we recall the basic definition of a Poisson structure and give a formal derivation of (5). The material is standard but not always easy to find in textbooks. Let V ⊆ L 2 (R d , µ; R ^N ) be a Banach space of squared-integrable functions over a domain Ω ⊆ R ^d with values in R ^N . We define functional derivatives of a function in C 1 (V ) as an element of L 2 (Ω, µ; R ^N ) according (2) with the pairing given by the L 2 scalar product,

$$\langle w, v \rangle = \int_{\Omega} w(x) \cdot v(x) d\mu(x),$$ 77

for all v ∈ V and w ∈ L 2 (Ω, µ; R ^N ), where µ is a measure on Ω, e.g., the Lebesgue measure.

We consider in particular the consequences of the Leibniz identity on the structure of the bracket. Let α : C∞(V ) × C∞(V ) → C∞(V ) be a generic bi-linear map satisfying the symmetry condition

$$\alpha(\mathcal{F}, \mathcal{G}) = \sigma \alpha(\mathcal{G}, \mathcal{F}), \quad \sigma^2 = 1,$$ (A.1)

and the Leibniz identity,

$$\alpha(\mathcal{F}, \mathcal{GH}) = \alpha(\mathcal{F}, \mathcal{G})\mathcal{H} + \mathcal{G}\alpha(\mathcal{F}, \mathcal{H}). \tag{A.2}$$

Hence, α is either symmetric or antisymmetric, depending on whether σ = 1 or σ = −1, and a derivation in all arguments. Both Poisson and metric brackets, cf. Section 2.1, are special cases of such a bilinear form.

Given the restriction functions Rx,j (u) = u^j (x), which are defined for functions u that are at least continuous, let

$$\mathscr{A}{ij}(u; x, x') := \alpha(\mathcal{R}{x,i}, \mathcal{R}_{x',j})(u),$$

and we have

$$\mathcal{A}{ij}(u; x, x') = \sigma \mathcal{A}{ji}(u; x', x).$$

We claim that a bilinear form α satisfying (A.1) and (A.2) is such that

1) it vanishes on constants, that is,

$$G(v) = a \in \mathbb{R} \text{ for all } v \implies \alpha(F, G) = 0 \text{ for all } F$$

2) it has the representation, for any F, G,

$$\alpha(\mathcal{F},\mathcal{G}) = \sum{i,j} \int{\Omega} \int_{\Omega} \frac{\delta \mathcal{F}(u)}{\delta ui}(x) \mathscr{A}{ij}(u;x,x') \frac{\delta \mathcal{F}(u)}{\delta u_j}(x') d\mu(x') d\mu(x).$$

In particular, 2) implies Eqs. (5).

Claim 1) follows from the bilinearity of the form and the Leibniz property. If G(v) = a is a constant function, for any pair of functions F and H, one has

$$a\alpha(\mathcal{F}, \mathcal{H}) = \alpha(\mathcal{F}, a\mathcal{H}) = \alpha(\mathcal{F}, \mathcal{G}\mathcal{H})$$ $$= \alpha(\mathcal{F}, \mathcal{G})\mathcal{H} + a\alpha(\mathcal{F}, \mathcal{H}),$$

from which one deduces α(F, G) = 0 as claimed.

Claim 2) requires Taylor formula: For any u, u⁰ ∈ V

$$\mathcal{F}(u) = \mathcal{F}(u0) + \theta{\mathcal{F}}(u_0, u)(u - u_0),$$

where, with v = u − u0,

$$\begin{split} \theta_{\mathcal{F}}(u_0, u)v &= \int_0^1 D\mathcal{F}\big((1-t)u_0 + tu\big)vdt \ &= \int0^1 \int\Omega \frac{\delta\mathcal{F}}{\delta u}\big((1-t)u_0 + tu\big) \cdot vd\mu(x)dt \ &= \sum_i \int0^1 \int\Omega \frac{\delta\mathcal{F}}{\delta u_i}\big((1-t)u_0 + tu\big)v_id\mu(x)dt. \end{split}$$

For a fixed u0, F(u0) and G(u0) are constants and using claim 1) we have

$$\begin{split} \alpha(\mathcal{F},\mathcal{G}) &= \alpha \big(\mathcal{F} - \mathcal{F}(u_0), \mathcal{G} - \mathcal{G}(u0)\big) \ &= \alpha \big(\theta{\mathcal{F}}(u_0,\cdot)(\cdot - u0), \theta{\mathcal{G}}(u_0,\cdot)(\cdot - u0)\big) \ &= \sum{i,j} \int_0^1 \int0^1 \int\Omega \int\Omega A{ij} dt ds d\mu(x) d\mu(x'), \end{split}$$

where we have formally exchanged the integrals and the bi-linear form and for brevity we have defined

$$A_{ij} = \alpha \left( \frac{\delta \mathcal{F}((1-t)u_0 + tu)}{\delta u_i} (x) (u(x) - u_0(x))_i, \frac{\delta \mathcal{G}((1-s)u_0 + su)}{\delta u_j} (x') (u(x') - u_0(x'))_j \right).$$

In the latter expression the first argument of α is the product of the functionals u 7→ δF((1 − t)u⁰ + tu)/δuⁱ |^x and u 7→ (u − u0)^j |x; analogously for the second argument. We can use Leibniz identity and evaluate at u = u⁰ with the result that

$$A_{ij} = \frac{\delta \mathcal{F}(u_0)}{\delta ui}(x) \alpha \left(\mathcal{R}{x,i}, \mathcal{R}_{x',j}\right) (u_0) \frac{\delta \mathcal{G}(u_0)}{\delta u_j}(x').$$

Therefore,

$$\alpha(\mathcal{F},\mathcal{G})(u0) = \sum{i,j} \int{\Omega} \int{\Omega} \frac{\delta \mathcal{F}(u_0)}{\delta ui}(x) \mathscr{A}{ij}(u_0,x,x') \frac{\delta \mathcal{G}(u_0)}{\delta u_j}(x') d\mu(x') d\mu(x),$$

and since the point u⁰ ∈ V is arbitrary this is claim 2). This argument however is purely formal: The restriction function Rx,j is defined only for functions that can be evaluated at a point x, e.g. continuous functions, thus excluding L p functions for any p. We have assumed that functional derivative exists and in exchanging the integral with the form α one needs some continuity in order to pass to the limit after approximating the integrals by finite sums. We have also freely exchanged the integration order.

Appendix B. On continuous bilinear forms on Hilbert spaces

Let H be a Hilbert space over the fields of real numbers and with scalar product (·, ·) and induced norm ∥ · ∥. If a : H × H → R is a continuous, positive bilinear form, where continuity means that there exists C > 0 such that

$$0 \le a(u, v) \le C||u||||v||,$$

for all u, v ∈ H, then one can find a bounded, symmetric, positive-definite linear operator A : H → H such that

$$a(u,v) = (u,Av),$$

for all u, v ∈ H. In order to find A one observes that, for any u fixed, ℓ(v) = a(u, v) is a bounded linear function from H → R. The Riesz representation theorem [52, Theorem 8.12] yields a unique element w ∈ H such that

$$\ell(v) = (w, v),$$

for all v ∈ H and ∥w∥ = sup{ℓ(v) : v ∈ H, ∥v∥ = 1} ≤ C∥u∥. Since w is unique, one can set w = Au, and A is a bounded linear operator on H. Then

$$a(u, v) = \ell(v) = (Au, v).$$

Positivity and symmetry follow from the positivity and symmetry of a.

Appendix C. A Lagrangian variational principle for Beltrami fields

In this appendix, we give a self-contained overview of the variational principle for Beltrami fields. This is the constant-pressure version of the variational principle for full MHD equilibria obtained by Kendall [75], and formulated in a modern language.

While the special case of linear Beltrami fields obey Woltjer's principle of least magnetic energy at constant magnetic helicity, cf. Section 2.2.3, general Beltrami fields minimize energy under a much stronger constraint.

On a bounded (not necessarily simply connected) domain Ω ⊂ R 3 , we fix a reference magnetic field B⁰ ∈ V , where V is the space of vector fields B ∈ [L 2 (Ω)]³ satisfying the conditions

div B = 0, in Ω, n · B = 0, on ∂Ω. (C.1)

For any Φ : Ω → Ω an element of the group Diff(Ω) of diffeomorphisms of the domain Ω, we define

$$B = \Phi_* B_0 = \frac{D\Phi B_0}{\det D\Phi} \circ \Phi^{-1}, \tag{C.2}$$

where DΦ is the Jacobian matrix of Φ (defined by (DΦ)ij = ∂xjΦi) and det DΦ ̸= 0 is its determinant. Then B is the push-forward of the fields B⁰ with the map Φ. A direct calculation show that

$$(\det D\Phi) \operatorname{div} B = \operatorname{div} B_0,$$

hence div B⁰ = 0 imply div B = 0. Analogously one can show that the boundary condition B⁰ · n = 0 on ∂Ω is preserved by the diffeomorphism as we have that if x = Φ(x0) and x⁰ ∈ ∂Ω, then x ∈ ∂Ω,

$$\begin{split} n(x) &= \frac{{}^tD\Phi^{-1}(x)n(x_0)}{|{}^tD\Phi^{-1}(x)n(x_0)|},\ n(x) \cdot B(x) &= \frac{B_0(x_0) \cdot n(x_0)}{\det D\Phi \cdot |{}^tD\Phi^{-1}(x)n(x_0)|}. \end{split}$$

(This can be proven by recalling that for a sufficiently regular domain, near a point x⁰ ∈ ∂Ω there is a function f such that f > 0 in Ω and f = 0 on ∂Ω; then n(x0) ∝ ∇f(x0) and this transforms like a 1-form under Φ.)

Therefore the push-forward formula (C.2) maps B⁰ ∈ V into B ∈ V .

Given B⁰ ∈ V , we define the entropy function on the group Diff(Ω) as the magnetic energy stored in B, that is

$$S(\Phi) = S(\Phi; B0) = \int{\Omega} \frac{|B|^2}{8\pi} dx, \qquad (C.3)$$

where B = Φ∗B0. The entropy depends parametrically on the initial field B0.

We can now state the variational principle for (24). For any B⁰ ∈ V fixed, if Φ is a critical point of (C.3), then B = Φ∗B⁰ ∈ V satisfies (24). More explicitly, this means that if

$$\frac{d}{d\varepsilon}\mathcal{S}(\Phi^{\varepsilon})\Big|_{\varepsilon=0} = 0,$$ (C.4)

for any curve ε 7→ Φ ^ε ∈ Diff(Ω) such that Φ^ε |ε=0 = Φ, then B = Φ∗B⁰ is a Beltrami field obtained by mapping the given field B⁰ by the action of the diffeomorphism Φ.

In order to prove the variational principle (C.4) let B^ε (x) = Φ^ε ^∗B0(x), and introduce the displacement field

$$\xi^{\varepsilon} = \partial_{\varepsilon} \Phi^{\varepsilon} \circ (\Phi^{\varepsilon})^{-1}$$ . (C.5)

The definition is equivalent to

$$\partial_{\varepsilon} \Phi^{\varepsilon}(x_0) = \xi^{\varepsilon}(x), \quad x = \Phi^{\varepsilon}(x_0).$$

Then, one obtains

$$\partial{\varepsilon} B^{\varepsilon} = \operatorname{curl}(\xi^{\varepsilon} \times B^{\varepsilon}), \quad B^{\varepsilon}|{\varepsilon=0} = B,$$ (C.6)

and we compute from (C.4),

$$\frac{d}{d\varepsilon} \mathcal{S}(\Phi^{\varepsilon}) \Big|{\varepsilon=0} = -\int{\Omega} \left[ \frac{1}{4\pi} (\operatorname{curl} B) \times B \right] \cdot \xi dx + \int_{\Omega} n \cdot [(\xi \times B) \times B] d\sigma = 0,$$

where ξ = ξ ε |ε=0. The boundary term vanishes due to the identity (ξ × B) × B = (B · ξ)B − B² ξ and the boundary condition B · n = 0, ξ · n = 0; the latter follows from the fact that Φ preserves the boundary, hence ξ|∂^Ω must be tangent to ∂Ω. Since the derivative of S(Φ^ε ) has to vanish for any curve Φ^ε ∈ Diff(Ω) and thus for every ξ, we deduce that B satisfies (24).

We remark that, since B is the push-forward of a known field B0, the field-line topology of B is the same as that of B0. Magnetic helicity is also preserved.

This is a variant of the variational principle (C.4) at the basis of the relaxation method of Chodura and Schl¨uter [18, 14, 19] and of the SIESTA code [21].

References

[1] P. J. Morrison, Bracket formulation for irreversible classical fields, Physics Letters A 100 (8) (1984) 423 – 427. doi:10.1016/0375-9601(84)90635-2.
[2] P. J. Morrison, Some observations regarding brackets and dissipation, Tech. Rep. PAM–228, University of California at Berkeley, available at arXiv:2403.14698v1 [mathph] 15 Mar 2024 (March 1984).
[3] P. J. Morrison, A paradigm for joined Hamiltonian and dissipative systems, Physica D: Nonlinear Phenomena 18 (1986) 410–419. doi:10.1016/0167-2789(86)90209-5.
[4] P. J. Morrison, Hamiltonian description of the ideal fluid, Rev. Mod. Phys. 70 (1998) 467–521. doi:10.1103/RevModPhys.70.467.
[5] J. P. Freidberg, Ideal MHD, Cambridge University Press, 2014. doi:10.1017/CBO9780511795046. [6] O. P. Bruno, P. Laurence, Existence of three-dimensional toroidal MHD equilibria with non-
constant pressure, Communications on Pure and Applied Mathematics 49 (7) (1996) 717–764. doi:10.1002/(SICI)1097-0312(199607)49:7717::AID-CPA33.0.CO;2-C. [7] F. Bauer, O. Betancourt, P. Garabedian, A Computational Method in Plasma Physics, Springer-
Verlag, 1978. doi:10.1007/978-3-642-85470-5. [8] S. P. Hirshman, J. C. Whitson, Steepest-descent moment method for three-dimensional magneto-
hydrodynamic equilibria, Physics of Fluids 26 (12) (1983) 3553–3568. doi:10.1063/1.864116. [9] L.-M. Imbert-G´erard, E. J. Paul, A. M. Wright, An Introduction to Stellarators: From Magnetic Fields to Symmetries and Optimization, Society for Industrial and Applied Mathematics, 2024.
doi:https://doi.org/10.1137/1.9781611978223. [10] H. Grad, H. Rubin, Hydromagnetic equilibria and force-free fields, Journal of Nuclear Energy (1954) 7 (3-4) (1958) 284–285.
[11] A. Reiman, H. Greenside, Calculation of three-dimensional MHD equilibria with islands and stochastic regions, Computer Physics Communications 43 (1) (1986) 157 – 167. doi:10.1016/ 0010-4655(86)90059-7. [12] H. K. Moffatt, Magnetostatic equilibria and analogous Euler flows of arbitrarily complex topol-
ogy. part 1. fundamentals, Journal of Fluid Mechanics 159 (1985) 359–378. doi:10.1017/ S0022112085003251.
[13] V. I. Arnold, B. A. Khesin, Topological Methods in Hydrodynamics, Springer-Verlag, 1998. [14] H. Moffatt, Some topological aspects of fluid dynamics, Journal of Fluid Mechanics 914 (2021) P1. doi:10.1017/jfm.2020.230.
[15] H. Alfv´en, Existence of electromagnetic-hydrodynamic waves, Nature 150 (1942) 405–406. doi: 10.1038/150405d0. [16] K. Harafuji, T. Hayashi, T. Sato, Computational study of three-dimensional magnetohydrody-
namic equilibria in toroidal helical systems, Journal of Computational Physics 81 (1) (1989) 169– 192. doi:10.1016/0021-9991(89)90069-7. [17] Y. Suzuki, N. Nakajima, K. Watanabe, et al., Development and application of HINT2 to helical
system plasmas, Nuclear Fusion 46 (11) (2006) L19. doi:10.1088/0029-5515/46/11/L01. [18] R. Chodura, A. Schlueter, A 3D code for MHD equilibrium and stability, Journal of Computational
Physics 41 (1) (1981) 68 – 88. doi:10.1016/0021-9991(81)90080-2. [19] T. Wiegelmann, T. Sakurai, Solar force-free magnetic fields, Living Reviews in Solar Physics 9 (1) (2012) 5. doi:10.1007/s41116-020-00027-4. [20] Y. Brenier, X. Duan, An integrable example of gradient flow based on optimal transport of
differential forms, Calculus of Variations and Partial Differential Equations 57 (5) (2018) 125. doi:10.1007/s00526-018-1370-6. [21] S. P. Hirshman, R. Sanchez, C. R. Cook, SIESTA: A scalable iterative equilibrium solver for
toroidal applications, Physics of Plasmas 18 (6), 062504 (2011). doi:10.1063/1.3597155. [22] B. Coquinot, P. J. Morrison, A general metriplectic framework with application to dissipative extended magnetohydrodynamics, Journal of Plasma Physics 86 (3) (2020) 835860302. doi:10.
1017/s0022377820000392. [23] P. J. Morrison, M. H. Updike, Inclusive curvature-like framework for describing dissipation: Metriplectic 4-bracket dynamics, Physical Review E 109 (2024) 045202.
[24] A. Zaidni, P. J. Morrison, S. Benjelloun, Thermodynamically consistent Cahn–Hilliard–Navier–Stokes equations using the metriplectic dynamics formalism, Physica
D 468 (2024) 134303. doi:https://doi.org/10.1016/j.physd.2024.134303. [25] N. Sato, P. J. Morrison, A collision operator for describing dissipation in noncanonical phase space, Fund. Plasma Phys. 10 (2024) 100054.
[26] A. Zaidni, P. J. Morrison, Metriplectic 4-bracket algorithm for constructing thermodynamically consistent dynamical systems, arXiv (2024) arXiv:2412.arXiv:\physics.comp-ph/2412. [27] M. W. Hirsch, S. Smale, R. L. Devaney, Differential Equations, Dynamical Systems, and an Intro-
duction to Chaos, third edition Edition, Elsevier, 2013. doi:10.1016/C2009-0-61160-0. [28] M. Materassi, Metriplectic algebra for dissipative fluids in Lagrangian formulation, Entropy 17 (3) (2015) 1329–1346. doi:10.3390/e17031329.
[29] M. Materassi, E. Tassi, Metriplectic framework for dissipative magneto-hydrodynamics, Physica D: Nonlinear Phenomena 241 (6) (2012) 729–734. doi:10.1016/j.physd.2011.12.013.
[30] M. Kraus, E. Hirvijoki, Metriplectic integrators for the Landau collision operator, Physics of
Plasmas 24 (10) (2017) 102311. doi:10.1063/1.4998610.
[31] W. Barham, P. J. Morrison, A. Zaidni, A thermodynamically consistent discretization of 1D thermal-fluid models using their metriplectic 4-bracket structure, Communications in Nonlinear Science and Numerical Simulation 145 (2025) 108683.
URL https://www.sciencedirect.com/science/article/pii/S1007570425000942
[32] G. K. Vallis, G. G. Carnevale, W. R. Young, Extremal energy properties and construction of stable
solutions of the Euler equations, J. Fluid Mech. 207 (1989) 133–152. doi:https://doi.org/10.
1017/S0022112089002533. [33] G. G. Carnevale, G. K. Vallis, Pseudo-advective relaxation to stable states of inviscid two-dimensional fluids, J. Fluid Mech. 213 (1990) 549–571. doi:https://doi.org/10.1017/
S0022112090002440. [34] G. Flierl, P. Morrison, Hamiltonian Dirac simulated annealing: Application to the calculation of vortex states, Physica D: Nonlinear Phenomena 240 (2) (2011) 212 – 232, "Nonlinear Excursions" Symposium and Volume in Physica D to honor Louis N. Howard's scientific career. doi:10.1016/
j.physd.2010.08.011. [35] Y. Chikasue, M. Furukawa, Adjustment of vorticity fields with specified values of Casimir invariants as initial condition for simulated annealing of an incompressible, ideal neutral fluid and its mhd in two dimensions, Journal of Fluid Mechanics 774 (2015) 443–459. doi:10.1017/jfm.2015.263. [36] Y. Chikasue, M. Furukawa, Simulated annealing applied to two-dimensional low-beta reduced
magnetohydrodynamics, Physics of Plasmas (1994-present) 22 (2) (2015). doi:10.1063/1.4913234. [37] M. Furukawa, P. J. Morrison, Simulated annealing for three-dimensional low-beta reduced MHD equilibria in cylindrical geometry, Plasma Physics and Controlled Fusion 59 (5) (2017) 054001. doi:10.1088/1361-6587/aa5863. [38] M. Furukawa, T. Watanabe, P. J. Morrison, et al., Calculation of large-aspect-ratio tokamak and toroidally-averaged stellarator equilibria of high-beta reduced magnetohydrodynamics via simu-
lated annealing, Physics of Plasmas 25 (8) (2018) 082506. doi:10.1063/1.5038043. [39] M. Furukawa, P. J. Morrison, Simulated annealing of reduced magnetohydrodynamic systems,
Reviews of Modern Plasma Physics 9 (1) (2025) 15. URL https://doi.org/10.1007/s41614-025-00185-8 [40] R. Brockett, Dynamical systems that sort lists, diagonalize matrices, and solve linear programming

problems, Linear Algebra and its Applications 146 (1991) 79–91. doi:10.1109/cdc.1988.194420.

[41] A. M. Bloch, R. W. Brockett, T. S. Ratiu, Completely integrable gradient flows, Communications in Mathematical Physics 147 (1) (1992) 57–74. doi:10.1007/bf02099528. [42] A. M. Bloch, P. J. Morrison, T. S. Ratiu, Flows in the normal and Kaehler metrics and triple
bracket generated metriplectic systems, in Recent Trends in Dynamical Systems, eds. A. Johann et al., Springer Proceedings in Mathematics & Statistics 35 (2013) 371–415.
[43] A. Hasegawa, Self-organization processes in continuous media, Advances in Physics 34 (1) (1985) 1–42. doi:10.1080/00018738500101721.
[44] Z. Yoshida, S. M. Mahajan, Variational principles and self-organization in two-fluid plasmas, Phys. Rev. Lett. 88 (2002) 095001. doi:10.1103/PhysRevLett.88.095001.
[45] L. Woltjer, On hydromagnetic equilibrium, Proceedings of the National Academy of Sciences 44 (9) (1958) 833–841. doi:10.1073/pnas.44.9.833. [46] J. B. Taylor, Relaxation of toroidal plasma and generation of reverse magnetic fields, Phys. Rev.
Lett. 33 (1974) 1139–1141. doi:10.1103/PhysRevLett.33.1139. [47] J. B. Taylor, Relaxation and magnetic reconnection in plasmas, Rev. Mod. Phys. 58 (1986) 741– 763. doi:10.1103/revmodphys.58.741.
[48] H. Qin, W. Liu, H. Li, et al., Woltjer-Taylor state without Taylor's conjecture: Plasma relaxation at all wavelengths, Phys. Rev. Lett. 109 (2012) 235001. doi:10.1103/PhysRevLett.109.235001.
- [49] M. F. Adams, E. Hirvijoki, M. G. Knepley, et al., Landau collision integral solver with adaptive mesh refinement on emerging architectures, SIAM Journal on Scientific Computing 39 (6) (2017) C452–C465. doi:10.1137/17M1118828.
- [50] F. Gay-Balmaz, D. D. Holm, Selective decay by Casimir dissipation in inviscid fluids, Nonlinearity 26 (2) (2013) 495. doi:10.1088/0951-7715/26/2/495. [51] F. Gay-Balmaz, D. D. Holm, A geometric theory of selective decay with applications in MHD,
- Nonlinearity 27 (8) (2014) 1747. doi:10.1088/0951-7715/27/8/1747. [52] J. K. Hunter, B. Nachtergaele, Applied Analysis, World Scientific Publishing Company, 2001.
- [53] J. E. Marsden, T. Ratiu, R. Abraham, Manifolds, Tensor Analysis, and Applications, third edition Edition, Springer, 2001. [54] P. J. Morrison, Poisson brackets for fluids and plasmas, Cont. Math. 88 (1982) 13–46. doi:

10.1063/1.33633.

NY, 2003.

[55] P. J. Morrison, S. Eliezer, Spontaneous symmetry breaking and neutral stability in the noncanonical Hamiltonian formalism, Phys. Rev. A 33 (1986) 4205–4214.
- [56] Z. Yoshida, T. Tokieda, P. J. Morrison, Rattleback: A model of how geometric singularity induces dynamic chirality, Phys. Lett. A 381 (2017) 2772–2777.
- [57] Z. Yoshida, P. J. Morrison, Deformation of Lie-Poisson algebras and chirality, J. Math. Phys. 61 (2020) 082901.
- [58] Z. Yoshida, P. J. Morrison, Hierarchical structure of noncanonical Hamiltonian systems, Physica Scripta 91 (2) (2016) 024001. doi:10.1088/0031-8949/91/2/024001. [59] M. E. Taylor, Partial Differential Equations III: Nonlinear Equations, Vol. 117 of Applied Mathematical Sciences, Springer New York, 2011. doi:10.1007/978-1-4419-7049-7.
- tational Physics 93 (1) (1991) 1 107. doi:10.1016/0021-9991(91)90074-u. [61] M. E. Taylor, Partial Differential Equations I: Basic Theory, Vol. 115 of Applied Mathematical Sciences, Springer New York, 2011. doi:10.1007/978-1-4419-7055-8.

[60] T. Takeda, S. Tokuda, Computation of MHD equilibrium of tokamak plasma, Journal of Compu-

[62] M. D. Kruskal, R. M. Kulsrud, Equilibrium of a magnetically confined plasma in a toroid, Physics of Fluids 1 (4) (1958) 265–274. doi:10.1063/1.1705884. [63] H. Grad, Some new variational properties of hydromagnetic equilibria, Physics of Fluids 7 (8) (1964) 1283–1292. doi:10.1063/1.1711373.
[64] L. L. LoDestro, L. D. Pearlstein, On the Grad-Shafranov equation as an eigenvalue problem, with implications for q solvers, Physics of Plasmas 1 (1) (1994) 90–95. doi:10.1063/1.870464. [65] A. Pataki, A. J. Cerfon, J. P. Freidberg, et al., A fast, high-order solver for the Grad-Shafranov
equation, Journal of Computational Physics 243 (2013) 28–45. doi:10.1016/j.jcp.2013.02.045. [66] Z. Yoshida, Y. Giga, Remarks on spectra of operator rot, Mathematische Zeitschrift 204 (1) (1990) 235–245. doi:10.1007/BF02570870. [67] T.-Z. Boulmezaoud, Y. Maday, T. Amari, On the linear force-free fields in bounded and unbounded

three-dimensional domains, ESAIM: M2AN 33 (2) (1999) 359–393. doi:10.1051/m2an:1999121. [68] R. L. Dewar, M. J. Hole, M. McGann, et al., Relaxed plasma equilibria and entropy-related plasma

self-organization principles, Entropy 10 (4) (2008) 621–634. doi:10.3390/e10040621. URL http://dx.doi.org/10.3390/e10040621 [69] S. R. Hudson, R. L. Dewar, G. Dennis, et al., Computation of multi-region relaxed magnetohy-
drodynamic equilibria, Physics of Plasmas 19 (11) (2012) 112502. arXiv:http://dx.doi.org/10. 1063/1.4765691, doi:10.1063/1.4765691. [70] A. M. Dixon, M. A. Berger, P. K. Browning, E. R. Priest, A generalization of the Woltjer minimum-
energy principle, Astron. Astrophys. 225 (1989) 156–166. [71] P. Laurence, M. Avellaneda, On Woltjer's variational principle for force-free fields, Journal of
- Mathematical Physics 32 (5) (1991) 1240–1253. doi:10.1063/1.529321. [72] H. K. Moffatt, The degree of knottedness of tangled vortex lines, Journal of Fluid Mechanics null (1969) 117–129. doi:10.1017/s0022112069000991. [73] V. I. Arnold, The asymptotic Hopf invariant and its applications, Springer Berlin Heidelberg,
- Berlin, Heidelberg, 2014, pp. 357–375. doi:10.1007/978-3-642-31031-7\_32. [74] T. Vogel, On the asymptotic linking number, Proc. of the American Mathematical Society 131 (7) (2003) 2289–2297. doi:10.1090/s0002-9939-02-06792-8.
- [75] P. C. Kendall, The variational formulation of the magneto-hydrostatic equations, Astrophys. J. 131 (1960) 681. [76] T. Amari, C. Boulbe, T. Z. Boulmezaoud, Computing Beltrami fields, SIAM Journal on Scientific
- Computing 31 (5) (2009) 3217–3254. doi:10.1137/070700942. [77] D. Malhotra, A. Cerfon, L.-M. Imbert-G´erard, et al., Taylor states in stellarators: A fast highorder boundary integral solver, Journal of Computational Physics 397 (2019) 108791. doi:10.
- 1016/j.jcp.2019.06.067. [78] D. Henry, Geometric Theory of Semilinear Parabolic Equations, Springer, Berlin, 1981.
- [79] R. Temam, Infinite-Dimensional Dynamical Systems in Mechanics and Physics, Springer-Verlag New York, 1998. [80] S. Wiggins, Introduction to Applied Nonlinear Dynamical Systems and Chaos, Springer New York,
- [81] V. Moretti, Analytical Mechanics, Springer Nature Switzerland, 2023. [82] J. La Salle, S. Lefschetz, Stability by Liapunov's Direct Method with Applications, Academic
- Press, 1961.
- [83] G. P. Beretta, A theorem on Lyapunov stability for dynamical systems and a conjecture on a
property of entropy, J. Math. Phys. 27 (1) (1986) 305–308. doi:10.1063/1.527390.
[84] B. Polyak, Gradient methods for the minimisation of functionals, USSR Computational Mathematics and Mathematical Physics 3 (4) (1963) 864–878. doi:10.1016/0041-5553(63)90382-3.
[85] S. Lojasiewicz, Sur les trajectoires du gradient d'une fonction analytique, in: Seminari di Geometria 1982–1983, Universita di Bologna, Istituto di Geometria, Dipartamento di Matematica, 1984, p. 115–117.
[86] J. Palis, W. Melo, Geometric Theory of Dynamical Systems, Springer New York, 1982. doi: 10.1007/978-1-4612-5703-5.
[87] H. Karimi, J. Nutini, M. Schmidt, Linear convergence of gradient and proximal-gradient methods under the Polyak- Lojasiewicz condition, in: Machine Learning and Knowledge Discovery in Databases, Springer International Publishing, Cham, 2016, pp. 795–811.
[88] J. Rauch, Hyperbolic Partial Differential Equations and Geometric Optics, Graduate studies in mathematics, American Mathematical Society, 2012. [89] C. Bressan, Metriplectic relaxation for calculating equilibria: theory and structure-preserving
discretization, Ph.D. thesis, Technische Universit¨at M¨unchen (2023). URL https://mediatum.ub.tum.de/1686142 [90] C. Bressan, M. Kraus, P. J. Morrison, et al., Relaxation to magnetohydrodynamics equilibria via collision brackets, Journal of Physics: Conference Series 1125 (2018) 012002. doi:10.1088/
1742-6596/1125/1/012002. [91] L. C. Evans, Partial Differential Equations, Graduate studies in mathematics, American Mathematical Society, Providence (R.I.), 1998.
[92] P. Monk, Finite Element Methods for Maxwell's equations, Oxford Un. Press, 2003.
[93] M. S. Alnaes, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg, C. N. Richardson, J. Ring, M. E. Rognes, G. N. Wells, The FEniCS project version 1.5, Archive of Numerical Software 3 (2015). doi:10.11588/ans.2015.100.20553.
[94] A. Logg, K.-A. Mardal, G. N. Wells, et al., Automated Solution of Differential Equations by the Finite Element Method, Springer, 2012. doi:10.1007/978-3-642-23099-8.
[95] M. S. Alnaes, A. Logg, K. B. Ølgaard, M. E. Rognes, G. N. Wells, Unified form language: A domain-specific language for weak formulations of partial differential equations, ACM Transactions on Mathematical Software 40 (2014). doi:10.1145/2566630.
[96] A. Logg, G. N. Wells, DOLFIN: Automated finite element computing, ACM Transactions on Mathematical Software 37 (2010). doi:10.1145/1731022.1731030.
[97] A. Logg, G. N. Wells, J. Hake, DOLFIN: a C++/Python finite element library, in: A. Logg, K.-A. Mardal, G. N. Wells (Eds.), Automated Solution of Differential Equations by the Finite Element Method, Vol. 84 of Lecture Notes in Computational Science and Engineering, Springer, 2012, Ch. 10.
[98] F. Herrnegger, in: Proceedings of 5th European Conference on Controlled Fusion and Plasma Physics, Grenoble, European Physical Society, 1972, p. 26.
[99] E. K. Maschke, Exact solutions of the MHD equilibrium equation for a toroidal plasma, Plasma Physics 15 (6) (1973) 535. doi:10.1088/0032-1028/15/6/006.
[100] P. J. Mc Carthy, Analytical solutions to the Grad-Shafranov equation for tokamak equilibrium with dissimilar source functions, Phys. Plasmas 6 (9) (1999) 3554–3560. doi:10.1063/1.873630.
[101] E. Zoni, Y. G¨u¸cl¨u, Solving hyperbolic-elliptic problems on singular mapped disk-like domains with the method of characteristics and spline finite elements, Journal of Computational Physics 398 (2019) 108889.
- URL https://www.sciencedirect.com/science/article/pii/S002199911930587X
[102] W. H. Yang, P. A. Sturrock, S. K. Antiochos, Force-free magnetic fields: The magneto-frictional method, Astrophys. J. 309 (1986) 383–391.
[103] J. A. Klimchuk, P. A. Sturrock, Three-dimensional force-free magnetic fields and flare energy buildup (1992).
[104] G. Valori, B. Kliem, M. Fuhrmann, Magnetofrictional extrapolations of Low and Lou's force-free equilibria, Solar Physics 245 (2) (2007) 263–285. doi:10.1007/s11207-007-9046-y.
[105] G. Valori, B. Kliem, T. T¨or¨ok, V. S. Titov, Testing magnetofrictional extrapolation with the Titov-D´emoulin model of solar active regions, A&A 519 (Sep. 2010). doi:10.1051/0004-6361/201014416. [106] K. Hu, Y.-J. Lee, J. Xu, Helicity-conservative finite element discretization for incompressible MHD
systems, Journal of Computational Physics 436 (2021) 110284. [107] M. He, P. E. Farrell, K. Hu, B. D. Andrews, Topology-preserving discretization for the magneto-
frictional equations arising in the Parker conjecture (2025). doi:https://doi.org/10.48550/ arXiv.2501.11654.
[108] S. Candelaresi, D. Pontin, G. Hornig, Mimetic methods for Lagrangian relaxation of magnetic fields, SIAM Journal on Scientific Computing 36 (6) (2014) B952–B968. doi:10.1137/140967404.
[109] Y. Guo, C. Xia, R. Keppens, G. Valori, Magneto-frictional modeling of coronal nonlinear forcefree fields. I. Testing with analytic solutions, The Astrophysical Journal 828 (2) (2016) 82. doi: 10.3847/0004-637X/828/2/82.
[110] M. Padilla, O. Gross, F. Kn¨oppel, A. Chern, U. Pinkall, P. Schr¨oder, Filament based plasma, ACM Trans. Graph. 41 (4) (jul 2022). doi:10.1145/3528223.3530102.
[111] O. Gross, U. Pinkall, P. Schr¨oder, Plasma knots, Physics Letters A 480 (2023) 128986.
[112] A. R. Yeates, On the limitations of magneto-frictional relaxation, Geophysical & Astrophysical Fluid Dynamics 116 (4) (2022) 305–320. doi:10.1080/03091929.2021.2021197. URL https://doi.org/10.1080/03091929.2021.2021197

bressanmetriplecticrelaxationequilibria2025.md 275 KB

História Raw

status: converted

Metriplectic relaxation to equilibria

Abstract

Contents

1. Introduction

1.1. Overview of equilibrium calculations

1.2. Overview of metriplectic relaxation

2. Metriplectic dynamics and variational principles for equilibria of fluids and plasmas

2.1. Metriplectic dynamics

2.2. Examples of equilibrium problems

2.2.1. Reduced Euler equations

2.2.2. Axisymmetric MHD equilibria

2.2.3. Beltrami fields

3. On the relaxation of metriplectic systems

3.1. Finite-dimensional systems: Lyapunov stability

3.2. Finite-dimensional systems: the Polyak– Lojasiewicz inequality

3.3. Infinite-dimensional systems: tentative generalization of the PL inequality

4. Two examples: metric double brackets and projectors

4.1. Metric double brackets

4.2. Projector-based metric bracket

5. Collision-like metric brackets

5.1. General construction of collision-like brackets

5.2. Collision-like brackets based on div–grad operators

5.3. Collision-like brackets based on curl–curl operators

5.4. Application to the reduced Euler equations

5.5. Application to Grad-Shafranov equilibria

6. Diffusion-like metric brackets

6.1. General construction of the brackets

6.2. Diffusion-like brackets based on div–grad operators

6.3. Diffusion-like brackets based on curl–curl operators

6.4. Application to nonlinear Beltrami fields

7. Conclusions

Appendix A. Bilinear forms and Leibniz identity

Appendix B. On continuous bilinear forms on Hilbert spaces

Appendix C. A Lagrangian variational principle for Beltrami fields

References

bressanmetriplecticrelaxationequilibria2025.md 275 KB História Raw

status: converted

Metriplectic relaxation to equilibria

Abstract

Contents

1. Introduction

1.1. Overview of equilibrium calculations

1.2. Overview of metriplectic relaxation

2. Metriplectic dynamics and variational principles for equilibria of fluids and plasmas

2.1. Metriplectic dynamics

2.2. Examples of equilibrium problems

2.2.1. Reduced Euler equations

2.2.2. Axisymmetric MHD equilibria

2.2.3. Beltrami fields

3. On the relaxation of metriplectic systems

3.1. Finite-dimensional systems: Lyapunov stability

3.2. Finite-dimensional systems: the Polyak– Lojasiewicz inequality

3.3. Infinite-dimensional systems: tentative generalization of the PL inequality

4. Two examples: metric double brackets and projectors

4.1. Metric double brackets

4.2. Projector-based metric bracket

5. Collision-like metric brackets

5.1. General construction of collision-like brackets

5.2. Collision-like brackets based on div–grad operators

5.3. Collision-like brackets based on curl–curl operators

5.4. Application to the reduced Euler equations

5.5. Application to Grad-Shafranov equilibria

6. Diffusion-like metric brackets

6.1. General construction of the brackets

6.2. Diffusion-like brackets based on div–grad operators

6.3. Diffusion-like brackets based on curl–curl operators

6.4. Application to nonlinear Beltrami fields

7. Conclusions

Appendix A. Bilinear forms and Leibniz identity

Appendix B. On continuous bilinear forms on Hilbert spaces

Appendix C. A Lagrangian variational principle for Beltrami fields

References

bressanmetriplecticrelaxationequilibria2025.md 275 KB

História Raw