Analytical Solution of the Unsimplified Incompressible Navier-Stokes Equations for Arbitrary Cauchy Data by Second-Curl Poisson Self-Similarity and Fixed-Point Triquartic Decoupling Amid Modified Cole-Hopf Transformation

1. Introduction

Let $\Omega \subseteq {\mathbb{R}}^3$ be simply connected, $(\rho ,\nu )\in {\mathbb{R}}_{>0}^2$, $t\in {\mathbb{R}}_{\ge 0}$, $\mathcal{F}\in \{\mathbf{v},\mathbf{f}\}:\left[(\mathbf{r}\in \Omega )\times t\right]\to {\mathbb{R}}^3$, and $P:\left[(\mathbf{r}\in \Omega )\times t\right]\to \mathbb{R}$. Then, the incompressible Navier-Stokes momentum (1) and continuity (2) equations describe the velocity $\mathbf{v}\left(\mathbf{r},t\right)$ and pressure $P\left(\mathbf{r},t\right)$ of an incompressible Newtonian fluid of density $\rho$ and kinematic viscosity $\nu$ [1]. External forces $\mathbf{f}\left(\mathbf{r},t\right)$ are given while $\mathbf{v}\left(\mathbf{r},t\right)$ and $P\left(\mathbf{r},t\right)$ are sought.

$${\displaystyle \frac{D\mathbf{v}}{Dt}}={\displaystyle \frac{\partial \mathbf{v}}{\partial t}}+\left(\mathbf{v}·\nabla \right)\mathbf{v}=\mathbf{f}-{\displaystyle \frac{\nabla P}{\rho }}+\nu {\nabla }^2\mathbf{v}$$

(1)

$$\nabla ·\mathbf{v}=0$$

(2)

Though the powerful generality of (1-2) exemplified in both their virtually limitless degree of freedom to replicate incompressible fluid motion and their elegant derivation from first principles of Newtonian mechanics has cemented (1-2) as the central mathematical figure of fluid mechanics, they remain analytically unsolved for general $\{\mathbf{v}\left(\mathbf{r},t\right),P\left(\mathbf{r},t\right),\mathbf{f}\left(\mathbf{r},t\right)\}$. Analytical solutions are only known for very simple flow problems which exhibit negligible effects of nonlinear advective term coupling and are derived by imposing many simplifying assumptions.

2. Literary review

Analytical closure of (1) and (2) is an active and ongoing field of research spanning two centuries. Due to the emphasis of this work on a method of dealing with a model to increase its tractability rather than questioning the model itself, the literary record is divided into five eras spanning roughly half a century apiece with each focusing on the prevailing philosophy regarding the study of (1) and (2).

2.1. The Discovery Era (c. 1750-1850)

The discovery of the Navier-Stokes equations is considered to have occurred in five independent, effectively contemporaneous discoveries spread over the course of less than only three decades of the early 19th century. The discoverers in chronological order were Claude-Louis Navier in 1822, Augustin-Louis Cauchy in 1823, Siméon Denis Poisson in 1829, Adhémar Jean Claude Barré de Saint-Venant in 1837, and George Gabriel Stokes in 1845 [2,3,4]. The first attempt at analytically solving a specific case of (1-4) was most likely the Hagen-Poiseuille flow [5] which sought to model viscous laminar flow through a round pipe and resulted in a parabolic velocity profile taking as parameters the kinematic viscosity, pipe radius, pipe length, and pressure differential across the pipe length. Attempts to expand Hagen-Poiseuille flow into the more general context of turbulent pipe flow were roadblocked due to the mathematical difficulty increasing at an alarming rate as more information was included in analyses. The original purpose behind these discoveries was to contrive a way for the Euler equations written earlier in 1755 [6] to account for viscosity, which had included all terms of (1) and (2) except the viscous term $\nu {\nabla }^2\mathbf{v}\left(\mathbf{r},t\right)$. The Euler equations had become well known to mathematicians and physicists in the late 18th to early 19th centuries, and it is said that the Euler equations were the first historically documented set of partial differential equations. Even the metaphysical idea of describing the movement of something continually yielding under applied shear stress by the unyielding, shear-resistant truth of mathematical coherence was highly innovative for its time, and this idea and the equations born thereof remain an innovation of Euler that unfortunately seems to be one of his lesser known achievements by the public. Despite the historic significance of the Euler equations as the Charlemagne one could write on an index card that fathered the epistemological dynasty of PDEs as mathematical encodings of cosmic phenomena, they fall short of providing a complete model of Newtonian fluids since they contain no viscous term to slow the fluid to a stop upon cessation of its motive forces. With the Navier-Stokes equations reaching their standardized form in the mid-19th century, the idea of an analytical interpretation of fluid mechanics as the Newtonian laws of motion applied to substances lacking shear modulus was rendered sufficiently encoded in mathematical language to be counted a field worthy of deep mathematical study.

2.2. The Parametric Era (c. 1850-1920)

The most notable findings of the late 19th century regarding the mathematics of fluid mechanics were likely that of Osborne Reynolds, who realized that injecting a thin line of dye into a glass water pipe allowed the empirical observation of flow patterns when the observation of a free surface was impossible or insufficient for the desired level of understanding. When the average magnitude of the inertial forces within the flow was more than 2200 times greater than that of the viscous forces, the dye trail ceased to be a straight line and began to swirl to and fro until it diffused into the surrounding fluid. Regardless of the pipe dimensions, this laminar-to-turbulent transition consistently began to occur at an inertia-viscosity ratio of roughly 2200 and was completed on the order of 5000 [7]. Various figures are given in the literature since the transition is understood to occur in a certain parametric range as opposed to a known exact theoretical value. The idea of using a dimensionless inertia-viscosity ratio to quantify how turbulent a flow is caught on surprisingly quickly as the Reynolds number and was observed to have different critical value ranges for different flow geometries, and even for different dimensions of some geometries, since the inertia calculation requires the setting of a length scale. For example, flow over a cuboid body with breadth at least ten times less than its length or width (a flat plate in the literature) has a critical Reynolds number of approximately $5.0\times {10}^5$ [8], nearly 230 times the lower end of the critical Reynolds number regime of an ideal round pipe. The problem of ideal pipe flow became a classic introductory problem in fluid mechanics around this time, with the Moody chart [9] plotting the relative roughness of a round pipe wall against both the Reynolds number and the Darcy-Weisbach friction factor. The empirical emphasis in the study of fluid mechanics persists today with engineers and scientists being taught to use the Moody chart in the normative course of their undergraduate studies. The Darcy-Weisbach friction factor is used in the Darcy-Weisbach equation, which continues to serve as an effective model of laminar or turbulent flow through a pipe with a specified internal surface roughness which affects the range of the critical Reynolds number regime. The eponymous friction factor is and was normally found through the complementary yet notoriously implicit Colebrook-White equation, which was recently solved in closed form using the Lambert W function [10]. Many similar thermofluidic correlations such as Zukauskas’ correlation [11] and the Blasius boundary layer approximation [12] were discovered in this period. The invention of the airplane at the turn of the 20th century made the pursuit of solutions to the Navier-Stokes equations a highly desirable endeavor even for incompressible flow, though the required mathematics remained poorly understood. The bulk of fluid mechanical knowledge dating to this period was gathered experimentally through wind-tunnel tests and channels, using smoke to highlight turbulent patterns in gases, and dye for liquids building from the Reynolds pipe [13]. Mockups of entire aircraft were constructed to ensure airflow over wings was sufficiently laminar to ensure fuel efficiency and safety from aerodynamic flutter or stall. To assist aerospace engineers in this process, the National Aeronautic and Administration (NACA) compiled their well-known four-digit and five-digit [14] series of airfoil cross-sections mathematically conceived of as polynomial splines fitted to measured fluid mechanical data. Empirical formulae such as the lift and drag equations compressed complicated fluid dependencies into a single parameter known as a lift (or drag) coefficient while the reasons why they took on the values they did remained a mystery. The mystique surrounding such parameters, together with the Reynolds number and a menagerie of other named dimensionless quantities such as the Mach and Prandtl numbers, was sufficiently high that it was commonly taught by mentors that one should not waste time meditating on such things. The parameters worked sufficiently well for their applications, and the requirement of a general and analytical solution to the Navier-Stokes equations to explain said parameters largely served as an “elephant in the room.”

2.3. The Analog Era (c. 1920-1960)

A multitude of lessons in aerodynamics were learned through aerial combat in both World Wars. Higher airspeeds and altitudes meant higher stakes regarding personnel safety and mission reliability as well as higher Reynolds and Mach numbers. Propeller-driven civilian airliners were becoming commonplace at similar speeds and altitudes. Greater amounts of turbulence made it necessary to search for more sophisticated models of fluid mechanics than burying the swirls in fine-tuned parameters. With the war effort having improved punched-card computing technology, the first direct numerical simulation (DNS) algorithms were published for implementation into the analog computers of the day, not only in aerodynamics but also in weather forecasting [15] when Lorenz discovered the iconic Lorenz attractor while formulating a meteorological model [16], raising questions of how chaos theory could be exploited. Both the incompressible and compressible Navier-Stokes equations became capable of being studied only in theoretical contexts. In 1958, Olga Ladyzhenskaya published a proof [17,18] that smooth and globally defined solutions to the incompressible Navier-Stokes equations always exist in two space dimensions, showing also an inequality that a given velocity field must satisfy to qualify as a solution. Ladyzhensakya’s inequality was also written for three-dimensional flow, with the insightful observation that exponents on two instances of the Euclidean norm became ${\displaystyle \frac{1}{4}}$ and ${\displaystyle \frac{3}{4}}$ in three-dimensional flow rather than them both equaling ${\displaystyle \frac{1}{2}}$ in two-dimensional flow. It was later discovered that Ladyzhensakya’s inequality was a special case of the Gagliardo-Nirenberg interpolation inequality [19], used more generally in functional analysis to help define Sobolev spaces of elliptic PDE solutions. Consequently, much of the subsequent effort to prove desirable information about the Navier-Stokes equations such as existence, smoothness, and uniqueness has centered around the study of its Sobolev solution spaces, while the exponential mismatch in Ladyzhenskaya’s inequality being often cited as one of the largest epistemological roadblocks. Perhaps the most significant closed-form solution breakthrough in this era was the solution of the viscous Burgers equation by Julian Cole [20] and Eberhard Hopf [21]. The Burgers equation was introduced by Jan Burgers [22] as a simplified analytical model of compressible flow containing only one nonlinear term and no pressure or force terms, sufficient to begin analytical attacks on the Navier-Stokes equations, but in no way able to achieve even a two-dimensional solution. This start bore rather early fruit in what would become known as the Cole-Hopf transformation in which an arbitrary change of variables was introduced into a class of second-order PDEs of which the viscous Burgers equation was a special reduced case. The change of variables was then chosen so as to force reduction to a heat equation, which has had well-known analytical solutions since the time of Fourier. The reduction to a heat equation allowed the solution thereof to be substituted directly into the change of variables that caused the heat equation to arise, and the inversion of said mapping allowed initial and boundary conditions to be swiftly and smoothly reconciled between those of the heat equation and those of the viscous Burgers equation. Even though there was relatively little analytical study of the Navier-Stokes equations in industry settings during this period, the prevailing kingship of the parametric interpretation of fluid mechanics was insufficient to stop academia and industrial dissenters alike from pursuing analytical solutions, and some of the best-known analytical solutions of the Navier-Stokes equations date to this time. The best-known are the Landau-Squire jet [23,24] describing a steady jet of fluid colliding with a flat perpendicular obstruction, the Taylor-Green vortex [25] describing a tessellated pattern of vortices decaying to constant velocity, and the Gromeka-Arnold-Beltrami-Childress flow [26] describing highly chaotic flow trajectories with relatively simple formulae under an inviscid assumption. Many classifications of flows were grouped and named during this period, such as Beltrami flows whose Lamb vector $(\nabla \times \mathbf{v})\times \mathbf{v}$ [27] vanishes everywhere. Potential flows were widely studied as well during this time, whose velocities exist as the gradient of a scalar field, but cannot describe vortex stretching, vortex shedding, or any other energy trade routes between off-diagonal advective terms. Though potential flows thus fall short of describing turbulence, they have been pivotal in studying laminar flow around and through various geometries under inviscid or irrotational approximations since vector fields that are both irrotational and incompressible or are effectively inviscid are analytically describable as solutions to the relatively simple and well-understood Laplace equation. Attempting to generalize this framework to rotational, turbulence-bearing flows resulted in a dramatic increase in energy trade between velocity components and their partial derivatives and thus became much more difficult to address. The framework in which the rotational and irrotational parts of the flow would be separated from one another was nonetheless proven useful in specific flow problems and became known as the Helmholtz decomposition [28].

2.4. The Statistical Era (c. 1960-2000)

Ladyzhenskaya, Burgers, Cole, Hopf and the others from the previous era would see their labors bear fruit during the age of digitalization brought about by the invention of the both the transistor and the integrated circuit. The triumph of human spaceflight, exemplified in its highest degree in achieving physical human presence on the Moon, required such a strong drive of innovation that computing technology improved drastically enough to enable viable development of the first statistical models of turbulence, the mathematical backbone of modern CFD methods. Much of the knowledge that turbulence is a dissipative phenomenon at that time and before was intuitive, evidenced by Lewis Fry Richardson writing the rhyme “Little whirls have smaller swirls, with differing velocities; little whirls have smaller swirls, and so on to viscosity” in the 1930s [29]. The idea behind this intuition was given statistical formalism by Andrey Kolmogorov [30], who discovered that the size of the smallest eddies in a turbulent flow were predictably proportional to the fluid density and viscosity. The length scale spectra formed when controlling for different flow variables became known as the Kolmogorov microscales, and the flow variable of turbulent kinetic energy (TKE) became an important variable in computational studies since it had become more widely known from Kolmogorov and his contemporaries that the kinetic energy released as heat by the dissipation of eddies was a significant phenomenon, and many began to theorize that the thermodynamic purpose of turbulence is to serve as a kinetic heatsink. With the innovations in computing technology in full swing in the 1970s and 1980s coupled with the unique closed-form solution to the second-order nonlinear flow problem of the Burgers equation having been successfully found, the aerospace industry theorized that transonic airfoils could be redesigned to combat shock waves aft of their upper cambers. Airfoil optimization had become an important concern for marketability since the normative cruising speed of a commercial jet airliner had stabilized at roughly 560 mph, sufficiently close to the sound barrier to cause transonic shock waves. New aircraft designs featured supercritical airfoils in light of deeper study into transonic phenomena. In an attempt to account for more information and push for greater generality in analytical methods, the far more complicated problem of solving the Burgers equation for arbitrary forcing had become paramount since pressure and density would then be connected to a continuity equation and an equation of state, enabling a full theory of one-dimensional compressible flow. Attempting to use the Cole-Hopf transformation resulted in a Schrodinger-like equation [31,32] which likewise remained unsolved, and the search exhausted nearly all known methods until the works of Sophus Lie in group-theoretic topology underwent a two-decade renaissance. The concepts of Lie point and group symmetries were exploited to simplify various classes of forced Burgers equations, and a flurry of research ensued [33] presenting closed-form solutions to large categories of forced Burgers equations. Towards the end of this renaissance, it was discovered [34] that the unsolved variable-coefficient Riccati equation, a first-order ODE quadratic in its unknown [35], was the fundamental obstruction to achieving general forced Burgers solution, and said solution became contingent on the solubility of the Riccati equation hidden within the Lie symmetries. It was in this time that the more well-known numerical solution algorithms such as k-$ϵ$ [36], k-$\omega$ [37], LES [38], RANS [39], SIMPLE [40], and PISO [41] ubiquitous in CFD software today were formulated, enabled by the pioneering work of Bram van Leer [42], Peter Lax [43], and Sergei Godunov [44] building on the work of Richard Courant, Kurt Friedrichs, and Hans Lewy [45] together with John von Neumann [46].

2.5. The Millennium Era (c. 2000-Present)

The advent of CFD and the optimization of reusable spacecraft in the 1980s and 1990s raised awareness of the vitality of the Navier-Stokes equations to aerospace engineering even further. Researchers such as Antony Jameson developed custom software for highly efficient aerospace applications [47] used at NASA while figures such as Stan Osher and Ami Harten [48] contributed to their widespread implementation. In 2000, the Clay Mathematics Institute announced seven open problems in pure and applied mathematics dubbed the Millennium Problems due to their importance in furthering mathematical knowledge in the 21st century. As of this writing, only one solution has been achieved, namely the Poincaré conjecture by Grigori Perelman [49] who continued the work of Richard Hamilton [50] on the same problem. The question of whether the Navier-Stokes equations possess globally smooth $\mathbf{v}\left(\mathbf{r},t\right)$ and $P\left(\mathbf{r},t\right)$ on ${\mathbb{R}}^3$ or ${\mathbb{T}}^3$, regardless of the value of ${\mathbf{v}}_0\left(\mathbf{r}\right)$ or $P_0\left(\mathbf{r}\right)$ so long as both of the latter are well-defined to a specified degree, is listed as a Millennium Problem [51] and is referred to by several names such as Navier-Stokes continuity and smoothness, Navier-Stokes existence and uniqueness, Navier-Stokes existence and smoothness, or simply the Navier-Stokes problem. This proclamation by the CMI drastically increased public awareness of the equations while stimulating many to attempt the problem. Perhaps the closest anyone has come to a full resolution of the Navier-Stokes Millennium Problem has been a 2015 work by Terence Tao [52] proving that (1-4) are capable of “blowing up” (a term used often in the literature to refer to the value of $\mathbf{v}\left(\mathbf{r},t\right)$ or $P\left(\mathbf{r},t\right)$ diverging or otherwise becoming undefined) in finite time when the fine scale of turbulent structure is averaged to analyze its bulk motion while ensuring that the averaged $\mathbf{v}\left(\mathbf{r},t\right)$ obeys conservation of energy, suggesting that resolution of the Navier-Stokes Millennium Problem requires closer handling of the nonlinear terms of (1-3), possibly with no averaging methods at all. One three-dimensional, time-evolving solution was found through the study of periodicity constraints on the velocity and pressure fields was found by Matteo Antuono in 2020 [53] while several solutions for slightly looser conditions were found by a generalized version of the well-known PDE method of separation of variables and tabulated in one of the many compendiums of Andrei Polyanin et al [54].

2.6. Summary of the Literature

The Navier-Stokes equations, though chaotic in their behavior, do nonetheless possess a particular mathematical structure and appear when numerically simulated to conserve kinetic energy and to a surprisingly great extent have been shown to obey many correlative laws. Many common threads are observed running through groups of such correlations, pointing to deeper connections of which known correlative laws are merely specific examples, as the presence of a spring on one mountain and a few other similar-tasting ones on adjacent hills points to an aquifer system. These facts about the Navier-Stokes equations, taken together, point to a very large and deeply hidden mathematical structure underlying the cacophony of approximative methods which characterize the modern study of fluid mechanics. This work seeks to reveal the incompressible wing of this hidden structure by analytically solving the Navier-Stokes equations from first principles with no quarter given, whether to their inherent difficulty or their reputation as such.

3. Non-Mathematical Reasons for Solution Unknownness

It is beyond dispute that analytical solutions of the unsimplified Navier-Stokes equations for arbitrary flow geometry, initial conditions, and boundary conditions are considered to be mathematically uncharted territory. With the finding of such solutions being a highly daunting task with surprisingly few attempts at full analytical closure relative to the expected impact of such a solution justifying the human effort required of its derivation, the question of why knowledge remains minimal becomes just as important as the mathematical strategy a proposed solution would utilize.

3.1. Economic dependence on Navier-Stokes as an unsolvable algorithm

In almost every practical scenario encountered either academically or industrially involving fluid mechanics, the Navier-Stokes equations are treated numerically rather than analytically using specialized high-performance computing infrastructure. While the advent of high-performance computing systems has spawned a great deal of advancements spanning and connecting diverse fields such as aerospace, meteorology, oceanography, and even cosmology which have added a tremendous amount of value to economies worldwide, this prevailing state of affairs has led to a dynamic in which specialized software-based knowledge has added more value than deep mathematical knowledge regarding the same subject. While these methods have proven extremely useful for fast turnaround times in engineering contexts, the methods themselves support a billion-dollar industry since high-performance computing infrastructure is needed to implement them. Specialists of this field with sufficient knowledge of the Navier-Stokes equations thus understand them as what is to engineering as SHA-256 is to cryptocurrency: the Navier-Stokes equations are seen as too profitable to solve. Even though the act of solving them would unlock even more profitable endeavors and create entirely new industries and thousands of jobs, the fact that an economic shift due to the rapid change in simulation technology would come first is too deep a thought to bear for most prospective human solvers, especially when considering their need to first acquire sufficient mathematical knowledge to understand the problem and its scope, then develop novel methods to solve the problem, and legibly show that their proposed solution is indeed valid, all while working possibly the same role that enabled acquisition of mathematical and fluid mechanical knowledge. Even positions sufficiently close to fluids or aerospace engineering departments could be discouraged from attempting to solve, since conflicts of interest would likely arise. Since the chief objective of most in salaried positions is to earn a living rather than to disrupt their respective industries or institutional roles even when doing so would profoundly advance them, it is rare indeed for such a one to attempt solution, though in most cases only these parties have necessary knowledge to do so. The resultant situation is then perpetuated by these mutually conflicting dependencies, barring an unlikely external intervention.

3.2. Straw Men Are Lazy: Lack of Action Falsely Implied by Lack of Necessity

The question of why the markedly human trait of sheer intellectual curiosity has not resulted in at least one large-scale analytical attack on the Navier-Stokes equations, sans uncertainty-driven fear of changing the status quo at a rate possibly higher than that to which the status quo is capable of adapting without major structural change, finds its most plausible explanation in the equally human phenomenon of socioeconomically justified inaction. Although necessity is indeed the mother of innovation, the lack of immediate, perceived, or popularly foreseeable necessity of the completion of a task is insufficient grounds for willful lack of said action, the discouragement of efforts to act, or the dismissal thereof. In a word, the reality that humanity does not need to solve (1) and (2) analytically does not itself imply that it should not do so. One might even claim that humanity refusing to solve (1) and (2) is unethical since the technological development enabled by a solution could increase the quality of life for many.

3.3. Status Loss Phobia Due to Deeper Philosophical Implications

While the relatively low degree of public knowledge of the Navier-Stokes equations is not the same in those domains that depend on them from the academic to the design engineer, they arise often in contexts which effectively form buffer zones between the natural sciences and the metaphysical. Horace Lamb is famously quoted [55] as having stated, “I am an old man now, and when I die and go to Heaven there are two matters on which I hope for enlightenment. One is quantum electrodynamics and the other is the turbulent motion of fluids. And about the former I am rather more optimistic.” Lamb’s quote, which is often misattributed to Werner Heisenberg, was articulating the prevalent idea that chaos is utterly incomprehensible to humans despite its ubiquity in nature, implying that collecting exhaustive theoretical and practical knowledge of something as chaotic as a turbulent fluid is a task so arduous that one would require divine intervention in order to complete it. Because of the metaphysical significance of deterministically predicting turbulent flow in that there would be an exact, strict, hardsided formula replicating fluid chaos from first principles and thus doing away with the idea of true randomness in the universe, a certain fraction of internal discouragement from attacking the Navier-Stokes equations in pursuit of two explicit, independent, and deterministic formulae giving the velocity and pressure of virtually any liquid at all has likely come from the awareness of the implications of the postmodern taboo on cosmic determinism when the discovery of a general Navier-Stokes solution would result in a philosophical investigation into the problem of randomness that could veer into the metaphysical or even the theological realms. Such a discussion could easily place discoverers at odds with the zeitgeists of nihilism and moral relativism which have demanded high conformity in the developed world. Individual and mass forebodings of philosophical fallout both in academia and the inevitable spread in some form to the industry with the name of a person or group attached have almost certainly played a role in the small number of attempts at analytical closure, especially when compared to the number of professionals with access to sufficient tooling, talent, and knowledge to study the equations. It is highly unlikely that a prospective human solver would willfully and knowingly expend a lifetime of effort risking their public image for a goal whose effect is readily achievable in relevant industry contexts, even though a solution would drastically heighten GDP growth per unit time within those contexts, unless said solver presents their results in such a way as to spur discourse more thoughtful and innovative than ad hominem suppression of theories deemed threatening to enshrined assumptions.

4. Structure of Proposed Navier-Stokes Solution

4.1. Foundation as Three Postulates

The metamathematical and philosophical foundation of the entire solution approach to (1) and (2) is outlined in Postulates 1-3 and is fundamentally rooted in the school of mathematical Platonism.

Postulate 1. 

All quantities in the solution to the unsimplified incompressible Navier-Stokes equations (1) and (2) are capable of being expressed without approximations or numerically computed components, using the loose and informal definition of a “closed form” solution as having every part thereof being fully symbolic, fully nonlinear, and fully evaluable under arbitrary input data.

Postulate 2. 

The mathematical essence of the solution to the unsimplified incompressible Navier-Stokes equations (1) and (2) is a set of two separate and explicit expressions capable of independent investigation and computation, coupled only in their initial and boundary data, one for velocity as a real Cartesian vector field $\mathbf{v}\left(\mathbf{r},t\right)$ of three space dimensions ${\left|\right|\mathbf{r}\left|\right|}^2=x^2+y^2+z^2$ and time t, and the other a real scalar field $P\left(\mathbf{r},t\right)$ also of three space dimensions ${\left|\right|\mathbf{r}\left|\right|}^2=x^2+y^2+z^2$ and time t.

Postulate 3. 

The solution $\left\{\mathbf{v}\left(\mathbf{r},t\right),P\left(\mathbf{r},t\right)\right\}$ to the unsimplified incompressible Navier-Stokes equations (1) and (2) remains mathematically regular, regardless of the presence or absence of manifested chaos or turbulence, or of the degree of chaotic severity.

4.2. Irrotational Flow by Helmholtz Decomposed Potential Flow

The solution of the Navier-Stokes equations proposed herein is presented in three main parts. The first part begins by imposing a Helmholtz decomposition on the velocity field $\mathbf{v}\left(\mathbf{r},t\right)$, separating it into an irrotational part ${\mathbf{v}}_{\mathrm{irr}}\left(\mathbf{r},t\right)$ and a rotational part ${\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)$ and observing the well-known result that irrotational flow requires that ${\mathbf{v}}_{\mathrm{irr}}\left(\mathbf{r},t\right)$ by the gradient of a scalar field. Combining this conservative form of ${\mathbf{v}}_{\mathrm{irr}}\left(\mathbf{r},t\right)$ with incompressibility results in a harmonic ${\mathbf{v}}_{\mathrm{irr}}\left(\mathbf{r},t\right)$, that is, a ${\mathbf{v}}_{\mathrm{irr}}\left(\mathbf{r},t\right)$ satisfying a Laplace equation whose solutions are also well-known. The Dirichlet boundary condition in the potential scalar required by the Laplace equation solution is converted to a boundary condition in ${\mathbf{v}}_{\mathrm{irr}}\left(\mathbf{r},t\right)$ by inversion of the gradient operator. The boundary condition in ${\mathbf{v}}_{\mathrm{irr}}\left(\mathbf{r},t\right)$ is left without conversion to $\mathbf{v}\left(\mathbf{r},t\right)$ until the rotational flow is found to maintain clarity of methodology in later steps.

4.3. Rotational Flow by Vector-Triquartic Modified Cole-Hopf Methods

The second part is the largest and most technically complex since it seeks a general rotational velocity field ${\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)$ and therefore requires direct dealings with the nonlinear advective terms. To find the ${\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)$, the curl of the momentum equations is taken twice, eliminating both ${\mathbf{v}}_{\mathrm{irr}}\left(\mathbf{r},t\right)$ and $\nabla P\left(\mathbf{r},t\right)$ and leaving only that velocity information which is encoded in ${\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)$. Because incompressibility requires that $\nabla ·\mathbf{v}\left(\mathbf{r},t\right)=0$, the identity $\nabla \times \left(\nabla \times \mathbf{v}\left(\mathbf{r},t\right)\right)=\nabla \left(\nabla ·\mathbf{v}\left(\mathbf{r},t\right)\right)-{\nabla }^2\mathbf{v}\left(\mathbf{r},t\right)$ reduces to $\nabla \times \left(\nabla \times \mathbf{v}\left(\mathbf{r},t\right)\right)=-{\nabla }^2\mathbf{v}\left(\mathbf{r},t\right)$, rendering the entire momentum equation a Laplace equation whose unknown is itself the sum of all terms of the momentum equations sans $\nabla P\left(\mathbf{r},t\right)$. Solving this Laplace equation effectively sets the momentum equations equal to a boundary integral containing a boundary condition in the image of the momentum equations, showing that the rotational Navier-Stokes momentum equations possess a hidden forcing term containing their own boundary conditions that occupies what would normally be the place of $\nabla P\left(\mathbf{r},t\right)$. These reduced momentum equations, due to their importance in addressing self-similarity in Navier-Stokes solutions with respect to boundary conditions, are termed the pseudo-depressurized momentum equations. Since addressing their solution requires entirely new mathematical methods, the Cole-Hopf transformation is considered not as the rote recollection of a historically significant result, but as a specific example of a more general process of finding a change of variables reducing an unsolvable PDE problem to a solvable problem such that the solution to the more general PDE once counted unsolvable is recovered in closed form from the reduced problem. The Cole-Hopf transformation was engineered to satisfy the Burgers equation $u_t+u{u}_x=\nu u_{xx}$ by reduction to a solvable heat equation ${\psi }_t(x,t)=\nu {\psi }_{xx}(x,t)$ via $u(x,t)=-2\nu {(ln\left(\psi (x,t)\right))}_x$ not by attempting various substitution until a suitable transformation was discovered, but by introducing an arbitrary transformation $u(x,t)=\Theta \left(\psi \right(x,t\left)\right)$ that, when substituted into the Burgers equation, yielded $\Theta \left(\psi \right)=-2\nu {(ln\left(\psi (x,t)\right))}_x$ and thus $u(x,t)=-2\nu {(ln\left(\psi (x,t)\right))}_x$ when $\varphi (x,t)$ was defined as a solution of the heat equation and the remaining terms treated as an ODE in $\Theta \left(\psi \right(x,t\left)\right)$ and solved for the same. The underlying rhetoric of the derivation of the Cole-Hopf transformation sought to answer the question “What does $\Theta \left(\psi \right(x,t\left)\right)$ need to be in order to solve $u_t(x,t)+u(x,t)u_x(x,t)=\nu u_{xx}(x,t)$ if $\psi (x,t)$ is assumed to satisfy ${\psi }_t(x,t)=\nu {\psi }_{xx}(x,t)$ whose closed-form solution is known?”, as opposed to the subtly less mathematically closed question being asked in this work, “What $u(x,t)=\Theta \left(\varphi \right(x,t\left)\right)$ can be substituted into the PDE $F(u_t(x,t),u_x(x,t),\nu u_{xx}(x,t),u(x,t))=0$ to find $u(x,t)$ in terms of $\psi (x,t)$ if the latter satisfies a PDE $G({\psi }_t(x,t),\nu {\psi }_{xx}(x,t))=0$ with a known closed-form solution?” The essence of this important nuance in the derivation philosophy of the Cole-Hopf transformation is that said transformation may, in principle, be generalized to attack very broad classes of second-order PDEs and PDE systems.

4.3.1. The Modified Cole-Hopf (MCH) Transformation Framework

Because the aforementioned strategy of modifying the Cole-Hopf transformation to accommodate a much broader palette of second-order PDE solutions is used also in this subtask of solving the rotational incompressible Navier-Stokes equations, it is important to formalize this modified Cole-Hopf (MCH) transformation framework for tackling stubborn second-order PDEs lest any confusion arise regarding which methods are used for which solutions to which PDEs. Definitions 1-6, established by deduction from the features of the Cole-Hopf transformation, are used in this work whenever MCH transformations arise.

Definition 1. 

Let a modified Cole-Hopf transformation (or MCH transformation) be a change of variables which maps a second-order PDE $N\left[u\left(\mathbf{r},t\right),\mathbf{v}\left(\mathbf{r},t\right)·\nabla u\left(\mathbf{r},t\right),\nu {\nabla }^{2}u\left(\mathbf{r},t\right),\right.$$\left.F\left(\mathbf{r},t\right)\right]=0$ to a form $L\left[v\left(\mathbf{r},t\right),\mathbf{v}\left(\mathbf{r},t\right)·\right.$$\left.\nabla v\left(\mathbf{r},t\right),\nu {\nabla }^{2}v\left(\mathbf{r},t\right),F\left(\mathbf{r},t\right)\right]=0$ whose general closed-form solution is known, allowing the solution of the former to be derived from the solution of the latter in the initial and boundary conditions of the former, and the mapping between the two PDEs is determined by the given structure of the former and the prescribed structure of the latter.

Definition 2. 

Let the subject of an MCH transformation be the PDE $N[u\left(\mathbf{r},t\right),\mathbf{v}\left(\mathbf{r},t\right)·\nabla u\left(\mathbf{r},t\right),$$\nu {\nabla }^{2}u\left(\mathbf{r},t\right),$$F\left(\mathbf{r},t\right)]=0$ to be solved via the MCH transformation. In the case of the Cole-Hopf transformation, the MCH subject is the Burgers equation. In the case of this work, the MCH subject is the pseudo-depressurized momentum equations.

Definition 3. 

Let the kernel of an MCH transformation be the PDE $L[v\left(\mathbf{r},t\right),\mathbf{v}\left(\mathbf{r},t\right)·\nabla v\left(\mathbf{r},t\right),$$\nu {\nabla }^{2}v\left(\mathbf{r},t\right),$$F\left(\mathbf{r},t\right)]=0$ to which the MCH transformation seeks to reduce its subject PDE so as to render the subject PDE solvable.

Definition 4. 

Let the residue of an MCH transformation be the expression remaining after the MCH subject terms forming the MCH kernel are canceled by substituting an arbitrary mapping $u=f\left(v\right)$, enabling the solution $v\left(\mathbf{r},t\right)=N^{-1}\left(0\right)$ of the MCH kernel to be expressed as $u\left(\mathbf{r},t\right)=f\left(L^{-1}\left(0\right)\right)$.

Definition 5. 

Let the nomenclature MCH lift refer to the mapping f in Definition 4.

Definition 6. 

Let the reconciliation of an of an MCH transformation be a set of mappings $\{v_0\left(\mathbf{r}\right)=f^{-1}\left(u_0\left(\mathbf{r}\right)\right),v_{\Omega }\left(\mathbf{r},t\right)=f^{-1}\left(u_{\Omega }\left(\mathbf{r},t\right)\right)\}$ between initial data $v_0$ and boundary data $v_{\Omega }\left(\mathbf{r},t\right)$ in the kernel solution $v\left(\mathbf{r},t\right)$ to that $\{u_0\left(\mathbf{r}\right),u_{\Omega }\left(\mathbf{r},t\right)\}$ in $u\left(\mathbf{r},t\right)$ such that the MCH kernel solution $v\left(\mathbf{r},t\right)=L^{-1}\left(0\right)$ might have its initial and boundary conditions expressed in terms of $\{u_0\left(\mathbf{r}\right),u_{\Omega }\left(\mathbf{r},t\right)\}$ since it is used to express $u\left(\mathbf{r},t\right)$ via $u\left(\mathbf{r},t\right)=f\left(L^{-1}\left(0\right)\right)$.

4.3.2. Origin of Turbulence Identified as Quartic Root Multiplicity Collapse in MCH Lift

When the pseudo-depressurized momentum equation are treated as a vector-valued MCH subject and the three MCH kernels assumed to be heat equations of the form ${\Phi }_t\left(\mathbf{r},t\right)=\nu {\nabla }^2\Phi \left(\mathbf{r},t\right)$, the vector Laplacian terms and time derivative terms are absorbed into the kernel, leaving the nonlinear terms equal to the pseudo-depressurized forcing terms composed entirely of known information in the categories of velocity boundary data and the forcing terms $\mathbf{f}\left(\mathbf{r},t\right)$. Writing the nonlinear advective terms as $({\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)·\nabla )u_{\mathrm{rot}}\left(\mathbf{r},t\right)\widehat{ı}+({\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)·\nabla )v_{\mathrm{rot}}\left(\mathbf{r},t\right)\widehat{j}+({\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)·\nabla )w_{\mathrm{rot}}\left(\mathbf{r},t\right)\widehat{k}$ instead of $({\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)·\nabla ){\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)$ allows ${\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)$ to be treated as a fixed-point multiplier accounting for all influences of all parts of ${\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)$ on each of its components, in turn allowing the three MCH residues to be solved for $\{u_{\mathrm{rot}}({\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right),\mathbf{r},t),v_{\mathrm{rot}}({\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right),\mathbf{r},t),w_{\mathrm{rot}}($${\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right),\mathbf{r},t\left)\right\}$ respectively. The aforementioned solution establishes a coupled set of three quadratic polynomials in $\left\{u_{\mathrm{rot}}\left(\mathbf{r},t\right),v_{\mathrm{rot}}\left(\mathbf{r},t\right),w_{\mathrm{rot}}\left(\mathbf{r},t\right)\right\}$ when considered as $\left\{u_{\mathrm{rot}}(u_{\mathrm{rot}}\left(\mathbf{r},t\right),v_{\mathrm{rot}}\left(\mathbf{r},t\right),w_{\mathrm{rot}}\left(\mathbf{r},t\right),\mathbf{r},t),\right.$ $\left.v_{\mathrm{rot}}(u_{\mathrm{rot}}\left(\mathbf{r},t\right),v_{\mathrm{rot}}\left(\mathbf{r},t\right),w_{\mathrm{rot}}\left(\mathbf{r},t\right),\mathbf{r},t),w_{\mathrm{rot}}(u_{\mathrm{rot}}\left(\mathbf{r},t\right),v_{\mathrm{rot}}\left(\mathbf{r},t\right),w_{\mathrm{rot}}\left(\mathbf{r},t\right),\mathbf{r},t)\right\}$. An invocation of Bézout’s theorem implies that each velocity component is formed from the roots of a quartic polynomial equation with real-valued coefficients, suggesting that the number of unique ${\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)$ (and therefore $\mathbf{v}\left(\mathbf{r},t\right)$ since ${\mathbf{v}}_{\mathrm{irr}}\left(\mathbf{r},t\right)$ is not multivalued) per set of total velocity initial and boundary data $\{{\mathbf{v}}_0\left(\mathbf{r}\right),{\mathbf{v}}_{\Omega }\left(\mathbf{r},t\right)\}$ is at least one and at most $2·2^3=16$. Substituting $\alpha \left(\mathbf{r},t\right)=u_{\mathrm{rot}}\left(\mathbf{r},t\right){\left(w_{\mathrm{rot}}\left(\mathbf{r},t\right)\right)}^{-1}$ and $\beta \left(\mathbf{r},t\right)=v_{\mathrm{rot}}\left(\mathbf{r},t\right){\left(w_{\mathrm{rot}}\left(\mathbf{r},t\right)\right)}^{-1}$ into the system for $w_{\mathrm{rot}}\left(\mathbf{r},t\right)\ne 0$ reduces the system to two coupled quadratics $\left\{\alpha \left(\beta \left(\mathbf{r},t\right),\mathbf{r},t\right),\beta \left(\alpha \left(\mathbf{r},t\right),\mathbf{r},t\right)\right\}$ which is far more tractable, as substituting one into the other results in a quartic polynomial in $\beta \left(\mathbf{r},t\right)$, whose roots are well understood to be capable of closed-form expression. Substituting $\alpha \left(\mathbf{r},t\right)$ as a set of quartic roots into the remaining polynomial results in a quadratic in $\beta \left(\mathbf{r},t\right)$. The substitution defining $\{\alpha \left(\mathbf{r},t\right),\beta \left(\mathbf{r},t\right)\}$ coupled with one more sign choice resulting from a square of $w_{\mathrm{rot}}\left(\mathbf{r},t\right)$ arising in the substitution of $\{\alpha \left(\mathbf{r},t\right),\beta \left(\mathbf{r},t\right)\}$ is then used to reconstruct the MCH lifts $\left\{u_{\mathrm{rot}}\left(\mathbf{r},t\right),v_{\mathrm{rot}}\left(\mathbf{r},t\right),w_{\mathrm{rot}}\left(\mathbf{r},t\right)\right\}$ as quartic roots with coefficients formed from the MCH kernels $\{{\Phi }_x\left(\mathbf{r},t\right),{\Phi }_y\left(\mathbf{r},t\right),{\Phi }_z\left(\mathbf{r},t\right)\}$, their first derivatives, and the forcing terms of the pseudo-depressurized momentum equations. The MCH reconciliation begins its derivation at the point marking completion of the MCH residue derivation and the commencement of the MCH lift derivation. Considering the MCH residue as a set of first-order PDEs in $\{{\Phi }_x\left(\mathbf{r},t\right),{\Phi }_y\left(\mathbf{r},t\right),{\Phi }_z\left(\mathbf{r},t\right)\}$ resembling an eikonal equation with $u_{\mathrm{rot}}\left(\mathbf{r},t\right),v_{\mathrm{rot}}\left(\mathbf{r},t\right),w_{\mathrm{rot}}\left(\mathbf{r},t\right)$ as variable coefficients permits solution by the method of characteristics. Evaluating the newly minted $\{{\Phi }_x\left({\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right),\mathbf{r},t\right),{\Phi }_y$$\left.{\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right),\mathbf{r},t\right),{\Phi }_z\left({\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right),\mathbf{r},t\right)\}$ at $t=0$ and on $\mathbf{r}\in \partial \Omega \left(t\right)$ yields all required MCH reconciliations. The edge case $w_{\mathrm{rot}}\left(\mathbf{r},t\right)=0$ not covered by the quartic solution during the MCH lift derivation, corresponding also to two-dimensional flow, is covered by solving the coupled quadratic system in $\{u_{\mathrm{rot}}\left(\mathbf{r},t\right),v_{\mathrm{rot}}\left(\mathbf{r},t\right)\}$ resulting from assuming $w_{\mathrm{rot}}\left(\mathbf{r},t\right)=0$ at the same point at which the MCH reconciliation derivation was begun. Since the case of $w_{\mathrm{rot}}\left(\mathbf{r},t\right)=0$ arises also in situations of two-dimensional flow, the nature of the unevenness in the norm exponents of Ladyzhenskaya’s inequality in three dimensions as opposed to two is hypothesized to involve this quartic-quadratic collapse. An additional edge case is covered wherein $v_{\mathrm{rot}}\left(\mathbf{r},t\right)=w_{\mathrm{rot}}\left(\mathbf{r},t\right)=0$ in which $u_{\mathrm{rot}}\left(\mathbf{r},t\right)$ becomes merely the square root of a Poisson integral divided by the logarithmic gradient of a heat equation solution. The final remaining edge case of an MCH lift in case of unbounded flow, which is a setup with practical relevance only to the Millennium Problem, is addressed in a separate section through Lambert W functions since the derivation chain of the MCH lift forks at the question of bounded or unbounded flow, with the quartic root MCH lift holding only in the presence of boundary conditions. Because quartic equations degenerate into linear, quadratic or cubic equations under certain coefficient symmetry classes while gaining multiplicity of their roots, and all quartic coefficients must vanish in order to produce ${\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)=\mathbf{0}$ when considering irrotational flow as the “state of least turbulence” while not always being wholly laminar, the quartic roots losing multiplicity due to more information being required to encode more complicated aspects of the flow is the most likely mathematical phenomenon to blame for the complexity of turbulence. The algebraic structures of quartic root coefficient symmetries are therefore the most probable means of identifying the exact theoretical edge of turbulence within the solution space of the incompressible Navier-Stokes equations (1) and (2).

4.4. Pressure Recovery by Pressure-Poisson Formulation

The third and final part of the solution to (1) and (2) is the recovery of the pressure field $P\left(\mathbf{r},t\right)$. Taking the divergence of the original momentum equations (1) transforms $\nabla P\left(\mathbf{r},t\right)$ into ${\nabla }^{2}P\left(\mathbf{r},t\right)$, allowing solution for $P\left(\mathbf{r},t\right)$ as a Poisson equation given a Dirichlet boundary condition in $P\left(\mathbf{r},t\right)$, a technique standard in the literature [62]. Neumann or mixed boundary conditions are also admissible if necessary or desired due to the flexibility of the Poisson equation. The final result is a closed-form expression for $P\left(\mathbf{r},t\right)$ when the independent $\mathbf{v}\left(\mathbf{r},t\right)$ is substituted. The two independent, analytically closed expressions $\left\{\mathbf{v}\left(\mathbf{r},t\right),P\left(\mathbf{r},t\right)\right\}$ constitute the general analytical solution of the incompressible Navier-Stokes equations (1) and (2) for a bounded domain.

4.5. Unbounded Flows and the Millennium Problem

In addition to a brief discourse on the implementation of the solution in practical contexts, the remainder of this work takes the form of a deep investigation of the mathematical structure of the velocity and pressure fields. The conditions on $\mathbf{v}\left(\mathbf{r},t\right)$ and $P\left(\mathbf{r},t\right)$ that guarantee their global continuity and smoothness for the purpose of resolving the Clay Millennium Problem amount to global $C^{\infty }$ differentiability and global square-integrability assuming $C^{\infty }$ initial conditions that decay at infinity are evaluated at last with the written form of the general solution in hand. This intersection between solutions from an engineering perspective and from a pure mathematics perspective identifies previously hidden links in the mathematical structure of both. The solution process begins with Claim 1 covering solution for simply connected $\Omega \subset \mathbb{R}$.

Claim 1. 

Let $\mathbf{v}:\left[(\mathbf{r}\in \Omega )\times t\in {\mathbb{R}}_{\ge 0}\right]\to {\mathbb{R}}^3$, $P:\left[(\mathbf{r}\in \Omega )\times t\in {\mathbb{R}}_{\ge 0}\right]\to \mathbb{R}$ and $\mathbf{f}:\left[(\mathbf{r}\in \Omega )\times t\in {\mathbb{R}}_{\ge 0}\right]\to {\mathbb{R}}^3$ on a simply connected $\Omega \subseteq {\mathbb{R}}^3$. Then, there exists a pair $\{\mathbf{v}\left(\mathbf{r},t\right),P\left(\mathbf{r},t\right)\}$ globally satisfying (1) and (2) on Ω for arbitrary initial velocity ${\mathbf{v}}_0\left(\mathbf{r}\right)=\mathbf{v}\left(\mathbf{r},0\right)$ and boundary velocity ${\mathbf{v}}_{\Omega }\left(\mathbf{r},t\right)=\mathbf{v}\left(\mathbf{r}\in \partial \Omega ,t\right)$ under standard regularity assumptions.

5. Solving for the Irrotational Part of the Velocity Field

The derivation of the solution to (1) and (2) is formally begun by seeking a general form of an irrotational and incompressible velocity field (Claim 2.) Though this task is largely known and performed regularly in practical situations, its derivation is nonetheless included for the dual purpose of establishing philosophical and notational context in preparation for the rotational flow derivation that follows, and to display utmost mathematical clarity in disseminating the written solution to (1) and (2).

Claim 2. 

Let $\mathbf{v}:\left[(\mathbf{r}\in \Omega )\times t\in {\mathbb{R}}_{\ge 0}\right]\to {\mathbb{R}}^3$ on a simply connected $\Omega \subseteq {\mathbb{R}}^3$. Then, there exists an irrotational yet solenoidal velocity field ${\mathbf{v}}_{\mathrm{rot}}(\mathbf{r},t)$ satisfying (1) and (2) and constructible from initial ${\mathbf{v}}_0\left(\mathbf{r}\right)=\mathbf{v}(\mathbf{r},0)$ or boundary ${\mathbf{v}}_{\Omega }(\mathbf{r},t)=\mathbf{v}(\mathbf{r}\in \Omega ,t)$ data in $\mathbf{v}(\mathbf{r},t)$ using appropriate irrotational velocity initial and boundary data $\{{\mathbf{v}}_{{\mathrm{irr}}_0}\left(\mathbf{r}\right),{\mathbf{v}}_{{\mathrm{irr}}_{\Omega }}\left(\mathbf{r}\right)\}$ arising from the derivation of ${\mathbf{v}}_{\mathrm{irr}}(\mathbf{r},t)$.

5.1. Invocation of the Helmholtz Decomposition

It is well understood [28] that any $\mathbf{v}\left(\mathbf{r},t\right)$ satisfying (1) and (2) may be separated into an irrotational part ${\mathbf{v}}_{\mathrm{irr}}:(\mathbf{r}\in \Omega \times t\in {\mathbb{R}}_{\ge 0})\to {\mathbb{R}}^3=-\nabla \varphi \left(\mathbf{r},t\right)$ for $\varphi :(\mathbf{r}\in {\mathbb{R}}^3\times t\in {\mathbb{R}}_{\ge 0})\to \mathbb{R}$ so that $\nabla \times {\mathbf{v}}_{\mathrm{irr}}\left(\mathbf{r},t\right)=\mathbf{0}$ everywhere, and a rotational part ${\mathbf{v}}_{\mathrm{rot}}:(\mathbf{r}\in \Omega \times t\in {\mathbb{R}}_{\ge 0})\to {\mathbb{R}}^3|\nabla \times {\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)\ne \mathbf{0}$ via the Helmholtz decomposition (3).

$$\mathbf{v}\left(\mathbf{r},t\right)={\mathbf{v}}_{\mathrm{irr}}\left(\mathbf{r},t\right)+{\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)$$

(3)

5.2. Derivation of the Potential Flow Scalar

It is known [56] that flows which are both incompressible and irrotational are termed potential flows and obey the conservative relationship ${\mathbf{v}}_{\mathrm{irr}}:(\mathbf{r}\in \Omega \times t\in {\mathbb{R}}_{\ge 0})\to {\mathbb{R}}^3=\nabla \varphi \left(\mathbf{r},t\right)$ with the potential flow scalar $\varphi \left(\mathbf{r},t\right)$ also obeying the Laplace equation ${\nabla }^2\varphi \left(\mathbf{r},t\right)=0$. Solving said Laplace equation yields the potential flow scalar (4). Substituting $\varphi \left(\mathbf{r},t\right)$ into $\nabla \varphi \left(\mathbf{r},t\right)$ recovers ${\mathbf{v}}_{\mathrm{irr}}\left(\mathbf{r},t\right)$ (5) in terms of the boundary value ${\varphi }_{\Omega }\left(\mathbf{r},t\right)$ of the potential flow scalar.

$${\nabla }^2\varphi \left(\mathbf{r},t\right)=0\Rightarrow \varphi \left(\mathbf{r},t\right)={\oint }_{\partial \Omega }{\displaystyle \frac{{\varphi }_{\Omega }\left({\mathbf{r}}^{\prime },t\right)}{-2\pi }}{\displaystyle \frac{\partial G_{\mathrm{PD}}\left(\mathbf{r},{\mathbf{r}}^{\prime },t\right)}{\partial {\widehat{N}}_{{\mathbf{r}}^{\prime }}}}d{S}_{{\mathbf{r}}^{\prime }}$$

(4)

$${\mathbf{v}}_{\mathrm{irr}}\left(\mathbf{r},t\right)=-\nabla \varphi \left(\mathbf{r},t\right)=\nabla {\oint }_{\partial \Omega }{\displaystyle \frac{{\varphi }_{\Omega }\left({\mathbf{r}}^{\prime },t\right)}{2\pi }}{\displaystyle \frac{\partial G_{\mathrm{PD}}\left(\mathbf{r},{\mathbf{r}}^{\prime },t\right)}{\partial {\widehat{N}}_{{\mathbf{r}}^{\prime }}}}d{S}_{{\mathbf{r}}^{\prime }}$$

(5)

Inverting the gradient is possible on simply connected $\Omega$. To reconcile Cauchy data on the boundary in ${\mathbf{v}}_{\mathrm{irr}}\left(\mathbf{r},t\right)$ with that of $\varphi \left(\mathbf{r},t\right)$, the inverse gradient of the definition of $\varphi \left(\mathbf{r},t\right)$ is taken and the result $\varphi \left(\mathbf{r},t\right)={\nabla }^{-1}{\mathbf{v}}_{\mathrm{irr}}\left(\mathbf{r},t\right)$ evaluated on the boundary (6). The resultant expression of ${\varphi }_{\Omega }\left(\mathbf{r},t\right)$ in terms of ${\mathbf{v}}_{{\Omega }_{\mathrm{irr}}}\left({\mathbf{r}}^{\prime },t\right)$ is substituted into (5) to recover the irrotational velocity field (7) expressed in its own boundary conditions.

$${\mathbf{v}}_{\mathrm{irr}}\left(\mathbf{r},t\right)=\nabla \varphi \left(\mathbf{r},t\right)\Rightarrow \varphi \left(\mathbf{r},t\right)={\nabla }^{-1}{\mathbf{v}}_{\mathrm{irr}}\left(\mathbf{r},t\right)\Rightarrow {\varphi }_{\Omega }\left(\mathbf{r},t\right)={\nabla }^{-1}{\mathbf{v}}_{{\Omega }_{\mathrm{irr}}}\left(\mathbf{r},t\right)$$

(6)

$${\mathbf{v}}_{\mathrm{irr}}\left(\mathbf{r},t\right)=\nabla {\oint }_{\partial \Omega }{\displaystyle \frac{{\nabla }^{-1}\left({\mathbf{v}}_{{\Omega }_{\mathrm{irr}}}\left({\mathbf{r}}^{\prime },t\right)\right)}{2\pi }}{\displaystyle \frac{\partial G_{\mathrm{PD}}\left(\mathbf{r},{\mathbf{r}}^{\prime },t\right)}{\partial {\widehat{N}}_{{\mathbf{r}}^{\prime }}}}d{S}_{{\mathbf{r}}^{\prime }}$$

(7)

5.3. Inversion of the Helmholtz Decomposition by Projection Operators

Recovering the total velocity field $\mathbf{v}\left(\mathbf{r},t\right)$ from ${\mathbf{v}}_{\mathrm{irr}}\left(\mathbf{r},t\right)$ (7) and ${\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)$ (91) when the latter two are written in terms of irrotational velocity initial and boundary conditions $\left\{{\mathbf{v}}_{0_{\mathrm{irr}}}\left(\mathbf{r}\right),{\mathbf{v}}_{{\Omega }_{\mathrm{irr}}}\left(\mathbf{r},t\right)\right\}$ and rotational velocity initial and boundary conditions $\left\{{\mathbf{v}}_{0_{\mathrm{rot}}}\left(\mathbf{r}\right),{\mathbf{v}}_{{\Omega }_{\mathrm{rot}}}\left(\mathbf{r},t\right)\right\}$ respectively is only possible if $\{{\mathbf{v}}_{\mathrm{irr}}\left(\mathbf{r},t\right),{\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)\}$ are written in terms of total velocity initial and boundary data $\left\{{\mathbf{v}}_0\left(\mathbf{r}\right),{\mathbf{v}}_{\Omega }\left(\mathbf{r},t\right)\right\}$. The relationship between the velocity field and its irrotational and rotational parts as exemplified in the Helmholtz decomposition (3) may be inverted in the sense of extracting the irrotational and rotational parts of the velocity field from the total velocity field. The expressions of the irrotational part (8) and rotational part (9) of the velocity field as functions of the total velocity field are standard components of the derivation of the Helmholtz decomposition itself. The definition of the primed gradient operator ${\nabla }^{\prime }$ is given by (10) while $\widehat{\mathbf{n}}$ is the unit outward normal vector to $\partial \Omega$.

$${\mathbf{v}}_{\mathrm{irr}}\left(\mathbf{r},t\right)=\nabla {\oint }_{\partial \Omega }\widehat{\mathbf{n}}·{\displaystyle \frac{\mathbf{v}\left({\mathbf{r}}^{\prime },t\right)}{4\pi |\mathbf{r}-{\mathbf{r}}^{\prime }|}}d{S}_{{\mathbf{r}}^{\prime }}-\nabla {\int }_{\Omega }{\displaystyle \frac{{\nabla }^{\prime }·\mathbf{v}\left({\mathbf{r}}^{\prime },t\right)}{4\pi |\mathbf{r}-{\mathbf{r}}^{\prime }|}}d{V}_{{\mathbf{r}}^{\prime }}$$

(8)

$${\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)=\nabla \times {\int }_{\Omega }{\displaystyle \frac{{\nabla }^{\prime }\times \mathbf{v}\left({\mathbf{r}}^{\prime },t\right)}{4\pi |\mathbf{r}-{\mathbf{r}}^{\prime }|}}d{V}_{{\mathbf{r}}^{\prime }}-\nabla \times {\oint }_{\partial \Omega }\widehat{\mathbf{n}}\times {\displaystyle \frac{\mathbf{v}\left({\mathbf{r}}^{\prime },t\right)}{4\pi |\mathbf{r}-{\mathbf{r}}^{\prime }|}}d{S}_{{\mathbf{r}}^{\prime }}$$

(9)

$$\nabla \left({\displaystyle \frac{1}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}}\right)=-{\nabla }^{\prime }\left({\displaystyle \frac{1}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}}\right)$$

(10)

Evaluating (8) and (9) at $t=0$ gives expressions for ${\mathbf{v}}_{0_{\mathrm{irr}}}$ (11) and ${\mathbf{v}}_{0_{\mathrm{rot}}}$ (12) respectively.

$${\mathbf{v}}_{0_{\mathrm{irr}}}\left(\mathbf{r}\right)=\nabla {\oint }_{\partial \Omega }\widehat{\mathbf{n}}·{\displaystyle \frac{{\mathbf{v}}_0\left({\mathbf{r}}^{\prime }\right)}{4\pi |\mathbf{r}-{\mathbf{r}}^{\prime }|}}d{S}_{{\mathbf{r}}^{\prime }}-\nabla {\int }_{\Omega }{\displaystyle \frac{{\nabla }^{\prime }·{\mathbf{v}}_0\left({\mathbf{r}}^{\prime }\right)}{4\pi |\mathbf{r}-{\mathbf{r}}^{\prime }|}}d{V}_{{\mathbf{r}}^{\prime }}$$

(11)

$${\mathbf{v}}_{0_{\mathrm{rot}}}\left(\mathbf{r}\right)=\nabla \times {\int }_{\Omega }{\displaystyle \frac{{\nabla }^{\prime }\times {\mathbf{v}}_0\left({\mathbf{r}}^{\prime }\right)}{4\pi |\mathbf{r}-{\mathbf{r}}^{\prime }|}}d{V}_{{\mathbf{r}}^{\prime }}-\nabla \times {\oint }_{\partial \Omega }\widehat{\mathbf{n}}\times {\displaystyle \frac{{\mathbf{v}}_0\left({\mathbf{r}}^{\prime }\right)}{4\pi |\mathbf{r}-{\mathbf{r}}^{\prime }|}}d{S}_{{\mathbf{r}}^{\prime }}$$

(12)

Evaluating (8) and (9) on $\mathbf{r}\in \partial \Omega$ gives expressions for ${\mathbf{v}}_{{\Omega }_{\mathrm{irr}}}$ (13) and ${\mathbf{v}}_{{\Omega }_{\mathrm{rot}}}$ (14) respectively.

$${\mathbf{v}}_{{\Omega }_{\mathrm{irr}}}\left(\mathbf{r},t\right)=\nabla {\oint }_{\partial \Omega }\widehat{\mathbf{n}}·{\displaystyle \frac{{\mathbf{v}}_{\Omega }\left({\mathbf{r}}^{\prime },t\right)}{4\pi |\mathbf{r}-{\mathbf{r}}^{\prime }|}}d{S}_{{\mathbf{r}}^{\prime }}-\nabla {\int }_{\Omega }{\displaystyle \frac{{\nabla }^{\prime }·{\mathbf{v}}_{\Omega }\left({\mathbf{r}}^{\prime },t\right)}{4\pi |\mathbf{r}-{\mathbf{r}}^{\prime }|}}d{V}_{{\mathbf{r}}^{\prime }}$$

(13)

$${\mathbf{v}}_{{\Omega }_{\mathrm{rot}}}\left(\mathbf{r},t\right)=\nabla \times {\int }_{\Omega }{\displaystyle \frac{{\nabla }^{\prime }\times {\mathbf{v}}_{\Omega }\left({\mathbf{r}}^{\prime },t\right)}{4\pi |\mathbf{r}-{\mathbf{r}}^{\prime }|}}d{V}_{{\mathbf{r}}^{\prime }}-\nabla \times {\oint }_{\partial \Omega }\widehat{\mathbf{n}}\times {\displaystyle \frac{{\mathbf{v}}_{\Omega }\left({\mathbf{r}}^{\prime },t\right)}{4\pi |\mathbf{r}-{\mathbf{r}}^{\prime }|}}d{S}_{{\mathbf{r}}^{\prime }}$$

(14)

5.4. Final Statement of Irrotational Velocity Field

Substituting (13) into (7) yields ${\mathbf{v}}_{\mathrm{irr}}\left(\mathbf{r},t\right)$ written in terms of boundary data in $\mathbf{v}\left(\mathbf{r},t\right)$ (15). Because ${\mathbf{v}}_{\mathrm{irr}}\left(\mathbf{r},t\right)$ is irrotational by construction as the gradient of a scalar function $\varphi \left(\mathbf{r},t\right)$ and depends only on arbitrary ${\mathbf{v}}_{\Omega }\left(\mathbf{r},t\right)$ which is not necessarily irrotational but is rather the Dirichlet boundary condition of $\mathbf{v}\left(\mathbf{r},t\right)$ as a whole, Claim 2 stands with its result summarized in Lemma 1.

$${\mathbf{v}}_{\mathrm{irr}}\left(\mathbf{r},t\right)=\nabla {\oint }_{\partial \Omega }\left({\oint }_{\partial \Omega }\widehat{\mathbf{n}}·{\displaystyle \frac{{\mathbf{v}}_{\Omega }\left({\mathbf{r}}^{\prime \prime },t\right)}{8{\pi }^2|{\mathbf{r}}^{\prime }-{\mathbf{r}}^{\prime \prime }|}}d{S}_{{\mathbf{r}}^{\prime \prime }}-{\int }_{\Omega }{\displaystyle \frac{{\nabla }^{\prime }·{\mathbf{v}}_{\Omega }\left({\mathbf{r}}^{\prime \prime },t\right)}{8{\pi }^2|{\mathbf{r}}^{\prime }-{\mathbf{r}}^{\prime \prime }|}}d{V}_{{\mathbf{r}}^{\prime \prime }}\right){\displaystyle \frac{\partial G_{\mathrm{PD}}\left(\mathbf{r},{\mathbf{r}}^{\prime },t\right)}{\partial {\widehat{N}}_{{\mathbf{r}}^{\prime }}}}d{S}_{{\mathbf{r}}^{\prime }}$$

(15)

Lemma 1. 

Let $\Omega \subseteq {\mathbb{R}}^3$ be simply connected and $\mathbf{v}:\left[(\mathbf{r}\in \Omega )\times t\in {\mathbb{R}}_{\ge 0}\right]\to {\mathbb{R}}^3$. Then, there exists an irrotational velocity field ${\mathbf{v}}_{\mathrm{irr}}\left(\mathbf{r},t\right)$ (15) satisfying (1) and (2) and constructible from Dirichlet boundary data ${\mathbf{v}}_{\Omega }(\mathbf{r},t)=\mathbf{v}(\mathbf{r}\in \Omega ,t)$ of $\mathbf{v}(\mathbf{r},t)$ by inverting the Helmholtz decomposition and evaluating the result on the boundary $\partial \Omega$ (13).

6. Solving for the Rotational Part of the Velocity Field

Almost all of the unknown structure of the Navier-Stokes equations lies in the rotational part of the velocity field. The nonlinear interactions between the advective terms of the momentum equations (1) are nearly restricted out of existence under a vanishing velocity curl, and if left to themselves under no simplifications at all, are necessary preconditions to turbulence. Therefore, the bulk of the attention given, both in this work (Claim 3) and in its predecessors seeking the closure of the incompressible Navier-Stokes equations, is and has been to the rotational part of the velocity field.

Claim 3. 

Let $\mathbf{v}:\left[(\mathbf{r}\in \Omega )\times t\in {\mathbb{R}}_{\ge 0}\right]\to {\mathbb{R}}^3$ and $\mathbf{f}:\left[(\mathbf{r}\in \Omega )\times t\in {\mathbb{R}}_{\ge 0}\right]\to {\mathbb{R}}^3$ on $\Omega \subseteq {\mathbb{R}}^3$. Then, there exists a rotational yet solenoidal velocity field ${\mathbf{v}}_{\mathrm{rot}}(\mathbf{r},t)$ satisfying (1) and (2) and constructible from initial data ${\mathbf{v}}_0\left(\mathbf{r}\right)=\mathbf{v}(\mathbf{r},0)$ or boundary data ${\mathbf{v}}_{\Omega }(\mathbf{r},t)=\mathbf{v}(\mathbf{r}\in \Omega ,t)$ of $\mathbf{v}(\mathbf{r},t)$ via (12) and (14) respectively in conjunction with appropriate rotational velocity initial and boundary data $\{{\mathbf{v}}_{{\mathrm{rot}}_0}\left(\mathbf{r}\right),{\mathbf{v}}_{{\mathrm{rot}}_{\Omega }}\left(\mathbf{r}\right)\}$ arising from the derivation of ${\mathbf{v}}_{\mathrm{rot}}(\mathbf{r},t)$, from which the total $\mathbf{v}(\mathbf{r},t)$ may be recovered by substitution into the Helmholtz decomposition (3) along with ${\mathbf{v}}_{\mathrm{irr}}\left(\mathbf{r},t\right)$ (Theorem 1.)

6.1. Second-Curl Reduction of Momentum Equations

Reduction of the well-known vector field identity $\nabla \times \left(\nabla \times \mathbf{v}\right)=\nabla \left(\nabla ·\mathbf{v}\right)-{\nabla }^2\mathbf{v}$ to $\nabla \times \left(\nabla \times \mathbf{v}\right)=-{\nabla }^2\mathbf{v}$ is not only enabled, but mandated, by incompressibility which requires that $\nabla ·\mathbf{v}=0$. Applying said reduced identity to (1) results in the fourth-order PDE (16) sans pressure per the additional vector field identities $\left\{\nabla \times \left(\nabla P\right)=\mathbf{0}\wedge \nabla \times \mathbf{0}=\mathbf{0}\right\}\forall P:{\mathbb{R}}^3\to \mathbb{R}$.

$$\nabla \times \left(\nabla \times \left({\displaystyle \frac{\partial \mathbf{v}\left(\mathbf{r},t\right)}{\partial t}}\right)\right)+\nabla \times \left(\nabla \times \left(\left(\mathbf{v}\left(\mathbf{r},t\right)·\nabla \right)\mathbf{v}\left(\mathbf{r},t\right)\right)\right)=\nabla \times \left(\nabla \times \mathbf{f}\left(\mathbf{r},t\right)\right)+$$

$$\nabla \times \left(\nabla \times \left(\nu {\nabla }^2\mathbf{v}\left(\mathbf{r},t\right)\right)\right)$$

$$\Rightarrow {\displaystyle \frac{\partial {\nabla }^2{\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)}{\partial t}}+{\nabla }^2\left(\left({\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)·\nabla \right){\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)\right)={\nabla }^2\mathbf{f}\left(\mathbf{r},t\right)+\nu {\nabla }^4{\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)$$

(16)

The distributive property of addition borne by the vector Laplacian as a linear differential operator permits the vector Laplacians of (16) to be factored, yielding a vector Poisson equation (17) whose unknown is itself the terms of the momentum equations sans the external force and pressure terms. The vector Laplacian of the external force terms is treated as the forcing term of this Poisson equation to fully encapsulate its time dependency in the interior of $\Omega$ as a given quantity. Treating ${\nabla }^2\mathbf{f}\left(\mathbf{r},t\right)$ as part of the Poisson unknown would cause the solution of (17) to account only for the boundary values of ${\nabla }^2\mathbf{f}\left(\mathbf{r},t\right)$ as opposed to all of its values.

$${\nabla }^2\left({\displaystyle \frac{\partial {\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)}{\partial t}}+\left({\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)·\nabla \right){\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)-\nu {\nabla }^2{\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)\right)={\nabla }^2\mathbf{f}\left(\mathbf{r},t\right)$$

(17)

The result of solving the vector Poisson equation (17) using the well-studied Green’s function convolution method [63] is a three-dimensional forced Burgers equation (18) whose forcing term involves its own initial and boundary conditions in addition to the Navier-Stokes external force term. This recursive, highly nonlinear relationship between the velocity and its own Cauchy data could shed light not only on the self-similarity observed in turbulence [57], but also on what additional deep structure must be added to a three-dimensional Burgers equation in order to render it an effective Navier-Stokes equation. Since this form of the Burgers equation does not yet possess a standard name in the literature, both the discussion of the matter in the mathematical community and the nomenclatural clarity demanded of a Navier-Stokes solution is enabled by christening (18) the pseudo-depressurized form of the momentum equations in light of their absorption of pressure effects into their velocity Cauchy data.

$${\displaystyle \frac{\partial {\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)}{\partial t}}+\left({\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)·\nabla \right){\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)-\nu {\nabla }^2{\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)={\int }_{\Omega }{\displaystyle \frac{{\nabla }^2\mathbf{f}\left({\mathbf{r}}^{\prime },t\right)}{4\pi }}{G}_{\mathrm{PD}}(\mathbf{r},{\mathbf{r}}^{\prime })d{V}_{{\mathbf{r}}^{\prime }}+{\displaystyle \frac{1}{2\pi }}{\oint }_{\partial \Omega }$$

$$\left(\nu {\nabla }^2{\mathbf{v}}_{{\mathrm{rot}}_{\Omega }}\left({\mathbf{r}}^{\prime },t\right)-{\displaystyle \frac{\partial {\mathbf{v}}_{{\mathrm{rot}}_{\Omega }}\left({\mathbf{r}}^{\prime },t\right)}{\partial t}}-\left({\mathbf{v}}_{{\mathrm{rot}}_{\Omega }}\left({\mathbf{r}}^{\prime },t\right)·\nabla \right){\mathbf{v}}_{{\mathrm{rot}}_{\Omega }}\left({\mathbf{r}}^{\prime },t\right)\right){\displaystyle \frac{\partial G_{\mathrm{PD}}(\mathbf{r},{\mathbf{r}}^{\prime })}{\partial {\widehat{N}}_{{\mathbf{r}}^{\prime }}}}d{S}_{{\mathbf{r}}^{\prime }}$$

(18)

6.2. Solving the Pseudo-Depressurized Momentum Equations by MCH

Close observation of the patterns in the advective terms $\left(\mathbf{v}\left(\mathbf{r},t\right)·\nabla \right)\mathbf{v}\left(\mathbf{r},t\right)$ of the pseudo-depressurized momentum equations (18) allows said terms to be rewritten as $\left(\mathbf{v}\left(\mathbf{r},t\right)·\nabla \right)v_i\left(\mathbf{r},t\right)$ for the ith momentum equation (19) and (20) and (21) without loss of generality since the notation $\mathbf{v}\left(\mathbf{r},t\right)$ is nothing more than a shorthand for the $v_i\left(\mathbf{r},t\right)$ as vector field components. The expected outcome of this step is to find closed form $v_i\left(\mathbf{r},t\right)$ taking $\mathbf{v}\left(\mathbf{r},t\right)$ as a parameter such that the components $\{v_1,v_2,v_3\}=\{u,v,w\}$ of $\mathbf{v}\left(\mathbf{r},t\right)$ might be expanded, enabling at least an algebraic solution for explicit $\{u,v,w\}$ independent of one another in the next step. This treatment of the advective terms has some literary precedent in analytical [58] and numerical [59] studies alike which term said methods “fixed-point” solution methods. The functions in this subsection have their independent variables omitted for clarity unless so-called “dummy variables” are needed for integration.

$${\displaystyle \frac{\partial u_{\mathrm{rot}}\left(\mathbf{r},t\right)}{\partial t}}+\left({\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)·\nabla \right)u_{\mathrm{rot}}\left(\mathbf{r},t\right)-\nu {\nabla }^2{u}_{\mathrm{rot}}\left(\mathbf{r},t\right)={\displaystyle \frac{1}{4\pi }}{\int }_{\Omega }{\nabla }^2{f}_1\left({\mathbf{r}}^{\prime },t\right)G_{\mathrm{PD}}(\mathbf{r},{\mathbf{r}}^{\prime })d{V}_{{\mathbf{r}}^{\prime }}+{\displaystyle \frac{1}{2\pi }}{\oint }_{\partial \Omega }$$

$$\left(\nu {\nabla }^2{u}_{{\mathrm{rot}}_{\Omega }}\left({\mathbf{r}}^{\prime },t\right)-{\displaystyle \frac{\partial u_{{\mathrm{rot}}_{\Omega }}\left({\mathbf{r}}^{\prime },t\right)}{\partial t}}-\left({\mathbf{v}}_{{\mathrm{rot}}_{\Omega }}\left({\mathbf{r}}^{\prime },t\right)·\nabla \right)u_{{\mathrm{rot}}_{\Omega }}\left({\mathbf{r}}^{\prime },t\right)\right){\displaystyle \frac{\partial G_{\mathrm{PD}}(\mathbf{r},{\mathbf{r}}^{\prime })}{\partial {\widehat{N}}_{{\mathbf{r}}^{\prime }}}}d{S}_{{\mathbf{r}}^{\prime }}{\widehat{x}}_1$$

(19)

$${\displaystyle \frac{\partial v_{\mathrm{rot}}\left(\mathbf{r},t\right)}{\partial t}}+\left({\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)·\nabla \right)v_{\mathrm{rot}}\left(\mathbf{r},t\right)-\nu {\nabla }^2{v}_{\mathrm{rot}}\left(\mathbf{r},t\right)={\displaystyle \frac{1}{4\pi }}{\int }_{\Omega }{\nabla }^2{f}_2\left({\mathbf{r}}^{\prime },t\right)G_{\mathrm{PD}}(\mathbf{r},{\mathbf{r}}^{\prime })d{V}_{{\mathbf{r}}^{\prime }}+{\displaystyle \frac{1}{2\pi }}{\oint }_{\partial \Omega }$$

$$\left(\nu {\nabla }^2{v}_{{\mathrm{rot}}_{\Omega }}\left({\mathbf{r}}^{\prime },t\right)-{\displaystyle \frac{\partial v_{{\mathrm{rot}}_{\Omega }}\left({\mathbf{r}}^{\prime },t\right)}{\partial t}}-\left({\mathbf{v}}_{{\mathrm{rot}}_{\Omega }}\left({\mathbf{r}}^{\prime },t\right)·\nabla \right)v_{{\mathrm{rot}}_{\Omega }}\left({\mathbf{r}}^{\prime },t\right)\right){\displaystyle \frac{\partial G_{\mathrm{PD}}(\mathbf{r},{\mathbf{r}}^{\prime })}{\partial {\widehat{N}}_{{\mathbf{r}}^{\prime }}}}d{S}_{{\mathbf{r}}^{\prime }}{\widehat{x}}_2$$

(20)

$${\displaystyle \frac{\partial w_{\mathrm{rot}}\left(\mathbf{r},t\right)}{\partial t}}+\left({\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)·\nabla \right)w_{\mathrm{rot}}\left(\mathbf{r},t\right)-\nu {\nabla }^2{w}_{\mathrm{rot}}\left(\mathbf{r},t\right)={\displaystyle \frac{1}{4\pi }}{\int }_{\Omega }{\nabla }^2{f}_3\left({\mathbf{r}}^{\prime },t\right)G_{\mathrm{PD}}(\mathbf{r},{\mathbf{r}}^{\prime })d{V}_{{\mathbf{r}}^{\prime }}+{\displaystyle \frac{1}{2\pi }}{\oint }_{\partial \Omega }$$

$$\left(\nu {\nabla }^2{w}_{{\mathrm{rot}}_{\Omega }}\left({\mathbf{r}}^{\prime },t\right)-{\displaystyle \frac{\partial w_{{\mathrm{rot}}_{\Omega }}\left({\mathbf{r}}^{\prime },t\right)}{\partial t}}-\left({\mathbf{v}}_{{\mathrm{rot}}_{\Omega }}\left({\mathbf{r}}^{\prime },t\right)·\nabla \right)w_{{\mathrm{rot}}_{\Omega }}\left({\mathbf{r}}^{\prime },t\right)\right){\displaystyle \frac{\partial G_{\mathrm{PD}}(\mathbf{r},{\mathbf{r}}^{\prime })}{\partial {\widehat{N}}_{{\mathbf{r}}^{\prime }}}}d{S}_{{\mathbf{r}}^{\prime }}{\widehat{x}}_3$$

(21)

The pseudo-depressurized momentum equations with augmented advective term notation (19) and (20) and (21) are now staged for an implicit fixed-point analytical solution by an MCH method per Definitions 1-5. Commencing with a substitution of $u_{\mathrm{rot}}\left(\mathbf{r},t\right)={\Theta }_x^{\prime }\left({\Phi }_x\left(\mathbf{r},t\right)\right)$, $v_{\mathrm{rot}}\left(\mathbf{r},t\right)={\Theta }_y^{\prime }\left({\Phi }_y\left(\mathbf{r},t\right)\right)$, and $w_{\mathrm{rot}}\left(\mathbf{r},t\right)={\Theta }_z^{\prime }\left({\Phi }_z\left(\mathbf{r},t\right)\right)$ into (19), (20), and (21) respectively while leaving the ${\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)$ in $\left({\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)·\nabla \right)v_{i_{\mathrm{rot}}}\left(\mathbf{r},t\right)$ as a fixed point generates a system of ODEs in ${\Theta }_i\left({\Phi }_i\left(\mathbf{r},t\right)\right)$ and of PDEs in ${\Phi }_i\left(\mathbf{r},t\right)$. Terms involving initial and boundary conditions, which are counted as given information, are collectively treated as known forcing terms akin to the treatment of $\mathbf{f}\left(\mathbf{r},t\right)$.

$${\Theta }_1^{\prime }\left({\Phi }_1\right){\displaystyle \frac{\partial {\Phi }_1}{\partial t}}+{\Theta }_1^{\prime }\left({\Phi }_1\right)\left(\left(\mathbf{v}·\nabla \right){\Phi }_1\right)-{\Theta }_1^{\prime }\left({\Phi }_1\right)\nu {\nabla }^2{\Phi }_1-\nu {\Theta }_1^{\prime \prime }\left({\Phi }_1\right){\left|\left|\nabla {\Phi }_1\right|\right|}^2={\int }_{\Omega }{\displaystyle \frac{{\nabla }^2{f}_1\left({\mathbf{r}}^{\prime },t\right)}{4\pi }}{G}_{\mathrm{PD}}(\mathbf{r},{\mathbf{r}}^{\prime })d{V}_{{\mathbf{r}}^{\prime }}$$

$$+{\oint }_{\partial \Omega }\left(\nu {\nabla }^2{u}_{{\mathrm{rot}}_{\Omega }}\left({\mathbf{r}}^{\prime },t\right)-{\displaystyle \frac{\partial u_{{\mathrm{rot}}_{\Omega }}\left({\mathbf{r}}^{\prime },t\right)}{\partial t}}-\left({\mathbf{v}}_{{\mathrm{rot}}_{\Omega }}\left({\mathbf{r}}^{\prime },t\right)·\nabla \right)u_{{\mathrm{rot}}_{\Omega }}\left({\mathbf{r}}^{\prime },t\right)\right){\displaystyle \frac{\frac{\partial G_{\mathrm{PD}}(\mathbf{r},{\mathbf{r}}^{\prime })}{\partial {\widehat{N}}_{{\mathbf{r}}^{\prime }}}}{2\pi }}d{S}_{{\mathbf{r}}^{\prime }}{\widehat{x}}_1$$

(22)

$${\Theta }_2^{\prime }\left({\Phi }_2\right){\displaystyle \frac{\partial {\Phi }_2}{\partial t}}+{\Theta }_2^{\prime }\left({\Phi }_2\right)\left(\left(\mathbf{v}·\nabla \right){\Phi }_2\right)-{\Theta }_2^{\prime }\left({\Phi }_2\right)\nu {\nabla }^2{\Phi }_2-\nu {\Theta }_2^{\prime \prime }\left({\Phi }_2\right){\left|\left|\nabla {\Phi }_2\right|\right|}^2={\int }_{\Omega }{\displaystyle \frac{{\nabla }^2{f}_2\left({\mathbf{r}}^{\prime },t\right)}{4\pi }}{G}_{\mathrm{PD}}(\mathbf{r},{\mathbf{r}}^{\prime })d{V}_{{\mathbf{r}}^{\prime }}$$

$$+{\oint }_{\partial \Omega }\left(\nu {\nabla }^2{v}_{{\mathrm{rot}}_{\Omega }}\left({\mathbf{r}}^{\prime },t\right)-{\displaystyle \frac{\partial v_{{\mathrm{rot}}_{\Omega }}\left({\mathbf{r}}^{\prime },t\right)}{\partial t}}-\left({\mathbf{v}}_{{\mathrm{rot}}_{\Omega }}\left({\mathbf{r}}^{\prime },t\right)·\nabla \right)v_{{\mathrm{rot}}_{\Omega }}\left({\mathbf{r}}^{\prime },t\right)\right){\displaystyle \frac{\frac{\partial G_{\mathrm{PD}}(\mathbf{r},{\mathbf{r}}^{\prime })}{\partial {\widehat{N}}_{{\mathbf{r}}^{\prime }}}}{2\pi }}d{S}_{{\mathbf{r}}^{\prime }}{\widehat{x}}_2$$

(23)

$${\Theta }_3^{\prime }\left({\Phi }_3\right){\displaystyle \frac{\partial {\Phi }_3}{\partial t}}+{\Theta }_3^{\prime }\left({\Phi }_3\right)\left(\left(\mathbf{v}·\nabla \right){\Phi }_3\right)-{\Theta }_3^{\prime }\left({\Phi }_3\right)\nu {\nabla }^2{\Phi }_3-\nu {\Theta }_3^{\prime \prime }\left({\Phi }_3\right){\left|\left|\nabla {\Phi }_3\right|\right|}^2={\int }_{\Omega }{\displaystyle \frac{{\nabla }^2{f}_3\left({\mathbf{r}}^{\prime },t\right)}{4\pi }}{G}_{\mathrm{PD}}(\mathbf{r},{\mathbf{r}}^{\prime })d{V}_{{\mathbf{r}}^{\prime }}$$

$$+{\oint }_{\partial \Omega }\left(\nu {\nabla }^2{w}_{{\mathrm{rot}}_{\Omega }}\left({\mathbf{r}}^{\prime },t\right)-{\displaystyle \frac{\partial w_{{\mathrm{rot}}_{\Omega }}\left({\mathbf{r}}^{\prime },t\right)}{\partial t}}-\left({\mathbf{v}}_{{\mathrm{rot}}_{\Omega }}\left({\mathbf{r}}^{\prime },t\right)·\nabla \right)w_{{\mathrm{rot}}_{\Omega }}\left({\mathbf{r}}^{\prime },t\right)\right){\displaystyle \frac{\frac{\partial G_{\mathrm{PD}}(\mathbf{r},{\mathbf{r}}^{\prime })}{\partial {\widehat{N}}_{{\mathbf{r}}^{\prime }}}}{2\pi }}d{S}_{{\mathbf{r}}^{\prime }}{\widehat{x}}_3$$

(24)

6.3. Defining the MCH Kernel as a Vector of Heat Equation Solutions

If the MCH kernels ${\Phi }_k\left(\mathbf{r},t\right)$ are constrained such that the componentwise differences between their time derivative and Laplacian terms in (22), (23), (24) vanish, then the ${\Phi }_k\left(\mathbf{r},t\right)$ are defined as solutions of heat equations (25) while the remaining terms (26)–(28) of (22)–(24) respectively form the MCH residues, establishing relationships between (25) and the rotational velocity field by ${\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)={\Theta }_1\left({\Phi }_1\left(\mathbf{r},t\right)\right){\widehat{x}}_1+{\Theta }_2\left({\Phi }_2\left(\mathbf{r},t\right)\right){\widehat{x}}_2+{\Theta }_3\left({\Phi }_3\left(\mathbf{r},t\right)\right){\widehat{x}}_3$, providing an equal number of equations and unknowns necessary to close the system (25) and (26) and (27) and (28) with the known methods of solving heat equations beginning the cascade. The functions $B_1\left(\mathbf{r},t\right)$, $B_2\left(\mathbf{r},t\right)$, and $B_3\left(\mathbf{r},t\right)$ are used from this point forward to abbreviate the RHS of (26) and (27) and (28) which consist entirely of Poisson domain and boundary integrals.

$$\left[\begin{array}{c}{\Theta }_1^{\prime }\left({\Phi }_1\left(\mathbf{r},t\right)\right){\displaystyle \frac{\partial {\Phi }_1\left(\mathbf{r},t\right)}{\partial t}}=\nu {\Theta }_1^{\prime }\left({\Phi }_1\left(\mathbf{r},t\right)\right){\nabla }^2{\Phi }_1\left(\mathbf{r},t\right){\widehat{x}}_1\\ {\Theta }_2^{\prime }\left({\Phi }_2\left(\mathbf{r},t\right)\right){\displaystyle \frac{\partial {\Phi }_2\left(\mathbf{r},t\right)}{\partial t}}=\nu {\Theta }_2^{\prime }\left({\Phi }_2\left(\mathbf{r},t\right)\right){\nabla }^2{\Phi }_2\left(\mathbf{r},t\right){\widehat{x}}_2\\ {\Theta }_3^{\prime }\left({\Phi }_3\left(\mathbf{r},t\right)\right){\displaystyle \frac{\partial {\Phi }_3\left(\mathbf{r},t\right)}{\partial t}}=\nu {\Theta }_3^{\prime }\left({\Phi }_3\left(\mathbf{r},t\right)\right){\nabla }^2{\Phi }_3\left(\mathbf{r},t\right){\widehat{x}}_3\end{array}\right]\Rightarrow \left[\begin{array}{c}{\displaystyle \frac{\partial {\Phi }_1\left(\mathbf{r},t\right)}{\partial t}}=\nu {\nabla }^2{\Phi }_1\left(\mathbf{r},t\right){\widehat{x}}_1\\ {\displaystyle \frac{\partial {\Phi }_2\left(\mathbf{r},t\right)}{\partial t}}=\nu {\nabla }^2{\Phi }_2\left(\mathbf{r},t\right){\widehat{x}}_2\\ {\displaystyle \frac{\partial {\Phi }_3\left(\mathbf{r},t\right)}{\partial t}}=\nu {\nabla }^2{\Phi }_3\left(\mathbf{r},t\right){\widehat{x}}_3\end{array}\right]$$

(25)

$${\Theta }_1^{\prime }\left({\Phi }_1\left(\mathbf{r},t\right)\right)\left(\left({\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)·\nabla \right){\Phi }_1\left(\mathbf{r},t\right)\right)-\nu {\Theta }_1^{\prime \prime }\left({\Phi }_1\left(\mathbf{r},t\right)\right){\left|\left|\nabla {\Phi }_1\left(\mathbf{r},t\right)\right|\right|}^2={\int }_{\Omega }{\displaystyle \frac{{\nabla }^2{f}_1\left({\mathbf{r}}^{\prime },t\right)d{V}_{{\mathbf{r}}^{\prime }}}{4\pi {\left(G_{\mathrm{PD}}(\mathbf{r},{\mathbf{r}}^{\prime })\right)}^{-1}}}$$

$$+{\oint }_{\partial \Omega }{\displaystyle \frac{\nu {\nabla }^2{u}_{{\mathrm{rot}}_{\Omega }}\left({\mathbf{r}}^{\prime },t\right)-\frac{\partial u_{{\mathrm{rot}}_{\Omega }}\left({\mathbf{r}}^{\prime },t\right)}{\partial t}-\left({\mathbf{v}}_{{\mathrm{rot}}_{\Omega }}\left({\mathbf{r}}^{\prime },t\right)·\nabla \right)u_{{\mathrm{rot}}_{\Omega }}\left({\mathbf{r}}^{\prime },t\right)}{2\pi {\left(\frac{\partial G_{\mathrm{PD}}(\mathbf{r},{\mathbf{r}}^{\prime })}{\partial {\widehat{N}}_{{\mathbf{r}}^{\prime }}}\right)}^{-1}}}d{S}_{{\mathbf{r}}^{\prime }}=B_1\left(\mathbf{r},t\right){\widehat{x}}_1$$

(26)

$${\Theta }_2^{\prime }\left({\Phi }_2\left(\mathbf{r},t\right)\right)\left(\left({\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)·\nabla \right){\Phi }_2\left(\mathbf{r},t\right)\right)-\nu {\Theta }_2^{\prime \prime }\left({\Phi }_2\left(\mathbf{r},t\right)\right){\left|\left|\nabla {\Phi }_2\left(\mathbf{r},t\right)\right|\right|}^2={\int }_{\Omega }{\displaystyle \frac{{\nabla }^2{f}_2\left({\mathbf{r}}^{\prime },t\right)d{V}_{{\mathbf{r}}^{\prime }}}{4\pi {\left(G_{\mathrm{PD}}(\mathbf{r},{\mathbf{r}}^{\prime })\right)}^{-1}}}$$

$$+{\oint }_{\partial \Omega }{\displaystyle \frac{\nu {\nabla }^2{v}_{{\mathrm{rot}}_{\Omega }}\left({\mathbf{r}}^{\prime },t\right)-\frac{\partial v_{{\mathrm{rot}}_{\Omega }}\left({\mathbf{r}}^{\prime },t\right)}{\partial t}-\left({\mathbf{v}}_{{\mathrm{rot}}_{\Omega }}\left({\mathbf{r}}^{\prime },t\right)·\nabla \right)v_{{\mathrm{rot}}_{\Omega }}\left({\mathbf{r}}^{\prime },t\right)}{2\pi {\left(\frac{\partial G_{\mathrm{PD}}(\mathbf{r},{\mathbf{r}}^{\prime })}{\partial {\widehat{N}}_{{\mathbf{r}}^{\prime }}}\right)}^{-1}}}d{S}_{{\mathbf{r}}^{\prime }}=B_2\left(\mathbf{r},t\right){\widehat{x}}_2$$

(27)

$${\Theta }_3^{\prime }\left({\Phi }_3\left(\mathbf{r},t\right)\right)\left(\left({\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)·\nabla \right){\Phi }_3\left(\mathbf{r},t\right)\right)-\nu {\Theta }_3^{\prime \prime }\left({\Phi }_3\left(\mathbf{r},t\right)\right){\left|\left|\nabla {\Phi }_3\left(\mathbf{r},t\right)\right|\right|}^2={\int }_{\Omega }{\displaystyle \frac{{\nabla }^2{f}_3\left({\mathbf{r}}^{\prime },t\right)d{V}_{{\mathbf{r}}^{\prime }}}{4\pi {\left(G_{\mathrm{PD}}(\mathbf{r},{\mathbf{r}}^{\prime })\right)}^{-1}}}$$

$$+{\oint }_{\partial \Omega }{\displaystyle \frac{\nu {\nabla }^2{w}_{{\mathrm{rot}}_{\Omega }}\left({\mathbf{r}}^{\prime },t\right)-\frac{\partial w_{{\mathrm{rot}}_{\Omega }}\left({\mathbf{r}}^{\prime },t\right)}{\partial t}-\left({\mathbf{v}}_{{\mathrm{rot}}_{\Omega }}\left({\mathbf{r}}^{\prime },t\right)·\nabla \right)w_{{\mathrm{rot}}_{\Omega }}\left({\mathbf{r}}^{\prime },t\right)}{2\pi {\left(\frac{\partial G_{\mathrm{PD}}(\mathbf{r},{\mathbf{r}}^{\prime })}{\partial {\widehat{N}}_{{\mathbf{r}}^{\prime }}}\right)}^{-1}}}d{S}_{{\mathbf{r}}^{\prime }}=B_3\left(\mathbf{r},t\right){\widehat{x}}_3$$

(28)

Solving the three heat equations (25) is easily achievable with standard Green’s function convolution methods [62]. $G_{\mathrm{HD}}\left(\mathbf{r},{\mathbf{r}}^{\prime },t\right)$ is the Green’s function associated with the heat equation under Dirichlet boundary conditions. The initial conditions ${\Phi }_{k_0}\left(\mathbf{r}\right)$ and Dirichlet boundary conditions ${\Phi }_{k_{\Omega \left(t\right)}}$ in the heat equation solutions ${\Phi }_k\left(\mathbf{r},t\right)$ for $k\in {\mathbb{N}}_{[1,3]}$ must eventually be written in terms of their velocity field counterparts $\{u,v,w\}$ to ensure analytical closure via MCH reconciliation.

$${\Phi }_k\left(\mathbf{r},t\right)={\int }_{\Omega }{\Phi }_{k_0}\left({\mathbf{r}}^{\prime }\right)G_{\mathrm{HD}}\left(\mathbf{r},{\mathbf{r}}^{\prime },t\right)d{V}_{{\mathbf{r}}^{\prime }}+\nu {\int }_0^t{\oint }_{\partial \Omega }{\Phi }_{k_{\Omega }}\left({\mathbf{r}}^{\prime },\tau \right){\displaystyle \frac{\partial G_{\mathrm{HD}}\left(\mathbf{r},{\mathbf{r}}^{\prime },t-\tau \right)}{\partial {\widehat{N}}_{{\mathbf{r}}^{\prime }}}}d{S}_{{\mathbf{r}}^{\prime }}d\tau$$

(29)

6.4. Deriving the MCH Lift as Roots of a Triquadratic Algebraic System

Treating the MCH residues (26) and (27) and (28) as a set of three second-order linear ODEs in ${\Theta }_k\left({\Phi }_k\left(\mathbf{r},t\right)\right)$ respectively allows for their solution (30) and (31) and (32). Because $u_{\mathrm{rot}}\left(\mathbf{r},t\right)={\Theta }_1^{\prime }\left({\Phi }_1\left(\mathbf{r},t\right)\right)$, $v_{\mathrm{rot}}={\Theta }_2^{\prime }\left({\Phi }_2\left(\mathbf{r},t\right)\right)$, and $w_{\mathrm{rot}}={\Theta }_3^{\prime }\left({\Phi }_3\left(\mathbf{r},t\right)\right)$, a relationship between ${\Phi }_k\left(\mathbf{r},t\right)$, $u_{\mathrm{rot}}\left(\mathbf{r},t\right)$, $v_{\mathrm{rot}}\left(\mathbf{r},t\right)$, and $w_{\mathrm{rot}}\left(\mathbf{r},t\right)$ is derived (33) which serves as an implicit MCH lift. The constants of integration $C_{ij}|(i,j)\in \{0,1,2\}\times \{0,1\}$ are free to choose after the precedent of the original Cole-Hopf transformation since their value does not affect the ability of (30)–(32) to satisfy (26)–(28) respectively and thus serve as an MCH lift, and it is practical to choose $C_{ij}$ such that the resultant expressions are as simple as possible while continuing to satisfy (26)–(28) respectively. To minimize nonlinearity in the MCH lifts while preserving the nonlinearity of ${\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)$, it is chosen that $C_{ij}=0$ for bounded flows. For the unique case of completely unbounded flow suffusing all of space with fluid, which must be addressed for the Millennium Problem, a simultaneous vanishing of $C_{ij}$ and $B_k\left(\mathbf{r},t\right)$ would require that ${\mathbf{v}}_{\mathrm{rot}}=\mathbf{0}$, which is both insufficiently general and physically unrealistic for any case. For the unbounded Millennium Problem case, $C_{i0}=1\wedge C_{i1}=0$ is the simplest choice of $C_{ij}$, which causes (30)–(32) to reduce to (33). Since unboundedness is the lone exception to the rule that almost all practical flow problems involve a boundary in at least some respect, the case of bounded flow is assumed from this point forth and unto the final statement of the total velocity field, while the unbounded case will be dealt with in its own section specifically addressing the Millennium Problem.

$$u_{\mathrm{rot}}\left(\mathbf{r},t\right)={\Theta }_1\left({\Phi }_1\left(\mathbf{r},t\right)\right)={\displaystyle \frac{\nu C_{00}{\left|\left|\nabla {\Phi }_1\left(\mathbf{r},t\right)\right|\right|}^2}{\left({\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)·\nabla \right){\Phi }_1\left(\mathbf{r},t\right)}}exp\left({\displaystyle \frac{\left(\left({\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)·\nabla \right){\Phi }_1\left(\mathbf{r},t\right)\right){\Phi }_1\left(\mathbf{r},t\right)}{\nu {\left|\left|\nabla {\Phi }_1\left(\mathbf{r},t\right)\right|\right|}^2}}\right)-$$

$${\displaystyle \frac{B_1\left(\mathbf{r},t\right){\Phi }_1\left(\mathbf{r},t\right)}{\left(\mathbf{v}·\nabla \right){\Phi }_1\left(\mathbf{r},t\right)}}+C_{01}{\widehat{x}}_1$$

(30)

$$v_{\mathrm{rot}}\left(\mathbf{r},t\right)={\Theta }_2\left({\Phi }_2\left(\mathbf{r},t\right)\right)={\displaystyle \frac{\nu C_{10}{\left|\left|\nabla {\Phi }_2\left(\mathbf{r},t\right)\right|\right|}^2}{\left({\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)·\nabla \right){\Phi }_2\left(\mathbf{r},t\right)}}exp\left({\displaystyle \frac{\left(\left({\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)·\nabla \right){\Phi }_2\left(\mathbf{r},t\right)\right){\Phi }_2\left(\mathbf{r},t\right)}{\nu {\left|\left|\nabla {\Phi }_2\left(\mathbf{r},t\right)\right|\right|}^2}}\right)-$$

$${\displaystyle \frac{B_2\left(\mathbf{r},t\right){\Phi }_2\left(\mathbf{r},t\right)}{\left(\mathbf{v}·\nabla \right){\Phi }_2\left(\mathbf{r},t\right)}}+C_{11}{\widehat{x}}_2$$

(31)

$$w_{\mathrm{rot}}\left(\mathbf{r},t\right)={\Theta }_3\left({\Phi }_3\left(\mathbf{r},t\right)\right)={\displaystyle \frac{\nu C_{20}{\left|\left|\nabla {\Phi }_3\left(\mathbf{r},t\right)\right|\right|}^2}{\left({\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)·\nabla \right){\Phi }_3\left(\mathbf{r},t\right)}}exp\left({\displaystyle \frac{\left(\left({\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)·\nabla \right){\Phi }_3\left(\mathbf{r},t\right)\right){\Phi }_3\left(\mathbf{r},t\right)}{\nu {\left|\left|\nabla {\Phi }_3\left(\mathbf{r},t\right)\right|\right|}^2}}\right)-$$

$${\displaystyle \frac{B_3\left(\mathbf{r},t\right){\Phi }_3\left(\mathbf{r},t\right)}{\left(\mathbf{v}·\nabla \right){\Phi }_3\left(\mathbf{r},t\right)}}+C_{21}{\widehat{x}}_3$$

(32)

$$u_{\mathrm{rot}}\left(\mathbf{r},t\right)=-{\displaystyle \frac{B_1\left(\mathbf{r},t\right){\Phi }_1\left(\mathbf{r},t\right)}{\left({\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)·\nabla \right){\Phi }_1\left(\mathbf{r},t\right)}}{\widehat{x}}_1,v_{\mathrm{rot}}\left(\mathbf{r},t\right)=-{\displaystyle \frac{B_2\left(\mathbf{r},t\right){\Phi }_2\left(\mathbf{r},t\right)}{\left({\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)·\nabla \right){\Phi }_2\left(\mathbf{r},t\right)}}{\widehat{x}}_2,$$

$$w_{\mathrm{rot}}\left(\mathbf{r},t\right)=-{\displaystyle \frac{B_3\left(\mathbf{r},t\right){\Phi }_3\left(\mathbf{r},t\right)}{\left({\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)·\nabla \right){\Phi }_3\left(\mathbf{r},t\right)}}{\widehat{x}}_3$$

(33)

Expanding the notation ${\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)$ as $u_{\mathrm{rot}}\left(\mathbf{r},t\right){\widehat{x}}_1+v_{\mathrm{rot}}\left(\mathbf{r},t\right){\widehat{x}}_2+w_{\mathrm{rot}}\left(\mathbf{r},t\right){\widehat{x}}_3$ allows (33) to be rewritten (34)–(36) solely in terms of $u_{\mathrm{rot}}\left(\mathbf{r},t\right)$, $v_{\mathrm{rot}}\left(\mathbf{r},t\right)$, and $w_{\mathrm{rot}}\left(\mathbf{r},t\right)$. This rendition of the MCH residue constitutes a nonlinear system of three algebraic equations in the three velocity components and is therefore neither overdetermined nor underdetermined. The notation ${\nabla }_m{\Phi }_k$, denoting the mth component of $\nabla {\Phi }_k$, is used to retain generality in the chosen coordinate system.

$$u_{\mathrm{rot}}\left(\mathbf{r},t\right)={\displaystyle \frac{-B_1\left(\mathbf{r},t\right){\Phi }_1\left(\mathbf{r},t\right)}{u_{\mathrm{rot}}\left(\mathbf{r},t\right){\nabla }_1{\Phi }_1\left(\mathbf{r},t\right)+v_{\mathrm{rot}}\left(\mathbf{r},t\right){\nabla }_2{\Phi }_1\left(\mathbf{r},t\right)+w_{\mathrm{rot}}\left(\mathbf{r},t\right){\nabla }_3{\Phi }_1\left(\mathbf{r},t\right)}}{\widehat{x}}_1$$

(34)

$$v_{\mathrm{rot}}\left(\mathbf{r},t\right)={\displaystyle \frac{-B_2\left(\mathbf{r},t\right){\Phi }_2\left(\mathbf{r},t\right)}{u_{\mathrm{rot}}\left(\mathbf{r},t\right){\nabla }_1{\Phi }_2\left(\mathbf{r},t\right)+v_{\mathrm{rot}}\left(\mathbf{r},t\right){\nabla }_2{\Phi }_2\left(\mathbf{r},t\right)+w_{\mathrm{rot}}\left(\mathbf{r},t\right){\nabla }_3{\Phi }_2\left(\mathbf{r},t\right)}}{\widehat{x}}_2$$

(35)

$$w_{\mathrm{rot}}\left(\mathbf{r},t\right)={\displaystyle \frac{-B_3\left(\mathbf{r},t\right){\Phi }_3\left(\mathbf{r},t\right)}{u_{\mathrm{rot}}\left(\mathbf{r},t\right){\nabla }_1{\Phi }_3\left(\mathbf{r},t\right)+v_{\mathrm{rot}}\left(\mathbf{r},t\right){\nabla }_2{\Phi }_3\left(\mathbf{r},t\right)+w_{\mathrm{rot}}\left(\mathbf{r},t\right){\nabla }_3{\Phi }_3\left(\mathbf{r},t\right)}}{\widehat{x}}_3$$

(36)

Multiplying (34)–(36) by the denominators of their RHS recasts the system (34)–(36) as a quadratically nonlinear algebraic system (37)–(39). The task of finding the MCH lifts, namely expressions of $u_{\mathrm{rot}}\left(\mathbf{r},t\right)$, $v_{\mathrm{rot}}\left(\mathbf{r},t\right)$, and $w_{\mathrm{rot}}\left(\mathbf{r},t\right)$ in terms of ${\Phi }_k\left(\mathbf{r},t\right)|k\in {\mathbb{N}}_{[1,3]}$, from the MCH residues has been reduced to solving (37)–(39).

$$u_{\mathrm{rot}}^2\left(\mathbf{r},t\right){\nabla }_1{\Phi }_1\left(\mathbf{r},t\right)+u_{\mathrm{rot}}\left(\mathbf{r},t\right)v_{\mathrm{rot}}\left(\mathbf{r},t\right){\nabla }_2{\Phi }_1\left(\mathbf{r},t\right)+u_{\mathrm{rot}}\left(\mathbf{r},t\right)w_{\mathrm{rot}}\left(\mathbf{r},t\right){\nabla }_3{\Phi }_1\left(\mathbf{r},t\right)=$$

$$-B_1\left(\mathbf{r},t\right){\Phi }_1\left(\mathbf{r},t\right){\widehat{x}}_1$$

(37)

$$u_{\mathrm{rot}}\left(\mathbf{r},t\right)v_{\mathrm{rot}}\left(\mathbf{r},t\right){\nabla }_1{\Phi }_2\left(\mathbf{r},t\right)+v_{\mathrm{rot}}^2\left(\mathbf{r},t\right){\nabla }_2{\Phi }_2\left(\mathbf{r},t\right)+v_{\mathrm{rot}}\left(\mathbf{r},t\right)w_{\mathrm{rot}}\left(\mathbf{r},t\right){\nabla }_3{\Phi }_2\left(\mathbf{r},t\right)=$$

$$-B_2\left(\mathbf{r},t\right){\Phi }_2\left(\mathbf{r},t\right){\widehat{x}}_2$$

(38)

$$u_{\mathrm{rot}}\left(\mathbf{r},t\right)w_{\mathrm{rot}}\left(\mathbf{r},t\right){\nabla }_1{\Phi }_3\left(\mathbf{r},t\right)+v_{\mathrm{rot}}\left(\mathbf{r},t\right)w_{\mathrm{rot}}\left(\mathbf{r},t\right){\nabla }_2{\Phi }_3\left(\mathbf{r},t\right)+w_{\mathrm{rot}}^2\left(\mathbf{r},t\right){\nabla }_3{\Phi }_3\left(\mathbf{r},t\right)=$$

$$-B_3\left(\mathbf{r},t\right){\Phi }_3\left(\mathbf{r},t\right){\widehat{x}}_3$$

(39)

Substituting (41) into the algebraic system of equations (37)–(39) governing the MCH lifts reduces the latter to the simpler biquadratic system (41) and (42) with one constraint equation (43) while the coefficients of (37)–(39) are abbreviated (44) for clarity.

$$\alpha \left(\mathbf{r},t\right)={\displaystyle \frac{u_{\mathrm{rot}}\left(\mathbf{r},t\right)}{w_{\mathrm{rot}}\left(\mathbf{r},t\right)}},\beta \left(\mathbf{r},t\right)={\displaystyle \frac{v_{\mathrm{rot}}\left(\mathbf{r},t\right)}{w_{\mathrm{rot}}\left(\mathbf{r},t\right)}}$$

(40)

$$w_{\mathrm{rot}}^2\left(\mathbf{r},t\right)\left(a_1{\alpha }^2\left(\mathbf{r},t\right)+a_2\alpha \left(\mathbf{r},t\right)\beta \left(\mathbf{r},t\right)+a_3\alpha \left(\mathbf{r},t\right)\right)=a_4$$

(41)

$$w_{\mathrm{rot}}^2\left(\mathbf{r},t\right)\left(a_5\alpha \left(\mathbf{r},t\right)\beta \left(\mathbf{r},t\right)+a_6{\beta }^2\left(\mathbf{r},t\right)+a_7\beta \left(\mathbf{r},t\right)\right)=a_8$$

(42)

$$w_{\mathrm{rot}}^2\left(\mathbf{r},t\right)\left(a_9\alpha \left(\mathbf{r},t\right)+a_{10}\beta \left(\mathbf{r},t\right)+a_{11}\right)=a_{12}$$

(43)

$$\left[\begin{array}{cccc}{a}_1& a_2& a_3& a_4\\ a_5& a_6& a_7& a_8\\ a_9& a_{10}& a_{11}& a_{12}\end{array}\right]=\left[\begin{array}{cccc}{\nabla }_1{\Phi }_1\left(\mathbf{r},t\right)& {\nabla }_2{\Phi }_1\left(\mathbf{r},t\right)& {\nabla }_3{\Phi }_1\left(\mathbf{r},t\right)& -B_1\left(\mathbf{r},t\right){\Phi }_1\left(\mathbf{r},t\right)\\ {\nabla }_1{\Phi }_2\left(\mathbf{r},t\right)& {\nabla }_2{\Phi }_2\left(\mathbf{r},t\right)& {\nabla }_3{\Phi }_2\left(\mathbf{r},t\right)& -B_2\left(\mathbf{r},t\right){\Phi }_2\left(\mathbf{r},t\right)\\ {\nabla }_1{\Phi }_3\left(\mathbf{r},t\right)& {\nabla }_2{\Phi }_3\left(\mathbf{r},t\right)& {\nabla }_3{\Phi }_3\left(\mathbf{r},t\right)& -B_3\left(\mathbf{r},t\right){\Phi }_3\left(\mathbf{r},t\right)\end{array}\right]$$

(44)

Substituting $w_{\mathrm{rot}}^2=a_{12}{\left(a_9\alpha \left(\mathbf{r},t\right)+a_{10}\beta \left(\mathbf{r},t\right)+a_{11}\right)}^{-1}$ into (41) and (42) results in two coupled quadratic polynomials (45) and (46) in $\alpha \left(\mathbf{r},t\right)$ and $\beta \left(\mathbf{r},t\right)$. According to Bézout’s theorem [60], the maximum number of total unique roots of a coupled system in the class of (45) and (46) is four, implying that there exist up to two unique values of $\alpha \left(\mathbf{r},t\right)$ and $\beta \left(\mathbf{r},t\right)$ for each set of initial and boundary conditions in ${\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)$.

$$a_{12}{a}_1{\alpha }^2\left(\mathbf{r},t\right)+a_{12}{a}_2\alpha \left(\mathbf{r},t\right)\beta \left(\mathbf{r},t\right)+\left(a_{12}{a}_3-a_4{a}_9\right)\alpha \left(\mathbf{r},t\right)-a_4{a}_{10}\beta \left(\mathbf{r},t\right)=a_4{a}_{11}$$

(45)

$$a_{12}{a}_5\alpha \left(\mathbf{r},t\right)\beta \left(\mathbf{r},t\right)+a_{12}{a}_6{\beta }^2\left(\mathbf{r},t\right)+\left(a_{12}{a}_7-a_8{a}_{10}\right)\beta \left(\mathbf{r},t\right)-a_8{a}_9\alpha \left(\mathbf{r},t\right)=a_8{a}_{11}$$

(46)

Solving (45) for $\alpha \left(\mathbf{r},t\right)$ (47) expresses $\alpha \left(\mathbf{r},t\right)$ in terms of $\beta \left(\mathbf{r},t\right)$.

$$\alpha \left(\mathbf{r},t\right)={\displaystyle \frac{a_4{a}_9-a_{12}{a}_3-a_{12}{a}_2\beta \left(\mathbf{r},t\right)}{a_{12}{a}_1}}\pm$$

$${\displaystyle \frac{\sqrt{{\left(a_{12}{a}_3-a_4{a}_9+a_{12}{a}_2\beta \left(\mathbf{r},t\right)\right)}^2+4a_{12}{a}_1\left(a_4{a}_{11}+a_4{a}_{10}\beta \left(\mathbf{r},t\right)\right)}}{a_{12}{a}_1}}$$

(47)

Substituting (47) into (46) results in the quartic equation (48) whose abbreviated coefficients (49) contain layered abbreviations (50) due to their immense typographical size.

$$A\beta {\left(\mathbf{r},t\right)}^4+B\beta {\left(\mathbf{r},t\right)}^3+C\beta {\left(\mathbf{r},t\right)}^2+D\beta \left(\mathbf{r},t\right)+E=0$$

(48)

$$\left\{\begin{array}{c}A={\gamma }_1{\delta }_1^2-{\gamma }_3{\delta }_1{\delta }_3\hfill \\ B=2{\gamma }_1{\delta }_1{\delta }_2+{\gamma }_4{\delta }_3^2+2{\gamma }_1{\delta }_1{\delta }_4-{\gamma }_3{\delta }_2{\delta }_3-{\gamma }_2{\delta }_1{\delta }_3-{\gamma }_3{\delta }_3{\delta }_4\hfill \\ C=2{\gamma }_1{\delta }_1{\delta }_5+{\gamma }_1{\delta }_2^2+2{\gamma }_1{\delta }_2{\delta }_4+{\gamma }_1{\delta }_4^2+{\gamma }_5{\delta }_3^2-{\gamma }_2{\delta }_2{\delta }_3-{\gamma }_2{\delta }_3{\delta }_4-{\gamma }_3{\delta }_3{\delta }_5\hfill \\ D=2{\gamma }_1{\delta }_2{\delta }_5+2{\gamma }_1{\delta }_4{\delta }_5-{\gamma }_2{\delta }_3{\delta }_5\hfill \\ E={\gamma }_1{\delta }_5^2\hfill \end{array}\right.$$

(49)

$$\left\{\begin{array}{c}{\gamma }_1=a_{12}{a}_1\hfill \\ {\gamma }_2=a_{12}{a}_2\hfill \\ {\gamma }_3=a_{12}{a}_3-a_4{a}_9\hfill \\ {\gamma }_4=-a_4{a}_{10}\hfill \\ {\gamma }_5=-a_4{a}_{11}\hfill \end{array}\right.\left\{\begin{array}{c}{\delta }_1=a_{12}{a}_6\hfill \\ {\delta }_2=a_{12}{a}_5\hfill \\ {\delta }_3=a_{12}{a}_7-a_8{a}_{10}\hfill \\ {\delta }_4=-a_8{a}_9\hfill \\ {\delta }_5=-a_8{a}_{11}\hfill \end{array}\right.$$

(50)

The four roots of the general quartic equation (48) are given by () where the subscripts $\kappa$ and $\sigma$ indicate causal connections in sign choice. Since the coefficients of the quartic (48) are all real-valued due to the heat and Poisson equation solutions which form all coefficients inheriting the real-valuedness of their Cauchy data, which depends ultimately on that of rotational velocity field Cauchy data, the quartic (48) must possess at least one real root. Therefore, regardless of the coefficient values (49) and (50), there always exists at least one real $\alpha \left(\mathbf{r},t\right)$, with the differences between any other real $\alpha \left(\mathbf{r},t\right)$ consisting only in the sign choices exemplified in (51). In situations where vanishing coefficients cause (51) to degenerate into polynomials of lower degree, closed-form solutions for $\alpha \left(\mathbf{r},t\right)$ are known by the cubic and quadratic formulae, with the only difference between degenerate cases and the full quartic case being the number of unique real $\alpha \left(\mathbf{r},t\right)$ per set of initial and boundary conditions in ${\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)$.

$$\beta \left(\mathbf{r},t\right)={\pm }_{\sigma }\sqrt{{\displaystyle \frac{a}{3}}+{\displaystyle \frac{\zeta }{2}}-{\displaystyle \frac{\frac{a^2}{6}+2c}{12\zeta }}}{\pm }_{\kappa }\sqrt{-\left({\displaystyle \frac{5a}{6}}+{\displaystyle \frac{\zeta }{2}}+{\displaystyle \frac{\frac{a^2}{12}+c}{6\zeta }}{\pm }_{\sigma }{\displaystyle \frac{b}{\sqrt{\frac{4a}{3}+2\zeta +\frac{\frac{a^2}{6}+2c}{3\zeta }}}}\right)}-{\displaystyle \frac{B}{4A}}|$$

$$\{\zeta =\sqrt[3]{\sqrt{{\displaystyle \frac{1}{4}}{\left({\displaystyle \frac{ac}{3}}-{\displaystyle \frac{a^3}{108}}-{\displaystyle \frac{b^2}{8}}\right)}^2+{\displaystyle \frac{1}{27}}{\left(-{\displaystyle \frac{a^2}{12}}-c\right)}^3}-{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{ac}{3}}-{\displaystyle \frac{a^3}{108}}-{\displaystyle \frac{b^2}{8}}\right)},.$$

$$.a={\displaystyle \frac{C}{A}}-{\displaystyle \frac{3B^2}{8A^2}},b={\displaystyle \frac{B^3}{8A^3}}-{\displaystyle \frac{BC}{2A^2}}+{\displaystyle \frac{D}{A}},c={\displaystyle \frac{C{B}^2}{16A^3}}+{\displaystyle \frac{E}{A}}-{\displaystyle \frac{BD}{4A^2}}-{\displaystyle \frac{3B^4}{256A^4}}\}$$

(51)

The full expression of $\alpha \left(\mathbf{r},t\right)$ in terms of (44) is recovered by substituting (51) into (47). Substituting the definitions (40) of $\alpha \left(\mathbf{r},t\right)$ and $\beta \left(\mathbf{r},t\right)$ into one another (52) and the result into (43) enables direct recovery of the MCH lift (53) of the ${\widehat{x}}_1$-component $u_{\mathrm{rot}}\left(\mathbf{r},t\right)$ of the velocity field from $\alpha \left(\mathbf{r},t\right)$, $\beta \left(\mathbf{r},t\right)$, and (44).

$$\alpha \left(\mathbf{r},t\right)={\displaystyle \frac{u_{\mathrm{rot}}\left(\mathbf{r},t\right)}{w_{\mathrm{rot}}\left(\mathbf{r},t\right)}}\wedge \beta \left(\mathbf{r},t\right)={\displaystyle \frac{v_{\mathrm{rot}}\left(\mathbf{r},t\right)}{w_{\mathrm{rot}}\left(\mathbf{r},t\right)}}\Rightarrow {\displaystyle \frac{u_{\mathrm{rot}}\left(\mathbf{r},t\right)}{\alpha \left(\mathbf{r},t\right)}}=w_{\mathrm{rot}}\left(\mathbf{r},t\right)={\displaystyle \frac{v_{\mathrm{rot}}\left(\mathbf{r},t\right)}{\beta \left(\mathbf{r},t\right)}}$$

(52)

$$u_{\mathrm{rot}}\left(\mathbf{r},t\right)=\pm \sqrt{{\displaystyle \frac{{\nabla }_x{\Phi }_x\left(\mathbf{r},t\right)}{{\nabla }_y{\Phi }_x\left(\mathbf{r},t\right)+{\nabla }_z{\Phi }_x\left(\mathbf{r},t\right)\alpha \left(\mathbf{r},t\right)-B_x\left(\mathbf{r},t\right)\beta \left(\mathbf{r},t\right)}}}{\widehat{x}}_1$$

(53)

The MCH lifts for $v_{\mathrm{rot}}\left(\mathbf{r},t\right)$ (54) and $w_{\mathrm{rot}}\left(\mathbf{r},t\right)$ (55) are obtained from the expression (53) for $u_{\mathrm{rot}}\left(\mathbf{r},t\right)$ by substituting (53) into each of (40) and solving for $v_{\mathrm{rot}}\left(\mathbf{r},t\right)$ and $w_{\mathrm{rot}}\left(\mathbf{r},t\right)$ respectively.

$$v_{\mathrm{rot}}\left(\mathbf{r},t\right)={\displaystyle \frac{\alpha \left(\mathbf{r},t\right)u_{\mathrm{rot}}\left(\mathbf{r},t\right)}{\beta \left(\mathbf{r},t\right)}}=\pm \sqrt{{\displaystyle \frac{{\alpha }^2\left(\mathbf{r},t\right){\nabla }_x{\Phi }_x\left(\mathbf{r},t\right)/{\beta }^2\left(\mathbf{r},t\right)}{{\nabla }_y{\Phi }_x\left(\mathbf{r},t\right)+{\nabla }_z{\Phi }_x\left(\mathbf{r},t\right)\alpha \left(\mathbf{r},t\right)-B_x\left(\mathbf{r},t\right)\beta \left(\mathbf{r},t\right)}}}{\widehat{x}}_2$$

(54)

$$w_{\mathrm{rot}}\left(\mathbf{r},t\right)={\displaystyle \frac{\alpha \left(\mathbf{r},t\right)u_{\mathrm{rot}}\left(\mathbf{r},t\right)}{{\beta }^2\left(\mathbf{r},t\right)}}=\pm \sqrt{{\displaystyle \frac{{\alpha }^2\left(\mathbf{r},t\right){\nabla }_x{\Phi }_x\left(\mathbf{r},t\right)/{\beta }^4\left(\mathbf{r},t\right)}{{\nabla }_y{\Phi }_x\left(\mathbf{r},t\right)+{\nabla }_z{\Phi }_x\left(\mathbf{r},t\right)\alpha \left(\mathbf{r},t\right)-B_x\left(\mathbf{r},t\right)\beta \left(\mathbf{r},t\right)}}}{\widehat{x}}_3$$

(55)

The finding (53)–(55) of all three components of ${\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)$ establishes the MCH lift (56) of ${\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)$ in vector notation.

$${\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)=\pm {\displaystyle \frac{1{\widehat{x}}_1+\frac{\alpha \left(\mathbf{r},t\right)}{\beta \left(\mathbf{r},t\right)}{\widehat{x}}_2+\frac{\alpha \left(\mathbf{r},t\right)}{{\beta }^2\left(\mathbf{r},t\right)}{\widehat{x}}_3}{{\left(\frac{{\nabla }_x{\Phi }_x\left(\mathbf{r},t\right)}{{\nabla }_y{\Phi }_x\left(\mathbf{r},t\right)+{\nabla }_z{\Phi }_x\left(\mathbf{r},t\right)\alpha \left(\mathbf{r},t\right)-B_x\left(\mathbf{r},t\right)\beta \left(\mathbf{r},t\right)}\right)}^{-1/2}}}$$

(56)

6.4.1. MCH Lift Edge Case: One Irrotational or Zero Velocity Component

For $w_{\mathrm{rot}}\left(\mathbf{r},t\right)\ne 0$ and $a_{12}\ne 0$, the solution (53)–(55) of the triquadratic MCH lift system is complete and bijective for all $(u_{\mathrm{rot}}\left(\mathbf{r},t\right),v_{\mathrm{rot}}\left(\mathbf{r},t\right),w_{\mathrm{rot}}\left(\mathbf{r},t\right))\in {\mathbb{R}}^3$ except $w_{\mathrm{rot}}\left(\mathbf{r},t\right)=0$. Accounting for the branch $w_{\mathrm{rot}}\left(\mathbf{r},t\right)=0$ is therefore required to maintain the internal self-consistency of the MCH lift solution (53)–(55) and therefore for the consistency of the solution to (1) and (2) at large. Assuming $w_{\mathrm{rot}}\left(\mathbf{r},t\right)=0$ from the point (41)–(43) in the proof causes (41)–(43) to degenerate into a biquadratic MCH kernel system (57) and (58).

$$u_{\mathrm{rot}}^2\left(\mathbf{r},t\right){\nabla }_1{\Phi }_1\left(\mathbf{r},t\right)+u_{\mathrm{rot}}\left(\mathbf{r},t\right)v_{\mathrm{rot}}\left(\mathbf{r},t\right){\nabla }_2{\Phi }_1\left(\mathbf{r},t\right)=-B_1\left(\mathbf{r},t\right){\Phi }_1\left(\mathbf{r},t\right){\widehat{x}}_1$$

(57)

$$u_{\mathrm{rot}}\left(\mathbf{r},t\right)v_{\mathrm{rot}}\left(\mathbf{r},t\right){\nabla }_1{\Phi }_2\left(\mathbf{r},t\right)+v_{\mathrm{rot}}^2\left(\mathbf{r},t\right){\nabla }_2{\Phi }_2\left(\mathbf{r},t\right)=-B_2\left(\mathbf{r},t\right){\Phi }_2\left(\mathbf{r},t\right){\widehat{x}}_2$$

(58)

Transcribing the ${\Phi }_k\left(\mathbf{r},t\right)$ notation in (57) and (58) to the $a_k|k\in {\mathbb{N}}_{[1,12]}$ notation (44) increases the legibility of the polynomial algebraic manipulations necessary to recover ${\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)$ from ${\Phi }_k\left(\mathbf{r},t\right)$.

$$u_{\mathrm{rot}}^2\left(\mathbf{r},t\right)a_1+u_{\mathrm{rot}}\left(\mathbf{r},t\right)v_{\mathrm{rot}}\left(\mathbf{r},t\right)a_2=a_4,u_{\mathrm{rot}}\left(\mathbf{r},t\right)v_{\mathrm{rot}}\left(\mathbf{r},t\right)a_5+v_{\mathrm{rot}}{\left(\mathbf{r},t\right)}^2{a}_6=a_8$$

(59)

Using $r\left(\mathbf{r},t\right)$ to denote the magnitude of $u\left(\mathbf{r},t\right)$ in units of $v\left(\mathbf{r},t\right)$ permits both $u\left(\mathbf{r},t\right)$ and $v\left(\mathbf{r},t\right)$ to be found in terms of $r\left(\mathbf{r},t\right)$ (60) and (61). Substituting (60) into (61) yields a quadratic equation in $r\left(\mathbf{r},t\right)$ with solution (62).

$$u_{\mathrm{rot}}\left(\mathbf{r},t\right)=r\left(\mathbf{r},t\right)v_{\mathrm{rot}}\left(\mathbf{r},t\right)$$

(60)

$$v_{\mathrm{rot}}\left(\mathbf{r},t\right)=\sqrt{{\displaystyle \frac{a_8}{a_{5}r\left(\mathbf{r},t\right)+a_6}}}$$

(61)

$$r\left(\mathbf{r},t\right)={\displaystyle \frac{a_4{a}_5-a_2{a}_8\pm \sqrt{{\left(a_4{a}_5-a_2{a}_8\right)}^2+4a_1{a}_4{a}_6{a}_8}}{a_1{a}_8}}$$

(62)

Substituting (62) into both (60) and (61) yields $u_{\mathrm{rot}}\left(\mathbf{r},t\right)$ (63) and $v_{\mathrm{rot}}\left(\mathbf{r},t\right)$ (64) in terms of ${\Phi }_k\left(\mathbf{r},t\right)$, forming the bidimensionally rotational velocity field as an MCH lift, ready for its MCH reconciliation between the heat equations and the velocity field in the next section. Since the choice of sign appears within expressions appearing multiple times in the same expression, the subscript $ϵ$ is included to emphasize the causal commonality of sign choices across the aforementioned expressions. There are therefore at most two ${\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)=u_{\mathrm{rot}}\left(\mathbf{r},t\right){\widehat{x}}_1+v_{\mathrm{rot}}\left(\mathbf{r},t\right){\widehat{x}}_2$ which may be formed from the same set of coefficients (44).

$$u_{\mathrm{rot}}\left(\mathbf{r},t\right)={\displaystyle \frac{a_4{a}_5-a_2{a}_8{\pm }_{ϵ}\sqrt{{\left(a_4{a}_5-a_2{a}_8\right)}^2+4a_1{a}_4{a}_6{a}_8}}{a_1{a}_8{\left(\frac{a_8}{a_5\left(\frac{a_4{a}_5-a_2{a}_8{\pm }_{ϵ}\sqrt{{\left(a_4{a}_5-a_2{a}_8\right)}^2+4a_1{a}_4{a}_6{a}_8}}{a_1{a}_8}\right)+a_6}\right)}^{-1/2}}}{\widehat{x}}_1$$

(63)

$$v_{\mathrm{rot}}\left(\mathbf{r},t\right)=\sqrt{{\displaystyle \frac{a_8}{a_5\left(\frac{a_4{a}_5-a_2{a}_8{\pm }_{ϵ}\sqrt{{\left(a_4{a}_5-a_2{a}_8\right)}^2+4a_1{a}_4{a}_6{a}_8}}{a_1{a}_8}\right)+a_6}}}{\widehat{x}}_2$$

(64)

The situation (63) and (64) arises also as a direct consequence of assuming two-dimensional flow. The markedly simpler polynomial unwrapping of the MCH lift structure as two coupled quadratic equations (57) and (58) as opposed to three coupled quartics (37)–(39), whose explicit root expressions require only one component-rational substitution (62) as opposed to two (40), decreases the possible multivaluedness of $v_{\mathrm{rot}}\left(\mathbf{r},t\right)$ and consequently $v\left(\mathbf{r},t\right)$. It is well-known that there exists a generally lesser degree of chaoticity in two-dimensional flows as opposed to three-dimensional flows due to the impossibility of vortex stretching, suggesting the possibility that sufficient vortex stretching to produce nontrivial $w_{\mathrm{rot}}\left(\mathbf{r},t\right)$ arises when the quartic MCH lifts lose multiplicity of at least two of their roots. At any rate, the Euclidean norm exponent lateralization in the three-dimensional case of the Ladyzhenskaya inequality and the equal exponents in its two-dimensional case may find its explanation by this quartic-quadratic collapse encoding dimensional degeneracy.

6.4.2. MCH Lift Edge Case: Two Irrotational or Zero Velocity Components

The additional edge case of (37)–(39) where $v_{\mathrm{rot}}\left(\mathbf{r},t\right)=w_{\mathrm{rot}}\left(\mathbf{r},t\right)=0$ causes the former to reduce to one quadratic equation in $u_{\mathrm{rot}}\left(\mathbf{r},t\right)$ with the straightforward square root solution (65). Similar to (63) and (64), there are at most two ${\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)=u_{\mathrm{rot}}\left(\mathbf{r},t\right){\widehat{x}}_1$ which may be formed from the same set of coefficients (44). This edge case bears one limitation, namely that $u_{\mathrm{rot}}\left(\mathbf{r},t\right)\in \mathbb{C}\forall B_1\left(\mathbf{r},t\right)>0$, implying that (44) is most reliable in situations of boundary advection, although such situations almost always involve all three rotational velocity components being nonzero. The most consistent interpretation of (44) as regards the question of global existence of a real-valued velocity field on $\Omega \subset {\mathbb{R}}^3$ is the assignment of a value to ${\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)$ when $v_{\mathrm{rot}}\left(\mathbf{r},t\right)$ and $u_{\mathrm{rot}}\left(\mathbf{r},t\right)$ simultaneously change signs. In such a situation on an unbounded domain, $B_1\left(\mathbf{r},t\right)=0$ and therefore ${\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)=\mathbf{0}$. Since ${\mathbf{v}}_{\mathrm{irr}}\left(\mathbf{r},t\right)$ only involves boundary values of $\mathbf{v}\left(\mathbf{r},t\right)$, specifying $\Omega ={\mathbb{R}}^3\wedge v_{\mathrm{rot}}\left(\mathbf{r},t\right)=w_{\mathrm{rot}}\left(\mathbf{r},t\right)=0$ ultimately causes $\mathbf{v}\left(\mathbf{r},t\right)=\mathbf{0}$.

$$u_{\mathrm{rot}}^2\left(\mathbf{r},t\right){\nabla }_1{\Phi }_1\left(\mathbf{r},t\right)=-B_1\left(\mathbf{r},t\right){\Phi }_1\left(\mathbf{r},t\right){\widehat{x}}_1\Rightarrow u_{\mathrm{rot}}\left(\mathbf{r},t\right)=\pm \sqrt{{\displaystyle \frac{-B_1\left(\mathbf{r},t\right){\Phi }_1\left(\mathbf{r},t\right)}{{\nabla }_1{\Phi }_1\left(\mathbf{r},t\right)}}}{\widehat{x}}_1$$

(65)

6.5. Deriving the MCH Reconciliation by Inversion of the MCH Residue

Rewriting (37)–(39) in vector form (66)–(68) emphasizes the fact that the MCH residue system, when framed as an MCH reconciliation whose unknowns are ${\Phi }_k\left(\mathbf{r},t\right)|k\in \{1,2,3\}$ rather than as an MCH lift whose unknowns are the components of ${\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)$, is decoupled in ${\Phi }_k\left(\mathbf{r},t\right)$ and therefore solvable in parallel.

$$\left({\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)·\nabla \right){\Phi }_1\left(\mathbf{r},t\right)=-{\displaystyle \frac{B_1\left(\mathbf{r},t\right)}{u_{\mathrm{rot}}\left(\mathbf{r},t\right)}}{\Phi }_1\left(\mathbf{r},t\right){\widehat{x}}_1$$

(66)

$$\left({\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)·\nabla \right){\Phi }_2\left(\mathbf{r},t\right)=-{\displaystyle \frac{B_2\left(\mathbf{r},t\right)}{v_{\mathrm{rot}}\left(\mathbf{r},t\right)}}{\Phi }_2\left(\mathbf{r},t\right){\widehat{x}}_2$$

(67)

$$\left({\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)·\nabla \right){\Phi }_3\left(\mathbf{r},t\right)=-{\displaystyle \frac{B_3\left(\mathbf{r},t\right)}{w_{\mathrm{rot}}\left(\mathbf{r},t\right)}}{\Phi }_3\left(\mathbf{r},t\right){\widehat{x}}_3$$

(68)

Applying the change of variables ${\Phi }_k\left(\mathbf{r},t\right)=exp\left({\Xi }_k\left(\mathbf{r},t\right)\right)$ to (66)–(68) cancels ${\Phi }_k\left(\mathbf{r},t\right)$ from the RHS of the MCH residues, leaving the system of first-order linear PDEs (69)–(71) coupled only through their mutual mentioning of ${\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)$ and which lack undifferentiated instances of their unknowns.

$$\left({\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)·\nabla \right){\Xi }_1\left(\mathbf{r},t\right)=-{\displaystyle \frac{B_1\left(\mathbf{r},t\right)}{u_{\mathrm{rot}}\left(\mathbf{r},t\right)}}{\widehat{x}}_1$$

(69)

$$\left({\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)·\nabla \right){\Xi }_2\left(\mathbf{r},t\right)=-{\displaystyle \frac{B_2\left(\mathbf{r},t\right)}{v_{\mathrm{rot}}\left(\mathbf{r},t\right)}}{\widehat{x}}_2$$

(70)

$$\left({\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)·\nabla \right){\Xi }_3\left(\mathbf{r},t\right)=-{\displaystyle \frac{B_3\left(\mathbf{r},t\right)}{w_{\mathrm{rot}}\left(\mathbf{r},t\right)}}{\widehat{x}}_3$$

(71)

The system (69)–(71) may be solved using the method of characteristics. Defining the stream function $V_{t^{\prime }}\left(\mathbf{r}\right)$ of the rotational velocity ${\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)$ by the first-order ODE (72) with solution (73) defines characteristic curves of the system.

$${\displaystyle \frac{d{V}_{t^{\prime }}\left(\mathbf{r}\right)}{d{t}^{\prime }}}={\mathbf{v}}_{\mathrm{rot}}\left(V_{t^{\prime }}\left(\mathbf{r}\right),t\right),V_0\left(\mathbf{r}\right)=\mathbf{r}$$

(72)

$$t^{\prime }={\int }_{\mathbf{r}}^{V_{t^{\prime }}\left(\mathbf{r}\right)}{\displaystyle \frac{d\Lambda }{{\mathbf{v}}_{\mathrm{rot}}\left(\Lambda ,t\right)}}$$

(73)

Though (72) and (73) defines characteristic curves, they are insufficiently imbued with information regarding $\Omega$ to be used in practical scenarios requiring a specific $\Omega$. Specifying a backward exit point $T(\Omega ,\mathbf{r},t)$ (74) at which an existing characteristic passing through the specified $(\mathbf{r},t)$ intersects $\Omega$ (74), while substituting (73) gives $T(\Omega ,\mathbf{r},t)$ in terms of the rotational velocity field ${\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)$ (75).

$$T(\Omega ,\mathbf{r},t)=sup\left\{t^{\prime }\ge 0:V_{-s}\left(\mathbf{r}\right)\in \Omega \forall s\in [0,t^{\prime }]\right\}$$

(74)

$$T(\Omega ,\mathbf{r},t)=sup\left\{{\int }_{\mathbf{r}}^{V_{t^{\prime }}\left(\mathbf{r}\right)}{\displaystyle \frac{d\Lambda }{{\mathbf{v}}_{\mathrm{rot}}\left(\Lambda ,t\right)}}\ge 0:V_{-s}\left(\mathbf{r}\right)\in \Omega \forall s\in \left[0,{\int }_{\mathbf{r}}^{V_{t^{\prime }}\left(\mathbf{r}\right)}{\displaystyle \frac{d\Lambda }{{\mathbf{v}}_{\mathrm{rot}}\left(\Lambda ,t\right)}}\right]\right\}$$

(75)

Integrating the RHS of (69) and (70) and (71) along the path of each characteristic from $s=0$ to $s=T(\Omega ,\mathbf{r},t)$ enables all streamlines of ${\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)$ on a simply connected $\Omega$ to be accounted for when calculating ${\Xi }_k\left(\mathbf{r},t\right)$.

$${\Xi }_1\left(\mathbf{r},t\right)={\int }_0^{T(\Omega ,\mathbf{r},t)}-{\displaystyle \frac{B_1\left(V_{-s}\left(\mathbf{r}\right),t\right)}{u_{\mathrm{rot}}\left(V_{-s}\left(\mathbf{r}\right),t\right)}}ds$$

(76)

$${\Xi }_2\left(\mathbf{r},t\right)={\int }_0^{T(\Omega ,\mathbf{r},t)}-{\displaystyle \frac{B_2\left(V_{-s}\left(\mathbf{r}\right),t\right)}{v_{\mathrm{rot}}\left(V_{-s}\left(\mathbf{r}\right),t\right)}}ds$$

(77)

$${\Xi }_3\left(\mathbf{r},t\right)={\int }_0^{T(\Omega ,\mathbf{r},t)}-{\displaystyle \frac{B_3\left(V_{-s}\left(\mathbf{r}\right),t\right)}{w_{\mathrm{rot}}\left(V_{-s}\left(\mathbf{r}\right),t\right)}}ds$$

(78)

Substituting (76)–(78) into the change of variables ${\Phi }_k\left(\mathbf{r},t\right)=exp\left({\Xi }_k\left(\mathbf{r},t\right)\right)$ which gave rise to ${\Xi }_k\left(\mathbf{r},t\right)$ yields expressions (79)–(81) of ${\Phi }_k\left(\mathbf{r},t\right)$ all in terms of ${\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)$ which serve as the template of the MCH reconciliation.

$${\Phi }_1\left(\mathbf{r},t\right)=exp\left({\int }_0^{T(\Omega )}-{\displaystyle \frac{B_1\left(V_{-s}\left(\mathbf{r}\right),t\right)}{u_{\mathrm{rot}}\left(V_{-s}\left(\mathbf{r}\right),t\right)}}ds\right)$$

(79)

$${\Phi }_2\left(\mathbf{r},t\right)=exp\left({\int }_0^{T(\Omega )}-{\displaystyle \frac{B_2\left(V_{-s}\left(\mathbf{r}\right),t\right)}{v_{\mathrm{rot}}\left(V_{-s}\left(\mathbf{r}\right),t\right)}}ds\right)$$

(80)

$${\Phi }_3\left(\mathbf{r},t\right)=exp\left({\int }_0^{T(\Omega )}-{\displaystyle \frac{B_3\left(V_{-s}\left(\mathbf{r}\right),t\right)}{w_{\mathrm{rot}}\left(V_{-s}\left(\mathbf{r}\right),t\right)}}ds\right)$$

(81)

Evaluating (79)–(81) at $t=0$ and at $\mathbf{r}\in \partial \Omega$ yield the initial conditions ${\Phi }_{k_0}\left(\mathbf{r}\right)$ of ${\Phi }_k\left(\mathbf{r},t\right)$ in terms of those of ${\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)$ (82)–(84) and the boundary conditions ${\Phi }_{k_{\Omega }}\left(\mathbf{r},t\right)$ of ${\Phi }_k\left(\mathbf{r},t\right)$ also in terms of those of ${\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)$ (85)–(87) respectively.

$${\Phi }_{1_0}\left(\mathbf{r},t\right)=exp\left({\int }_0^{T(\Omega ,\mathbf{r},0)}-{\displaystyle \frac{B_{1_0}\left(V_{-s}\left(\mathbf{r}\right)\right)}{u_{{\mathrm{rot}}_0}\left(V_{-s}\left(\mathbf{r}\right)\right)}}ds\right)$$

(82)

$${\Phi }_{2_0}\left(\mathbf{r},t\right)=exp\left({\int }_0^{T(\Omega ,\mathbf{r},0)}-{\displaystyle \frac{B_{2_0}\left(V_{-s}\left(\mathbf{r}\right)\right)}{v_{{\mathrm{rot}}_0}\left(V_{-s}\left(\mathbf{r}\right)\right)}}ds\right)$$

(83)

$${\Phi }_{3_0}\left(\mathbf{r},t\right)=exp\left({\int }_0^{T(\Omega ,\mathbf{r},0)}-{\displaystyle \frac{B_{3_0}\left(V_{-s}\left(\mathbf{r}\right)\right)}{w_{{\mathrm{rot}}_0}\left(V_{-s}\left(\mathbf{r}\right)\right)}}ds\right)$$

(84)

$${\Phi }_{1_{\Omega }}\left(\mathbf{r},t\right)=exp\left({\int }_0^{T(\partial \Omega ,\mathbf{r},t)}-{\displaystyle \frac{B_{1_{\Omega }}\left(V_{-s}\left(\mathbf{r}\right)\in \partial \Omega ,t\right)}{u_{{\mathrm{rot}}_{\Omega }}\left(V_{-s}\left(\mathbf{r}\right)\in \partial \Omega ,t\right)}}ds\right)$$

(85)

$${\Phi }_{2_{\Omega }}\left(\mathbf{r},t\right)=exp\left({\int }_0^{T(\partial \Omega ,\mathbf{r},t)}-{\displaystyle \frac{B_{2_{\Omega }}\left(V_{-s}\left(\mathbf{r}\right)\in \partial \Omega ,t\right)}{v_{{\mathrm{rot}}_{\Omega }}\left(V_{-s}\left(\mathbf{r}\right)\in \partial \Omega ,t\right)}}ds\right)$$

(86)

$${\Phi }_{3_{\Omega }}\left(\mathbf{r},t\right)=exp\left({\int }_0^{T(\partial \Omega ,\mathbf{r},t)}-{\displaystyle \frac{B_{3_{\Omega }}\left(V_{-s}\left(\mathbf{r}\right)\in \partial \Omega ,t\right)}{w_{{\mathrm{rot}}_{\Omega }}\left(V_{-s}\left(\mathbf{r}\right)\in \partial \Omega ,t\right)}}ds\right)$$

(87)

Substituting (82)–(29) for $k=1$, $k=2$, and $k=3$ respectively generates formulae for bounded $\{{\Phi }_1\left(\mathbf{r},t\right),{\Phi }_2\left(\mathbf{r},t\right),{\Phi }_3\left(\mathbf{r},t\right)\}$ (88)–(90) which rely on $\{{\mathbf{v}}_{{\mathrm{rot}}_0}\left(\mathbf{r}\right),{\mathbf{v}}_{{\mathrm{rot}}_{\Omega }}\left(\mathbf{r},t\right)\}$ rather than $\{{\Phi }_0\left(\mathbf{r}\right),{\Phi }_{\Omega }\left(\mathbf{r},t\right)\}$. The Cauchy data of the heat equation solutions and of the velocity field is now reconciled such that no heat equation Cauchy data appears in the final expression of the rotational velocity field without being expressed in terms of velocity field Cauchy data.

$${\Phi }_1\left(\mathbf{r},t\right)={\int }_{\Omega }exp\left({\int }_0^{T(\Omega ,{\mathbf{r}}^{\prime },0)}-{\displaystyle \frac{B_{1_0}\left(V_{-s}\left({\mathbf{r}}^{\prime }\right)\right)}{u_{{\mathrm{rot}}_0}\left(V_{-s}\left({\mathbf{r}}^{\prime }\right)\right)}}ds\right)G_{\mathrm{HD}}\left(\mathbf{r},{\mathbf{r}}^{\prime },t\right)d{V}_{{\mathbf{r}}^{\prime }}+$$

$$\nu {\int }_0^t{\oint }_{\partial \Omega }exp\left({\int }_0^{T(\partial \Omega ,{\mathbf{r}}^{\prime },\tau )}-{\displaystyle \frac{B_{1_{\Omega }}\left(V_{-s}\left({\mathbf{r}}^{\prime }\right)\in \partial \Omega ,\tau \right)}{u_{{\mathrm{rot}}_{\Omega }}\left(V_{-s}\left({\mathbf{r}}^{\prime }\right)\in \partial \Omega ,\tau \right)}}ds\right){\displaystyle \frac{\partial G_{\mathrm{HD}}\left(\mathbf{r},{\mathbf{r}}^{\prime },t-\tau \right)}{\partial {\widehat{N}}_{{\mathbf{r}}^{\prime }}}}d{S}_{{\mathbf{r}}^{\prime }}d\tau$$

(88)

$${\Phi }_2\left(\mathbf{r},t\right)={\int }_{\Omega }exp\left({\int }_0^{T(\Omega ,{\mathbf{r}}^{\prime },0)}-{\displaystyle \frac{B_{2_0}\left(V_{-s}\left({\mathbf{r}}^{\prime }\right)\right)}{v_{{\mathrm{rot}}_0}\left(V_{-s}\left({\mathbf{r}}^{\prime }\right)\right)}}ds\right)G_{\mathrm{HD}}\left(\mathbf{r},{\mathbf{r}}^{\prime },t\right)d{V}_{{\mathbf{r}}^{\prime }}+$$

$$\nu {\int }_0^t{\oint }_{\partial \Omega }exp\left({\int }_0^{T(\partial \Omega ,{\mathbf{r}}^{\prime },\tau )}-{\displaystyle \frac{B_{2_{\Omega }}\left(V_{-s}\left({\mathbf{r}}^{\prime }\right)\in \partial \Omega ,\tau \right)}{v_{{\mathrm{rot}}_{\Omega }}\left(V_{-s}\left({\mathbf{r}}^{\prime }\right)\in \partial \Omega ,\tau \right)}}ds\right){\displaystyle \frac{\partial G_{\mathrm{HD}}\left(\mathbf{r},{\mathbf{r}}^{\prime },t-\tau \right)}{\partial {\widehat{N}}_{{\mathbf{r}}^{\prime }}}}d{S}_{{\mathbf{r}}^{\prime }}d\tau$$

(89)

$${\Phi }_3\left(\mathbf{r},t\right)={\int }_{\Omega }exp\left({\int }_0^{T(\Omega ,{\mathbf{r}}^{\prime },0)}-{\displaystyle \frac{B_{3_0}\left(V_{-s}\left({\mathbf{r}}^{\prime }\right)\right)}{w_{{\mathrm{rot}}_0}\left(V_{-s}\left({\mathbf{r}}^{\prime }\right)\right)}}ds\right)G_{\mathrm{HD}}\left(\mathbf{r},{\mathbf{r}}^{\prime },t\right)d{V}_{{\mathbf{r}}^{\prime }}+$$

$$\nu {\int }_0^t{\oint }_{\partial \Omega }exp\left({\int }_0^{T(\partial \Omega ,{\mathbf{r}}^{\prime },\tau )}-{\displaystyle \frac{B_{3_{\Omega }}\left(V_{-s}\left({\mathbf{r}}^{\prime }\right)\in \partial \Omega ,\tau \right)}{w_{{\mathrm{rot}}_{\Omega }}\left(V_{-s}\left({\mathbf{r}}^{\prime }\right)\in \partial \Omega ,\tau \right)}}ds\right){\displaystyle \frac{\partial G_{\mathrm{HD}}\left(\mathbf{r},{\mathbf{r}}^{\prime },t-\tau \right)}{\partial {\widehat{N}}_{{\mathbf{r}}^{\prime }}}}d{S}_{{\mathbf{r}}^{\prime }}d\tau$$

(90)

6.6. Final Statement of Rotational Velocity Field as MCH Subject

The rotational velocity field ${\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)$ is expressed (91) where the ${\Phi }_k\left(\mathbf{r},t\right)$ are given by (94), in which $\{{\mathbf{v}}_{k_{{\mathrm{rot}}_0}}\left(\mathbf{r},t\right),{\mathbf{v}}_{k_{{\mathrm{rot}}_{\Omega }}}\left(\mathbf{r},t\right)\}$ is the kth component of $\{{\mathbf{v}}_{{\mathrm{rot}}_0}\left(\mathbf{r},t\right),{\mathbf{v}}_{{\mathrm{rot}}_{\Omega }}\left(\mathbf{r},t\right)\}$ and $\{\alpha \left(\mathbf{r},t\right),\beta \left(\mathbf{r},t\right)\}$ are given as solutions of the biquadratic system (92) and (93). The backward exit times and characteristic definitions for the initial conditions and boundary conditions of (94) are found by evaluating the expressions (74) and (75) at $t=0$ (95) and (96) and on $\mathbf{r}\in \partial \Omega$ (97) and (98) respectively. The boundary self-similarity terms $\{B_1\left(\mathbf{r},t\right),B_2\left(\mathbf{r},t\right),B_3\left(\mathbf{r},t\right)\}$ are given in vector form (99) for brevity.

$${\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)=\pm {\displaystyle \frac{1{\widehat{x}}_1+\frac{\alpha \left(\mathbf{r},t\right)}{\beta \left(\mathbf{r},t\right)}{\widehat{x}}_2+\frac{\alpha \left(\mathbf{r},t\right)}{{\beta }^2\left(\mathbf{r},t\right)}{\widehat{x}}_3}{{\left(\frac{{\nabla }_x{\Phi }_x\left(\mathbf{r},t\right)}{{\nabla }_y{\Phi }_x\left(\mathbf{r},t\right)+{\nabla }_z{\Phi }_x\left(\mathbf{r},t\right)\alpha \left(\mathbf{r},t\right)-B_x\left(\mathbf{r},t\right)\beta \left(\mathbf{r},t\right)}\right)}^{-1/2}}}$$

(91)

$${\nabla }_1{\Phi }_1\left(\mathbf{r},t\right){\nabla }_3{\Phi }_2\left(\mathbf{r},t\right){\alpha }^2\left(\mathbf{r},t\right)+\left({\nabla }_1{\Phi }_1\left(\mathbf{r},t\right){\nabla }_2{\Phi }_2\left(\mathbf{r},t\right)-{\nabla }_1{\Phi }_1\left(\mathbf{r},t\right)B_2\left(\mathbf{r},t\right)\beta \left(\mathbf{r},t\right)-\right.$$

$$\left.{\nabla }_1{\Phi }_2\left(\mathbf{r},t\right){\nabla }_3{\Phi }_1\left(\mathbf{r},t\right)\right)\alpha \left(\mathbf{r},t\right)+{\nabla }_1{\Phi }_2\left(\mathbf{r},t\right)B_1\left(\mathbf{r},t\right)\beta \left(\mathbf{r},t\right)-{\nabla }_1{\Phi }_2\left(\mathbf{r},t\right){\nabla }_2{\Phi }_1\left(\mathbf{r},t\right)=0$$

(92)

$${\nabla }_1{\Phi }_1\left(\mathbf{r},t\right)B_3\left(\mathbf{r},t\right){\beta }^2\left(\mathbf{r},t\right)-\left({\nabla }_1{\Phi }_1\left(\mathbf{r},t\right){\nabla }_2{\Phi }_3\left(\mathbf{r},t\right)+{\nabla }_1{\Phi }_1\left(\mathbf{r},t\right){\nabla }_3{\Phi }_3\left(\mathbf{r},t\right)\alpha \left(\mathbf{r},t\right)+\right.$$

$$\left.{\nabla }_1{\Phi }_3\left(\mathbf{r},t\right)B_1\left(\mathbf{r},t\right)\right)\beta \left(\mathbf{r},t\right)+{\nabla }_1{\Phi }_3\left(\mathbf{r},t\right){\nabla }_3{\Phi }_1\left(\mathbf{r},t\right)\alpha \left(\mathbf{r},t\right)+{\nabla }_1{\Phi }_3\left(\mathbf{r},t\right){\nabla }_2{\Phi }_1\left(\mathbf{r},t\right)=0$$

(93)

$${\Phi }_k\left(\mathbf{r},t\right)={\int }_{\Omega }exp\left({\int }_0^{T(\Omega ,{\mathbf{r}}^{\prime },0)}-{\displaystyle \frac{B_{k_0}\left(V_{-s}\left({\mathbf{r}}^{\prime }\right)\right)}{{\mathbf{v}}_{k_{{\mathrm{rot}}_0}}\left(V_{-s}\left({\mathbf{r}}^{\prime }\right)\right)}}ds\right)G_{\mathrm{HD}}\left(\mathbf{r},{\mathbf{r}}^{\prime },t\right)d{V}_{{\mathbf{r}}^{\prime }}+$$

$$\nu {\int }_0^t{\oint }_{\partial \Omega }exp\left({\int }_0^{T(\partial \Omega ,{\mathbf{r}}^{\prime },\tau )}-{\displaystyle \frac{B_{k_{\Omega }}\left(V_{-s}\left({\mathbf{r}}^{\prime }\right)\in \partial \Omega ,\tau \right)}{{\mathbf{v}}_{k_{{\mathrm{rot}}_{\Omega }}}\left(V_{-s}\left({\mathbf{r}}^{\prime }\right)\in \partial \Omega ,\tau \right)}}ds\right){\displaystyle \frac{\partial G_{\mathrm{HD}}\left(\mathbf{r},{\mathbf{r}}^{\prime },t-\tau \right)}{\partial {\widehat{N}}_{{\mathbf{r}}^{\prime }}}}d{S}_{{\mathbf{r}}^{\prime }}d\tau$$

(94)

$$T(\Omega ,\mathbf{r},0)=sup\left\{{\int }_{\mathbf{r}}^{V_{t^{\prime }}\left(\mathbf{r}\right)}{\displaystyle \frac{d\Lambda }{{\mathbf{v}}_{{\mathrm{rot}}_0}\left(\Lambda \right)}}\ge 0:V_{-s}\left(\mathbf{r}\right)\in \Omega \forall s\in \left[0,{\int }_{\mathbf{r}}^{V_{t^{\prime }}\left(\mathbf{r}\right)}{\displaystyle \frac{d\Lambda }{{\mathbf{v}}_{{\mathrm{rot}}_0}\left(\Lambda \right)}}\right]\right\}$$

(95)

$$t_{1_0}^{\prime }={\int }_{\mathbf{r}}^{V_{t^{\prime }}\left(\mathbf{r}\right)}{\displaystyle \frac{d\Lambda }{{\mathbf{v}}_{{\mathrm{rot}}_0}\left(\Lambda \right)}}$$

(96)

$$T(\partial \Omega ,\mathbf{r},t)=sup\left\{{\int }_{\mathbf{r}}^{V_{t^{\prime }}\left(\mathbf{r}\right)}{\displaystyle \frac{d\Lambda }{{\mathbf{v}}_{{\mathrm{rot}}_{\Omega }}\left(\Lambda ,t\right)}}\ge 0:V_{-s}\left(\mathbf{r}\right)\in \partial \Omega \forall s\in \left[0,{\int }_{\mathbf{r}}^{V_{t^{\prime }}\left(\mathbf{r}\right)}{\displaystyle \frac{d\Lambda }{{\mathbf{v}}_{{\mathrm{rot}}_{\Omega }}\left({\Lambda }_1,t\right)}}\right]\right\}$$

(97)

$$t_{1_{\Omega }}^{\prime }={\int }_{\mathbf{r}}^{V_{t^{\prime }}\left(\mathbf{r}\right)}{\displaystyle \frac{d\Lambda }{{\mathbf{v}}_{{\mathrm{rot}}_{\Omega }}\left(\Lambda ,t\right)}}$$

(98)

$$B_1\left(\mathbf{r},t\right){\widehat{x}}_1+B_2\left(\mathbf{r},t\right){\widehat{x}}_2+B_3\left(\mathbf{r},t\right){\widehat{x}}_3={\int }_{\Omega \left(t\right)}{\displaystyle \frac{{\nabla }^2\mathbf{f}\left({\mathbf{r}}^{\prime },t\right)}{4\pi }}{G}_p(\mathbf{r},{\mathbf{r}}^{\prime })d{V}_{{\mathbf{r}}^{\prime }}+{\displaystyle \frac{1}{2\pi }}{\oint }_{\partial \Omega \left(t\right)}$$

$$\left(\nu {\nabla }^2{\mathbf{v}}_{{\mathrm{rot}}_{\Omega \left(t\right)}}\left({\mathbf{r}}^{\prime },t\right)-{\displaystyle \frac{\partial {\mathbf{v}}_{{\mathrm{rot}}_{\Omega \left(t\right)}}\left({\mathbf{r}}^{\prime },t\right)}{\partial t}}-\left({\mathbf{v}}_{{\mathrm{rot}}_{\Omega \left(t\right)}}\left({\mathbf{r}}^{\prime },t\right)·\nabla \right){\mathbf{v}}_{{\mathrm{rot}}_{\Omega \left(t\right)}}\left({\mathbf{r}}^{\prime },t\right)\right){\displaystyle \frac{\partial G_{\mathrm{PD}}(\mathbf{r},{\mathbf{r}}^{\prime })}{\partial {\widehat{N}}_{{\mathbf{r}}^{\prime }}}}{S}_{{\mathbf{r}}^{\prime }}$$

(99)

Since a rotational velocity field ${\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)$ was found (91)–(99) that satisfies (1) and (2) by its MCH-based construction and may therefore be used in combination with an irrotational velocity field ${\mathbf{v}}_{\mathrm{irr}}\left(\mathbf{r},t\right)$ (Lemma 1) to reconstruct a total velocity field $\mathbf{v}\left(\mathbf{r},t\right)$ satisfying (1) and (2), Claim 3 stands with Lemma 2 specifying the result of its proof.

Lemma 2. 

If $\mathbf{v}:\left[(\mathbf{r}\in \Omega )\times t\in {\mathbb{R}}_{\ge 0}\right]\to {\mathbb{R}}^3$ satisfy (1) and (2) given $\Omega \subseteq {\mathbb{R}}^3$ and $\mathbf{f}:\left[(\mathbf{r}\in \Omega )\times t\in {\mathbb{R}}_{\ge 0}\right]\to {\mathbb{R}}^3$, then there exists a rotational yet solenoidal velocity field ${\mathbf{v}}_{\mathrm{rot}}(\mathbf{r},t)$ (91)–(99) constructible from initial data ${\mathbf{v}}_0\left(\mathbf{r}\right)=\mathbf{v}(\mathbf{r},0)$ and boundary data ${\mathbf{v}}_{\Omega }(\mathbf{r},t)=\mathbf{v}(\mathbf{r}\in \Omega ,t)$ of $\mathbf{v}(\mathbf{r},t)$, from which the total $\mathbf{v}(\mathbf{r},t)$ is recovered through the Helmholtz decomposition (3) and ${\mathbf{v}}_{\mathrm{irr}}\left(\mathbf{r},t\right)$ (Lemma 1.)

7. Final Assembly of the Velocity Field

The irrotational part ${\mathbf{v}}_{\mathrm{irr}}\left(\mathbf{r},t\right)$ (Lemma 1) and rotational part ${\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)$ (Lemma 2) of the velocity field $\mathbf{v}\left(\mathbf{r},t\right)$ are now in hand (15) and (91)–(99). Per the Helmholtz decomposition (3), the sum of ${\mathbf{v}}_{\mathrm{irr}}\left(\mathbf{r},t\right)$ (15) and ${\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)$ (91) equals $\mathbf{v}\left(\mathbf{r},t\right)$ (100), where the scalars $\alpha \left(\mathbf{r},t\right)$ and $\beta \left(\mathbf{r},t\right)$ are determined by solving the biquadratic system (101) and (102). The heat equation solutions as MCH kernels ${\Phi }_k\left(\mathbf{r},t\right)|k=\{1,2,3\}$ required to solve (101) and (102) are given by (103). Substituting (12) and (14) into (95)–(99) expresses the equally necessary depressurized momentum boundary terms $B_k\left(\mathbf{r},t\right)|k=\{1,2,3\}$, initial and boundary backward exit times $\left\{T\right(\Omega ,\mathbf{r},0),T(\partial \Omega ,\mathbf{r},t\left)\right\}$ (104) and (106), and characteristic definitions $\{V_{t_0^{\prime }}\left(\mathbf{r}\right),V_{t_{\Omega }^{\prime }}\left(\mathbf{r},t\right)\}$ (105) and (107) on the boundary and at the initial respectively, wholly in terms of ${\mathbf{v}}_{\Omega }\left(\mathbf{r},t\right)$ (108). Substituting (108) and (103) into (101) and (102) and the results into (100) yields the closed-form velocity field $\mathbf{v}\left(\mathbf{r},t\right)$ of the incompressible Navier-Stokes equations (1) and (2) on simply connected $\Omega \subseteq {\mathbb{R}}^3$.

$$\mathbf{v}\left(\mathbf{r},t\right)=\nabla {\oint }_{\partial \Omega }\left({\oint }_{\partial \Omega }\widehat{\mathbf{n}}·{\displaystyle \frac{{\mathbf{v}}_{\Omega }\left({\mathbf{r}}^{\prime \prime },t\right)}{8{\pi }^2|{\mathbf{r}}^{\prime }-{\mathbf{r}}^{\prime \prime }|}}d{S}_{{\mathbf{r}}^{\prime \prime }}-{\int }_{\Omega }{\displaystyle \frac{{\nabla }^{\prime }·{\mathbf{v}}_{\Omega }\left({\mathbf{r}}^{\prime \prime },t\right)}{8{\pi }^2|{\mathbf{r}}^{\prime }-{\mathbf{r}}^{\prime \prime }|}}d{V}_{{\mathbf{r}}^{\prime \prime }}\right){\displaystyle \frac{\partial G_{\mathrm{PD}}\left(\mathbf{r},{\mathbf{r}}^{\prime },t\right)}{\partial {\widehat{N}}_{{\mathbf{r}}^{\prime }}}}d{S}_{{\mathbf{r}}^{\prime }}$$

$$\pm \left(1{\widehat{x}}_1+{\displaystyle \frac{\alpha \left(\mathbf{r},t\right)}{\beta \left(\mathbf{r},t\right)}}{\widehat{x}}_2+{\displaystyle \frac{\alpha \left(\mathbf{r},t\right)}{{\beta }^2\left(\mathbf{r},t\right)}}{\widehat{x}}_3\right)\sqrt{{\displaystyle \frac{{\nabla }_x{\Phi }_x\left(\mathbf{r},t\right)}{{\nabla }_y{\Phi }_x\left(\mathbf{r},t\right)+{\nabla }_z{\Phi }_x\left(\mathbf{r},t\right)\alpha \left(\mathbf{r},t\right)-B_x\left(\mathbf{r},t\right)\beta \left(\mathbf{r},t\right)}}}$$

(100)

$${\nabla }_1{\Phi }_1\left(\mathbf{r},t\right){\nabla }_3{\Phi }_2\left(\mathbf{r},t\right){\alpha }^2\left(\mathbf{r},t\right)+\left({\nabla }_1{\Phi }_1\left(\mathbf{r},t\right){\nabla }_2{\Phi }_2\left(\mathbf{r},t\right)-{\nabla }_1{\Phi }_1\left(\mathbf{r},t\right)B_2\left(\mathbf{r},t\right)\beta \left(\mathbf{r},t\right)-\right.$$

$$\left.{\nabla }_1{\Phi }_2\left(\mathbf{r},t\right){\nabla }_3{\Phi }_1\left(\mathbf{r},t\right)\right)\alpha \left(\mathbf{r},t\right)+{\nabla }_1{\Phi }_2\left(\mathbf{r},t\right)B_1\left(\mathbf{r},t\right)\beta \left(\mathbf{r},t\right)-{\nabla }_1{\Phi }_2\left(\mathbf{r},t\right){\nabla }_2{\Phi }_1\left(\mathbf{r},t\right)=0$$

(101)

$${\nabla }_1{\Phi }_1\left(\mathbf{r},t\right)B_3\left(\mathbf{r},t\right){\beta }^2\left(\mathbf{r},t\right)-\left({\nabla }_1{\Phi }_1\left(\mathbf{r},t\right){\nabla }_2{\Phi }_3\left(\mathbf{r},t\right)+{\nabla }_1{\Phi }_1\left(\mathbf{r},t\right){\nabla }_3{\Phi }_3\left(\mathbf{r},t\right)\alpha \left(\mathbf{r},t\right)+\right.$$

$$\left.{\nabla }_1{\Phi }_3\left(\mathbf{r},t\right)B_1\left(\mathbf{r},t\right)\right)\beta \left(\mathbf{r},t\right)+{\nabla }_1{\Phi }_3\left(\mathbf{r},t\right){\nabla }_3{\Phi }_1\left(\mathbf{r},t\right)\alpha \left(\mathbf{r},t\right)+{\nabla }_1{\Phi }_3\left(\mathbf{r},t\right){\nabla }_2{\Phi }_1\left(\mathbf{r},t\right)=0$$

(102)

$${\Phi }_k\left(\mathbf{r},t\right)={\int }_{\Omega }exp\left({\int }_0^{T(\Omega ,{\mathbf{r}}^{\prime },0)}-B_{k_0}\left(V_{-s}\left({\mathbf{r}}^{\prime }\right)\right)\left({\left(\nabla \times {\int }_{\Omega }{\displaystyle \frac{{\nabla }^{\prime }\times {\mathbf{v}}_0\left(V_{-s}\left({\mathbf{r}}^{\prime \prime }\right)\right)}{4\pi |V_{-s}\left({\mathbf{r}}^{\prime }\right)-V_{-s}\left({\mathbf{r}}^{\prime \prime }\right)|}}d{V}_{V_{-s}\left({\mathbf{r}}^{\prime \prime }\right)}\right)}_k-\right.\right.$$

$$\left.{\left.{\left(\nabla \times {\oint }_{\partial \Omega }\widehat{\mathbf{n}}\times {\displaystyle \frac{{\mathbf{v}}_0\left(V_{-s}\left({\mathbf{r}}^{\prime \prime }\right)\right)}{4\pi |V_{-s}\left({\mathbf{r}}^{\prime }\right)-V_{-s}\left({\mathbf{r}}^{\prime \prime }\right)|}}d{S}_{V_{-s}\left({\mathbf{r}}^{\prime \prime }\right)}\right)}_k\right)}^{-1}ds\right)G_{\mathrm{HD}}\left(\mathbf{r},{\mathbf{r}}^{\prime },t\right)d{V}_{{\mathbf{r}}^{\prime }}+$$

$$\nu {\int }_0^t{\oint }_{\partial \Omega }exp\left({\int }_0^{T(\partial \Omega ,{\mathbf{r}}^{\prime },\tau )}-B_{k_{\Omega }}\left(V_{-s}\left({\mathbf{r}}^{\prime }\right)\in \partial \Omega ,\tau \right)\left({\left(\nabla \times {\int }_{\Omega }{\displaystyle \frac{{\nabla }^{\prime }\times {\mathbf{v}}_{\Omega }\left(V_{-s}\left({\mathbf{r}}^{\prime \prime }\right)\in \partial \Omega ,\tau \right)}{4\pi |V_{-s}\left({\mathbf{r}}^{\prime }\right)-V_{-s}\left({\mathbf{r}}^{\prime \prime }\right)|}}d{V}_{V_{-s}\left({\mathbf{r}}^{\prime \prime }\right)}\right)}_k-\right.\right.$$

$$\left.{\left.{\left(\nabla \times {\oint }_{\partial \Omega }\widehat{\mathbf{n}}\times {\displaystyle \frac{{\mathbf{v}}_{\Omega }\left(V_{-s}\left({\mathbf{r}}^{\prime \prime }\right)\in \partial \Omega ,\tau \right)}{4\pi |V_{-s}\left({\mathbf{r}}^{\prime }\right)-V_{-s}\left({\mathbf{r}}^{\prime \prime }\right)|}}d{S}_{V_{-s}\left({\mathbf{r}}^{\prime \prime }\right)}\right)}_k\right)}^{-1}ds\right){\displaystyle \frac{\partial G_{\mathrm{HD}}\left(\mathbf{r},{\mathbf{r}}^{\prime },t-\tau \right)}{\partial {\widehat{N}}_{{\mathbf{r}}^{\prime }}}}d{S}_{{\mathbf{r}}^{\prime }}d\tau$$

(103)

$$T(\partial \Omega ,\mathbf{r},t)=sup\left\{{\int }_{\mathbf{r}}^{V_{t_{\Omega }^{\prime }}\left(\mathbf{r}\right)}{\left(\nabla \times {\int }_{\Omega }{\displaystyle \frac{{\nabla }^{\prime }\times {\mathbf{v}}_{\Omega }\left({\Lambda }^{\prime },t\right)}{4\pi |\Lambda -{\Lambda }^{\prime }|}}d{V}_{{\Lambda }^{\prime }}-\nabla \times {\oint }_{\partial \Omega }\widehat{\mathbf{n}}\times {\displaystyle \frac{{\mathbf{v}}_{\Omega }\left({\Lambda }^{\prime },t\right)}{4\pi |\Lambda -{\Lambda }^{\prime }|}}d{S}_{{\Lambda }^{\prime }}\right)}^{-1}\right.$$

$$d\Lambda \ge 0:V_{-s}\left(\mathbf{r}\right)\in \partial \Omega \forall s\in \left[0,\left({\int }_{\mathbf{r}}^{V_{t_{\Omega }^{\prime }}\left(\mathbf{r}\right)}\left(\nabla \times {\int }_{\Omega }{\displaystyle \frac{{\nabla }^{\prime }\times {\mathbf{v}}_{\Omega }\left({\Lambda }^{\prime },t\right)}{4\pi |\Lambda -{\Lambda }^{\prime }|}}d{V}_{{\Lambda }^{\prime }}-\right.\right.\right.$$

$$\left.\left.\left.{\left.\nabla \times {\oint }_{\partial \Omega }\widehat{\mathbf{n}}\times {\displaystyle \frac{{\mathbf{v}}_{\Omega }\left({\Lambda }^{\prime },t\right)}{4\pi |\Lambda -{\Lambda }^{\prime }|}}d{S}_{{\Lambda }^{\prime }}\right)}^{-1}d\Lambda \right)\right]\right\}$$

(104)

$$t_{\Omega }^{\prime }={\int }_{\mathbf{r}}^{V_{t_{\Omega }^{\prime }}\left(\mathbf{r}\right)}{\left(\nabla \times {\int }_{\Omega }{\displaystyle \frac{{\nabla }^{\prime }\times {\mathbf{v}}_{\Omega }\left({\Lambda }^{\prime },t\right)}{4\pi |\Lambda -{\Lambda }^{\prime }|}}d{V}_{{\Lambda }^{\prime }}-\nabla \times {\oint }_{\partial \Omega }\widehat{\mathbf{n}}\times {\displaystyle \frac{{\mathbf{v}}_{\Omega }\left({\Lambda }^{\prime },t\right)}{4\pi |\Lambda -{\Lambda }^{\prime }|}}d{S}_{{\Lambda }^{\prime }}\right)}^{-1}d\Lambda$$

(105)

$$T(\Omega ,\mathbf{r},0)=sup\left\{{\int }_{\mathbf{r}}^{V_{t_0^{\prime }}\left(\mathbf{r}\right)}{\left(\nabla \times {\int }_{\Omega }{\displaystyle \frac{{\nabla }^{\prime }\times {\mathbf{v}}_0\left({\Lambda }^{\prime }\right)}{4\pi |\Lambda -{\Lambda }^{\prime }|}}d{V}_{{\Lambda }^{\prime }}-\nabla \times {\oint }_{\partial \Omega }\widehat{\mathbf{n}}\times {\displaystyle \frac{{\mathbf{v}}_0\left({\Lambda }^{\prime }\right)}{4\pi |\Lambda -{\Lambda }^{\prime }|}}d{S}_{{\Lambda }^{\prime }}\right)}^{-1}\right.$$

$$d\Lambda \ge 0:V_{-s}\left(\mathbf{r}\right)\in \partial \Omega \forall s\in \left[0,\left({\int }_{\mathbf{r}}^{V_{t_0^{\prime }}\left(\mathbf{r}\right)}\left(\nabla \times {\int }_{\Omega }{\displaystyle \frac{{\nabla }^{\prime }\times {\mathbf{v}}_0\left({\Lambda }^{\prime }\right)}{4\pi |\Lambda -{\Lambda }^{\prime }|}}d{V}_{{\Lambda }^{\prime }}-\right.\right.\right.$$

$$\left.\left.\left.{\left.\nabla \times {\oint }_{\partial \Omega }\widehat{\mathbf{n}}\times {\displaystyle \frac{{\mathbf{v}}_0\left({\Lambda }^{\prime }\right)}{4\pi |\Lambda -{\Lambda }^{\prime }|}}d{S}_{{\Lambda }^{\prime }}\right)}^{-1}d\Lambda \right)\right]\right\}$$

(106)

$$t_0^{\prime }={\int }_{\mathbf{r}}^{V_{t_0^{\prime }}\left(\mathbf{r}\right)}{\left(\nabla \times {\int }_{\Omega }{\displaystyle \frac{{\nabla }^{\prime }\times {\mathbf{v}}_0\left({\Lambda }^{\prime }\right)}{4\pi |\Lambda -{\Lambda }^{\prime }|}}d{V}_{{\Lambda }^{\prime }}-\nabla \times {\oint }_{\partial \Omega }\widehat{\mathbf{n}}\times {\displaystyle \frac{{\mathbf{v}}_0\left({\Lambda }^{\prime }\right)}{4\pi |\Lambda -{\Lambda }^{\prime }|}}d{S}_{{\Lambda }^{\prime }}\right)}^{-1}d\Lambda$$

(107)

$$B_1{\widehat{x}}_1+B_2{\widehat{x}}_2+B_3{\widehat{x}}_3={\int }_{\Omega }{\displaystyle \frac{{\nabla }^2\mathbf{f}\left({\mathbf{r}}^{\prime },t\right)}{4\pi }}{G}_{\mathrm{PD}}(\mathbf{r},{\mathbf{r}}^{\prime })d{V}_{{\mathbf{r}}^{\prime }}+{\oint }_{\partial \Omega }\left(\nu {\nabla }^2\left(\nabla \times {\int }_{\Omega }{\displaystyle \frac{{\nabla }^{\prime }\times {\mathbf{v}}_{\Omega }\left({\mathbf{r}}^{\prime \prime },t\right)}{8{\pi }^2\left|{\mathbf{r}}^{\prime }-{\mathbf{r}}^{\prime \prime }\right|}}d{V}_{{\mathbf{r}}^{\prime \prime }}-\right.\right.$$

$$\left.\nabla \times {\oint }_{\partial \Omega }\widehat{\mathbf{n}}\times {\displaystyle \frac{{\mathbf{v}}_{\Omega }\left({\mathbf{r}}^{\prime \prime },t\right)}{8{\pi }^2|{\mathbf{r}}^{\prime }-{\mathbf{r}}^{\prime \prime }|}}d{S}_{{\mathbf{r}}^{\prime \prime }}\right)-{\displaystyle \frac{\partial }{\partial t}}\left(\nabla \times {\int }_{\Omega }{\displaystyle \frac{{\nabla }^{\prime }\times {\mathbf{v}}_{\Omega }\left({\mathbf{r}}^{\prime \prime },t\right)}{8{\pi }^2\left|{\mathbf{r}}^{\prime }-{\mathbf{r}}^{\prime \prime }\right|}}d{V}_{{\mathbf{r}}^{\prime \prime }}-\nabla \times {\oint }_{\partial \Omega }\widehat{\mathbf{n}}\times {\displaystyle \frac{{\mathbf{v}}_{\Omega }\left({\mathbf{r}}^{\prime \prime },t\right)}{8{\pi }^2|{\mathbf{r}}^{\prime }-{\mathbf{r}}^{\prime \prime }|}}d{S}_{{\mathbf{r}}^{\prime \prime }}\right)$$

$$-\left(\left(\nabla \times {\int }_{\Omega }{\displaystyle \frac{{\nabla }^{\prime }\times {\mathbf{v}}_{\Omega }\left({\mathbf{r}}^{\prime \prime },t\right)}{8{\pi }^2\left|{\mathbf{r}}^{\prime }-{\mathbf{r}}^{\prime \prime }\right|}}d{V}_{{\mathbf{r}}^{\prime \prime }}-\nabla \times {\oint }_{\partial \Omega }\widehat{\mathbf{n}}\times {\displaystyle \frac{{\mathbf{v}}_{\Omega }\left({\mathbf{r}}^{\prime \prime },t\right)}{8{\pi }^2|{\mathbf{r}}^{\prime }-{\mathbf{r}}^{\prime \prime }|}}d{S}_{{\mathbf{r}}^{\prime \prime }}\right)·\nabla \right)\left(\nabla \times {\int }_{\Omega }{\displaystyle \frac{{\nabla }^{\prime }\times {\mathbf{v}}_{\Omega }\left({\mathbf{r}}^{\prime \prime },t\right)}{8{\pi }^2\left|{\mathbf{r}}^{\prime }-{\mathbf{r}}^{\prime \prime }\right|}}d{V}_{{\mathbf{r}}^{\prime \prime }}\right.$$

$$\left.\left.-\nabla \times {\oint }_{\partial \Omega }\widehat{\mathbf{n}}\times {\displaystyle \frac{{\mathbf{v}}_{\Omega }\left({\mathbf{r}}^{\prime \prime },t\right)}{8{\pi }^2|{\mathbf{r}}^{\prime }-{\mathbf{r}}^{\prime \prime }|}}d{S}_{{\mathbf{r}}^{\prime \prime }}\right)\right){\displaystyle \frac{\partial G_{\mathrm{PD}}(\mathbf{r},{\mathbf{r}}^{\prime })}{\partial {\widehat{N}}_{{\mathbf{r}}^{\prime }}}}d{S}_{{\mathbf{r}}^{\prime }}$$

(108)

8. Recovery of the Pressure Field

The velocity field $\mathbf{v}\left(\mathbf{r},t\right)$ is now known (100)–(108) with no direct mentioning of the pressure field $P\left(\mathbf{r},t\right)$. The pressure field is therefore capable of being entirely recovered from the velocity field. Taking the divergence of (1) and isolating the pressure gradient $\nabla P\left(\mathbf{r},t\right)$ cancels the viscous term by the identity $\nu \nabla ·\left({\nabla }^2\mathbf{v}\right)=\nu {\nabla }^2\left(\nabla ·\mathbf{v}\right)=0\forall \mathbf{v}:{\mathbb{R}}^3\to {\mathbb{R}}^3$ and establishes a Poisson equation (109) in $P\left(\mathbf{r},t\right)$ by the identity $\nabla ·\left(\nabla P\right)={\nabla }^{2}P\forall P:{\mathbb{R}}^3\to \mathbb{R}$. This pressure-Poisson formulation of (1) is standard in the literature [62].

$$\nabla ·\left(\nabla P\left(\mathbf{r},t\right)\right)=\rho \nabla ·\mathbf{f}\left(\mathbf{r},t\right)-\rho \nabla ·{\displaystyle \frac{D\mathbf{v}\left(\mathbf{r},t\right)}{Dt}}-\nabla ·\left(\nu {\nabla }^2\mathbf{v}\left(\mathbf{r},t\right)\right)={\nabla }^{2}P\left(\mathbf{r},t\right)\Rightarrow$$

$${\nabla }^{2}P\left(\mathbf{r},t\right)=\rho \nabla ·\mathbf{f}\left(\mathbf{r},t\right)-\rho \nabla ·{\displaystyle \frac{D\mathbf{v}\left(\mathbf{r},t\right)}{Dt}}$$

(109)

Solving the pressure-Poisson equation (109) using a standard Green’s function convolution method [63] for Dirichlet boundary data defines the scalar pressure field (110) as a function of the velocity field. The only requirement of this arrangement is the existence of a Dirichlet boundary condition in pressure $P_{\Omega }\left(\mathbf{r},t\right)$ which may be modified to accommodate a Neumann or mixed boundary condition if needed. The Dirichlet case is shown here for brevity, where $G_{\mathrm{PD}}\left(\mathbf{r},{\mathbf{r}}^{\prime }\right)$ is the Green’s function associated with the Poisson equation in three dimensions under Dirichlet boundary conditions.

$$P\left(\mathbf{r},t\right)={\int }_{\Omega }{\displaystyle \frac{\rho \nabla ·\mathbf{f}\left({\mathbf{r}}^{\prime },t\right)-\rho \nabla ·\frac{D\mathbf{v}\left({\mathbf{r}}^{\prime },t\right)}{Dt}}{4\pi {\left(G_{\mathrm{PD}}\left(\mathbf{r},{\mathbf{r}}^{\prime }\right)\right)}^{-1}}}d{V}_{{\mathbf{r}}^{\prime }}-{\oint }_{\partial \Omega }{\displaystyle \frac{P_{\Omega }\left({\mathbf{r}}^{\prime },t\right)}{2\pi }}{\displaystyle \frac{\partial G_{\mathrm{PD}}\left(\mathbf{r},{\mathbf{r}}^{\prime }\right)}{\partial {\widehat{N}}_{{\mathbf{r}}^{\prime }}}}d{S}_{{\mathbf{r}}^{\prime }}$$

(110)

Because both the total vector velocity field $\mathbf{v}\left(\mathbf{r},t\right)$ (100)–(108) and the scalar pressure field $P\left(\mathbf{r},t\right)$ (110) are now known as functions of $(x_1,x_2,x_3,t)$ without neglecting any nonlinear effects, the equation set (100)–(108) and (110) together constitute the explicit general solution of the incompressible Navier-Stokes equations (1) and (2) for arbitrary flow geometry $\Omega \subseteq {\mathbb{R}}^3$ and arbitrary Dirichlet initial and boundary conditions $\{{\mathbf{v}}_0\left(\mathbf{r}\right),{\mathbf{v}}_{\Omega }\left(\mathbf{r},t\right),P_{\Omega }\left(\mathbf{r},t\right)\}$. Since the components (53)–(55) of ${\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)$ contain the same $\alpha \left(\mathbf{r},t\right)$ and $\beta \left(\mathbf{r},t\right)$, which possess a combined maximum of eight unique real values per $\left\{{\mathbf{v}}_{{\mathrm{rot}}_{\Omega \left(t\right)}}\left(\mathbf{r},t\right),{\mathbf{v}}_{{\mathrm{rot}}_0}\left(\mathbf{r}\right)\right\}$ due to the fact that they are the roots of a biquadratic algebraic system, and the only multi-valuedness of ${\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)$ apart from that of $\alpha \left(\mathbf{r},t\right)$ and $\beta \left(\mathbf{r},t\right)$ consists in its choice of square root shared across each of its three components, the maximum number of unique ${\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)$ for each $\left\{{\mathbf{v}}_{{\mathrm{rot}}_{\Omega \left(t\right)}}\left(\mathbf{r},t\right),{\mathbf{v}}_{{\mathrm{rot}}_0}\left(\mathbf{r}\right)\right\}$ equals $2·2^3=16$. Since ${\mathbf{v}}_{\mathrm{irr}}\left(\mathbf{r},t\right)$ lacks multi-valuedness, the number of unique bounded $\mathbf{v}\left(\mathbf{r},t\right)$ that (1) and (2) are mathematically capable of generating from the same set of initial and boundary conditions remains sixteen. The only determining factor is the number of real zero-multiplicity roots of each quartic. All of these results constitute proof of Claim 1 with the relevant results summarized in Theorem 1.

Theorem 1. 

Let $\Omega \subset {\mathbb{R}}^3$ be simply connected, $\mathbf{v}:\mathbf{r}\in \Omega \times t\in {\mathbb{R}}_{\ge 0}\to {\mathbb{R}}^3$, $P:\mathbf{r}\in \Omega \times t\in {\mathbb{R}}_{\ge 0}\to \mathbb{R}$, $\mathbf{f}:\mathbf{r}\in \Omega \times t\in {\mathbb{R}}_{\ge 0}\to {\mathbb{R}}^3$, and ${\mathbf{v}}_0\left(\mathbf{r}\right)=\mathbf{v}\left(\mathbf{r},0\right)$. Then, $\mathbf{v}\left(\mathbf{r},t\right)$ (100)–(108) and $P\left(\mathbf{r},t\right)$ (110) satisfy both (1) and (2) such that the number of unique bounded $\mathbf{v}\left(\mathbf{r},t\right)$ which (1) and (2) are capable of generating from the same initial and boundary condition set $\{{\mathbf{v}}_0\left(\mathbf{r}\right),{\mathbf{v}}_{\Omega }\left(\mathbf{r},t\right),P_{\Omega }\left(\mathbf{r},t\right)\}$ is at least one and at most sixteen, while the number of unique $P\left(\mathbf{r},t\right)$ equals the number of unique $\mathbf{v}\left(\mathbf{r},t\right)$.

9. The Navier-Stokes Existence and Smoothness Problem

The official statement regarding the open question of Navier-Stokes existence and smoothness as a Millennium Problem was issued by the Clay Mathematics Institute in 2000 [51] and asks for a smooth and globally defined solution to (1) and (2) for all unbounded space and for all time given zero external forces and initial velocity and pressure fields which decay sufficiently rapidly that they are rendered square integrable on all space. A choice is given between $\Omega ={\mathbb{R}}^3$ or $\Omega ={\mathbb{T}}^3$. Since formulating $\mathbf{v}\left(\mathbf{r},t\right)$ for unbounded $\Omega$ is best suited to ${\mathbb{R}}^3$ due to the added complications of periodic boundary conditions on ${\mathbb{T}}^3$, this work addresses the Clay problem in the topology of ${\mathbb{R}}^3$. The reason that the Clay problem prescribes $\mathbf{f}=\mathbf{0}$ lies in the physical meaning of the Clay problem when considered as an ocean of fluid infinitely wide and infinitely deep. Due to the Seeliger-Neumann paradox of an infinite amount of mass having infinite gravitational potential while canceling said potential due to perfectly even mass distribution [64], the focus is removed from the external force field and placed rightly on the velocity field by mandating that forces cancel. The prescribed space of initial velocity fields ${\mathbf{v}}_0\left(\mathbf{r}\right)$ are chosen to be those which are infinitely differentiable, compactly supported, and square integrable because such functional analytic properties are expected to encapsulate physically relevant encodings of an initial velocity field. The content of Problem A within the official Clay problem statement for $\Omega ={\mathbb{R}}^3$ may be distilled into Conditions 1-3 while the claim of this work as regards the Millennium Problem is outlined in Claim 4.

Condition 1. 

The velocity $\mathbf{v}\left(\mathbf{r},t\right)$ and pressure $P\left(\mathbf{r},t\right)$ must satisfy the Navier-Stokes equations (1) and (2) throughout three dimensions of unbounded space for all time. Formally, $\mathbf{v}:\mathbf{r}\in {\mathbb{R}}^3\times t\in {\mathbb{R}}_{\ge 0}\to {\mathbb{R}}^3\wedge P:\mathbf{r}\in {\mathbb{R}}^3\times t\in {\mathbb{R}}_{\ge 0}\to \mathbb{R}\Rightarrow {\displaystyle \frac{\partial \mathbf{v}}{\partial t}}+\left(\mathbf{v}·\nabla \right)\mathbf{v}=\nu {\nabla }^2\mathbf{v}-{\displaystyle \frac{\nabla P}{\rho }}\wedge \nabla ·\mathbf{v}=0$.

Condition 2. 

The velocity $\mathbf{v}\left(\mathbf{r},t\right)$ and pressure $P\left(\mathbf{r},t\right)$ must remain infinitely differentiable everywhere for all time. Formally, $\{\mathbf{v}\in C^{\infty }\left({\mathbb{R}}^3,(0,t]\right)\wedge P\in C^{\infty }\left({\mathbb{R}}^3,(0,t]\right)\}\forall \{t,{\mathbf{v}}_0\in C_c^{\infty }\left({\mathbb{R}}^3\right)\wedge \nabla ·{\mathbf{v}}_0=0\}$.

Condition 3. 

The velocity $\mathbf{v}\left(\mathbf{r},t\right)$ and pressure $P\left(\mathbf{r},t\right)$ must remain globally square integrable for all time. Formally, $\{\mathbf{v}\in L^2\left({\mathbb{R}}^3,(0,t]\right)\wedge P\in L^2\left({\mathbb{R}}^3,(0,t]\right)\}$$\forall \{t,{\mathbf{v}}_0\in L^2\left({\mathbb{R}}^3\right)\}$.

Claim 4. 

Let $\mathbf{v}:\mathbf{r}\in {\mathbb{R}}^3\times t\in {\mathbb{R}}_{\ge 0}\to {\mathbb{R}}^3$, $P:\mathbf{r}\in {\mathbb{R}}^3\times t\in {\mathbb{R}}_{\ge 0}\to \mathbb{R}$, and ${\mathbf{v}}_0\left(\mathbf{r}\right)=\mathbf{v}\left(\mathbf{r},0\right)$. Then, $\mathbf{v}\left(\mathbf{r},t\right)$ and $P\left(\mathbf{r},t\right)$ satisfy both (1) and (2) such that $\{\mathbf{v}\left(\mathbf{r},t\right)\in \{C^{\infty }\left({\mathbb{R}}^3,(0,t]\right)\cap L^2\left({\mathbb{R}}^3,(0,t]\right)\}\wedge P\left(\mathbf{r},t\right)\in \{C^{\infty }\left({\mathbb{R}}^3,(0,t]\right)\cap L^2\left({\mathbb{R}}^3,(0,t]\right)\}\forall \{t,{\mathbf{v}}_0\left(\mathbf{r}\right)\in \{C_c^{\infty }\left({\mathbb{R}}^3\right)\cap L^2\left({\mathbb{R}}^3\right)\}\wedge \nabla ·{\mathbf{v}}_0\left(\mathbf{r}\right)=0\}$.

9.1. Proofs of Global Existence of $\mathbf{v}\left(\mathbf{r},t\right)$ and $P\left(\mathbf{r},t\right)$ on $\Omega ={\mathbb{R}}^3$ for All $t\in {\mathbb{R}}_{\ge 0}$

The first step to proving Claim 4 is establishing existence and uniqueness of $\mathbf{v}\left(\mathbf{r},t\right)$ and $P\left(\mathbf{r},t\right)$ in compliance with Condition 1. The subtask of establishing existence and uniqueness must be completed before smoothness and square integrability are addressed because smoothness and square integrability presuppose both existence and some degree of uniqueness of the functions to be analyzed. The aim of this subtask is summarized in Claim 4.1.

Claim 4.1. 

Let $\mathbf{v}:\mathbf{r}\in {\mathbb{R}}^3\times t\in {\mathbb{R}}_{\ge 0}\to {\mathbb{R}}^3$, $P:\mathbf{r}\in {\mathbb{R}}^3\times t\in {\mathbb{R}}_{\ge 0}\to \mathbb{R}$, and ${\mathbf{v}}_0\left(\mathbf{r}\right)=\mathbf{v}\left(\mathbf{r},0\right)\in \{C_c^{\infty }\left({\mathbb{R}}^3\right)\cap L^2\left({\mathbb{R}}^3\right)\}$. Then, there exists a unique pair of $\mathbf{v}\left(\mathbf{r},t\right)$ and $P\left(\mathbf{r},t\right)$ satisfying both (1) and (2) everywhere on $\Omega ={\mathbb{R}}^3$ for all time $t\in {\mathbb{R}}_{\ge 0}$.

9.1.1. Proof of Global Existence of Velocity Field $\mathbf{v}\left(\mathbf{r},t\right)$ on $\Omega ={\mathbb{R}}^3$ for All t

The global existence and uniqueness of the velocity field is addressed first since $P\left(\mathbf{r},t\right)$ presupposes $\mathbf{v}\left(\mathbf{r},t\right)$ (109). As stated in Section 6.4, the unboundedness of $\Omega$ implies that all boundary conditions $\{{\mathbf{v}}_{\Omega },{\mathbf{v}}_{{\mathrm{irr}}_{\Omega }},{\mathbf{v}}_{{\mathrm{rot}}_{\Omega }},P_{\Omega }\}$ vanish with the further implication that $B_1{\widehat{x}}_1+B_2{\widehat{x}}_2+B_3{\widehat{x}}_3=\mathbf{0}$ everywhere. Modifying (30)–(32) to reflect $\left\{C_{k1}=0\wedge C_{k0}=1\right\}|k\in {\mathbb{N}}_{[1,3]}$ and therefore avoid a trivial $\mathbf{v}\left(\mathbf{r},t\right)$ as a consequence of $B_1{\widehat{x}}_1+B_2{\widehat{x}}_2+B_3{\widehat{x}}_3=\mathbf{0}$ reduces (30)–(32) to (111)–(113) respectively. The aim of handling the global existence and uniqueness of the velocity field is summarized in Claim 4.1.1.

Claim 4.1.1. 

Let $\mathbf{v}:\mathbf{r}\in {\mathbb{R}}^3\times t\in {\mathbb{R}}_{\ge 0}\to {\mathbb{R}}^3$ and ${\mathbf{v}}_0\left(\mathbf{r}\right)=\mathbf{v}\left(\mathbf{r},0\right)\in \{C_c^{\infty }\left({\mathbb{R}}^3\right)\cap L^2\left({\mathbb{R}}^3\right)\}$. Then, there exists a unique $\mathbf{v}\left(\mathbf{r},t\right)$ which satisfies both (1) and (2) everywhere on $\Omega ={\mathbb{R}}^3$ for all time $t\in {\mathbb{R}}_{\ge 0}$.

The Helmholtz decomposition (3) remains unique for unbounded $\Omega$ so long as velocity decays appropriately, which it must in physically sensible unbounded flows, and especially so for compactly supported ${\mathbf{v}}_0\left(\mathbf{r}\right)$. Since the irrotational velocity field depends on boundary conditions and inherits even its time dependency through boundary conditions, and all boundary conditions vanish for an unbounded domain, it is clear that $\mathbf{v}\left(\mathbf{r},t\right)={\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)$ for $\Omega ={\mathbb{R}}^3$. Applying $\mathbf{v}\left(\mathbf{r},t\right)={\mathbf{v}}_{\mathrm{rot}}\left(\mathbf{r},t\right)$ to (30)–(32) yields (111)–(113) respectively.

$$u\left(\mathbf{r},t\right)=u_{\mathrm{rot}}\left(\mathbf{r},t\right)={\displaystyle \frac{\nu {\left|\left|\nabla {\Phi }_1\left(\mathbf{r},t\right)\right|\right|}^2}{\left(\mathbf{v}\left(\mathbf{r},t\right)·\nabla \right){\Phi }_1\left(\mathbf{r},t\right)}}exp\left({\displaystyle \frac{\left(\left(\mathbf{v}\left(\mathbf{r},t\right)·\nabla \right){\Phi }_1\left(\mathbf{r},t\right)\right){\Phi }_1\left(\mathbf{r},t\right)}{\nu {\left|\left|\nabla {\Phi }_1\left(\mathbf{r},t\right)\right|\right|}^2}}\right){\widehat{x}}_1$$

(111)

$$v\left(\mathbf{r},t\right)=v_{\mathrm{rot}}\left(\mathbf{r},t\right)={\displaystyle \frac{\nu {\left|\left|\nabla {\Phi }_2\left(\mathbf{r},t\right)\right|\right|}^2}{\left(\mathbf{v}\left(\mathbf{r},t\right)·\nabla \right){\Phi }_2\left(\mathbf{r},t\right)}}exp\left({\displaystyle \frac{\left(\left(\mathbf{v}\left(\mathbf{r},t\right)·\nabla \right){\Phi }_2\left(\mathbf{r},t\right)\right){\Phi }_2\left(\mathbf{r},t\right)}{\nu {\left|\left|\nabla {\Phi }_2\left(\mathbf{r},t\right)\right|\right|}^2}}\right){\widehat{x}}_2$$

(112)

$$w\left(\mathbf{r},t\right)=w_{\mathrm{rot}}\left(\mathbf{r},t\right)={\displaystyle \frac{\nu {\left|\left|\nabla {\Phi }_3\left(\mathbf{r},t\right)\right|\right|}^2}{\left(\mathbf{v}\left(\mathbf{r},t\right)·\nabla \right){\Phi }_3\left(\mathbf{r},t\right)}}exp\left({\displaystyle \frac{\left(\left(\mathbf{v}\left(\mathbf{r},t\right)·\nabla \right){\Phi }_3\left(\mathbf{r},t\right)\right){\Phi }_3\left(\mathbf{r},t\right)}{\nu {\left|\left|\nabla {\Phi }_3\left(\mathbf{r},t\right)\right|\right|}^2}}\right){\widehat{x}}_3$$

(113)

Multiplying (111)–(113) by the reciprocals of their exponential prefactors and solving the resultant expressions for $\left(\mathbf{v}\left(\mathbf{r},t\right)·\nabla \right){\Phi }_k\left(\mathbf{r},t\right)$ via the principal branch $W_{n\in \{0,1\}}$ of the Lambert W function yields (114)–(116) respectively.

$$u\left(\mathbf{r},t\right){\nabla }_1{\Phi }_1\left(\mathbf{r},t\right)+v\left(\mathbf{r},t\right){\nabla }_2{\Phi }_1\left(\mathbf{r},t\right)+w\left(\mathbf{r},t\right){\nabla }_3{\Phi }_1\left(\mathbf{r},t\right)=$$

$${\displaystyle \frac{\nu {\left|\left|\nabla {\Phi }_1\left(\mathbf{r},t\right)\right|\right|}^2}{u\left(\mathbf{r},t\right)}}{\mathrm{W}}_n\left({\displaystyle \frac{{\Phi }_1\left(\mathbf{r},t\right)}{u\left(\mathbf{r},t\right)}}\right){\widehat{x}}_1$$

(114)

$$u\left(\mathbf{r},t\right){\nabla }_1{\Phi }_2\left(\mathbf{r},t\right)+v\left(\mathbf{r},t\right){\nabla }_2{\Phi }_2\left(\mathbf{r},t\right)+w\left(\mathbf{r},t\right){\nabla }_3{\Phi }_2\left(\mathbf{r},t\right)=$$

$${\displaystyle \frac{\nu {\left|\left|\nabla {\Phi }_2\left(\mathbf{r},t\right)\right|\right|}^2}{v\left(\mathbf{r},t\right)}}{\mathrm{W}}_n\left({\displaystyle \frac{{\Phi }_2\left(\mathbf{r},t\right)}{v\left(\mathbf{r},t\right)}}\right){\widehat{x}}_2$$

(115)

$$u\left(\mathbf{r},t\right){\nabla }_1{\Phi }_3\left(\mathbf{r},t\right)+v\left(\mathbf{r},t\right){\nabla }_2{\Phi }_3\left(\mathbf{r},t\right)+w\left(\mathbf{r},t\right){\nabla }_3{\Phi }_3\left(\mathbf{r},t\right)=$$

$${\displaystyle \frac{\nu {\left|\left|\nabla {\Phi }_3\left(\mathbf{r},t\right)\right|\right|}^2}{w\left(\mathbf{r},t\right)}}{\mathrm{W}}_n\left({\displaystyle \frac{{\Phi }_3\left(\mathbf{r},t\right)}{w\left(\mathbf{r},t\right)}}\right){\widehat{x}}_3$$

(116)

Introducing the invertible substitutions (117)–(119) places (114)–(116) into a form (120)–(122) similar to (111)–(113) yet distinct so as to assist in the next step regarding the handling of the heavy transcendental nonlinearities.

$$u\left(\mathbf{r},t\right)={\displaystyle \frac{a_5}{{\omega }_1\left(\mathbf{r},t\right)exp\left({\omega }_1\left(\mathbf{r},t\right)\right)}}\iff {\omega }_1\left(\mathbf{r},t\right)=W_0\left({\displaystyle \frac{a_5}{u\left(\mathbf{r},t\right)}}\right)$$

(117)

$$v\left(\mathbf{r},t\right)={\displaystyle \frac{a_{10}}{{\omega }_2\left(\mathbf{r},t\right)exp\left({\omega }_2\left(\mathbf{r},t\right)\right)}}\iff {\omega }_2\left(\mathbf{r},t\right)=W_0\left({\displaystyle \frac{a_{10}}{u\left(\mathbf{r},t\right)}}\right)$$

(118)

$$w\left(\mathbf{r},t\right)={\displaystyle \frac{a_{15}}{{\omega }_3\left(\mathbf{r},t\right)exp\left({\omega }_3\left(\mathbf{r},t\right)\right)}}\iff {\omega }_3\left(\mathbf{r},t\right)=W_0\left({\displaystyle \frac{a_{15}}{u\left(\mathbf{r},t\right)}}\right)$$

(119)

$${\displaystyle \frac{a_1{a}_5}{{\omega }_1\left(\mathbf{r},t\right)exp\left({\omega }_1\left(\mathbf{r},t\right)\right)}}+{\displaystyle \frac{a_2{a}_{10}}{{\omega }_2\left(\mathbf{r},t\right)exp\left({\omega }_2\left(\mathbf{r},t\right)\right)}}+{\displaystyle \frac{a_3{a}_{15}}{{\omega }_3\left(\mathbf{r},t\right)exp\left({\omega }_3\left(\mathbf{r},t\right)\right)}}=$$

$${\displaystyle \frac{a_4}{a_5}}{\omega }_1^2\left(\mathbf{r},t\right)exp\left({\omega }_1\left(\mathbf{r},t\right)\right)$$

(120)

$${\displaystyle \frac{a_6{a}_5}{{\omega }_1\left(\mathbf{r},t\right)exp\left({\omega }_1\left(\mathbf{r},t\right)\right)}}+{\displaystyle \frac{a_7{a}_{10}}{{\omega }_2\left(\mathbf{r},t\right)exp\left({\omega }_2\left(\mathbf{r},t\right)\right)}}+{\displaystyle \frac{a_8{a}_{15}}{{\omega }_3\left(\mathbf{r},t\right)exp\left({\omega }_3\left(\mathbf{r},t\right)\right)}}=$$

$${\displaystyle \frac{a_9}{a_{10}}}{\omega }_2^2\left(\mathbf{r},t\right)exp\left({\omega }_2\left(\mathbf{r},t\right)\right)$$

(121)

$${\displaystyle \frac{a_{11}{a}_5}{{\omega }_1\left(\mathbf{r},t\right)exp\left({\omega }_1\left(\mathbf{r},t\right)\right)}}+{\displaystyle \frac{a_{12}{a}_{10}}{{\omega }_2\left(\mathbf{r},t\right)exp\left({\omega }_2\left(\mathbf{r},t\right)\right)}}+{\displaystyle \frac{a_{13}{a}_{15}}{{\omega }_3\left(\mathbf{r},t\right)exp\left({\omega }_3\left(\mathbf{r},t\right)\right)}}=$$

$${\displaystyle \frac{a_{14}}{a_{15}}}{\omega }_3^2\left(\mathbf{r},t\right)exp\left({\omega }_3\left(\mathbf{r},t\right)\right)$$

(122)

Careful observation of patterns in the coupling points of (120)–(122) allow for their specification as a nonlinear matrix system of three coupled algebraic transcendental equations (123).

$$\left[\begin{array}{ccc}{a}_1& a_2& a_3\\ a_6& a_7& a_8\\ a_{11}& a_{12}& a_{13}\end{array}\right]\left[\begin{array}{c}{\displaystyle \frac{a_5}{{\omega }_1\left(\mathbf{r},t\right)exp\left({\omega }_1\left(\mathbf{r},t\right)\right)}}\\ {\displaystyle \frac{a_{10}}{{\omega }_2\left(\mathbf{r},t\right)exp\left({\omega }_2\left(\mathbf{r},t\right)\right)}}\\ {\displaystyle \frac{a_{15}}{{\omega }_3\left(\mathbf{r},t\right)exp\left({\omega }_3\left(\mathbf{r},t\right)\right)}}\end{array}\right]=\left[\begin{array}{c}{\displaystyle \frac{a_4}{a_5}}{\omega }_1\left(\mathbf{r},t\right)exp\left({\omega }_1\left(\mathbf{r},t\right)\right)\\ {\displaystyle \frac{a_9}{a_{10}}}{\omega }_2\left(\mathbf{r},t\right)exp\left({\omega }_2\left(\mathbf{r},t\right)\right)\\ {\displaystyle \frac{a_{14}}{a_{15}}}{\omega }_3\left(\mathbf{r},t\right)exp\left({\omega }_3\left(\mathbf{r},t\right)\right)\end{array}\right]$$

(123)

Defining symbols for the coefficient matrices (124) and unknown vectors (125) assist in writing (123) as a vector-valued form of an exponential-transcendental algebraic equation (126) which, if decoupled, would be capable of solution by the Lambert W function.

$$\mathbf{A}:=\left[\begin{array}{ccc}{a}_1& a_2& a_3\\ a_6& a_7& a_8\\ a_{11}& a_{12}& a_{13}\end{array}\right],\mathbf{B}:=\mathrm{diag}\left(a_5,a_{10},a_{15}\right),\mathbf{C}:=\mathrm{diag}\left({\displaystyle \frac{a_4}{a_5}},{\displaystyle \frac{a_9}{a_{10}}},{\displaystyle \frac{a_{14}}{a_{15}}}\right)$$

(124)

$$\kappa \left(\omega \left(\mathbf{r},t\right)\right)=\sum _{i=1}^3{\displaystyle \frac{exp\left(-{\omega }_i\left(\mathbf{r},t\right)\right)}{{\omega }_i\left(\mathbf{r},t\right)}},\eta \left(\omega \left(\mathbf{r},t\right)\right)=\sum _{i=1}^3{\omega }_i^2\left(\mathbf{r},t\right)exp\left({\omega }_i\left(\mathbf{r},t\right)\right){\widehat{x}}_i$$

(125)

$$\mathbf{A}\mathbf{B}\kappa \left(\omega \left(\mathbf{r},t\right)\right)=\mathbf{C}\eta \left(\omega \left(\mathbf{r},t\right)\right)$$

(126)

Transcendental problems similar to (126) have been addressed in the literature. Because such problems are canonically of the type of nonlinearity which has been considered tractable by the Lambert W function while their multivariate coupling hampers such effort, there has been research into multivariate versions of the Lambert W function [65]. While not many closed forms of multivariate Lambert W functions yet exist, the existence and uniqueness of solutions in the abstract may still be reliably proven. It is possible that explicit forms of multivariate Lambert W functions require highly general hypergeometric functions such as the Weierstrass elliptic function or the Meijer G function. For the purposes of addressing existence and uniqueness of $\mathbf{v}$, the inverse of (126) is considered as the trivariate Lambert W function $W_{\left(3\right)}$ (127) whose subscript k denotes the vector component ${\omega }_k$ of its value. The velocity field $\mathbf{v}\left(\mathbf{r},t\right)$ is then recovered (128) by substituting ${\omega }_k$ into (117)–(119).

$$\omega \left(\mathbf{r},t\right)=F^{-1}\left(\mathbf{A},\mathbf{B},\mathbf{C}\right)|F\left(\omega \left(\mathbf{r},t\right)\right)=\mathbf{A}\mathbf{B}\kappa \left(\omega \left(\mathbf{r},t\right)\right)=\mathbf{C}\eta \left(\omega \left(\mathbf{r},t\right)\right)=0$$

(127)

$$\mathbf{v}\left(\mathbf{r},t\right)=\sum _{k=1}^3\left[{\displaystyle \frac{a_{5k}exp\left(-W_{\left(3\right)}{\left(\mathbf{A},\mathbf{B},\mathbf{C}\right)}_k\right)}{W_{\left(3\right)}{\left(\mathbf{A},\mathbf{B},\mathbf{C}\right)}_k}}\right]{\widehat{x}}_k$$

(128)

To prove existence and uniqueness of the multivariate Lambert W-function in order to ensure its establishment as a valid component of an explicit expression of an unbounded velocity field (128), the method of linear homotopy is used [67]. The method of linear homotopy works by first establishing existence (and uniqueness if necessary) of a known subcase of the problem to be addressed, and then the existence of a sufficiently smooth mapping between the known subcase and the fully general solution is proven, which is designed to preserve essential functional analytic properties such as uniqueness so that the general solution may also have its existence and uniqueness proven. An edge case of (114)–(116) solvable in closed form is that which corresponds to $\mathbf{A}=\mathrm{diag}\left(\mathbf{A}\right)$, that is, when all nondiagonal components of $\mathbf{A}$ vanish (129)–(131), enabling the Lambert W function to be used in recovering the MCH lifts from the heat equation MCH kernels to a reduced version of the velocity field ${\mathbf{v}}_{\mathrm{loc}}\left(\mathbf{r},t\right)$. This decoupled case of (114)–(116) is used as the initial condition of the homotopic mapping $\mathcal{H}:T\in {\mathbb{R}}_{[0,1]}\to \mathbb{R}$ where $\mathcal{H}\left(0\right)={\mathbf{v}}_{\mathrm{loc}}\left(\mathbf{r},t\right)$ and $\mathcal{H}\left(1\right)=\mathbf{v}\left(\mathbf{r},t\right)$. The solution to (114)–(116) is given by (132) while its equivalent ${\omega }_{\mathrm{loc}}\left(\mathbf{r},t\right)$ with $\omega \left(\mathbf{r},t\right)$ is given by (133).

$${\displaystyle \frac{{\nabla }_1{\Phi }_1\left(\mathbf{r},t\right)}{\nu {\left|\left|\nabla {\Phi }_1\left(\mathbf{r},t\right)\right|\right|}^2}}{u}_{\mathrm{rot}}^2\left(\mathbf{r},t\right)={\mathrm{W}}_0\left({\displaystyle \frac{{\Phi }_1\left(\mathbf{r},t\right)}{u_{\mathrm{rot}}\left(\mathbf{r},t\right)}}\right){\widehat{x}}_1$$

(129)

$${\displaystyle \frac{{\nabla }_2{\Phi }_2\left(\mathbf{r},t\right)}{\nu {\left|\left|\nabla {\Phi }_2\left(\mathbf{r},t\right)\right|\right|}^2}}{v}_{\mathrm{rot}}^2\left(\mathbf{r},t\right)={\mathrm{W}}_0\left({\displaystyle \frac{{\Phi }_2\left(\mathbf{r},t\right)}{v_{\mathrm{rot}}\left(\mathbf{r},t\right)}}\right){\widehat{x}}_2$$

(130)

$${\displaystyle \frac{{\nabla }_3{\Phi }_3\left(\mathbf{r},t\right)}{\nu {\left|\left|\nabla {\Phi }_3\left(\mathbf{r},t\right)\right|\right|}^2}}{w}_{\mathrm{rot}}^2\left(\mathbf{r},t\right)={\mathrm{W}}_0\left({\displaystyle \frac{{\Phi }_3\left(\mathbf{r},t\right)}{w_{\mathrm{rot}}\left(\mathbf{r},t\right)}}\right){\widehat{x}}_3$$

(131)

$${\mathbf{v}}_{\mathrm{loc}}\left(\mathbf{r},t\right)=\sum _{i=1}^3\left[\pm \sqrt{{\displaystyle \frac{3}{2}}\sqrt{{\displaystyle \frac{{\nabla }_i{\Phi }_i\left(\mathbf{r},t\right)}{\nu {\left|\left|\nabla {\Phi }_i\left(\mathbf{r},t\right)\right|\right|}^2}}}{\mathrm{W}}_0\left({\displaystyle \frac{2}{3}}\sqrt[3]{{\displaystyle \frac{{\Phi }_i^2\left(\mathbf{r},t\right){\nabla }_i{\Phi }_i\left(\mathbf{r},t\right)}{\nu {\left|\left|\nabla {\Phi }_i\left(\mathbf{r},t\right)\right|\right|}^2}}}\right)}{\widehat{x}}_i\right]$$

(132)

$${\omega }_{\mathrm{loc}}\left(\mathbf{r},t\right)=\sum _{i=1}^3\left[{\displaystyle \frac{a_{5i}}{{\mathrm{W}}_0{\left(\pm \sqrt{\frac{3}{2}\sqrt{\frac{{\nabla }_i{\Phi }_i\left(\mathbf{r},t\right)}{\nu {\left|\left|\nabla {\Phi }_i\left(\mathbf{r},t\right)\right|\right|}^2}}{\mathrm{W}}_0\left(\frac{2}{3}\sqrt[3]{\frac{{\Phi }_i^2\left(\mathbf{r},t\right){\nabla }_i{\Phi }_i\left(\mathbf{r},t\right)}{\nu {\left|\left|\nabla {\Phi }_i\left(\mathbf{r},t\right)\right|\right|}^2}}\right)}\right)}^{-1}}}\right]{\widehat{x}}_i$$

(133)

The linear homotopy mapping is defined (134) so as to render $\mathcal{H}\left(0\right)$ as representing ${\omega }_{\mathrm{loc}}\left(\mathbf{r},t\right)$ (133) and $\mathcal{H}\left(1\right)$ as representing $\omega \left(\mathbf{r},t\right)$ (127).

$$\mathcal{H}\left(T\right)=\mathrm{diag}\left(a_1,a_7,a_{13}\right)+T\left(\mathbf{A}-\mathrm{diag}\left(a_1,a_7,a_{13}\right)\right)$$

(134)

The Jacobian matrix $D_{\mathbf{v}\left(\mathbf{r},t\right)}F(\mathbf{v}\left(\mathbf{r},t\right),T)$ of (127) is taken with respect to $\mathbf{v}\left(\mathbf{r},t\right)$ (135), taking care to include (134) in its formulation.

$$D_{\omega \left(\mathbf{r},t\right)}F(\omega \left(\mathbf{r},t\right),T)=\mathbf{C}\mathrm{diag}\left({\displaystyle \frac{\partial \left({\omega }_1^2\left(\mathbf{r},t\right)exp\left({\omega }_1\left(\mathbf{r},t\right)\right)\right)}{\partial {\omega }_1\left(\mathbf{r},t\right)}},{\displaystyle \frac{\partial \left({\omega }_2^2\left(\mathbf{r},t\right)exp\left({\omega }_2\left(\mathbf{r},t\right)\right)\right)}{\partial {\omega }_2\left(\mathbf{r},t\right)}},\right.$$

$$\left.{\displaystyle \frac{\partial \left({\omega }_3^2\left(\mathbf{r},t\right)exp\left({\omega }_3\left(\mathbf{r},t\right)\right)\right)}{\partial {\omega }_3\left(\mathbf{r},t\right)}}\right)-\mathcal{H}\left(T\right)\mathbf{B}\mathrm{diag}\left({\displaystyle \frac{\partial \left({({\omega }_1^2\left(\mathbf{r},t\right)exp\left({\omega }_1\left(\mathbf{r},t\right)\right))}^{-1}\right)}{\partial {\omega }_1\left(\mathbf{r},t\right)}},\right.$$

$$\left.{\displaystyle \frac{\partial \left({({\omega }_2^2\left(\mathbf{r},t\right)exp\left({\omega }_2\left(\mathbf{r},t\right)\right))}^{-1}\right)}{\partial {\mathbf{v}}_2\left(\mathbf{r},t\right)}},{\displaystyle \frac{\partial \left({({\omega }_3^2\left(\mathbf{r},t\right)exp\left({\omega }_3\left(\mathbf{r},t\right)\right))}^{-1}\right)}{\partial {\omega }_3\left(\mathbf{r},t\right)}}\right)$$

(135)

Evaluating the first partial derivatives in the diagonal matrix multiplying $\mathbf{C}$ (136) in (135) and also in the diagonal matrix multiplying $\mathbf{B}$ (137) in (135) reveal a natural exponential constraint that ${\omega }_i\left(\mathbf{r},t\right)>0$. Substituting these evaluations of the system Jacobian matrix entries into the system Jacobian matrix itself (135) yields constraints on the diagonal and off-diagonal terms, specifically that the diagonal terms must be both positive and nonzero everywhere (138), while the off-diagonal terms (139) need only be nonnegative everywhere.

$${\displaystyle \frac{\partial \left({\omega }_i^2\left(\mathbf{r},t\right)exp\left({\omega }_i\left(\mathbf{r},t\right)\right)\right)}{\partial {\omega }_i\left(\mathbf{r},t\right)}}={\omega }_i\left(\mathbf{r},t\right)exp\left({\omega }_i\left(\mathbf{r},t\right)\right)\left({\omega }_i\left(\mathbf{r},t\right)+2\right)>0\Rightarrow {\omega }_i\left(\mathbf{r},t\right)>0$$

(136)

$${\displaystyle \frac{\partial \left({\left({\omega }_i^2\left(\mathbf{r},t\right)exp\left({\omega }_i\left(\mathbf{r},t\right)\right)\right)}^{-1}\right)}{\partial {\omega }_i\left(\mathbf{r},t\right)}}={\displaystyle \frac{-1-{\omega }_i\left(\mathbf{r},t\right)}{{\omega }_i^2\left(\mathbf{r},t\right)exp\left({\omega }_i\left(\mathbf{r},t\right)\right)}}<0\Rightarrow {\omega }_i\left(\mathbf{r},t\right)>-1$$

(137)

$${\left[D_{\omega \left(\mathbf{r},t\right)}F(\omega \left(\mathbf{r},t\right),T)\right]}_{ii}=C_i{\displaystyle \frac{\partial \left({\omega }_i^2\left(\mathbf{r},t\right)exp\left({\omega }_i\left(\mathbf{r},t\right)\right)\right)}{\partial {\omega }_i\left(\mathbf{r},t\right)}}-a_{ii}{b}_i{\displaystyle \frac{\partial \left({\left({\omega }_i^2\left(\mathbf{r},t\right)exp\left({\omega }_i\left(\mathbf{r},t\right)\right)\right)}^{-1}\right)}{\partial {\omega }_i\left(\mathbf{r},t\right)}}$$

$$=C_i{\omega }_i\left(\mathbf{r},t\right)exp\left({\omega }_i\left(\mathbf{r},t\right)\right)\left({\omega }_i\left(\mathbf{r},t\right)+2\right)+a_{ii}{b}_i{\displaystyle \frac{1+{\omega }_i\left(\mathbf{r},t\right)}{{\omega }_i^2\left(\mathbf{r},t\right)exp\left({\omega }_i\left(\mathbf{r},t\right)\right)}}>0$$

(138)

$${\left[D_{\omega \left(\mathbf{r},t\right)}F(\omega ,T)\right]}_{ij}=-a_{ij}{b}_i{\displaystyle \frac{\partial \left({\left({\omega }_i^2\left(\mathbf{r},t\right)exp\left({\omega }_i\left(\mathbf{r},t\right)\right)\right)}^{-1}\right)}{\partial {\omega }_i\left(\mathbf{r},t\right)}}=$$

$$a_{ij}{b}_i{\displaystyle \frac{1+{\omega }_i\left(\mathbf{r},t\right)}{{\omega }_i^2\left(\mathbf{r},t\right)exp\left({\omega }_i\left(\mathbf{r},t\right)\right)}}\ge 0$$

(139)

Invoking the real-analytic Implicit Function Theorem [68] regarding the nonzero partial derivative terms in the system Jacobian (135) reveals a unique real-analytic mapping $\omega \left(T\right)$ (140) such that $\omega \left(T\right)>0$ for small T. This mapping may be extended over arbitrary T.

$$\omega \left(T\right)|F\left(\omega \left(T\right),T\right)=0\forall \{T\in [0,ϵ],ϵ\in {\mathbb{R}}_{>0}\}$$

(140)

Let $0<\underline{\omega }\left(\mathbf{r},t\right)\le \omega \left(\mathbf{r},t\right)\le \overline{\omega }\left(\mathbf{r},t\right)$ and define a priori$\mathcal{B}={[\underline{\omega }\left(\mathbf{r},t\right),\overline{\omega }\left(\mathbf{r},t\right)]}^3$. Then, the mappings $\mathfrak{A}\left(\omega \left(\mathbf{r},t\right)\right)$ (141) and $\mathfrak{B}\left(\omega \left(\mathbf{r},t\right)\right)$ (142) are continuous on $\mathcal{B}$ with minimum ${\mathfrak{A}}^{*}$ (143) and maximum ${\mathfrak{B}}_{*}$ (144) on $\mathcal{B}$ provided that $\{0<\mathcal{H}(T)<\infty \}\forall T$. Reasonable choices of $\underline{\omega }\left(\mathbf{r},t\right)$ and $\overline{\omega }\left(\mathbf{r},t\right)$ which actualize this positive finite homotopy mapping constraint are (145) and (146) respectively, where ${\omega }_i\left(\mathbf{r},t\right)$ is the ith component of $\omega \left(\mathbf{r},t\right)$.

$$\mathfrak{A}\left(\omega \left(\mathbf{r},t\right)\right):=\underset{i}{min}\left({\displaystyle \frac{C_{6i-5}{\omega }_i\left(\mathbf{r},t\right)\left({\omega }_i\left(\mathbf{r},t\right)+2\right)}{exp\left({\omega }_i\left(\mathbf{r},t\right)\right)}}\right)$$

(141)

$$\mathfrak{B}\left(\omega \left(\mathbf{r},t\right)\right):=\underset{i}{max}\left(\sum _{i\ne j}{a}_{ij}{b}_j\left|{\displaystyle \frac{-1-{\omega }_i\left(\mathbf{r},t\right)}{{\omega }_i^2\left(\mathbf{r},t\right)exp\left({\omega }_i\left(\mathbf{r},t\right)\right)}}\right|\right)$$

(142)

$${\mathfrak{A}}_{*}:=\underset{\omega \left(\mathbf{r},t\right)\in \mathcal{B}}{min}\left(\mathfrak{A}\left(\omega \left(\mathbf{r},t\right)\right)\right)>0$$

(143)

$${\mathfrak{B}}^{*}:=\underset{\omega \left(\mathbf{r},t\right)\in \mathcal{B}}{max}\left(\mathfrak{B}\left(\omega \left(\mathbf{r},t\right)\right)\right)<\infty$$

(144)

$$\underline{\omega }\left(\mathbf{r},t\right)={\displaystyle \frac{1}{2}}\underset{i}{min}\left({\omega }_{i_{\mathrm{loc}}}\left(\mathbf{r},t\right)\right)$$

(145)

$$\overline{\omega }\left(\mathbf{r},t\right)=2\underset{i}{max}\left({\omega }_{i_{\mathrm{loc}}}\left(\mathbf{r},t\right)\right)$$

(146)

Comparing (147) the bounds of both the diagonal and off-diagonal system Jacobian terms (138) and (139) with the minimum ${\mathfrak{A}}^{*}$ (143) and maximum ${\mathfrak{B}}_{*}$ (144) of the mappings (141) and (142) respectively on $\mathcal{B}$ while observing that the identity of $\mathcal{H}\left(T\right)$ as a convex combination of both itself and $\mathrm{diag}(a_1,a_7,a_{13})$ (134) implies that the dominant entries of $\mathbf{A}$ are its row diagonals (148), effectively ensuring Gershgorin nonsingularity of the system Jacobian (135). The row diagonal dominance of $\mathcal{H}\left(T\right)$ implies a similar dominance property for the system Jacobian (135) on $\mathcal{B}$. Since Gershgorin nonsingular matrices possess nonnegative inverses, and both the lower and upper bounds (143) and (144) of $\mathcal{B}$ are arbitrary in their magnitude, (140) is globally injective on $\mathcal{B}$.

$$\left\{{\left[D_{\omega \left(\mathbf{r},t\right)}F(\omega \left(\mathbf{r},t\right),T)\right]}_{ii}\ge {\mathfrak{A}}_{*}\wedge \sum _{i\ne j}{\left[D_{\omega \left(\mathbf{r},t\right)}F(\omega \left(\mathbf{r},t\right),T)\right]}_{ij}\le {\mathfrak{B}}^{*}\right\}\forall T\in [0,1]$$

(147)

$${\mathfrak{A}}_{*}>{\mathfrak{B}}^{*}$$

(148)

If $\mathcal{T}$ is defined by (149), then, $\mathcal{T}$ is a nonempty set because (140) implies $[0,ϵ]\subset \mathcal{T}$. If ${\tau }_0\in \mathcal{T}$, then (148) predicts a nonzero system Jacobian determinant on the branch $T=T_0$ which would be analytically extended via (140), implying an open $\mathcal{T}$. But if by the nonemptiness of $\mathcal{T}$ there exists some $T=T_n\in \mathcal{T}$ and $T_n\to T_{*}$ for some $T_{*}\in \mathcal{T}$, then $\omega \left(\mathbf{r},t\right)$ corresponding to $T_n$ lies in compact $\mathcal{B}$ and therefore admits a convergent subsequence $\omega \left(\mathbf{r},t\right):T=T_n\to \omega \left(\mathbf{r},t\right):T=T_{*}$. Combining the facts that the continuity of (140) yields $F(\omega \left(\mathbf{r},t\right):T=T_{*},T=T_{*})=0$ and that the system Jacobian is Gershgorin nonsingular implies a unique $\mathbf{v}\left(\mathbf{r},t\right):T=T_{*}$ on $\mathcal{B}$ [66]. The whole space of $\mathcal{T}$ is covered because of the arbitrary choice of $T_n$ and $T_{*}$. Therefore, $\mathcal{T}$ is a closed set, and because the only nonempty set that is simultaneously open and closed is a set with exactly one element, there exists exactly one $\omega \left(\mathbf{r},t\right)$ for all $0\le T\le 1$. Because $\omega \left(\mathbf{r},t\right)$ is unique and mapped to (128) by (129)–(131), the full $\mathbf{v}\left(\mathbf{r},t\right)$ (128) at $\mathcal{H}\left(1\right)$ and general $\mathbf{A}$ is also unique.

$$\mathcal{T}:=\{T\in [0,1]:\exists !{\mathbf{v}}^{\mathrm{T}}\left(\mathbf{r},t\right)\in \mathcal{B}\}$$

(149)

The final analysis is that there exists a unique linear homotopy $\mathcal{H}\left(T\right)$ between ${\mathbf{v}}_{\mathrm{loc}}\left(\mathbf{r},t\right)$ using the standard principal branch of the Lambert W function and $\mathbf{v}\left(\mathbf{r},t\right)$ using the trivariate Lambert W function $W_{\left(3\right)}$, proving existence and uniqueness of the full velocity field $\mathbf{v}\left(\mathbf{r},t\right)$ (128) satisfying (1) and (2) on ${\mathbb{R}}^3$ in a multivariate Lambert W functional context. Both the general unbounded MCH lift system (114)–(116) and its decoupled, locally solvable case ${\mathbf{v}}_{\mathrm{loc}}\left(\mathbf{r},t\right)$ (132) from which the $\mathbf{v}\left(\mathbf{r},t\right)$ was homotopized now possess known analytical solutions. The next step in proving global existence and uniqueness of $\mathbf{v}\left(\mathbf{r},t\right)$ is to investigate the behavior of the MCH kernels and reconciliations. The fact that the ${\Phi }_k\left(\mathbf{r},t\right)$ satisfy unbounded heat equations does not itself imply that ${\Phi }_k\left(\mathbf{r},t\right)\ge 0$ and $\left|\right|{\nabla }_k{\Phi }_k\left(\mathbf{r},t\right)\left|\right|>0$, both of which must be met in order to rule out blowup in ${\mathbf{v}}_{\mathrm{loc}}\left(\mathbf{r},t\right)$ (132) and therefore in $\mathbf{v}\left(\mathbf{r},t\right)$ (128). For said conditions to be met, the initial condition ${\Phi }_{k_0}$ of the heat equation must be nonnegative everywhere. Solving the MCH kernel heat equations (25) on $\Omega ={\mathbb{R}}^3$ which implies that ${\Phi }_{{\Omega }_k}(\mathbf{r},t)=0$ everywhere for $k\in \{1,2,3\}$. The heat equation solutions (29) reduce to their unbounded case (150) which is well understood [69]. Even when ${\Phi }_{k_0}\left(\mathbf{r}\right)>0$ everywhere, global existence and uniqueness of a solution to the heat equation on ${\mathbb{R}}^3$ is guaranteed so long as ${\Phi }_{k_0}\left(\mathbf{r}\right)$ has subquadratic exponential growth $\left|\right|{\Phi }_{k_0}\left|\right|\le Cexp\left({\lambda \left|\right|\mathbf{r}\left|\right|}^2\right)\forall \{C\in \mathbb{R},\lambda \in {\mathbb{R}}_{>0}\}$ at infinity [70]. The combined conditions of positivity and subquadratic growth yield a global existence and uniqueness condition (151) for ${\Phi }_k\left(\mathbf{r},t\right)$ insofar as said existence depends on ${\Phi }_{k_0}\left(\mathbf{r},t\right)$. The condition (151) is a very loose condition since the initial velocity field in the Millennium Problem is posed not as growing, but as decaying to zero as $\left|\mathbf{r}\right|\to \infty$.

$${\displaystyle \frac{\partial {\Phi }_k\left(\mathbf{r},t\right)}{\partial t}}=\nu {\nabla }^2{\Phi }_k\left(\mathbf{r},t\right)\Rightarrow {\Phi }_k\left(\mathbf{r},t\right)={\int }_{{\mathbb{R}}^3}{\displaystyle \frac{{\Phi }_{k_0}\left({\mathbf{r}}^{\prime }\right)}{8{\left(\pi \nu t\right)}^{3/2}}}exp\left(-{\displaystyle \frac{{\left|\left|\mathbf{r}-{\mathbf{r}}^{\prime }\right|\right|}^2}{4\nu t}}\right)d{V}_{{\mathbf{r}}^{\prime }}$$

(150)

$$0\le \left|\right|{\Phi }_{k_0}\left(\mathbf{r}\right)\left|\right|\le Cexp\left({\lambda \left|\right|\mathbf{r}\left|\right|}^2\right)\forall \{C\in \mathbb{R},\lambda \in {\mathbb{R}}_{>0}\}\mathrm{as}\mathbf{r}\to \infty$$

(151)

The MCH reconciliation formulae (106) and (107) adapting the heat equation initial condition ${\Phi }_{k_0}\left(\mathbf{r}\right)$ to the velocity field initial condition ${\mathbf{v}}_0\left(\mathbf{r}\right)$ reduce to transport invariance in the ${\Phi }_k\left(\mathbf{r}\right)$-oriented recasting of the algebraic MCH lift system into the first-order PDE system (66)–(68) for the MCH reconciliation when $\Omega ={\mathbb{R}}^3\Rightarrow B_k\left(\mathbf{r},t\right)=0|k\in \{1,2,3\}$ (152). Solving (152) by the method of characteristics produces a system similar in form to (72) and (73) prescribing a flow map $V_{t^{\prime }}\left(\mathbf{r}\right)$ encoding characteristics $V_{t^{\prime }}\left(\mathbf{r}\right)$ of ${\mathbf{v}}_0\left(\mathbf{r}\right)$ (153) and a connection to ${\Phi }_{k_0}\left(\mathbf{r}\right)$ via the finite and nonnegative Césaro average orbit integral (154) for arbitrary $h\in C^{\infty }$ which is square integrable so long as ${\mathbf{v}}_0\left(\mathbf{r}\right)$ is. Because ${\mathbf{v}}_0\left(\mathbf{r}\right)\in L^2$, choosing $h=\left|\right|{\mathbf{v}}_0\left(\mathbf{r}\right){\left|\right|}^2$ ensures that (154) becomes physically and analytically meaningful (155) while using the same ${\mathbf{v}}_0\left(\mathbf{r}\right)\in \{C_c^{\infty }\left({\mathbb{R}}^3\right)\cap L^2\left({\mathbb{R}}^3\right)\}$.

$$\left({\mathbf{v}}_0\left(\mathbf{r}\right)·\nabla \right){\Phi }_{k_0}\left(\mathbf{r}\right)=0$$

(152)

$$\left\{{\displaystyle \frac{d{V}_{t^{\prime }}\left(\mathbf{r}\right)}{d{t}^{\prime }}}={\mathbf{v}}_0\left(V_{t^{\prime }}\left(\mathbf{r}\right)\right)\wedge V_0\left(\mathbf{r}\right)=\mathbf{r}\right\}\Rightarrow t^{\prime }={\int }_{\mathbf{r}}^{V_{t^{\prime }}\left(\mathbf{r}\right)}{\displaystyle \frac{d\mathbf{s}}{{\mathbf{v}}_0\left(\mathbf{s}\right)}}$$

(153)

$${\Phi }_{k_0}\left(\mathbf{r}\right)=\underset{\Gamma \to 0}{lim}\left({\int }_{-\Gamma }^{\Gamma }{\displaystyle \frac{h\left(V_{t^{\prime }}\left(\mathbf{r}\right)\right)d{t}^{\prime }}{\Gamma }}\right)$$

(154)

$${\Phi }_{k_0}\left(\mathbf{r}\right)=\underset{\Gamma \to 0}{lim}\left({\int }_{-\Gamma }^{\Gamma }{\displaystyle \frac{\left|\right|{\mathbf{v}}_0\left(V_{t^{\prime }}\left(\mathbf{r}\right)\right){\left|\right|}^{2}d{t}^{\prime }}{\Gamma }}\right)$$

(155)

Since ${\mathbf{v}}_0\left(\mathbf{r}\right)$ is compactly supported, ${\Phi }_{k_0}\left(\mathbf{r}\right)$ becomes strictly positive everywhere while also lying in $\{C_c^{\infty }\left({\mathbb{R}}^3\right)\cap L^2\left({\mathbb{R}}^3\right)\}$. This strict positivity meets the condition (151) guaranteeing global regularity of ${\Phi }_k\left(\mathbf{r},t\right)$ for the required ${\mathbf{v}}_0\left(\mathbf{r}\right)\in L^2\left({\mathbb{R}}^3\right)$. The heat equation solutions ${\Phi }_k\left(\mathbf{r},t\right)$ (156) for unbounded flow on $\Omega ={\mathbb{R}}^3$ and arbitrary ${\mathbf{v}}_0\left(\mathbf{r}\right)\in \{C_c^{\infty }\left({\mathbb{R}}^3\right)\cap L^2\left({\mathbb{R}}^3\right)\}$ are given (156) by substituting (155) into (150).

$${\Phi }_k\left(\mathbf{r},t\right)={\int }_{{\mathbb{R}}^3}{\displaystyle \frac{{lim}_{\Gamma \to 0}\left({\int }_{-\Gamma }^{\Gamma }\frac{\left|\right|{\mathbf{v}}_0\left(V_{t^{\prime }}\left(\mathbf{r}\right)\right){\left|\right|}^{2}d{t}^{\prime }}{\Gamma }\right)}{8{\left(\pi \nu t\right)}^{3/2}}}exp\left(-{\displaystyle \frac{{\left|\left|\mathbf{r}-{\mathbf{r}}^{\prime }\right|\right|}^2}{4\nu t}}\right)d{V}_{{\mathbf{r}}^{\prime }},t^{\prime }={\int }_{\mathbf{r}}^{V_{t^{\prime }}\left(\mathbf{r}\right)}{\displaystyle \frac{d\mathbf{s}}{{\mathbf{v}}_0\left(\mathbf{s}\right)}}$$

(156)

Because any prescribed ${\mathbf{v}}_0\left(\mathbf{r}\right)\in \{C_c^{\infty }\left({\mathbb{R}}^3\right)\cap L^2\left({\mathbb{R}}^3\right)\}$ generates a set of globally defined ${\Phi }_k\left(\mathbf{r},t\right)>0\forall \{\left(\mathbf{r},t\right)\in \left({\mathbb{R}}^3\times t\in {\mathbb{R}}_{\ge 0}\right)\}$ (156) compliant with (151) so as to generate a unique $\mathbf{v}\left(\mathbf{r},t\right)$ (128) which also exists globally, the velocity field (128) meets Condition 1. The result is stated in Lemma 2.1.

Lemma 2.1. 

Let $\mathbf{v}:\left(\mathbf{r}\in {\mathbb{R}}^3\times t\in {\mathbb{R}}_{\ge 0}\right)\to {\mathbb{R}}^3$ and $\mathbf{v}(\mathbf{r},0)={\mathbf{v}}_0\left(\mathbf{r}\right)\in \{C_c^{\infty }\left({\mathbb{R}}^3\right)\cap L^2\left({\mathbb{R}}^3\right)\}$. Then, there exists a unique velocity field $\mathbf{v}(\mathbf{r},t)$ (128) satisfying (1) and (2) everywhere on $\Omega ={\mathbb{R}}^3$ for all time $0\le t<\infty$ in compliance with Condition 1.

9.1.2. Proof of Global Existence of Pressure Field $P\left(\mathbf{r},t\right)$ on $\Omega ={\mathbb{R}}^3$ for All t

Claim 4.1.2. 

Let $P:\left(\mathbf{r}\in {\mathbb{R}}^3\times t\in {\mathbb{R}}_{\ge 0}\right)\to \mathbb{R}$. Then, there exists a unique pressure field $P(\mathbf{r},t)$ (157) satisfying (1) and (2) everywhere on $\Omega ={\mathbb{R}}^3$ for all time $0\le t<\infty$ given a corresponding unique velocity field $\mathbf{v}\left(\mathbf{r},t\right)$ (Lemma 2.1) which satisfies (1) and (2) together with $P\left(\mathbf{r},t\right)$ in compliance with Condition 1.

The pressure-Poisson formulation (109) and solution (110) is readily extensible to an unbounded case by setting $\Omega ={\mathbb{R}}^3$ [63], implying that $P_{\Omega }(\mathbf{r},t)=0$ everywhere, and setting $\mathbf{f}(\mathbf{r},t)=\mathbf{0}$ everywhere (157). Since $\mathbf{v}\left(\mathbf{r},t\right)$ both exists and is unique everywhere, and $P(\mathbf{r},t)$ depends on $\mathbf{v}\left(\mathbf{r},t\right)$, it follows from Lemma 2.1 that $P(\mathbf{r},t)$ exists and is unique everywhere (Lemma 2.2).

$$P\left(\mathbf{r},t\right)=-{\displaystyle \frac{\rho }{4\pi }}{\int }_{{\mathbb{R}}^3}{\displaystyle \frac{\nabla ·\frac{D\mathbf{v}\left({\mathbf{r}}^{\prime },t\right)}{Dt}}{\left|\left|\mathbf{r}-{\mathbf{r}}^{\prime }\right|\right|}}d{V}_{{\mathbf{r}}^{\prime }}$$

(157)

Lemma 2.2. 

Let $P:\left(\mathbf{r}\in {\mathbb{R}}^3\times t\in {\mathbb{R}}_{\ge 0}\right)\to \mathbb{R}$. Then, there exists a unique pressure field $P(\mathbf{r},t)$ (157) satisfying (1) and (2) everywhere on $\Omega ={\mathbb{R}}^3$ for all time $0\le t<\infty$ whose uniqueness is implied by that of its corresponding velocity field $\mathbf{v}\left(\mathbf{r},t\right)$ (Lemma 2.1) which satisfies (1) and (2) together with $P(\mathbf{r},t)$ in compliance with Condition 1.

9.1.3. Conclusion of Existence

Since there always exists at least one unique pair of a velocity field (128) and a corresponding pressure field (157) together satisfying (1) and (2) on $\Omega ={\mathbb{R}}^3$ given both $\mathbf{f}\left(\mathbf{r},t\right)=\mathbf{0}$ and ${\mathbf{v}}_0\left(\mathbf{r}\right)\in L^2\left({\mathbb{R}}^3\right)$ and thus obeying Condition 1, Claim 4.1 stands with its result summarized in Lemma 3.

Lemma 3. 

Let $\mathbf{v}:\left(\mathbf{r}\in {\mathbb{R}}^3\times t\in {\mathbb{R}}_{\ge 0}\right)\to {\mathbb{R}}^3$, $P:(\mathbf{r}\in {\mathbb{R}}^3\times t\in {\mathbb{R}}_{\ge 0})\to \mathbb{R}$, $\mathbf{f}:(\mathbf{r}\in {\mathbb{R}}^3\times t\in {\mathbb{R}}_{\ge 0})\to {\mathbb{R}}^3$, and $\mathbf{v}(\mathbf{r},0)={\mathbf{v}}_0\left(\mathbf{r}\right)\in \{C_c^{\infty }\left({\mathbb{R}}^3\right)\cap L^2\left({\mathbb{R}}^3\right)\}$. Then, there exists a unique velocity field $\mathbf{v}(\mathbf{r},t)$ (128) and pressure field $P(\mathbf{r},t)$ (157) together satisfying (1) and (2) everywhere on all $\Omega ={\mathbb{R}}^3$ for all time $t\in {\mathbb{R}}_{\ge 0}$ under $\mathbf{f}\left(\mathbf{r},t\right)=\mathbf{0}$, obeying Condition 1.

9.2. Proofs of global smoothness of $\mathbf{v}\left(\mathbf{r},t\right)$ and $P\left(\mathbf{r},t\right)$ on $\Omega ={\mathbb{R}}^3$ for all t

Claim 4.2. 

Let $\mathbf{v}:\left(\mathbf{r}\in {\mathbb{R}}^3\times t\in {\mathbb{R}}_{\ge 0}\right)\to {\mathbb{R}}^3$ and $P:(\mathbf{r}\in {\mathbb{R}}^3\times t\in {\mathbb{R}}_{\ge 0})\to \mathbb{R}$, $\mathbf{f}:(\mathbf{r}\in {\mathbb{R}}^3\times t\in {\mathbb{R}}_{\ge 0})\to {\mathbb{R}}^3$, and $\mathbf{v}(\mathbf{r},0)={\mathbf{v}}_0\left(\mathbf{r}\right)\in \{C_c^{\infty }\left({\mathbb{R}}^3\right)\cap L^2\left({\mathbb{R}}^3\right)\}$. Then, the unique velocity field $\mathbf{v}(\mathbf{r},t)$ (Lemma 2.1) and pressure field $P(\mathbf{r},t)$ (Lemma 2.2) together satisfying (1) and (2) everywhere on all $\Omega ={\mathbb{R}}^3$ for all time $t\in {\mathbb{R}}_{\ge 0}$ under $\mathbf{f}\left(\mathbf{r},t\right)=\mathbf{0}$ in compliance with Condition 1 (Lemma 3) are also both $C^{\infty }$ differentiable everywhere on $\Omega ={\mathbb{R}}^3$ for all time $t\in {\mathbb{R}}_{\ge 0}$ in further compliance with Condition 2.

9.2.1. Proof of Global Smoothness of Velocity Field $\mathbf{v}\left(\mathbf{r},t\right)$ on $\Omega ={\mathbb{R}}^3$ for All t

Claim 4.2.1. 

Let $\mathbf{v}:\left(\mathbf{r}\in {\mathbb{R}}^3\times t\in {\mathbb{R}}_{\ge 0}\right)\to {\mathbb{R}}^3$ and $\mathbf{v}(\mathbf{r},0)={\mathbf{v}}_0\left(\mathbf{r}\right)\in L^2\left({\mathbb{R}}^3\right)$. Then, the velocity field $\mathbf{v}(\mathbf{r},t)$ satisfying (1) and (2) everywhere on $\Omega ={\mathbb{R}}^3$ for all time $0\le t<\infty$ in compliance with Condition 1 (Lemma 2.1) is also $C^{\infty }$ everywhere on $\Omega ={\mathbb{R}}^3$ for all time $t\in {\mathbb{R}}_{\ge 0}$ in further compliance with Condition 2.

It is known that the parabolicity of heat equation solutions causes their time evolution to be $C^{\infty }$ smooth everywhere for all time, regardless of the smoothness of their initial conditions [70]. Formally, if ${\Phi }_{k_0}:\mathbf{r}\in {\mathbb{R}}^3\to \mathbb{R}$, then ${\Phi }_k\left(\mathbf{r},t\right)\in C^{\infty }\left({\mathbb{R}}^3(0,t]\right)\forall t\in {\mathbb{R}}_{>0}$. The heat equation initial condition ${\Phi }_{k_0}\left(\mathbf{r}\right)$ has nonetheless been proven to both be unique per (151) and be $C^{\infty }$ smooth for all ${\mathbf{v}}_0\left(\mathbf{r}\right)\in C_c^{\infty }\left({\mathbb{R}}^3\right)$ (155). It therefore follows that ${\Phi }_k\left(\mathbf{r},t\right)\in C^{\infty }\left({\mathbb{R}}^3[0,t]\right)\forall t\in {\mathbb{R}}_{\ge 0}$. The entire coefficient matrix $\mathbf{A}$ (124) is made up entirely of gradient components of heat equation solutions ${\Phi }_k\left(\mathbf{r},t\right)$ which not only make up all entries of both $\mathbf{B}$ and $\mathbf{C}$ as well, but also have been proven to lie within $\{C^{\infty }\left(\mathbb{R}(0,t]\right)\cap L^2\left(\mathbb{R}(0,t]\right)\}\forall t\in {\mathbb{R}}_{\ge 0}$. Since quotients of $C^{\infty }$ functions are also $C^{\infty }$ provided the denominators never vanish, and ${\Phi }_k\left(\mathbf{r},t\right)>0$, it is clear that $a_{ij}\in \{C^{\infty }\left(\mathbb{R}(0,t]\right)\cap L^2\left(\mathbb{R}(0,t]\right)\}\forall \{i\in \{1,2,3\},j\in \{1,2,3\},t\in {\mathbb{R}}_{\ge 0}\}\Rightarrow \{\mathcal{P}\in \{C^{\infty }\left(\mathbb{R}(0,t]\right)\cap L^2\left(\mathbb{R}(0,t]\right)\}\forall \{i\in \{1,2,3\},j\in \{1,2,3\},t\in {\mathbb{R}}_{\ge 0}\}\}\forall \mathcal{P}\in \{b_j,c_j\}$. $C^{\infty }$ smoothness of (128) follows from that of $\mathbf{A}$, $\mathbf{B}$, and $\mathbf{C}$ if there exists a nonempty open set $\mathfrak{U}\subset {\mathbb{R}}^3$ such that for all $p\in \mathfrak{U}$, the following four conditions hold: positivity of the a-coefficients not in $\mathbf{A}$ (158), the diagonal dominance of $\mathbf{A}$ for some ${\delta }_{\mathrm{diag}}>0$ (159), an a priori positive bound $\mathcal{B}$ on $\omega (\mathbf{r},t)$ (160) such that $\omega \left(p\right)\in \mathcal{B}$ for some $0<\underline{\omega }(\mathbf{r},t)<\overline{\omega }(\mathbf{r},t)<\infty$ independent of p, and a uniform Jacobian gap (161) with a uniform diagonal dominance gap on $\mathcal{B}$ (162). These conditions hold by the reasoning in the existence and uniqueness sections.

$$\mathcal{A}>0\forall \mathcal{A}\in \{a_4\left(p\right),a_5\left(p\right),a_9\left(p\right),a_{1}0\left(p\right),a_{1}4\left(p\right),a_{1}5\left(p\right)\}$$

(158)

$$a_{ii}\left(p\right)\ge {\delta }_{\mathrm{diag}}+\sum _{i\ne j}{a}_{ij}\left(p\right)|i\in \{1,2,3\}$$

(159)

$$\mathcal{B}={\left[\underline{\omega }(\mathbf{r},t),\overline{\omega }(\mathbf{r},t)\right]}^3\subset {(0,\infty )}^3$$

(160)

$$D_{\omega }\left(\omega (\mathbf{r},t),p\right)={\displaystyle \frac{\partial F(\mathfrak{U},p)}{\partial \mathfrak{U}}}=\mathbf{C}\left(p\right)\mathrm{diag}\left({\kappa }^{\prime }\left({\omega }_i(\mathbf{r},t)\right)\right)-\mathbf{A}\left(p\right)\mathbf{B}\left(p\right)\mathrm{diag}\left({\eta }^{\prime }\left(\omega (\mathbf{r},t)\right)\right)\left|\right\{{\kappa }^{\prime }\left(\omega (\mathbf{r},t)\right)=$$

$$\left.\omega (\mathbf{r},t)exp\left(\omega (\mathbf{r},t)\right)\left(\omega (\mathbf{r},t)+2\right)>0,{\eta }^{\prime }\left(\omega (\mathbf{r},t)\right)={\displaystyle \frac{-(1+\omega (\mathbf{r},t\left)\right)}{\omega (\mathbf{r},t)exp\left(\omega \right(\mathbf{r},t\left)\right)}}<0\right\}$$

(161)

$$\underset{p\in \mathfrak{U}}{inf}\left(\underset{u\in \mathcal{B}}{inf}\left(\underset{i}{min}{\left(D_{\omega (\mathbf{r},t)}\left(\omega (\mathbf{r},t),p\right)\right)}_{ii}\right)\right)>\underset{p\in \mathfrak{U}}{sup}\left(\underset{u\in \mathcal{B}}{sup}\left(\underset{i}{max}{\left(\sum _{i\ne j}{D}_{\omega (\mathbf{r},t)}\left(\omega (\mathbf{r},t),p\right)\right)}_{ij}\right)\right)$$

(162)

Since $a_{ij}\in C^{\infty }$, $\kappa \in C^{\infty }$ on $(0,\infty )$, and $\eta \in C^{\infty }$ also on $(0,\infty )$, it must also be that $F\in C^{\infty }$ on $(\omega ,p)$. By (162), the system Jacobian is invertible for all $(\omega ,p)\in \mathcal{B}\times \mathfrak{U}$ due to the property of Gershgorin row dominance with positive diagonals regarding the system Jacobian. Taking these two facts together with an invocation of (140) yields $\omega \left(p\right)\in C^1\left(\mathfrak{U}\right)$, while iterating differentiation to arbitrary order leads to $\omega \left(p\right)\in C^{\infty }\left(\mathfrak{U}\right)$. Since $\mathbf{v}(\mathbf{r},t)$ is a composition of both $\omega (\mathbf{r},t)$ and $a_{5i}|i\in \{1,2,3\}$, which lie in $C^{\infty }(0,\infty )$, it must be that $\mathbf{v}(\mathbf{r},t)\in C^{\infty }\left({\mathbb{R}}^3(0,t]\right)\forall t\in {\mathbb{R}}_{\ge 0}$ (163).

$$\omega (\mathbf{r},t)\in C^{\infty }(\mathfrak{U};{\mathbb{R}}^3)\wedge \mathbf{v}(\mathbf{r},t)={\displaystyle \frac{a_5{\widehat{x}}_1}{{\omega }_1(\mathbf{r},t)exp\left({\omega }_1(\mathbf{r},t)\right)}}+{\displaystyle \frac{a_{10}{\widehat{x}}_2}{{\omega }_2(\mathbf{r},t)exp\left({\omega }_2(\mathbf{r},t)\right)}}$$

$$+{\displaystyle \frac{a_{15}{\widehat{x}}_3}{{\omega }_3(\mathbf{r},t)exp\left({\omega }_3(\mathbf{r},t)\right)}}\Rightarrow \mathbf{v}(\mathbf{r},t):\mathfrak{U}\to {\mathbb{R}}^3\Rightarrow \mathbf{v}(\mathbf{r},t)\in C^{\infty }\left({\mathbb{R}}^3(0,t]\right)\forall t\in {\mathbb{R}}_{\ge 0}$$

(163)

Since the implication of the fact that $\mathbf{A}$, $\mathbf{B}$, and $\mathbf{C}$ have unanimously $C^{\infty }$ entries is that $\mathbf{v}(\mathbf{r},t)\in C^{\infty }\left({\mathbb{R}}^3(0,t]\right)\forall t\in {\mathbb{R}}_{\ge 0}$, Claim 4.2.1 stands with its result summarized in Lemma 3.1.

Lemma 3.1. 

Let $\mathbf{v}:\left(\mathbf{r}\in {\mathbb{R}}^3\times t\in {\mathbb{R}}_{\ge 0}\right)\to {\mathbb{R}}^3$ and $\mathbf{v}(\mathbf{r},0)={\mathbf{v}}_0\left(\mathbf{r}\right)\in \{C_c^{\infty }\left({\mathbb{R}}^3\right)\cap L^2\left({\mathbb{R}}^3\right)\}$. Then, the velocity field $\mathbf{v}(\mathbf{r},t)$ (128) satisfying (1) and (2) everywhere on $\Omega ={\mathbb{R}}^3$ for all time $t\in {\mathbb{R}}_{\ge 0}$ in compliance with Condition 1 (Lemma 2.1) is also $C^{\infty }$ differentiable everywhere on $\Omega ={\mathbb{R}}^3$ for all time $t\in {\mathbb{R}}_{\ge 0}$ in further compliance with Condition 2.

9.2.2. Proof of Global Smoothness of Pressure Field $P\left(\mathbf{r},t\right)$ on $\Omega ={\mathbb{R}}^3$ for All t

Claim 4.2.2. 

Let $P:\left(\mathbf{r}\in {\mathbb{R}}^3\times t\in {\mathbb{R}}_{\ge 0}\right)\to \mathbb{R}$. Then, the unique pressure field $P(\mathbf{r},t)$ satisfying (1) and (2) everywhere on $\Omega ={\mathbb{R}}^3$ for all time $0\le t<\infty$ given a corresponding velocity field $\mathbf{v}\left(\mathbf{r},t\right)$ (Lemma 2.2) in compliance with Condition 1 is also $C^{\infty }$ differentiable everywhere on $\Omega ={\mathbb{R}}^3$ for all time $t\in {\mathbb{R}}_{\ge 0}$ in further compliance with Condition 2 under appropriate Poisson decay restrictions.

The pressure field is calculated from the velocity field using an unbounded Poisson integral (157) which inherits its elliptic regularity from its source term, namely the divergence of the material derivative (164) of the velocity field $\mathbf{v}(\mathbf{r},t)$ (128). The elliptic regularity of (164) may be deduced from that of (128) by understanding that because the sum or product of two $C^{\infty }$ smooth functions also lies in $C^{\infty }$, and because the first derivative of a $C^{\infty }$ function also lies in $C^{\infty }$, and finally because a $C^{\infty }$ Poisson source term yields at least a $C_{\mathrm{loc}}^{\infty }$ solution, it must be that at least $\nabla ·{\displaystyle \frac{D\mathbf{v}(\mathbf{r},t)}{Dt}}\in C_{\mathrm{loc}}^{\infty }$. To ensure the stronger condition of global infinite smoothness $\nabla ·{\displaystyle \frac{D\mathbf{v}(\mathbf{r},t)}{Dt}}\in C^{\infty }$, the Poisson source term (164) must decay as $\left|\mathbf{r}\right|\to \infty$ sufficiently fast (165). This condition (165) that the divergence of the material derivative decay strictly faster at infinity than an inverse cubic polynomial is perhaps the strongest condition encountered thus far in this work for decay of the velocity field, but it is nonetheless required to ensure $\mathbf{v}(\mathbf{r},t)\in C^{\infty }$, and finds its manifestation in ultimately restricting the allowable decay rate of ${\mathbf{v}}_0\left(\mathbf{r}\right)$.

$$\nabla ·{\displaystyle \frac{D\mathbf{v}(\mathbf{r},t)}{Dt}}={\displaystyle \frac{\partial }{\partial x_1}}\left({\displaystyle \frac{\partial u(\mathbf{r},t)}{\partial t}}+u(\mathbf{r},t){\displaystyle \frac{\partial u(\mathbf{r},t)}{\partial x_1}}+v(\mathbf{r},t){\displaystyle \frac{\partial u(\mathbf{r},t)}{\partial x_2}}+w(\mathbf{r},t){\displaystyle \frac{\partial u(\mathbf{r},t)}{\partial x_3}}\right)+$$

$${\displaystyle \frac{\partial }{\partial x_2}}\left({\displaystyle \frac{\partial v(\mathbf{r},t)}{\partial t}}+u(\mathbf{r},t){\displaystyle \frac{\partial v(\mathbf{r},t)}{\partial x_1}}+v(\mathbf{r},t){\displaystyle \frac{\partial v(\mathbf{r},t)}{\partial x_2}}+w(\mathbf{r},t){\displaystyle \frac{\partial v(\mathbf{r},t)}{\partial x_3}}\right)+$$

$${\displaystyle \frac{\partial }{\partial x_3}}\left({\displaystyle \frac{\partial w(\mathbf{r},t)}{\partial t}}+u(\mathbf{r},t){\displaystyle \frac{\partial w(\mathbf{r},t)}{\partial x_1}}+v(\mathbf{r},t){\displaystyle \frac{\partial w(\mathbf{r},t)}{\partial x_2}}+w(\mathbf{r},t){\displaystyle \frac{\partial w(\mathbf{r},t)}{\partial x_3}}\right)$$

(164)

$$\nabla ·{\displaystyle \frac{D\mathbf{v}(\mathbf{r},t)}{Dt}}<O\left({\left|\mathbf{r}\right|}^{-3-ϵ}\right)|ϵ>0$$

(165)

Since the velocity field is $C^{\infty }$ smooth, and the pressure field requires that the divergence of the material derivative of the velocity field decay at infinity strictly faster than the cubic polynomial in order to lie in $C^{\infty }$, which ultimately quantifies a minimum allowable decay rate at infinity for the initial velocity field in order to produce an acceptably decaying divergence of the material derivative of the velocity field, Claim 4.2.2 stands with (165) as the Poisson decay condition. The result is summarized in Lemma 3.2.

Lemma 3.2. 

Let $P:\left(\mathbf{r}\in {\mathbb{R}}^3\times t\in {\mathbb{R}}_{\ge 0}\right)\to \mathbb{R}$. Then, the unique pressure field $P(\mathbf{r},t)$ (157) satisfying (1) and (2) everywhere on $\Omega ={\mathbb{R}}^3$ for all time $0\le t<\infty$ given a corresponding velocity field $\mathbf{v}\left(\mathbf{r},t\right)$ (128) (Lemma 2.2) in compliance with Condition 1 is also $C^{\infty }$ differentiable everywhere on $\Omega ={\mathbb{R}}^3$ for all time $t\in {\mathbb{R}}_{\ge 0}$ in further compliance with Condition 2 so long as $\nabla ·{\displaystyle \frac{D\mathbf{v}(\mathbf{r},t)}{Dt}}$ decays strictly faster than the inverse cubic polynomial as $\left|\mathbf{r}\right|\to \infty$.

9.2.3. Conclusion of Smoothness

Because it has been proven that there indeed exists a velocity field $\mathbf{v}(\mathbf{r},t)\in C^{\infty }\left({\mathbb{R}}^3(0,t]\right)$ and pressure field $P(\mathbf{r},t)\in C^{\infty }\left({\mathbb{R}}^3(0,t]\right)$ satisfying (1) and (2) for all $t\in {\mathbb{R}}_{\ge 0}$ given $\mathbf{v}(\mathbf{r},0)={\mathbf{v}}_0\left(\mathbf{r}\right)\in \{C_c^{\infty }\left({\mathbb{R}}^3\right)\cap L^2\left({\mathbb{R}}^3\right)\}$, Claim 4.2 stands with its result summarized in Lemma 4.

Lemma 4. 

Let $\mathbf{v}:\left(\mathbf{r}\in {\mathbb{R}}^3\times t\in {\mathbb{R}}_{\ge 0}\right)\to {\mathbb{R}}^3$, $P:(\mathbf{r}\in {\mathbb{R}}^3\times t\in {\mathbb{R}}_{\ge 0})\to \mathbb{R}$, $\mathbf{f}:(\mathbf{r}\in {\mathbb{R}}^3\times t\in {\mathbb{R}}_{\ge 0})\to {\mathbb{R}}^3$, and $\mathbf{v}(\mathbf{r},0)={\mathbf{v}}_0\left(\mathbf{r}\right)\in \{C_c^{\infty }\left({\mathbb{R}}^3\right)\cap L^2\left({\mathbb{R}}^3\right)\}$. Then, both the unique velocity field $\mathbf{v}(\mathbf{r},t)$ (128) (Lemma 2.1) and corresponding unique pressure field $P(\mathbf{r},t)$ (157) (Lemma 2.2) together satisfying (1) and (2) everywhere on all $\Omega ={\mathbb{R}}^3$ for all time $0\le t<\infty$ under $\mathbf{f}\left(\mathbf{r},t\right)=\mathbf{0}$ and $\mathbf{v}(\mathbf{r},0)={\mathbf{v}}_0\left(\mathbf{r}\right)\in L^2\left({\mathbb{R}}^3\right)$ (Lemma 3) are also both $C^{\infty }$ differentiable everywhere on $\Omega ={\mathbb{R}}^3$ for all time $t\in {\mathbb{R}}_{\ge 0}$ in compliance with Condition 2 so long as $\nabla ·{\displaystyle \frac{D\mathbf{v}(\mathbf{r},t)}{Dt}}$ decays strictly faster than the inverse cubic polynomial as $\left|\mathbf{r}\right|\to \infty$.

9.3. Proofs of Square Integrability of $\mathbf{v}\left(\mathbf{r},t\right)$ and $P\left(\mathbf{r},t\right)$ on $\Omega ={\mathbb{R}}^3$ for all t

The question of whether the Navier-Stokes equations possess bounded kinetic energy per unit volume may be answered by ensuring global square integrability of the velocity and pressure fields. The formal goal of this final section addressing Condition 3 is outlined in Claim 4.3.

Claim 4.3. 

Let $\mathbf{v}:\left(\mathbf{r}\in {\mathbb{R}}^3\times t\in {\mathbb{R}}_{\ge 0}\right)\to {\mathbb{R}}^3$, $P:(\mathbf{r}\in {\mathbb{R}}^3\times t\in {\mathbb{R}}_{\ge 0})\to \mathbb{R}$, $\mathbf{f}:(\mathbf{r}\in {\mathbb{R}}^3\times t\in {\mathbb{R}}_{\ge 0})\to {\mathbb{R}}^3$, and $\mathbf{v}(\mathbf{r},0)={\mathbf{v}}_0\left(\mathbf{r}\right)\in L^2\left({\mathbb{R}}^3\right)$. Then, both the unique velocity field $\mathbf{v}(\mathbf{r},t)$ (128) and corresponding unique pressure field $P(\mathbf{r},t)$ (157) together satisfying (1) and (2) everywhere on all $\Omega ={\mathbb{R}}^3$ for all time $0\le t<\infty$ under $\mathbf{f}\left(\mathbf{r},t\right)=\mathbf{0}$ and $\mathbf{v}(\mathbf{r},0)={\mathbf{v}}_0\left(\mathbf{r}\right)\in \{C_c^{\infty }\left({\mathbb{R}}^3\right)\cap L^2\left({\mathbb{R}}^3\right)\}$ (Lemma 3) which are also both $C^{\infty }$ differentiable everywhere under the same conditions (Lemma 4) in compliance with both Conditions 1 and 2, are rendered $L^2$ square integrable on $\Omega ={\mathbb{R}}^3$ for all time $t\in {\mathbb{R}}_{\ge 0}$ in further compliance with Condition 3.

9.3.1. Proof of Global Square Integrability of Velocity Field $\mathbf{v}\left(\mathbf{r},t\right)$ on $\Omega ={\mathbb{R}}^3$ for All t

As with the existence, uniqueness, and smoothness proofs, the square integrability of the velocity field is proven first because the pressure inherits many of its functional analytical properties from the velocity field. The specific goal of proving that the velocity field is square integrable on $\Omega ={\mathbb{R}}^3$ for all time $t\in {\mathbb{R}}_{\ge 0}$ is outlined in Claim 4.3.1.

Claim 4.3.1. 

Let $\mathbf{v}:\left(\mathbf{r}\in {\mathbb{R}}^3\times t\in {\mathbb{R}}_{\ge 0}\right)\to {\mathbb{R}}^3$, $\mathbf{f}:(\mathbf{r}\in {\mathbb{R}}^3\times t\in {\mathbb{R}}_{\ge 0})\to {\mathbb{R}}^3$, and $\mathbf{v}(\mathbf{r},0)={\mathbf{v}}_0\left(\mathbf{r}\right)\in L^2\left({\mathbb{R}}^3\right)$. Then, the unique velocity field $\mathbf{v}(\mathbf{r},t)$ satisfying (1) and (2) everywhere on all $\Omega ={\mathbb{R}}^3$ for all time $0\le t<\infty$ under $\mathbf{f}\left(\mathbf{r},t\right)=\mathbf{0}$ and $\mathbf{v}(\mathbf{r},0)={\mathbf{v}}_0\left(\mathbf{r}\right)\in \{C_c^{\infty }\left({\mathbb{R}}^3\right)\cap L^2\left({\mathbb{R}}^3\right)\}$ (Lemma 2.1) which is also both $C^{\infty }$ differentiable everywhere under the same conditions (Lemma 3.1) in compliance with both Conditions 1 and 2, is rendered $L^2$ square integrable on $\Omega ={\mathbb{R}}^3$ for all time $t\in {\mathbb{R}}_{\ge 0}$ in further compliance with Condition 3.

It is known that heat equations possess solutions that remain square integrable for all time so long as their initial conditions are square integrable [70] due to the $L^2$-contractivity of the heat semigroup. Formally, if ${\Phi }_{k_0}\left(\mathbf{r}\right)\in L^2\left({\mathbb{R}}^3\right)$, then ${\Phi }_k\left(\mathbf{r},t\right)\in L^2\left({\mathbb{R}}^3(0,t]\right)\forall t\in {\mathbb{R}}_{\ge 0}$. Since ${\Phi }_{k_0}\left(\mathbf{r}\right)$ indeed lies in $\{C_c^{\infty }\left({\mathbb{R}}^3\right)\cap L^2\left({\mathbb{R}}^3\right)\}$ (155) for ${\mathbf{v}}_0\left(\mathbf{r}\right)\in \{C_c^{\infty }\left({\mathbb{R}}^3\right)\cap L^2\left({\mathbb{R}}^3\right)\}$, it follows that ${\Phi }_k\left(\mathbf{r},t\right)\in L^2\left({\mathbb{R}}^3(0,t]\right)\forall \{t\in {\mathbb{R}}_{\ge 0}\wedge k\in \{1,2,3\}$. Because all entries $\{a_{ij},b_j,c_j\}$ of the coefficient matrices $\{\mathbf{A},\mathbf{B},\mathbf{C}\}$ are either ${\Phi }_k\left(\mathbf{r},t\right)$ or gradient components thereof, coupled with the reality that the gradient of a square integrable solution to the heat equation is also square integrable, it is implied that $a_{ij}\in \left\{L^2\left(\mathbb{R}(0,t]\right)\right\}\forall \{i\in \{1,2,3\},j\in \{1,2,3\},t\in {\mathbb{R}}_{\ge 0}\}\Rightarrow \{\mathcal{P}\in \left\{L^2\left(\mathbb{R}(0,t]\right)\right\}\forall \{i\in \{1,2,3\},j\in \{1,2,3\},t\in {\mathbb{R}}_{\ge 0}\}\}\forall \mathcal{P}\in \{b_j,c_j\}$. Square integrability of (128) follows from that of $\mathbf{A}$, $\mathbf{B}$, and $\mathbf{C}$ if, similar to the smoothness argument, there exists a nonempty open set $\mathfrak{U}\subset {\mathbb{R}}^3$ such the following three conditions hold for all $p\in \mathfrak{U}$: mutual positivity and essential boundedness of the a-coefficients not in $\mathbf{A}$ (166) and the nonnegativity of the a-coefficients in $\mathbf{A}$ (167), and the row diagonal dominance of $\mathbf{A}$ for some ${\delta }_{\mathrm{diag}}>0$ (168).

$$\mathcal{A}>0\wedge \underline{a}\le \mathcal{A}\le \overline{a}|\mathcal{A}\in \{a_4,a_5,a_9,a_{10},a_{14},a_{15}\}$$

(166)

$$\mathbf{A}\left(p\right)=\left[\begin{array}{ccc}{a}_1& a_2& a_3\\ a_6& a_7& a_8\\ a_{11}& a_{12}& a_{13}\end{array}\right]\left(p\right)\ge 0$$

(167)

$$a_{ii}\left(p\right)\ge {\delta }_{\mathrm{diag}}+\sum _{i\ne j}{a}_{ij}\left(p\right)|i\in \{1,2,3\}$$

(168)

The fact that $\kappa$ is strictly decreasing on the interval $(0,\infty )$ and $\eta$ is strictly increasing on the same interval establishes the inequalities (169) and (170) whose notation is given by (171).

$$c_i\kappa \left({\omega }_j\right)\le \left(\sum _j{a}_{ij}{b}_j\right)\eta \left({\omega }_{\min }\right)$$

(169)

$$c_i\kappa \left({\omega }_{\min }\right)\ge a_{ii}{b}_i\eta \left({\omega }_{\max }\right)$$

(170)

$$b_j=\left\{a_5,a_{10},a_{15}\right\},c_i\in \left\{{\displaystyle \frac{a_4}{a_5}},{\displaystyle \frac{a_9}{a_{10}}},{\displaystyle \frac{a_{14}}{a_{15}}}\right\}$$

(171)

If ${\omega }_{\min }$ and ${\omega }_{\max }$ are taken as the unique positive solutions of (172), then the constant parameters of (172) obey the three constraints $\underline{c}=\underline{a}/\overline{a}$, $\overline{c}=\overline{a}/\underline{a}$, and ${\overline{a}}_B=3{\left(\overline{a}\right)}^2$ with a fourth constraint ${\underline{a}}_D={\delta }_{\mathrm{diag}}\underline{a}$ implied by the row diagonal dominance (168). The monotonicity of $\kappa$ and $\eta$ yields $0<{\omega }_{\min }\left(p\right)\le {\omega }_i\left(p\right)\le {\omega }_{\max }\left(p\right)<\infty$ uniformly in p for all $i\in \{1,2,3\}$, and therefore for almost all $p\in {\mathbb{R}}^3$, there exists a unique $\omega \left(p\right)\in {[{\omega }_{\min },{\omega }_{\max }]}^3$. Since ${\omega }_{\min }\left(p\right)$ and ${\omega }_{\max }\left(p\right)$ have both been shown to exist and be bounded away from zero and infinity, the constant $M_i$ defining the $L^2$ bound (173)–(175) for the ith component of $\omega \left(\mathbf{r},t\right)$ is derived (176). Using the substitutions (117)–(119) to convert respective (173)–(175) into equivalent statements for $\mathbf{v}\left(\mathbf{r},t\right)$ yield criteria (177)–(179) which firmly seat $\mathbf{v}\left(\mathbf{r},t\right)$ in the $L^2$ function space.

$$\underline{c}\kappa \left({\omega }_{\min }\right)={\overline{a}}_B\eta \left({\omega }_{\min }\right),\overline{c}\kappa \left({\omega }_{\max }\right)={\underline{a}}_D\eta \left({\omega }_{\max }\right)$$

(172)

$$\left|\right|u\left(p\right){\left|\right|}_{L^2}={\left|\left|{\displaystyle \frac{a_5\left(p\right)}{{\omega }_1\left(p\right)exp\left({\omega }_1\left(p\right)\right)}}\right|\right|}_{L^2}\le M_1\left|\right|a_5\left(p\right){\left|\right|}_{L^2}$$

(173)

$$\left|\right|v\left(p\right){\left|\right|}_{L^2}={\left|\left|{\displaystyle \frac{a_{10}\left(p\right)}{{\omega }_2\left(p\right)exp\left({\omega }_2\left(p\right)\right)}}\right|\right|}_{L^2}\le M_2\left|\right|a_{10}\left(p\right){\left|\right|}_{L^2}$$

(174)

$$\left|\right|w\left(p\right){\left|\right|}_{L^2}={\left|\left|{\displaystyle \frac{a_{15}\left(p\right)}{{\omega }_3\left(p\right)exp\left({\omega }_3\left(p\right)\right)}}\right|\right|}_{L^2}\le M_3\left|\right|a_{15}\left(p\right){\left|\right|}_{L^2}$$

(175)

$$M_i=\underset{\omega \left(\mathbf{r},t\right)\in [{\omega }_{\min }\left(\mathbf{r},t\right),{\omega }_{\max }\left(\mathbf{r},t\right)]}{max}\left({\displaystyle \frac{1}{\omega \left(\mathbf{r},t\right)exp\left(\omega \left(\mathbf{r},t\right)\right)}}\right)<\infty$$

(176)

$$\left|\right|u\left(\mathbf{r},t\right){\left|\right|}_{L^2}\le M_1\left|\right|a_5{\left|\right|}_{L^2}$$

(177)

$$\left|\right|v\left(\mathbf{r},t\right){\left|\right|}_{L^2}\le M_2\left|\right|a_{10}{\left|\right|}_{L^2}$$

(178)

$$\left|\right|w\left(\mathbf{r},t\right){\left|\right|}_{L^2}\le M_3\left|\right|a_{15}{\left|\right|}_{L^2}$$

(179)

Since the coefficient matrices $\mathbf{A}$, $\mathbf{B}$, and $\mathbf{C}$ are $L^2$ square integrable, it has been proven that the velocity field is also $L^2$ square integrable. Due to the unique and bounded existence of a $\mathbf{v}\left(\mathbf{r},t\right)\in L^2\left({\mathbb{R}}^3(0,t]\right)$, Claim 4.3.1 stands with the result summarized in Lemma 4.1.

Lemma 4.1. 

Let $\mathbf{v}:\left(\mathbf{r}\in {\mathbb{R}}^3\times t\in {\mathbb{R}}_{\ge 0}\right)\to {\mathbb{R}}^3$, $\mathbf{f}:(\mathbf{r}\in {\mathbb{R}}^3\times t\in {\mathbb{R}}_{\ge 0})\to {\mathbb{R}}^3$, and $\mathbf{v}(\mathbf{r},0)={\mathbf{v}}_0\left(\mathbf{r}\right)\in L^2\left({\mathbb{R}}^3\right)$. Then, the unique velocity field $\mathbf{v}(\mathbf{r},t)$ (128) satisfying (1) and (2) everywhere on all $\Omega ={\mathbb{R}}^3$ for all time $0\le t<\infty$ under $\mathbf{f}\left(\mathbf{r},t\right)=\mathbf{0}$ and $\mathbf{v}(\mathbf{r},0)={\mathbf{v}}_0\left(\mathbf{r}\right)\in \{C_c^{\infty }\left({\mathbb{R}}^3\right)\cap L^2\left({\mathbb{R}}^3\right)\}$ (Lemma 2.1) which is also both $C^{\infty }$ differentiable everywhere under the same conditions (Lemma 3.1) in compliance with both Conditions 1 and 2, is rendered $L^2$ square integrable on $\Omega ={\mathbb{R}}^3$ for all time $t\in {\mathbb{R}}_{\ge 0}$ in further compliance with Condition 3.

9.3.2. Proof of global square integrability of pressure field $P\left(\mathbf{r},t\right)$ on $\Omega ={\mathbb{R}}^3$ for all t

Claim 4.3.2. 

Let $P:(\mathbf{r}\in {\mathbb{R}}^3\times t\in {\mathbb{R}}_{\ge 0})\to \mathbb{R}$, $\mathbf{f}:(\mathbf{r}\in {\mathbb{R}}^3\times t\in {\mathbb{R}}_{\ge 0})\to {\mathbb{R}}^3$, and $\mathbf{v}(\mathbf{r},0)={\mathbf{v}}_0\left(\mathbf{r}\right)\in L^2\left({\mathbb{R}}^3\right)$. Then, the unique pressure field $P(\mathbf{r},t)$ satisfying (1) and (2) everywhere on all $\Omega ={\mathbb{R}}^3$ for all time $t\in {\mathbb{R}}_{\ge 0}$ under $\mathbf{f}\left(\mathbf{r},t\right)=\mathbf{0}$ and $\mathbf{v}(\mathbf{r},0)={\mathbf{v}}_0\left(\mathbf{r}\right)\in \{C_c^{\infty }\left({\mathbb{R}}^3\right)\cap L^2\left({\mathbb{R}}^3\right)\}$ (Lemma 2.2) which are also both $C^{\infty }$ differentiable everywhere under the same conditions (Lemma 3.2) in compliance with both Conditions 1 and 2, is rendered $L^2$ square integrable on $\Omega ={\mathbb{R}}^3$ for all time $0\le t<\infty$ in further compliance with Condition 3.

The calculation of the pressure field from the velocity field using an unbounded Poisson integral (157) with a square integrable source term must be traced through the functional analytic filters of the Poisson equation to guarantee generation of a square integrable Poisson solution for pressure. The source term, namely the divergence of the material derivative (164) of the velocity field $\mathbf{v}(\mathbf{r},t)$ (128), must lie in the Sobolev space $H^{-2}\left({\mathbb{R}}^3(0,t]\right)$ because its Fourier transform must decay strictly faster than the positive branch of the square root function as $\left|\mathbf{r}\right|\to 0$ in order for $P(\mathbf{r},t)\in L^2\left({\mathbb{R}}^3(0,t]\right)$. Since the Poisson equation does not directly deal with time but rather takes time as a parameter, the reality that $\mathbf{v}(\mathbf{r},t)$ already maintains its existence, uniqueness, and smoothness at all times under stated assumptions regarding ${\mathbf{v}}_0\left(\mathbf{r}\right)$ implies that (164) shares the same properties and time may therefore be considered a completed part of the subtask of proving that $P(\mathbf{r},t)\in L^2\left({\mathbb{R}}^3(0,t]\right)$. It is understood that if $\mathbf{v}(\mathbf{r},t)\in L^2\left({\mathbb{R}}^3(0,t]\right)$ (Lemma 4.1), then (164) lies in the Sobolev space $H^{-1}$ in the distributional sense, while taking the divergence of an $H^{-1}$ function yields an $H^{-2}$ function. Therefore, the velocity field $\mathbf{v}(\mathbf{r},t)$ (128) generates a $P(\mathbf{r},t)\in L^2\left({\mathbb{R}}^3(0,t]\right)$ for all $t\in {\mathbb{R}}_{\ge 0}$. Since there exists a $P(\mathbf{r},t)\in L^2\left({\mathbb{R}}^3(0,t]\right)$ paired with a matching $\mathbf{v}(\mathbf{r},t)\in L^2\left({\mathbb{R}}^3(0,t]\right)$, Claim 4.3.2 stands with the result summarized in Lemma 4.2.

Lemma 4.2. 

Let $P:(\mathbf{r}\in {\mathbb{R}}^3\times t\in {\mathbb{R}}_{\ge 0})\to \mathbb{R}$, $\mathbf{f}:(\mathbf{r}\in {\mathbb{R}}^3\times t\in {\mathbb{R}}_{\ge 0})\to {\mathbb{R}}^3$, and $\mathbf{v}(\mathbf{r},0)={\mathbf{v}}_0\left(\mathbf{r}\right)\in L^2\left({\mathbb{R}}^3\right)$. Then, the unique pressure field $P(\mathbf{r},t)$ (157) satisfying (1) and (2) everywhere on all $\Omega ={\mathbb{R}}^3$ for all time $t\in {\mathbb{R}}_{\ge 0}$ under $\mathbf{f}\left(\mathbf{r},t\right)=\mathbf{0}$ and $\mathbf{v}(\mathbf{r},0)={\mathbf{v}}_0\left(\mathbf{r}\right)\in \{C_c^{\infty }\left({\mathbb{R}}^3\right)\cap L^2\left({\mathbb{R}}^3\right)\}$ (Lemma 2.2) which are also both $C^{\infty }$ differentiable everywhere under the same conditions (Lemma 3.2) in compliance with both Conditions 1 and 2, is rendered $L^2$ square integrable on $\Omega ={\mathbb{R}}^3$ for all time $t\in {\mathbb{R}}_{\ge 0}$ in further compliance with Condition 3.

9.3.3. Conclusion of Square Integrability

Because it has been proven that there indeed exists a velocity field $\mathbf{v}(\mathbf{r},t)\in L^2\left({\mathbb{R}}^3(0,t]\right)$ and pressure field $P(\mathbf{r},t)\in L^2\left({\mathbb{R}}^3(0,t]\right)$ satisfying (1) and (2) for all $t\in {\mathbb{R}}_{\ge 0}$ given $\mathbf{v}(\mathbf{r},0)={\mathbf{v}}_0\left(\mathbf{r}\right)\in \{C_c^{\infty }\left({\mathbb{R}}^3\right)\cap L^2\left({\mathbb{R}}^3\right)\}$, Claim 4.3 stands with its result summarized in Lemma 5.

Lemma 5. 

Let $\mathbf{v}:\left(\mathbf{r}\in {\mathbb{R}}^3\times t\in {\mathbb{R}}_{\ge 0}\right)\to {\mathbb{R}}^3$, $P:(\mathbf{r}\in {\mathbb{R}}^3\times t\in {\mathbb{R}}_{\ge 0})\to \mathbb{R}$, $\mathbf{f}:(\mathbf{r}\in {\mathbb{R}}^3\times t\in {\mathbb{R}}_{\ge 0})\to {\mathbb{R}}^3$, and $\mathbf{v}(\mathbf{r},0)={\mathbf{v}}_0\left(\mathbf{r}\right)\in \{C_c^{\infty }\left({\mathbb{R}}^3\right)\cap L^2\left({\mathbb{R}}^3\right)\}$. Then, both the unique velocity field $\mathbf{v}(\mathbf{r},t)$ (128) and corresponding unique pressure field $P(\mathbf{r},t)$ (157) together satisfying (1) and (2) everywhere on all $\Omega ={\mathbb{R}}^3$ for all time $t\in {\mathbb{R}}_{\ge 0}$ under $\mathbf{f}\left(\mathbf{r},t\right)=\mathbf{0}$ (Lemma 3) which are also both $C^{\infty }$ differentiable everywhere under the same conditions (Lemma 4) in compliance with both Conditions 1 and 2, lie furthermore in the function space $L^2$ and are therefore square integrable on $\Omega ={\mathbb{R}}^3$ for all time $t\in {\mathbb{R}}_{\ge 0}$ in further compliance with Condition 3.

9.4. Conclusion of Navier-Stokes Existence and Smoothness

The intersection of the conclusions of Lemmae 3, 4, and 5 yield the statement that $\mathfrak{v}\in \{C^{\infty }\left({\mathbb{R}}^3[0,t]\right)\cap L_2\left({\mathbb{R}}^3[0,t]\right)\}\forall \{\mathfrak{v}\in \{\mathbf{v}\left(\mathbf{r},t\right),P\left(\mathbf{r},t\right)\}\wedge t\in {\mathbb{R}}_{\ge 0}\}$ when $\Omega \left(t\right)={\mathbb{R}}^3$ and ${\mathbf{v}}_0\in \{C_c^{\infty }\left({\mathbb{R}}^3\right)\cap L_2\left({\mathbb{R}}^3\right)\}$ decays sufficiently quickly to satisfy (165). Because such a $\left\{\mathbf{v}\left(\mathbf{r},t\right),P\left(\mathbf{r},t\right)\right\}$ pair meets Conditions 1-3 together, Claim 4 stands with the result stated in Theorem 2. Because both the aforementioned velocity field $\mathbf{v}\left(\mathbf{r},t\right)$ and its corresponding pressure field $P\left(\mathbf{r},t\right)$ satisfying unsimplified (1) and (2) are indeed uniquely defined, infinitely differentiable, and square integrable everywhere for arbitrary $t\in {\mathbb{R}}_{\ge 0}$ and ${\mathbf{v}}_0\left(\mathbf{r}\right)\in \{C_c^{\infty }\left({\mathbb{R}}^3\right)\cap L^2\left({\mathbb{R}}^3\right)\}$, the findings of their functional analytic properties herein constitute a resolution of the Navier-Stokes existence and smoothness problem insofar as it is defined by the Clay Mathematics Institute.

Theorem 2. 

Let $\mathbf{v}:\left(\mathbf{r}\in {\mathbb{R}}^3\times t\in {\mathbb{R}}_{\ge 0}\right)\to {\mathbb{R}}^3$, $P:(\mathbf{r}\in {\mathbb{R}}^3\times t\in {\mathbb{R}}_{\ge 0})\to \mathbb{R}$, $\mathbf{f}:(\mathbf{r}\in {\mathbb{R}}^3\times t\in {\mathbb{R}}_{\ge 0})\to {\mathbb{R}}^3$, and $\mathbf{v}(\mathbf{r},0)={\mathbf{v}}_0\left(\mathbf{r}\right)\in L^2\left({\mathbb{R}}^3\right)$. Then, there exists both a unique velocity field $\mathbf{v}(\mathbf{r},t)\in \{C^{\infty }\left({\mathbb{R}}^3[0,t]\right)\cap L_2\left({\mathbb{R}}^3[0,t]\right)\}\forall t\in {\mathbb{R}}_{\ge 0}$ (128) and a corresponding unique pressure field $P(\mathbf{r},t)\in \{C^{\infty }\left({\mathbb{R}}^3[0,t]\right)\cap L_2\left({\mathbb{R}}^3[0,t]\right)\}\forall t\in {\mathbb{R}}_{\ge 0}$ (157) together satisfying (1) and (2) everywhere on all $\Omega ={\mathbb{R}}^3$ for all time $t\in {\mathbb{R}}_{\ge 0}$ under $\mathbf{f}\left(\mathbf{r},t\right)=\mathbf{0}$ and a $\mathbf{v}(\mathbf{r},0)={\mathbf{v}}_0\left(\mathbf{r}\right)\in \{C_c^{\infty }\left({\mathbb{R}}^3\right)\cap L^2\left({\mathbb{R}}^3\right)\}$, wherein the global smoothness of $P(\mathbf{r},t)$ is conditional on $\nabla ·{\displaystyle \frac{D\mathbf{v}(\mathbf{r},t)}{Dt}}$ decaying strictly faster than the inverse cubic polynomial as $\left|\mathbf{r}\right|\to \infty$ and is otherwise locally smooth (Lemmae 3, 4, 5.)

10. Conclusions

Through a novel multi-layered cascade of known PDE solutions interlaced with MCH transformations ending in a fixed-point overturning of all nine advective terms unwrapping the velocity field into quartic-root deterministic chaos, a highly complex yet equally practical solution of the incompressible Navier-Stokes equations was successfully derived (Theorem 1). This solution was shown to be smooth and globally defined for all time, satisfying requested Clay criteria for the Millennium Problem due to the functional analytic properties of the solution being inherited ultimately from the robust heat equation and Lambert W function (Theorem 2), identifying the commonalities and differences between the Navier-Stokes equations known to engineers and those known to pure mathematicians. Innovative use of vector calculus identities and handling of pressure in the rotational part of the velocity allowed the ultimate separation of the final explicit velocity and pressure fields for their independent use and also in their study, particularly the ways in which both quantities chaotically influence themselves and one another. The use of heat equations as MCH kernels and Poisson equations when handling pressure-related phenomena allowed unbounded flows to be treated separately from bounded flows in the Millennium Problem context, while bounded flows were also given appropriate treatment. Due to the deterministic nature of the velocity and pressure formulae, crisp details of familiar phenomena associated with turbulence such as vortex shedding, chaotic yet physically consistent interactions between pressure and velocity, and large mixing scales are now able to be studied in close detail with the power of a personal computer. Because this general, flexible, and comprehensive analytical solution enables a direct and immediate increase from CFD benchmarks by a factor of several million without any numerical methods at all, this work is expected to bring significant value to the CFD software industry and would thus spread fairly rapidly in some form to virtually any research and development setting in industry or academia involving fluid mechanics. It is entirely possible that new corporate ventures and new technologies will arise due to the effective democratization of Newtonian liquid turbulence simulation, while many other unsolved problems in mathematics, physics, and engineering will be seen to possess solutions implied from those enclosed within this work due to the construction and formalization of novel methods of tackling nonlinear PDEs. These ripple effects are projected to be compounded by the underlying philosophies of nonlinear chaos being shaken, since these results show clearly that fluid turbulence is in no wise a wholly unpredictable, lawless, or truly random phenomenon, but merely an orderly dynamical system too complex for human intuition to fully grasp without mathematical aid. The complicated, intuitively unpredictable nature of turbulence has historically led many a philosopher and poet to speak and write of fluids as having lives of their own, with apparent stability in some situations and violent outbursts in others with often unknown root causes being seen as signs of an inbuilt ability to respond to perceived changes in their environment in a similar manner to an animal. This observation is somewhat sensible in a scientific context insofar as the reality that liquid water is a fundamental constituent of all known life forms, comprising 80 per cent of mammalian physical constitution, and those sciences essential to human life as hemodynamics and viscoelastic tissue deformation depend on at least a generalized or reduced form of (1-2). Indeed, the most common root cause of death in the developed world is neither cancer nor physical accidents nor infectious disease, but atherosclerosis, the progressive obstruction of blood vessels by lipid deposits which eventually starves critical tissues of oxygen. Outside of the human body, atmospheric disturbances such as tropical cyclones and tornadoes are some of the most destructive events encountered in daily life on Earth, and when accounting also for volcanic eruptions and tsunamis, the Navier-Stokes equations are the theoretical physical entities which draw the largest portion of the narrow line between life and death. The track record of storms, floods, heart attacks, strokes, drownings, poisonings, sinking ships, and aviation disasters has frightened civilizations scattered across the planet since time immemorial, and the visceral human fear of the winds and waves becoming roof-tearing gales and wall-melting tsunamis remains largely untempered even amid the triumphant rise of human aviation, spaceflight, and meteorological forecasting. Indeed, since the time it would take a human to solve (1-2) has long been thought to exceed a normative human life expectancy, surely at least one researcher has perished in the struggle to harness the last frontier of classical physics. It is clear then, that if the now freshly slain biblical Leviathan according to the mathematical had been counted among the living creatures, then the apex predator of all known land animals would be neither Felis tigris nor Panthera leo nor even Homo sapiens, but Imperator tumultus.