### A canonical Hamiltonian formulation of the Navier–Stokes problem

1. Introduction Given the title of this paper, it is incumbent on the authors to assure the reader that we do not claim to have done the impossible. A viscous fluid is, after all, a non-Hamiltonian system (Millikan Reference Millikan1929; Finlayson Reference Finlayson1972a , Reference Finlayson b ). There is no action integral for which

## 1. Introduction

Given the title of this paper, it is incumbent on the authors to assure the reader that we do not claim to have done the impossible. A viscous fluid is, after all, a non-Hamiltonian system (Millikan Reference Millikan1929; Finlayson Reference Finlayson1972*a*,Reference Finlayson*b*). There is no action integral for which Hamilton’s principle (Hamilton Reference Hamilton1833, Reference Hamilton1834, Reference Hamilton1835) yields the Navier–Stokes equations (Stokes Reference Stokes1845; Anderson, Tannehill & Pletcher Reference Anderson, Tannehill and Pletcher1984; Anderson Reference Anderson1995; Batchelor Reference Batchelor2000; White Reference White2006; Pozrikidis Reference Pozrikidis2009; Cengel & Cimbala Reference Cengel and Cimbala2018) in their usual form (Millikan Reference Millikan1929; Finlayson Reference Finlayson1972*a*,Reference Finlayson*b*), and we do not claim otherwise. Remarkably, however, a Hamiltonian formulation can still be found by considering a mathematically equivalent higher-order problem, as we will now demonstrate via simple example.

### 1.1. A motivating example

Consider the first-order initial-value problem

(1.1*a*,*b*)begin{equation} dot{v}={-}v, quad v(0)=1, end{equation}

with unique solution $v(t)={rm e}^{-t}$. Here, $v(t)$ can be interpreted as the velocity of a lumped mass moving in a viscous medium in one dimension with linear damping. Like the traditional Navier–Stokes equations (Stokes Reference Stokes1845; Anderson *et al.* Reference Anderson, Tannehill and Pletcher1984; Anderson Reference Anderson1995; Batchelor Reference Batchelor2000; White Reference White2006; Pozrikidis Reference Pozrikidis2009; Cengel & Cimbala Reference Cengel and Cimbala2018), this too is an intrinsically non-Hamiltonian problem, in that there is no action $mathcal {S}$ for which Hamilton’s principle ($delta mathcal {S}=0$) yields the governing equation $dot {v}=-v$. Yet, if we simply differentiate both sides of the equation ($ddot {v}=-dot {v}$), use the original equation to write $dot {v}=-v$ and apply the additional initial condition $dot {v}(0)=-v(0)=-1$, we arrive at the mathematically equivalent second-order problem

(1.2*a*–*c*)begin{equation} ddot{v}=v, quad v(0)=1, quad dot{v}(0)={-}1, end{equation}

which has the same unique solution $v(t)={rm e}^{-t}$ but which is Hamiltonian – not in the sense that the total mechanical energy is conserved, but in the sense that it has mathematically Hamiltonian structure.

As first observed by Sanders (Reference Sanders2021, Reference Sanders2022, Reference Sanders2023*a*,Reference Sanders*b*) and Sanders & Inman (Reference Sanders and Inman2023), the associated action can be obtained by writing the original equation in standard form ($mathcal {R}equiv dot {v}+v=0$), squaring the residual $mathcal {R}$ and integrating over time

(1.3)

begin{align}

mathcal{S}^{*}[v]=inttext{d}tbig(tfrac{1}{2}mathcal{R}^{2}big)=int

text{d}tbig[tfrac{1}{2}big(dot{v}^{2}+2vdot{v}+v^{2}big)big]
siminttext{d}tbig[tfrac{1}{2}big(dot{v}^{2}+v^{2}big)big],

end{align}

where we have used the fact that $2vdot {v}=text {d}(v^{2})/text {d}t$ is a total time derivative and can therefore be excluded from the action without changing the resulting Euler–Lagrange equation (Lanczos Reference Lanczos1970). This is the so-called ‘time-averaged principle of least squares’ (Sanders Reference Sanders2021, Reference Sanders2022, Reference Sanders2023*a*,Reference Sanders*b*; Sanders & Inman Reference Sanders and Inman2023): since $mathcal {R}=0$ is a local minimum of $mathcal {R}^{2}$, it is also a local minimum of $int text {d}t(mathcal {R}^{2})$. Varying $v$, the first variation of $mathcal {S}^{*}$ is

(1.4)begin{equation} deltamathcal{S}^{*}=inttext{d}tleft[dot{v}deltadot{v}+vdelta vright]=inttext{d}tleft[(-ddot{v}+v)delta vright]+left[dot{v}delta vright]_{t_{1}}^{t_{2}}, end{equation}

yielding the second-order equation $ddot {v}=v$ and revealing the canonically conjugate ‘momentum’ ${rm pi} equiv dot {v}$. Here, and in what follows, we will use the symbol ${rm pi}$ for canonically conjugate momenta, as is customary in Hamiltonian field theory, in order to avoid later confusion with the pressure field $p$. Since the mathematical constant $3.14159ldots$ does not appear in the present work, there will be no ambiguity.

The corresponding Hamiltonian is obtained via the Legendre transform

(1.5)

begin{equation}

H^{*}[v,{rm pi}]={rm pi}dot{v}-tfrac{1}{2}big(dot{v}^{2}+v^{2}big)=tfrac{1}{2}big({rm pi}^{2}-v^{2}big).

end{equation}

Notably, this Hamiltonian has nothing to do with the total mechanical energy of the system, although it is a conserved quantity. In fact, $H^{*}=0$ for the actual motion satisfying ${rm pi} equiv dot {v}=-v$. We note in passing that Liouville’s theorem is satisfied, as the motion occurs along the line ${rm pi} =-v$, so that the phase-space volume, being always zero, is conserved. Hamilton’s equations

(1.6*a*,*b*)begin{equation} dot{v}=frac{partial H^{*}}{partial{rm pi}}, quad dot{rm pi}={-}frac{partial H^{*}}{partial v}, end{equation}

are mathematically equivalent to the second-order problem $ddot {v}=v$ and therefore also mathematically equivalent to the original, first-order problem.

The associated Hamilton–Jacobi equation (Hamilton Reference Hamilton1833, Reference Hamilton1834, Reference Hamilton1835; Jacobi Reference Jacobi1837, Reference Jacobi1842–1843; Whittaker Reference Whittaker1904; Lanczos Reference Lanczos1970) is

(1.7)begin{equation} frac{1}{2}left(frac{partial{S}^{*}}{partial v}right)^{2}-frac{1}{2}v^{2}+frac{partial{S}^{*}}{partial t}=0, end{equation}

where Hamilton’s principal function ${S}^{*}={S}^{*}(v,t)$ serves as the generating function for a canonical transformation to a new coordinate $phi$ which is constant and equal to its initial value. Although this is almost identical to the Hamilton–Jacobi equation for the simple harmonic oscillator – the only difference being the sign in front of $(1/2)v^{2}$ – the usual separable solution of the form ${S}^{*}(v,t)=W(v)+T(t)$, where $W(v)$ and $T(t)$ are functions of *v* and *t* (respectively), does not work, as the reader may check.

Instead, let us use a trial solution of the form

(1.8)begin{equation} {S}^{*}(v,t)=F(t)v+tfrac{1}{2}v^{2}+f(t), end{equation}

where $F(t)$ and $f(t)$ are as yet undetermined functions of $t$. This trial solution was chosen to cancel the term $(1/2)v^{2}$ from the equation. Substituting our trial solution into the Hamilton–Jacobi equation, we find that

(1.9)begin{equation} tfrac{1}{2}[F(t)]^{2}+[F(t)+F'(t)]v+f'(t)=0. end{equation}

In order for this equation to hold for all $v$, we must have the following:

(1.10)begin{equation} F(t)+F'(t)=0 quad Rightarrow quad F(t)=alpha {rm e}^{{-}t}, end{equation}

where $alpha$ is a constant of integration which will be used to transform to the new coordinate, and

(1.11)begin{equation} tfrac{1}{2}[F(t)]^{2}+f'(t)=0 quad Rightarrow quad f(t)=tfrac{1}{4}alpha^{2}{rm e}^{{-}2t}+gamma, end{equation}

where $gamma$ is another constant of integration which is simply additive and can therefore be discarded.

In this way, we have that

(1.12)begin{equation} {S}^{*}(v,t;alpha)=alpha {rm e}^{{-}t}v+tfrac{1}{2}v^{2}+tfrac{1}{4}alpha^{2}{rm e}^{{-}2t}. end{equation}

With one constant of integration ($alpha$) to match the single degree of freedom ($v$), this is a complete solution to the Hamilton–Jacobi equation. The new coordinate $phi$ (which is constant and equal to its initial value) is obtained via the canonical transformation

(1.13)begin{equation} phi=frac{partial{S}^{*}}{partialalpha}={rm e}^{{-}t}v+frac{1}{2}alpha {rm e}^{{-}2t}. end{equation}

The numerical value of $alpha$ is in turn obtained via the canonical relation

(1.14)begin{equation} {rm pi}=frac{partial{S}^{*}}{partial v}=alpha {rm e}^{{-}t}+v, end{equation}

which, evaluated at $t=0$, gives $alpha =-2$ (recall that ${rm pi} =dot {v}$, and $dot {v}(0)=-v(0)=-1$). Using the fact that the new coordinate $phi$ is equal to its initial value, we have that

(1.15)begin{equation} {rm e}^{{-}t}v-{rm e}^{{-}2t}=v(0)-1=0, end{equation}

giving the correct solution $v(t)={rm e}^{-t}$.

In summary, by doubling the order of the governing equation and supplying additional auxiliary conditions, we made a non-Hamiltonian problem into a Hamiltonian one (Sanders Reference Sanders2021, Reference Sanders2022, Reference Sanders2023*a*,Reference Sanders*b*; Sanders & Inman Reference Sanders and Inman2023). Furthermore, this simple example demonstrates that the method correctly gives the solution to the original, non-Hamiltonian problem. Indeed, it would appear that every non-Hamiltonian problem belongs to an equivalence class of problems with the same solution, and within each such equivalence class there are Hamiltonian variants. The remainder of this paper applies that concept to the isotropic Navier–Stokes problem (Stokes Reference Stokes1845; Anderson *et al.* Reference Anderson, Tannehill and Pletcher1984; Anderson Reference Anderson1995; Batchelor Reference Batchelor2000; White Reference White2006; Pozrikidis Reference Pozrikidis2009; Cengel & Cimbala Reference Cengel and Cimbala2018).

### 1.2. The Navier–Stokes problem

The incompressible Navier–Stokes equations (Stokes Reference Stokes1845; Anderson *et al.* Reference Anderson, Tannehill and Pletcher1984; Anderson Reference Anderson1995; Batchelor Reference Batchelor2000; White Reference White2006; Pozrikidis Reference Pozrikidis2009; Cengel & Cimbala Reference Cengel and Cimbala2018) are given by

(1.16)$$begin{gather} rho dot{u}_{i} + rho u_{i,j}u_{j} + p_{,i} – mu u_{i,jj} – rho b_{i} = 0, end{gather}$$

(1.17)$$begin{gather}u_{i,i} = 0, end{gather}$$

where $rho$ is the constant and uniform density, $u_{i}=u_{i}(x_{j},t)$ is the velocity field, $p=p(x_{j},t)$ is the pressure field, $b_{i}=b_{i}(x_{j},t)$ is the body force field, subscript Roman indices label Euclidean tensor components ($i,j=1,2,3$), the $x_{j}$ are Eulerian spatial coordinates, $t$ is time, $mu$ is the dynamic viscosity, a dot over a symbol denotes a partial time derivative ($dot {u}_{i}=partial {u}_{i}/partial t$), a comma in a subscript indicates a spatial gradient ($p_{,i}=partial p/partial x_{i}$) and we employ the Einstein summation convention on repeated subscript indices. To be clear, the notation $u_{i}(x_{j},t)$ signifies that each component of the velocity field is a function of all three spatial coordinates $(x_{1},x_{2},x_{3})$ and time $t$. It has the same meaning as other common notations, such as $u_{i}(boldsymbol {x},t)$ and $boldsymbol {u}(boldsymbol {x},t)$. Likewise for all other field quantities. In the case of a uniform gravitational field, $b_{i}=g_{i}$ coincides with the local acceleration due to gravity; however, in what follows, we make no assumptions about the functional form of $b_{i}(x_{j},t)$: it is completely arbitrary. There are four unknown field quantities: $u_{i}(x_{j},t)$ and $p(x_{j},t)$.

We seek, ultimately, a functional

(1.18)begin{equation} H^{*}=H^{*}[u_{i},p,{rm pi}_{j},{rm pi}_{4};t], end{equation}

where (${rm pi} _{i},{rm pi} _{4}$) are suitable ‘momenta’ conjugate to the field quantities ($u_{i},p$), such that Hamilton’s canonical equations

(1.19*a*,*b*)$$begin{gather} dot{u}_{i}=frac{delta H^{*}}{delta{rm pi}_{i}}, quaddot{p}=frac{delta H^{*}}{delta{rm pi}_{4}}, end{gather}$$

(1.20*a*,*b*)$$begin{gather}dot{rm pi}_{i}={-}frac{delta H^{*}}{delta u_{i}}, quaddot{rm pi}_{4}={-}frac{delta H^{*}}{delta p}, end{gather}$$

constitute a mathematically equivalent second-order formulation of the problem, where $delta H^{*}/delta u_{i}$, $delta H^{*}/delta p$, $delta H^{*}/delta {rm pi}_{i}$ and $delta H^{*}/delta {rm pi}_{4}$ are the Volterra (Reference Volterra1930) functional derivatives of $H^{*}$ with respect to the field quantities and the conjugate momenta. We will find that this is generally possible for a compressible fluid. For an incompressible fluid, the equation $dot {p}=delta H^{*}/delta {rm pi}_{4}$ will need to be replaced by the incompressibility condition $u_{i,i}=0$, consistent with the well-known result that the pressure usually serves as a Lagrange multiplier for the incompressibility constraint (Lanczos Reference Lanczos1970; Badin & Crisciani Reference Badin and Crisciani2018) (refer to p. 361 of Lanczos (Reference Lanczos1970) and pp. 137 and 141 of Badin & Crisciani Reference Badin and Crisciani2018).

The remainder of this paper is organized as follows. Section 2 gives a comprehensive overview of the relevant literature to date. Sections 3 and 4 contain the main results of the present work, culminating in a conserved Hamiltonian functional $H^{*}$ satisfying Hamilton’s equations (1.19*a*,*b*) and (1.20*a*,*b*) for the mathematically equivalent second-order problem, along with the accompanying Hamilton–Jacobi equation (Hamilton Reference Hamilton1833, Reference Hamilton1834, Reference Hamilton1835; Jacobi Reference Jacobi1837, Reference Jacobi1842–1843; Whittaker Reference Whittaker1904; Lanczos Reference Lanczos1970). Section 5 contains a discussion of the physical interpretation of the second-order formulation. Section 6 presents a brief case study in the form of one-dimensional flow over an infinite, flat plate. Finally, § 7 concludes the paper with a few closing remarks and an outline of how the present formulation can aid in resolving the question of existence and uniqueness of solutions to the Navier–Stokes problem.

By the end of the paper, we will have achieved precisely what the title promises: a canonical Hamiltonian formulation of the problem, opening new avenues toward resolution of one of the most famous unsolved problems in mathematics.

## 2. Literature review

The field of analytical mechanics, with foundations planted in Hamilton’s principle of stationary action (Hamilton Reference Hamilton1833, Reference Hamilton1834, Reference Hamilton1835) or d’Alembert’s principle of virtual work (d’Alembert Reference d’Alembert1743), has been vital to the development of both classical and quantum physics since the eighteenth century. This approach is versatile and helpful to the physical understanding of the problem in question, and the foundation, structure and utility of Hamiltonian formalism is well documented (Becker Reference Becker1954; Taylor Reference Taylor2005; Hamill Reference Hamill2014; Bohn Reference Bohn2018; Cline Reference Cline2023; Fowler Reference Fowler2023). The supporting mathematics of the calculus of variations as well as symplectic and differential geometry can also be found in many excellent sources (Arnold Reference Arnold1989; Berndt Reference Berndt2001; Hall Reference Hall2003; Boas Reference Boas2006; Stone & Goldbart Reference Stone and Goldbart2009; Gelfand & Fomin Reference Gelfand and Fomin2012; Arfken Reference Arfken2013; Needham Reference Needham2021). It is therefore no surprise that researchers have been applying analytical formalism to classical fluids dating back to the time of Lagrange (Lagrange Reference Lagrange1811; Lichtenstein Reference Lichtenstein1929; Morrison Reference Morrison1998, Reference Morrison2006; Berdichevsky Reference Berdichevsky2009*a*,Reference Berdichevsky*b*; dell’Isola & Gavrilyuk Reference dell’Isola and Gavrilyuk2011; Badin & Crisciani Reference Badin and Crisciani2018; Bedford Reference Bedford2021).

The task of obtaining solutions to the governing equations of fluid flow represents one of the most challenging problems in science and engineering. In most cases, the mathematical formulation is expressed as an initial-boundary-value problem: a set of coupled, nonlinear partial differential equations, which are to be solved subject to various initial and boundary conditions. The degree of complication of the governing equations depends on the type of the fluid. For a viscous fluid where the transport phenomena of friction and thermal conduction are included, the governing equations are called the Navier–Stokes equations (Stokes Reference Stokes1845; Anderson *et al.* Reference Anderson, Tannehill and Pletcher1984; Anderson Reference Anderson1995; Batchelor Reference Batchelor2000; White Reference White2006; Pozrikidis Reference Pozrikidis2009; Cengel & Cimbala Reference Cengel and Cimbala2018). The Navier–Stokes equations are derived by applying fundamental physical principles – conservation of mass, conservation of momentum and conservation of energy – to a viscous fluid, and the derivation can be found in any fluid mechanics textbook (Anderson *et al.* Reference Anderson, Tannehill and Pletcher1984; Anderson Reference Anderson1995; Batchelor Reference Batchelor2000; White Reference White2006; Pozrikidis Reference Pozrikidis2009; Cengel & Cimbala Reference Cengel and Cimbala2018). It is important to recognize that the Navier–Stokes equations as they are known today were not developed solely by Navier and Stokes; indeed, Poisson, Cauchy and others were also heavily involved in their development (Darrigol Reference Darrigol2002). As far as the present authors are aware, to date there is still no firm answer to the question of whether or not there always exist unique, smooth, non-singular solutions to the three-dimensional Navier–Stokes equations (Lemarie-Rieusset Reference Lemarie-Rieusset2018), and this constitutes one of the most famous unsolved problems in mathematics.

The application of analytical mechanics (Goldstein Reference Goldstein1980; Arnold Reference Arnold1989; Fetter & Walecka Reference Fetter and Walecka2003; Gelfand & Fomin Reference Gelfand and Fomin2012) to the field of fluid mechanics (Lanczos Reference Lanczos1970) has recently seen a resurgence in interest (Salmon Reference Salmon1983, Reference Salmon1988; Brenier Reference Brenier2017; Giga, Kirshtein & Liu Reference Giga, Kirshtein and Liu2018; Mottaghi, Gabbai & Benaroya Reference Mottaghi, Gabbai and Benaroya2019; Taroco, Blanco & Feijoo Reference Taroco, Blanco and Feijoo2020; Bedford Reference Bedford2021; Mavroeidis & Athanassoulis Reference Mavroeidis and Athanassoulis2022) after a long history. In the absence of non-conservative forces, an inviscid fluid is a Hamiltonian system, and so the classical Hamiltonian theory applies. Serrin (Reference Serrin1959), Benjamin (Reference Benjamin1984) and Holm, Marsden & Ratiu (Reference Holm, Marsden and Ratiu1986) have all described variational and Hamiltonian formulations of incompressible, inviscid fluid flow. Roberts (Reference Roberts1972) presented a Hamiltonian dynamic for weakly interacting vortices. This research obtained the canonical equations of Hamiltonian dynamics for a set of two well-separated vortex rings by setting up a Hamiltonian to define the set. Olver (Reference Olver1982) showed that the Euler equations of inviscid and incompressible fluid flow can be put into Hamiltonian form. Benjamin & Olver (Reference Benjamin and Olver1982) investigated the Hamiltonian structure of the water waves problem. They examined the symmetry groups of this problem, finding that Hamiltonian analysis enables the solution of conservative elements of the problem. However, the study also acknowledged that further study is needed to identify the physical significance of the mathematical results. Maddocks & Pego (Reference Maddocks and Pego1995) presented a novel Hamiltonian formulation of ideal fluid flow expressed in material coordinates. Their Hamiltonian formulation arises from a general approach for constrained systems that is not restricted to problems in fluid mechanics. Rather, it is widely applicable for obtaining unconstrained Hamiltonian dynamical systems from Lagrangian field equations that are subject to pointwise constraints. More recently, Arnold (Reference Arnold2014) also studied the Hamiltonian nature of the ideal Euler equations.

Viscous forces are non-conservative, which presents a fundamental challenge when applying Hamilton’s principle to viscous fluids (Millikan Reference Millikan1929; Finlayson Reference Finlayson1972*a*,Reference Finlayson*b*; Lemarie-Rieusset Reference Lemarie-Rieusset2018). Indeed, it is a well-known theorem (first proven by Millikan Reference Millikan1929) that the Navier–Stokes equations in their usual form cannot be derived from a classical action principle (Millikan Reference Millikan1929; Finlayson Reference Finlayson1972*a*,Reference Finlayson*b*). Millikan (Reference Millikan1929) summarizes his main result as follows:

It is impossible to derive the equations of steady motion of a viscous, incompressible fluid from a variation principle involving as Lagrangian function an expression in the velocity components and their first-order space derivatives, unless conditions are imposed on these velocity components such that all of the terms $vu_{,2}$, $wu_{,3}$, $wv_{,3}$, $uv_{,1}$, $uw_{,1}$, $vw_{,2}$ disappear from their positions in the Navier–Stokes equations (Millikan Reference Millikan1929).

(It should be noted that the six terms referred to above come from the convective acceleration $u_{i,j}u_{j}$, and Millikan (Reference Millikan1929) uses the notation $u=u_{1}$, $v=u_{2}$ and $w=u_{3}$.) More generally, it has been shown that the existence of variational formulations is related to self-adjointness of the system with respect to a standard duality relation, a property that all non-conservative systems lack (Vainberg Reference Vainberg1964). Within the last 80 years, many alternative methods have been developed in an attempt to circumvent the non-self-adjointness of dissipative systems (Prigogine & Glansdorff Reference Prigogine and Glansdorff1965; Biot Reference Biot1970; Finlayson Reference Finlayson1972*a*; Lebon & Lambermont Reference Lebon and Lambermont1973; Tonti Reference Tonti1973; Magri Reference Magri1974; Telega Reference Telega1979; Tonti Reference Tonti1984; Filippov Reference Filippov1989; Sieniutycz Reference Sieniutycz2000; Robinson Reference Robinson2001; Galley Reference Galley2013; Kim, Dargush & Lee Reference Kim, Dargush and Lee2016; Mottaghi *et al.* Reference Mottaghi, Gabbai and Benaroya2019; Taroco *et al.* Reference Taroco, Blanco and Feijoo2020; Bersani & Caressa Reference Bersani and Caressa2021; Junker & Balzani Reference Junker and Balzani2021). The mathematical study of alternative variational methods as applied to the Navier–Stokes equations in particular remains an ongoing endeavour (Oseledets Reference Oseledets1989; Vujanovic & Jones Reference Vujanovic and Jones1989; Doering & Gibbon Reference Doering and Gibbon1995; Fukagawa & Fujitani Reference Fukagawa and Fujitani2012; Jones Reference Jones2015; Gay-Balmaz & Yoshimura Reference Gay-Balmaz and Yoshimura2017; Hieber, Robinson & Shibata Reference Hieber, Robinson and Shibata2017; Hochgerner Reference Hochgerner2018; Gay-Balmaz & Yoshimura Reference Gay-Balmaz and Yoshimura2019*a*,Reference Gay-Balmaz and Yoshimura*b*; Rashad *et al.* Reference Rashad, Califano, Brugnoli, Schuller and Stramigioli2021; Gonzalez & Taha Reference Gonzalez and Taha2022; Taha & Gonzalez Reference Taha and Gonzalez2022; Sanders Reference Sanders2023*b*).

Oseledets (Reference Oseledets1989) attempted to express the Navier–Stokes equations using Hamiltonian formalism. He was able to formalize the incompressible Euler equation but stated that his result is not valid for a compressible fluid. More recent attempts, such as Fukagawa & Fujitani (Reference Fukagawa and Fujitani2012), Jones (Reference Jones2015) and Gay-Balmaz & Yoshimura (Reference Gay-Balmaz and Yoshimura2017, Reference Gay-Balmaz and Yoshimura2019*a*,Reference Gay-Balmaz and Yoshimura*b*), have enforced dissipation using a non-holonomic constraint on the entropy. Hochgerner (Reference Hochgerner2018) attempted to obtain a Hamiltonian interacting particle system that could accurately model the fluid dynamics. His research separated the dynamics into slow (deterministic) and fast (stochastic) components to capture fine-scale effects. The study was able to derive the Navier–Stokes equation from a stochastic Hamiltonian system but ignored the stress tensor, was unable to separate configuration and momentum variables and did not establish energy conservation or dissipation.

Rashad *et al.* (Reference Rashad, Califano, Brugnoli, Schuller and Stramigioli2021) modelled the incompressible Navier–Stokes equations in so-called ‘port-Hamiltonian’ framework rather than the standard Hamiltonian framework. Their model used vector calculus instead of exterior calculus to minimize the number of operators. While the main goal of this research was increasing the interest of computational researchers in using vector calculus, they also demonstrated that vector calculus can help in the formulation of individual subsystems of Navier–Stokes equations and boundary ports of the model.

Gonzalez & Taha (Reference Gonzalez and Taha2022), Taha & Gonzalez (Reference Taha and Gonzalez2022) and Taha, Gonzalez & Shorbagy (Reference Taha, Gonzalez and Shorbagy2023) have recently applied Gauss’s principle of least constraint (Gauss Reference Gauss1829) to the Navier–Stokes problem. Using Gauss’s principle, Taha *et al.* (Reference Taha, Gonzalez and Shorbagy2023) have shown that, for an incompressible fluid, the magnitude of the pressure gradient is minimum over the domain, which they term the principle of minimum pressure gradient (PMPG). When applied to an inviscid fluid in two dimensions, the PMPG provides a closure condition based in first principles that could be considered a generalization of the Kutta condition to smooth geometries. It should be noted that Gauss’s principle (Gauss Reference Gauss1829) is fundamentally different from Hamilton’s principle (Hamilton Reference Hamilton1833, Reference Hamilton1834, Reference Hamilton1835). Whereas the Hamiltonian framework involves an invariant action integral and employs variations in the coordinates (or, in continuum mechanics, the field quantities), Gauss’s principle employs variations in the accelerations. As a result, the framework of Gauss’s principle does not lead to canonical transformations.

Particularly relevant to the present work, Sanders (Reference Sanders2021, Reference Sanders2022, Reference Sanders2023*a*,Reference Sanders*b*) has shown that the higher-order dynamics are ‘intrinsically variational’, in the sense that higher-derivative versions of the classical equations of motion can be derived from a minimum-action principle even for dissipative systems, thus allowing inherently non-Hamiltonian problems to be treated as though they are Hamiltonian. This discovery has already led to two applications: the direct modal analysis of damped dynamical systems (Sanders Reference Sanders2022) and, subsequently, a new and more efficient algorithm for computing a damped system’s resonant frequencies (Sanders & Inman Reference Sanders and Inman2023). Higher-derivative theories had been studied before in the realm of quantum gravity physics (Pais & Uhlenbeck Reference Pais and Uhlenbeck1950; Van den Berg & VanderVorst Reference Van den Berg and VanderVorst2002; Kalies & VanderVorst Reference Kalies and VanderVorst2004; Bender & Mannheim Reference Bender and Mannheim2008; Smilga Reference Smilga2009; Mostafazadeh Reference Mostafazadeh2010; Baleanu *et al.* Reference Baleanu, Petras, Asad and Velasco2012; Chen *et al.* Reference Chen, Fasiello, Lim and Tolley2013) but until now they have not been applied to classical fluids. While the Navier–Stokes equations, in their standard form, may be unsuited to Hamiltonian formalism (Millikan Reference Millikan1929; Finlayson Reference Finlayson1972*a*,Reference Finlayson*b*; Doering & Gibbon Reference Doering and Gibbon1995; Hieber *et al.* Reference Hieber, Robinson and Shibata2017; Lemarie-Rieusset Reference Lemarie-Rieusset2018), it will be shown here that the higher-order dynamics can be used to restate the problem in a form consistent with Hamiltonian and Hamilton–Jacobi formalism.

In conclusion, although the body of research surrounding the Navier–Stokes equations is extensive, it would appear that no canonical Hamiltonian formulation of the Navier–Stokes problem has been found to date. That is what the present work aims to achieve.

## 3. Lagrangian formulation of the problem

Although we are primarily interested in the incompressible form of the equations given by (1.16) and (1.17), here we will begin with the compressible form of the equations, with the understanding that we will eventually take the incompressible limit. For the compressible case, the linear momentum balance and continuity equations are given by

(3.1)$$begin{gather} mathcal{R}_{i}[u_{j},p,rho;x_{k},t]equiv rho dot{u}_{i} + rho u_{i,j}u_{j} + p_{,i} – mu u_{i,jj} – (mu+lambda)u_{j,ji} – rho b_{i} = 0, end{gather}$$

(3.2)$$begin{gather}mathcal{R}_{4}[u_{j},rho]equiv dot{rho} + rho_{,i}u_{i} + rho u_{i,i} = 0, end{gather}$$

where $rho =rho (x_{j},t)$ is the density field (now one of the unknown field quantities along with $u_{i}$ and $p$), and $lambda$ is an additional viscosity coefficient which, under Stokes’s (Stokes Reference Stokes1845) hypothesis, is related to $mu$ as $lambda =-2mu /3$, ensuring that the mechanical pressure agrees with the thermodynamic pressure. Henceforth, we will assume that all quantities have been suitably non-dimensionalized. The non-dimensional (constant) viscosities in (3.1) and (3.2) may be regarded as inverse Reynolds numbers, and the non-dimensional pressure may be considered to be normalized by the inertial scale $rho _0 U^2$, with $rho _0$ and $U$ appropriate density and velocity scales. As we will see later, starting from the compressible form of the equations will allow us to treat the pressure as a dynamical field variable alongside the velocities, rather than simply a Lagrange multiplier. Crucially, this will reveal in no uncertain terms what becomes of the momentum conjugate to the pressure (which will be identified later) in the incompressible limit.

In general, (3.1) and (3.2) would be appended with the energy equation, which introduces additional thermodynamic variables, such as temperature and enthalpy or entropy. Two of the thermodynamic variables are designated as ‘primary’, and equations of state are required to relate the remaining variables to the primary variables. Typically, pressure and temperature are chosen as the primary variables, and the equation of state for the density, for example, is expressed as $rho =rho (p,T)$. The conservation equations along with the equation of state constitute six equations for the six unknowns fields $(u_i,p,T,rho )$. Henceforth in the present work, we will take the temperature to be constant, though we intend to consider variations in temperature in future work.

An incompressible flow is one for which the material derivative of the density vanishes, i.e. $text {d}rho /text {d}t=dot {rho } + rho _{,i}u_{i}=0$, and this condition serves as an equation of state. It is usually also assumed, for the sake of simplicity, that the density is both constant and uniform, further reducing the equation of state $rho =rho (p,T)$ to specification of $rho =rho _0$ as a system parameter. Consequently, (3.2) reduces to $rho u_{i,i}=0$ and the energy equation is decoupled from the system. Accordingly, in the incompressible limit there are only four unknown field quantities $(u_{i},p)$ and the momentum balance and continuity equations suffice for the governing field equations.

We pause here to remark that all four field equations (3.1), (3.2) are first order in time with respect to the field quantities $u_{i}$ and $rho$. This will be important shortly, when we double the order of the equations. It should also be noted, as mentioned previously, that the first-order problem described above is inherently non-Hamiltonian, in that there is no action $mathcal {S}$ for which Hamilton’s principle ($delta mathcal {S}=0$) yields the first-order field equations (Millikan Reference Millikan1929; Finlayson Reference Finlayson1972*a*,Reference Finlayson*b*). Finally, we note that in the incompressible limit, $mathcal {R}_{4}$ becomes independent of $dot {rho }$ and is no longer first order in time.

### 3.1. Second-order formulation

Although the first-order formulation of the problem is intrinsically non-Hamiltonian (Millikan Reference Millikan1929; Finlayson Reference Finlayson1972*a*,Reference Finlayson*b*), nevertheless a Hamiltonian for the system may be found by considering a second-order formulation. Following Sanders (Reference Sanders2023*b*), we observe that the actual motion of the fluid corresponds to the particular fields $(u_{i},p,rho )$ for which the following action achieves a local minimum:

(3.3)begin{equation} mathcal{S}^{*}[u_{j},p,rho]=inttext{d}^{4}x(tfrac{1}{2}mathcal{R}_{i}mathcal{R}_{i}+tfrac{1}{2}mathcal{R}_{4}mathcal{R}_{4}), end{equation}

where $text {d}^{4}x=text {d}kern 0.06em x_{1},text {d}kern 0.06em x_{2},text {d}kern 0.06em x_{3},text {d}t$, and the integral is carried out over both the control volume $mathcal {V}$ occupied by the fluid ($x_{j}in mathcal {V}$) and the time period of interest ($tin [t_{1},t_{2}]$). It must be emphasized that this action contains no new physics. Again, this is simply the principle of least squares (Finlayson Reference Finlayson1972*b*) averaged over the space–time occupied by the fluid. The entire physics of the problem are already contained in the residuals ($mathcal {R}_{i}$, $mathcal {R}_{4}$).

Without an equation of state relating $rho$ to $p$, the problem is underconstrained with five unknown field quantities and only four dynamical field equations. Anticipating the case of incompressible flow, where the density is constant and the four field quantities are $u_{i}$ and $p$, henceforth we will assume an equation of state of the form $rho =hat {rho }(p)$, with $hat {rho }$ a known function determined either from first principles or empirically. In this way, the density field may be eliminated in favour of the pressure field, and the field equations assume the following form:

(3.4)$$begin{gather} mathcal{R}_{i}[u_{j},p;x_{k},t]equiv hat{rho} dot{u}_{i} + hat{rho} u_{i,j}u_{j} + p_{,i} – mu u_{i,jj} – (mu+lambda)u_{j,ji} – hat{rho} b_{i} = 0, end{gather}$$

(3.5)$$begin{gather}mathcal{R}_{4}[u_{j},p]equiv hat{rho}’dot{p} + hat{rho}’p_{,i}u_{i} + hat{rho} u_{i,i} = 0, end{gather}$$

where $hat {rho }’=text {d}hat {rho }/text {d}p$. We note that, under equilibrium conditions, $hat {rho }’$ is related to the speed of sound $c$ and the bulk modulus $K$ as $hat {rho }’=1/c^{2}=rho /K$ (for incompressible fluids, $hat {rho }’equiv 0$ and the speed of sound and bulk modulus are both infinite). Having specified $hat {rho }(p)$, $mu$, $lambda$ and $b_{i}(x_{j},t)$, and having prescribed appropriate auxiliary conditions (initial and boundary conditions), one seeks the four field quantities $(u_{i},p)$ satisfying the governing field equations and the auxiliary conditions. To recover the case of incompressible flow, we will eventually take $hat {rho }’equiv 0$.

We pause here to note that, even though the residuals ($mathcal {R}_{i}$, $mathcal {R}_{4}$) vanish for the actual motion, they are not trivially zero. That is, the residuals only vanish for the particular fields $(u_{i},p)$ which satisfy the first-order field equations (3.4) and (3.5); they do not vanish for every conceivable $(u_{i},p)$. Thus it is not appropriate to take $mathcal {R}_{i}equiv 0$, $mathcal {R}_{4}equiv 0$. We will return to this point later when we discuss the Hamiltonian formulation of the problem.

For now, we note that the action $mathcal {S}^{*}=mathcal {S}^{*}[u_{i},p]$ defines a Lagrangian

(3.6)begin{equation} L^{*}[u_{i},p;t]=int text{d}^{3}x(mathcal{L}^{*}) , end{equation}

where the integral is carried out over the volume $mathcal {V}$ only ($text {d}^{3}x=text {d}kern 0.06em x_{1},text {d}kern 0.06em x_{2},text {d}kern 0.06em x_{3}$), with Lagrangian density

(3.7)begin{equation} mathcal{L}^{*}[u_{j},p;x_{k},t]=tfrac{1}{2}mathcal{R}_{i}mathcal{R}_{i}+tfrac{1}{2}mathcal{R}_{4}mathcal{R}_{4}. end{equation}

Because the residuals ($mathcal {R}_{i}, mathcal {R}_{4}$) have been non-dimensionalized, the Lagrangian density is also dimensionless. Once again, even though the Lagrangian vanishes for the actual motion, it is not trivially zero, and it is not appropriate to take $L^{*}equiv 0$.

As noted above, the actual motion of the fluid corresponds to the particular fields $(u_{i},p)$ for which $mathcal {S}^{*}$ achieves a local minimum. To obtain the Euler–Lagrange equations, the conjugate momenta and the natural auxiliary conditions, we insist that $mathcal {S}^{*}$ not vary to first order ($delta mathcal {S}^{*}=0$) under small variations in the fields $(delta u_{i},delta p)$. Evaluating $delta mathcal {S}^{*}$, integrating by parts, and collecting like terms, we find that

(3.8*a*)begin{align} deltamathcal{S}^{*}&=inttext{d}^{4}xleft{left[-frac{partial}{partial t} left(hat{rho}mathcal{R}_{i}right) -frac{partial}{partial x_{j}} left(hat{rho} mathcal{R}_{i}u_{j}right)+ hat{rho} mathcal{R}_{j}u_{j,i}- mu mathcal{R}_{i,jj} right.right. end{align}

(3.8*b*)begin{align} &quad – left.(mu+lambda)mathcal{R}_{j,ij}+ hat{rho}’mathcal{R}_{4}p_{,i}-frac{partial}{partial x_{i}}left( hat{rho} mathcal{R}_{4}right)right]delta{u}_{i} end{align}

(3.8*c*)begin{align} &quad +left[vphantom{frac{partial}{partial t}}hat{rho}’mathcal{R}_{i} dot{u}_{i}+ hat{rho}’ mathcal{R}_{i}u_{i,j}u_{j}- mathcal{R}_{i,i}- hat{rho}’ mathcal{R}_{i}b_{i} + hat{rho}”mathcal{R}_{4}dot{p}right. end{align}

(3.8*d*)begin{align} &quad -left.left.frac{partial}{partial t}left( hat{rho}’mathcal{R}_{4}right) + hat{rho}”mathcal{R}_{4}p_{,i}u_{i} – frac{partial}{partial x_{i}}left( hat{rho}’mathcal{R}_{4}u_{i}right)+ hat{rho}’ mathcal{R}_{4}u_{i,i}right] delta p right} end{align}

(3.8*e*)begin{align} &quad +inttext{d}^{3}xleft[hat{rho}mathcal{R}_{i} delta{u}_{i}+hat{rho}’mathcal{R}_{4}delta{p}right]_{t_{1}}^ {t_{2}} end{align}

(3.8*f*)begin{align} &quad +inttext{d}^{2}x,text{d}tleft{left[hat{rho} mathcal{R}_{i}u_{j}n_{j}+mu mathcal{R}_{i,j}n_{j}+ (mu+lambda)mathcal{R}_{j,i}n_{j}+hat{rho} mathcal{R}_{4}n_{i}right]delta u_{i}right. end{align}

(3.8*g*)begin{align} &quad +left.left[- mu mathcal{R}_{i}n_{j}- (mu+lambda)mathcal{R}_{j}n_{i}right]delta u_{i,j}+left[mathcal{R}_{i}n_{i}+hat{rho}’mathcal{R}_{4}u_{i}n_{i}right]delta pright}, end{align}

where the purely volumetric integral ($text {d}^{3}x$) is carried out over $mathcal {V}$, the surface integral ($text {d}^{2}x$) is carried out over the boundary $partial mathcal {V}$ and $n_{i}$ is the unit outward normal vector to $partial mathcal {V}$. Note that, because we are using Eulerian coordinates $x_{j}$, the volume element $text {d}^{3}x$ is not to be varied.

The Euler–Lagrange equations (which hold for all $x_{j}in mathcal {V}$) may be read directly from the space–time ($text {d}^{4}x$) integral

(3.9)begin{align} delta{u}_{i}: quad &{-}frac{partial}{partial t} left(hat{rho}mathcal{R}_{i}right) -frac{partial}{partial x_{j}} left(hat{rho} mathcal{R}_{i}u_{j}right)+ hat{rho} mathcal{R}_{j}u_{j,i}- mu mathcal{R}_{i,jj} nonumber\ &- (mu+lambda)mathcal{R}_{j,ij}+ hat{rho}’mathcal{R}_{4}p_{,i}-frac{partial}{partial x_{i}}left( hat{rho} mathcal{R}_{4}right)=0, end{align}

(3.10)begin{align} delta p:quad & hat{rho}’mathcal{R}_{i} dot{u}_{i}+ hat{rho}’ mathcal{R}_{i}u_{i,j}u_{j}- mathcal{R}_{i,i}- hat{rho}’ mathcal{R}_{i}b_{i} + hat{rho}”mathcal{R}_{4}dot{p} nonumber\ &-frac{partial}{partial t}left( hat{rho}’mathcal{R}_{4}right) + hat{rho}”mathcal{R}_{4}p_{,i}u_{i} – frac{partial}{partial x_{i}}left(hat{rho}’mathcal{R}_{4}u_{i}right)+ hat{rho}’ mathcal{R}_{4}u_{i,i} =0. end{align}

It should be noted that all four Euler–Lagrange equations (3.9), (3.10) are second order in time, as they involve time derivatives of the residuals. By doubling the order of the equations, we have put the problem in Hamiltonian form, consistent with the general result of Sanders (Reference Sanders2023*b*). We also note that all four Euler–Lagrange equations of the second-order formulation are automatically satisfied by the solution to the first-order formulation (i.e. the actual motion), for which $mathcal {R}_{i}=0$ and $mathcal {R}_{4}=0$ everywhere and at all times.

Corresponding to each of the four field quantities is a canonically conjugate ‘momentum’ field, which can be read from (3.8*e*). The momenta conjugate to the velocities $u_{i}$ are

(3.11)begin{equation} {rm pi}_{i}equiv hat{rho}mathcal{R}_{i}, end{equation}

and the momentum conjugate to the pressure $p$ is

(3.12)begin{equation} {rm pi}_{4}equiv hat{rho}’mathcal{R}_{4}. end{equation}

In the forthcoming Hamiltonian formulation, the conjugate momenta will be used to eliminate the (partial) time derivatives $(dot {u}_{i},dot {p})$ of the field quantities from the Hamiltonian. In general, Hamilton’s principle would insist that the variations $(delta u_{i},delta p)$ vanish at the endpoints $t=t_{1}$ and $t=t_{2}$ to ensure that the purely volumetric integral (3.8*e*) vanishes identically. Interestingly, for the actual motion ($mathcal {R}_{i}=0,mathcal {R}_{4}=0$), the volumetric integral (3.8*e*) already vanishes even without taking $(delta u_{i},delta p)$ to vanish at $t_{1}$ and $t_{2}$. We interpret this to mean that the actual motion is the natural evolution of the second-order formulation (Sanders Reference Sanders2023*b*).

Although the conjugate momenta (${rm pi} _{i}$, ${rm pi} _{4}$) do not coincide with conventional linear or angular momenta, there is nonetheless a curious mathematical connection between the conjugate momenta and the linear momentum density $P_{i}=rho u_{i}$, which we will see presently from the natural boundary conditions. These are read directly from the surface ($text {d}^{2}x$) integral

(3.13)$$begin{gather} delta u_{i}:quad hat{rho} mathcal{R}_{i}u_{j}n_{j}+mu mathcal{R}_{i,j}n_{j}+ (mu+lambda)mathcal{R}_{j,i}n_{j}+hat{rho} mathcal{R}_{4}n_{i}=0, end{gather}$$

(3.14)$$begin{gather}delta u_{i,j}:quad – mu mathcal{R}_{i}n_{j}- (mu+lambda)mathcal{R}_{j}n_{i}=0, end{gather}$$

(3.15)$$begin{gather}delta p:quadmathcal{R}_{i}n_{i}+hat{rho}’mathcal{R}_{4}u_{i}n_{i}=0. end{gather}$$

This last condition, (3.15), establishes a connection between the new conjugate momenta and the conventional linear momenta. Multiplying (3.15) by $hat {rho }$, and noting that $hat {rho }u_{i}=rho u_{i}= P_{i}$, we find that

(3.16)begin{equation} ({rm pi}_{i}+{rm pi}_{4}P_{i})n_{i}=0. end{equation}

Evidently, boundary condition (3.15) states that the flux of the vector $varPi _{i}equiv {rm pi}_{i}+{rm pi} _{4}P_{i}$ through the boundary $partial mathcal {V}$ should vanish. It is interesting that this new vector $varPi _{i}$ contains both old and new momenta, with ${rm pi} _{4}$ ‘carried’ (i.e. given direction) by $P_{i}$. The actual physical meaning of these natural boundary conditions is less clear and may require further investigation.

### 3.2. Equivalence of the first- and second-order formulations

The first- and second-order formulations are mathematically equivalent, in the sense that imposing identical auxiliary conditions on the two formulations will yield identical solutions $(u_{i},p)$. In other words, with identical auxiliary conditions, $(u_{i},p)$ is a solution to the first-order formulation if and only if the same $(u_{i},p)$ is a solution to the second-order formulation.

The proof is straightforward. Consider the two formulations separately, and impose on the second-order formulation identical auxiliary conditions to those of the first-order formulation. In particular, just like the simple example given in § 1.1, the second-order formulation requires additional auxiliary conditions over and above those applied to the first-order formulation. These include initial conditions making $mathcal {R}_{i}(x_{k},0)=0$ and $mathcal {R}_{4}(x_{k},0)=0$ for all $x_{k}in mathcal {V}cup partial mathcal {V}$, along with boundary conditions making $mathcal {R}_{i}(x_{k},t)=0$, $mathcal {R}_{i,j}(x_{k},t)=0$ and $mathcal {R}_{4}(x_{k},t)=0$ for all $x_{k}in partial mathcal {V}$ and all times $t$. By supposition, the auxiliary conditions applied to the two formulations are identical, so it suffices to show that $(u_{i},p)$ satisfies the governing field equations (3.4), (3.5) of the first-order formulation ($mathcal {R}_{i}=0$ and $mathcal {R}_{4}=0$) everywhere in $mathcal {V}$ and at all times if and only if $(u_{i},p)$ satisfies the Euler–Lagrange equations (3.9), (3.10) of the second-order formulation everywhere in $mathcal {V}$ and at all times.

Suppose first that $(u_{i},p)$ satisfies the governing field equations (3.4), (3.5) of the first-order formulation everywhere in $mathcal {V}$ and at all times. Then $mathcal {R}_{i}=0$, $mathcal {R}_{4}=0$, and $(u_{i},p)$ is a trivial solution to the Euler–Lagrange equations (3.9), (3.10) of the second-order formulation.

Conversely, suppose that $(u_{i},p)$ satisfies the Euler–Lagrange equations (3.9), (3.10) of the second-order formulation everywhere in $mathcal {V}$ and at all times. We note that $(mathcal {R}_{i}=0,mathcal {R}_{4}=0)$ constitutes an equilibrium solution of the Euler–Lagrange equations (3.9), (3.10). Thus, because the initial conditions are chosen such that $mathcal {R}_{i}(x_{k},0)=0$ and $mathcal {R}_{4}(x_{k},0)=0$ for all $x_{k}in mathcal {V}cup partial mathcal {V}$, and because the boundary conditions are chosen such that $mathcal {R}_{i}(x_{k},t)=0$, $mathcal {R}_{i,j}(x_{k},t)=0$ and $mathcal {R}_{4}(x_{k},t)=0$ for all $x_{k}in partial mathcal {V}$ and all times $t$, then $mathcal {R}_{i}$ and $mathcal {R}_{4}$ will remain identically zero everywhere in $mathcal {V}$ for all future times. Thus, $(u_{i},p)$ satisfies the governing field equations (3.4), (3.5) of the first-order formulation everywhere in $mathcal {V}$ and at all times. This completes the proof, and we have established that the two formulations are equivalent.

## 4. Hamiltonian formulation of the problem

We are now ready to proceed with the Hamiltonian formulation of the problem. The Lagrangian $L^{*}$ has a corresponding Hamiltonian

(4.1)begin{equation} H^{*}=int text{d}^{3}x(mathcal{H}^{*}) , end{equation}

with the Hamiltonian density $mathcal {H}^{*}$ obtained from the Lagrangian density $mathcal {L}^{*}$ via the Legendre transform

(4.2)begin{equation} mathcal{H}^{*}={rm pi}_{i}dot{u}_{i}+{rm pi}_{4}dot{p}-mathcal{L}^{*}={rm pi}_{i}dot{u}_{i}+{rm pi}_{4}dot{p}-tfrac{1}{2}mathcal{R}_{i}mathcal{R}_{i}-tfrac{1}{2}mathcal{R}_{4}mathcal{R}_{4}. end{equation}

Again, this $H^{*}$ has nothing to do with the total mechanical energy of the system, although it is a conserved quantity, since $H^{*}=0$ for the actual motion – just as in the example of § 1.1. In order to write down Hamilton’s equations, we must express $mathcal {H}^{*}$ in terms of the fields and the conjugate momenta, eliminating $dot {u}_{i}$ and $dot {p}$.

Observe that $mathcal {R}_{i}={rm pi} _{i}/hat {rho }$, and ignoring for the moment the incompressible limit, we may write $mathcal {R}_{4}={rm pi} _{4}/hat {rho }’$ ($hat {rho }’neq 0$). In this way, using the functional expressions for the residuals given by (3.4) and (3.5), we find that

(4.3)$$begin{gather} dot{u}_{i}=frac{1}{(hat{rho})^{2}}{{rm pi}_{i}}-frac{1}{hat{rho}}left(hat{rho}u_{i,j}u_{j} +{ p_{,i} – mu u_{i,jj} – (mu+lambda)u_{j,ji}} – hat{rho}b_{i}right); end{gather}$$

(4.4)$$begin{gather}dot{p}=frac{1}{(hat{rho}’)^{2}}{{rm pi}_{4}}-frac{1}{hat{rho}’}left({hat{rho}’p_{,i}u_{i} + hat{rho}u_{i,i}}right), quad hat{rho}’neq0; end{gather}$$

and

(4.5)begin{align} mathcal{H}^{*}[u_{i},p,{rm pi}_{j},{rm pi}_{4};x_{k},t]&=frac{1}{2}frac{1}{(hat{rho})^{2}}{{rm pi}_{i}{rm pi}_{i}}-frac{1}{hat{rho}}(hat{rho}u_{i,j}u_{j} +{p_{,i} – mu u_{i,jj} – (mu+lambda)u_{j,ji}} – hat{rho}b_{i}){rm pi}_{i}nonumber\ &quad +frac{1}{2}frac{1}{(hat{rho}’)^{2}}{{rm pi}_{4}{rm pi}_{4}}-frac{1}{hat{rho}’}({hat{rho}’p_{,i}u_{i} + hat{rho}u_{i,i}}){rm pi}_{4}, quad hat{rho}’neq0. end{align}

Hamilton’s equations (Hamilton Reference Hamilton1834, Reference Hamilton1835) are as follows:

(4.6*a*,*b*)$$begin{gather} dot{u}_{i}=frac{delta H^{*}}{delta{rm pi}_{i}}, quaddot{p}=frac{delta H^{*}}{delta{rm pi}_{4}}, end{gather}$$

(4.7*a*,*b*)$$begin{gather}dot{rm pi}_{i}={-}frac{delta H^{*}}{delta u_{i}}, quad dot{rm pi}_{4}={-}frac{delta H^{*}}{delta p}, end{gather}$$

where $delta H^{*}/delta u_{i}$, $delta H^{*}/delta p$, $delta H^{*}/delta {rm pi}_{i}$ and $delta H^{*}/delta {rm pi}_{4}$ are the Volterra (Volterra Reference Volterra1930) functional derivatives of $H^{*}$ with respect to the field quantities and the conjugate momenta. Equations (4.6*a*,*b*) reproduce (4.3) and (4.4), respectively. Equations (4.7*a*,*b*) in turn reproduce the Euler–Lagrange equations (3.9), (3.10) of the second-order formulation.

We return now to our previous observation concerning the vanishing of the residuals. While $H^{*}$ vanishes for the particular fields $(u_{i},p)$ that satisfy the governing field equations (3.4), (3.5) of the first-order formulation, it does not vanish for every conceivable $(u_{i},p)$. The latter would imply, according to (4.6*a*,*b*), that $dot {u}_{i}equiv 0$ and $dot {p}equiv 0$, which is not generally the case. This observation, and the fact that (4.7*a*,*b*) faithfully reproduce the Euler–Lagrange equations (3.9), (3.10) of the second-order formulation, confirm that the Hamiltonian formulation described above is, in fact, a legitimate reformulation of the problem. In the following section, we develop the Hamilton–Jacobi theory as it relates to the present formulation, the goal being to find a canonical transformation to a new set of fields ($phi _{i}, phi _{4}$) and conjugate momenta for which the Hamiltonian does vanish identically.

To obtain the Hamiltonian for incompressible flow, we set $hat {rho }’equiv 0$ from the beginning, in which case $mathcal {R}_{4}$ reduces to $hat {rho } u_{i,i}$ and ${rm pi} _{4}$ vanishes identically, consistent with the fact that $mathcal {R}_{4}$ becomes independent of $dot {p}$. The Hamiltonian density in turn reduces to

(4.8)begin{equation} mathcal{H}^{*}={rm pi}_{i}dot{u}_{i}-tfrac{1}{2}mathcal{R}_{i}mathcal{R}_{i}, quad hat{rho}’equiv0, end{equation}

or, in terms of the conjugate momenta,

(4.9)begin{equation} mathcal{H}^{*}[u_{i},p,{rm pi}_{j};x_{k},t]=frac{1}{2}frac{1}{{rho}^{2}}{{rm pi}_{i}{rm pi}_{i}}-frac{1}{rho}({rho}u_{i,j}u_{j} + p_{,i} – mu u_{i,jj} – {rho}b_{i}){rm pi}_{i}, quad hat{rho}’equiv0,end{equation}

where the density $hat {rho }=rho$ is a constant and we have used the fact that $u_{i,i}=0$. Hamilton’s equations $dot {u}_{i}=delta H^{*}/delta {rm pi}_{i}$, $dot {{rm pi} }_{i}=-delta H^{*}/delta u_{i}$ and $0equiv dot {{rm pi} }_{4}=-delta H^{*}/delta p$ still apply (and reproduce the corresponding equations in the incompressible limit), but $dot {p}=delta H^{*}/delta {rm pi}_{4}$ must be replaced by the constraint that $u_{i,i}=0$. That the incompressibility condition should take the place of the pressure equation $dot {p}=delta H^{*}/delta {rm pi}_{4}$ is consistent with the well-known result that the pressure usually serves as the Lagrange multiplier for the incompressibility constraint (Lanczos Reference Lanczos1970; Badin & Crisciani Reference Badin and Crisciani2018).

### 4.1. Hamilton–Jacobi equation

One of the most significant aspects of the Hamiltonian formalism is that it leads to the transformation theory of Hamilton (Reference Hamilton1833, Reference Hamilton1834, Reference Hamilton1835) and Jacobi (Reference Jacobi1837, Reference Jacobi1842–1843) (see Whittaker (Reference Whittaker1904) and Lanczos (Reference Lanczos1970)), celebrated both for unifying particle mechanics with wave optics (Hamilton Reference Hamilton1833) and for its relationship to the Schrödinger equation of quantum mechanics (Schrödinger Reference Schrödinger1926*a*,Reference Schrödinger*b*). Here, we will obtain a Hamilton–Jacobi equation representing the Navier–Stokes problem.

In the context of discrete mechanics, Hamilton’s principal function is obtained as the solution to the Hamilton–Jacobi equation, which is in turn defined by the functional form of the Hamiltonian. Hamilton’s principal function provides the generating function for a canonical transformation to a new set of generalized coordinates and conjugate momenta for which the Hamiltonian vanishes identically, in which case Hamilton’s equations do, in fact, become trivial. The new coordinates and their conjugate momenta are simply equal to their initial values.

In the present context, the role of Hamilton’s principal function is played by a characteristic functional ${S}^{*}={S}^{*}[u_{i},p,t]$ (not to be confused with the action $mathcal {S}^{*}$, although they are related; see Appendix A), which is the solution to the following Hamilton–Jacobi equation:

(4.10)begin{equation} {H}^{*}left[u_{i},p,frac{delta {S}^{*}}{delta u_{j}},frac{delta {S}^{*}}{delta p};tright]+frac{partial {S}^{*}}{partial t}=0, end{equation}

where $delta {S}^{*}/delta u_{i}$ and $delta {S}^{*}/delta p$ are the Volterra (Volterra Reference Volterra1930) functional derivatives of ${S}^{*}$ with respect to the field quantities. Interested readers will find the derivation of (4.10) in Appendix A. Henceforth, we will refer to ${S}^{*}$ as ‘Hamilton’s principal functional’. Substituting for the conjugate momenta in (4.5), we obtain for the Hamilton–Jacobi equation

(4.11)begin{align} &int text{d}^{3}xleft[frac{1}{2}frac{1}{(hat{rho})^{2}}frac{delta {S}^{*}}{delta u_{i}}frac{delta {S}^{*}}{delta u_{i}}-frac{1}{hat{rho}}left(hat{rho}u_{i,j}u_{j} +{ p_{,i} – mu u_{i,jj} – (mu+lambda)u_{j,ji}} – hat{rho}b_{i}right)frac{delta {S}^{*}}{delta u_{i}}right. nonumber\ &quad left.+,frac{1}{2}frac{1}{(hat{rho}’)^{2}}frac{delta {S}^{*}}{delta p}frac{delta {S}^{*}}{delta p}-frac{1}{hat{rho}’}left({hat{rho}’p_{,i}u_{i} + hat{rho}u_{i,i}}right)frac{delta {S}^{*}}{delta p}right]+frac{partial {S}^{*}}{partial t}=0, quad hat{rho}’neq0. end{align}

In contrast to the four original field equations – (3.1) and (3.2) – the Hamilton–Jacobi equation (4.11) is a single equation in Hamilton’s principal functional ${S}^{*}$. This constitutes an equivalent formulation of the problem, as a complete and non-trivial solution to (4.11) is tantamount to an integration of Hamilton’s equations (4.6*a*,*b*) and (4.7*a*,*b*) (note that it is not appropriate to take ${S}^{*}equiv 0$ for the same reason that it is not appropriate to take $H^{*}equiv 0$). In this way, we have reduced the problem of finding four separate field quantities to that of finding a single functional in those field quantities. One need only deduce (or even guess) the general form of ${S}^{*}$ in order to solve the problem. If an analytical expression for ${S}^{*}$ can be obtained, it will lead via canonical transformation to a new set of fields ($phi _{i}, phi _{4}$) and conjugate momenta which are simply equal to their initial values, giving analytical expressions for the four original fields $(u_{i},p)$.

The case of incompressible flow requires care, as $hat {rho }’equiv 0$ and $hat {rho }’$ appears in the denominators of terms in (4.11). Even so, the Hamiltonian formulation remains well posed in the incompressible limit. Recall that, with $hat {rho }’equiv 0$, $mathcal {R}_{4}$ reduces to $hat {rho } u_{i,i}$, ${rm pi} _{4}$ vanishes identically and the Hamiltonian density reduces to (4.9). Substituting for the conjugate momenta ${rm pi} _{i}$ in (4.9), the corresponding Hamilton–Jacobi equation is

(4.12)begin{align} int text{d}^{3}xleft[frac{1}{2}frac{1}{{rho}^{2}}frac{delta {S}^{*}}{delta u_{i}}frac{delta {S}^{*}}{delta u_{i}}-frac{1}{rho}({rho}u_{i,j}u_{j} + p_{,i} – mu u_{i,jj} – {rho}b_{i})frac{delta {S}^{*}}{delta u_{i}}right]+frac{partial {S}^{*}}{partial t}=0, quad hat{rho}’equiv0, end{align}

with $delta {S}^{*}/delta p=0$, since again ${rm pi} _{4}$ vanishes identically for incompressible flow. Here, the merit of starting from the compressible form of the equations becomes fully evident, as it would not necessarily have been clear that $delta {S}^{*}/delta p$ should vanish in the incompressible limit without knowing that in general ${rm pi} _{4}=hat {rho }’mathcal {R}_{4}$. This is the form of the Hamilton–Jacobi equation as it relates to the traditional Navier–Stokes problem. In this case, the pressure is determined last of all, and is whatever it needs to be to enforce the incompressibility constraint $u_{i,i}=0$ (again consistent with the role of pressure as Lagrange multiplier Lanczos Reference Lanczos1970; Badin & Crisciani Reference Badin and Crisciani2018).

It must be acknowledged that the Hamilton–Jacobi equation developed above (either (4.11) for the compressible case or (4.12) in the incompressible limit) contains Volterra (Reference Volterra1930) functional derivatives and is thus by no means trivial to solve. Indeed, it appears that solving such equations is itself a long-standing open problem in mathematics, which has received very little attention since the first half of the twentieth century (Michal Reference Michal1926; Jordan & Pauli Reference Jordan and Pauli1928; Levy Reference Levy1951; Tatarskii Reference Tatarskii1961; Syavavko & Mel’nichak Reference Syavavko and Mel’nichak1974; Dieudonne Reference Dieudonne1981; Koval’chik Reference Koval’chik1993). Nevertheless, if a rigorous theory of such equations can be developed, the present formulation of the Navier–Stokes problem might be solved as one special case. The present authors submit that such an endeavour is worthwhile and merits further study.

We conclude this section by remarking that, in the inviscid limit ($mu =lambda =0$), all of the preceding formalism remains perfectly well posed. In that limit, the present approach yields a mathematically equivalent second-order formulation of the inviscid Euler equations, as one would expect. Interested readers will find the full details in Appendix B.

## 5. Discussion

In this section we provide some qualitative interpretations of the developments of § 3. More specifically, we investigate the incompressible form (via constant, uniform density) of the Euler–Lagrange equations (3.9) and (3.10) when the residuals $mathcal {R}_i$ and $mathcal {R}_4$ are substituted.

Our motivation is again the simple example of § 1.1 for which the first-order non-Hamiltonian system $dot {v}=-v$ was converted to the second-order Hamiltonian system $ddot {v}=v$ by (manual) elimination of the non-conservative ‘damping’ term $dot {v}$ (see Sanders (Reference Sanders2022) for a similar result for the damped harmonic oscillator converting from second-order to fourth-order dynamics). Sanders (Reference Sanders2023*b*) showed that the elimination process is ‘automated’ by the definition of the action in the first integral of (1.3), which is generalized to the action in (3.3) for our current continuum dynamics problem containing fields.

First consider the pressure equation (3.10) and corresponding natural boundary condition (3.15), which take the following incompressible forms:

(5.1)begin{equation} {-}mathcal{R}_{i,i}=0 quad forall x_jinmathcal{V}, quadtext{subject to}quad mathcal{R}_in_i = 0 quad forall x_jinpartialmathcal{V}. end{equation}

This higher-order field equation is simply the divergence of the residual $mathcal {R}_i$. Upon substituting for $mathcal {R}_i$ from (3.1) and subsequently imposing the incompressible continuity condition $u_{i,i}=mathcal {R}_4/rho =0$ from (3.2), we obtain

(5.2)begin{equation} p_{,ii} ={-}[rho u_ju_{i,j}]_{,i} + rho b_{i,i}, end{equation}

which is a Poisson equation for the pressure. The boundary condition is a Neumann type requiring the specification of the normal pressure gradient, $n_ip_{,i} = p_{,n} equiv f(x_j,t)$, on the boundary, where

(5.3)begin{equation} f(x_j,t) ={-}n_i[rhodot{u}_i+rho u_ju_{i,j}-mu u_{i,jj}-rho b_i]. end{equation}

Equation (5.2) and boundary condition (5.3) evolve the pressure in a manner that ensures the velocity field is solenoidal. This is a well-known pressure–velocity-based formulation commonly used in the numerical solution of incompressible flows (e.g. Ferziger & Peric Reference Ferziger and Peric2002; Pozrikidis Reference Pozrikidis2009).

Next, we consider the velocity equations (3.9) which, at present, have a more elusive physical interpretation. Here, we instead begin with the natural boundary conditions (3.13) and (3.14), which are due to the $delta u_i$ and $delta u_{i,j}$ variations. The incompressible versions of these equations are

(5.4*a*,*b*)begin{equation} rho mathcal{R}_{i}u_{j}n_{j}+mu mathcal{R}_{i,j}n_{j}=0 quadtext{and}quad – mu mathcal{R}_{i}n_{j}=0 quad forall x_jinpartialmathcal{V}. end{equation}

The boundary conditions involving the residual $mathcal {R}_i$ are those compatible with the first-order Navier–Stokes equations, such as the no-slip and no-penetration conditions. Indeed, if we specify the velocity vector of the actual motion on the boundary, then $mathcal {R}_iequiv 0$ there. Note that the pressure of the actual motion on the boundary will be known from the simultaneous solution of (5.2).

However, the gradient terms $mathcal {R}_{i,j}$ will introduce up to third-order spatial derivatives that must be specified. These represent the additional boundary conditions that must accompany the higher-order governing equation, which will be seen shortly to be second order in time and fourth order in space. Again, recall the example of § 1, whereby the system (1.2*a*–*c*) must be appended with a second (initial) condition specifying the (time) derivative of the coordinate $v(t)$. In the present context, these boundary conditions are ostensibly tantamount to specification of the viscous stress on the boundary by way of velocity gradients.

In general, the conditions at a boundary require two transition relations (Batchelor Reference Batchelor2000; White Reference White2006) to ultimately describe the momentum transport. Mathematically speaking, these conditions are the jump in velocity (momentum intensity) and the jump in stress (momentum flux). Under ordinary physical circumstances the velocity and stress are assumed to be continuous. However, this is one particular form of the transition relations, and there are familiar examples to which they do not apply. For example, at a liquid–gas interface the stress relation is modified to account for a non-zero jump in the normal stress that is balanced by a force due to surface tension (the tangential stress component usually still taken to be continuous). Similarly, in the event that molecular slip occurs, the typical transition relation gives an expression for the slip velocity (e.g. Thompson & Troian Reference Thompson and Troian1997; Thalakkottor & Mohseni Reference Thalakkottor and Mohseni2016). In the case of energy transport, analogous conditions are needed regarding jumps in temperature (intensity) and heat flow (flux), which are recognized as the concept of thermal contact resistance.

We now turn our attention to the Euler–Lagrange equations (3.9), which upon imposing incompressibility and expanding derivatives of product terms yields

(5.5)begin{equation} rhodot{mathcal{R}}_i + rho u_jmathcal{R}_{i,j} = rho mathcal{R}_{j}u_{j,i}- mu mathcal{R}_{i,jj} quad forall x_jinmathcal{V}. end{equation}

The left-hand side is the material derivative of the residual $mathcal {R}_i$. Our purpose here is to observe which terms from the first-order Navier–Stokes equation are ‘eliminated’ in the higher-order formulation. Specifically, we are interested in the non-conservative viscous terms; while the body force $b_i$ could also be non-conservative, we will not concern ourselves with this possibility. Direct substitution of $mathcal {R}_i$ into (5.5) generates many terms, but it is found that only one is cancelled: the viscous Laplacian of the (time derivative of the) velocity, namely $mu dot {u}_{i,jj}$. This term mutually appears from the $rho dot {mathcal {R}}_i$ and $-mu mathcal {R}_{i,jj}$ terms in (5.5). To maintain notional clarity, we write the residual as

(5.6)begin{equation} mathcal{R}_i = rhodot{u}_i – mu u_{i,kk} + tilde{mathcal{R}}_i, end{equation}

where index $k$ has been used to avoid confusion with gradient operators in (5.5) having index $j$, and $tilde {mathcal {R}}_i = rho u_ku_{i,k}-rho b_i$ are the remaining terms in the residual. Substituting the above into the first and last terms of (5.5), cancelling the aforementioned $mu dot {u}_{i,jj}$ term and then dividing out by the density gives

(5.7)begin{equation} rhoddot{u}_i+dot{tilde{mathcal{R}}}_i+u_jmathcal{R}_{i,j} = mathcal{R}_ju_{j,i} – nu[-mu u_{i,kkjj} + tilde{mathcal{R}}_{i,jj} ], end{equation}

where $nu =mu /rho$ is the kinematic viscosity (recall that all variables are non-dimensional). We see that this equation is second order in time and fourth order in space. Viscous terms still appear in the equation including second- and third-order spatial derivatives. Nevertheless, the technique detailed by Sanders (Reference Sanders2023*b*) and employed here evidently ensures that (5.7) has a corresponding Hamiltonian structure.

## 6. Case study

We can explore how this method can be applied by considering a simplified example with a known field solution. In looking at the variety of cases in which the Navier–Stokes equations have a known analytical solution, the simplest are those involving steady flows. While the Euler–Lagrange equations (3.9), (3.10) can be written for these cases, the corresponding Hamilton–Jacobi equation is trivial because for steady flows the fields are already equal to their initial values.

It is therefore worthwhile to examine the simplest unsteady flows, which should result in a non-trivial Hamilton–Jacobi equation. Indeed, there exists a class of flows for which the Navier–Stokes equations take the same simplified form: those in which the flow is incompressible and unidirectional (Batchelor Reference Batchelor2000). This class of problems include both of Stokes’s flows (Stokes Reference Stokes1851), in which a semi-infinite fluid is influenced by a boundary moving in its own plane. In the first of these cases, the boundary is impulsively started and in the second, the boundary oscillates. We can also include developing flow in a channel or pipe. The only difference between these flows results from initial and boundary conditions, but the Navier–Stokes equations and therefore the present Hamilton–Jacobi equation take the same form.

Here we will examine the case in which there is motion only in the $x_1$ direction, and the velocities take the form ${u_{i}} = {u_{1}(x_{2},t),0,0}$. In the absence of a body force, our pressure gradient in the $x_1$ direction is solely a function of time and the pressure gradients in the $x_2$ and $x_3$ directions are zero. There are thus only two unknown field quantities: $u_{1}(x_{2},t)$ and $p(x_{1},t)$, where $p$ is linear in $x_{1}$. The field equation of primary interest is

(6.1)begin{equation} mathcal{R}_{1} equiv rho dot{u}_{1} + p_{,1} – mu u_{1,22} = 0, end{equation}

and the remaining field equations are satisfied automatically by the assumed form of the fields. Following the procedure described above, the momenta conjugate to $u_{1}$ and $p$ are given by

(6.2*a*,*b*)begin{equation} {rm pi}_{1}equiv rhomathcal{R}_{1}, quad {rm pi}_{4}equiv0.end{equation}

This results in a Hamiltonian density given by

(6.3)begin{equation} mathcal{H}^{*}=frac{1}{2}frac{1}{{rho}^{2}}{rm pi}_{1}{rm pi}_{1}-frac{1}{rho}(p_{,1}-mu u_{1,22}){rm pi}_{1}. end{equation}

Hamilton’s principal functional ${S}^{*}={S}^{*}[u_{1},p,t]$ can be expressed as an integral over $x_{2}$ only, since the other spatial coordinates do not appear and may be integrated out. In this way, we may write the Hamilton–Jacobi equation as follows:

(6.4)begin{equation} int text{d}kern 0.06em x_2left[frac{1}{2}frac{1}{{rho}^{2}}frac{delta {S}^{*}}{delta u_{1}}frac{delta {S}^{*}}{delta u_{1}}-frac{1}{rho}(p_{,1}-mu u_{1,22})frac{delta {S}^{*}}{delta u_{1}}right]+frac{partial {S}^{*}}{partial t}=0, end{equation}

with $delta {S}^{*}/delta p=0$. The solution to (6.4) would provide a canonical transformation to a new set of coordinates, giving analytical expressions for $(u_{1},p)$.

Despite knowing the analytical solution for these fields in this particular example, the present authors have not been able to solve this Hamilton–Jacobi equation, since again the solution of such equations is itself an open problem (Michal Reference Michal1926; Jordan & Pauli Reference Jordan and Pauli1928; Levy Reference Levy1951; Tatarskii Reference Tatarskii1961; Syavavko & Mel’nichak Reference Syavavko and Mel’nichak1974; Dieudonne Reference Dieudonne1981; Koval’chik Reference Koval’chik1993). This example therefore appears to be a good place to start for tackling the general problem. Another interesting example to consider might be a two-dimensional Taylor–Green vortex such as that considered by Wu, Ma & Zhou (Reference Wu, Ma and Zhou2006).

## 7. Conclusion

This paper has presented a novel Hamiltonian formulation of the isotropic Navier–Stokes problem for both compressible and incompressible fluids. This canonical formulation opens several previously unexplored avenues toward a final resolution of the problem, which we briefly describe below.

Perhaps the most obvious route would be to solve the Hamilton–Jacobi equation – either (4.11) for the compressible case or (4.12) for the incompressible case – for Hamilton’s principal functional ${S}^{*}[u_{i},p,t]$ directly. If a complete solution for ${S}^{*}$ can be found, it will lead via canonical transformation to a new set of fields which are equal to their initial values, thereby giving analytical expressions for the original velocity and pressure fields. Alternatively, if one can simply establish based on emerging analytical techniques that a complete solution to this Hamilton–Jacobi equation does (or does not) exist under the usual assumptions, that will also settle the question of existence of solutions.

An alternative strategy might be to investigate the corresponding Lagrangian formulation based on the action $mathcal {S}^{*}$ as given by (3.3). Because the first- and second-order formulations are mathematically equivalent (recall the proof in § 3.2), $mathcal {S}^{*}$ must have as many local minima as there are solutions to the traditional, first-order formulation. Intuitively, it seems as though it ought to be possible to determine – or at least to establish bounds on – the number of critical points an action has based on the form of the Lagrangian (Van den Berg & VanderVorst Reference Van den Berg and VanderVorst2002; Kalies & VanderVorst Reference Kalies and VanderVorst2004). If one can establish that, under the usual assumptions, $mathcal {S}^{*}$ always has exactly one local minimum, or else demonstrate that there are cases where it fails to achieve a local minimum, that too will resolve the question of existence and uniqueness.

By no means is either of the above programs trivial. As pointed out in § 4.1, solving equations containing Volterra (Reference Volterra1930) functional derivatives is itself a long-standing open problem in mathematics, which has received very little attention since the first half of the twentieth century (Michal Reference Michal1926; Jordan & Pauli Reference Jordan and Pauli1928; Levy Reference Levy1951; Tatarskii Reference Tatarskii1961; Syavavko & Mel’nichak Reference Syavavko and Mel’nichak1974; Dieudonne Reference Dieudonne1981; Koval’chik Reference Koval’chik1993). One might even go so far as to call it a ‘forgotten’ open problem (as did one of the reviewers of the present paper, who generously drew our attention to Jordan & Pauli Reference Jordan and Pauli1928; Levy Reference Levy1951; Tatarskii Reference Tatarskii1961; Syavavko & Mel’nichak Reference Syavavko and Mel’nichak1974; Dieudonne Reference Dieudonne1981; Koval’chik Reference Koval’chik1993). We see the lack of work on such equations as a challenge, yes, but at the same time we also see it as a significant opportunity for advancing the field of analytical continuum mechanics. Perhaps, despite an apparent increase in complexity, a rigorous theory of such equations can be developed after all, in which case the present formulation of the Navier–Stokes problem might be solved as one example. We submit that, at the very least, such an endeavour merits further study, which we intend to continue in future work.

Finally, it is worth noting that the techniques employed here are by no means specific to the Navier–Stokes problem, nor are they restricted to the field of classical mechanics. The suitably averaged principle of least squares (Sanders Reference Sanders2021, Reference Sanders2022, Reference Sanders2023*a*,Reference Sanders*b*; Sanders & Inman Reference Sanders and Inman2023) can be applied to any traditionally non-Hamiltonian dynamical system in order to formulate a mathematically equivalent higher-order Hamiltonian system. It is believed that this fundamental result will also find uses in other branches of pure and applied mathematics.