A canonical Hamiltonian formulation of the Navier-Stokes problem

This paper presents a novel Hamiltonian formulation of the isotropic Navier-Stokes problem based on a minimum-action principle derived from the principle of least squares. This formulation uses the velocities u i ( x j , t ) and pressure p ( x j , t ) as the field quantities to be varied, along with canoni-cally conjugate momenta deduced from the analysis. From these, a conserved Hamiltonian functional H ∗ satisfying Hamilton’s canonical equations is constructed, and the associated Hamilton-Jacobi equation is formulated for both compressible and incompressible flows. This Hamilton-Jacobi equation reduces the problem of finding four separate field quantities ( u i , p ) to that of finding a single scalar functional in those fields—Hamilton’s principal functional S ∗ [ t, u i , p ] . Moreover, the transformation theory of Hamilton and Jacobi now provides a prescribed recipe for solving the Navier-Stokes problem: Find S ∗ . If an analytical expression for S ∗ can be obtained, it will lead via canonical transformation to a new set of fields which are simply equal to their initial values, giving analytical expressions for the original velocity and pressure fields. Failing that, if one can only show that a complete solution to this Hamilton-Jacobi equation does or does not exist, that will also resolve the question of existence of solutions. The method employed here is not specific to the Navier-Stokes problem or even to classical mechanics, and can be applied to any traditionally non-Hamiltonian problem.


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 [1][2][3].There is no action integral for which Hamilton's principle [4][5][6] yields the Navier-Stokes equations [7][8][9][10][11][12][13] in their usual form [1][2][3], 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.

A motivating example
Consider the first-order initial-value problem with unique solution v(t) = 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 [7][8][9][10][11][12][13], this too is an intrinsically non-Hamiltonian problem, in that there is no action S for which Hamilton's principle (δS = 0) yields the governing equation v = −v.Yet if we simply differentiate both sides of the equation (v = − v), use the original equation to write v = −v, and apply the additional initial condition v(0) = −v(0) = −1, we arrive at the mathematically equivalent second-order problem which has the same unique solution v(t) = 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 [14][15][16][17][18], the associated action can be obtained by writing the original equation in standard form (R ≡ v + v = 0), squaring the residual R, and integrating over time: where we have used the fact that 2v v = d(v 2 )/dt is a total time derivative and can therefore be excluded from the action without changing the resulting Euler-Lagrange equation [19].This is the so-called "time-averaged principle of least squares" [14][15][16][17][18]: since R = 0 is a local minimum of R 2 , it is also a local minimum of ´dt(R 2 ).Varying v, the first variation of S * is δS * = ˆdt vδ v + vδv = ˆdt (−v + v)δv + vδv yielding the second-order equation v = v and revealing the canonically conjugate "momentum" π ≡ v.Here and in what follows, we will use the symbol π 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... does not appear in the present work, there will be no ambiguity.
The corresponding Hamiltonian is obtained via the Legendre transform 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 π ≡ v = −v.
We note in passing that Liouville's theorem is satisfied, as the motion occurs along the line π = −v, so that the phase-space volume, being always zero, is conserved.Hamilton's equations are mathematically equivalent to the second-order problem v = v and therefore also mathematically equivalent to the original, first-order problem.The associated Hamilton-Jacobi equation [4][5][6][19][20][21][22] is where Hamilton's principal function S * = S * (v, t) serves as the generating function for a canonical transformation to a new coordinate φ 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) does not work, as the reader may check.
Instead, let us use a trial solution of the form 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 In order for this equation to hold for all v, we must have the following: where α is a constant of integration which will be used to transform to the new coordinate, and where γ is another constant of integration which is simply additive and can therefore be discarded.In this way, we have that With one constant of integration (α) to match the single degree of freedom (v), this is a complete solution to the Hamilton-Jacobi equation.The new coordinate φ (which is constant and equal to its initial value) is obtained via the canonical transformation The numerical value of α is in turn obtained via the canonical relation which, evaluated at t = 0, gives α = −2 (recall that π = v, and v(0) = −v(0) = −1).Using the fact that the new coordinate φ is equal to its initial value, we have that giving the correct solution v(t) = 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 [14][15][16][17][18].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 [7][8][9][10][11][12][13].

The Navier-Stokes problem
The incompressible Navier-Stokes equations [7][8][9][10][11][12][13] are given by where ρ is the constant and uniform density, 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, µ is the dynamic viscosity, a dot over a symbol denotes a partial time derivative ( ui = ∂u i /∂t), a comma in a subscript indicates a spatial gradient (p ,i = ∂p/∂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 (x, t) and u(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 where (π i , π 4 ) are suitable "momenta" conjugate to the field quantities (u i , p), such that Hamilton's canonical equations constitute a mathematically equivalent second-order formulation of the problem, where δH * /δu i , δH * /δp, δH * /δπ i , and δH * /δπ 4 are the Volterra [23] 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 ṗ = δH * /δπ 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 [19,24] (refer to p. 361 of [19] and pp.137 and 141 of [24]).
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 (19) and (20) for the mathematically equivalent second-order problem, along with the accompanying Hamilton-Jacobi equation [4][5][6][19][20][21][22]. 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, Section 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.

Literature review
The field of analytical mechanics, with foundations planted in Hamilton's principle of stationary action [4][5][6] or d'Alembert's principle of virtual work [25], 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 [26][27][28][29][30][31].The supporting mathematics of the calculus of variations as well as symplectic and differential geometry can also be found in many excellent sources [32][33][34][35][36][37][38][39].It is therefore no surprise that researchers have been applying analytical formalism to classical fluids dating back to the time of Lagrange [24,[40][41][42][43][44][45][46][47].
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 [7][8][9][10][11][12][13]. 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 [8][9][10][11][12][13].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 [48].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, nonsingular solutions to the three-dimensional Navier-Stokes equations [49], and this constitutes one of the most famous unsolved problems in mathematics.
The application of analytical mechanics [32,37,50,51] to the field of fluid mechanics [19] has recently seen a resurgence in interest [47,[52][53][54][55][56][57][58] after a long history.In the absence of nonconservative forces, an inviscid fluid is a Hamiltonian system, and so the classical Hamiltonian theory applies.Serrin [59], Benjamin [60], and Holm et al. [61] have all described variational and Hamiltonian formulations of incompressible, inviscid fluid flow.Roberts [62] 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 [63] showed that the Euler equations of inviscid and incompressible fluid flow can be put into Hamiltonian form.Benjamin and Olver [64] 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 and Pego [65] 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 [66] 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 [1][2][3]49].Indeed, it is a well known theorem (first proven by Millikan [1]) that the Navier-Stokes equations in their usual form cannot be derived from a classical action principle [1][2][3].Millikan [1] 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 [1].
Oseledets [82] 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 and Fujitani [85], Jones [86], and Gay-Balmaz and Yoshimura [88,90,91], have enforced dissipation using a non-holonomic constraint on the entropy.Hochgerner [89] attempted to obtain a Hamiltonian interacting particle system that could accurately model 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. [92] modeled 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 and Taha [93], Taha and Gonzalez [94], and Taha, Gonzalez, and Shorbagy [95] have recently applied Gauss's principle of least constraint [96] to the Navier-Stokes problem.Using Gauss's principle, Taha, Gonzalez, and Shorbagy [95] 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 [96] is fundamentally different from Hamilton's principle [4][5][6].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 [14][15][16][17][18] has shown that 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 [15] and subsequently a new and more efficient algorithm for computing a damped system's resonant frequencies [17].Higher-derivative theories had been studied before in the realm of quantum gravity physics [97][98][99][100][101][102][103][104] 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 [1-3, 49, 84, 87], it will be shown here that 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.

Lagrangian formulation of the problem
Although we are primarily interested in the incompressible form of the equations given by ( 16) and ( 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 where ρ = ρ(x j , t) is the density field (now one of the unknown field quantities along with u i and p), and λ is an additional viscosity coefficient which, under Stokes's [7] hypothesis, is related to µ as λ = −2µ/3, ensuring that the mechanical pressure agrees with the thermodynamic pressure.Henceforth we will assume that all quantities have been suitably non-dimensionalized.The nondimensional (constant) viscosities in ( 21) and ( 22) may be regarded as inverse Reynolds numbers, and the non-dimensional pressure may be considered to be normalized by the inertial scale ρ 0 U 2 , with ρ 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, ( 21) and ( 22) 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 ρ = ρ(p, T ).The conservation equations along with the equation of state constitute six equations for the six unknowns fields (u i , p, T, ρ).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., dρ/dt = ρ + ρ ,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 ρ = ρ(p, T ) to specification of ρ = ρ 0 as a system parameter.Consequently, ( 22) reduces to ρ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 ( 21), ( 22) are first-order in time with respect to the field quantities u i and ρ.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 S for which Hamilton's principle (δS = 0) yields the first-order field equations [1][2][3].Finally, we note that in the incompressible limit, R 4 becomes independent of ρ and is no longer first-order in time.

Second-order formulation
Although the first-order formulation of the problem is intrinsically non-Hamiltonian [1][2][3], nevertheless a Hamiltonian for the system may be found by considering a second-order formulation.Following Sanders [18], we observe that the actual motion of the fluid corresponds to the particular fields (u i , p, ρ) for which the following action achieves a local minimum: where d 4 x = dx 1 dx 2 dx 3 dt, and the integral is carried out over both the control volume V occupied by the fluid (x j ∈ V) and the time period of interest (t ∈ [t 1 , t 2 ]).It must be emphasized that this action contains no new physics.Again, this is simply the principle of least squares [3] averaged over the spacetime occupied by the fluid.The entire physics of the problem are already contained in the residuals (R i , R 4 ).Without an equation of state relating ρ 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 ρ = ρ(p), with ρ a known function determined either from first principles or empirically.In this way, the density field may be eliminated in favor of the pressure field, and the field equations assume the following form: where ρ′ = dρ/dp.We note that, under equilibrium conditions, ρ′ is related to the speed of sound c and the bulk modulus K as ρ′ = 1/c 2 = ρ/K (for incompressible fluids, ρ′ ≡ 0 and the speed of sound and bulk modulus are both infinite).Having specified ρ(p), µ, λ, 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 ρ′ ≡ 0.
We pause here to note that, even though the residuals (R i , 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 ( 24) and (25); they do not vanish for every conceivable (u i , p).Thus it is not appropriate to take R i ≡ 0, R 4 ≡ 0. We will return to this point later when we discuss the Hamiltonian formulation of the problem.
For now, we note that the action where the integral is carried out over the volume V only (d 3 x = dx 1 dx 2 dx 3 ), with Lagrangian density Because the residuals (R i , 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 * ≡ 0.
As noted above, the actual motion of the fluid corresponds to the particular fields (u i , p) for which S * achieves a local minimum.To obtain the Euler-Lagrange equations, the conjugate momenta, and the natural auxiliary conditions, we insist that S * not vary to first order (δS * = 0) under small variations in the fields (δu i , δp).Evaluating δS * , integrating by parts, and collecting like terms, we find that where the purely volumetric integral (d 3 x) is carried out over V, the surface integral (d 2 x) is carried out over the boundary ∂V, and n i is the unit outward normal vector to ∂V.Note that, because we are using Eulerian coordinates x j , the volume element d 3 x is not to be varied.The Euler-Lagrange equations (which hold for all x j ∈ V) may be read directly from the spacetime (d 4 x) integral: δp It should be noted that all four Euler-Lagrange equations ( 29), ( 30) 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 [18].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 R i = 0 and 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 (28e).The momenta conjugate to the velocities u i are and the momentum conjugate to the pressure p is In the forthcoming Hamiltonian formulation, the conjugate momenta will be used to eliminate the (partial) time derivatives ( ui , ṗ) of the field quantities from the Hamiltonian.In general, Hamilton's principle would insist that the variations (δu i , δp) vanish at the endpoints t = t 1 and t = t 2 to ensure that the purely volumetric integral (28e) vanishes identically.Interestingly, for the actual motion (R i = 0, R 4 = 0), the volumetric integral (28e) already vanishes even without taking (δu i , δ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 [18].Although the conjugate momenta (π i , π 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 = ρu i , which we will see presently from the natural boundary conditions.These are read directly from the surface (d 2 x) integral: This last condition, (35), establishes a connection between the new conjugate momenta and the conventional linear momenta.Multiplying (35) by ρ, and noting that ρu i = ρu i = P i , we find that Evidently, boundary condition (35) states that the flux of the vector Π i ≡ π i + π 4 P i through the boundary ∂V should vanish.It is interesting that this new vector Π i contains both old and new momenta, with π 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.

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 Section 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 R i (x k , 0) = 0 and R 4 (x k , 0) = 0 for all x k ∈ V ∪ ∂V, along with boundary conditions making R i (x k , t) = 0, R i,j (x k , t) = 0, and R 4 (x k , t) = 0 for all x k ∈ ∂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 ( 24), ( 25) of the first-order formulation (R i = 0 and R 4 = 0) everywhere in V and at all times if and only if (u i , p) satisfies the Euler-Lagrange equations ( 29), (30) of the second-order formulation everywhere in V and at all times.Suppose first that (u i , p) satisfies the governing field equations ( 24), ( 25) of the first-order formulation everywhere in V and at all times.Then R i = 0, R 4 = 0, and (u i , p) is a trivial solution to the Euler-Lagrange equations ( 29), (30) of the second-order formulation.
Conversely, suppose that (u i , p) satisfies the Euler-Lagrange equations ( 29), (30) of the secondorder formulation everywhere in V and at all times.We note that (R i = 0, R 4 = 0) constitutes an equilibrium solution of the Euler-Lagrange equations ( 29), (30).Thus, because the initial conditions are chosen such that R i (x k , 0) = 0 and R 4 (x k , 0) = 0 for all x k ∈ V ∪ ∂V, and because the boundary conditions are chosen such that R i (x k , t) = 0, R i,j (x k , t) = 0, and R 4 (x k , t) = 0 for all x k ∈ ∂V and all times t, then R i and R 4 will remain identically zero everywhere in V for all future times.Thus, (u i , p) satisfies the governing field equations ( 24), ( 25) of the first-order formulation everywhere in V and at all times.This completes the proof, and we have established that the two formulations are equivalent.

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 with the Hamiltonian density H * obtained from the Lagrangian density L * via the Legendre transform 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 Section 1.1.
In order to write down Hamilton's equations, we must express H * in terms of the fields and the conjugate momenta, eliminating ui and ṗ.
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 ( 24), (25) of the first-order formulation, it does not vanish for every conceivable (u i , p).The latter would imply, according to Equations (42), that ui ≡ 0 and ṗ ≡ 0, which is not generally the case.This observation, and the fact that Equations ( 43) faithfully reproduce the Euler-Lagrange equations ( 29), (30) 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 (φ i , φ 4 ) and conjugate momenta for which the Hamiltonian does vanish identically.
To obtain the Hamiltonian for incompressible flow, we set ρ′ ≡ 0 from the beginning, in which case R 4 reduces to ρu i,i and π 4 vanishes identically, consistent with the fact that R 4 becomes independent of ṗ.The Hamiltonian density in turn reduces to or, in terms of the conjugate momenta, where the density ρ = ρ is a constant and we have used the fact that u i,i = 0. Hamilton's equations ui = δH * /δπ i , πi = −δH * /δu i , and 0 ≡ π4 = −δH * /δp still apply (and reproduce the corresponding equations in the incompressible limit), but ṗ = δH * /δπ 4 must be replaced by the constraint that u i,i = 0.That the incompressibility condition should take the place of the pressure equation ṗ = δH * /δπ 4 is consistent with the well known result that the pressure usually serves as the Lagrange multiplier for the incompressibility constraint [19,24].

Hamilton-Jacobi equation
One of the most significant aspects of the Hamiltonian formalism is that it leads to the transformation theory of Hamilton and Jacobi [4][5][6][19][20][21][22], celebrated both for unifying particle mechanics with wave optics [4] and for its relationship to the Schrödinger equation of quantum mechanics [105,106].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 S * , although they are related; see Appendix A), which is the solution to the following Hamilton-Jacobi equation: where δS * /δu i and δS * /δp are the Volterra [23] functional derivatives of S * with respect to the field quantities.Interested readers will find the derivation of (46) in Appendix A. Henceforth we will refer to S * as "Hamilton's principal functional."Substituting for the conjugate momenta in (41), we obtain for the Hamilton-Jacobi equation In contrast to the four original field equations-( 21) and ( 22)-the Hamilton-Jacobi equation ( 47) is a single equation in Hamilton's principal functional S * .This constitutes an equivalent formulation of the problem, as a complete and nontrivial solution to ( 47) is tantamount to an integration of Hamilton's equations ( 42) and ( 43) (note that it is not appropriate to take S * ≡ 0 for the same reason that it is not appropriate to take H * ≡ 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 (φ i , φ 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 ρ′ ≡ 0 and ρ′ appears in the denominators of terms in (47).Even so, the Hamiltonian formulation remains well posed in the incompressible limit.Recall that, with ρ′ ≡ 0, R 4 reduces to ρu i,i , π 4 vanishes identically, and the Hamiltonian density reduces to (45).Substituting for the conjugate momenta π i in (45), the corresponding Hamilton-Jacobi equation is with δS * /δp = 0, since again π 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 δS * /δp should vanish in the incompressible limit without knowing that in general π 4 = ρ′ 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 [19,24]).
It must be acknowledged that the Hamilton-Jacobi equation developed above (either (47) for the compressible case or (48) in the incompressible limit) contains Volterra [23] 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 [107][108][109][110][111][112][113].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 endeavor is worthwhile and merits further study.
We conclude this section by remarking that, in the inviscid limit (µ = λ = 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.

Discussion
In this section we provide some qualitative interpretations of the developments of Section 3.More specifically, we investigate the incompressible form (via constant, uniform density) of the Euler-Lagrange equations ( 29) and (30) when the residuals R i and R 4 are substituted.
Our motivation is again the simple example of Section 1.1 for which the first-order non-Hamiltonian system v = −v was converted to the second-order Hamiltonian system v = v by (manual) elimination of the non-conservative 'damping' term v (see [15] for a similar result for the damped harmonic oscillator converting from second-to fourth-order dynamics).Sanders [18] showed that the elimination process is 'automated' by the definition of the action in the first integral of (3), which is generalized to the action in (23) for our current continuum dynamics problem containing fields.
First consider the pressure equation ( 30) and corresponding natural boundary condition (35), which take the following incompressible forms: This higher-order field equation is simply the divergence of the residual R i .Upon substituting for R i from ( 21) and subsequently imposing the incompressible continuity condition u i,i = R 4 /ρ = 0 from ( 22), we obtain: which is a Poisson equation for the pressure.The boundary condition is a Neumann type requiring the specification of the normal pressure gradient, n i p ,i = p ,n ≡ f (x j , t), on the boundary: Equation ( 50) and boundary condition (51) 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.[12,114]).Next, we consider the velocity equations ( 29) which, at present, have a more elusive physical interpretation.Here, we instead begin with the natural boundary conditions (33) and (34), which are due to the δu i and δu i,j variations.The incompressible versions of these equations are: ρR i u j n j + µR i,j n j = 0 and − µR i n j = 0 ∀x j ∈ ∂V.
The boundary conditions involving the residual 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 R i ≡ 0 there.Note that the pressure of the actual motion on the boundary will be known from the simultaneous solution of (50).However, the gradient terms 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 higherorder governing equation, which will be seen shortly to be second-order in time and fourth-order in space.Again, recall the example of Section 1, whereby the system (2) 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 [10,11] 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.[115,116]).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 (29), which upon imposing incompressibility and expanding derivatives of product terms yields: The left-hand side is the material derivative of the residual 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 R i into (53) generates many terms, but it is found that only one is canceled: the viscous Laplacian of the (time derivative of the) velocity, namely µ ui,jj .This term mutually appears from the ρ Ṙi and −µR i,jj terms in (53).To maintain notional clarity, we write the residual as: where index k has been used to avoid confusion with gradient operators in (53) having index j, and Ri = ρu k u i,k − ρb i are the remaining terms in the residual.Substituting the above into the first and last terms of ( 53), canceling the aforementioned µ ui,jj term, and then dividing out by the density gives: where ν = µ/ρ 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 [18] and employed here evidently ensures that (55) has a corresponding Hamiltonian structure.

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 ( 29), ( 30) 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 nontrivial 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 [10].This class of problems include both of Stokes's flows [117], in which a semiinfinite 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 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 This results in a Hamiltonian density given by: 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: with δS * /δp = 0.The solution to (59) 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 [107][108][109][110][111][112][113].This example therefore appears to be a good place to start for tackling the general problem.Another interesting example to consider might be a 2D Taylor-Green vortex such as that considered by Wu et al. [118].

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 (47) for the compressible case or (48) 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 S * as given by (23).Because the first-and second-order formulations are mathematically equivalent (recall the proof in Section 3.2), 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 [98,99].If one can establish that, under the usual assumptions, 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 Section 4.1, solving equations containing Volterra [23] 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 [107][108][109][110][111][112][113].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 References [108][109][110][111][112][113]).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 endeavor 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 suitablyaveraged principle of least squares [14][15][16][17][18] 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.
where t is the current (variable) time, t 0 is an arbitrary initial time, and L * denotes the Lagrangian (26) evaluated for fields satisfying the second-order Euler-Lagrange equations ( 29) and ( 30)not necessarily the first-order field equations ( 24) and (25).Crucially, because the fields do not necessarily satisfy the first-order field equations, it is not appropriate to take S * ≡ 0. We imagine that the time integral has already been carried out, so that the functional S * may be regarded as an integral over V, that is, S * [u i , p, t] = ˆd3 x (s * ) , (A.2) where s * is the Lagrangian density ( 27) evaluated for fields satisfying the second-order Euler-Lagrange equations (which will be denoted L * ) integrated over time from t 0 to t.
Starting from (A.1) and evaluating the first variation of S * as we did with the action S * in Section 3, we find that δS * = ˆd3 x π i δu i + π 4 δp where we have used the fact that the second-order Euler-Lagrange equations ( 29) and ( 30) are satisfied by definition and we have also enforced the natural boundary conditions.But, by the definition of Volterra [23]  as claimed in (46).

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 µ = λ = 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 [63].
The vanishing of viscous forces affects the residuals of the momentum balance equations: R i [u j , p, ρ; x k , t] ≡ ρ ui + ρu i,j u j + p ,i − ρb i = 0, (B.1) where we are using a bar to distinguish the inviscid residuals (R i ) from their more general forms given in (21).Applying the equation of state ρ = ρ(p), we obtain R i [u j , p; x k , t] ≡ ρ ui + ρu i,j u j + p ,i − ρb i = 0. (B.2) The residual 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 where in the latter case we have δS * /δp = 0, since π 4 vanishes identically for incompressible flow.
The second-order formulation described above is fundamentally different from-though again still mathematically equivalent to-the classical Hamiltonian formulation of the first-order Euler equations [63].
functional derivatives, we have that In this way, we may identify the conjugate momenta with the functional derivatives of S * :Here the integral is simply the Hamiltonian H * , with the conjugate momenta replaced by the functional derivatives in accordance with (A.5).Hence we arrive at the Hamilton-Jacobi equation ρR i u j + ρR j u j,i + ρ′ R 4 p ,i −δp : ρ′ R i ui + ρ′ R i u i,j u j − R i,i − ρ′ R i b i + ρ′′ R 4 ṗ ′ R 4 ) + ρ′′ R 4 p ,i u i − ∂ ∂x i (ρ ′ R 4 u i ) + ρ′ R 4 u i,i = 0. (B.5)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π i ≡ ρR i , (B.6)and the momentum π 4 conjugate to the pressure p remains as it was in the main text.For compressible flows, the resulting Hamiltonian density takes the formH * [u i , p, π j , π 4 ; x k , t] = 1 2 1 (ρ) 2 π i π i − 1 ρ ρu i,j u j + p ,i − ρb i π i ′ ) 2 π 4 π 4 −1 ρ′ ρ′ p ,i u i + ρu i,i π 4 , ρ′ = 0.
4 , (B.3) and the associated Euler-Lagrange equations are as follows: