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 Finlaysonb). 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 Finlayson1972a,Reference Finlaysonb), 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
\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
\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 Sanders2023a,Reference Sandersb) 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
\begin{align}
\mathcal{S}^{*}[v]=\int\text{d}t\big(\tfrac{1}{2}\mathcal{R}^{2}\big)=\int
\text{d}t\big[\tfrac{1}{2}\big(\dot{v}^{2}+2v\dot{v}+v^{2}\big)\big]
\sim\int\text{d}t\big[\tfrac{1}{2}\big(\dot{v}^{2}+v^{2}\big)\big],
\end{align}
where we have used the fact that 
$2v\dot {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 Sanders2023a,Reference Sandersb; 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
\begin{equation} \delta\mathcal{S}^{*}=\int\text{d}t\left[\dot{v}\delta\dot{v}+v\delta v\right]=\int\text{d}t\left[(-\ddot{v}+v)\delta v\right]+\left[\dot{v}\delta v\right]_{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.14159\ldots$ does not appear in the present work, there will be no ambiguity.
The corresponding Hamiltonian is obtained via the Legendre transform
\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
\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
\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
\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
\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:
\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
\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
\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
\begin{equation} \phi=\frac{\partial{S}^{*}}{\partial\alpha}={\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
\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
\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 Sanders2023a,Reference Sandersb; 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
$$\begin{gather} \rho \dot{u}_{i} + \rho u_{i,j}u_{j} + p_{,i} - \mu u_{i,jj} - \rho b_{i} = 0, \end{gather}$$
$$\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
\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
$$\begin{gather} \dot{u}_{i}=\frac{\delta H^{*}}{\delta{\rm \pi}_{i}}, \quad\dot{p}=\frac{\delta H^{*}}{\delta{\rm \pi}_{4}}, \end{gather}$$
$$\begin{gather}\dot{\rm \pi}_{i}={-}\frac{\delta H^{*}}{\delta u_{i}}, \quad\dot{\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.19a,b) and (1.20a,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 Berdichevsky2009a,Reference Berdichevskyb; 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 Finlayson1972a,Reference Finlaysonb; 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 Finlayson1972a,Reference Finlaysonb). 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 Finlayson1972a; 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 Yoshimura2019a,Reference Gay-Balmaz and Yoshimurab; Rashad et al. Reference Rashad, Califano, Brugnoli, Schuller and Stramigioli2021; Gonzalez & Taha Reference Gonzalez and Taha2022; Taha & Gonzalez Reference Taha and Gonzalez2022; Sanders Reference Sanders2023b).
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 Yoshimura2019a,Reference Gay-Balmaz and Yoshimurab), 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 Sanders2023a,Reference Sandersb) 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 Finlayson1972a,Reference Finlaysonb; 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
$$\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}$$
$$\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 =-2\mu /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 Finlayson1972a,Reference Finlaysonb). 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 Finlayson1972a,Reference Finlaysonb), nevertheless a Hamiltonian for the system may be found by considering a second-order formulation. Following Sanders (Reference Sanders2023b), 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:
\begin{equation} \mathcal{S}^{*}[u_{j},p,\rho]=\int\text{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 (
$t\in [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 Finlayson1972b) 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:
$$\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}$$
$$\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
\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
\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
\begin{align} \delta\mathcal{S}^{*}&=\int\text{d}^{4}x\left\{\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}
\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}
\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}
\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}
\begin{align} &\quad +\int\text{d}^{3}x\left[\hat{\rho}\mathcal{R}_{i} \delta{u}_{i}+\hat{\rho}'\mathcal{R}_{4}\delta{p}\right]_{t_{1}}^ {t_{2}} \end{align}
\begin{align} &\quad +\int\text{d}^{2}x\,\text{d}t\left\{\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}
\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 p\right\}, \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
\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}
\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 Sanders2023b). 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.8e). The momenta conjugate to the velocities 
$u_{i}$ are
\begin{equation} {\rm \pi}_{i}\equiv \hat{\rho}\mathcal{R}_{i}, \end{equation}
and the momentum conjugate to the pressure 
$p$ is
\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.8e) vanishes identically. Interestingly, for the actual motion (
$\mathcal {R}_{i}=0,\mathcal {R}_{4}=0$), the volumetric integral (3.8e) 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 Sanders2023b).
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
$$\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}$$
$$\begin{gather}\delta u_{i,j}:\quad - \mu \mathcal{R}_{i}n_{j}- (\mu+\lambda)\mathcal{R}_{j}n_{i}=0, \end{gather}$$
$$\begin{gather}\delta p:\quad\mathcal{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
\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
\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
\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
$$\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}$$
$$\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
\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:
$$\begin{gather} \dot{u}_{i}=\frac{\delta H^{*}}{\delta{\rm \pi}_{i}}, \quad\dot{p}=\frac{\delta H^{*}}{\delta{\rm \pi}_{4}}, \end{gather}$$
$$\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.6a,b) reproduce (4.3) and (4.4), respectively. Equations (4.7a,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.6a,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.7a,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
\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,
\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 
$0\equiv \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ödinger1926a,Reference Schrödingerb). 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:
\begin{equation} {H}^{*}\left[u_{i},p,\frac{\delta {S}^{*}}{\delta u_{j}},\frac{\delta {S}^{*}}{\delta p};t\right]+\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
\begin{align} &\int \text{d}^{3}x\left[\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.6a,b) and (4.7a,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
\begin{align} \int \text{d}^{3}x\left[\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 Sanders2023b) 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:
\begin{equation} {-}\mathcal{R}_{i,i}=0 \quad \forall x_j\in\mathcal{V}, \quad\text{subject to}\quad \mathcal{R}_in_i = 0 \quad \forall x_j\in\partial\mathcal{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
\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
\begin{equation} f(x_j,t) ={-}n_i[\rho\dot{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
\begin{equation} \rho \mathcal{R}_{i}u_{j}n_{j}+\mu \mathcal{R}_{i,j}n_{j}=0 \quad\text{and}\quad - \mu \mathcal{R}_{i}n_{j}=0 \quad \forall x_j\in\partial\mathcal{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}_i\equiv 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.2a–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
\begin{equation} \rho\dot{\mathcal{R}}_i + \rho u_j\mathcal{R}_{i,j} = \rho \mathcal{R}_{j}u_{j,i}- \mu \mathcal{R}_{i,jj} \quad \forall x_j\in\mathcal{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
\begin{equation} \mathcal{R}_i = \rho\dot{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
\begin{equation} \rho\ddot{u}_i+\dot{\tilde{\mathcal{R}}}_i+u_j\mathcal{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 Sanders2023b) 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
\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
\begin{equation} {\rm \pi}_{1}\equiv \rho\mathcal{R}_{1}, \quad {\rm \pi}_{4}\equiv0.\end{equation}This results in a Hamiltonian density given by
\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:
\begin{equation} \int \text{d}\kern 0.06em x_2\left[\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 Sanders2023a,Reference Sandersb; 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.
Acknowledgments
The authors wish to thank Maggie Sanders for posing thought-provoking questions during the development of this paper. We would also like to thank the reviewers for taking the time to provide genuinely constructive feedback and suggestions, which have improved the quality of the paper immeasurably. One reviewer in particular provided not only a number of important technical clarifications, but also a list of 35 additional references to augment the literature review. We express our profound gratitude for the reviewers’ efforts.
Funding
The authors declare that no funds, grants or other support were received during the preparation of this manuscript.
Declaration of interests
The authors report no conflict of interest.
Data availability
Any data appearing in the present work are available from the corresponding author upon reasonable request.
Author contributions
J.W.S. conceived the original idea for the work, wrote the Abstract, Introduction, Lagrangian formulation, Hamiltonian formulation, Conclusion and Appendix, and contributed to the Literature Review. A.C.D. wrote the Discussion, contributed to the Analysis and Literature Review and provided expertise in the area of fluid mechanics. N.J.W. wrote the Case Study, contributed to the Literature Review and provided expertise in the area of fluid mechanics. G.A.E. contributed to the Literature Review and provided expertise in the area of fluid mechanics. K.L.S. organized and wrote the Literature Review and provided expertise in the area of fluid mechanics. J.C.B. contributed to the Introduction and Literature Review, and provided expertise in the areas of analytical mechanics and variational calculus. All authors contributed equally in checking the accuracy of the work, discussing the results together as a group on multiple occasions, drawing conclusions and proofreading the manuscript.
Appendix A. Derivation of the Hamilton–Jacobi equation
In what follows, it is important to distinguish between two sets of solutions: solutions to the second-order Euler–Lagrange equations (3.9) and (3.10), and solutions to the first-order field equations (3.4) and (3.5). The latter are a subset of the former, but not vice versa.
We define Hamilton's principal functional 
${S}^{*}[u_{i},p,t]$ as
\begin{equation} {S}^{*}[u_{i},p,t]\equiv\int_{t_{0}}^{t}{\rm d}t(\tilde{L}^{*}), \end{equation}
where 
$t$ is the current (variable) time, 
$t_{0}$ is an arbitrary initial time and 
$\tilde {L}^{*}$ denotes the Lagrangian (3.6) evaluated for fields satisfying the second-order Euler–Lagrange equations (3.9) and (3.10) – not necessarily the first-order field equations (3.4) and (3.5). Crucially, because the fields do not necessarily satisfy the first-order field equations, it is not appropriate to take 
${S}^{*}\equiv 0$. We imagine that the time integral has already been carried out, so that the functional 
${S}^{*}$ may be regarded as an integral over 
$\mathcal {V}$, that is
\begin{equation} {S}^{*}[u_{i},p,t]=\int\text{d}^{3}x(s^{*}), \end{equation}
where 
$s^{*}$ is the Lagrangian density (3.7) evaluated for fields satisfying the second-order Euler–Lagrange equations (which will be denoted 
$\tilde {\mathcal {L}}^{*}$) integrated over time from 
$t_{0}$ to 
$t$.
Starting from (A1) and evaluating the first variation of 
${S}^{*}$ as we did with the action 
$\mathcal {S}^{*}$ in § 3, we find that
\begin{equation} \delta{S}^{*}=\int\text{d}^{3}x\big[{\rm \pi}_{i}\delta{u}_{i}+{\rm \pi}_{4}\delta{p}\big]_{t_{0}}^{t}, \end{equation}where we have used the fact that the second-order Euler–Lagrange equations (3.9) and (3.10) are satisfied by definition and we have also enforced the natural boundary conditions. But, by the definition of Volterra (Reference Volterra1930) functional derivatives, we have that
\begin{equation} \delta{S}^{*}=\int\text{d}^{3}x\left[\frac{\delta{S}^{*}}{\delta u_{i}}\delta{u}_{i}+\frac{\delta{S}^{*}}{\delta p}\delta{p}\right]. \end{equation}
In this way, we may identify the conjugate momenta with the functional derivatives of 
${S}^{*}$:
\begin{equation} {\rm \pi}_{i}=\frac{\delta{S}^{*}}{\delta u_{i}}, \quad {\rm \pi}_{4}=\frac{\delta{S}^{*}}{\delta p}. \end{equation}Now, starting from (A1) and evaluating the time derivative, we find that
\begin{equation} \frac{\text{d}{S}^{*}}{\text{d}t}=\tilde{L}^{*}=\int\text{d}^{3}x\big(\tilde{\mathcal{L}}^{*}\big). \end{equation}But, by the chain rule,
\begin{equation} \frac{\text{d}{S}^{*}}{\text{d}t}=\frac{\partial{S}^{*}}{\partial t}+\int\text{d}^{3}x\left[\frac{\delta{S}^{*}}{\delta u_{i}}\dot{u}_{i}+\frac{\delta{S}^{*}}{\delta p}\dot{p}\right]. \end{equation}In this way, we find that
\begin{equation} \int\text{d}^{3}x\left[\frac{\delta{S}^{*}}{\delta u_{i}}\dot{u}_{i}+\frac{\delta{S}^{*}}{\delta p}\dot{p}-\tilde{\mathcal{L}}^{*}\right]+\frac{\partial{S}^{*}}{\partial t}=0. \end{equation}
Here, the integral is simply the Hamiltonian 
$H^{*}$, with the conjugate momenta replaced by the functional derivatives in accordance with (A5a,b). Hence, we arrive at the Hamilton–Jacobi equation
\begin{equation} {H}^{*}\left[u_{i},p,\frac{\delta {S}^{*}}{\delta u_{j}},\frac{\delta {S}^{*}}{\delta p};t\right]+\frac{\partial {S}^{*}}{\partial t}=0, \end{equation}as claimed in (4.10).
Appendix B. Inviscid flow
In this appendix, we consider what becomes of the present formulation in the inviscid limit, i.e. upon setting the viscosities 
$\mu =\lambda =0$. While this is a relatively simple exercise, it is of interest by virtue of its relationship to the classical Euler equations, which of course are already conservative and so do not require a higher-order formulation to put them in Hamiltonian form (Olver Reference Olver1982).
The vanishing of viscous forces affects the residuals of the momentum balance equations
\begin{equation} \bar{\mathcal{R}}_{i}[u_{j},p,\rho;x_{k},t]\equiv \rho \dot{u}_{i} + \rho u_{i,j}u_{j} + p_{,i} - \rho b_{i} = 0, \end{equation}
where we are using a bar to distinguish the inviscid residuals (
$\bar {\mathcal {R}}_{i}$) from their more general forms given in (3.1). Applying the equation of state 
$\rho =\hat {\rho }(p)$, we obtain
\begin{equation} \bar{\mathcal{R}}_{i}[u_{j},p;x_{k},t]\equiv \hat{\rho} \dot{u}_{i} + \hat{\rho} u_{i,j}u_{j} + p_{,i} - \hat{\rho} b_{i} = 0. \end{equation}
The residual 
$\mathcal {R}_{4}$ of the mass balance is unaffected by viscous effects and will remain as it was in the main text.
The Lagrangian density is now
\begin{equation} \bar{\mathcal{L}}^{*}[u_{j},p;x_{k},t]=\tfrac{1}{2}\bar{\mathcal{R}}_{i}\bar{\mathcal{R}}_{i}+\tfrac{1}{2}\mathcal{R}_{4}\mathcal{R}_{4}, \end{equation}and the associated Euler–Lagrange equations are as follows:
\begin{equation} \delta{u}_{i}: \quad -\frac{\partial}{\partial t} (\hat{\rho}\bar{\mathcal{R}}_{i}) -\frac{\partial}{\partial x_{j}} (\hat{\rho} \bar{\mathcal{R}}_{i}u_{j})+ \hat{\rho} \bar{\mathcal{R}}_{j}u_{j,i} + \hat{\rho}'\mathcal{R}_{4}p_{,i}-\frac{\partial}{\partial x_{i}}( \hat{\rho} \mathcal{R}_{4})=0, \end{equation}
\begin{align} \delta p:\quad & \hat{\rho}'\bar{\mathcal{R}}_{i} \dot{u}_{i}+ \hat{\rho}' \bar{\mathcal{R}}_{i}u_{i,j}u_{j}- \bar{\mathcal{R}}_{i,i}- \hat{\rho}' \bar{\mathcal{R}}_{i}b_{i} + \hat{\rho}''\mathcal{R}_{4}\dot{p} \nonumber\\ &-\frac{\partial}{\partial t}( \hat{\rho}'\mathcal{R}_{4}) + \hat{\rho}''\mathcal{R}_{4}p_{,i}u_{i} - \frac{\partial}{\partial x_{i}}(\hat{\rho}'\mathcal{R}_{4}u_{i})+ \hat{\rho}' \mathcal{R}_{4}u_{i,i} =0. \end{align}Again, all four Euler–Lagrange equations are second order in time, as they involve time derivatives of the residuals. This is a mathematically equivalent second-order formulation of the inviscid Euler equations, as claimed in the main text.
The momenta conjugate to the velocities 
$u_{i}$ are now
\begin{equation} \bar{\rm \pi}_{i}\equiv \hat{\rho}\bar{\mathcal{R}}_{i}, \end{equation}
and the momentum 
${\rm \pi} _{4}$ conjugate to the pressure 
$p$ remains as it was in the main text. For compressible flows, the resulting Hamiltonian density takes the form
\begin{align} \bar{\mathcal{H}}^{*}[u_{i},p,\bar{\rm \pi}_{j},{\rm \pi}_{4};x_{k},t]&= \frac{1}{2}\frac{1}{(\hat{\rho})^{2}}{\bar{\rm \pi}_{i}\bar{\rm \pi}_{i}}- \frac{1}{\hat{\rho}}(\hat{\rho}u_{i,j}u_{j} + p_{,i} - \hat{\rho}b_{i})\bar{\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}Once again, the higher-order Hamiltonian has nothing to do with the total mechanical energy of the system, but rather vanishes for the actual motion (just as it does in the case of viscous flows). With this Hamiltonian density, Hamilton's equations assume the usual canonical form
$$\begin{gather} \dot{u}_{i}=\frac{\delta \bar{H}^{*}}{\delta\bar{\rm \pi}_{i}}, \quad \dot{p}= \frac{\delta \bar{H}^{*}}{\delta{\rm \pi}_{4}}, \end{gather}$$
$$\begin{gather}\dot{\bar{\rm \pi}}_{i}={-}\frac{\delta \bar{H}^{*}}{\delta u_{i}}, \quad \dot{\rm \pi}_{4}={-}\frac{\delta \bar{H}^{*}}{\delta p}, \end{gather}$$
where 
$\bar {H}^{*}$ denotes 
$\bar {\mathcal {H}}^{*}$ integrated over the volume 
$\mathcal {V}$. Equations (B9a,b) reproduce the Euler–Lagrange equations (B4) and (B5). In the case of incompressible flows, the Hamiltonian density reduces to
\begin{equation} \bar{\mathcal{H}}^{*}[u_{i},p,\bar{\rm \pi}_{j};x_{k},t]=\frac{1}{2}\frac{1}{{\rho}^{2}}{ \bar{\rm \pi}_{i}\bar{\rm \pi}_{i}}-\frac{1}{\rho}({\rho}u_{i,j}u_{j} + p_{,i} - {\rho}b_{i})\bar{\rm \pi}_{i}, \quad \hat{\rho}'\equiv0.\end{equation}
Hamilton's equations 
$\dot {u}_{i}=\delta \bar {H}^{*}/\delta \bar {{\rm \pi} }_{i}$, 
$\dot {\bar {{\rm \pi} }}_{i}=-\delta \bar {H}^{*}/\delta u_{i}$ and 
$0\equiv \dot {{\rm \pi} }_{4}=-\delta \bar {H}^{*}/\delta p$ still apply, but 
$\dot {p}=\delta \bar {H}^{*}/\delta {\rm \pi}_{4}$ must be replaced by the constraint that 
$u_{i,i}=0$, once again consistent with the fact that the pressure serves as Lagrange multiplier for the incompressibility constraint (Lanczos Reference Lanczos1970; Badin & Crisciani Reference Badin and Crisciani2018).
The Hamilton–Jacobi equation
\begin{equation} \bar{H}^{*}\left[u_{i},p,\frac{\delta \bar{S}^{*}}{\delta u_{j}},\frac{\delta \bar{S}^{*}}{\delta p};t\right]+\frac{\partial \bar{S}^{*}}{\partial t}=0 \end{equation}assumes the following forms for compressible and incompressible flows, respectively:
\begin{align} &\int \text{d}^{3}x\left[\frac{1}{2}\frac{1}{(\hat{\rho})^{2}}\frac{\delta \bar{S}^{*}}{\delta u_{i}}\frac{\delta \bar{S}^{*}}{\delta u_{i}}-\frac{1}{\hat{\rho}}(\hat{\rho}u_{i,j}u_{j} + p_{,i} - \hat{\rho}b_{i})\frac{\delta \bar{S}^{*}}{\delta u_{i}}\right. \nonumber\\ &\quad \left.+\,\frac{1}{2}\frac{1}{(\hat{\rho}')^{2}}\frac{\delta \bar{S}^{*}}{\delta p}\frac{\delta \bar{S}^{*}}{\delta p}-\frac{1}{\hat{\rho}'}({\hat{\rho}'p_{,i}u_{i} + \hat{\rho}u_{i,i}})\frac{\delta \bar{S}^{*}}{\delta p}\right]+\frac{\partial \bar{S}^{*}}{\partial t}=0, \quad \hat{\rho}'\neq0, \end{align}and
\begin{equation} \int \text{d}^{3}x\left[\frac{1}{2}\frac{1}{{\rho}^{2}}\frac{\delta \bar{S}^{*}}{\delta u_{i}}\frac{\delta \bar{S}^{*}}{\delta u_{i}}-\frac{1}{\rho}({\rho}u_{i,j}u_{j} + p_{,i} - {\rho}b_{i})\frac{\delta \bar{S}^{*}}{\delta u_{i}}\right]+\frac{\partial \bar{S}^{*}}{\partial t}=0, \quad \hat{\rho}'\equiv0, \end{equation}
where in the latter case we have 
$\delta \bar {S}^{*}/\delta p=0$, since 
${\rm \pi} _{4}$ vanishes identically for incompressible flow.
The second-order formulation described above is fundamentally different from – although again still mathematically equivalent to – the classical Hamiltonian formulation of the first-order Euler equations (Olver Reference Olver1982).
Open access
