# Hamiltonian hydrodynamic reductions of one-dimensional Vlasov equations

Rayan Oufar, Cristel Chandre

► To cite this version:

Rayan Oufar, Cristel Chandre. Hamiltonian hydrodynamic reductions of one-dimensional Vlasov equations. 2026. hal-05465405

HAL Id: hal-05465405

https://hal.science/hal-05465405v1

Preprint submitted on 19 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Hamiltonian hydrodynamic reductions of one-dimensional Vlasov equations

Rayan Oufar$^{a)}$ and Cristel Chandre$^{b)}$

CNRS, Aix Marseille Univ, I2M, Marseille, France

We investigate Hamiltonian fluid reductions of the one-dimensional Vlasov-Poisson
equation. Our approach utilizes the hydrodynamic Poisson bracket framework, which
allows us to systematically identify fundamental normal variables derived from the
analysis of the Casimir invariants of the resulting Poisson bracket. This framework
is then applied to analyze several well-established Hamiltonian closures of the one-
dimensional Vlasov equation, including the multi-delta distribution and the waterbag
models. Our key finding is that all of these seemingly distinct closures consistently
lead to the formulation of a unified form of parametric closures: When expressed in
terms of the identified normal variables, the parameterization across all these closures
is revealed to be polynomial and of the same degree. All these parametric closures
are uniquely generated from one of the moments, called $\mu_2$, a cubic polynomial in the
normal variables. This result establishes a structural connection between these differ-
ent physical models, offering a path toward a more unified and simplified description
of the one-dimensional Vlasov-Poisson dynamics through its reduced hydrodynamic
forms with an arbitrary number of fluid variables.

<footnote>a) Electronic mail: rayan.oufar@univ-amu.fr</footnote>

<footnote>b) Electronic mail: cristel.chandre@cnrs.fr</footnote>
# I. INTRODUCTION

The one-dimensional Vlasov-Poisson system describes the dynamics of a collisionless plasma through its phase-space distribution function $f(t, x, p)$, which gives the density of particles at time $t$, position $x$, and momentum $p$. It is written as

$$ \frac{\partial f}{\partial t} = -v \frac{\partial f}{\partial x} - qE(t, x) \frac{\partial f}{\partial p}, \qquad (1) $$

where $v = p/m$ denotes the particle velocity and $q$ is the electric charge of the particle species. The self-consistent electrostatic field $E(t, x)$ follows from Poisson's equation

$$ \frac{\partial E}{\partial x} = q \left( \int f(t, x, p) \, dp - n_0 \right), \qquad (2) $$

with $n_0$ denoting the density of a uniform neutralizing background that ensures overall charge neutrality. In what follows, we assume $q = m = 1$. Although the kinetic description (1)-(2) is complete, it evolves a function over the full phase space $(x, p)$, which can be computationally and analytically demanding, especially in higher dimensions, and provide more information than needed. A standard approach toward model reduction is to introduce fluid moments of the distribution $f$

$$ P_n(t, x) = \int p^n f(t, x, p) \, dp, \qquad (3) $$

for $n = 0, 1, \dots$, under the assumption that $f$ decays sufficiently rapidly as $|p| \to \infty$. These moments correspond to macroscopic plasma quantities: $P_0$ is the particle density, $P_1$ the momentum density, $P_2$ is related to pressure and temperature, and $P_3$ to the heat flux. By evolving a reduced set of variables in configuration space, one obtains a more tractable model with clearer and more intuitive physical interpretations.

Each dynamical equation for $P_n$, e.g., mass conservation, momentum balance, energy balance, involves a higher-order moment $P_{n+1}$. The term in the Vlasov equation responsible for this shift in moments is the streaming term $-v \partial f / \partial x$ (see also Ref. [1]). Thus, the hierarchy does not naturally close and, without additional modeling assumptions, remains untractable as an infinite ladder of partial differential equations.

To obtain a closed fluid model, one must prescribe a closure relation expressing higher-order moments in terms of lower-order ones [2–13]. Common closure strategies include:

* Isothermal closure: Pressure proportional to density, modeling a fixed-temperature regime.
* Adiabatic closure: A polytropic relation such as $P_2 \propto P_0^\gamma$, appropriate for collisionless compressions.

* Cold plasma closure: Neglecting higher-order moments ($P_{n \ge 2} = 0$), suitable for mono-energetic beams.

* Kinetic-based closures: Derived via asymptotic methods such as Chapman–Enskog expansions for weakly collisional regimes.

* Closures based on an ansatz on the shape of the distribution function $f$, such as the waterbags [14, 15], the multi-delta distributions [16, 17], or a generalized Maxwellian distribution [9, 18].

* Hamiltonian closures [19–21]: Ensuring the Hamiltonian structure for the truncated system, retaining the Hamiltonian character of the Vlasov equation.

The chosen closure controls the level of kinetic information retained and must balance physical accuracy, mathematical well-posedness, and computational feasibility. When selected appropriately, moment closures capture important macroscopic behavior—including wave propagation and pressure effects—while reducing the complexity inherent to full kinetic simulations. Usually these closures amount to finding suitable functions, aka, equations of state, of the form

$$P_{n \ge N} = P_{n \ge N}(P_0, P_1, \dots, P_{N-1}), \quad (4)$$

in order to truncate the infinite ladder to the first $N$ fluid moments. Such equations of state are often selected for their practicality and their modeling of specific conditions of the plasma. Given the explicit form of the bracket, $N-3$ functions $P_{n \ge N}$ of $N$ variables should be specified.

In Refs. [19–21], a strategy to find relevant closures was developed based on the Hamiltonian structure of the parent kinetic equation. By solving directly the set of nonlinear partial differential equations for the closure functions resulting from the Jacobi identity, Hamiltonian closure functions were identified and studied for $N \le 4$. In particular, in Ref. [21], it was shown that these closure functions for $N = 4$ can take very intricate forms, and sometimes cannot be formulated in an explicit way. Instead, it was shown that a parametric expression of the same cumbersome closure can express it in a quite simple form. These
parametric forms for the moments are

$$P_n = P_n(\nu_1, \dots, \nu_N),$$

for $n \ge 0$, and with as many carefully selected parameters $\nu_i$ as the number of field variables. In particular, for $n = 0, \dots, N-1$, the first $N$ relations can be locally inverted $\nu_i = \nu_i(P_0, \dots, P_{N-1})$ and inserted in the remaining $N-3$ equations $P_{n \ge N}$ to define (locally) an equation of state of the form (4). Of course, this inversion should always be possible in principle (to ensure the Jacobi identity), but often cannot be written explicitly. In the following, we demonstrate that defining closures parametrically significantly simplifies the reduction: We show that the closure is generated by a single function $\mu_2$ of $N-2$ variables, instead of the initial $N-3$ functions $P_{n \ge N}$ of $N$ variables. Furthermore, by examining all known Hamiltonian closures of the one-dimensional Vlasov-Poisson equation, we show that they each admit a polynomial parameterization in $\nu_k$ uniquely determined by a single cubic polynomial $\mu_2$ in $N-2$ normal variables. This underlying structure reveals a deep connection between seemingly distinct closures, providing a unified framework for describing the one-dimensional Vlasov-Poisson dynamics through reduced hydrodynamic models with an arbitrary number of fluid variables.

In Sec. II, we recall some known elements on the Hamiltonian formulation of the Vlasov equation and its fluid reductions. In Sec. III, we define what hydrodynamic Poisson brackets for fluid and plasmas are, and provide some general properties. In Sec. IV, we investigate all known Hamiltonian closures in the light of the properties of the hydrodynamic brackets.

## II. PRELIMINARIES

The Vlasov equation possesses a Hamiltonian structure [22, 23], in the sense that it can be derived from a Poisson bracket $\{\cdot, \cdot\}$ and a Hamiltonian $H$ as $\dot{f} = \{f, H\}$.

We remind that a Poisson bracket on an algebra $\mathcal{A}$ is a bilinear operator $\{\cdot, \cdot\} : \mathcal{A} \times \mathcal{A} \to \mathcal{A}$ satisfying the following properties for all $F, G, H \in \mathcal{A}$:

- antisymmetry: $\{F, G\} = -\{G, F\}$,

- Jacobi identity: $\{F, \{G, H\}\} + \{G, \{H, F\}\} + \{H, \{F, G\}\} = 0$,

- Leibniz rule: $\{F, GH\} = \{F, G\}H + G\{F, H\}$.
The Hamiltonian associated with the Vlasov-Poisson equation is defined as

$$H[f] = \iint f \frac{p^2}{2} dx dp + \frac{1}{2} \int E^2 dx,$$

where $E[f]$ is given by Eq. (2) and the Poisson bracket

$$\{F, G\} = \iint f \left( \frac{\partial}{\partial x} \frac{\delta F}{\delta f} \frac{\partial}{\partial p} \frac{\delta G}{\delta f} - \frac{\partial}{\partial p} \frac{\delta F}{\delta f} \frac{\partial}{\partial x} \frac{\delta G}{\delta f} \right) dx dp. \quad (5)$$

Using the moments (3), the Poisson bracket reduces to the Kupershmidt-Manin bracket [24, 25]

$$\{F, G\} = \int m P_{n+m-1} (\partial_x F_n G_m - F_m \partial_x G_n) dx, \quad (6)$$

with implicit summation over repeated indices $n$ and $m$, and where $F_n$ denotes the functional derivative of the observable $F$ with respect to $P_n$. Using an integration by parts, this Poisson bracket is rewritten as

$$\{F, G\} = \int (\partial_x F_n \alpha_{nm}(\mathbf{P}) G_m + F_n \beta_{nm}(\mathbf{P}, \partial_x \mathbf{P}) G_m) \, dx, \quad (7)$$

where

$$\alpha_{nm}(\mathbf{P}) = (n + m)P_{n+m-1}, \quad (8a)$$

$$\beta_{nm}(\mathbf{P}, \partial_x \mathbf{P}) = n \partial_x P_{n+m-1}. \quad (8b)$$

We notice that $\alpha$ is symmetric and $\partial_x \alpha = \beta + \beta^t$, to ensure that the resulting bracket is antisymmetric. If the change of variables $\mathbf{P} \mapsto \mathbf{Q} = \mathbf{Q}(\mathbf{P})$, i.e., a change of variables which only depend on the values of the original field variables (and not, e.g., of its derivatives) is performed, the form (7) of the Poisson bracket remains unchanged with $\alpha$ and $\beta$ replaced by $\bar{\alpha}$ and $\bar{\beta}$ given by

$$\bar{\alpha}_{kl}(\mathbf{Q}) = \frac{\partial Q_k}{\partial P_n} \alpha_{nm} \frac{\partial Q_l}{\partial P_m}, \quad (9a)$$

$$\bar{\beta}_{kl}(\mathbf{Q}, \partial_x \mathbf{Q}) = \partial_x \frac{\partial Q_k}{\partial P_n} \alpha_{nm} \frac{\partial Q_l}{\partial P_m} + \frac{\partial Q_k}{\partial P_n} \beta_{nm} \frac{\partial Q_l}{\partial P_m}. \quad (9b)$$

In particular, we notice that $\bar{\alpha}$ is still a function of the values of the field variables (and not of their derivatives), and that if $\beta$ is linear in $\partial_x \mathbf{P}$, then $\bar{\beta}$ is a function of $\mathbf{Q}$ and $\partial_x \mathbf{Q}$, linear in $\partial_x \mathbf{Q}$. We also notice that the matrix $\bar{\alpha}$ is still symmetric and that it satisfies $\partial_x \bar{\alpha} = \bar{\beta} + \bar{\beta}^t$. The invariance of the Poisson brackets of the form (7) led to the definition of hydrodynamic
brackets [26–29]. Given this invariance, we explore the Hamiltonian closures within this family of hydrodynamic Poisson brackets.

Truncating Poisson brackets of the form (7), i.e., restricting the sum over $n$ and $m$ from 0 to $N-1$, and imposing a closure of the form (4), still preserve the shape of the bracket, in the sense that $\alpha$ remains only a function of the values of the field variables, and $\beta$ remains a function of the values and first derivatives of these field variables, linear in the first derivatives. However, the Jacobi identity identity is most likely broken for a generic equation of state of the form (4). As a result, the dynamics generated by truncated brackets of the form (7) is no longer Hamiltonian. The central question in finding Hamiltonian closures if for which Eq. (4), the bracket (7) is a Poisson bracket?

### III. HYDRODYNAMIC POISSON BRACKETS FOR FLUID AND PLASMAS

Here we consider the algebra $\mathcal{A}$ of functionals $F$ of the $N$ fields $\mathbf{u}(x) = (u_1(x), u_2(x), \dots, u_N(x))$. In what follows, we assume that the field variables $u_n$, and all their derivatives vanish at infinity. We also assume implicit summation over repeated indices.

#### A. Definition and properties of hydrodynamic Poisson brackets

Following the above discussion and in particular Eq. (7), we define a hydrodynamic bracket [26–28] as

$$ \{F, G\} = \int (\partial_x F_n \alpha_{nm}(\mathbf{u}) G_m + F_n \beta_{nm}(\mathbf{u}, \partial_x \mathbf{u}) G_m) \, dx, \quad (10) $$

where $F_n = \delta F / \delta u_n$ is the functional derivatives of $F$ with respect to the field variable $u_n$, and the implicit summation over $n$ and $m$ is from 1 to $N$. For these brackets, $\alpha(\mathbf{u})$ is symmetric, i.e., $\alpha = \alpha^t$, and only depends on the values of the functions $u_n(x)$ and $\beta_{nm}(\mathbf{u}, \partial_x \mathbf{u}) = \beta_{nmk}(\mathbf{u}) \partial_x u_k$, for all $n, m = 1, \dots, N$, i.e., $\beta$ is linear in $\partial_x \mathbf{u}$. These brackets are antisymmetric if and only if

$$ \partial_x \alpha_{nm} = \beta_{nm} + \beta_{mn}, $$

and they obviously satisfy the Leibniz rule. The truncated bracket (7) in the moments $(P_0, \dots, P_{N-1})$ is a hydrodynamic bracket for all $N$.
Remark: Using an integration by parts, a hydrodynamic bracket can be mapped to:

$$ \{F, G\} = \int a_{nm}(\mathbf{u}) (F_n \partial_x G_m - G_n \partial_x F_m) \, dx, \qquad (11) $$

where $\alpha_{nm} = a_{nm} + a_{mn}$ and $\beta_{nm} = \partial_x a_{mn}$, as it is for the bracket (6). The inconvenience of working with the brackets of the form (11) is the presence of a gauge invariance: By adding a constant antisymmetric matrix to $a$, the expression of the bracket, and in particular, $\alpha$ and $\beta$, does not change.

Under some non-degeneracy condition, the Poisson brackets (10) can be locally mapped to a canonical form where $\alpha$ is constant and $\beta = 0$. This result can be viewed as the “equivalent” of the Lie-Darboux theorem for finite dimensional Hamiltonian systems. It is formalized using the following theorem [26]:

**Theorem III.1.** (see Ref. [26]) For the bracket (10) with non-degenerate $\alpha$ to satisfy the Jacobi identity, a necessary and sufficient condition is that there exists a locally invertible change of field variables $\mathbf{u} \mapsto \mathbf{v} = \mathbf{v}(\mathbf{u})$ such that the bracket expressed in the new coordinates $\mathbf{v}$ is of hydrodynamic type with a constant metric $g$ and no $\beta$-part, i.e., in the field variables $\mathbf{v}$, bracket (10) becomes

$$ \{\overline{F}, \overline{G}\} = \int \partial_x \overline{F}_n g_{nm} \overline{G}_m \, dx, \qquad (12) $$

where $\overline{F}[\mathbf{v}] = F[\mathbf{u}]$ and $\overline{F}_n = \delta\overline{F}/\delta\nu_n$.

**Lemma III.1.** The elements of the bracket (10), namely $\alpha_{nm}$ and $\beta_{nm}$ are given by

$$ \alpha_{nm}(\mathbf{v}) = \frac{\partial u_n}{\partial \nu_k} g_{kl} \frac{\partial u_m}{\partial \nu_l}, \qquad (13) $$

$$ \beta_{nm}(\mathbf{v}, \partial_x \mathbf{v}) = \partial_x \frac{\partial u_n}{\partial \nu_k} g_{kl} \frac{\partial u_m}{\partial \nu_l}, \qquad (14) $$

*Proof.* This is a consequence of the identities (9) since $\beta = 0$ in the variables $\mathbf{v}$. $\square$

A Casimir invariant is an observable that lies in the center of $\mathcal{A}$, i.e, $\{C, F\} = 0$ for all $F \in \mathcal{A}$. As a consequence of Theorem III.1, a non-degenerate Poisson bracket of hydrodynamic type (10) has as many independent Casimir invariants as the number of field variables, and they are of the form

$$ C_n[\mathbf{u}] = \int \nu_n(\mathbf{u}) \, dx, $$

for $n = 1, \dots, N$.
**Corollary III.2.** Non-degenerate Poisson brackets of hydrodynamic type are classified by the signature $(l, q)$ of the metric $g$, where $q$ is the number of negative eigenvalues and $l$ is the number of positive ones. This signature is invariant by a change of coordinates of the type $\nu \mapsto \eta = \eta(\nu)$.

*Proof.* By Sylvester's law of inertia, two congruent metrics have the same signature, and since under a coordinate change $\nu \mapsto \eta = \eta(\nu)$, we have $\bar{g} = J^t g J$ (see Eq. (9a), where $J$ is the Jacobian of the change of variables, $\bar{g}$ and $g$ have the same signature. $\square$

By extension, we define the signature of a Poisson bracket of hydrodynamic type (10) as the signature of the associated metric $g$.

*Remark:* Using an appropriate rescaling and reshuffling of the field variables, a Poisson bracket of hydrodynamic type with signature $(l, q)$, can be cast into the form (12) with

$$g = \begin{pmatrix} \mathbb{I}_l & 0 \\ 0 & -\mathbb{I}_q \end{pmatrix},$$

where $\mathbb{I}_l$ is the $l \times l$ identity matrix.

**Definition III.1 (Kinetic hydrodynamic bracket).** A hydrodynamic bracket of the form (10) comes from a kinetic model when there exists a local change of variables $u \mapsto P(u)$ that maps it to a truncated Kupershmidt-Manin bracket (7) (i.e., truncated in the sense that the summations are restricted to $n, m = 0, \dots, N-1$).

**Corollary III.3.** A non-degenerate hydrodynamic Poisson bracket (10) coming from a kinetic model has a signature $(l, q)$ with $l, q \ge 1$.

*Proof.* By focusing on the subalgebra of observables $F[P_0, P_1]$ for the bracket (7), the matrix $\alpha$ becomes

$$\alpha = \begin{pmatrix} 0 & P_0 \\ P_0 & 2P_1 \end{pmatrix},$$

which has a negative determinant. Henceforth, it has a positive and a negative eigenvalue. $\square$

This corollary leads to the definition of partially decoupled brackets, which has the purpose of partially isolating $(P_0, P_1)$ from the rest.
## B. Partial decoupling of hydrodynamic Poisson brackets

**Definition III.2** (Partially decoupled bracket). A Poisson bracket of hydrodynamic type (10) is said to be partially decoupled if it has the following form:

$$ \{F, G\} = \int \left( \partial_x F_u G_\rho - F_\rho \partial_x G_u + \partial_x w_k \left( F_u \frac{G_k}{\rho} - \frac{F_k}{\rho} G_u \right) \right. \\ \left. + \partial_x \frac{F_k}{\rho} \alpha_{kl}(\mathbf{w}) \frac{G_l}{\rho} + \frac{F_k}{\rho} \beta_{kl}(\mathbf{w}, \partial_x \mathbf{w}) \frac{G_l}{\rho} \right) dx, \quad (15) $$

where the implicit sums over k, l are from 1 to N − 2 and $F_k = \delta F / \delta w_k$.

Here the field variables are $(\rho, u, w_1, \dots, w_{N-2})$. We notice that partially decoupled brackets are invariant under any change of field variables $\mathbf{w} \mapsto \mathbf{z}(\mathbf{w})$.

**Remark:** The set of functionals $F[\rho, u]$ and the set of functionals $F[\rho, \mathbf{w}]$ form subalgebras of the bracket (15). This remark leads to the definition of the microscopic part of a partially decoupled hydrodynamic bracket with $N − 2$ field variables:

$$ \{F, G\}_m = \int \left( \partial_x F_k \alpha_{kl}(\mathbf{w}) G_l + F_k \beta_{kl}(\mathbf{w}, \partial_x \mathbf{w}) G_l \right) dx, \quad (16) $$

i.e., where the macroscopic degree of freedom represented by the fluid density and fluid velocity $(\rho, u)$ has been eliminated.

**Proposition III.1.** The bracket (15) with $N$ field variables is a Poisson bracket if and only if the microscopic bracket (16) with $N − 2$ field variables is a Poisson bracket.

**Proof.** Suppose first that the bracket (15) is Poisson. Since the set of functionals $F[\rho, \mathbf{w}]$ is a subalgebra, it follows that the induced bracket (16) is a Poisson bracket.

Conversely, suppose that bracket (16) is a Poisson bracket. By Theorem III.1 there exists a local change of variables $\mathbf{w} \mapsto \mathbf{v}(\mathbf{w})$ such that $\{F, G\}_m = \int g_{kl} \partial_x F_k G_l dx$, where $F_k = \delta F / \delta v_k$.

In the variables $(\rho, \psi, \mathbf{v})$, the bracket $\{F, G\} = \int (\partial_x F_\psi G_\rho - F_\rho \partial_x G_\psi + \partial_x F_k g_{kl} G_l) dx$ is of Poisson type since the coefficients do not depend on the field variables. Under the following changes of variables $(\rho, \psi, \mathbf{v}) \mapsto (\rho, u, \eta)$ where $u = \mathbf{v} \cdot g^{-1} \mathbf{v} / (2\rho) + \psi$ and $\eta_k = \nu_k / \rho$ and then $\mathbf{v} \mapsto \mathbf{w}(\mathbf{v})$, this bracket is mapped to bracket (10). □

**Proposition III.2.** Any hydrodynamic bracket (10) with signature $(l, q)$ where $l, q \ge 1$ can be cast into a partially decoupled hydrodynamic bracket. In particular, any hydrodynamic
bracket coming from a kinetic equation, such as Vlasov, can be cast into a partially decoupled hydrodynamic bracket.

*Proof.* We consider a hydrodynamic Poisson bracket (10) with signature $(l, q)$ with $l, q \ge 1$. Using Corollary III.3, it is rewritten in the $\nu$ variables as

$$ \{F, G\} = \int (\partial_x F_1 G_1 - \partial_x F_2 G_2 + \varepsilon_k \partial_x F_k G_k) \, dx, $$

where $\varepsilon_k = \pm 1$ for $k = 3, \dots, N$. Under the following change of variables, $\rho = (\nu_1 + \nu_2)/\sqrt{2}$, $u = (\nu_1 - \nu_2)/\sqrt{2} + (\sum_{k \ge 3} \varepsilon_k \nu_k^2)/(2\rho)$ and $\bar{\nu}_k = \nu_k/\rho$, it takes the partially decoupled form (15). $\square$

It follows that a non-degenerate hydrodynamic bracket of kinetic origin can be cast in the following form

$$ \{F, G\} = \int \left( \partial_x F_u G_\rho - F_\rho \partial_x G_u + \partial_x \nu_k \left( F_u \frac{G_k}{\rho} - \frac{F_k}{\rho} G_u \right) + \partial_x \frac{F_k}{\rho} g_{kl} \frac{G_l}{\rho} \right) dx, \quad (17) $$

where $F_k = \delta F / \delta \nu_k$ and $g$ is a constant symmetric matrix.

**Lemma III.2.** A non-degenerate partially decoupled Poisson bracket of hydrodynamic type (10) has $N-2$ Casimir invariants of the form

$$ C[\rho, \nu] = \int \rho \nu \, dx, $$

in addition to the two Casimir invariants

$$ C[\rho] = \int \rho \, dx, \quad (18) $$

$$ C[\rho, u, \nu] = \int \left( u - \frac{\rho}{2} \nu \cdot g^{-1} \nu \right) dx, \quad (19) $$

where $\nu = \nu(\mathbf{w})$ are the coordinates which flatten the metric of bracket (16) into the metric $g$ as in Eq. (17).

*Proof.* A Casimir invariant $C = \int f(\mathbf{w}) \, dx$ of the bracket (16) yields a Casimir invariant $C = \int \rho f(\mathbf{w}) \, dx$ of the bracket (15). The conservation of the Casimir invariant $C[\rho, u, \nu]$ can be verified by a direct and straightforward calculation of $\{C, \rho\}$, $\{C, u\}$ and $\{C, \nu\}$. It then follows from the Leibniz rule that $\{C, F\} = 0$ for all $F$. $\square$
### C. Centered moments

To write a hydrodynamic bracket coming from a kinetic model into a partially decoupled hydrodynamic bracket, a strategy used in Refs. [19–21, 30] is to move to the moments centered around the velocity, i.e.,

$$S_n(t, x) = \frac{1}{\rho^{n+1}} \int (v-u)^n f(t, x, v) \, dv, \qquad (20)$$

where the fluid velocity is given by

$$u = \frac{1}{\rho} \int v f(t, x, v) \, dv.$$

In this way, we notice that $S_0 = 1$ and $S_1 = 0$. The centered moments $\mathbf{S}$ are obtained from the moments $\mathbf{P}$ as

$$S_n = \frac{1}{\rho^{n+1}} \sum_{k=0}^{n} \binom{n}{k} (-u)^{n-k} P_k,$$

and inverted as

$$P_n = \sum_{k=0}^{n} \binom{n}{k} \rho^{k+1} u^{n-k} S_k.$$

In particular, we notice that $S_n$ depends on $(P_0, \dots, P_n)$, and vice versa, $P_n$ depends on $(\rho, u, S_2, \dots, S_n)$. Using the change of variables $(P_0, P_1, \dots, P_{N-1}) \mapsto (\rho, u, S_2, \dots, S_{N-1})$, the bracket becomes a partially decoupled hydrodynamic bracket of the form (15) with

$$\alpha_{nm}(\mathbf{S}) = (n+m)S_{n+m-1} - m(n+1)S_n S_{m-1} - n(m+1)S_m S_{n-1},$$

$$\beta_{nm}(\mathbf{S}, \partial_x \mathbf{S}) = n \partial_x S_{n+m-1} - n S_{n-1} \partial_x S_m - n(m+1) S_m \partial_x S_{n-1} - n m S_{m-1} \partial_x S_n,$$

where the closure is simplified to

$$S_{n \ge N} = S_{n \ge N}(S_2, \dots, S_{N-1}).$$

In particular, we notice that these closures do not depend on $\rho$ and $u$.

Following Lemma III.2, another interesting set of centered variables is obtained by centering the moments not with respect to $u$ but with respect to $\psi(x) \equiv u - \rho\mu_1$ such that $C = \int \psi(x) \, dx$ given by Eq. (19) is a Casimir invariant of the bracket. Imposing that $C$ is a Casimir invariant defines the variable $\mu_1$. For the bracket expressed in the $\nu$-variables, it takes the form

$$\mu_1(x) = \frac{1}{2} \nu \cdot g^{-1} \nu, \qquad (21)$$
which is a homogeneous quadratic polynomial in the $N-2$ variables $\nu_k$. We also define $\mu_n$ for $n \ge 2$ as

$$ \mu_n(t, x) = \frac{1}{\rho^{n+1}} \int (v - \psi)^n f(t, x, v) \, dv. \qquad (22) $$

These variables are linked to the ones used in Ref. [31]. The $\nu$-variables are linked to the $S$ variables in the following way:

$$ S_n = \sum_{k=0}^{n} \binom{n}{k} (-\mu_1)^{n-k} \mu_k. $$

Next we perform the change of variables $(\rho, u, S_2, \dots, S_{N-1}) \mapsto (\rho, u, \mu_1, \dots, \mu_{N-2})$. In Appendix A, we show that the bracket expressed in these variables remains of the type (15) with

$$ \alpha_{nm}(\boldsymbol{\mu}) = (n+m)\mu_{n+m-1} - m\mu_{m-1}\gamma_n - n\mu_{n-1}\gamma_m, \qquad (23a) $$

$$ \beta_{nm}(\boldsymbol{\mu}, \partial_x \boldsymbol{\mu}) = n \partial_x \mu_{n+m-1} - n \partial_x \mu_{n-1} \gamma_m - m \mu_{m-1} \partial_x \gamma_n, \qquad (23b) $$

where the explicit expression of $\gamma_n$ depends on the specific closure. In the variables $\nu_k$, it becomes

$$ \gamma_n = (n+1)\mu_n - \nu_k \frac{\partial \mu_n}{\partial \nu_k}. \qquad (24) $$

If the parametric expression of $\mu_n$ is a homogeneous polynomial of degree $n+1$, we have that $\gamma_n = 0$. As a consequence, the equations for the elements $\alpha$ and $\beta$ in $\mu$ variables are identical as the equations in the $P$ variables (see Eq. (8)). This is the case for some of the cases examined below, such as the multi-delta closure and Burby's closure [31]. However, there exist cases where $\gamma_n \ne 0$, such as the waterbags.

In what follows, we use the notation for the field variables and the closure functions as follows:

* $P_n$: moments of $f$ in velocity given by Eq. (3),

* $S_n$: centered moments of $f$ around the fluid velocity $u$, given by Eq. (20),

* $\mu_n$: centered moments of $f$ around $\psi = u - \rho\mu_1$, given by Eq. (22),

* Greek letters such as $\nu_k$, $\eta_k$ or $\xi_k$ for the field variables which flatten $\alpha$ into a constant metric $g$ (and consequently no $\beta$-part).
Also, in the rest of the article, $N$ always designates the number of field variables of the main fluid model (i.e., including $\rho$ and fluid velocity).

**Proposition III.3.** The closure, i.e., $\mu_n = \mu_n(\nu_1, \dots, \nu_{N-2})$ for $n \ge 1$, is generated by a single moment $\mu_2$.

*Proof.* We assume that the hydrodynamic bracket (15) with elements $\alpha$ and $\beta$ given by Eq. (23) is a Poisson bracket. The change of variables from $\mu$ to $\nu$ ensures that

$$ \frac{\partial \mu_n}{\partial \nu_k} g_{kl} \frac{\partial \mu_m}{\partial \nu_l} = (n+m)\mu_{n+m-1} - m\mu_{m-1}\gamma_n - n\mu_{n-1}\gamma_m, \\ \partial_x \frac{\partial \mu_n}{\partial \nu_k} g_{kl} \frac{\partial \mu_m}{\partial \nu_l} = n \partial_x \mu_{n+m-1} - n \partial_x \mu_{n-1} \gamma_m - m \mu_{m-1} \partial_x \gamma_n. $$

For $m=2$, the first identity becomes

$$ \mu_{n+1} = \frac{1}{n+2} \left[ \frac{\partial \mu_n}{\partial \nu_k} g_{kl} \frac{\partial \mu_2}{\partial \nu_l} + 2\mu_1 \gamma_n + n \mu_{n-1} \gamma_2 \right], $$

where $\gamma_n$ is given by Eq. (24) and $\mu_1$ given by Eq. (21). Since $\gamma_n$ only depends on $\mu_n$ and its first derivatives, this recurrence relation ensures that all $\mu_{n \ge 3}$ are uniquely determined by $\mu_2$. $\square$

# IV. CLASSIFICATION OF KNOWN HAMILTONIAN FLUID CLOSURES

In this section, we investigate all the known Hamiltonian closures in the light of the results listed in Sec. III.

## A. The multi-delta distribution

We consider the Vlasov distribution (as considered in Refs. [16, 17, 32–38]) given by $N = 2M$ fields $a_k(x)$ and $v_k(x)$ for $k = 1, \dots, M$:

$$ f(t, x, v) = \sum_{k=1}^{M} a_k(t, x) \delta(v - v_k(t, x)), \qquad (25) $$

where $M$ is the number of Dirac delta functions. From the Poisson bracket (5), the Poisson bracket associated with the distribution (25) is given by:

$$ \{F, G\} = \int (\partial_x F_{a_k} G_{v_k} - \partial_x G_{a_k} F_{v_k}) \, dx. \qquad (26) $$
We notice that the above bracket is readily of hydrodynamic form (10) and its metric is flat, i.e., $\beta = 0$ and $\alpha$ is given by

$$ \alpha = \begin{pmatrix} 0 & \mathbb{I}_M \\ \mathbb{I}_M & 0 \end{pmatrix}. $$

Consequently, the signature of the bracket (26) is $(M, M)$. The moments $\mathbf{P}$ are given by

$$ P_n = \sum_{k=1}^{M} a_k v_k^n. $$

We notice that the moment $P_n$ is a homogeneous polynomial of degree $n + 1$ in $N = 2M$ variables. Next, we proceed with the partial decoupling by considering the fluid density and fluid velocity. We consider the following variables:

$$ \begin{aligned} \rho &= \sum_{l=1}^{M} a_l, \\ u &= \frac{1}{\rho} \sum_{l=1}^{M} a_l v_l, \\ \xi_k &= \frac{a_k}{\rho}, \\ \eta_k &= \frac{v_k - v_1}{\rho}, \end{aligned} $$

for $k = 2, \dots, M$. This change of variables is invertible, and its inverse is given by

$$ \begin{aligned} a_1 &= \rho \left( 1 - \sum_{l=2}^{M} \xi_l \right), \\ v_1 &= u - \rho \sum_{l=2}^{M} \xi_l \eta_l, \\ a_k &= \rho \xi_k, \\ v_k &= u - \rho \sum_{l=2}^{M} \xi_l \eta_l + \rho \eta_k, \end{aligned} $$

for $k = 2, \dots, N-1$. In the normal variables $\boldsymbol{\nu} = (\boldsymbol{\xi}, \boldsymbol{\eta})$, the bracket becomes of the form (17) with

$$ g = \begin{pmatrix} 0 & \mathbb{I}_{M-1} \\ \mathbb{I}_{M-1} & 0 \end{pmatrix}. $$

As a consequence, the signature of the microscopic bracket (16) is $(M - 1, M - 1)$.
The centered fluid moments $S$ are expressed as,

$$ 
\begin{aligned}
S_n &= \frac{1}{\rho^{n+1}} \sum_{i=1}^{M} a_i (v_i - u)^n, \\
&= \left(1 - \sum_{k=2}^{M} \xi_k\right) \left(-\sum_{k=2}^{M} \eta_k \xi_k\right)^n + \sum_{k=2}^{M} \xi_k \left(\eta_k - \sum_{l=2}^{M} \xi_l \eta_l\right)^n.
\end{aligned}
$$

Remark: $S_n$ is a polynomial of degree $2n$ in the $2M-2$ normal variables $(\xi, \eta)$.

Next we compute the $\mu$-variables by centering the moments around $\psi = u - \rho\mu_1$ where $\mu_1$ is given by Eq. (21) as $\mu_1 \equiv \sum_{k=2}^{M} \xi_k \eta_k$, i.e.,

$$ 
\begin{aligned}
\mu_n &= \frac{1}{\rho^{n+1}} \sum_{k=1}^{M} a_k (v_k - u + \rho\mu_1)^n, \\
&= \sum_{k=2}^{M} \xi_k \eta_k^n,
\end{aligned}
$$

since $v_1 - u + \rho\mu_1 = 0$. Given that $\mu_n$ is a homogeneous polynomial of degree $n+1$, the functions $\gamma_n$ given by Eq. (24) vanish. Consequently, under the changes of variables $(\rho, u, \xi_2, \dots, \xi_M, \eta_2, \dots, \eta_M) \to (\rho, u, \mu_1, \dots, \mu_{N-2})$, one gets the partially decoupled hydrodynamic bracket of the form (15) with

$$ 
\begin{aligned}
\alpha_{nm}(\boldsymbol{\mu}) &= (n+m)\mu_{n+m-1}, \\
\beta_{nm}(\boldsymbol{\mu}, \partial_x \boldsymbol{\mu}) &= m\partial_x \mu_{n+m-1},
\end{aligned}
$$

for $n, m = 1, \dots, N-2$ (see also Appendix A).

We notice that the expression of the moments $\mu_n$ in the $2M-2$ variables $(\xi, \eta)$ is the same as the moments $P_n$ in the $2M$ variables $(\mathbf{a}, \mathbf{v})$.

The cubic polynomial generating the closure is

$$ \mu_2 = \sum_{k=2}^{M} \xi_k \eta_k^2. $$

## B. The waterbag distribution

We consider the Vlasov distributions given by the waterbag ansatz (first introduced in Refs. [14, 15, 39, 40]), i.e., a piecewise constant function in the velocity domain,

$$ f(t, x, v) = \sum_{n=1}^{N} a_n \Theta(v - v_n(t, x)), $$
where $a_n$ are constant, $\Theta$ is the Heaviside distribution, and $v_n(t, x)$ are the contour velocities of the $N-1$ bags of height $\sigma_k = \sum_{n=1}^k a_n$ for $k = 1, \dots, N-1$. To ensure a compact support, we impose the condition $\sigma_N = \sum_{n=1}^N a_n = 0$.

In the algebra of functionals of $v = (v_1, \dots, v_N)$, the Poisson bracket is given by [41–44]

$$ \{F, G\}_{\text{wb}} = - \int \sum_{n=1}^{N} \frac{1}{a_n} \partial_x F_n G_n \, dx, $$

where $F_n = \delta F / \delta v_n$. We notice that this bracket is a hydrodynamic bracket (10) with $\beta = 0$ and $\alpha_{nm} = -\delta_{nm} / a_n$. As a consequence, $N$ independent Casimir invariants in the variables $v_n$ are given by

$$ C_n = \int v_n \, dx. \qquad (27) $$

We observe that the moments (3), given by

$$ P_n = - \frac{1}{n+1} \sum_{n=1}^{N} a_n v_n^{n+1}, $$

are homogeneous polynomials of degree $n+1$ in the $N$ variables $v_n$.

Next, we cast this Poisson bracket into a partially decoupled form (15) by introducing the normal variables :

$$ \begin{aligned} \rho &= - \sum_{n=1}^{N} a_n v_n, \\ u &= - \frac{1}{2\rho} \sum_{n=1}^{N} a_n v_n^2, \\ \nu_k &= \frac{1}{\rho} \sum_{l=1}^{k} \sigma_l (v_{l+1} - v_l), \end{aligned} $$

for $k = 1, \dots, N-2$. We notice that $\nu_k$ is the cumulative normalized particle density up to the $k$-th bag and, using the definition of $\nu_k$, we have $\nu_{N-1} = 1$ and $\nu_0 = 0$. This change of variables is invertible and its inverse is given by

$$ v_1 = u + \frac{\rho}{2} \sum_{k=2}^{N} a_k \left( \sum_{l=1}^{k-1} \frac{\nu_l - \nu_{l-1}}{\sigma_l} \right)^2, \qquad (28a) $$

$$ v_k = v_1 + \rho \sum_{l=1}^{k-1} \frac{\nu_l - \nu_{l-1}}{\sigma_l}, \qquad (28b) $$
for $k = 2, \dots, N$. By performing the change of variables $(v_1, \dots, v_N) \mapsto (\rho, u, \nu_1, \dots, \nu_{N-2})$, the Poisson bracket becomes a partially decoupled hydrodynamic bracket (17) with

$$g_{nm} = -\delta_{nm} \frac{\sigma_n \sigma_{n+1}}{a_{n+1}}.$$

Since $\sigma_n$ corresponds to the height of the n-th waterbag, the sign of $g_{nm}$ is the opposite of the sign of $a_{n+1}$. Consequently, the signature of $\{\cdot, \cdot\}_{\text{wb}}$ is $(N-k, k)$ where $k \in \{1, \dots, N-1\}$ is the number of positive a's. This means that the signature of the bracket is directly linked to the predefined number of increasing vs decreasing bags in the distribution. Since the heights of the bags are constant, this signature is a characteristic of the shape of the chosen initial distribution.

The centered fluids moments $S$ are given by [30]

$$S_n = -\frac{1}{(n+1)\rho^{n+1}} \sum_{k=1}^{N} a_k (v_k - u)^{n+1}, \\ = -\frac{1}{(n+1)} \sum_{k=1}^{N} a_k \eta_k^{n+1}(\nu_1, \dots, \nu_{N-2}),$$

where, using Eq. (28), the expressions of $\eta_k$ are given by

$$\eta_1 = \frac{1}{2} \sum_{k=2}^{N} a_k \left( \sum_{l=1}^{k-1} \frac{\nu_l - \nu_{l-1}}{\sigma_l} \right)^2, \\ \eta_k = \eta_1 + \sum_{l=1}^{k-1} \frac{\nu_l - \nu_{l-1}}{\sigma_l},$$

for $k = 2, \dots, N$, where $\nu_0 = 0$ and $\nu_{N-1} = 1$.

*Remark: $S_n$ is a polynomial of degree $2n$ in the normal variables $\nu$. We notice that for this closure, the polynomial is not homogeneous since*

$$S_n(\boldsymbol{\nu} = 0) = \frac{1 + (-1)^n}{(n+1)2^{n+1}a_N^n}.$$

Instead of centering the moments around $u$, we center the moments around $\psi$ defined by $\psi = u - \rho\mu_1$. The first moment $\mu_1$ is given by Eq. (21):

$$\mu_1 = \frac{1}{2} \sum_{k=1}^{N-2} \frac{\nu_k^2}{\lambda_k},$$

where $\lambda_k = -\sigma_k \sigma_{k+1}/a_{k+1}$. The explicit expression of $\psi$ as a function of $v_n$ becomes

$$\psi = \frac{v_N}{2} + \frac{\rho}{2a_N} = -\frac{1}{2a_N} \sum_{n=1}^{N-1} a_n v_n,$$
where, as expected, it can be shown to be a Casimir invariant given Eq. (27). The moments $\mu_n$ for $n \ge 2$ as functions of the variables $\nu_k$ are given by

$$ \begin{aligned} \mu_n &= -\frac{1}{n+1} \frac{1}{\rho^{n+1}} \sum_{k=1}^{N} a_k (\nu_k - \psi)^{n+1}, \\ &= \frac{(-1)^n}{n+1} \left( \sum_{k=1}^{N-1} a_k \left( \frac{1}{2a_N} + \sum_{l=k}^{N-1} \frac{\nu_l - \nu_{l-1}}{\sigma_l} \right)^{n+1} + \frac{1}{2^{n+1}a_N^n} \right). \end{aligned} $$

From this expression, we notice that $\mu_n$ is a polynomial of degree $n+1$ in the $N-2$ variables $\nu_k$. In the variables $\mu$ the brackets takes the form (15) with

$$ \alpha_{nm} = (n+m)\mu_{n+m-1} - m\mu_{m-1}\Lambda^n - n\mu_{n-1}\Lambda^m + 2\Lambda nm\mu_{n-1}\mu_{m-1}, $$

$$ \beta_{nm} = n\partial_x\mu_{n+m-1} - n\partial_x\mu_{n-1}\Lambda^m + 2\Lambda nm\partial_x\mu_{n-1}\mu_{m-1}, $$

where $\Lambda = -1/(2a_N)$. The inhomogeneous character of the polynomial $\mu_n$ is expressed by a non-zero $\gamma_n$ where $\gamma_n$ is given by Eq. (24)

$$ \gamma_n = \Lambda^n - n\Lambda\mu_{n-1}, $$

whereas this quantity is zero for all other known closures.

Although the coefficients $\alpha$ and $\beta$ depend only on the constant $\Lambda$, the realization of the moment variables $\mu$ as functions of the normal variables $\nu$ depends explicitly on the signs of the waterbag heights. Consequently, two waterbag distributions with identical $\Lambda$ but different signatures give rise to non-isomorphic Poisson structures, despite having formally identical reduced brackets in the $\mu$-variables.

The cubic polynomial generating the closure is

$$ \mu_2 = \frac{1}{3} \left( \sum_{k=1}^{N-1} a_k \left( \frac{1}{2a_N} + \sum_{l=k}^{N-1} \frac{\nu_l - \nu_{l-1}}{\sigma_l} \right)^3 + \frac{1}{8a_N^2} \right). $$

### C. Burby's Hamiltonian closure (see Ref. [31])

In this section we consider the Hamiltonian closure obtained in Ref. [31]. In the one-dimensional case, this closure corresponds to considering a Vlasov distribution of the form

$$ f(t, x, v) = \sum_{k=0}^{N-2} \frac{(-1)^k}{k!} \mu_k(t, x) \delta^{(k)}(v - \psi(t, x)), $$
where $\delta^{(k)}$ is the k-th derivative of the Dirac delta function defined as

$$ \int \delta^{(k)}(v)\phi(v) \, dv = (-1)^k \phi^{(k)}(0), $$

for an arbitrary function $\phi \in C^\infty(\mathbb{R})$. The moments $P_n$ are given by

$$ P_n = \sum_{k=0}^{n} \binom{n}{k} \mu_k \psi^{n-k}. $$

The centered moments $S_n$ are given by

$$ S_n = \sum_{k=0}^{n} \binom{n}{k} \bar{\mu}_k (-\bar{\mu}_1)^{n-k}, $$

where $\bar{\mu}_k = \mu_k / \rho^{k+1}$. As expected, we have $S_0 = \bar{\mu}_0 = 1$ and $S_1 = 0$. Therefore the centered moments $S_n$ are only functions of $\bar{\mu}_k$ for $k = 1, \dots, N-2$.

*Remark:* If $\bar{\mu}_n$ is a polynomial of degree $n+1$ for $n \ge 1$, then $S_n$ is a polynomial of degree $2n$. Therefore we notice the advantage of working with the variables $\bar{\mu}_n$ rather than the centered variables $S_n$, to obtain closure functions of the lowest possible degree.

We note that the first fluid moments are given by

$$ \rho = \mu_0, \\ u = \frac{\mu_1}{\rho} + \psi, $$

In the variables $(\rho, \psi, \mu_1, \dots, \mu_{N-2})$, the Poisson bracket becomes

$$ \{F, G\} = \int (\partial_x F\psi G_\rho - F_\rho \partial_x G_\psi + \partial_x F_k \alpha_{kl}(\boldsymbol{\mu}) G_l + F_k \beta_{kl}(\boldsymbol{\mu}, \partial_x \boldsymbol{\mu}) G_l) \, dx, \quad (29) $$

where the implicit summations over $k$ and $l$ run from 1 to $N-2$, and

$$ \alpha_{kl}(\boldsymbol{\mu}) = (k+l)\mu_{k+l-1}, \quad (30a) $$

$$ \beta_{kl}(\boldsymbol{\mu}, \partial_x \boldsymbol{\mu}) = k \partial_x \mu_{k+l-1}, \quad (30b) $$

and $\mu_k = 0$ for $k \ge N-1$. In other terms, the bracket is of hydrodynamic type and the elements $\alpha$ and $\beta$ are upper anti-triangular. Contrary to the examples of the multi-delta distribution and the waterbags, the expression of the Poisson bracket has an explicit dependence on the field variables, so computing the signature of the bracket in the original variables is not trivial.
In order to proceed, we consider the following change of variables: $(\rho, \psi, \mu_1, \dots, \mu_{N-2}) \mapsto (\rho, u, \bar{\mu}_1, \dots, \bar{\mu}_{N-2})$ defined by

$$u = \frac{\mu_1}{\rho} + \psi, \\ \bar{\mu}_k = \frac{\mu_k}{\rho^{k+1}},$$

for $k = 1, \dots, N-2$. This leads to the partial decoupling of the hydrodynamic bracket in the form given by Eq. (15), where the expressions of the elements $\bar{\alpha}$ and $\bar{\beta}$ satisfy Eq. (30) with $\bar{\mu}$ instead of $\mu$. In this way, the Poisson bracket is already in the reduced variables $\bar{\mu}$.

# 1. Expression of the Casimir invariants

Since the bracket (29) satisfies the Jacobi identity, it has $N$ Casimir invariants. Two obvious Casimir invariants are

$$C[\rho] = \int \rho \, dx, \\ C[\psi] = \int \psi \, dx.$$

The other $N-2$ Casimir invariants related to the specific shape of the elements $\alpha$ and $\beta$ are less trivial, but are of the form

$$C[\rho, \bar{\mu}] = \int \rho \nu(\bar{\mu}) \, dx,$$

as mentioned in Lemma III.2. The goal of this section to obtain their explicit expressions as functions of $\bar{\mu}$. In what follows, we drop the bars above $\mu$, and the scaling with $\rho$ is implicitly assumed.

Finding the expressions of the Casimir invariants amounts to finding the change of variables $\mu \mapsto \nu = \nu(\mu)$, or alternatively to determine the parametrization $\mu = \mu(\nu)$ which maps the bracket (15) into the bracket (17) and invert it. The equations to be solved are obtained from Lemma III.1 and Eq. (30) for $m \equiv N-2$:

$$\frac{\partial \mu_k}{\partial \nu_i} g_{ij} \frac{\partial \mu_l}{\partial \nu_j} = (k+l) \mu_{k+l-1}, \qquad (31a)$$

$$\partial_x \frac{\partial \mu_k}{\partial \nu_i} g_{ij} \frac{\partial \mu_l}{\partial \nu_j} = k \partial_x \mu_{k+l-1}, \qquad (31b)$$
for $k, l = 1, \dots, m$, and where $\mu_k = 0$ for $k \ge m + 1$. The implicit summations over $i$ and $j$ are taken from 1 to $m$. First, we assume the matrix $g$ to be the following $m \times m$ antidiagonal matrix:

$$ g = \begin{pmatrix} 0 & \cdots & 0 & 1 \\ \vdots & \ddots & 1 & 0 \\ 0 & 1 & \cdots & \vdots \\ 1 & 0 & \cdots & 0 \end{pmatrix}. \qquad (32) $$

The identities (31) to be solved simplify to:

$$ \frac{\partial \mu_k}{\partial \nu_j} \frac{\partial \mu_l}{\partial \nu_{m+1-j}} = (k+l) \mu_{k+l-1}, $$

$$ \partial_x \frac{\partial \mu_k}{\partial \nu_j} \frac{\partial \mu_l}{\partial \nu_{m+1-j}} = k \partial_x \mu_{k+l-1}. $$

The goal is to find the solution of these coupled partial differential equations for an arbitrary $m$.

**Lemma IV.1.** The functions $\mu_n^{(m)}$ of the $m + 1 - n$ last variables $\nu$, i.e.,

$$ \mu_n^{(m)} = \mu_n^{(m)}(\nu_n, \dots, \nu_m), $$

for $n = 1, \dots, m$, satisfy the following relations

$$ \frac{\partial \mu_k^{(m)}}{\partial \nu_j} \frac{\partial \mu_l^{(m)}}{\partial \nu_{m+1-j}} = 0, $$

$$ \partial_x \frac{\partial \mu_k^{(m)}}{\partial \nu_j} \frac{\partial \mu_l^{(m)}}{\partial \nu_{m+1-j}} = 0, $$

for $k + l > m + 1$.

*Proof.* The derivative $\partial \mu_k^{(m)} / \partial \nu_j$ is nonzero if $j \ge k$, and $\partial \mu_l^{(m)} / \partial \nu_{m+1-j}$ is nonzero if $m + 1 - j \ge l$. The condition $m + 1 - j \ge l$ together with $k + l > m + 1$ leads to $j < k$. So there is no possible index $j$ such that both partial derivatives are non-zero. $\square$

We look for solutions to Eq. (31) of this form. In what follows, we consider functions $\mu_n$ for different values of $m$, i.e., we define the family $\mathcal{M}$ of centered moments $\mu_n^{(m)}$ for $n = 1, \dots, m$ and $m = 1, \dots, N - 2$. Here $m$ is referred to as the level of the moment(s), and therefore $\mu_n^{(m)}$ is the $n$-th centered moment of level $m$.
**Lemma IV.2.** Let $\mathcal{M}$ be a family of centered moments $\mu_n^{(m)}(\nu)$ for $n = 1, \dots, m$ and $m = 1, \dots, M$, such that $\mu_n^{(m)}(0) = 0$. The two proposals are equivalent:

(i) $\mathcal{M}$ satisfies the following differential identities:

$$ \frac{\partial}{\partial \nu_k} [\mu_n^{(m)}(\nu_n, \dots, \nu_m)] = n \mu_{n-1}^{(k-1)}(\nu_{m-k+n}, \dots, \nu_m), \quad (33) $$

for all $k = n, \dots, m$, with $\mu_0^{(l)}(x_0, \dots, x_l) = x_0$ for all levels $l$.

(ii) $\mathcal{M}$ satisfies the following recursive relations:

$$ \mu_n^{(m)}(\nu_n, \dots, \nu_m) = \begin{cases} \frac{\nu_m^{m+1}}{m+1} & \text{if } n = m, \\ \sum_{k=0}^{n} \binom{n}{k} \nu_m^{n-k} \mu_k^{(m-n-1)}(\nu_{k+n}, \dots, \nu_{m-1}) & \text{if } n < m. \end{cases} $$

We notice that the identity (i) states that the partial derivatives of the $m$-th level moments $\mu_n^{(m)}$ are given by elements of $\mathcal{M}$ of lower levels, $\mu^{(k-1)}$ for $k = n, \dots, m$. In a similar way, we notice that the identity (ii) states that the elements of level $m$, namely, $\mu_n^{(m)}$ for $n < m$, are polynomials of degree $n$ in the last variable $\nu_m$ whose coefficients depend on the other elements of $\mathcal{M}$ of lower levels, $\mu^{(m-n-1)}$.

*Proof.* See Appendix B 1. $\square$

The induction identities (ii) are used to compute the analytic expressions of the moments $\mu_n^{(n)}$. The Python code and the expression of the first moments for levels $m \le 5$ are given in Appendix C.

**Corollary IV.1.** The centered moments $\mu_n^{(m)}$ of $\mathcal{M}$ are homogeneous polynomials of degree $n + 1$ for $n = 0, \dots, m$.

*Proof.* This follows by induction from the identity (ii) in Lemma IV.2. If $\mu_k^{(m-n-1)}$ is a homogeneous polynomial of degree $k + 1$ for all $k = 0, \dots, n$ and $n = 0, \dots, m - 1$, then $\nu_m^{n-k} \mu_k^{(m-n-1)}$ and, consequently, $\mu_n^{(m)}$ are homogeneous polynomials of degree $n + 1$. $\square$

**Proposition IV.1.** Let the centered moments $\mu_n^{(m)}$ satisfying the hypotheses of Lemma IV.2. Then the family of moments $\mu_n^{(m)}$ is a solution of Eq. (31).

*Proof.* See Appendix B 2. $\square$
**Proposition IV.2.** The moments $\mu_n^{(m)}$, for $n = 0, \dots, m$, are given by

$$ \mu_n^{(m)}(\nu_n, \dots, \nu_m) = \frac{1}{n+1} \sum_{\substack{n \le i_1, i_2, \dots, i_{n+1} \le m \\ i_1 + \dots + i_{n+1} = n(m+1)}} \nu_{i_1} \nu_{i_2} \dots \nu_{i_{n+1}}. \quad (34) $$

*Proof.* We prove Eq. (34) by induction. The first two cases $n = 0$ and $n = 1$ are straightforward. We assume Eq. (34) holds for $n - 1$, i.e.,

$$ \mu_{n-1}^{(m)}(\nu_{n-1}, \dots, \nu_m) = \frac{1}{n} \sum_{\substack{n-1 \le i_1, i_2, \dots, i_n \le m \\ i_1 + \dots + i_n = (n-1)(m+1)}} \nu_{i_1} \dots \nu_{i_n}. $$

We define the candidate formula for $n$ as

$$ F_n^{(m)}(\nu_n, \dots, \nu_m) = \frac{1}{n+1} \sum_{\substack{n \le i_1, i_2, \dots, i_{n+1} \le m \\ i_1 + \dots + i_{n+1} = n(m+1)}} \nu_{i_1} \dots \nu_{i_{n+1}}. $$

For all $k \in \{n, \dots, m\}$, we compute the partial derivative of $F_n^{(m)}$ with respect to $\nu_k$:

$$ \frac{\partial F_n^{(m)}}{\partial \nu_k} = \sum_{i_1 + \dots + i_n = n(m+1) - k} \nu_{i_1} \dots \nu_{i_n}. $$

By reindexing

$$ j_{\ell} = i_{\ell} - (m - k + 1), \quad \text{for } \ell = 1, \dots, n, $$

we rewrite the sum as

$$ \frac{\partial F_n^{(m)}}{\partial \nu_k} = n \mu_{n-1}^{(k-1)}(\nu_{m-k+n}, \dots, \nu_m), $$

which matches the identity (33). Since both $F_n^{(m)}$ and $\mu_n^{(m)}$ vanish at the origin, they must coincide: $\mu_n^{(m)}(\nu_n, \dots, \nu_m) = F_n^{(m)}(\nu_n, \dots, \nu_m)$. $\square$

**Proposition IV.3.** The parameterizations $\mu_n^{(m)} = \mu_n^{(m)}(\nu)$ are locally invertible for all $m$, and are in fact globally invertible when $m$ is even.

*Proof.* We proceed by inverting the anti-triangular system starting from $\mu_m^{(m)}$. We notice that the relation (ii) in Lemma IV.2 can be rewritten as

$$ \mu_n^{(m)} = \nu_n \nu_m^n + \chi_n^{(m)}(\nu_{n+1}, \dots, \nu_m), $$

where

$$ \chi_n^{(m)}(\nu_{n+1}, \dots, \nu_m) = \sum_{k=1}^{n} \binom{n}{k} \nu_m^{n-k} \mu_k^{(m-n-1)}(\nu_{k+n}, \dots, \nu_{m-1}), $$
First, we consider $m$ even. For $m \ge 1$, we have:

$$ \nu_m(\mu_m) = \operatorname{sgn}(\mu_m^{(m)})((m+1)|\mu_m^{(m)}|)^{1/(m+1)}. $$

We now move upwards in the anti-triangular system using the induction:

$$ \nu_n = (\operatorname{sgn}(\mu_n^{(m)}))^{n} ((m+1)|\mu_n^{(m)}|)^{-n/(m+1)} (\mu_n^{(m)} - \bar{\chi}_n^{(m)}(\mu_{n+1}^{(m)}, \dots, \mu_m^{(m)})), $$

where we have inserted the inversion relation from $m$ to $n-1$ such that

$$ \chi_n^{(m)}(\nu_{n+1}, \dots, \nu_m) = \bar{\chi}_n^{(m)}(\mu_{n+1}^{(m)}, \dots, \mu_m^{(m)}). $$

Using the same argument the relations (ii) in Lemma IV.2 can also be locally inverted for $m$ odd and $\mu_n^{(m)} > 0$. When $m$ is odd and $\mu_n^{(m)} < 0$, we consider the following variables:

$$ \bar{\mu}_n^{(m)}(\nu_n, \dots, \nu_m) = (-1)^n \mu_n^{(m)}(\nu_n, \dots, \nu_m). $$

If the polynomials $\bar{\mu}_n^{(m)}$ are solutions of Eq. (31) with $g$ given by Eq. (32), then the polynomials $\mu_n^{(m)}$ are solution of Eq. (31) with $-g$. We consider the same expressions for the inverse as for $\mu_n^{(m)}$ and replace $\mu_n^{(m)}$ by $(-1)^n \bar{\mu}_n^{(m)}$. $\square$

**Proposition IV.4.** The signature of the bracket (29) is $(\lfloor N/2 \rfloor, \lfloor N/2 \rfloor)$ or $(\lceil N/2 \rceil, \lfloor N/2 \rfloor)$.

*Proof.* This follows from a direct computation of the eigenvalues of $g$ given by Eq. (32). If $\mu_m < 0$, the metric is $-g$, so the eigenvalues are opposite to the ones for $\mu_m > 0$, where $m = N-2$. $\square$

For $N$ even, we notice that the signature is the same as for the multi-delta closure of Sec. IV A.

The cubic polynomial generating the closure is

$$ \mu_2(\nu_2, \dots, \nu_{N-2}) = \frac{1}{3} \sum_{\substack{2 \le i,j,k \le N-2 \\ i+j+k=2(N-1)}} \nu_i \nu_j \nu_k. $$

### D. Four-field Hamiltonian closure of Ref. [21]

We consider a fluid model obtained by using the first $N=4$ fluid moments derived in Ref. [21]. We recall that the Hamiltonian closure was not obtained from a closed form for the
distribution function, but by solving directly the Jacobi identity. We consider the variables $(\rho, u, S_2, S_3)$. The bracket has the form (15) with the $2 \times 2$ matrices $\alpha$ and $\beta$ given by

$$ \alpha = \begin{pmatrix} 4S_3 & 5S_4 - 9S_2^2 \\ 5S_4 - 9S_2^2 & 6S_5 - 24S_2S_3 \end{pmatrix}, $$

and

$$ \beta = \partial_x \begin{pmatrix} 2S_3 & 2S_4 - 3S_2^2 \\ 3S_4 - 6S_2^2 & 3S_5 - 12S_2S_3 \end{pmatrix}. $$

In Ref. [20], it was shown that the Jacobi identity leads to the following constraints on the closure functions $S_4$ and $S_5$:

$$ \begin{aligned} 6S_5 &= 12S_2S_3 + 4S_3 \frac{\partial S_4}{\partial S_2} + (5S_4 - 9S_2^2) \frac{\partial S_4}{\partial S_3}, \\ \frac{\partial S_5}{\partial S_2} &= 4S_3 + \frac{\partial S_4}{\partial S_3} \left( \frac{\partial S_4}{\partial S_2} - 3S_2 \right), \\ \frac{\partial S_5}{\partial S_3} &= \frac{\partial S_4}{\partial S_2} + \left( \frac{\partial S_4}{\partial S_3} \right)^2. \end{aligned} $$

Instead of solving explicitly these coupled partial differential equations, it is wiser to give solutions in parametric form:

$$ S_n = S_n(\Gamma_2, \Gamma_3), $$

for $n = 2, 3, 4, 5$. Solving for the constraints yields:

$$ \begin{aligned} S_2 &= \Gamma_2^3 + \Gamma_2(\kappa - \Gamma_2)\Gamma_3^2, \\ S_3 &= \Gamma_2\Gamma_3(\kappa - \Gamma_2)(3\Gamma_2^2 + (\kappa - 2\Gamma_2)\Gamma_3^2), \\ S_4 &= \frac{9\kappa}{5}\Gamma_2^5 + 6\Gamma_2^3(\kappa - \Gamma_2)^2\Gamma_3^2 \\ &\quad + \Gamma_2(\kappa - \Gamma_2)(\kappa^2 - 3\Gamma_2(\kappa - \Gamma_2))\Gamma_3^4, \\ S_5 &= 9\kappa\Gamma_2^5(k - \Gamma_2)\Gamma_3 + 10\Gamma_2^3(\kappa - \Gamma_2)^3\Gamma_3^3 \\ &\quad + \Gamma_2(\kappa - \Gamma_2)(\kappa - 2\Gamma_2)(\kappa^2 - 2\kappa\Gamma_2 + 2\Gamma_2^2)\Gamma_3^5, \end{aligned} $$

where $\kappa$ is a free parameter. We verify that, as expected, $S_n$ is a polynomial of degree $2n$ in the variables $(\Gamma_2, \Gamma_3)$, and these expressions verify the recurrence relation:

$$ S_{n+1} = \frac{3n}{n+2} S_2 S_{n-1} + \frac{1}{n+2} \left( \frac{\partial S_n}{\partial \Gamma_2} \frac{\partial S_2}{\partial \Gamma_3} + \frac{\partial S_n}{\partial \Gamma_3} \frac{\partial S_2}{\partial \Gamma_2} \right). $$
Under the change of variables $(\rho, u, S_2, S_3) \mapsto (\rho, u, \Gamma_2, \Gamma_3)$, the bracket becomes of the form (17) with

$$g = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix}.$$

Therefore the signature of the bracket is $(1, 1)$. Once again, centering with respect to $\mu_1 = \Gamma_2\Gamma_3$ simplifies the parametric expression of the moments:

$$\mu_1 = \Gamma_2\Gamma_3,$$

$$\mu_2 = \Gamma_2^3 + \kappa \Gamma_2 \Gamma_3^2,$$

$$\mu_3 = \kappa \Gamma_2 \Gamma_3 (3\Gamma_2^2 + \kappa \Gamma_3^2),$$

$$\mu_4 = \kappa \left( \frac{9\Gamma_2^5}{5} + 6\kappa\Gamma_2^3\Gamma_3^2 + \kappa^2\Gamma_2\Gamma_3^4 \right),$$

$$\mu_5 = \kappa^2 \Gamma_2 \Gamma_3 (9\Gamma_2^4 + 10\kappa\Gamma_2^2\Gamma_3^2 + \kappa^2\Gamma_3^4),$$

where we notice that $\mu_n$ is a homogeneous polynomial of degree $n+1$ in the two variables $(\Gamma_2, \Gamma_3)$. The parametric closure satisfies the following recursion relation:

$$\mu_{n+1} = \frac{1}{n+2} \left( \frac{\partial \mu_n}{\partial \Gamma_2} \frac{\partial \mu_2}{\partial \Gamma_3} + \frac{\partial \mu_n}{\partial \Gamma_3} \frac{\partial \mu_2}{\partial \Gamma_2} \right),$$

which is the same recursion relation as in Sec. IV C for the case $N=4$. In addition, we notice that for $\kappa=0$, we recover the case of Sec. IV C, with in particular, $\mu_n = 0$ for $n \ge 4$. However, contrary to the case in Sec. IV C, Casimir invariants in the variables $\mu_n$ are not expressed in a convenient and simple format, given that they require an inversion of a fourth-order polynomial (see Ref. [21] for more details).

From the above recursion relation, we deduce that all the information on the closure is contained in the function $\mu_2$, a homogeneous cubic polynomial. For $N=4$, the following generators of the closures are found:

* 4-field closure (Sec. IV D): $\mu_2 = \Gamma_2^3 + \kappa \Gamma_2 \Gamma_3^2,$

* Burby's closure (Sec. IV C): $\mu_2 = \Gamma_2^3/3,$

* multi-delta closure (Sec. IV A): $\mu_2 = \Gamma_2 \Gamma_3^2.$

Therefore the closure found in Ref. [21] is a linear interpolation between the multi-delta closure and Burby's closure [31]. As in the waterbag case, this closure has the advantage of
having a free parameter $\kappa$ to adjust depending on the kinetic features to be reproduced (see Ref. [21] for more detail).

Given that the closure is composed of homogeneous polynomials $\mu_n$ of degree $n + 1$, the change of variables $(\rho, u, \Gamma_2, \Gamma_3) \mapsto (\rho, u, \mu_1, \mu_2)$ yields a bracket of the form (15) with $\alpha_{nm}(\boldsymbol{\mu}) = (n + m)\mu_{n+m-1}$ and $\beta_{nm}(\boldsymbol{\mu}, \partial_x \boldsymbol{\mu}) = m\partial_x \mu_{n+m-1}$.

## CONCLUSIONS

In this article, we addressed the problem of finding closure functions for the fluid moments $P_n$, i.e., of the form $P_{n \ge N} = P_{n \ge N}(P_0, \dots, P_{N-1})$ by first finding suitable centered moments $\mu_n$ such that the closure reduces to

$$ \mu_{n \ge N-1} = \mu_{n \ge N-1}(\mu_1, \dots, \mu_{N-2}), $$

i.e., closure functions of $N-2$ variables. These closure functions are constrained by the requirement that the resulting dynamical systems should respect the Hamiltonian structure of the parent model, namely the one-dimensional Vlasov equation. We have shown that all these closure functions, when expressed parametrically with well-chosen normal variables $\nu_k$, are generated by a single moment $\mu_2$.

We found that all known hydrodynamic Hamiltonian closures for the one-dimensional Vlasov equation, as listed in Sec. III, share a fundamental property: There exists a convenient parameterization of the variables $\mu_{n < N-1}$ and the closure functions $\mu_{n \ge N-1}$, i.e.,

$$ \mu_n = \mu_n(\nu_1, \dots, \nu_{N-2}), $$

where $\mu_n$ is a polynomial of degree $n + 1$ in the normal variables $\nu_k$. It was shown that the bracket in the $\boldsymbol{\mu}$ variables and the change of variables from $(\rho, u, \mu_1, \dots, \mu_{N-2})$ to $(P_0, P_1, \dots, P_{N-1})$ involve the specification of $\mu_n$ for $n = 1, \dots, 2N - 3$ (by including the closure functions $\mu_n$ appearing explicitly in the bracket). The specification of these polynomials of increasing degree from 2 to $2N - 2$ determines completely the hydrodynamic Hamiltonian closure. We have demonstrated that a given closure is uniquely determined by a single moment, $\mu_2$. For all known closures, this moment is expressed as a cubic polynomial in the $N-2$ variables $\nu_k$. Expressed in the normal variables, the Hamiltonian is given by:

$$ H[\rho, u, \nu_1, \dots, \nu_{N-2}] = \frac{1}{2} \int [\rho u^2 + \rho^3 (\mu_2(\nu_1, \dots, \nu_{N-2}) - \mu_1(\nu_1, \dots, \nu_{N-2})^2) + E^2] \, dx, $$
where $\mu_1$ is a homogeneous quadratic polynomial, defined by Eq. (21). Within this framework, the internal energy density $\rho^2(\mu_2 - \mu_1^2)/2$ is given by the difference between a cubic polynomial in the $N-2$ variables $\nu_k$ and the square of a homogeneous quadratic polynomial in $\nu_k$. Consequently, the specification of the polynomial $\mu_2$ fully characterizes the Hamiltonian system, i.e., Poisson bracket and Hamiltonian. Despite defining a Hamiltonian system, any cubic polynomial does not necessarily lead to a model linked to a reduction of the Vlasov equation. However, following up on a suggestion made in Ref. [31], this family of solutions can be a good starting point for data-driven closures, using sparse identification of non-linear dynamics (SINdy [45]) or sparse physics-informed discovery of empirical relations (SPIDER [46]).

In addition to providing natural variables in which the closures get the simple polynomial form, these normal variables $\nu_k$ greatly simplify the expression of the Poisson bracket:

$$ \{F, G\} = \sum_k \epsilon_k \int \partial_x \frac{\delta F}{\delta \nu_k} \frac{\delta G}{\delta \nu_k} dx, $$

which is very convenient from a computational point of view since it allows for a split of the Poisson bracket as a direct sum of Poisson brackets. Each Poisson bracket individually preserves all the Casimir invariants. The implementation of a split-integrator (by splitting the Poisson bracket, not necessarily the Hamiltonian) ensures the conservation of these quantities.

## ACKNOWLEDGMENTS

The authors would like to thank useful conversations with Josh Burby, Phil Morrison, Maxime Perin and Brad Shadwick. This work has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Programme (Grant Agreement No 101052200 — EUROfusion). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them.
# DATA AVAILABILITY STATEMENT

This manuscript describes an entirely theoretical approach and does not report on any primary data. Data sharing is therefore not applicable.

# AUTHOR CONTRIBUTION STATEMENT

**R. Oufar:** Methodology, Formal analysis, Writing - Original Draft. **C. Chandre:** Conceptualization, Methodology, Formal analysis, Software, Writing - Original Draft.

## Appendix A: Expression of the bracket in the $\mu$-variables

We perform the change of variables $(P_0, P_1, \dots, P_{N-1}) \mapsto (\rho, \psi, \mu_1, \dots, \mu_{N-2})$, given by

$$P_n = \sum_{k=0}^{n} \binom{n}{k} \mu_k \rho^{k+1} \psi^{n-k}, \qquad (A1)$$

obtained from

$$P_n = \int (v - \psi + \psi)^n f(t, x, v) \, dv.$$

The relation (A1) is inverted:

$$\mu_n = \frac{1}{\nu_0^{n+1}} \sum_{k=0}^{n} \binom{n}{k} P_k(-\nu_1)^{n-k}. \qquad (A2)$$

Following Lemma III.1, the moments $P_n$ satisfy the following relations in the variables $\nu$:

$$(n+m)P_{n+m-1} = \frac{\partial P_n}{\partial \nu_k} g_{kl} \frac{\partial P_m}{\partial \nu_l}, \qquad (A3a)$$

$$n \partial_x P_{n+m-1} = \partial_x \frac{\partial P_n}{\partial \nu_k} g_{kl} \frac{\partial P_m}{\partial \nu_l}, \qquad (A3b)$$

for $n, m = 0, \dots, N-1$, and where $g$ is the constant metric. We assume that $g$ takes the block form

$$g = \begin{pmatrix} g_0 & \mathbb{O} \\ \mathbb{O}^t & \bar{g} \end{pmatrix},$$

where $\mathbb{O}$ is the $2 \times (N-2)$ zero matrix, $\bar{g}$ is a $(N-2) \times (N-2)$ constant matrix, and

$$g_0 = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix},$$
to account for the partial decoupling in the variables $\nu_0 = \rho$ and $\nu_1 = \psi$, i.e., $\rho$ is canonically conjugate to $\psi$. Using Eq. (A3a) for $m = 0$, we have

$$ \frac{\partial P_n}{\partial \nu_1} = n P_{n-1}. \qquad (A4) $$

By differentiating Eq. (A1) with respect to $\nu_1$, and using the equation above, we show that

$$ \frac{\partial \mu_n}{\partial \nu_1} = 0. $$

From Eq. (A3b) with $n = 1$, we have

$$ \left( g_{kl} \frac{\partial^2 P_1}{\partial \nu_k \partial \nu_l} - \delta_{kl} \right) \frac{\partial P_m}{\partial \nu_l} = 0, $$

for all $m = 0, \dots, N - 1$. Therefore, $P_1$ is a quadratic polynomial in $\nu$:

$$ P_1 = \frac{1}{2} \nu \cdot g^{-1} \nu + \mathbf{h} \cdot \nu + k, $$

where $\mathbf{h}$ and $k$ are constants. Using a gauge on $\nu$ (arbitrary translation in $\nu$ by a constant vector), we can always assume that $\mathbf{h} = \mathbf{0}$. From Eq. (A3a) with $n = m = 1$, we deduce that $k = 0$. As a consequence, for all closures, we can always find variables $\nu$ such that $P_1$ is a homogeneous quadratic polynomial in these variables, i.e.,

$$ P_1 = \frac{1}{2} \nu \cdot g^{-1} \nu. $$

Ensuring that $P_1$ is a homogeneous quadratic polynomial is fixing the gauge on $\nu$. Since $P_1 = \nu_0\nu_1 + \nu_0^2\mu_1$, the centered moment $\mu_1$ can be rewritten as

$$ \mu_1 = \frac{1}{2} \frac{\nu_k}{\nu_0} (\bar{g}^{-1})_{kl} \frac{\nu_l}{\nu_0}, $$

where the implicit summation over $k, l$ is from 2 to $N - 1$. Next, from Eq. (A2), we compute the partial derivative of $\mu_n$ with respect to $\nu_0$:

$$ \nu_0 \frac{\partial \mu_n}{\partial \nu_0} = -(n+1)\mu_n + \frac{1}{\nu_0^{n+1}} \sum_{k=0}^{n} \binom{n}{k} \nu_0 \frac{\partial P_k}{\partial \nu_0} (-\nu_1)^{n-k}, \qquad (A5) $$

and using Eq. (A3a) with $m = 1$:

$$ \nu_0 \frac{\partial P_k}{\partial \nu_0} = (k+1)P_k - \nu_1 \frac{\partial P_k}{\partial \nu_1} - \nu_l \frac{\partial P_k}{\partial \nu_l}, $$

where the implicit summation over $l$ is from 2 to $N - 1$. The operator $\nu_l \partial/\partial \nu_l$ commutes with the sum over $k$ in the right-hand side of Eq. (A5). As for the operator $\nu_1 \partial/\partial \nu_1$, we
use Eq. (A4), to show that the term of the right-hand side of Eq. (A5) involving $\nu_1\partial/\partial\nu_1$ cancels out with the term involving $(k+1)P_k$. Consequently, we have

$$ \nu_0 \frac{\partial \mu_n}{\partial \nu_0} = -\nu_k \frac{\partial \mu_n}{\partial \nu_k}, $$

where the implicit sum over $k$ is from 2 to $N-1$. Thus leads to the solution:

$$ \mu_n(\nu_0, \nu_2, \dots, \nu_{N-1}) = \bar{\mu}_n \left( \frac{\nu_2}{\nu_0}, \dots, \frac{\nu_{N-1}}{\nu_0} \right). $$

In other words, the centered moments $\mu_n$ for $n \ge 1$ are functions of the rescaled normal variables $\nu_k/\nu_0$ for $k = 2, \dots, N-1$. Next we compute

$$ \nu_0^2 \frac{\partial \mu_n}{\partial \nu_k} \frac{\partial \mu_m}{\partial \nu_l} = \frac{1}{\nu_0^{n+m}} \sum_{q=0}^{n} \sum_{p=0}^{m} \binom{n}{q} \binom{m}{p} \frac{\partial P_q}{\partial \nu_k} \frac{\partial P_p}{\partial \nu_l} (-\nu_1)^{n+m-q-p}, $$

for $n, m \ge 2$ and where the implicit sums over $k$ and $l$ are from 2 to $N-1$. Using Eq. (A3a), we have

$$ \frac{\partial P_n}{\partial \nu_k} g_{kl} \frac{\partial P_m}{\partial \nu_l} = \frac{\partial P_n}{\partial \nu_k} \bar{g}_{kl} \frac{\partial P_m}{\partial \nu_l} + m P_{m-1} \frac{\partial P_n}{\partial \nu_0} + n P_{n-1} \frac{\partial P_m}{\partial \nu_0}, $$

where we have used Eq. (A4). Using the Vandermonde identity,

$$ \sum_{q=0}^{s} \binom{n}{q} \binom{m}{s-q} = \binom{n+m}{s}, $$

the expression simplifies to

$$ \frac{\partial \bar{\mu}_n}{\partial \bar{\nu}_k} \frac{\partial \bar{\mu}_m}{\partial \bar{\nu}_l} = (n+m) \bar{\mu}_{n+m-1} - n \bar{\nu}_{n-1} \gamma_m - m \bar{\nu}_{m-1} \gamma_n, $$

where $\gamma_n$ is given by Eq. (24) and $\bar{\nu} = \nu/\nu_0$. Using similar calculations, we have

$$ \partial_x \frac{\partial \bar{\mu}_n}{\partial \bar{\nu}_k} \frac{\partial \bar{\mu}_m}{\partial \bar{\nu}_l} = n \partial_x \bar{\mu}_{n+m-1} - n \gamma_m \partial_x \nu_{n-1} - m \nu_{m-1} \partial_x \gamma_n. $$

As a consequence, the bracket in the $\bar{\mu}$ variables is a hydrodynamic bracket of the form (15) with elements

$$ \bar{\alpha}_{nm} = (n+m)\bar{\mu}_{n+m-1} - n\bar{\nu}_{n-1}\gamma_m - m\bar{\nu}_{m-1}\gamma_n, \qquad (A6) $$

$$ \bar{\beta}_{nm} = n \partial_x \bar{\mu}_{n+m-1} - n \gamma_m \partial_x \nu_{n-1} - m \nu_{m-1} \partial_x \gamma_n. \qquad (A7) $$
# Appendix B: Proofs of statements in Sec. IV C

## 1. Proof of Lemma IV.2

Since we are looking for solutions of the type defined in Lemma IV.1, $\mu_n^{(m)}$ only depends on $m+1-n$ variables, e.g., $(\nu_n, \dots, \nu_m)$, only the mention of the index of the last variable is needed, e.g., without ambiguity, $\mu_n^{(m)}(\dots, x_l)$ is thus a shorthand for $\mu_n^{(m)}(x_{n-m+l}, \dots, x_l)$. Consequently, the identity (i) can be rewritten as

$$ \partial_k \mu_n^{(m)}(\dots, x_{m+1-n}) = n \mu_{n-1}^{(n+k-2)}(\dots, x_{m+1-n}), \quad (B1) $$

where $\partial_k$ is the partial derivative of $\mu_n^{(m)}$ with respect to the k-th variable (here $x_k$). In particular, we notice that the derivative with respect to the last variable, i.e., for $k = m+1-n$, the variable dependencies do not change between the left and the right-hand side, and

$$ \partial_{m+1-n} \mu_n^{(m)}(x_1, \dots, x_{m+1-n}) = n \mu_{n-1}^{(m-1)}(x_1, \dots, x_{m+1-n}). $$

• (i) implies (ii): We assume that (ii) is satisfied at all lower levels than $m$ (i.e., for $\mu_n^{(k)}$ with $k < m$), and we show that (ii) is satisfied at level $m$ using (i). From Eq. (B1), we compute the k-th partial derivative of $\mu_n^{(m)}$ with respect to $\nu_m$:

$$ \frac{\partial^k}{\partial \nu_m^k} [\mu_n^{(m)}(\dots, \nu_m)] = n(n-1)\dots(n-k+1)\mu_{n-k}^{(m-k)}(\dots, \nu_m), $$

for $1 \le k \le n$, and

$$ \frac{\partial^{n+1}}{\partial \nu_m^{n+1}} [\mu_n^{(m)}(\dots, \nu_m)] = n! \delta_{mn}, $$

$$ \frac{\partial^k}{\partial \nu_m^k} [\mu_n^{(m)}(\dots, \nu_m)] = 0, $$

for $k > n + 1$, where we have used $\mu_0^{(l)}(\dots, x_l) = x_0$. Consequently, $\mu_n^{(m)}$ is a polynomial in the variable $\nu_m$ of degree $n$ if $n < m$, and degree $m+1$ if $n = m$. By writing the polynomial using a series expansion around $\nu_m = 0$, we have for $m > n$,

$$ \begin{aligned} \mu_n^{(m)}(\dots, \nu_m) &= \sum_{k=0}^{n} \frac{\nu_m^{n-k}}{(n-k)!} \frac{\partial^{n-k} \mu_n^{(m)}}{\partial \nu_m^{n-k}}(\dots, \nu_{m-1}, 0), \\ &= \sum_{k=0}^{n} \binom{n}{k} \nu_m^{n-k} \mu_k^{(m-n+k)}(\dots, \nu_{m-1}, 0). \end{aligned} $$
From (ii) at lower levels, i.e., for $k < n$, we have by taking $\nu_m = 0$:

$$ \mu_k^{(m-n+k)}(\dots, \nu_{m-1}, 0) = \mu_k^{(m-n-1)}(\dots, \nu_{m-1}). $$

In order to prove that this identity is also true for $k = n$, we differentiate with respect to $\nu_l$ for $l = n, \dots, m-1$, and show that

$$ \frac{\partial}{\partial \nu_l} \mu_n^{(m)}(\dots, \nu_{m-1}, 0) = \frac{\partial}{\partial \nu_l} \mu_n^{(m-n-1)}(\dots, \nu_{m-1}), $$

since this identity is equivalent to

$$ \mu_{n-1}^{(l-1)}(\dots, \nu_{m-1}, 0) = \mu_{n-1}^{(l-n-1)}(\dots, \nu_{m-1}), $$

which is true for all $l = n, \dots, m-1$, given (ii) at lower levels. Assuming that $\mu_n^{(m)}(0) = 0$, we deduce that

$$ \mu_n^{(m)}(\dots, \nu_{m-1}, 0) = \mu_n^{(m-n-1)}(\dots, \nu_{m-1}). $$

The identity (ii) follows for $n < m$. As for $n = m$, since

$$ \frac{\partial}{\partial \nu_m} [\mu_m^{(m)}(\nu_m)] = m \mu_{m-1}^{(m-1)}(\nu_m), $$

and $\nu_0^{(0)}(\nu_m) = 0$, we deduce that $\nu_m^{(m)}(\nu_m) = \nu_m^{m+1}/(m+1)$ by induction.

• (ii) implies (i):

We assume that (ii) is satisfied for $\mu_n^{(p)}$ where $p = 1, \dots, m-1$ (and $n = 1, \dots, p$). For $k = m$, the partial derivative of $\mu_n^{(m)}$ with respect to $\nu_m$ becomes

$$ \begin{aligned} \frac{\partial \mu_n^{(m)}}{\partial \nu_m}(\dots, \nu_m) &= \sum_{k=0}^{n-1} \binom{n}{k} (n-k) \nu_m^{n-k-1} \mu_k^{(m-n-1)}(\dots, \nu_{m-1}), \\ &= n \sum_{k=0}^{n-1} \binom{n-1}{k} \nu_m^{n-k-1} \mu_k^{(m-n-1)}(\dots, \nu_{m-1}), \\ &= n \mu_{n-1}^{(m-1)}(\dots, \nu_m), \end{aligned} $$

where we have used $(n-k)\binom{n}{k} = n\binom{n-1}{k}$. Next, for $k < m$, the partial derivative of $\mu_n^{(m)}$ with respect to $\nu_l$ is given by

$$ \begin{aligned} \frac{\partial \mu_n^{(m)}}{\partial \nu_l}(\dots, \nu_m) &= \sum_{k=0}^{n} \binom{n}{k} \nu_m^{n-k} \frac{\partial}{\partial \nu_l} [\mu_k^{(m-n-1)}(\dots, \nu_{m-1})], \\ &= \sum_{k=0}^{n} k \binom{n}{k} \nu_m^{n-k} \mu_{k-1}^{(l-n-1)}(\dots, \nu_{m-1}), \\ &= n \mu_{n-1}^{(l-1)}(\dots, \nu_m), \end{aligned} $$

where we have used (i) at strictly lower orders in the form of Eq. (B1), and $k\binom{n}{k} = n\binom{n-1}{k-1}$.
## 2. Proof of Proposition IV.1

Given the functional dependencies of $\mu_k^{(m)}$, we already know that these identities are satisfied for $k+l > m+1$ at any level $m$ (see Lemma IV.1).

• Proof of Eq. (31a):

From (i) of Lemma IV.2, we obtain

$$ \frac{\partial \mu_k^{(m)}}{\partial \nu_j} \frac{\partial \mu_l^{(m)}}{\partial \nu_{m+1-j}} = kl \sum_{j=k}^{m-l+1} \mu_{k-1}^{(j-1)}(\dots, \nu_m) \mu_{l-1}^{(m-j)}(\dots, \nu_m). $$

We extend the summation to all $j$ by assuming that $\mu_k^{(j)} = 0$ if $k > j$ or $j < 0$. Below the index $j$ is implicitly summed. By using (ii) of Lemma IV.2, the summation is further expanded to

$$ \begin{aligned} \frac{\partial \mu_k^{(m)}}{\partial \nu_j} \frac{\partial \mu_l^{(m)}}{\partial \nu_{m+1-j}} &= kl \sum_{k', l'} \binom{k-1}{k'} \nu_m^{k-1-k'} \mu_{k'}^{(j-k-1)}(\dots, \nu_{m-1}) \binom{l-1}{l'} \nu_m^{l-1-l'} \mu_{l'}^{(m-j-l)}(\dots, \nu_{m-1}), \\ &= kl \sum_{k', l'} \binom{k-1}{k'} \binom{l-1}{l'} \nu_m^{k-1-k'} \nu_m^{l-1-l'} \mu_{k'}^{(j-k-1)}(\dots, \nu_{m-1}) \mu_{l'}^{(m-j-l)}(\dots, \nu_{m-1}). \end{aligned} $$

Next we use (ii) at levels strictly smaller than $m$:

$$ \mu_{k'}^{(j-k-1)} \mu_{l'}^{(m-j-l)} = \frac{k' + l' + 2}{(k' + 1)(l' + 1)} \mu_{k' + l' + 1}^{(m-k-l)}. $$

Using $k \binom{k-1}{k'} = (k' + 1) \binom{k}{k'+1}$, we have

$$ \frac{\partial \mu_k^{(m)}}{\partial \nu_j} \frac{\partial \mu_l^{(m)}}{\partial \nu_{m+1-j}} = \sum_{k'} \sum_s \binom{k}{k' + 1} \binom{l}{s - k' + 1} \nu_m^{k+l-s-2} (s+2) \mu_{s+1}^{m-k-l}. $$

We simplify the product of the binomial coefficients using the Vandermonde identity

$$ \sum_{k'} \binom{k}{k' + 1} \binom{l}{s - k' + 1} = \binom{k + l}{s + 2}, $$

to get

$$ \frac{\partial \mu_k^{(m)}}{\partial \nu_j} \frac{\partial \mu_l^{(m)}}{\partial \nu_{m+1-j}} = \sum_s (s+2) \binom{k+l}{s+2} \nu_m^{k+l-s-2} \mu_{s+1}^{(m-k-l)}, $$

which is equal to $(k+l)\mu_{k+l-1}^{(m)}$ since $(s+2)\binom{k+l}{s+2} = (k+l)\binom{k+l-1}{s+1}$.

To initiate the induction, we prove that the identity (31a) is trivially satisfied at level $m=1$ where there is only one non-zero polynomial, $\nu_1^{(1)}(\nu_1) = \nu_1^2/2$.

• Proof of Eq. (31b):
This proof follows the same lines as for Eq. (31a). Proving Eq. (31b) at level m is equivalent to proving:

$$ \frac{\partial^2 \mu_k^{(m)}}{\partial \nu_n \partial \nu_j} \frac{\partial \mu_l^{(m)}}{\partial \nu_{m+1-j}} = k \frac{\partial \mu_{k+l-1}^{(m)}}{\partial \nu_n}. $$

We make use of the identities (i) of Lemma IV.2, from which we deduce:

$$ \begin{aligned} \frac{\partial}{\partial \nu_n} [\mu_k^{(m)}(\dots, \nu_m)] &= k \mu_{k-1}^{(n-1)}(\dots, \nu_m), \\ \frac{\partial^2}{\partial \nu_n \partial \nu_j} [\mu_k^{(m)}(\dots, \nu_m)] &= k(k-1) \mu_{k-2}^{(n-2+j-m)}(\dots, \nu_m), \\ \frac{\partial}{\partial \nu_{m+1-j}} [\mu_l^{(m)}(\dots, \nu_m)] &= l \mu_{l-1}^{(m-j)}(\dots, \nu_m). \end{aligned} $$

Furthermore, we use (ii) from Lemma IV.2, i.e.,

$$ \mu_{k-2}^{(n-2+j-m)} \mu_{l-1}^{(m-j)} = \frac{k+l-1}{(k-1)l} \mu_{k+l-2}^{(n-1)}. $$

Consequently,

$$ \begin{aligned} \frac{\partial^2 \mu_k^{(m)}}{\partial \nu_n \partial \nu_j} \frac{\partial \mu_l^{(m)}}{\partial \nu_{m+1-j}} &= k(k+l-1)\mu_{k+l-2}^{(n-1)}, \\ &= k \frac{\partial}{\partial \nu_n} [\mu_{k+l-1}^{(m)}(\dots, \nu_m)], \end{aligned} $$

where we have used (i) of Lemma IV.2.

To initiate the recursion, we also consider the lowest level m = 1 where the identity (31b) is trivially satisfied.

## Appendix C: Python code for the computation of Casimir invariants associated with Burby's closure [31]

```python
import sympy as sp
```

```python
def compute mu(n, l, *args):
    if len(args) == 0 or l < 0 or l > n:
        return 0
    if len(args) != n - l + 1:
        raise ValueError(f"Dimension of mu({n},{l}) should be {n - l + 1}
                        and not {len(args)}")
    mu = args[0] * args[-1]**l if len(args)>=2 else args[0]**(l + 1) / (l + 1)
    for k in range(1, l + 1):
```
```
mu += sp.binomial(1, k) * args[-1]**(1 - k)\n
    * compute Mu(n - 1 - 1, k, *args[k:-1])
return mu
```

def compute_vec Mu(n, *args):

```
mu = [None] * n

for l in range(n):
    mu[l] = sp.simplify(compute Mu(n, l + 1, *args[1:]))

return mu
```

def inverse Mu(n, *args):

```
x_ = spalletx(sp. symbols(f"y1: {n+1}"))
x_[-1] = ((n + 1) * args[-1])**sp. Rational(1, n + 1)

for 1 in range(2, n + 1):
    x_[-1] = args[-1] - sp. simplify(compute Mu(n, n - 1 + 1, *x_[-1]))\n
    - x_[-1] * x_[-1]**(n - 1 + 1))

x_[-1] += x_[-1]**(1 - 1 - n)

return [eq for row in sp. simplify(x_). tolist() for eq in row]
```

def print_var (vars,**kwargs):

```
output, name = kwargs.get("output", "normal"), kwargs.get("name", "x")
vars = [vars] if isinstance (vars, sp. Basic) else vars

for var in vars:
    if output == "normal":

```
m = 3
x = sp. symbols(f"x1: {m + 1}")
print_var(compute_vec Mu(m, *x))
print_var(compute_vec Mu(m, *x), output="latex", name=r"\nu")
print_var(inverse mu(m, *x))
print_var(inverse mu(m, *x), output="latex", name=r"\mu")
```

The first levels $m$ are given below:

• $m=2$:

$$ \mu_1^{(2)} = \nu_1 \nu_2, $$

$$ \mu_2^{(2)} = \frac{\nu_2^3}{3}. $$
• $m = 3$:

$$\mu_1^{(3)} = \nu_1 \nu_3 + \frac{\nu_2^2}{2},$$

$$\mu_2^{(3)} = \nu_2 \nu_3^2,$$

$$\mu_3^{(3)} = \frac{\nu_3^4}{4}.$$

• $m = 4$:

$$\mu_1^{(4)} = \nu_1 \nu_4 + \nu_2 \nu_3,$$

$$\mu_2^{(4)} = \nu_2 \nu_4^2 + \nu_3^2 \nu_4,$$

$$\mu_3^{(4)} = \nu_3 \nu_4^3,$$

$$\mu_4^{(4)} = \frac{\nu_4^5}{5}.$$

• $m = 5$:

$$\mu_1^{(5)} = \nu_1 \nu_5 + \nu_2 \nu_4 + \frac{\nu_3^2}{2},$$

$$\mu_2^{(5)} = \nu_2 \nu_5^2 + 2 \nu_3 \nu_4 \nu_5 + \frac{\nu_4^3}{3},$$

$$\mu_3^{(5)} = \nu_3 \nu_5^3 + \frac{3}{2} \nu_4^2 \nu_5^2,$$

$$\mu_4^{(5)} = \nu_4 \nu_5^4,$$

$$\mu_5^{(5)} = \frac{\nu_5^6}{6}.$$

Their inversion provide the explicit expression of the density of the Casimir invariants $C = \int \rho \nu_n(\mu) dx$:

• $m = 2$:

$$\nu_1^{(2)} = \frac{3^{\frac{2}{3}} \mu_1}{3^{\sqrt{3}} \mu_2},$$

$$\nu_2^{(2)} = \sqrt[3]{3} \sqrt[3]{\mu_2}.$$

• $m = 3$:

$$\nu_1^{(3)} = \frac{\sqrt{2} (8 \mu_1 \mu_3 - \mu_2^2)}{16 \mu_3^{\frac{5}{4}}},$$

$$\nu_2^{(3)} = \frac{\mu_2}{2\sqrt{\mu_3}},$$

$$\nu_3^{(3)} = \sqrt{2} \sqrt[4]{\mu_3}.$$
• m = 4:

$$ \nu_1^{(4)} = \frac{5^{\frac{4}{5}} (25\mu_1\mu_4^2 - \mu_3(5\mu_2\mu_4 - \mu_3^2))}{125\mu_4^{\frac{11}{5}}}, $$

$$ \nu_2^{(4)} = \frac{5^{\frac{3}{5}} (5\mu_2\mu_4 - \mu_3^2)}{25\mu_4^{\frac{7}{5}}}, $$

$$ \nu_3^{(4)} = \frac{5^{\frac{2}{5}}\mu_3}{5\mu_4^{\frac{3}{5}}}, $$

$$ \nu_4^{(4)} = \sqrt[5]{5\sqrt[5]{\mu_4}}. $$

• m = 5:

$$ \nu_1^{(5)} = \frac{6^{\frac{5}{6}} (5184\mu_1\mu_5^3 - 864\mu_2\mu_4\mu_5^2 - 432\mu_3^2\mu_5^2 + 504\mu_3\mu_4^2\mu_5 - 91\mu_4^4)}{31104\mu_5^{\frac{19}{6}}}, $$

$$ \nu_2^{(5)} = \frac{6^{\frac{2}{3}} (27\mu_2\mu_5^2 - \mu_4(9\mu_3\mu_5 - 2\mu_4^2))}{162\mu_5^{\frac{7}{3}}}, $$

$$ \nu_3^{(5)} = \frac{\sqrt{6}(4\mu_3\mu_5 - \mu_4^2)}{24\mu_5^{\frac{3}{2}}}, $$

$$ \nu_4^{(5)} = \frac{\sqrt[3]{6}\mu_4}{6\mu_5^{\frac{2}{3}}}, $$

$$ \nu_5^{(5)} = \sqrt[6]{6\sqrt[6]{\mu_5}}. $$

REFERENCES

¹R. D. Hazeltine and F. L. Waelbroeck, *The Framework of Plasma Physics*, Frontiers in Physics (Westview Press, Boulder, Colorado, 2004).

²S. I. Braginskii, *Reviews of Plasma Physics* **1**, 205 (1965).

³E. Ott and R. N. Sudan, *The Physics of Fluids* **12**, 2388 (1969).

⁴R. D. Hazeltine, M. Kotschenreuther, and P. J. Morrison, *The Physics of Fluids* **28**, 2466 (1985).

⁵G. W. Hammett and F. W. Perkins, *Phys. Rev. Lett.* **64**, 3019 (1990).

⁶G. W. Hammett, M. A. Beer, W. Dorland, S. C. Cowley, and S. A. Smith, *Plasma Physics and Controlled Fusion* **35**, 973 (1993).

⁷H. Sugama, T.-H. Watanabe, and W. Horton, *Physics of Plasmas* **10**, 726 (2003).

⁸T. Passot and P. L. Sulem, *Nonlinear Processes in Geophysics* **11**, 245 (2004).
⁹B. A. Shadwick, G. M. Tarkenton, and E. H. Esarey, Phys. Rev. Lett. **93**, 175002 (2004).

¹⁰B. A. Shadwick, G. M. Tarkenton, E. Esarey, and C. B. Schroeder, Physics of Plasmas **12**, 056710 (2005).

¹¹P. Goswami, T. Passot, and P. L. Sulem, Physics of Plasmas **12**, 102109 (2005).

¹²Y. Sarazin, G. Dif-Pradalier, D. Zarzoso, X. Garbet, P. Ghendrih, and V. Grandgirard, Plasma Physics and Controlled Fusion **51**, 115003 (2009).

¹³B. Shadwick, G. Tarkenton, E. Esarey, and F. M. Lee, Communications in Nonlinear Science and Numerical Simulation **17**, 2153 (2012).

¹⁴D. C. Depackh, Journal of Electronics and Control **13**, 417 (1962).

¹⁵P. Bertrand and M. R. Feix, Physics Letters A **28**, 68 (1968).

¹⁶L. Gosse, S. Jin, and X. Li, Mathematical Models and Methods in Applied Sciences **13**, 1689 (2003).

¹⁷R. O. Fox, Journal of Computational Physics **228**, 7771 (2009).

¹⁸L. de Guillebon and C. Chandre, Physics Letters A **376**, 3172 (2012).

¹⁹M. Perin, C. Chandre, P. J. Morrison, and E. Tassi, Annals of Physics **348**, 50 (2014).

²⁰M. Perin, C. Chandre, P. J. Morrison, and E. Tassi, Journal of Physics A: Mathematical and Theoretical **48**, 275501 (2015).

²¹C. Chandre and B. A. Shadwick, Physics of Plasmas **29**, 102101 (2022).

²²P. J. Morrison, AIP Conference Proceedings **88**, 13 (1982).

²³P. J. Morrison, Reviews of Modern Physics **70**, 467 (1998).

²⁴B. A. Kupershmidt and Y. I. Manin, Functional Analysis and Its Applications **12**, 20 (1978).

²⁵B. A. Kupershmidt, Physics Letters A **121**, 167 (1987).

²⁶B. A. Dubrovin and S. P. Novikov, Soviet Mathematics. Doklady **27**, 665 (1983).

²⁷O. I. Mokhov, Functional Analysis and Its Applications **22**, 336 (1988).

²⁸B. A. Dubrovin and S. P. Novikov, Russian Mathematical Surveys **44**, 35 (1989).

²⁹O. I. Mokhov, Functional Analysis and Its Applications **42**, 33 (2008).

³⁰M. Perin, C. Chandre, P. J. Morrison, and E. Tassi, Physics of Plasmas **22**, 092309 (2015).

³¹J. W. Burby, Scientific Reports **13**, 18286 (2023).

³²C. Yuan and R. O. Fox, Journal of Computational Physics **230**, 8216 (2011).

³³C. Chalons, D. Kah, and M. Massot, Communications in Mathematical Sciences **10**, 1241 (2012).
<footnote n="34">³⁴Y. Cheng and J. A. Rossmanith, Journal of Computational and Applied Mathematics **262**, 384 (2014).</footnote>

<footnote n="35">³⁵A. Inglebert, A. Ghizzo, T. Reveille, D. Del Sarto, P. Bertrand, and F. Califano, Euro-physics Letters **95**, 45002 (2011).</footnote>

<footnote n="36">³⁶A. Inglebert, A. Ghizzo, T. Reveille, D. Del Sarto, P. Bertrand, and F. Califano, Plasma Physics and Controlled Fusion **54**, 085004 (2012).</footnote>

<footnote n="37">³⁷A. Inglebert, A. Ghizzo, T. Reveille, P. Bertrand, and F. Califano, Physics of Plasmas **19**, 122109 (2012).</footnote>

<footnote n="38">³⁸M. Antoine, A. Ghizzo, D. D. Sarto, and E. Deriaz, Journal of Plasma Physics **91**, E131 (2025).</footnote>

<footnote n="39">³⁹P. Bertrand, G. Baumann, and M. R. Feix, Physics Letters A **29**, 489 (1969).</footnote>

<footnote n="40">⁴⁰P. Bertrand and M. R. Feix, Plasma Physics **18**, 655 (1976).</footnote>

<footnote n="41">⁴¹L. Yu, Journal of Physics A: Mathematical and General **33**, 8127 (2000).</footnote>

<footnote n="42">⁴²A. A. Chesnokov and M. V. Pavlov, Acta Applicandae Mathematicae **122**, 367 (2012).</footnote>

<footnote n="43">⁴³P. J. Morrison and G. I. Hagstrom, Continuum Hamiltonian Hopf bifurcation I, in *Non-linear Physical Systems* (2013) Chap. 12, pp. 247–281.</footnote>

<footnote n="44">⁴⁴G. I. Hagstrom and P. J. Morrison, Continuum Hamiltonian Hopf bifurcation II, in *Non-linear Physical Systems* (2013) Chap. 13, pp. 283–310.</footnote>

<footnote n="45">⁴⁵S. L. Brunton and J. N. Kutz, *Data-Driven Science and Engineering: Machine Learning, Dynamical Systems, and Control*, 2nd ed. (Cambridge University Press, 2022).</footnote>

<footnote n="46">⁴⁶D. R. Gurevich, M. R. Golden, P. A. Reinbold, and R. O. Grigoriev, Journal of Fluid Mechanics **996**, A25 (2024).</footnote>