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 ui(xj, t)and pressure p(xj, t)as the field quantities to be varied, along with canoni- cally conjugate momenta deduced from the analysis. From these, a conserved Hamiltonian func- tional 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 ( ui,p) to that of finding a single scalar functional in those fields—Hamilton’s principal functional S∗[ui, p, t]. 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. 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 [1–3]. There is no action integral for which Hamilton’s principle [4–6] yields the Navier-Stokes equations [7–13] in their usual form [1–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. 1 1.1 A motivating example Consider the first-order initial-value problem ˙v=−v, v (0) = 1 , (1) 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–13], this too is an intrinsically non-Hamiltonian problem, in that there is no actionSfor 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 ¨v=v, v (0) = 1 ,˙v(0) =−1, (2) which has the same unique solution v(t) =e−tbut which isHamiltonian—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–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: S∗[v] =ˆ dt1 2R2 =ˆ dt1 2 ˙v2+ 2v˙v+v2 ∼ˆ dt1 2 ˙v2+v2 , (3) where we have used the fact that 2v˙v=d(v2)/dtis 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–18]: since R= 0is a local minimum of R2, it is also a local minimum of´ dt(R2). Varying v, the first variation of S∗is δS∗=ˆ dt ˙vδ˙v+vδv =ˆ dt (−¨v+v)δv + ˙vδvt2 t1, (4) yielding the second-order equation ¨v=vand 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 H∗[v, π] =π˙v−1 2 ˙v2+v2 =1 2 π2−v2 . (5) Notably, this Hamiltonian has nothing to do with the total mechanical energy of the system, al- though it isa 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 ˙v=∂H∗ ∂π,˙π=−∂H∗ ∂v(6) 2 are mathematically equivalent to the second-order problem ¨v=vand therefore also mathemati- cally equivalent to the original, first-order problem. The associated Hamilton-Jacobi equation [4–6, 19–22] is 1 2∂S∗ ∂v2 −1 2v2+∂S∗ ∂t= 0, (7) 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)v2—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 S∗(v, t) =F(t)v+1 2v2+f(t), (8) where F(t)andf(t)are as yet undetermined functions of t. This trial solution was chosen to cancel the term (1/2)v2from the equation. Substituting our trial solution into the Hamilton-Jacobi equation, we find that 1 2[F(t)]2+ [F(t) +F′(t)]v+f′(t) = 0 . (9) In order for this equation to hold for all v, we must have the following: F(t) +F′(t) = 0 ⇒ F(t) =αe−t, (10) where αis a constant of integration which will be used to transform to the new coordinate, and 1 2[F(t)]2+f′(t) = 0 ⇒ f(t) =1 4α2e−2t+γ, (11) where γis another constant of integration which is simply additive and can therefore be discarded. In this way, we have that S∗(v, t;α) =αe−tv+1 2v2+1 4α2e−2t. (12) 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 ϕ=∂S∗ ∂α=e−tv+1 2αe−2t. (13) The numerical value of αis in turn obtained via the canonical relation π=∂S∗ ∂v=αe−t+v, (14) 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 e−tv−e−2t=v(0)−1 = 0 , (15) 3 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–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–13]. 1.2 The Navier-Stokes problem The incompressible Navier-Stokes equations [7–13] are given by ρ˙ui+ρui,juj+p,i−µui,jj−ρbi= 0, (16) ui,i= 0, (17) where ρis the constant and uniform density, ui=ui(xj, t)is the velocity field, p=p(xj, t)is the pressure field, bi=bi(xj, t)is the body force field, subscript Roman indices label Euclidean tensor components ( i, j= 1,2,3), the xjare Eulerian spatial coordinates, tis time, µis the dynamic viscosity, a dot over a symbol denotes a partial time derivative ( ˙ui=∂ui/∂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 ui(xj, t)signifies that each component of the velocity field is a function of all three spatial coordinates (x1, x2, x3)and time t. It has the same meaning as other common notations, such as ui(x, t)andu(x, t). Likewise for all other field quantities. In the case of a uniform gravitational field, bi=gicoincides with the local acceleration due to gravity; however, in what follows, we make no assumptions about the functional form of bi(xj, t): it is completely arbitrary. There are four unknown field quantities: ui(xj, t)andp(xj, t). We seek, ultimately, a functional H∗=H∗[ui, p, π j, π4;t], (18) where ( πi, π4) are suitable “momenta” conjugate to the field quantities ( ui, p), such that Hamilton’s canonical equations ˙ui=δH∗ δπi, ˙p=δH∗ δπ4, (19) ˙πi=−δH∗ δui, ˙π4=−δH∗ δp, (20) constitute a mathematically equivalent second-order formulation of the problem, where δH∗/δui, δH∗/δp,δH∗/δπi, and δH∗/δπ 4are the V olterra [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 ˙p=δH∗/δπ 4will need to be replaced by the incompressibility condition ui,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]). 4 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, cul- minating 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–6, 19–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. 2 Literature review The field of analytical mechanics, with foundations planted in Hamilton’s principle of stationary action [4–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–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–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–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 formu- lation is expressed as an initial-boundary-value problem: a set of coupled, nonlinear partial dif- ferential 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–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–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–58] 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 [59], Benjamin [60], and Holm et al. [61] have all described variational and Hamiltonian formulations of incompressible, inviscid fluid flow. Roberts [62] presented a Hamil- tonian 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 5 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–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–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 im- posed on these velocity components such that all of the terms vu,2,wu,3,wv,3,uv,1, uw,1,vw,2disappear from their positions in the Navier-Stokes equations [1]. (It should be noted that the six terms referred to above come from the convective acceleration ui,juj, and Millikan [1] uses the notation u=u1,v=u2, and w=u3.) 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 [67]. Within the last 80 years, many alternative methods have been developed in an attempt to circumvent the non-self-adjointness of dissipative systems [2, 56, 57, 68–81]. The mathematical study of alternative variational methods as applied to the Navier-Stokes equations in particular remains an ongoing endeavor [18, 82–94]. 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 vec- tor 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. 6 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 generaliza- tion of the Kutta condition to smooth geometries. It should be noted that Gauss’s principle [96] is fundamentally different from Hamilton’s principle [4–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–18] has shown that higher-order dy- namics 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 sys- tem’s resonant frequencies [17]. Higher-derivative theories had been studied before in the realm of quantum gravity physics [97–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 for- malism [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 exten- sive, 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 (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 Ri[uj, p, ρ;xk, t]≡ρ˙ui+ρui,juj+p,i−µui,jj−(µ+λ)uj,ji−ρbi= 0, (21) R4[uj, ρ]≡˙ρ+ρ,iui+ρui,i= 0, (22) where ρ=ρ(xj, t)is the density field (now one of the unknown field quantities along with uiand 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 non- dimensional (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 ρ0U2, with ρ0andUappropriate 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 7 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 ad- ditional thermodynamic variables, such as temperature and enthalpy or entropy. Two of the ther- modynamic 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 equa- tions for the six unknowns fields (ui, 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= ˙ρ+ρ,iui= 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 ρ=ρ0as a system parameter. Consequently, (22) reduces to ρui,i= 0 and the energy equation is decoupled from the system. Accordingly, in the incompressible limit there are only four unknown field quantities (ui, 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 re- spect to the field quantities uiandρ. 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 Sfor which Hamilton’s principle (δS= 0) yields the first-order field equations [1–3]. Finally, we note that in the incompressible limit,R4becomes independent of ˙ρ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 [1–3], nev- ertheless 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 (ui, p, ρ)for which the following action achieves a local minimum: S∗[uj, p, ρ] =ˆ d4x1 2RiRi+1 2R4R4 , (23) where d4x=dx1dx2dx3dt, and the integral is carried out over both the control volume Voccupied by the fluid ( xj∈ V) and the time period of interest ( t∈[t1, t2]). 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 ( Ri,R4). Without an equation of state relating ρtop, 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 uiandp, 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: Ri[uj, p;xk, t]≡ˆρ˙ui+ ˆρui,juj+p,i−µui,jj−(µ+λ)uj,ji−ˆρbi= 0, (24) 8 R4[uj, p]≡ˆρ′˙p+ ˆρ′p,iui+ ˆρui,i= 0, (25) where ˆρ′=dˆρ/dp. We note that, under equilibrium conditions, ˆρ′is related to the speed of sound cand the bulk modulus Kasˆρ′= 1/c2=ρ/K (for incompressible fluids, ˆρ′≡0and the speed of sound and bulk modulus are both infinite). Having specified ˆρ(p),µ,λ, and bi(xj, t), and having prescribed appropriate auxiliary conditions (initial and boundary conditions), one seeks the four field quantities (ui, 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 ( Ri,R4) vanish for the actual motion, they are nottrivially zero. That is, the residuals only vanish for the particular fields (ui, p)which satisfy the first-order field equations (24) and (25); they do not vanish for every conceivable (ui, p). Thus it is not appropriate to take Ri≡0,R4≡0. We will return to this point later when we discuss the Hamiltonian formulation of the problem. For now, we note that the action S∗=S∗[ui, p]defines a Lagrangian L∗[ui, p;t] =ˆ d3x(L∗), (26) where the integral is carried out over the volume Vonly (d3x=dx1dx2dx3), with Lagrangian density L∗[uj, p;xk, t] =1 2RiRi+1 2R4R4. (27) Because the residuals ( Ri,R4) 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 (ui, p)for which S∗achieves a local minimum. To obtain the Euler-Lagrange equations, the conjugate mo- menta, and the natural auxiliary conditions, we insist that S∗not vary to first order ( δS∗= 0) under small variations in the fields (δui, δp). Evaluating δS∗, integrating by parts, and collecting like terms, we find that δS∗=ˆ d4x −∂ ∂t(ˆρRi)−∂ ∂xj(ˆρRiuj) + ˆρRjuj,i−µRi,jj (28a) −(µ+λ)Rj,ij+ ˆρ′R4p,i−∂ ∂xi(ˆρR4) δui (28b) + ˆρ′Ri˙ui+ ˆρ′Riui,juj− R i,i−ˆρ′Ribi+ ˆρ′′R4˙p (28c) −∂ ∂t(ˆρ′R4) + ˆρ′′R4p,iui−∂ ∂xi(ˆρ′R4ui) + ˆρ′R4ui,i δp (28d) +ˆ d3x ˆρRiδui+ ˆρ′R4δpt2 t1(28e) +ˆ d2xdt ˆρRiujnj+µRi,jnj+ (µ+λ)Rj,inj+ ˆρR4ni δui (28f)  −µRinj−(µ+λ)Rjni δui,j+ Rini+ ˆρ′R4uini δp , (28g) 9 where the purely volumetric integral (d3x) is carried out over V, the surface integral (d2x) is carried out over the boundary ∂V, and niis the unit outward normal vector to ∂V. Note that, because we are using Eulerian coordinates xj, the volume element d3xis not to be varied. The Euler-Lagrange equations (which hold for all xj∈ V ) may be read directly from the spacetime (d4x) integral: δui:−∂ ∂t(ˆρRi)−∂ ∂xj(ˆρRiuj) + ˆρRjuj,i−µRi,jj −(µ+λ)Rj,ij+ ˆρ′R4p,i−∂ ∂xi(ˆρR4) = 0 , (29) δp: ˆρ′Ri˙ui+ ˆρ′Riui,juj− R i,i−ˆρ′Ribi+ ˆρ′′R4˙p −∂ ∂t(ˆρ′R4) + ˆρ′′R4p,iui−∂ ∂xi(ˆρ′R4ui) + ˆρ′R4ui,i= 0. (30) 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 Ri= 0andR4= 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 uiare πi≡ˆρRi, (31) and the momentum conjugate to the pressure pis π4≡ˆρ′R4. (32) In the forthcoming Hamiltonian formulation, the conjugate momenta will be used to eliminate the (partial) time derivatives ( ˙ui,˙p)of the field quantities from the Hamiltonian. In general, Hamilton’s principle would insist that the variations (δui, δp)vanish at the endpoints t=t1andt=t2to ensure that the purely volumetric integral (28e) vanishes identically. Interestingly, for the actual motion ( Ri= 0,R4= 0), the volumetric integral (28e) already vanishes even without taking (δui, δp)to vanish at t1andt2. 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 Pi=ρui, which we will see presently from the natural boundary conditions. These are read directly from the surface (d2x) integral: δui: ˆρRiujnj+µRi,jnj+ (µ+λ)Rj,inj+ ˆρR4ni= 0, (33) δui,j:−µRinj−(µ+λ)Rjni= 0, (34) δp:Rini+ ˆρ′R4uini= 0. (35) 10 This last condition, (35), establishes a connection between the new conjugate momenta and the conventional linear momenta. Multiplying (35) by ˆρ, and noting that ˆρui=ρui=Pi, we find that (πi+π4Pi)ni= 0. (36) Evidently, boundary condition (35) states that the flux of the vector Πi≡πi+π4Pithrough the boundary ∂Vshould vanish. It is interesting that this new vector Πicontains both old and new momenta, with π4“carried” ( i.e., given direction) by Pi. 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 (ui, p). In other words, with identical auxiliary conditions, (ui, p)is a solution to the first-order formulation if and only if the same (ui, 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 Ri(xk,0) = 0 andR4(xk,0) = 0 for all xk∈ V ∪ ∂V, along with boundary conditions making Ri(xk, t) = 0 ,Ri,j(xk, t) = 0 , andR4(xk, t) = 0 for all xk∈∂Vand all times t. By supposition, the auxiliary conditions applied to the two formulations are identical, so it suffices to show that (ui, p)satisfies the governing field equations (24), (25) of the first-order formulation ( Ri= 0andR4= 0) everywhere in Vand at all times if and only if (ui, p)satisfies the Euler-Lagrange equations (29), (30) of the second-order formulation everywhere in Vand at all times. Suppose first that (ui, p)satisfies the governing field equations (24), (25) of the first-order formulation everywhere in Vand at all times. Then Ri= 0,R4= 0, and (ui, p)is a trivial solution to the Euler-Lagrange equations (29), (30) of the second-order formulation. Conversely, suppose that (ui, p)satisfies the Euler-Lagrange equations (29), (30) of the second- order formulation everywhere in Vand at all times. We note that (Ri= 0,R4= 0) constitutes an equilibrium solution of the Euler-Lagrange equations (29), (30). Thus, because the initial condi- tions are chosen such that Ri(xk,0) = 0 andR4(xk,0) = 0 for all xk∈ V ∪ ∂V, and because the boundary conditions are chosen such that Ri(xk, t) = 0 ,Ri,j(xk, t) = 0 , andR4(xk, t) = 0 for all xk∈∂Vand all times t, thenRiandR4will remain identically zero everywhere in Vfor all future times. Thus, (ui, p)satisfies the governing field equations (24), (25) of the first-order formulation everywhere in Vand 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 H∗=ˆ d3x(H∗), (37) 11 with the Hamiltonian density H∗obtained from the Lagrangian density L∗via the Legendre trans- form H∗=πi˙ui+π4˙p− L∗=πi˙ui+π4˙p−1 2RiRi−1 2R4R4. (38) Again, this H∗has nothing to do with the total mechanical energy of the system, although it isa 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 ˙uiand˙p. Observe that Ri=πi/ˆρ, and ignoring for the moment the incompressible limit, we may write R4=π4/ˆρ′(ˆρ′̸= 0). In this way, using the functional expressions for the residuals given by (24) and (25), we find that ˙ui=1 (ˆρ)2πi−1 ˆρ ˆρui,juj+p,i−µui,jj−(µ+λ)uj,ji−ˆρbi ; (39) ˙p=1 (ˆρ′)2π4−1 ˆρ′ ˆρ′p,iui+ ˆρui,i ,ˆρ′̸= 0; (40) and H∗[ui, p, π j, π4;xk, t] =1 21 (ˆρ)2πiπi−1 ˆρ ˆρui,juj+p,i−µui,jj−(µ+λ)uj,ji−ˆρbi πi +1 21 (ˆρ′)2π4π4−1 ˆρ′ ˆρ′p,iui+ ˆρui,i π4,ˆρ′̸= 0. (41) Hamilton’s equations [5, 6] are as follows: ˙ui=δH∗ δπi, ˙p=δH∗ δπ4, (42) ˙πi=−δH∗ δui, ˙π4=−δH∗ δp, (43) where δH∗/δui,δH∗/δp,δH∗/δπi, and δH∗/δπ 4are the V olterra [23] functional derivatives of H∗with respect to the field quantities and the conjugate momenta. Equations (42) reproduce (39) and (40), respectively. Equations (43) in turn reproduce the Euler-Lagrange equations (29), (30) 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 (ui, p)that satisfy the governing field equations (24), (25) of the first-order formulation, it does not vanish for every conceivable (ui, p). The latter would imply, according to Equations (42), that ˙ui≡0and˙p≡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. 12 To obtain the Hamiltonian for incompressible flow, we set ˆρ′≡0from the beginning, in which caseR4reduces to ˆρui,iandπ4vanishes identically, consistent with the fact that R4becomes independent of ˙p. The Hamiltonian density in turn reduces to H∗=πi˙ui−1 2RiRi,ˆρ′≡0, (44) or, in terms of the conjugate momenta, H∗[ui, p, π j;xk, t] =1 21 ρ2πiπi−1 ρ ρui,juj+p,i−µui,jj−ρbi πi,ˆρ′≡0, (45) where the density ˆρ=ρis a constant and we have used the fact that ui,i= 0. Hamilton’s equa- tions ˙ui=δH∗/δπi,˙πi=−δH∗/δui, and 0≡˙π4=−δH∗/δpstill apply (and reproduce the corresponding equations in the incompressible limit), but ˙p=δH∗/δπ 4must be replaced by the constraint that ui,i= 0. That the incompressibility condition should take the place of the pressure equation ˙p=δH∗/δπ 4is consistent with the well known result that the pressure usually serves as the Lagrange multiplier for the incompressibility constraint [19, 24]. 4.1 Hamilton-Jacobi equation One of the most significant aspects of the Hamiltonian formalism is that it leads to the transforma- tion theory of Hamilton and Jacobi [4–6, 19–22], celebrated both for unifying particle mechanics with wave optics [4] and for its relationship to the Schr ¨odinger equation of quantum mechan- ics [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∗[ui, 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: H∗ ui, p,δS∗ δuj,δS∗ δp;t +∂S∗ ∂t= 0, (46) where δS∗/δuiandδS∗/δpare the V olterra [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 ˆ d3x1 21 (ˆρ)2δS∗ δuiδS∗ δui−1 ˆρ ˆρui,juj+p,i−µui,jj−(µ+λ)uj,ji−ˆρbiδS∗ δui +1 21 (ˆρ′)2δS∗ δpδS∗ δp−1 ˆρ′ ˆρ′p,iui+ ˆρui,iδS∗ δp +∂S∗ ∂t= 0,ˆρ′̸= 0. (47) 13 In contrast to the four original field equations—(21) and (22)—the Hamilton-Jacobi equation (47) is asingle equation in Hamilton’s principal functional S∗. This constitutes an equivalent formula- tion 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∗≡0for 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 (ui, p). The case of incompressible flow requires care, as ˆρ′≡0andˆρ′appears in the denominators of terms in (47). Even so, the Hamiltonian formulation remains well posed in the incompressible limit. Recall that, with ˆρ′≡0,R4reduces to ˆρui,i,π4vanishes identically, and the Hamiltonian density reduces to (45). Substituting for the conjugate momenta πiin (45), the corresponding Hamilton-Jacobi equation is ˆ d3x1 21 ρ2δS∗ δuiδS∗ δui−1 ρ ρui,juj+p,i−µui,jj−ρbiδS∗ δui +∂S∗ ∂t= 0,ˆρ′≡0, (48) with δS∗/δp= 0, since again π4vanishes 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∗/δpshould vanish in the incompressible limit without knowing that in general π4= ˆρ′R4. 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 ui,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 V olterra [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–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. 5 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 RiandR4are substituted. Our motivation is again the simple example of Section 1.1 for which the first-order non- Hamiltonian system ˙v=−vwas converted to the second-order Hamiltonian system ¨v=vby 14 (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 inte- gral 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: − R i,i= 0∀xj∈ V,subject to Rini= 0∀xj∈∂V. (49) This higher-order field equation is simply the divergence of the residual Ri. Upon substituting for Rifrom (21) and subsequently imposing the incompressible continuity condition ui,i=R4/ρ= 0 from (22), we obtain: p,ii=−[ρujui,j],i+ρbi,i, (50) which is a Poisson equation for the pressure. The boundary condition is a Neumann type requiring the specification of the normal pressure gradient, nip,i=p,n≡f(xj, t), on the boundary: f(xj, t) =−ni ρ˙ui+ρujui,j−µui,jj−ρbi . (51) 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 δuiandδui,jvariations. The incompressible versions of these equations are: ρRiujnj+µRi,jnj= 0 and−µRinj= 0∀xj∈∂V. (52) The boundary conditions involving the residual Riare 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 Ri≡0there. Note that the pressure of the actual motion on the boundary will be known from the simultaneous solution of (50). However, the gradient termsRi,jwill 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 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 15 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 incom- pressibility and expanding derivatives of product terms yields: ρ˙Ri+ρujRi,j=ρRjuj,i−µRi,jj∀xj∈ V. (53) The left-hand side is the material derivative of the residual Ri. Our purpose here is to observe which terms from the first-order Navier-Stokes equation are ‘eliminated’ in the higher-order for- mulation. Specifically, we are interested in the non-conservative viscous terms; while the body force bicould also be non-conservative, we will not concern ourselves with this possibility. Direct substitution of Riinto (53) generates many terms, but it is found that only one is canceled: the vis- cous Laplacian of the (time derivative of the) velocity, namely µ˙ui,jj. This term mutually appears from the ρ˙Riand−µRi,jjterms in (53). To maintain notional clarity, we write the residual as: Ri=ρ˙ui−µui,kk+˜Ri, (54) where index khas been used to avoid confusion with gradient operators in (53) having index j, and˜Ri=ρukui,k−ρbiare the remaining terms in the residual. Substituting the above into the first and last terms of (53), canceling the aforementioned µ˙ui,jjterm, and then dividing out by the density gives: ρ¨ui+˙˜Ri+ujRi,j=Rjuj,i−ν −µui,kkjj +˜Ri,jj , (55) 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. 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 (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 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. 16 Here we will examine the case in which there is motion only in the x1direction, and the velocities take the form {ui}={u1(x2, t),0,0}. In the absence of a body force, our pressure gradient in the x1direction is solely a function of time and the pressure gradients in the x2andx3 directions are zero. There are thus only two unknown field quantities: u1(x2, t)andp(x1, t), where pis linear in x1. The field equation of primary interest is R1≡ρ˙u1+p,1−µu1,22= 0, (56) and the remaining field equations are satisfied automatically by the assumed form of the fields. Following the procedure described above, the momenta conjugate to u1andpare given by π1≡ρR1, π 4≡0. (57) This results in a Hamiltonian density given by: H∗=1 21 ρ2π1π1−1 ρ p,1−µu1,22 π1. (58) Hamilton’s principal functional S∗=S∗[u1, p, t]can be expressed as an integral over x2only, 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: ˆ dx21 21 ρ2δS∗ δu1δS∗ δu1−1 ρ p,1−µu1,22δS∗ δu1 +∂S∗ ∂t= 0, (59) withδS∗/δp= 0. The solution to (59) would provide a canonical transformation to a new set of coordinates, giving analytical expressions for (u1, 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–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]. 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 pre- viously 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∗[ui, p, t]directly. If a complete solution for S∗can be found, it will lead via canonical trans- formation 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. 17 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 mathemat- ically 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 assump- tions,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 V olterra [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 cen- tury [107–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