Terry Moschandreou
Intermediate and Senior Mathematics, Thames Valley District School Board, London, Canada.
DOI: 10.4236/apm.2026.163009   PDF    HTML   XML   84 Downloads   436 Views   Citations

Abstract

This paper presents a framework for solving the Navier-Stokes equations on the 3D torus $T^3$ . In this work, a review is made and a sequel follows with regard to the author’s work in 2021 in [1] where the Geometric Calculus method was used to obtain Equation (6) in that reference. Developing an algebraic and differential procedure to rewrite the Navier-Stokes equations in a form that renders solutions in a particular i-th direction of fluid flow, the expressions solved for give finite time singularities in the form of the LambertW function on the surface of the ball embedded in ${\mathbb{T}}^3={\left[0,1\right]}^3$ . In this paper, it is argued that these singularities for a single $u_i$ direction can be shifted outwards by using fixed point methods via nth compositions of LambertW functions. In this paper, the solution to Equation (6) is provided, which proves that the triple product $u_x{u}_y{u}_z$ has both singular and no finite time blowup. We develop a geometric-analytic framework for the three-dimensional incompressible Navier-Stokes equations on the torus ${\mathbb{T}}^3$ that yields a periodic global attractor and clarifies the role of singular structures arising in nonlinear velocity interactions. Building on a componentwise reformulation of the Navier-Stokes equations, the nonlinear transport term is resolved explicitly in physical space, revealing a rigid normal-form structure governed by compositions of the Lambert $W$ function. Singular solutions initially appear on embedded spherical manifolds through branch-point behavior of $W$ , but we show that these singularities are removable in finite space and can be shifted to infinity via fixed-point methods and iterated Lambert $W$ compositions and increasing Tori approaching ${ℝ}^3$ in the infinite limit. A key mechanism is the synchronization of local Lambert $W$ dynamics with global elliptic structure through a Weierstrass $\wp$ -$\zeta$ mapping. By exploiting homogeneity and quasi-periodicity, a drift-corrected Weierstrass zeta potential is constructed that is strictly periodic on ${\mathbb{T}}^3$ . The resulting velocity field is bounded, oscillatory in time, and invariant under the Navier-Stokes evolution, with secular growth absorbed into the pressure gauge. Elliptic degeneracy controls the singular-regular transition and ensures that all accessible singularities are isolated and integrable. The analysis identifies Lambert $W$ profiles as invariant manifolds of finite codimension and establishes the existence of a periodic attractor for the 3D Navier-Stokes equations on ${\mathbb{T}}^3$ , providing a concrete mechanism by which nonlinear singular behavior is regularized globally.

Keywords

Navier-Stokes Equations, Geometric Depletion, BKM Criterion, Weierstrass -Function, Sphere, Overlap, LambertW, Regular, Smoothness, Alignment, Non-Alignment, Elliptic, Synchronization, Normal Form, Vortex Stretching, Viscosity

Share and Cite:

1. Introduction

The three-dimensional incompressible Navier-Stokes equations on compact manifolds remain a central and unresolved problem in mathematical fluid mechanics. Despite the existence of global weak solutions and extensive partial regularity theory, the mechanisms governing singularity formation, nonlinear depletion, and long-time dynamics in three dimensions are not fully understood. In particular, the interaction between geometry (of embedded spheres in each Tori), nonlinearity, and analytic structure continues to motivate the development of alternative formulations capable of resolving the internal structure of the nonlinear term beyond energy estimates alone. Singularities for a single $u_i$ direction in the Navier-Stokes formulation can be shifted outwards by using fixed point methods via nth compositions of LambertW functions and increasing Tori. This idea has recently been published in [2] and is a continuation of the work in [3]. Beyond the classical velocity-vorticity framework, the relative alignment and non-alignment of velocity components themselves play a decisive role in shaping nonlinear dynamics. When velocity components or dominant modal contributions align locally, the quadratic nonlinearity may undergo partial or complete depletion, suppressing effective energy transfer across scales. This phenomenon is implicit in standard energy methods and is made explicit in Fourier-space analyses, where aligned velocity triads fail to sustain strong cascades [4]. Physically, such alignment corresponds to a reduction in transverse shear and relative motion between interacting flow structures, weakening vortex stretching and nonlinear transport.

Conversely, strong non-alignment between velocity components enhances nonlinear coupling and promotes energy and enstrophy transfer. This distinction is emphasized in classical turbulence theory, where orthogonal or transverse velocity interactions are responsible for sustained nonlinear dynamics and turbulent mixing [5]. From a geometric standpoint, non-alignment introduces additional degrees of freedom that amplify vortex stretching and destabilize purely advective transport. The geometric framework developed here makes this distinction explicit by decomposing the velocity field into radial, tangential, and transverse components on embedded spherical structures, allowing alignment and non-alignment effects to be analyzed directly in physical space rather than solely through spectral methods. The framework in [6] confirms that velocity fields near embedded spherical surfaces admit an orthogonal decomposition into radial, tangential, and binormal components; when radial and tangential components are of comparable magnitude, the binormal direction represents a genuine third degree of freedom rather than a higher-order correction, enabling strong vortex stretching.

The present work focuses on singular solutions supported on spherical surfaces embedded in the three-torus $T^3$ . Rather than treating singularities as obstructions to well-posedness, we show that they can be systematically analyzed, displaced, and ultimately regularized through a combination of fixed-point methods, elliptic structure, and functional composition. In particular, solutions arising from algebraic and differential reformulations of the Navier-Stokes equations naturally involve the Lambert $W$ function, whose branch structure captures implicit nonlinear inversion mechanisms inherent in characteristic-based formulations. Such Lambert $W$ -type singularities appear on the surface of embedded spheres in $T^3$ and reflect degeneracies in the nonlinear interaction geometry rather than genuine finite-time blowup.

A central observation of this paper is that singularities arising in a single velocity component can be shifted outward through iterated composition of Lambert $W$ functions. This idea, recently developed in [2] (Ch. 8), provides a mechanism by which singular behavior on a fixed torus can be transported to arbitrarily large spatial scales by embedding $T^3$ into expanding tori $T_{\alpha }^3={\left[-\alpha ,\alpha \right]}^3$ . As shown in [2], singularities confined to compact regions migrate to infinity as $\alpha \to \infty$ , suggesting that apparent finite-space singularities may correspond instead to infinite-space blowup. This perspective reframes the singularity question in terms of geometric transport rather than intrinsic breakdown.

The triple velocity product $u_x{u}_y{u}_z$ plays a distinguished role in this analysis. Solving the governing PDE derived in [1] reveals that this product admits both singular and non-blowup regimes depending on the accessibility of Lambert $W$ branch points. Crucially, while the Lambert $W$ function possesses a genuine singularity at $W=-1$ , the algebraic structure of the triple product ensures that this singularity does not correspond to finite-time blowup of physically relevant norms. Instead, it manifests as a removable singularity when embedded within an appropriate elliptic and periodic framework. The present work demonstrates that such singularities can be rendered dynamically harmless by mapping the Lambert $W$ evolution onto a Weierstrass elliptic lattice, where quasi-periodicity is neutralized through drift compensation.

This approach connects naturally with broader themes in modern Navier-Stokes analysis. Classical frameworks developed by Leray [7], Fujita-Kato [8], Kato [9], Cannone [10], and Koch-Tataru [11] establish existence, uniqueness, and stability within carefully chosen function spaces such as $L^2$ , $H^s$ , $BM{O}^{-1}$ , and ${\chi }^0$ . However, these approaches control the nonlinear term through a priori estimates and do not resolve its internal algebraic or geometric structure. As a result, they do not single out a preferred functional form for solutions. In contrast, the present work adopts a direct PDE and interaction-geometry perspective, resolving individual nonlinear channels and identifying Lambert $W$ profiles as invariant manifolds under Navier-Stokes evolution.

This philosophy aligns with recent efforts to understand Navier-Stokes dynamics through normal-form transformations, renormalization, and invariant-manifold theory, as well as geometric depletion mechanisms arising from alignment constraints. From this perspective, Lambert $W$ functions arise not as exotic special functions, but as structurally forced objects generated by implicit nonlinear inversion under Navier-Stokes scaling.

The main contribution of this paper is therefore twofold. First, we provide a rigorous derivation of Lambert $W$ -based solutions to the periodic Navier-Stokes equations on $T^3$ , demonstrating closure under nonlinear interactions and compatibility with elliptic pressure constraints. Second, we show that the singular structures inherent in these solutions can be regularized and displaced through composition and elliptic embedding, leading to the existence of a periodic attractor and excluding finite-time blowup on compact domains. By synchronizing internal Lambert $W$ dynamics with external Weierstrass lattice frequencies, we obtain a unique, drift-corrected periodic velocity that is bounded for all time.

More broadly, this work suggests that singularity formation in the three-dimensional Navier-Stokes equations may be governed as much by geometric transport and functional composition as by local amplification mechanisms. The interaction between Lambert $W$ inversion, elliptic degeneracy, and toroidal geometry provides a concrete framework in which singularities can be analyzed, classified, and displaced rather than merely bounded. This perspective opens a new avenue for understanding long-time dynamics, periodic attractors, and the role of implicit nonlinear structure in three-dimensional incompressible flows.

2. Method

The methodology adopts a direct PDE and interaction-geometry perspective, resolving individual nonlinear channels and identifying Lambert $W$ profiles as invariant manifolds under Navier-Stokes evolution. These normal forms can be shifted outward through iterative Lambert W compositions and increasing tori, which when combined with the Weierstrass framework, embeds into globally periodic, bounded dynamics. This provides the bridge from local interaction geometry to the global periodic attractor described in this work.

2.1. Direct PDE and Interaction-Geometry Perspective

The geometric discussion thus far emphasizes that nonlinear activity in the Navier-Stokes equations is governed not merely by vorticity, but by the interaction geometry of velocity components themselves. Alignment suppresses nonlinear transfer, while non-alignment activates it. This observation motivates a formulation in which nonlinear interactions are resolved componentwise and geometrically, rather than being absorbed into a priori functional estimates. The present section makes this perspective explicit at the level of the partial differential equations.

We begin by writing the velocity field in components,

$u=\left(u_x,u_y,u_z\right),$

So that the quadratic nonlinearity decomposes as

$\left(u\cdot \nabla \right)u={\displaystyle \sum }_{i,j\in \left\{x,y,z\right\}}{u}_i {\partial }_i{u}_j.$

Each term $u_i{\partial }_i{u}_j$ represents a distinct interaction channel, encoding the transport of the $j$ -th velocity component along the $i$ -th direction. In contrast to global energy methods, which estimate this sum as a whole, we analyze these channels individually, tracking how geometry, scaling, and functional form propagate through the dynamics.

To formalize this, define the interaction operators

${ℐ}_{ij}\left(u\right):=u_i{\partial }_i{u}_j.$

The nonlinear structure of the Navier-Stokes equations may then be viewed as the superposition of the actions of ${ℐ}_{ij}$ over all component pairs. This viewpoint naturally leads to the consideration of velocity profiles that are closed under all interaction channels, meaning that no new singular structures are generated when any ${ℐ}_{ij}$ acts on the profile.

We therefore consider the set

$\mathcal{S}:={\displaystyle \underset{i,j\in \left\{x,y,z\right\}}{\cap }\text{ker}}\left(\text{singular growth induced by}{ℐ}_{ij}\right),$

whose elements are velocity fields invariant, in a structural sense, under the full nonlinear interaction geometry. Being in $\mathcal{S}$ imposes severe rigidity: such profiles must be stable under multiplication, differentiation, and implicit inversion, while remaining compatible with Navier-Stokes scaling and incompressibility.

This rigidity is precisely what singles out the Lambert $W$ function. If

$f\left(x\right)=W\left(g\left(x\right)\right),$

then

$f^{\prime }\left(x\right)=\frac{g^{\prime }\left(x\right)}{g\left(x\right) \left[1+W\left(g\left(x\right)\right)\right]},$

so differentiation introduces no new transcendental structure beyond rational combinations of $W$ itself. Moreover, equations of the characteristic form

${\partial }_{t}u+u{\partial }_{z}u=\Phi \left(u\right),$

which arise naturally when a dominant transport direction is isolated (as in the geometric decomposition discussed in the Introduction), lead via the method of characteristics to implicit relations of the type

$u{\text{e}}^u=\text{known function}.$

The unique inversion of such relations is given by the Lambert $W$ function.

This observation is reinforced by scaling considerations. Under the Navier-Stokes scaling

$x\mapsto \lambda x,t\mapsto {\lambda }^{2}t,u\mapsto {\lambda }^{-1}u,$

Logarithmic corrections arise that destabilize purely polynomial or exponential ansätze. The Lambert $W$ function naturally absorbs these corrections, remaining invariant in form under renormalization. In this sense, Lambert $W$ profiles are not ad hoc solutions, but normal forms adapted to the intrinsic scaling of the equations.

Taken together, these considerations motivate the following structural claim.

Conjecture (Normal-Form Universality). Among all locally self-similar, renormalization-invariant velocity profiles that are closed under the full set of nonlinear interaction operators ${ℐ}_{ij}$ , the only admissible explicit or implicit solution class is generated by the Lambert $W$ function.

This conjecture is not an existence theorem but a classification principle. It asserts that Lambert $W$ profiles form invariant manifolds-of finite codimension-within the Navier-Stokes flow. From the interaction-geometry viewpoint introduced above, these manifolds represent configurations in which alignment and non-alignment effects balance in such a way that nonlinear growth is structurally constrained rather than merely estimated.

Crucially, this perspective explains why singular behavior appears naturally on embedded geometric structures, such as spheres in $T^3$ , as discussed later in this paper. The singularities are not artifacts of the method, but geometric markers of where interaction channels become degenerate. As shown in subsequent sections, these singularities can be shifted outward through iterative Lambert $W$ compositions and, when combined with the Weierstrass framework, embedded into globally periodic, bounded dynamics. This provides the bridge from local interaction geometry to the global periodic attractor described in this work.

It is claimed that Lambert $W$ profiles are invariant manifolds of finite codimension in the Navier-Stokes flow and there is closure under all interaction operators ${ℐ}_{ij}$ .

A main theorem which leads to the LambertW-$\wp$ closure solutions is as follows.

Theorem 1 (Existence of W function solutions to the 3D Navier-Stokes problem that are given by Equation (65): PDE) Assuming the Poisson equation,

$u_i^2{\nabla }^{2}P\left(x,t\right)=S\left(x,t\right)$ (1)

And the following PDE holds,

${\left(\frac{{\partial }^2{V}_z}{\partial z\partial t}\right)}^2=\frac{1}{\kappa \left(x,y,z,t\right)}{\left(\frac{\partial V_z}{\partial t}\right)}^3$ (2)

where $u_i=V_z+C$ , $C$ is a shift related to $\kappa :\left({ℝ}^3,R^{+}\right)\to ℝ$ where $V_z$ is at least in $C^0\left(x,t\right)$ with the pressure $P$ in $C^2\left(x,t\right)$ , then,

$\kappa \left(x,y,z,t\right)=\frac{2^{2/3}3}{36C_3}{\left[\frac{\partial }{\partial z}{\wp }^{-1}\left(\frac{\partial }{\partial t}W\left(-{\text{e}}^{\Xi -1}\right)\right)\right]}^{-2}.$ (3)

where $C_3\in R$ and $\wp$ is the WeierstrassP function, ${\wp }^{-1}$ is it’s inverse and $W$ is the LambertW function defined on an affine spatio-temporal space $\text{Ξ}=x-y+2z-t+C_6$ , with $C_6=O\left(\delta \right)$ . The expression for $\delta$ is connected to $\kappa$ as $\delta =\frac{1}{\sqrt{\kappa 1}}$ , where $\kappa 1$ is the derivative of $\kappa$ . In principle both $\kappa$ and $\delta$ are functions in general, where they become constants when we look at large spatial gradients in a given i-th direction, in particular for $V_z$ in the $z$ direction. Then, the solution problem can be defined by Equations (5) - (7) which can be shown to reduce to Equation (65). Moreover, the solutions given by $u_i$ or $V_z$ are periodic in both space and time by means of the Lambert to Weierstrass mapping.

Now to establish the grounds for these assertions made in the main theorem we are required to review the work in the proofs of the connection of the form of the Navier-Stokes equations to the PDEs in [1]-[3] and references therein. Since the work in these references were valid in the $L^p$ spaces for $p\in \left[2,3\right)$ , a note here is required. The Poisson equation was used to relate the velocities to the pressure terms. In order to remain in these spaces for $p$ (in particular it was found that $p>3/2$ is necessary), the following PDE was required in the definition of the Poisson equation,

$u_i^2{\nabla }^{2}P\left(x,t\right)=S\left(x,t\right)$ (4)

where $u_i$ are the components of the Navier-Stokes flow with $P$ representing the pressure and $S$ some arbitrary function used in the solution approach. To be specific the solution in the previous references made use of this PDE (Poisson equation). In [1] and [2], a geometric calculus approach was used to rewrite the Navier-Stokes equations in a general $u_i$ direction for $i=1,2,3$ . The PDEs defining the $u_i$ were possible to develop mainly due to the Poisson equation. Thus, we use the $\frac{\partial P}{\partial x_i}$ term in the governing PDE developed by Geometric Calculus approach. The transition to the PDE in Equations (5) - (7) is a result of adding to the original Navier-Stokes equations (after applying row operations [1]-[3], a pivot function $a\nabla \cdot a$ which is used as a place Holder in obtaining the solution of the Navier-Stokes equations. See Equation (6) in [1] where this pivot function has been added. In the subsequent parts of the paper there the pivot function is used in the following sense: If we have two operators $L_1=L_2$ and $L_1=L_3$ then necessarily $L_2=L_3$ so that the appearance of $L_1$ is not seen anymore. It is in this way that the PDE problem for the Navier-Stokes equations was set up. As mentioned, the following PDE captures the dynamics of the original Navier-Stokes equations, which now are written for the $u_i$ direction.

The 3D incompressible unsteady Navier-Stokes Equations (NSEs) in Cartesian coordinates may be listed for the velocity field,

$\rho \left(\frac{\partial }{\partial t^{\ast }}+u^{\ast j}{\nabla }_{\ast j}\right)u_i^{\ast }-\mu {\nabla }_{\ast }^2{u}_i^{\ast }+{\nabla }_{\ast i}{P}^{\ast }=\rho F_i^{\ast }$

where $\rho$ is constant density, $\mu$ is dynamic viscosity, $F^{\ast }=F^{\ast i}{e}_i$ are the body forces on the fluid. In some cases, it may be elected to reparametrize the components of the velocity vector, and pressure to $u={\left(u\right)}^i{e}_i$ , $P={\left(P\right)}^i{e}_i$ , coordinates $x_i$ and time $t$ according to the following form utilizing the non-dimensional quantity $\delta$ ($\delta \le 0$ ) :

$u_i^{\ast }=\frac{1}{\delta }{u}_i,P_i^{\ast }=\frac{1}{{\delta }^2}{P}_i,x_i^{\ast }=\delta x_i,t^{\ast }={\delta }^{2}t$

The Navier-Stokes equations above in $\ast$ variables are proven to be equivalent to the following PDE [1] [2] in non-star variables for the $u_i$ component and $x_i$ , $t$ variables, (here $u_i=V_z+\sqrt{\kappa 1}$ (for example $\kappa 1$ is a derivative of $\kappa$ in a preferred direction in space or time)):

$P=P_1+{\left(\frac{\partial V_z}{\partial t}\right)}^2{P}_2+\delta S\left(x,t\right){\left(\frac{\partial V_z}{\partial t}\right)}^2$ (5)

where $S={\left(V_z+\sqrt{\kappa 1}\right)}^2{\nabla }^{2}P$ ,

$\begin{array}{l}{P}_1=\frac{2\delta \left(\frac{\partial V_z}{\partial t}\right)\left(\sqrt{\kappa 1}{V}_z+\frac{V_z^2}{2}+\frac{\kappa 1}{2}\right)\rho \left(\frac{{\partial }^3{V}_z}{\partial t\partial z^2}\right)}{3}\\ +\frac{{\left(\frac{\partial V_z}{\partial t}\right)}^2\mu \left(\delta -1\right)\left(\frac{{\partial }^3{V}_z}{\partial x^2\partial z}\right)}{3}+\frac{{\left(\frac{\partial V_z}{\partial t}\right)}^2\mu \left(\delta -1\right)\left(\frac{{\partial }^3{V}_z}{\partial y^2\partial z}\right)}{3}\\ +\frac{{\left(\frac{\partial V_z}{\partial t}\right)}^2\mu \left(\delta -1\right)\left(\frac{{\partial }^3{V}_z}{\partial z^3}\right)}{3}+\frac{2\delta \left(\sqrt{\kappa 1}{V}_z+\frac{V{z}^2}{2}+\frac{\kappa 1}{2}\right)\rho {\left(\frac{{\partial }^2{V}_z}{\partial t\partial z}\right)}^2}{3}\\ +\frac{2\left(\frac{\left(-\delta +1\right){\left(\frac{\partial V_z}{\partial t}\right)}^2}{2}+\left(V_z+\sqrt{\kappa 1}\right)\left(\frac{\partial V_z}{\partial t}\right)\left(\frac{\partial V_z}{\partial z}\right)\delta -\frac{\delta \Phi \left(t\right)}{2}\right)\rho \left(\frac{{\partial }^2{V}_z}{\partial t\partial z}\right)}{3}\end{array}$ (6)

$\begin{array}{l}{P}_2:=\left(\frac{2}{3}+\left(\rho -\frac{2}{3}\right)\delta \right)\left(V_z+\sqrt{\kappa 1}\right)\left(\frac{{\partial }^2{V}_z}{\partial z^2}\right)-\frac{2\left(V_z+\sqrt{\kappa 1}\right)\left(\delta -1\right)\left(\frac{{\partial }^2{V}_z}{\partial x^2}\right)}{3}\\ -\frac{2\left(V_z+\sqrt{\kappa 1}\right)\left(\delta -1\right)\left(\frac{{\partial }^2{V}_z}{\partial y^2}\right)}{3}+\frac{\left(-\delta +1\right)\left(\frac{{\partial }^2{V}_z}{\partial x\partial z}\right)}{3}+\frac{\left(-\delta +1\right)\left(\frac{{\partial }^2{V}_z}{\partial y\partial z}\right)}{3}\\ +\frac{\left(1+\left(3\rho -1\right)\delta \right){\left(\frac{\partial V_z}{\partial z}\right)}^2}{3}-\frac{\delta }{3}+\frac{1}{3}\end{array}$ (7)

where $S$ is defined in Equation (1). Note the Laplacian for the pressure is written as an integral over an epsilon ball along each of the infinitely many branches of the WeierstrassP function. There are precisely three Laplacians in Equations (5) - (7), one is for the pressure and the other two are in terms of the velocity $V_z$ . The work of Rumer and Fet was used in [3] and reference within to write the Laplacians as integrals over epsilon balls. In [2], there, $x,y$ components $V_x,V_y$ vanish on a suitable manifold as shown in this paper and in [3] (where the space $J_{y_i}$ was defined with a calculation showing that an operator involving all three velocity components and their derivatives, $\mathcal{X}$ , is precisely zero on this space and we see it to be true on the boundary of an embedded ball in ${\mathbb{T}}^3={\left[0,1\right]}^3$ (see Section 4). It is important to define $\delta$ which is used in two ways in this paper. The first is that we assume alignment of two vector fields in general and then separately non-alignment. The expression for $\delta$ is that it is the negative reciprocal of $\kappa 1\left(x,y,z,t\right)$ ($\kappa 1$ is the derivative $\kappa \text{'}$ ) and the following eigen-type problem holds true:

${\left(\frac{{\partial }^2{V}_z}{\partial z\partial t}\right)}^2=\frac{1}{\kappa \left(x,y,z,t\right)}{\left(\frac{\partial V_z}{\partial t}\right)}^3$ (8)

The solution of this PDE gives a general form in terms of the WeierstrassP function (consistent with [2] and [3]). It is:

$V_z\left(x,y,z,t\right)=2^{\frac{2}{3}}\left({\displaystyle \int \wp \left({\displaystyle \int \frac{2^{\frac{2}{3}}\sqrt{3}}{6\sqrt{\kappa \left(x,y,z,t\right)2^{\frac{2}{3}}}}\text{d}z}+f_1\left(x,y,t\right);g_2,g_3\right)\text{d}t}\right)+f_2\left(x,y,z\right)$

In the PDE given by Equations (5) - (7), the expression $\left(\frac{{\partial }^2{V}_z}{\partial z\partial t}\right)$ appears at a few places. Substituting in Equations (5) - (7) introduces $\frac{1}{\kappa }$ at these places and then multiplying by $\kappa$ throughout aligns $\delta$ with $\kappa$ so that we have the PDE which gives the LambertW solution. It is a straightforward bookkeeping approach to see that this results as in Equation (65) in the later section where we have defined the representative governing equations. Attention must be given to the $P_1$ operator expression in Equation (6). Here the first term and sixth term expressions (the sixth part contains also $\Phi \left(t\right)$ which is set to zero) (in the sum of 6 parts) when added together vanish. If this is due to $\kappa =0$ then one can show that there is finite time blowup which must occur at a pole. (Substitute $\kappa =0$ into the $P_1$ operator part and solve this set to zero [2].) However, we exclude the finite blowup case since $\kappa \ne 0$ must hold in Equation (8). What remains in Equations (5) - (7) will lead to Equation (65) which solves as a LambertW solution.

The quantity $\kappa \left(x,t\right)=\kappa \left(x,y,z,t\right)$ is a refinement geometry parameter measuring the transverse concentration of vorticity relative to its local direction. Let,

$\omega \left(x,t\right)=\nabla \times u\left(x,t\right),\xi \left(x,t\right)=\frac{\omega \left(x,t\right)}{\left|\omega \left(x,t\right)\right|}$

denote the vorticity and vorticity direction field wherever $\omega \ne 0$ .

We define the transverse gradient operator

${\nabla }_{\perp }:=\left(I-\xi \otimes \xi \right)\nabla ,$

which projects spatial derivatives onto the plane orthogonal to the local vorticity direction.

Definition 1 (Refinement Geometry $\kappa$ ) The refinement geometry $\kappa \left(x,t\right)=1/{\kappa }_{geom}\left(x,t\right)$ is defined (up to universal constants) by

$\boxed{{\kappa }_{geom}\left(x,t\right)~\frac{{\left|{\nabla }_{\perp }\omega \left(x,t\right)\right|}^2}{{\left|\omega \left(x,t\right)\right|}^2}.}$

The quantity ${\kappa }_{geom}$ is scale-invariant under the Navier-Stokes scaling and diverges precisely when vorticity concentrates transversely.

Equivalently, if $r_{\perp }\left(x,t\right)$ denotes the local transverse filament radius associated with the vortex tube geometry, then

$\boxed{\kappa \left(x,t\right)~r_{\perp }{\left(x,t\right)}^2.}$

Thus:

• small $\kappa$ corresponds to thin, highly refined vortex filaments, strong transverse gradients,

• large $\kappa$ corresponds to thick vortex filaments, well spread vorticity,

• bounded $\kappa$ prevents transverse collapse.

2.1.1. Connection of $r_{\perp }~\sqrt{x_i}$ to $W$ Solution-Lambert $W$ Branch Geometry and Transverse Gradient Scaling

We consider the function

$W\text{}\left(-e^{-x-1}\right),$

where $W$ denotes either real branch $W_0$ or $W_{-1}$ of the Lambert $W$ function. The argument approaches the branch point

$z_{\ast }=-e^{-1}\text{as}x\to 0^{+}.$

A classical expansion of the Lambert $W$ function near this branch point yields

$W\left(z\right)=-1\pm \sqrt{2\left(ez+1\right)}+O\left(ez+1\right).$

Substituting $z=-e^{-x-1}$ gives

$ez+1=1-e^{-x}~x\left(x\to 0^{+}\right),$

and therefore

$W\left(-e^{-x-1}\right)+1~\pm \sqrt{2x}.$ (9)

Differentiating (9) with respect to $x$ yields

${\partial }_{x}W\left(-e^{-x-1}\right)~\frac{1}{\sqrt{x}}\text{as}x\to 0^{+}.$ (10)

This divergence reflects a square-root branch singularity. The function itself remains finite at $x=0$ , while its first derivative becomes unbounded.

Assume that $x$ measures squared transverse distance from a vortex tube axis, so that

$r_{\perp }~\sqrt{x}.$

Then

${\partial }_{\perp }~{\partial }_x~\frac{1}{\sqrt{x}}~\frac{1}{r_{\perp }}.$

Suppose the longitudinal velocity component has the form

$u_z\left(x\right)~W\left(-e^{-x-1}\right)+\text{smooth terms}.$

Using (10),

${\partial }_{\perp }{u}_z~\frac{1}{\sqrt{x}}~\frac{1}{r_{\perp }}.$

Define the geometric curvature parameter

${\kappa }_{\text{geom}}=\frac{{\left|{\nabla }_{\perp }\omega \right|}^2}{{\left|\omega \right|}^2}.$

The above scaling implies

${\kappa }_{\text{geom}}~\frac{1}{r_{\perp }^2}.$

Equivalently, defining

$\kappa :={\left(\kappa \right)}_{\text{geom}}^{-1},$

we obtain

$\boxed{\kappa ~r_{\perp }^2.}$

The Lambert $W$ structure encodes a loss of invertibility of the flow map in a single transverse direction, producing strong but integrable transverse gradients. This square-root behavior corresponds geometrically to a fold or caustic in the overlapping-sphere (or bent-tube) construction, not to a velocity blowup.

In particular:

• the velocity remains bounded,

• transverse gradients scale like $1/r_{\perp }$ ,

• the curvature parameter scales like $1/r_{\perp }^2$ ,

• and the singularity is purely geometric, arising from branch structure.

In Equation (8), through the transformations from nonstar to $\ast$ variables, keeping in mind that $\kappa$ must be transformed as $\kappa =\kappa \left(x^{\text{*}},y^{\text{*}},z^{\text{*}},t^{\text{*}}\right)$ , Equation (8) will be transformed so that $\delta$ disappears and everything is in $\ast$ variables. The vortex tube has as its envelope $\pm 1/\sqrt{x}$ where the physical vortex tube matches ± of the derivative of the $W$ solution as seen previously.

Finally, for a near bent vortex tube, one can have large longitudinal derivative of the $z$ component velocity as long as $\kappa ~1/r_{\perp }^2$ and it remains finite.

2.1.2. Vortex Stretching Control via $\kappa$ and Weierstrass $\wp$

Let $\kappa =\kappa \left(x,y,z,t\right)>0$ be a smooth function, and define

$U\left(x,y,z,t\right):={\displaystyle {\int }_0^z\frac{1}{\sqrt{\kappa \left(x,y,z^{\prime },t\right)}}\text{d}{z}^{\prime }}.$

Let $\wp \left(U;g_2,g_3\right)$ denote the Weierstrass $\wp$ -function with invariants $g_2,g_3\in ℝ$ , and let

${\wp }^{\prime }\left(U;g_2,g_3\right):=\frac{\text{d}}{\text{d}U}\wp \left(U;g_2,g_3\right)$

be its derivative.

Theorem 2 (Stretching Amplification via Small $\kappa$ ) Let

$U\left(x,y,z,t\right)={\displaystyle \int \frac{1}{\sqrt{\kappa \left(x,y,z,t\right)}}\text{d}z},$

where $\kappa \left(x,y,z,t\right)>0$ is smooth, and consider the Weierstrass elliptic function $\wp \left(U;g_2,g_3\right)$ with invariants $g_2,g_3\in ℝ$ such that the cubic

$4s^3-g_{2}s-g_3=0$

has three distinct real roots $e_1>e_2>e_3$ .

Let ${\wp }^{\prime }\left(U\right)$ denote the derivative of $\wp$ with respect to its argument. Then, the following hold:

1) (Pole-controlled growth) If $U\left(x,y,z,t\right)$ approaches a lattice point (pole) of $\wp$ , i.e.,

$U\to \omega \in \text{Λ}=\left\{2m{\omega }_1+2n{\omega }_2:m,n\in \mathbb{Z}\right\},$

then

$\left|{\wp }^{\prime }\left(U\right)\right|\to \infty ,$

and this growth is independent of the size of $\kappa$ .

2) (Root locations are bounded) The zeros of ${\wp }^{\prime }\left(U\right)$ occur precisely at the half-periods

$U\in \left\{{\omega }_1,{\omega }_2,{\omega }_3\right\},{\omega }_3={\omega }_1+{\omega }_2,$

where

$\wp \left({\omega }_i\right)=e_i.$

Hence, ${\wp }^{\prime }\left(U\right)=0$ only at these discrete points.

3) (Effect of small $\kappa$ ) If $\kappa \left(x,y,z,t\right)$ becomes small along a segment of a vortex tube, then

$\left|\frac{\partial U}{\partial z}\right|=\frac{1}{\sqrt{\kappa }}$

becomes large, causing $U$ to vary rapidly with respect to $z$ . This increases the likelihood that $U$ crosses a pole location in finite spatial distance.

4) (Amplification criterion) Large values of

$\left|{\wp }^{\prime }\left(U\left(x,y,z,t\right)\right)\right|$

occur only if $U$ is close to a pole of $\wp$ . Being near the midpoints between roots $e_1,e_2,e_3$ (i.e., values where $\wp \left(U\right)=e_i$ ) instead yields

${\wp }^{\prime }\left(U\right)\approx 0.$

5) (Conclusion) Small $\kappa$ alone does not directly make ${\wp }^{\prime }\left(U\right)$ large; rather, it acts through the mapping $U\left(z\right)$ , potentially driving $U$ toward a pole of $\wp$ . Thus, vortex stretching amplification in this formulation is a geometric effect of pole proximity, not merely of the magnitude of $\kappa$ .

2.1.3. Analysis of the Operator Reduced to Zero in $P_1$ in Equation (6) and $P_2$ in Equation (7)

We consider the addition of the first and sixth expressions previously in Equation (6) together with the very last term in Equation (7) and pressure is linear in space,

that is, after factoring a $\frac{\partial V_z}{\partial t}$ common term out. The following section provides a solution of $\kappa$ once $V_z$ is determined.

2.1.4. Exact Solution of the Extended PDE via the Weierstrass Ansatz for General Spatio-Temporal Pressure

We consider the PDE for the function $V_z\left(x,y,z,t\right)$ :

$\begin{array}{c}E\left[V_z\text{, }\kappa \text{, }P\right]=\frac{2}{3}\left(\frac{\partial u_z}{\partial t}\right)u_z^2\left(\frac{{\partial }^3{u}_z}{\partial t\partial z^2}\right)\delta \rho \\ +\left[\frac{2}{3}\left(\frac{1}{2}\left(1-\delta \right){\left(\frac{\partial u_z}{\partial t}\right)}^2+u_z\frac{\partial u_z}{\partial t}\frac{\partial u_z}{\partial z}\delta \right)\right]\rho \frac{{\partial }^2{u}_z}{\partial t\partial z}\\ -4u_z^2{\left(\frac{\partial u_z}{\partial t}\right)}^2\left(\frac{1}{2}{L}_p\delta \right)+\left(1-\delta \right){\left(\frac{\partial u_z}{\partial t}\right)}^2\end{array}$

where, $u_i=V_z+\sqrt{\kappa 1}$ , $V_z$ is the function to be determined, $\kappa 1\left(x,y,z,t\right)$ is a functional to be chosen ($\kappa 1$ is shown in Figure 1 in blue colour for the solution of Equation (8) with $v_3=V_z$ given by Equation (65)), and $P\left(x,y,z,t\right)$ is the pressure. Here, we have a general situation where pressure is coupled with the PDE and once $E=0$ , we are left with the remaining part of Equations (5) - (7) which leads to Equation (65). Here, there was a shift in $z$ for this $\kappa 1$ definition. It was shown that $z\to z+0.25$ approximately will lead to the $\kappa 1$ representations to be the same. Also if $\kappa 1=\kappa \text{'}>0$ then for $z_0>0$ , $\kappa >z_0>0$ for all $z>0$ . This is important since $\kappa$ must be greater than or equal to zero in the square root in the solution obtained. For example, if ${\nabla }^{2}P=0$ as a special case, we have the plot for $\kappa 1\left(z\right)$ in 1 in red (solve for $\kappa 1$ in $E$ above).

We have that the derivatives increase in order when the pressure is harmonic. This is precisely why we consider the derivative of $\kappa$ , for example wrt to $z$ , that is $\kappa 1$ . Here, we have taken the solution given in Equation (65) which is $V_z=\text{LambertW}\left(-\text{exp}\left(x-y+2z-t-1\right)\right)$ and substituted in the PDE $E\left[V_z,\kappa ,P\right]$ above and solved for $\kappa 1$ .

Figure 1. Special case of $\kappa 1\left(z\right)$ vs. $z$ for $\delta =1/\sqrt{\kappa 1}$ , where $\kappa 1={\kappa }^{\prime }$ , ${\kappa }^{\prime }=\frac{\partial \kappa }{\partial z}$ . The blue graph is based on Equation (8) and the red graph is based on equation defining $E\left[V_z,\kappa 1,P\right]$ . Changing the $\kappa 1$ due to a small shift horizontally and vertically matches them exactly and also the $\kappa \ge 0$ definitions.

2.1.5. Finding $\kappa$ Using the ${\wp }^{-1}$ Function

Now, the interesting finding is that if we solve Equation (8) for $\kappa$ and differentiate it wrt to $z$ , we have $\kappa 1$ exactly as in 1.

We investigate the relationship between a geometric field $\kappa \left(x,y,z,t\right)$ and a potential function $H\left(x,y,z,t\right)$ mediated by the inverse Weierstrass elliptic function ${\wp }^{-1}$ .

The governing integral equation is defined as:

${\displaystyle \int \frac{1}{6}\frac{2^{2/3}\sqrt{3}}{\sqrt{\kappa \left(x,y,z,t\right)2^{2/3}}}\text{d}z}+f_1\left(x,y,t\right)={\wp }^{-1}\left(\frac{\partial H\left(x,y,z,t\right)}{\partial t};g_2,g_3\right)$ (11)

where $g_2$ and $g_3$ are the Weierstrass invariants and $f_1$ is an arbitrary function of integration. Our objective is to determine the conditions under which the field $\kappa$ exhibits a Hölder singularity with exponent $\alpha =1/3$ at $y=0$ (taking $z=-2y$ ), particularly reproducing a profile characterized by rapid growth for $y<0$ and asymptotic growth for $y>0$ .

Differential Formulation and Inverse Function Properties

To analyze the local behavior of $\kappa$ , we differentiate both sides of the integral equation with respect to the spatial coordinate $z$ . Applying the Fundamental Theorem of Calculus to the left-hand side and the chain rule to the right-hand side yields:

$\frac{1}{\sqrt{12\kappa }}=\frac{\partial }{\partial z}\left[{\wp }^{-1}\left(H_t;g_2,g_3\right)\right]$ (12)

Utilizing the standard derivative for the inverse Weierstrass function, $\frac{\text{d}}{\text{d}u}{\wp }^{-1}\left(u\right)={\left(4u^3-g_{2}u-g_3\right)}^{-1/2}$ , where $u=H_t$ , we obtain:

$\frac{1}{\sqrt{12\kappa }}=\frac{H_{tz}}{\sqrt{4H_t^3-g_2{H}_t-g_3}}$ (13)

Squaring and isolating the field $\kappa$ provides the explicit algebraic relation:

$\kappa \left(x,y,z,t\right)=\frac{4H_t^3-g_2{H}_t-g_3}{12{\left(H_{tz}\right)}^2}$ (14)

This expression reveals that the regularity of $\kappa$ is strictly dependent on the behavior of $H_t$ relative to the roots of the cubic polynomial $P\left(u\right)=4u^3-g_{2}u-g_3$ .

Singularity Analysis and Hölder Exponents

A Hölder singularity of exponent $\alpha$ implies that as $z\to 0$ , $\kappa ~{\left|z\right|}^{\alpha }$ . We investigate the power-law behavior by assuming a local expansion for the time-derivative of $H$ near the origin: $H_t~U_0+c{\left|z\right|}^n$ , where $U_0$ is a root of the cubic polynomial.

Case I: Simple Root ($m=1$ )

If $U_0$ is a simple root such that $P\left(U_0\right)=0$ but $P^{\prime }\left(U_0\right)\ne 0$ , the numerator in the $\kappa$ expression is dominated by the linear term of the Taylor expansion:

$P\left(H_t\right)\approx P^{\prime }\left(U_0\right)\left(H_t-U_0\right)~{\left|z\right|}^n$ (15)

The denominator ${\left(H_{tz}\right)}^2$ scales as:

${\left(nc{z}^{n-1}\right)}^2~{\left|z\right|}^{2n-2}$ (16)

Combining these, the scaling for $\kappa$ is found to be:

$\kappa ~\frac{{\left|z\right|}^n}{{\left|z\right|}^{2n-2}}={\left|z\right|}^{2-n}$ (17)

To achieve the target Hölder exponent $\alpha =1/3$ , we solve $2-n=1/3$ , yielding $n=5/3$ .

Case II: Triple Root (Cusp Geometry, $m=3$ )

If the invariants are chosen such that the cubic has a triple root (which occurs when $g_2=g_3=0$ at $u=0$ ), the curve possesses a cusp singularity. In this regime:

$P\left(H_t\right)~{\left(H_t\right)}^3~{\left|z\right|}^{3n}$ (18)

The resulting $\kappa$ scaling becomes:

$\kappa ~\frac{{\left|z\right|}^{3n}}{{\left|z\right|}^{2n-2}}={\left|z\right|}^{n+2}$ (19)

To obtain $n+2=1/3$ , one would require $n=-5/3$ , implying a divergence in $H_t$ . Thus, the observed $\alpha =1/3$ regularity is mathematically consistent with a simple root geometry at the origin.

2.1.6. Reconstruction of Asymmetric Profiles

Numerical observations indicate a significant asymmetry in the $\kappa$ profile. To reproduce this, we introduce a piecewise coefficient for the spatial gradient of $H$ . We select $g_2=1$ and $g_3=0$ , establishing a simple root at $U_0=0$ . The potential function $H$ is constructed such that:

$\frac{\partial H}{\partial t}=\{\begin{array}{ll}{C}_L{\left(-z\right)}^{5/3}\hfill & z<0\hfill \\ C_R{z}^{5/3}\hfill & z>0\hfill \end{array}$ (20)

Substituting this into the derived formula for $\kappa$ , the leading-order behavior near the origin is:

$\kappa \approx \frac{-g_2{H}_t}{12H_{tz}^2}=\frac{-C{z}^{5/3}}{12{\left(\frac{5}{3}C{z}^{2/3}\right)}^2}\propto z^{1/3}$ (21)

The asymmetry is replicated by setting $C_L\gg C_R$ (e.g., $C_L=5.0$ , $C_R=0.5$ ). This ensures that for $z<0$ , the field enters the singularity with a higher amplitude, while for $z>0$ , the decay is significantly more gradual. The analysis demonstrates that a Hölder continuity of 1/3 in the $\kappa$ field is an emergent property of the 5/3 power-law growth of the underlying potential $H$ when evaluated near a simple root of the Weierstrass cubic. The inverse elliptic transformation acts as a non-linear operator that translates smooth fractional growth into singular geometric profiles.

2.1.7. Leading-Order Analysis at the Origin

As $z\to 0^{+}$ , we assume $H_t\to 0$ (approaching a root of $P\left(u\right)$ ) and $e^{-z}\to 1$ . Neglecting higher-order cubic terms, the ODE reduces to:

$H_{tz}^2\approx \frac{-g_2{H}_t}{12z^{1/3}}\Rightarrow \frac{\text{d}{H}_t}{\sqrt{H_t}}\approx \sqrt{\frac{-g_2}{12}}{z}^{-1/6}\text{d}z$ (22)

Integrating both sides, we obtain the scaling law for the potential:

$2\sqrt{H_t}\approx \sqrt{\frac{-g_2}{12}}\left(\frac{6}{5}{z}^{5/6}\right)\Rightarrow H_t\left(z\right)~\mathcal{C}{z}^{5/3}$ (23)

This confirms that the $z^{1/3}$ Hölder singularity in the geometric field $\kappa$ is uniquely generated by a $z^{5/3}$ scaling of the underlying potential $H_t$ .

2.1.8. Asymptotic Positive Saturation in the Far-Field Limit

In the context of the $H$ -solution derived from the inverse Weierstrass mapping, the phenomenon of asymptotic positive saturation represents a critical topological transition where the geometric field $\kappa$ stabilizes into a steady-state horizontal asymptote. Analytically, for the field to satisfy ${\lim }_{z\to \infty }\kappa \left(z\right)=L$ , the governing relation requires a specific equilibrium between the cubic potential and its spatial evolution:

$12L\approx \frac{4H_t^3-g_2{H}_t-g_3}{H_{tz}^2},\text{as}z\to \infty$ (24)

This condition implies that as $H_t$ diverges or approaches a non-singular trajectory, the growth of the spatial gradient $H_{tz}$ must be strictly slaved to the square root of the Weierstrass cubic, namely $H_{tz}~\sqrt{P\left(H_t\right)}$ . Physically, this saturation describes a localized regularization where the singular dissipation energy-characteristic of the $C^{1/3}$ Hölder singularity at the origin is dissipated into a uniform laminar geometry. This far-field stabilization prevents the re-emergence of spurious singularities and ensures the global boundedness of the regularized flow field.

2.1.9. Analysis of the LambertW Identity with a General $\kappa \left(x,y,z,t\right)$

We aim to find $\kappa$ with a different $H$ in particular the derivative of the LambertW function wrt to $t$ . We consider the equation upon differentiating the immediate above expression wrt to t and setting equal to the derivative of the $W$ function solution in terms of affine functions,

$\wp \left({\displaystyle \int \frac{1}{6}\frac{2^{2/3}\sqrt{3}}{\sqrt{\kappa \left(x,y,z,t\right)2^{2/3}}}\text{d}z}+f_1\left(x,y,t\right);g_2,g_3\right)=\frac{\partial }{\partial t}W\left(-e^{x-y+2z-t-1}\right)$ (25)

Setting equal to the derivative of $W$ wrt to t is justified if Equations (5) - (7) gives a $W$ solution which we do infact prove leads to Equation (65). Now, , where $G$ is the partial of $W$ wrt to t. Next, we invert the WeierstrassP function giving an elliptic integral, so that

${\displaystyle \int \frac{1}{6}\frac{2^{2/3}\sqrt{3}}{\sqrt{\kappa \left(x,y,z,t\right)2^{2/3}}}\text{d}z}+f_1\left(x,y,t\right)={\wp }^{-1}\left(\frac{\partial }{\partial t}W\left(-e^{x-y+2z-t-1}\right)\right)$ (26)

where ${\wp }^{-1}(\cdot )$ is an elliptic integral.

Differentiate w. r. t. $z$

Since $f_1$ is independent of $z$ , differentiating both sides with respect to $z$ gives

$\frac{\partial }{\partial z} {\wp }^{-1}\left(G\left(x,y,z,t\right)\right)=\frac{\frac{\partial G}{\partial z}}{{\wp }^{\prime }\left({\wp }^{-1}\left(G\left(x,y,z,t\right)\right)\right)}.$

where $G\left(x,y,z,t\right)=\frac{\partial }{\partial t}W\left(-e^{x-y+2z-t-1}\right)$ . We start with the equation:

${\displaystyle \int \frac{1}{6}\frac{2^{2/3}\sqrt{3}}{\sqrt{\kappa \left(x,y,z,t\right)2^{2/3}}}\text{d}z}+f_1\left(x,y,t\right)={\wp }^{-1}\left(\frac{\partial }{\partial t}W\left(-e^{x-y+2z-t-1}\right)\right),$ (27)

where $W$ is the Lambert W function. Differentiating both sides with respect to $z$ , we get:

$\frac{1}{6}\frac{2^{2/3}\sqrt{3}}{\sqrt{\kappa \left(x,y,z,t\right)2^{2/3}}}=\frac{\partial }{\partial z}\left[{\wp }^{-1}\left(\frac{\partial }{\partial t}W\left(-e^{x-y+2z-t-1}\right)\right)\right].$ (28)

Solving for $\kappa \left(x,y,z,t\right)$ , we square both sides:

$\kappa \left(x,y,z,t\right)=\frac{2^{2/3}3}{36}{\left[\frac{\partial }{\partial z}{\wp }^{-1}\left(\frac{\partial }{\partial t}W\left(-e^{x-y+2z-t-1}\right)\right)\right]}^{-2}.$ (29)

Thus, the solution for $\kappa$ in terms of $x,y,z,t$ is:

$\kappa \left(x,y,z,t\right)=\frac{2^{2/3}3}{36}{\left[\frac{\partial }{\partial z}{\wp }^{-1}\left(\frac{\partial }{\partial t}W\left(-e^{x-y+2z-t-1}\right)\right)\right]}^{-2}.$ (30)

2.1.10. Behavior of $\kappa$ under Large Lambertw Derivatives

If the spatial derivative of the LambertW term is large, then the coefficient $\kappa \left(x,y,z,t\right)$ becomes small. More precisely,

$\left|{\partial }_{z}W\right|\text{large}\Rightarrow \kappa \text{small}.$ (31)

In fact, $\kappa$ decays quadratically with respect to the derivative:

$\kappa ~\frac{1}{{\left({\partial }_z\left(\text{LambertW term}\right)\right)}^2}$ (32)

In Equations (5) - (7), there is a $\rho \left(V_z+\sqrt{\kappa 1}\right)\frac{{\partial }^2{V}_z}{\partial z^2}$ term which as a result of $\rho$ and steep partial gradients in $z$ based on the form of the $\wp$ function solution, dominates the term $\left(V_z+\sqrt{\kappa 1}\right){\nabla }^2{V}_z$ that is $\rho \left(V_z+\sqrt{\kappa 1}\right)\frac{{\partial }^2{V}_z}{\partial z^2}\gg \left(V_z+C\kappa 1\right){\nabla }^2{V}_z$ .

2.1.11. Exact Structure of $\kappa$

From the transformed equation, $\kappa$ takes the form

$\kappa \left(x,y,z,t\right)=C^{\ast }{\left({\partial }_z{\wp }^{-1}\left({\partial }_{t}W\left(-e^{x-y+2z-t-1}\right)\right)\right)}^{-2},$ (33)

where $C^{\text{*}}>0$ is a constant depending only on fixed parameters.

Thus, $\kappa$ is inversely proportional to the square of a $z$ -derivative.

$\Phi \left(x,y,z,t\right)={\wp }^{-1}\text{}\left(\frac{\partial }{\partial t}W\left(-e^{x-y+2z-t-1}\right)\right),$

where $\wp$ is the Weierstrass elliptic function with invariants $\left(g_2,g_3\right)$ , ${\wp }^{-1}$ denotes a local inverse branch, and $W$ is the principal Lambert $W$ -function. Define,

$\theta :=x-y+2z-t-1.$

We compute

${\partial }_z\text{Φ}\left(x,y,z,t\right),$

$\frac{\text{d}}{\text{d}w} {\wp }^{-1}\left(w\right)=\frac{1}{{\wp }^{\prime }\left({\wp }^{-1}\left(w\right)\right)},{\wp }^{\prime }{\left(z\right)}^2=4\wp {\left(z\right)}^3-g_2\wp \left(z\right)-g_3.$

Thus,

${\partial }_z\Phi =\frac{{\partial }_z\left({\partial }_{t}W\left(-e^{\theta }\right)\right)}{{\wp }^{\prime }\left({\wp }^{-1}\left({\partial }_{t}W\left(-e^{\theta }\right)\right)\right)}.$

Using

$\frac{\text{d}}{\text{d}s}W\left(s\right)=\frac{W\left(s\right)}{s\left(1+W\left(s\right)\right)},$

we obtain

${\partial }_{t}W\left(-e^{\theta }\right)=\frac{W\left(-e^{\theta }\right)}{1+W\left(-e^{\theta }\right)}.$

Differentiating with respect to $z$ ,

${\partial }_z{\partial }_{t}W\left(-e^{\theta }\right)=2{\partial }_{\theta }\left(\frac{W\left(-e^{\theta }\right)}{1+W\left(-e^{\theta }\right)}\right)=\frac{2W\left(-e^{\theta }\right)}{{\left(1+W\left(-e^{\theta }\right)\right)}^3}.$

Therefore,

$\begin{array}{l}{\partial }_z{\wp }^{-1}\left({\partial }_{t}W\left(-e^{x-y+2z-t-1}\right)\right)\\ =\frac{2W\left(-e^{\theta }\right)}{{\left(1+W\left(-e^{\theta }\right)\right)}^3{\wp }^{\prime }\left({\wp }^{-1}\left({\partial }_{t}W\left(-e^{\theta }\right)\right)\right)},\theta =x-y+2z-t-1.\end{array}$

Using the Weierstrass identity

${\wp }^{\prime }\left(u\right)=\sqrt{4\wp {\left(u\right)}^3-g_2\wp \left(u\right)-g_3},$

we may write equivalently

${\partial }_z\Phi =\frac{2W\left(-e^{\theta }\right)}{{\left(1+W\left(-e^{\theta }\right)\right)}^3\sqrt{4{\left({\partial }_{t}W\left(-e^{\theta }\right)\right)}^3-g_2\left({\partial }_{t}W\left(-e^{\theta }\right)\right)-g_3}}.$

2.1.12. Inverse of the Weierstrass $\wp$ -Function at Infinity

We consider the Weierstrass elliptic function $\wp \left(u;g_2,g_3\right)$ associated with the period lattice

$\text{Λ}=\left\{m{\omega }_1+n{\omega }_2:m,n\in \mathbb{Z}\right\}.$

Behavior near the pole

The function $\wp \left(u\right)$ has a double pole at $u=0$ (and at every lattice point $u\in \text{Λ}$ ). Its Laurent expansion near the origin is given by

$\wp \left(u\right)=\frac{1}{u^2}+O\left(u^2\right),u\to 0.$ (34)

Consequently,

$\wp \left(u\right)\to \infty \text{as}u\to 0\left(\mathrm{mod}\Lambda \right).$ (35)

Inverse of the $\wp$ -function

The inverse ${\wp }^{-1}$ is a multivalued function, since $\wp$ is even and doubly periodic:

$\wp \left(u\right)=\wp \left(-u\right).$

A local branch of the inverse is defined via the elliptic integral,

${\wp }^{-1}\left(z\right)={\displaystyle {\int }_{\infty }^z\frac{\text{d}s}{\sqrt{4s^3-g_{2}s-g_3}}}.$ (36)

Inverse at infinity

By definition of the elliptic integral, we obtain

${\wp }^{-1}\left(\infty \right)=\left[0\right].$ (37)

More precisely, since $\wp$ is doubly periodic,

${\wp }^{-1}\left(\infty \right)=\left\{ 0+\lambda :\lambda \in \text{Λ} \right\}.$ (38)

${\wp }^{-1}\left(\infty \right)=0\left(\mathrm{mod}\Lambda \right).$ (39)

2.1.13. Singular Limits

If a function $G\left(x,y,z,t\right)$ satisfies $\left|G\left(x,y,z,t\right)\right|\to \infty$ , then

${\wp }^{-1}\left(G\left(x,y,z,t\right)\right)\to 0,$

and therefore

$\frac{\partial }{\partial z}{\wp }^{-1}\left(G\left(x,y,z,t\right)\right)\to \infty .$

In applications where a coefficient $\kappa \left(x,y,z,t\right)$ depends inversely on this derivative squared, this implies

$\kappa \left(x,y,z,t\right)\to 0,$

revealing a singular-regular transition governed by the elliptic structure of the Weierstrass $\wp$ -function.

Lemma 1 (Elliptic degeneracy of $\kappa$ ) Let be a Weierstrass elliptic function with lattice $\Lambda$ , and let

$\kappa \left(x,y,z,t\right)=C{\left[{\partial }_z{\wp }^{-1}\left(G\left(x,y,z,t\right)\right)\right]}^{-2},C>0,$

where $G\left(x,y,z,t\right)$ is a smooth function satisfying

$\left|G\left(x,y,z,t\right)\right|\to \infty$ along infinitely many spatial curves,

then

$\kappa \left(x,y,z,t\right)\to 0$ at infinitely many spatial points.

Proof. The Weierstrass function $\wp \left(u\right)$ has a double pole at every lattice point $\lambda \in \text{Λ}$ , so that

$\wp \left(u\right)\to \infty \text{as}u\to \lambda .$

Consequently,

${\wp }^{-1}\left(z\right)\to \lambda \text{as}z\to \infty ,$

for each $\lambda \in \text{Λ}$ , implying that ${\wp }^{-1}$ has infinitely many branches.

Differentiating the inverse yields

$\frac{\text{d}}{\text{d}z}{\wp }^{-1}\left(z\right)=\frac{1}{{\wp }^{\prime }\left({\wp }^{-1}\left(z\right)\right)}.$

Near a lattice point,

${\wp }^{\prime }\left(u\right)~-\frac{2}{{\left(u-\lambda \right)}^3},$

so the derivative of the inverse diverges as $z\to \infty$ :

$\left|{\partial }_z{\wp }^{-1}\left(G\right)\right|\to \infty .$

Hence

$\kappa \left(x,y,z,t\right)\propto {\left({\partial }_z{\wp }^{-1}\left(G\right)\right)}^{-2}\to 0,$

at each such point. Since $\text{Λ}$ is infinite, this occurs at infinitely many spatial locations.

2.1.14. Elements of Equations (5) - (7) and Expressing Laplacians over $ϵ$ Balls

It works out that if ${\wp }^{\prime }=\kappa$ (if $\kappa =0$ then we solve for the zeros $e_i$ as usual), where ($\wp$ is the WeierstrassP function) this makes the derivative wrt to $z$ , ${\partial }_z\Phi$ , as the parameter scaling $\delta$ . In Equations (5) - (7), the term $M\left(ϵ\right)$ defined in [3] appears when we write the Laplacian as an integral on an epsilon ball. Here in [2] and [3] with references within from [Rumer and Fet], $M\left(ϵ\right)=10/{ϵ}^2$ which is infinite when $ϵ\to 0$ . Note that $M\left(ϵ\right)\cdot \kappa =O\left(1\right)$ hence $M\left(ϵ\right)={\kappa }^{-1}$ simplifying the Laplacian terms. The important observation is that we have an ${\left(\frac{\partial V_z}{\partial t}\right)}^3$ term in the PDE compared to an expression of operators multiplied by ${\left(\frac{\partial V_z}{\partial t}\right)}^2$ term. Next, $\delta =-\frac{1}{\sqrt{\kappa 1}}$ and ${\left(V_z+\sqrt{\kappa 1}\right)}^2{\left(\frac{{\partial }^2{V}_z}{\partial z\partial t}\right)}^2\delta =\frac{1}{3}\frac{1}{\kappa }{\left(\frac{\partial V_z}{\partial t}\right)}^3\delta$ as used in Equations (5) - (7) to obtain Equation (65). This has a factor of $\frac{1}{\sqrt{\kappa 1}}\frac{1}{\kappa }$ in the PDE. Multiply by this factor throughout the PDE. Most other occurrences in the PDE depend on $\delta$ . Now, Figure 2(a) and Figure 2(b) have the plots of $\wp$ versus $U={\displaystyle \int \frac{1}{6}\frac{2^{2/3}\sqrt{3}}{\sqrt{\kappa \left(x,y,z,t\right)2^{2/3}}}\text{d}z}+f_1\left(x,y,t\right)$ .

Figure 2. Plot of $\wp$ versus $U={\displaystyle \int \frac{1}{6}\frac{2^{2/3}\sqrt{3}}{\sqrt{\kappa \left(x,y,z,t\right)2^{2/3}}}\text{d}z}+f_1\left(x,y,t\right)$ for $m=100$ associated with the two invariants where $g_1=3m^2$ and $g_2=m^3$ .

2.1.15. Weierstrass $\wp$ Invariants, Plots, Inverse $\wp$ Functions Local Square-Root Behavior and Singuarities Due to Elliptic Degeneracy Prescribed at Critical Values

Here, in the two plots, Figure 2(a) and Figure 2(b), it is seen that with $m$ large the $P$ functions become uniformly narrower. As periods shrink, fixed spatial variations intersect more branch structure. So, we can trust that our function $\kappa$ will be accurate and the PDE of Equations (5) - (7) will be valid with increasing spatial density and variation. The multivalued (branched) structure of the WeierstrassP function allows us to choose a local inverse branch on which $\kappa \left(x,y,z,t\right)$ can be treated as effectively independent of the spatial variables to leading order, depending only on how the solution approaches a critical value of $\wp$ . Of course, time $t$ and its variation are still part of the LambertW structure.

Lemma 2. Let $\wp \left(u;g_2,g_3\right)$ be the Weierstrass elliptic function with invariants $g_2,g_3\in ℝ$ satisfying $g_2^3\ne 27g_3^2$ (non-degenerate case). Let ${\wp }^{-1}$ denote a local inverse branch of $\wp$ .

Define

$U\left(x,y,z,t\right):={\wp }^{-1}\left({\partial }_{t}W\left(-e^{x-y+2z-t-1}\right)\right),$

where $W$ is the principal branch of the Lambert $W$ function. Then, the spatial derivative ${\partial }_{z}U$ blows up if and only if

${\partial }_{t}W\left(-e^{x-y+2z-t-1}\right)\in \left\{e_1,e_2,e_3\right\},$

where $e_1,e_2,e_3$ are the three distinct real or complex roots of

$4w^3-g_{2}w-g_3=0.$

Proof. The Weierstrass $\wp$ -function satisfies the classical differential identity

${\left({\wp }^{\prime }\left(u\right)\right)}^2=4\wp {\left(u\right)}^3-g_2\wp \left(u\right)-g_3.$ (40)

Hence, ${\wp }^{\prime }\left(u\right)=0$ if and only if $\wp \left(u\right)$ equals one of the roots $e_1,e_2,e_3$ of the cubic polynomial on the right-hand side of (40).

Since ${\wp }^{-1}$ is locally the inverse of $\wp$ , the inverse function theorem implies

$\frac{\text{d}}{\text{d}w}{\wp }^{-1}\left(w\right)=\frac{1}{{\wp }^{\prime }\left(u\right)},w=\wp \left(u\right),$

whenever ${\wp }^{\prime }\left(u\right)\ne 0$ .

Differentiating $U$ with respect to $z$ yields

${\partial }_{z}U=\frac{1}{{\wp }^{\prime }\left(U\right)}{\partial }_z\left[{\partial }_{t}W\left(-e^{x-y+2z-t-1}\right)\right].$

The Lambert $W$ term is not smooth for all finite arguments and therefore ${\partial }_z{\partial }_{t}W(\cdot )$ can be infinite. Also, ${\partial }_{z}U$ can blow up when the denominator ${\wp }^{\prime }\left(U\right)$ vanishes.

By (40), this occurs if and only if

$\wp \left(U\right)\in \left\{e_1,e_2,e_3\right\},$

which, by the definition of $U$ , is equivalent to

${\partial }_{t}W\left(-e^{x-y+2z-t-1}\right)\in \left\{e_1,e_2,e_3\right\}.$

This establishes the equivalence and completes the proof.

Remark 1. Near such a critical value $e_i$ , the inverse function exhibits the local square-root behavior

${\wp }^{-1}\left(w\right)=u_i+\mathcal{O}\left(\sqrt{w-e_i}\right),$

so that ${\partial }_{z}U$ diverges algebraically. The singularity is therefore non-oscillatory and arises from elliptic degeneracy rather than from the Lambert $W$ forcing itself.

Lemma 3. Let $\wp \left(u\right)$ denote the Weierstrass $\wp$ -function with invariants $g_2,g_3\ne 0$ , and let $e_1,e_2,e_3$ be the roots of

$4w^3-g_{2}w-g_3=0.$

Suppose $e_1=-1$ (or any fixed root). Define the function

$u\left(x,y,z,t\right)={\wp }^{-1}\left({\partial }_{t}W\left(-e^{x-y+2z-t-1}\right)\right),$

where ${\wp }^{-1}$ is a local inverse branch of $\wp$ . Then, ${\partial }_{z}u$ blows up whenever

${\partial }_{t}W\left(-e^{x-y+2z-t-1}\right)=e_1=-1.$

Proof. 1) By definition, the derivative of the inverse function satisfies

$\frac{\partial }{\partial w}{\wp }^{-1}\left(w\right)=\frac{1}{{\wp }^{\prime }\left(u\right)},w=\wp \left(u\right).$

2) The critical points of $\wp$ are given by ${\wp }^{\prime }\left(u\right)=0$ . These occur precisely at $u$ such that $\wp \left(u\right)=e_i$ , $i=1,2,3$ .

3) Substituting $w={\partial }_{t}W\left(-e^{x-y+2z-t-1}\right)$ into the chain rule for ${\partial }_z$ , we have

${\partial }_{z}u=\frac{1}{{\wp }^{\prime }\left(u\right)} {\partial }_z\left[{\partial }_{t}W\left(-e^{x-y+2z-t-1}\right)\right].$

4) If ${\partial }_{t}W\left(-e^{x-y+2z-t-1}\right)=e_1=-1$ , then ${\wp }^{\prime }\left(u\right)=0$ and hence

${\partial }_{z}u=\frac{{\partial }_z\left({\partial }_{t}W\right)}{0}\to \infty .$

5) Therefore, the derivative along $z$ blows up at the set of points $\left(x,y,z,t\right)$ where the Lambert $W$ -driven forcing equals $e_1=-1$ .

2.1.16. Local Independence of Weierstrass Critical Points from Spacetime Variables

In this section, we prove that there is local independence of Weierstrass critical points from spacetime variables depending only on the elliptic invariants $\left(g_2,g_3\right)$ . We have the following theorem.

Theorem 3 (Independence of the spacetime parametrization) Let $\wp \left(u;g_2,g_3\right)$ be the Weierstrass elliptic function with invariants

$g_2,g_3\in ℝ,\Delta :=g_2^3-27g_3^2\ne 0,$

and let $e_1,e_2,e_3$ be the distinct real roots of

$4w^3-g_{2}w-g_3=0.$

Let $u_i\in ℂ$ be a half-period satisfying

$\wp \left(u_i\right)=e_i,{\wp }^{\prime }\left(u_i\right)=0.$

Let $w:\Omega \subset {ℝ}^3\times ℝ\to ℝ$ be a $C^1$ function, and define

$u\left(x,y,z,t\right):={\wp }^{-1}\text{}\left(w\left(x,y,z,t\right)\right),$

using a single-valued local inverse branch defined in a neighborhood of $e_i$ .

Then, there exists a neighborhood $U\subset \Omega$ such that:

${\wp }^{\prime }\left(u\left(x,y,z,t\right)\right)=0$ if and only if $w\left(x,y,z,t\right)=e_i$ for all $\left(x,y,z,t\right)\in U$ .

The condition ${\wp }^{\prime }\left(u\right)=0$ is independent of the spacetime variables $\left(x,y,z,t\right)$ and depends only on the elliptic invariants $\left(g_2,g_3\right)$ .

The derivative of the inverse map satisfies the local singular behavior

$\left|{\nabla }_{x,y,z,t}u\left(x,y,z,t\right)\right|~\frac{\left|{\nabla }_{x,y,z,t}w\left(x,y,z,t\right)\right|}{\sqrt{ \left|w\left(x,y,z,t\right)-e_i\right|}},$

and therefore diverges as $w\left(x,y,z,t\right)\to e_i$ .

Proof. Since $\text{Δ}\ne 0$ , the roots $e_1,e_2,e_3$ are distinct and the Weierstrass function satisfies the algebraic identity

${\left({\wp }^{\prime }\left(u\right)\right)}^2=4\left(\wp \left(u\right)-e_1\right)\left(\wp \left(u\right)-e_2\right)\left(\wp \left(u\right)-e_3\right).$

Thus, ${\wp }^{\prime }\left(u\right)=0$ if and only if $\wp \left(u\right)=e_i$ for some $i\in \left\{1,2,3\right\}$ .

Fix such an index $i$ . By analyticity of $\wp$ , there exists a neighborhood $V\subset ℂ$ of $u_i$ on which $\wp$ is locally invertible except at $u_i$ . The inverse map ${\wp }^{-1}$ is well-defined on a neighborhood of $e_i$ excluding the branch point.

Let $w\left(x,y,z,t\right)$ take values in this neighborhood and define

$u\left(x,y,z,t\right):={\wp }^{-1}\left(w\left(x,y,z,t\right)\right).$

Then

${\wp }^{\prime }\left(u\left(x,y,z,t\right)\right)=0\iff \wp \left(u\left(x,y,z,t\right)\right)=e_i\iff w\left(x,y,z,t\right)=e_i,$

which proves (1).

Statement (2) follows since the equation ${\wp }^{\prime }\left(u\right)=0$ depends only on the elliptic curve parameters $\left(g_2,g_3\right)$ and the variable $u$ , and not on how $u$ is parametrized by spacetime variables.

To prove (3), expand $\wp$ in a Taylor series at $u_i$ :

$\wp \left(u\right)=e_i+\frac{1}{2}{\wp }^{″}\left(u_i\right){\left(u-u_i\right)}^2+O\left({\left(u-u_i\right)}^3\right).$

Inverting this relation yields

$u-u_i~\sqrt{\frac{2\left(w-e_i\right)}{{\wp }^{″}\left(u_i\right)}}.$

Differentiation with respect to spacetime variables gives

$\nabla u~\frac{\nabla w}{\sqrt{w-e_i}},$

which diverges as $w\to e_i$ .

Remark 2. The theorem shows that any PDE solution whose structure includes a composition with ${\wp }^{-1}$ inherits a geometrically determined singularity whenever the composed quantity reaches a critical elliptic value $e_i$ . The singularity mechanism is independent of the spacetime parametrization and therefore invariant under coordinate changes.

2.2. The Pivot Function in Equation (6) [1] and the Question of Finite Time Blowup

As to the discussion of the pivot function $a\nabla \cdot a$ , this expression is periodic in ${\mathbb{T}}^3$ and hence it’s integral is zero. Also, inequalities have been used to show that this is confirmed and that the sum of the main operators to PNS system become zero which is consistent with the previous PDE given by Equations (5) - (7). It is noteworthy to mention that the addition of the pivot function in Equation (6) of reference [1] demands that we finally return to Equation (6) there to solve an $L_1=L_2$ PDE. This is the final step in the solution approach. Hence, it will be seen that a solution of the form for $u_x{u}_y{u}_z$ is given as:

$u_x{u}_y{u}_z=-\frac{W+1}{W-1}$ (41)

where $W$ is the LambertW function in particular:

$W=\text{LambertW}\left(-\text{exp}\left(-\Xi -1\right)\right)+C$

where $\text{Ξ}$ is an affine function $\text{Ξ}\left(x,t\right)$ . The interesting development of such a proven result is that it entails both singular and no finite time blowup on the surface of an imaginary surface, the sphere. This was obtained by solving Equation (6) in [1] directly which is shown in Appendix A. The solution of the PNS system is not smooth and has no finite time blowup since $W=1$ is not realized if a real solution branch is required as we are looking at explicitly real valued solutions to the real Navier-Stokes equations.

2.3. Iterating and Shifting $W$ Function Solutions by Fixed Point Methods and Application of Nth Order Compositions

Zeros of compositions of $W$ function solutions on the sphere inscribed in the torus need to be shifted outwards by applying higher order compositions in order to leave smooth solutions behind. It has been shown in [2] that using infinite order nth compositions via fixed point methods that the points on a sphere for a given torus can move outwards when we consider expanding tori: $T_{\alpha }^3={\left[-\alpha ,\alpha \right]}^3$ . The essence of such idea has been presented in [2]. It is based on the following idea, Let

$G_n\left(x\right)=\Psi \left(F_n\left(x\right)\right),$

where $\Psi$ is singular at 0 (for example $\Psi \left(s\right)=1/s$ , $\log s$ , or $W\left(s\right)$ at a branch point).

If $F_n\left(x_n\right)=0$ , then $G_n$ has a singularity at $x=x_n$ . Thus, zeros of $F_n$ correspond to singularities of $G_n$ .

Consider iterated expressions of the form

$F_n\left(x\right)={\Phi }^{\left(n\right)}\left(x\right)+1,\Phi \left(x\right)=W\left(-\left(1+2^n\right)W\left(-e^{-x-1}\right)\right).$

At a zero $x_n$ ,

$F_n\left(x_n\right)=0\iff {\Phi }^{\left(n\right)}\left(x_n\right)=-1.$

For the Lambert $W$ function,

$W\left(z\right)=-1$ if and only if $z=-e^{-1}$ ,

which is a branch-point singularity of the inverse mapping.

Hence:

$\begin{array}{l}\text{zero of }{F}_n\iff \text{branch point of the outer }W\\ \iff \text{singularity of the composed expression}.\end{array}$

The locations $x_n$ of the zeros satisfy

$x_n~\text{log}\left(1+2^n\right)\to +\infty .$

2.4. Derivation of the Periodic Velocity $V_z$ via Frequency Scaling (Lambert $W$ Included)

2.4.1. Vertical Velocity with Space-Time Dependent Lambert $W$

We define the vertical velocity using the Weierstrass kernel. The argument now explicitly includes the space-time dependent Lambert $W$ evolution:

$V_z\left(x,y,z,t\right)=\mathcal{C}{\displaystyle \int \wp \left(\sqrt{\omega }{u}_z\left(x,y,z,t\right)+f_1\left(x,y,t\right)\right)\text{d}t}+f_2\left(x,y,z\right),$ (42)

with the full space-time Lambert W contribution in $u_z$ :

$u_z\left(x,y,z,t\right)={\displaystyle \int \frac{2^{2/3}\sqrt{3}}{6\sqrt{2^{2/3}\kappa \left(x,y,z,t\right)}}\text{d}z}$ (43)

where $V_z=W\left(x,y,z,t\right)$ , $W$ is the Lambert $W$ evolution (nonlinear in space and time), and $f_1\left(x,y,t\right)$ will be linear in space and time to offset exactly one of the $\zeta$ functions. Now, $u_z$ fully encodes the space-time nonlinear Lambert W dynamics.

2.4.2. Small $f_1$ Linearization

Assuming $f_1\ll \sqrt{\omega }{u}_z$ , the addition formula can be linearized:

$\wp \left(\sqrt{\omega }{u}_z+f_1\right)\approx \wp \left(\sqrt{\omega }{u}_z\right)+f_1{\wp }^{\prime }\left(\sqrt{\omega }{u}_z\right)+\mathcal{O}\left(f_1^2\right)$ (44)

The first term contains the full space-time Lambert $W$ evolution. The second term is small, time-dependent, and captures the singular drift from $f_1$ .

We want to justify the linearization:

$\wp \left(\sqrt{\omega }{u}_z+f_1\right)\approx \wp \left(\sqrt{\omega }{u}_z\right)+f_1{\wp }^{\prime }\left(\sqrt{\omega }{u}_z\right)+\mathcal{O}\left(f_1^2\right)$

The exact addition formula for the Weierstrass $\wp$ function is

$\wp \left(u+v\right)=-\wp \left(u\right)-\wp \left(v\right)+\frac{1}{4}{\left(\frac{{\wp }^{\prime }\left(u\right)-{\wp }^{\prime }\left(v\right)}{\wp \left(u\right)-\wp \left(v\right)}\right)}^2$

For arbitrary $u,v$ . Here, we set

$u=\sqrt{\omega }{u}_z,v=f_1.$

For small $f_1$ , we can expand $\wp \left(u+f_1\right)$ as a Taylor series around $u$ :

$\wp \left(u+f_1\right)=\wp \left(u\right)+f_1{\wp }^{\prime }\left(u\right)+\frac{f_1^2}{2}{\wp }^{″}\left(u\right)+\cdots$

Neglecting higher-order terms gives the first-order linearization:

$\wp \left(u+f_1\right)\approx \wp \left(u\right)+f_1{\wp }^{\prime }\left(u\right).$

The addition formula also gives

$\wp \left(u+f_1\right)-\wp \left(u\right)=\frac{1}{4}{\left(\frac{{\wp }^{\prime }\left(u\right)-{\wp }^{\prime }\left(f_1\right)}{\wp \left(u\right)-\wp \left(f_1\right)}\right)}^2-\wp \left(u\right)-\wp \left(f_1\right).$

For $f_1\to 0$ , we can expand

$\wp \left(f_1\right)\approx \frac{1}{f_1^2}+\cdots ,{\wp }^{\prime }\left(f_1\right)\approx -\frac{2}{f_1^3}+\cdots$

Substituting into the addition formula and expanding in powers of $f_1$ , the leading term is

$\wp \left(u+f_1\right)-\wp \left(u\right)\approx f_1{\wp }^{\prime }\left(u\right),$

which agrees with the Taylor expansion.

Thus, for small $f_1$ , we have,

$\boxed{\wp \left(\sqrt{\omega }{u}_z+f_1\right)\approx \wp \left(\sqrt{\omega }{u}_z\right)+f_1{\wp }^{\prime }\left(\sqrt{\omega }{u}_z\right)+\mathcal{O}\left(f_1^2\right).}$

2.4.3. Order of the Linearized Term $f_1{\wp }^{\prime }\left(\sqrt{\omega }{u}_z\right)$

Consider the linearized term arising from the small $f_1$ expansion:

$f_1{\wp }^{\prime }\left(\sqrt{\omega }{u}_z\right).$

2.4.4. Definition of $f_1$ and Θ

We take

$f_1\le {\lambda }_1<\infty$

so that $f_1$ is bounded everywhere in the lattice and the derivative ${\wp }^{\prime }\left(\sqrt{\omega }{u}_z\right)$ depends on the nonlinear Lambert $W$ evolution, and is bounded (oscillatory), away from the poles of the $\wp$ function. Here, the $W$ and $\wp$ functions are related at their singularities. What is to be achieved in this is to define a drift function which together with the quasi-periodic WeierstrassZeta function $\zeta$ (as a difference) produce a periodic function with the poles of $\zeta$ included. When a lattice period goes to infinity, the elliptic structure degenerates, and the poles of the WeierstrassP, $\wp$ function recede to infinity, leaving behind a simpler (non-elliptic) function. We have:

• The Weierstrass $\zeta \left(z\right)$ function alone is quasi-periodic, not strictly periodic.

• The linear drift term $-\frac{\eta }{\Omega }\Theta$ cancels the quasi-periodicity, so that the combination

$\zeta \left(\Theta \right)-\frac{\eta }{\Omega }\Theta$

is strictly periodic.

• Even if the lattice periods become very large or tend to infinity, this cancellation ensures that a velocity $V_{\text{target}}$ periodic exits over the relevant spatial and temporal domain.

Along the synchronization boundary, $u_z~\Theta /\sqrt{\omega }$ , so that

${\wp }^{\prime }\left(\sqrt{\omega }{u}_z\right)~\mathcal{O}\left(1\right)\left(\text{bounded amplitude}\right).$

Hence, the product scales as

$f_1{\wp }^{\prime }\left(\sqrt{\omega }{u}_z\right)~O\left(1\right)$

Although $f_1$ is small in magnitude, its time integral carries a factor of $\omega t$ :

${\displaystyle \int \frac{\eta }{\Omega }{f}_1{\wp }^{\prime }\left(\sqrt{\omega }{u}_z\right)\text{d}t}=\frac{\eta }{\Omega }\left(-\omega \right)t=-\frac{\eta }{\Omega }\omega t$ (45)

${\displaystyle \int f_1{\wp }^{\prime }\left(\sqrt{\omega }{u}_z\right)\text{d}\Theta }~\Theta ,$

and can be canceled with the WeierstrassZeta term to ensure periodicity.

The integral becomes:

$\begin{array}{c}I=\mathcal{C}{\displaystyle \int \wp \left(\sqrt{\omega }{u}_z+f_1\right)\text{d}t}\\ \approx \mathcal{C}{\displaystyle \int \wp \left(\sqrt{\omega }{u}_z\left(x,y,z,t\right)\right)\text{d}t}+\mathcal{C}{\displaystyle \int f_1\left(x,y,t\right){\wp }^{\prime }\left(\sqrt{\omega }{u}_z\left(x,y,z,t\right)\right)\text{d}t}\end{array}$ (46)

1) First term: oscillatory contribution including space-time Lambert $W$ .

2) Second term: linear in $f_1$ , small time-dependent singular evolution.

2.4.5. Compensator $f_2$ and Target Velocity

To ensure strict periodicity:

$f_2\left(x,y,z\right)=V_{\text{target}}\left(x,y,z,t\right)-I,$ (47)

with

$V_{\text{target}}\left(x,y,z,t\right)=2^{2/3}\frac{\mathcal{A}}{\omega }\left[\zeta \left(\Theta \right)-\frac{\eta }{\Omega }\Theta \right],\Theta =k\cdot x-\omega t$ (48)

The drift from the small $f_1$ term is absorbed into $f_2$ or into the pressure gauge, leaving $u=\nabla V_z$ strictly periodic.

2.4.6. Temporal Derivative and Frequency Matching

Using $\frac{\text{d}}{\text{d}z}\zeta \left(z\right)=-\wp \left(z\right)$ :

$\frac{\partial V_{\text{target}}}{\partial t}=2^{2/3}\frac{\mathcal{A}}{\omega }\left[-\wp \left(\Theta \right)-\frac{\eta }{\Omega }\right]\left(-\omega \right)$ (49)

$=2^{2/3}\mathcal{A}\wp \left(\Theta \right)+2^{2/3}\mathcal{A}\frac{\eta }{\Omega }$ (50)

$\frac{\partial I}{\partial t}\approx \mathcal{C}\wp \left(\sqrt{\omega }{u}_z\left(x,y,z,t\right)\right)+\mathcal{C}{f}_1{\wp }^{\prime }\left(\sqrt{\omega }{u}_z\left(x,y,z,t\right)\right)$ (51)

$\approx 2^{2/3}\wp \left(\Theta \right)+\text{small}{f}_1\text{correction}$ (52)

The leading term matches perfectly. The small $f_1$ term is explicitly controlled.

2.4.7. Homogeneity Scaling at Synchronization Boundary

At $v\approx \Theta /\sqrt{\omega }$ :

$\wp \left(v\right)=\wp \left(\Theta /\sqrt{\omega }\right)=\omega \wp \left(\text{Θ}\right)$ (53)

This ensures synchronization between lattice frequency and internal nonlinear evolution.

2.4.8. Final Periodic Velocity Expression

$V_z\left(x,y,z,t\right)=\frac{2^{2/3}}{\omega }\left[\zeta \left(\Theta \right)-\frac{\eta }{\Omega }\Theta \right]$ (54)

$\Theta$ now implicitly encodes the full space-time Lambert $W$ evolution in $u_z$ and the small time-dependent singularity in $f_1$ . Solution is bounded, strictly periodic, and fully accounts for space-time nonlinear dynamics.

Observation: By incorporating the Lambert $W$ evolution directly into $u_z$ , the space-time structure of the solution is fully retained, while the small $f_1$ term allows controlled drift cancellation, resulting in a strictly periodic velocity field.

2.4.9. The Mapping Strategy: Lambert to Weierstrass

The proposed solution identifies that the local singular behavior, often described by branches of the Lambert $W$ function, can be regularized by embedding the temporal evolution into the argument of a Weierstrass $\wp$ function.

By setting the internal scaling constant $C^{\text{*}}=\sqrt{\omega }$ , we invoke the Homogeneity Property:

$\wp \left(\frac{z}{\sqrt{\omega }};\frac{g_2}{{\omega }^2},\frac{g_3}{{\omega }^3}\right)=\omega \wp \left(z;g_2,g_3\right)$ (55)

This allows the “speed” of the singular part to be perfectly synchronized with the frequency of the periodic wave, facilitating the cancellation of divergent terms. The Weierstrass Zeta function $\zeta \left(\Theta \right)$ , which serves as the base for the velocity $V_z$ , is quasi-periodic:

$\zeta \left(\Theta +\Omega \right)=\zeta \left(\Theta \right)+\eta$ (56)

To satisfy the Millennium problem requirement of a solution on a compact manifold (the 3-Torus ${\mathbb{T}}^3$ ), we utilize a drift-correction term. The final potential is defined as:

$V_z\left(\Theta \right)=\frac{2^{2/3}}{\omega }\left[\zeta \left(\Theta \right)-\frac{\eta }{\Omega }\Theta \right]$ (57)

As proven in the periodicity analysis, $V_z\left(\Theta +\Omega \right)=V_z\left(\Theta \right)$ . This ensures that the velocity field $u=\nabla V_z$ is strictly periodic and bounded.

2.4.10. Gauge Invariance and Temporal Stability

To match the temporal gradients of the integral $I$ and the target $V_{target}$ , we utilize the gauge symmetry $\Phi \to \Phi -Ct$ . By choosing $C=\frac{\eta \omega }{\Omega }$ , the secular linear growth in time is absorbed into the global pressure gauge. This ensures:

The fluctuations match: ${\partial }_t{V}_{\text{target}}={\partial }_{t}I=\mathcal{C}\wp \left(\Theta \right)$ . The spatial residue $f_2$ is stationary: $\frac{\partial f_2}{\partial t}=0$ . The pressure Poisson equation remains well-posed.

The mapping of the Navier-Stokes velocity to a drift-corrected Weierstrass Zeta lattice provides a candidate for global regularity. It bridges the gap between local non-linear singularities and global periodic stability, ensuring that the velocity field remains smooth and the energy remains finite for all time.

2.5. No Finite Time Blowup of $\left(W+1\right)/\left(W-1\right)$

We analyze whether the expression

$\frac{W+1}{W-1}$ (58)

indicates a finite-time blowup in the Navier-Stokes equations, where

$W=\text{LambertW}\left(-e^{-\Xi -1}\right),\Xi =x-y+2z-\omega t+C.$ (59)

2.5.1. Algebraic Singularity of the Ratio

The rational expression diverges only when

$W\to 1.$ (60)

At this value, the denominator vanishes and the ratio becomes unbounded. Thus, any apparent blowup signaled by this expression would require that the Lambert $W$ function attain the value $W=1$ .

2.5.2. Accessibility of $W=1$

However, the argument of the Lambert $W$ function in this construction is

$z=-e^{-\Xi -1},$ (61)

which is strictly negative for all real values of $\Xi$ . On the real branches of the Lambert $W$ function:

The principal branch $W_0\left(z\right)$ satisfies $W_0\left(z\right)\le 0$ for $z\in \left[-e^{-1},0\right)$ . The lower branch $W_{-1}\left(z\right)$ satisfies $W_{-1}\left(z\right)\le -1$ .

Therefore, the value $W=1$ is not attainable for any real $\text{Ξ}$ in this formulation. Consequently, the algebraic pole at $W=1$ is not dynamically accessible in physical space-time. The genuine singularity of the Lambert $W$ function occurs at the branch point

$W=-1,$ (62)

where the derivative $W^{\prime }\left(z\right)$ diverges. This is the only real singularity of $W$ and represents the location where velocity gradients could, in principle, become unbounded.

Notably, at this point, the ratio evaluates to

$\frac{W+1}{W-1}=\frac{0}{-2}=0,$ (63)

which is completely finite. Thus, the ratio $\frac{W+1}{W-1}$ does not detect the actual singular behavior of the Lambert $W$ function.

In the constructed solution, the Lambert $W$ representation is replaced by a Weierstrass-based velocity of the form

$V_z\left(\Theta \right)=\frac{2^{2/3}}{\omega }\left[\zeta \left(\Theta \right)-\frac{\eta }{\Omega }\Theta \right],\Theta =k\cdot x-\omega t.$ (64)

Under this mapping:

The Lambert branch point $W=-1$ is transported to poles of the Weierstrass Zeta function. These poles occur only at discrete lattice values of $\Theta$ . In physical space-time, they correspond to isolated, measure-zero sets.

The linear drift correction removes quasi-periodicity and prevents secular growth in time.

2.5.3. Implications for Finite-Time Blowup

Finite-time blowup in the Navier-Stokes equations requires divergence of a physical norm (such as velocity, vorticity, or enstrophy) over a nonzero spatial region at finite time. In the present construction:

The only reachable singularities are isolated lattice poles. These singularities are integrable and spatially localized. No velocity or energy norm diverges at finite time.

Therefore, the presence of the ratio $\frac{W+1}{W-1}$ does not imply finite-time blowup.

The apparent pole at $W=1$ in the expression $\frac{W+1}{W-1}$ is an algebraic artifact and is not reachable in the real-valued solution. The true Lambert $W$ singularity at $W=-1$ is mapped, via the Weierstrass formulation, to isolated lattice poles that do not produce finite-time blowup in the Navier-Stokes equations.

Finite-time blowup is a property of norms and dynamical accessibility, not of algebraic expressions alone.

2.6. Corollary (No Free Temporal Perturbations and Periodic Solutions to PNS)

Any additional term $H\left(x,y,z,t\right)$ satisfying

$\frac{\partial H}{\partial t}\ne 0$

cannot be absorbed into $f_2$ . It must either:

Be expressible as a gauge term $C\left(t\right)$ , or destroy stationarity and periodicity. The Lambert $W$ dynamics fixes the internal clock. The Weierstrass homogeneity fixes the scaling. The quasi-periodicity fixes the drift term. The stationarity condition eliminates all remaining freedom. Periodic regularization forces a unique potential.

Once the internal singular dynamics and the external lattice frequency are synchronized, the Navier-Stokes velocity admits exactly one stationary compensator and no spatiotemporal freedom remains. However, if $H\left(x,y,z,t\right)$ is strictly periodic in $t$ with the same period $T=2\pi /\omega$ as the lattice, then it can still be added without destroying the overall periodicity. That is because even though ${\partial }_{t}H\ne 0$ locally, the integral over one period cancels out, and the velocity field remains periodic:

${\tilde{V}}_z\left(x,y,z,t\right)=V_z\left(x,y,z,t\right)+H\left(x,y,z,t\right),H\left(x,y,z,t+T\right)=H\left(x,y,z,t\right)$

Then, by construction:

$\begin{array}{c}{\tilde{V}}_z\left(x,y,z,t+T\right)=V_z\left(x,y,z,t\right)+H\left(x,y,z,t+T\right)\\ =V_z\left(x,y,z,t\right)+H\left(x,y,z,t\right)\\ ={\tilde{V}}_z\left(x,y,z,t\right),\end{array}$

so the total velocity is still periodic, even though $H$ contributes a nonzero derivative ${\partial }_{t}H$ .

The exact drift cancellation mechanism for the Weierstrass Zeta function only applies to the base $V_z=I+f_2$ . Adding $H$ does not interfere with the cancellation as long as $H$ is periodic in time with the same period.

Lemma 4 (Preservation of Spatio-Temporal Periodicity) Let $V_z\left(x,y,z\right)$ be spatially periodic with periods $\left(L_x,L_y,L_z\right)$ , i.e.,

$\begin{array}{l}{V}_z\left(x+L_x,y,z\right)=V_z\left(x,y,z\right),\\ V_z\left(x,y+L_y,z\right)=V_z\left(x,y,z\right),\\ V_z\left(x,y,z+L_z\right)=V_z\left(x,y,z\right).\end{array}$

Let $H\left(x,y,z,t\right)$ be spatially periodic with the same spatial periods and temporally periodic with period $T$ , i.e.,

$\begin{array}{l}H\left(x+L_x,y,z,t\right)=H\left(x,y,z,t\right),\\ H\left(x,y+L_y,z,t\right)=H\left(x,y,z,t\right),\\ H\left(x,y,z+L_z,t\right)=H\left(x,y,z,t\right),\end{array}$

and

$H\left(x,y,z,t+T\right)=H\left(x,y,z,t\right).$

Define the combined velocity,

${\tilde{V}}_z\left(x,y,z,t\right):=V_z\left(x,y,z\right)+H\left(x,y,z,t\right).$

Then, ${\tilde{V}}_z$ is spatially periodic with periods $\left(L_x,L_y,L_z\right)$ and temporally periodic with period $T$ .

Explanation of Oscillatory Behavior and Compensation on the Real Number Line

The compensated Weierstrass zeta velocity is

$V_z\left(x\right)=C_{\text{amp}}\left[\zeta \left(kx\right)-\frac{\eta }{{\Omega }_1}kx\right].$

The raw Weierstrass zeta function $\zeta \left(kx\right)$ is quasi-periodic, meaning that it drifts linearly over the lattice:

$\zeta \left(x+2{\Omega }_1\right)=\zeta \left(x\right)+2\eta .$

The linear term $-\left(\eta /{\Omega }_1\right)kx$ exactly cancels this drift, producing a strictly periodic $V_z\left(x\right)$ . The derivative $\text{d}{V}_z/\text{d}x$ is nonzero only where $\zeta \left(kx\right)$ varies sharply; away from the singularities, the derivative is small, which is why $V_z$ appears nearly constant. When we superimpose a sinusoidal function

$V_z\left(x\right)+A\sin \left(\frac{\pi x}{{\Omega }_1}\right),$

the sinusoid is designed to vanish at the lattice points $x=n\cdot 2{\Omega }_1$ (the poles of $\zeta$ ), because

$\sin \left(\frac{\pi n2{\Omega }_1}{{\Omega }_1}\right)=\sin \left(2n\pi \right)=0.$

Therefore, the superposition does not perturb the velocity at the poles, and the resulting function retains exact periodicity.

See Appendix B for a calculation using Maple 2024 on how a linear term can offset a quasi periodic WeierstrassZeta function to create a sum that is periodic. So finally, we have that: We can always add the smooth parts of the WeierstrassZeta function together with the offset to $H$ which make it periodic and smooth in space except at the poles. We then add this result to a function that is pure periodic added to the negative of $H$ mod poles (with cancellation) which makes it purely periodic together with the poles.

2.7. On Adding a General Spatio-Temporal Function to $V_z\left(x,y,z,t\right)$

We consider a velocity component $V_z\left(x,y,z,t\right)$ constructed so as to cancel a singular temporal evolution through the introduction of a compensating term. In the original construction, one writes

$V_z=I+f_2,$

where $I$ is a temporally singular integral term and $f_2$ is chosen, so that

${\partial }_t{f}_2=0.$

This condition ensures that no new time dependence is introduced during the singular cancellation process.

The question arises when a further perturbation

$H\left(x,y,z,t\right)$

is added to the velocity component, forming

${\tilde{V}}_z=V_z+H.$

If $H$ depends on time, then in general ${\partial }_{t}H\ne 0$ . The purpose of this section is to explain why such a perturbation can nevertheless be absorbed without destroying periodicity or reintroducing secular growth, provided $H$ is periodic in time. The earlier requirement,

${\partial }_t{f}_2=0$

was not a general prohibition against time dependence, but a condition imposed at a specific stage of the construction. At that stage:

A singular temporal integral was being matched to a target potential, Exact annihilation of internal singular evolution was required, no additional time dependence could be introduced by the compensator.

Thus, $f_2=f_2\left(x,y,z\right)$ was required to be time-independent in order to complete the singular cancellation.

The addition of $H$ occurs after the singular cancellation and gauge fixing have already been performed. The new object

${\tilde{V}}_z=V_z+H$

is no longer part of the compensator used to cancel singular terms. Consequently, the requirement is no longer that

${\partial }_{t}H=0.$

Instead, the physically meaningful requirement becomes the absence of secular growth. The Navier-Stokes equations do not forbid time dependence itself; they forbid unbounded growth. The correct condition is therefore

${\displaystyle {\int }_0^T{\partial }_{t}H\left(x,y,z,t\right)\text{d}t}=0,$

where $T=\frac{2\pi }{\omega }$ is the temporal period.

Equivalently,

$H\left(x,y,z,t+T\right)=H\left(x,y,z,t\right).$

This condition replaces the earlier constraint ${\partial }_t{f}_2=0$ in the post-construction setting. Compute the time derivative of the modified velocity:

${\partial }_t{\tilde{V}}_z={\partial }_t{V}_z+{\partial }_{t}H.$

Integrating over one period yields

${\displaystyle {\int }_0^T{\partial }_t{\tilde{V}}_z\text{d}t}={\tilde{V}}_z\left(T\right)-{\tilde{V}}_z\left(0\right).$

Since both $V_z$ and $H$ are periodic with period $T$ ,

$V_z\left(T\right)=V_z\left(0\right),H\left(T\right)=H\left(0\right),$

and hence

${\tilde{V}}_z\left(T\right)-{\tilde{V}}_z\left(0\right)=0.$

Therefore, no secular drift is introduced.

2.7.1. Pressure Gauge Considerations

The pressure-velocity relation may be written schematically as

$\frac{p}{\rho }=-{\partial }_t\Phi -\frac{1}{2}{\left|\nabla \Phi \right|}^2+\cdots$

If $H$ is periodic in time, then ${\partial }_{t}H$ has zero mean over one period and contributes only oscillatory terms to the pressure. No secular pressure gradient arises, and the pressure gauge freedom

$\Phi \mapsto \Phi -C\left(t\right)$

remains intact. If instead

$H\left(x,y,z,t\right)=a\left(x,y,z\right)t+\widehat{H}\left(x,y,z,t\right),$

where $\widehat{H}$ is periodic in time, then

${\displaystyle {\int }_0^T{\partial }_{t}H\text{d}t}=a\left(x,y,z\right)T\ne 0.$

Such a term cannot be absorbed. It would reintroduce secular growth, destroy compactness on ${\mathbb{T}}^3$ , and violate bounded energy conditions.

2.7.2. Periodic and Oscillatory Attractor

It is therefore incorrect to require

${\partial }_{t}H=0.$

The correct condition is

$H\left(x,y,z,t+T\right)=H\left(x,y,z,t\right)$

or equivalently

${\displaystyle {\int }_0^T{\partial }_{t}H\left(x,y,z,t\right)\text{d}t}=0.$

Under the sole requirement of temporal periodicity, spatio-temporal smooth perturbations $H\left(x,y,z,t\right)$ may be absorbed into the vertical velocity component without destroying periodicity, boundedness, or energy control. The original condition ${\partial }_t{f}_2=0$ applies only at the singular cancellation stage and does not restrict later periodic perturbations. Now, from Equations (5) - (7), if the Laplacian of pressure is written as:

$\Delta P=\frac{H\left(x,y,z,t\right)}{G\left(x,y,z,t\right)},$

where G has singular support at z = z_0 for some z_0 < 0 and smooth for z_0 > 0 and the Laplacian of pressure is smooth due to the function G. Here, however, the function H is smooth and periodic in space and time. Then, from Equations (5) - (7), it can be shown that a LambertW solution occurs if $G$ is singular otherwise $\frac{\text{d}{V}_z}{\text{d}t}$ has a mean part $M$ plus an oscillatory part. From $V_z=H-\omega t$ , differentiating and using the fact that $\frac{\text{d}{V}_z}{\text{d}t}$ is purely oscillatory after choosing $\omega$ as $M$ , then

$\frac{\partial H}{\partial t}=\text{oscillatory}\text{term}+\omega$

Integrating in time:

$H\left(t\right)=\omega t+\text{oscillatory}\text{function}$

but recall $V_z=H-\omega t$ , so the linear term in $\omega t$ cancels. Therefore, after choosing $\omega$ to remove the mean, $V_z\left(x,y,z,t\right)$ is purely oscillatory in time, up to an additive constant. This proves that we have an attractor for the 3D Navier-Stokes equations on ${\mathbb{T}}^3$ that is periodic and in particular oscillatory if the pressure Laplacian is such as well.

3. Governing Equations

A unique representation for the PNS system

Consider the function representing the $y_3$ component of the Navier-Stokes equations,

$v_3=v_3\left(y_1,y_2,y_3,s\right)$

and the PDE

$\begin{array}{l}\begin{array}{c}M:=-\frac{\partial v_3}{\partial s}A\mu \left(\delta -1\right)-\frac{2}{3}{\mu }^2\left(\delta -1\right)\frac{\partial v_3}{\partial s}B\\ +\rho v_3\frac{{\partial }^2{v}_3}{\partial y_3^2}+\rho {\left(\frac{\partial v_3}{\partial y_3}\right)}^2+\frac{1}{3}\frac{{\partial }^2{v}_3}{\partial y_1\partial y_3}=0\end{array}\hfill \end{array}$ (65)

where $\rho$ is the density and $\mu$ is the dynamic viscosity of the fluid. There was a relabeling of $V_z$ to $v_3$ , $s=t$ and $y_i=x_i$ . The function $\delta$ is defined to be $1/\sqrt{\kappa 1\left(y_1,y_2,y_3,s\right)}$ and is related to the WeierstrassP function as described previously in this work. The problem becomes more uniform in space variation when the variants $g_2$ and $g_3$ are chosen large enough. So, $u=v/\delta$ where $u$ and $v$ are any two vector solutions to the representing PDE for PNS system in ${\mathbb{T}}^3$ [2]. We perform a substitution to reduce this PDE to an ODE in a single variable $\Xi$ , defined by

$\Xi :=C_1{y}_1+C_2{y}_2+C_3{y}_3+Cs+C_5+C_6,$ (66)

where $C$ absorbs the constant factor in front of $s$ .

Assume

$v_3\left(y_1,y_2,y_3,s\right)=V\left(\Xi \right)$

for some function $V$ of one variable. By the chain rule, the derivatives transform as follows:

$\frac{\partial v_3}{\partial s}=C V^{\prime }\left(\Xi \right),$

$\frac{\partial v_3}{\partial y_3}= C_3 V^{\prime }\left(\Xi \right),$

$\frac{{\partial }^2{v}_3}{\partial y_3^2}=C_3^2 V^{″}\left(\Xi \right),$

$\frac{{\partial }^2{v}_3}{\partial y_1\partial y_3}=C_1{C}_3 V^{″}\left(\Xi \right),$

where primes denote derivatives with respect to $\Xi$ .

Substitution into the PDE

Substituting these into (65), we obtain:

$\begin{array}{c}M=-\left(C{V}^{\prime }\right)A\mu \left(\delta -1\right)-\left(C{V}^{\prime }\right)\frac{2}{3}B{\mu }^2\left(\delta -1\right)\\ +\rho V\cdot \left(C_3^2{V}^{″}\right)+\rho {\left(C_3{V}^{\prime }\right)}^2+\frac{1}{3}\left(C_1{C}_3{V}^{″}\right).\end{array}$

Collecting like terms gives the ODE

$\left(\rho C_3^2 V+\frac{1}{3} C_1{C}_3\right)V^{″}+\rho C_3^2{\left(V^{\prime }\right)}^2-C\left(A\mu \left(\delta -1\right)+\frac{2}{3}B{\mu }^2\left(\delta -1\right)\right)V^{\prime }=0.$ (67)

Define

$\kappa :=\rho C_3^2,\alpha :=\frac{1}{3}{C}_1{C}_3,\beta :=C\left(A\mu \left(\delta -1\right)+\frac{2}{3}B{\mu }^2\left(\delta -1\right)\right).$

Then, (67) becomes

$\left(\kappa V+\alpha \right)V^{″}+\kappa {\left(V^{\prime }\right)}^2-\beta V^{\prime }=0.$ (68)

Equivalently,

$\left(\kappa V+\alpha \right)V^{″}=\beta V^{\prime }-\kappa {\left(V^{\prime }\right)}^2.$

This is the reduced second-order nonlinear ODE for $V\left(\Xi \right)$ . Any solution $V\left(\Xi \right)$ of this ODE yields a solution

$v_3\left(y_1,y_2,y_3,s\right)=V\left(C_1{y}_1+C_2{y}_2+C_3{y}_3+Cs+C_5+C_6\right)$

of the original PDE.

If $C_3=0$ , the ODE degenerates to a linear equation in $V$ :

$\alpha V^{″}-\beta V^{\prime }=0.$

If $\kappa \ne 0$ but $\kappa V+\alpha$ vanishes, the ODE has singular points. One can also define $W\left(\Xi \right):=V^{\prime }\left(\Xi \right)$ and reduce (68) to a first-order ODE in $W$ if desired.

Assume that

$\Xi :=C_1{y}_1+C_2{y}_2+C_3{y}_3+Cs+C_5+C_6$

and enforce the constraints

$C_2=-C_1,C_3=2C_1.$

Assume that $v_3$ depends on $\left(y_1,y_2,y_3,s\right)$ only through $\Xi$ :

$v_3\left(y_1,y_2,y_3,s\right)=f\left(\Xi \right).$

$\frac{\partial v_3}{\partial s}=f^{\prime }\left(\Xi \right)\frac{\partial \Xi }{\partial s}=f^{\prime }\left(\Xi \right)C{C}_5,$

$\frac{\partial v_3}{\partial y_1}=f^{\prime }\left(\Xi \right)C_1,\frac{\partial v_3}{\partial y_2}=f^{\prime }\left(\Xi \right)C_2=-f^{\prime }\left(\Xi \right)C_1,\frac{\partial v_3}{\partial y_3}=f^{\prime }\left(\Xi \right)C_3=2f^{\prime }\left(\Xi \right)C_1.$

Computing second derivatives,

$\frac{{\partial }^2{v}_3}{\partial y_3^2}=f^{″}\left(\Xi \right){\left(C_3\right)}^2=4f^{″}\left(\Xi \right)C_1^2,$

$\frac{{\partial }^2{v}_3}{\partial y_1\partial y_3}=f^{″}\left(\Xi \right)C_1{C}_3=2f^{″}\left(\Xi \right)C_1^2.$

Substituting into $M$ PDE,

Original PDE:

$\begin{array}{c}M=-\frac{\partial v_3}{\partial s}A\mu \left(\delta -1\right)-\frac{2}{3}{\mu }^2\left(\delta -1\right)B\frac{\partial v_3}{\partial s}+\\ \rho v_3\frac{{\partial }^2{v}_3}{\partial y_3^2}+\rho {\left(\frac{\partial v_3}{\partial y_3}\right)}^2+\frac{1}{3}\frac{{\partial }^2{v}_3}{\partial y_1\partial y_3}.\end{array}$

Substitute derivatives in terms of $f\left(\Xi \right)$ and constants:

$-\frac{\partial v_3}{\partial s}A\mu \left(\delta -1\right)-\frac{2}{3}{\mu }^2\left(\delta -1\right)B\frac{\partial v_3}{\partial s}=-f^{\prime }\left(\Xi \right)C{C}_5\left[A\mu \left(\delta -1\right)+\frac{2}{3}{\mu }^2\left(\delta -1\right)B\right],$

$\rho v_3\frac{{\partial }^2{v}_3}{\partial y_3^2}=4\rho f\left(\Xi \right)f^{″}\left(\Xi \right)C_1^2,$

$\rho {\left(\frac{\partial v_3}{\partial y_3}\right)}^2=4\rho f^{\prime }{\left(\Xi \right)}^2{C}_1^2,$

$\frac{1}{3}\frac{{\partial }^2{v}_3}{\partial y_1\partial y_3}=\frac{2}{3}{f}^{″}\left(\Xi \right)C_1^2.$

Combine into ODE in $f\left(\Xi \right)$ .

$\begin{array}{l}-f^{\prime }\left(\Xi \right)C{C}_5\left[A\mu \left(\delta -1\right)+\frac{2}{3}{\mu }^2\left(\delta -1\right)B\right]+4\rho f\left(\Xi \right)f^{″}\left(\Xi \right)C_1^2\\ +4\rho f^{\prime }{\left(\Xi \right)}^2{C}_1^2-\frac{2}{3}{f}^{″}\left(\Xi \right)C_1^2=0.\end{array}$

Factor $C_1^2$ where possible:

$\begin{array}{l}4\rho f\left(\Xi \right)f^{″}\left(\Xi \right)+4\rho f^{\prime }{\left(\Xi \right)}^2-\frac{2}{3}{f}^{″}\left(\Xi \right)\\ -\frac{C{C}_5}{C_1^2}\left[A\mu \left(\delta -1\right)+\frac{2}{3}{\mu }^2\left(\delta -1\right)B\right]f^{\prime }\left(\Xi \right)=0.\end{array}$

Define

$\kappa :=\frac{C{C}_5}{C_1^2}\left[A\mu \left(\delta -1\right)+\frac{2}{3}{\mu }^2\left(\delta -1\right)B\right].$

Then, the PDE reduces to the nonlinear second-order ODE:

$\left(4\rho f\left(\Xi \right)-\frac{2}{3}\right)f^{″}\left(\Xi \right)+4\rho {\left(f^{\prime }\left(\Xi \right)\right)}^2-\kappa f^{\prime }\left(\Xi \right)=0.$

Hence, we have the following representitive system to the 3D Navier-Stokes equations which solves the PNS system in the ball and sphere inscribed in $T={\left[0,1\right]}^3$ (it will be nonsmooth on the plane $y_1-y_2+2\ast y_3=0$ ),

$\left(4\rho \left(f\left(\Xi \right)+C\right)-\frac{2}{3}\right)\frac{{\text{d}}^{2}f}{\text{d}{\Xi }^2}+4\rho {\left(\frac{\text{d}f}{\text{d}\Xi }\right)}^2-\kappa \frac{\text{d}f}{\text{d}\Xi }=0.$ (69)

4. The Geometry of Spheres and Where the $W$ Function Loses Smoothness

For the following $y_2=y_1^2+1/6$ , $y_1=y_3^2+1/6$ and $y_3=y_2^2+1/6$

and,

$\begin{array}{l}{v}_2=-2y_1{v}_1\\ v_1=-2y_3{v}_3\\ v_3=-2y_2{v}_2\end{array}$

$\begin{array}{l}{\mathcal{J}}_3\left(y_3,s\right)=\{s\in {ℝ}^{+},y\in B\left(y_{c_i};R\right):2y_1{v}_1+v_2=0\text{and}A{y}_1+B{y}_2+C=0,\\ y_1,y_2\in I\times I\left(I\subset ℝ\right)\text{and}{y}_2=y_1^2+1/6\\ \text{and}{v}_3\left(y_1,y_2,y_3,s\right)\in C^0\left({\mathbb{T}}^3\right)\}\end{array}$ (70)

where $R$ is a sufficiently large positive number.

It can be calculated that on this subspace of solutions to PNS equations, the following is identically zero,

$\begin{array}{c}\mathcal{X}=\left(\left(\delta -1\right)v_1\frac{\partial v_3}{\partial s}+2\rho v_3\frac{\partial v_1}{\partial s}\right)\frac{{\partial }^2{v}_3}{\partial y_3\partial y_1}\\ +\left(\left(\delta -1\right)v_2\frac{\partial v_3}{\partial s}+2\rho v_3\frac{\partial v_2}{\partial s}\right)\frac{{\partial }^2{v}_3}{\partial y_3\partial y_2}\\ -\frac{\partial v_3}{\partial s}[v_3\frac{\partial v_3}{\partial y_1}\frac{{\partial }^2{v}_1}{\partial y_3\partial s}+v_3\frac{\partial v_3}{\partial y_2}\frac{{\partial }^2{v}_2}{\partial y_3\partial s}-\frac{\partial v_3}{\partial y_1}\frac{\partial v_3}{\partial y_3}\frac{\partial v_1}{\partial s}\\ -\frac{\partial v_3}{\partial y_2}\frac{\partial v_3}{\partial y_3}\frac{\partial v_2}{\partial s}+v_3\frac{\partial v_3}{\partial y_1}\frac{\partial v_1}{\partial s}\frac{{\partial }^2{v}_3}{\partial s\partial y_3}+v_3\frac{\partial v_3}{\partial y_2}\frac{\partial v_2}{\partial s}\frac{{\partial }^2{v}_3}{\partial s\partial y_3}]\\ =0\end{array}$ (71)

The expression $\mathcal{X}$ has been obtained in [2] (Ch. 8) as part of an integral form of the periodic Navier-Stokes equations. Continuing we have the other two spaces for the 3D Navier-Stokes equations,

$\begin{array}{l}{\mathcal{J}}_2\left(y_2,s\right)=\{s\in {ℝ}^{+},y\in B\left(y_{c_i};R\right):2y_3{v}_3+v_1=0\text{and}A{y}_3+B{y}_1+C=0,\\ y_1,y_3\in I\times I\left(I\subset ℝ\right)\text{and}{y}_1=y_3^2+1/6\text{and}\\ v_2\left(y_1,y_2,y_3,s\right)\in C^0\left({\mathbb{T}}^3\right)\}\end{array}$ (72)

$\begin{array}{l}{\mathcal{J}}_1\left(y_1,s\right)=\{s\in {ℝ}^{+},y\in B\left(y_{c_i};R\right):2y_2{v}_2+v_3=0\text{and}A{y}_2+B{y}_3+C=0,\\ y_2,y_3\in I\times I\left(I\subset ℝ\right)\text{and}{y}_3=y_2^2+1/6\\ \text{and}{v}_1\left(y_1,y_2,y_3,s\right)\in C^0\left({\mathbb{T}}^3\right)\}\end{array}$ (73)

Consider the quadratic equation defining a sphere:

$y_1^2+y_2^2+y_3^2-y_1-y_2-y_3+\frac{1}{2}=0.$ (74)

Complete the square for each coordinate:

$y_i^2-y_i={\left(y_i-\frac{1}{2}\right)}^2-\frac{1}{4},i=1,2,3.$

Substituting into (74):

${\displaystyle \sum }_{i=1}^3{\left(y_i-\frac{1}{2}\right)}^2-\frac{3}{4}+\frac{1}{2}=0\Rightarrow {\displaystyle \sum }_{i=1}^3{\left(y_i-\frac{1}{2}\right)}^2=\frac{1}{4}.$

Hence, the sphere has:

$\text{center}:C=\left(\frac{1}{2},\frac{1}{2},\frac{1}{2}\right),\text{radius}:r=\frac{1}{2}.$

It can be verified that,

$\mathcal{X}\equiv 0.$

Let $y_1,y_2,y_3$ satisfy

$y_2=y_1^2+\frac{1}{6},y_1=y_3^2+\frac{1}{6},y_3=y_2^2+\frac{1}{6},$

and lie on the sphere

$y_1^2+y_2^2+y_3^2-y_1-y_2-y_3+\frac{1}{2}=0.$

Because the three relations are cyclic and identical in form, a natural family of solutions is the diagonal one $y_1=y_2=y_3=y$ . Substituting into the iteration gives the scalar equation

$y=y^2+\frac{1}{6},$

or

$y^2-y+\frac{1}{6}=0.$

The quadratic has discriminant $\text{Δ}=1-4\cdot \frac{1}{6}=1-\frac{2}{3}=\frac{1}{3}$ , so

$y=\frac{1\pm \sqrt{1/3}}{2}=\frac{1\pm \frac{1}{\sqrt{3}}}{2}.$

Write the two roots as

$y^{(\pm )}=\frac{1\pm 1/\sqrt{3}}{2}.$

For these values,

$3y^2-3y+\frac{1}{2}=3\left(y^2-y+\frac{1}{6}\right)=0,$

so the sphere equation is satisfied. Thus, the points on the sphere solving the system are the two diagonal points

$\left(y_1,y_2,y_3\right)=\left(y^{(+)},y^{(+)},y^{(+)}\right)\text{and}\left(y_1,y_2,y_3\right)=\left(y^{(-)},y^{(-)},y^{(-)}\right),$

that is,

$\left(\frac{1+1/\sqrt{3}}{2},\frac{1+1/\sqrt{3}}{2},\frac{1+1/\sqrt{3}}{2}\right)\text{and}\left(\frac{1-1/\sqrt{3}}{2},\frac{1-1/\sqrt{3}}{2},\frac{1-1/\sqrt{3}}{2}\right).$

Numeric approximations:

$y^{(+)}\approx 0.788675,y^{(-)}\approx \mathrm{0.211325.}$

(These diagonal solutions satisfy both the cyclic relations and the sphere. See Figure 3. Non-diagonal solutions could in principle exist for the cyclic system, but the two symmetric diagonal points are the solutions that lie on the given sphere.)

As a side note we are required to check the continuity equation using the cyclic relations for the $v_i$ which becomes,

$\frac{\partial v_1}{\partial y_1}-2y_1\frac{\partial v_1}{\partial y_2}-\frac{1}{2y_3}\frac{\partial v_1}{\partial y_3}+\frac{v_1}{2y_3^2}=0$

and also to see if a similar LambertW solution for $v_1$ exists from this PDE, which can be proven to be the case. The following proves useful,

Let $W\left(z\right)$ denote the Lambert $W$ function, defined by

$W\left(z\right)e^{W\left(z\right)}=z.$

Then, for $x>0$ and $y>0$ , the following identity holds exactly:

$W\left(x\right)+W\left(y\right)=W\left(xy\left(\frac{1}{W\left(x\right)}+\frac{1}{W\left(y\right)}\right)\right)$ (I)

Figure 3. Diagonal solutions (points in black and red on the sphere) satisfy both the cyclic relations and the sphere.

Proof. Set $X=W\left(x\right)$ and $Y=W\left(y\right)$ . Then

$x=X{e}^X,y=Y{e}^Y,$

and

$xy\left(\frac{1}{X}+\frac{1}{Y}\right)=XY{e}^{X+Y}\frac{X+Y}{XY}=\left(X+Y\right)e^{X+Y}.$

By definition of $W$ , we have

$W\left(xy\left(\frac{1}{W\left(x\right)}+\frac{1}{W\left(y\right)}\right)\right)=W\left(\left(X+Y\right)e^{X+Y}\right)=X+Y=W\left(x\right)+W\left(y\right).$

As a consequence, this identity also holds asymptotically for $x\gg 1$ and $y\gg 1$ .

The solution of the continuity equation for $v_1$ is for general $F_1$ ,:

$v_1=F_1\left(y_1^2+y_2,y_3^2+y_1,s\right)F\left(y_3\right)y_3=\frac{e^{W\left(U\right)}}{e^{-W\left(y_3\right)}}$

so $v_1=e^{W\left(U\right)+W\left(y_3\right)}$ and $F\left(y_3\right)=W\left(y_3\right)/y_3$ , for $y_3\ne 0$ . We use the identity Equation (I) above to determine $v_1$ . Keeping in mind that the cyclic conditions are used, $y_2=y_1^2+1/6$ and $y_1=y_3^2+1/6$ which together on the sphere with a shift in order to be on ${\left[0,1\right]}^3$ gives $2y_1^2+2y_3^2=-2y_2^2+1/2$ . Substituting into the form for $v_1$ gives,

$v_1=e^{W\left(e^{-2y_2^2-s-1}\left(y_2{}^2+1/6\right)\left({\left(W\left(-e^{-2y_2^2-s-1}\right)\right)}^{-1}+{\left(W\left(y_2^2+1/6\right)\right)}^{-1}\right)\right)}$

and noting also that the constant 1/6 was cancelled by adding a constant $C_6$ to the time term $s$ in the general expression for $F_1$ . Finally, it is seen that $x=U=e^{-2y_2^2-s-1}>0$ and since $y=y_3=y_2^2+1/6>0$ the implicit formula for the sum $W\left(U\right)+W\left(y_2^2+1/6\right)$ can be expressed as a LambertW function as in Equation (I).

Changing back to $x_i$ variables we consider the intersection of the sphere

$x^2+y^2+z^2=\frac{1}{4}$ (75)

and the cubic surface

$xyz=\frac{1}{8d},$ (76)

for some $d>0$ . The cubic surface is obtained by solving the previous cyclic relations and transferring to $x$ variables through $\delta <0$ . To determine when real intersections exist, we apply the arithmetic-geometric mean inequality

$\frac{x^2+y^2+z^2}{3}\ge {\left(\left|xyz\right|\right)}^{2/3}.$ (77)

Substituting the defining equations of the surfaces gives

$\frac{1/4}{3}\ge {\left(\frac{1}{8d}\right)}^{2/3},$ (78)

Figure 4. Singularities of $v_3$ in red ring intersection regions of $xyz=1/\left(8d\right)$ and $x^2+y^2+z^2=1/4$ , $d\ge 3\sqrt{3}$ .

Figure 5. Singularities of $v_3$ in ring intersection regions of $\left(x-1/2\right)\left(y-1/2\right)\left(z-1/2\right)=1/\left(8d\right)$ , ${\left(x-1/2\right)}^2+{\left(y-1/2\right)}^2+{\left(z-1/2\right)}^2=1/4$ , and $x-y+2z=-1$ , $d\ge 3\sqrt{3}$ .

or equivalently

$\frac{1}{12}\ge {\left(8d\right)}^{-2/3}.$ (79)

Rearranging yields

${\left(8d\right)}^{2/3}\ge 12,$ (80)

and hence

$8d\ge {12}^{3/2}=12\sqrt{12}=24\sqrt{3}.$ (81)

Therefore, real solutions exist if and only if

$d\ge 3\sqrt{3}\approx \mathrm{5.196.}$ (82)

For $d\ge 3\sqrt{3}$ , the cubic surface $xyz=\frac{1}{8d}$ intersects the sphere $x^2+y^2+z^2=\frac{1}{4}$ along a one-dimensional curve in ${ℝ}^3$ . For $d<3\sqrt{3}$ , the two surfaces do not intersect. See Figure 4 where the solution to governing PDE $M=0$ earlier has singularities on the red rings in the plot for the sphere centered at 0. Next look at Figure 5 for the plot for the sphere centered at $x_c=\left(1/2,1/2,1/2\right)$ .

4.1. Geometric Interpretation of the Centered Cubic, Plane, and Sphere on [0, 1]³

We consider three geometric objects embedded in the unit cube

${\left[0,1\right]}^3\subset {ℝ}^3,$

all expressed relative to the distinguished center point

$c:=\left(\frac{1}{2},\frac{1}{2},\frac{1}{2}\right).$

4.1.1. Centered Coordinate System

Introduce shifted coordinates

$X=x-\frac{1}{2},Y=y-\frac{1}{2},Z=z-\frac{1}{2},$

so that the cube ${\left[0,1\right]}^3$ is recentered symmetrically about the origin

$\left(X,Y,Z\right)\in {\left[-\frac{1}{2},\frac{1}{2}\right]}^3.$

all three surfaces in the Maple program used are naturally expressed in these centered variables.

4.1.2. The Cubic Surface

The cubic surface is defined by

$\left(x-\frac{1}{2}\right)\left(y-\frac{1}{2}\right)\left(z-\frac{1}{2}\right)=\frac{1}{8d},$

or equivalently

$XYZ=\frac{1}{8d}.$

This is a translated and rescaled version of the standard cubic surface

$xyz=\text{constant},$

with the following properties:

• The coordinate planes $X=0$ , $Y=0$ , $Z=0$ act as asymptotic planes.

• The surface has eight connected components in ${ℝ}^3$ , though only portions lying inside ${\left[0,1\right]}^3$ are rendered.

• The parameter $d>0$ controls the distance of the surface from the coordinate planes; larger $d$ pushes the surface closer to the center.

This surface models a multiplicative constraint symmetric about the cube center.

4.1.3. The Plane

The plane plotted in Maple is

$\left(x-\frac{1}{2}\right)-\left(y-\frac{1}{2}\right)+2\left(z-\frac{1}{2}\right)=C_{\text{plane}},$

which simplifies to

$X-Y+2Z=C_{\text{plane}}.$

This is an affine plane with normal vector

$\nabla \Phi =\left(1,-1,2\right),$

and represents a linear constraint passing through or near the centered cube depending on the value of $C_{\text{plane}}$ .

Geometrically:

• The plane slices the cube obliquely.

• Its orientation is fixed; only its offset changes.

• It is positioned to examine tangency or near-tangency with the cubic surface.

4.1.4. The Sphere

The sphere is defined by

${\left(x-\frac{1}{2}\right)}^2+{\left(y-\frac{1}{2}\right)}^2+{\left(z-\frac{1}{2}\right)}^2=\frac{1}{4},$

or

$X^2+Y^2+Z^2={\left(\frac{1}{2}\right)}^2.$

Thus, the sphere:

• is centered exactly at $c=\left(\frac{1}{2},\frac{1}{2},\frac{1}{2}\right)$ ,

• has radius $r=\frac{1}{2}$ ,

• is the largest sphere fully contained inside the unit cube.

This sphere provides a natural geometric reference for symmetry and compactness.

4.1.5. Tangency Structure

At special parameter values, the plane may be tangent to:

• the cubic surface,

• the sphere,

• or both simultaneously.

Tangency occurs when the surface normal of the cubic,

$\nabla \left(XYZ\right)=\left(YZ,XZ,XY\right),$

aligns with the plane normal $\left(1,-1,2\right)$ at a point satisfying

$XYZ=\frac{1}{8d}.$

Because the surfaces are centered, this tangency structure is symmetric with respect to sign changes in $\left(X,Y,Z\right)$ , though only the portion inside ${\left[0,1\right]}^3$ is visible. The value of $C_{\text{plane}}=-1$ in the plot.

4.2. Cyclic Relations (or Inequalities) for the Exact Solution of PNS System

Consider a triple of real variables $y=\left(y_1,y_2,y_3\right)\in {ℝ}^3$ satisfying the following cyclic relations (or inequalities):

$y_2\ge y_1^2+\frac{1}{6},y_3\ge y_2^2+\frac{1}{6},y_1\ge y_3^2+\frac{1}{6}.$

These inequalities capture a natural cyclic dependency among the coordinates, reminiscent of iterative constraints that arise in discrete dynamical systems and geometric embeddings.

Sum over $i$

Consider the three cyclic relations (or inequalities)

$y_2=y_1^2+\frac{1}{6},y_3=y_2^2+\frac{1}{6},y_1=y_3^2+\frac{1}{6}.$

Summing the three identities gives

$y_1+y_2+y_3=y_1^2+y_2^2+y_3^2+\frac{1}{2},$

which rearranges to

$y_1^2+y_2^2+y_3^2-\left(y_1+y_2+y_3\right)+\frac{1}{2}=0.$

Completing the square yields the sphere equation

${\displaystyle \sum }_{i=1}^3{\left(y_i-\frac{1}{2}\right)}^2=\frac{1}{4}.$

If instead the three relations are replaced by the coordinate wise inequalities

$y_2\le y_1^2+\frac{1}{6},y_3\le y_2^2+\frac{1}{6},y_1\le y_3^2+\frac{1}{6},$

then summing yields

$y_1+y_2+y_3\le y_1^2+y_2^2+y_3^2+\frac{1}{2},$

so

${\displaystyle \sum }_{i=1}^3{\left(y_i-\frac{1}{2}\right)}^2\ge \frac{1}{4},$

i.e., every point satisfying the three “≤” inequalities lies on or outside the sphere (outside the open ball of radius 1/2 centered at $\left(\frac{1}{2},\frac{1}{2},\frac{1}{2}\right)$ ).

Conversely, if

$y_2\ge y_1^2+\frac{1}{6},y_3\ge y_2^2+\frac{1}{6},y_1\ge y_3^2+\frac{1}{6},$

then summing gives

${\displaystyle \sum }_{i=1}^3{\left(y_i-\frac{1}{2}\right)}^2\le \frac{1}{4},$

i.e., such points lie inside or on the closed ball of radius 1/2.

Remark

The implications above are one-way: the collection of three coordinatewise inequalities (all “≥” or all “≤”) implies the corresponding inclusion relative to the ball, but being inside (or outside) the ball does not force each coordinate inequality individually.

Let

$A:=-2y_1^2+\frac{1}{6},B:=-2y_3^2+\frac{1}{6},C:=-2y_2^2+\frac{1}{6},$

so the system of inequalities is

$v_3\le A{v}_2,v_2\le B{v}_1,v_1\le C{v}_3.$

Multiply the three inequalities (all real numbers), and set $P:=v_1{v}_2{v}_3$ . This gives

$P\le \left(ABC\right) P.$

Rearranging yields

$0\le \left(ABC-1\right)P.$

Thus, the product $\left(ABC-1\right)P$ must be nonnegative. Equivalently:

$\{\begin{array}{l}ABC>1\Rightarrow P\ge 0,\hfill \\ ABC<1\Rightarrow P\le 0,\hfill \\ ABC=1\Rightarrow \text{no information on}P\text{from the product inequality}.\hfill \end{array}$

Using the specific form of the coefficients, observe that

$A,B,C\le \frac{1}{6},$

because $-2y_i^2\le 0$ . Hence

$ABC\le {\left(\frac{1}{6}\right)}^3=\frac{1}{216}<1.$

Therefore $ABC<1$ always, and the second case above applies. Thus, we obtain the global implication

$v_1{v}_2{v}_3\le 0.$

If none of $v_1,v_2,v_3$ is zero, then the product is strictly negative, so an odd number of the $v_i$ must be negative (one or three of them).

5. Lambert-Velocity, the Level $E$ and Interior Smoothness

Consider the scalar velocity-type function

$v\left(y_1,y_2,y_3,s\right)=W\text{}\left(-e^{E\left(y,s\right)}\right)+c,$

where $W$ denotes the principal real branch of the Lambert $W$ -function (so $W:\left[-1/e,0\right]\to \left[-1,0\right]$ real-valued), $c\in ℝ$ is a constant, and the level function is now defined by the sign-changed formula

$E\left(y,s\right):=-y_1-y_2-y_3-s-1.$

The Lambert argument is $z\left(y,s\right)=-e^{E\left(y,s\right)}$ . Now:

$W$ is real-analytic on $\left(-1/e,0\right)$ and strictly increasing there; its only real singular point is $z=-1/e$ (corresponding to $E=-1$ ), where $W^{\prime }$ blows up with square-root type behaviour. Consequently $v$ is real-analytic at $\left(y,s\right)$ if and only if $E\left(y,s\right)<-1$ (equivalently $z\left(y,s\right)\in \left(-1/e,0\right)$ ). The locus $E=-1$ is precisely the potential non-smoothness locus of the Lambert piece.

We work under the cyclic inequalities on the $y$ -variables:

$y_2\ge y_1^2+\frac{1}{6},y_1\ge y_3^2+\frac{1}{6},y_3\ge y_2^2+\frac{1}{6}.$

These define a closed domain (the “ball”)

$\Omega :=\left\{\left(y_1,y_2,y_3\right):y_2\ge y_1^2+\frac{1}{6},y_1\ge y_3^2+\frac{1}{6},y_3\ge y_2^2+\frac{1}{6}\right\},$

and its boundary (the “sphere”) $\partial \Omega$ is the set where all three equalities hold.

It is useful to note the two real fixed points of $f\left(t\right)=t^2+\frac{1}{6}$ :

$\begin{array}{l}{t}^2-t+\frac{1}{6}=0\Rightarrow \alpha =\frac{1-\sqrt{1/3}}{2}\approx 0.211324865405,\\ \beta =\frac{1+\sqrt{1/3}}{2}\approx \mathrm{0.788675134485.}\end{array}$

Proposition 1 (Interior strict inequality for E) Assume the cyclic inequalities above hold on $\Omega$ . Suppose further that the boundary $\partial \Omega$ corresponds to the lower fixed point in the sense that the equalities on $\partial \Omega$ are realized by

$y_1=y_2=y_3=\alpha$ on $\partial \Omega$ .

Choose $s$ so that $E=-1$ on that boundary, i.e.,

$-3\alpha -s-1=-1,$

which is equivalent to $s=-3\alpha$ . Then, every strict interior point of $\Omega$ (a point where at least one of the three defining inequalities is strict) satisfies

$E\left(y,s\right)<-1.$

Consequently, the Lambert argument satisfies $z\left(y,s\right)\in \left(-1/e,0\right)$ strictly in the interior, and the composition $v\left(y,s\right)=W\left(-e^{E\left(y,s\right)}\right)+c$ is real-analytic (smooth) at every interior point.

Proof. First observe that any coordinate of a point in $\Omega$ must satisfy $t\ge \alpha$ . Indeed define $g\left(t\right)=t-\left(t^2+\frac{1}{6}\right)=-\left(t-\alpha \right)\left(t-\beta \right)$ . If $t<\alpha$ then $g\left(t\right)<0$ , i.e., $t<t^2+\frac{1}{6}$ , which contradicts the defining inequality $t\ge t^2+\frac{1}{6}$ . Thus, for every $\left(y_1,y_2,y_3\right)\in \Omega$ , we have

$y_1\ge \alpha ,y_2\ge \alpha ,y_3\ge \alpha ,$

and hence the sum $S\left(y\right):=y_1+y_2+y_3$ satisfies

$S\left(y\right)\ge 3\alpha \text{on}\Omega ,$

with equality only when $y_1=y_2=y_3=\alpha$ .

Now fix $s=-3\alpha$ so that on the boundary triple $\left(\alpha ,\alpha ,\alpha \right)$ we have

$E\left(\alpha ,\alpha ,\alpha ,s\right)=-3\alpha -s-1=-1.$

Take any strict interior point $y\in {\text{Ω}}^{\circ }$ . By strictness at least one coordinate is $>\alpha$ , so $S\left(y\right)>3\alpha$ . Evaluate $E$ at $y$ :

$E\left(y,s\right)=-S\left(y\right)-s-1=-S\left(y\right)+3\alpha -1=-1-\left(S\left(y\right)-3\alpha \right).$

Since $S\left(y\right)-3\alpha >0$ we conclude $E\left(y,s\right)<-1$ . This proves the proposition.

6. Alignment, Elliptic Scaling, and the Nonlinear Inertial Term

6.1. Aligned Decompositions and Elliptic Representations

Let $u\left(x,t\right)$ and $v\left(x,t\right)$ be velocity fields defined on ${ℝ}^3\times \left[0,T\right)$ . We say that $u$ and $v$ are aligned if there exists a scalar function $\text{Λ}\left(x,t\right)$ such that

$u\left(x,t\right)=\text{Λ}\left(x,t\right)v\left(x,t\right).$ (83)

We are particularly interested in the case where $\text{Λ}$ admits an elliptic representation of the form

$\Lambda \left(x,t\right)={\wp }^{-1}\left(G\left(x,t\right)\right)\text{or}\Lambda \left(x,t\right)=\zeta \left(\Theta \right)-\frac{\eta }{\Omega }\Theta ,$ (84)

where $\wp$ and $\zeta$ denote the Weierstrass elliptic and zeta functions, respectively, and $\Theta =k\cdot x-\omega t$ is an affine space-time phase. Such representations encode isolated, lattice-structured singularities with explicit algebraic growth rates.

A fundamental observation is that all singular behavior of $u$ is captured entirely by $\text{Λ}$ when $v$ remains smooth.

6.2. Effect of Alignment on the Nonlinear Inertial Term

The Navier-Stokes inertial term is given by

$\left(u\cdot \nabla \right)u.$

Substituting the aligned form (83), one obtains

$\left(u\cdot \nabla \right)u={\text{Λ}}^2\left(v\cdot \nabla \right)v+\text{Λ}\left(v\cdot \nabla \text{Λ}\right) v.$ (85)

This decomposition shows that:

The nonlinear term does not generate new directions, all singular amplification arises solely through derivatives of $\text{Λ}$ , and the geometry of the nonlinearity is degenerate under alignment. In particular, if $\text{Λ}$ is elliptic, then singularities of the inertial term are inherited from elliptic poles and are not dynamically generated.

6.3. Elliptic Control of Singularities

Let $\text{Λ}={\wp }^{-1}\left(G\right)$ . Then

$\nabla \text{Λ}=\frac{\nabla G}{{\wp }^{\prime }\left({\wp }^{-1}\left(G\right)\right)},$

and singularities occur only when ${\wp }^{\prime }\left({\wp }^{-1}\left(G\right)\right)=0$ , corresponding to the half-period values of the elliptic lattice. Thus, singularities are isolated in space-time, their locations are fixed geometrically, their growth rates are explicitly computable. The inertial term therefore transports regions where elliptic singularities occur rather than creating them.

6.4. Decomposition into Aligned and Non-Aligned Components

To analyze departures from alignment, we introduce the canonical decomposition

$u=\text{Λ}v+w,w\cdot v=0\text{pointwise}.$ (86)

Here,

$\text{Λ}=\frac{u\cdot v}{{\left|v\right|}^2},$

and $w$ measures the degree of non-alignment. Perfect alignment corresponds to $w\equiv 0$ .

6.5. Nonlinear Inertial Term under Non-Alignment

Substituting (86) into the inertial term yields

$\left(u\cdot \nabla \right)u=\left(\text{Λ}v\cdot \nabla \right)\left(\text{Λ}v\right)+\left(\text{Λ}v\cdot \nabla \right)w+\left(w\cdot \nabla \right)\left(\text{Λ}v\right)+\left(w\cdot \nabla \right)w.$ (87)

The first term coincides with the elliptically controlled aligned contribution (85). Every remaining term involves the non-aligned component $w$ . Consequently, any destabilizing nonlinear interaction necessarily arises through non-alignment.

6.6. Elliptic Thresholds for Non-Alignment

Elliptic representations of $\text{Λ}$ imply the pointwise bound

where denotes the elliptic pole lattice.

The most singular mixed term satisfies,

Thus, finite-time blowup can occur only if the non-aligned component satisfies

(88)

Elliptic alignment therefore imposes a quantitative threshold that non-alignment must exceed in order to trigger singular behavior.

6.7. Vorticity Stretching and Suppression of Blowup

The vorticity equation

${\partial }_t\omega +\left(u\cdot \nabla \right)\omega =\left(\omega \cdot \nabla \right)u+\nu \text{Δ}\omega$

shows that blowup is driven by vortex stretching.

Under the decomposition (86),

$\omega =\nabla \times \left(\text{Λ}v\right)+\nabla \times w.$

Stretching terms involving $\nabla \times \left(\text{Λ}v\right)$ inherit elliptic control, while genuinely destabilizing stretching requires coupling with $\nabla \times w$ . Hence, elliptic alignment suppresses vortex stretching unless non-alignment grows faster than elliptic pole scaling.

6.8. Invariant Elliptic Manifolds

Define the elliptic alignment manifold

${ℳ}_{\wp }=\left\{u=\text{Λ}v|\text{Λ}\text{elliptic},v\in C^{\infty }\right\}.$

This manifold is invariant under the Navier-Stokes nonlinearity in the sense that the inertial term is tangent to ${ℳ}_{\wp }$ . Non-alignment measures the distance of a solution from this manifold. Finite-time blowup corresponds to escape from ${ℳ}_{\wp }$ at a rate that exceeds elliptic control. Elliptic representations of the alignment scalar $\text{Λ}$ provide a geometric framework in which:

1) singularities are fixed geometric objects,

2) the nonlinear inertial term is kinematic rather than amplifying,

3) blowup requires rapid growth of non-alignment.

In this sense, finite-time singularity formation is reframed as a problem of geometric non-alignment rather than purely analytic instability.

Adding a perpendicular component breaks alignment

We have a vector field $u\in {ℝ}^3$ defined as

$u=\lambda \left(x,y,z,t\right)v,$

where $v\in {ℝ}^3$ is another vector field and $\lambda$ is a scalar function.

By definition, two vectors are aligned if one is a scalar multiple of the other. Here, $u$ and $v$ are aligned because $u=\lambda v$ . Now, define

$u_1=\lambda v+w$

where $w\ne 0$ is some vector field orthogonal to $u$ :

$u\cdot w=0.$

Two vectors $a$ and $b$ are aligned if

$b=\mu a$

for some scalar $\mu \in ℝ$ . If $u_1$ were aligned with $u$ , there would exist some $\mu$ such that

$u_1=\mu u.$

Plug in $u_1=u+w$ :

$u+w=\mu u\Rightarrow w=\left(\mu -1\right)u.$

To obtain a contradiction,

Take the dot product with $u$ on both sides:

$u\cdot w=\left(\mu -1\right)u\cdot u=\left(\mu -1\right){\Vert u\Vert }^2.$

By assumption, $u\cdot w=0$ , hence

$\left(\mu -1\right){\Vert u\Vert }^2=0\mu =1\left(\text{assuming}\text{u}\ne 0\right).$

Plug $\mu =1$ back:

$w=\left(\mu -1\right)u=0.$

Since $w\ne 0$ by construction, there is no scalar $\mu$ such that $u_1=\mu u$ . Therefore,

$\boxed{u_1\text{is}\text{not}\text{aligned}\text{with}u.}$

7. Perturbative Analysis with Non-Alignment Component $w$

We consider the scaled velocity decomposition with a non-aligned perturbation $w$ :

$u^{\text{*}}=\frac{1}{\delta }u+w,w\cdot u=0,$ (89)

where $u=\left(u_x,u_y,u_z\right)$ is the aligned component and $w=\left(w_x,w_y,w_z\right)$ is the non-alignment vector. We define,

$b=\frac{1}{\delta }\left(u_{x}i+u_{y}j+u_{z}k\right)+w,a=u_{z}b.$ (90)

7.1. Modified 2D-1D Composite Equations

Multiplying the first two Cartesian components of the scaled NSE by $i,j$ respectively, summing, and substituting Equation (89), we obtain the modified composite equation for $b$ including $w$ :

$\frac{\partial b}{\partial t}+\left(b\cdot \nabla \right)b-\frac{\mu }{\rho }{\nabla }^{2}b+\frac{1}{\rho }{\nabla }_{xy}P={\delta }^2{F}_T+\text{terms from }w.$ (91)

For the $z$ -component, multiply the scaled $z$ -NSE by $b$ and add to $b$ equation multiplied by $u_z$ :

$u_z\frac{\partial a}{\partial t}+a\cdot \nabla a-\frac{\mu }{\rho }{u}_z^2{\nabla }^{2}b+\frac{1}{\rho }{u}_z^2{\nabla }_{xy}P+u_{z}b\nabla \cdot a+{\delta }^2{u}_z^2{F}_T+\cdots =0,$ (92)

where all contributions from $w$ are explicitly included via $b$ (see Equation (90)).

7.2. Geometric Algebra Decomposition with $w$

Define the inertial vector:

$f=a\cdot \nabla a.$ (93)

Scalar part:

$\begin{array}{l}f\cdot \left(u_z^2\frac{\partial b}{\partial t}+u_{z}b\frac{\partial u_z}{\partial t}\right)+{\Vert f\Vert }^2+u_z\left(f\cdot b\right)\nabla \cdot a\\ =\frac{\mu }{\rho }{u}_z^{2}f\cdot {\nabla }^{2}b-\frac{1}{\rho }{u}_z^{2}f\cdot {\nabla }_{xy}P-\frac{1}{\rho }\left(f\cdot b\right)u_z\frac{\partial P}{\partial z}\\ +\frac{\mu }{\rho }{u}_z\left(f\cdot b\right){\nabla }^2{u}_z+f\cdot \text{Force terms}.\end{array}$ (94)

Vector part:

$\begin{array}{l}{u}_z^2\frac{\partial b}{\partial t}+a\frac{\partial u_z}{\partial t}+u_{z}b\nabla \cdot a\\ =\frac{\mu }{\rho }{u}_z^2{\nabla }^{2}b-\frac{1}{\rho }{u}_z^2{\nabla }_{xy}P-\frac{1}{\rho }{u}_{z}b\frac{\partial P}{\partial z}+\frac{\mu }{\rho }{u}_{z}b{\nabla }^2{u}_z+\text{Force terms}.\end{array}$ (95)

7.3. Taking the Divergence with $w$ Contributions

Let

$H=\frac{u_{z}b\cdot f}{{\partial }_t{u}_z}.$

Multiply divergence of Equation (95) by $H$ :

$\begin{array}{l}H\nabla u_z^2\cdot \frac{\partial b}{\partial t}+u_z^{2}H\frac{\partial }{\partial t}\left(\nabla \cdot b\right)+Ha\cdot \nabla \left({\partial }_t{u}_z\right)+H\nabla \cdot \left(b{u}_z\nabla \cdot a\right)\\ -\frac{\mu }{\rho }H\nabla \cdot \left(u_z^2{\nabla }^{2}b\right)+H\nabla \cdot \left(\frac{1}{\rho }{u}_z^2{\nabla }_{xy}P+\frac{1}{\rho }{u}_{z}b{\partial }_{z}P\right)\\ -\frac{\mu }{\rho }H\nabla \left(u_z{\nabla }^2{u}_z\right)\cdot b+H\nabla \cdot \left(\text{Forces}\right)\\ =\text{RHS from scalar eq}.\left(94\right).\end{array}$ (96)

7.4. Division by $u_z^{2}H$ and Omega Definitions with $w$

Define generalized operators [1] including $w$ :

$\begin{array}{l}{\Omega }_1=u_z^{-2}a\cdot \nabla \left({\partial }_t{u}_z\right),\\ {\Omega }_3=\frac{\mu }{\rho }{u}_z^{-2}\nabla \cdot \left(u_z^2{\nabla }^{2}b\right),\\ {\Omega }_4=u_z^{-2}\nabla \cdot \left(b{u}_z\nabla \cdot a\right),\\ {\Omega }_5=u_z^{-2}\nabla u_z^2\cdot {\partial }_{t}b,\\ {\Omega }_6=\frac{\mu }{\rho }{\partial }_t{u}_z{u}_z^{-2}{\nabla }^2{u}_z,\\ {\Omega }_7=u_z^{-2}\nabla \cdot \left(b{\rho }^{-1}{u}_z{\partial }_{z}P\right),\\ {\Omega }_8={\rho }^{-1}{u}_z^{-2}\left({\partial }_t{u}_z\right)\left({\partial }_{z}P\right),\\ Y=\nabla \cdot b,F=u_z^{-3}{\partial }_t{u}_z,G=u_z^{-2}{\left({\partial }_t{u}_z\right)}^2.\end{array}$

All terms now include contributions of $w$ via $b=u/\delta +w$ .

7.5. Final PDE in $u_z$ and $w$

The PDE corresponding to Equation (19) in [1] becomes:

$\begin{array}{l}{\partial }_t\left({\partial }_z{u}_z\right)+{\Omega }_5+\left(\frac{1}{\delta }-1\right)\frac{{\left({\partial }_t{u}_z\right)}^2}{u_z^2}-\frac{1}{u_z^2}\left(u_z{\partial }_z{u}_z{\partial }_t{u}_z\right)+\frac{\mu }{\rho }\frac{1-{\delta }^{-1}}{u_z^2}{\partial }_t{u}_z{\nabla }^2{u}_z\\ -\frac{1}{\rho u_z^2}{\partial }_t{u}_z{\partial }_{z}P+\frac{1}{u_z^2}{\displaystyle {\iint }_S\left(\frac{1}{\rho }{u}_z^2{\nabla }_{xy}P+u_{z}b\frac{1}{\rho }{\partial }_{z}P\right)\cdot n\text{d}S}+{\delta }^2{u}_z^{-2}{F}_T\cdot \nabla u_z^2\\ -{\delta }^3{u}_z^{-1}{\partial }_t{u}_z{\partial }_z{u}_z{F}_z+u_z^{-2}{F}_z{\partial }_t{u}_z\\ ={\displaystyle {\int }_U{\partial }_t{u}_z b\cdot \left(b\otimes \nabla u_z\right)\text{d}V},\end{array}$ (97)

where now $b=u/\delta +w$ , so all $w$ and derivative contributions are explicit in the nonlinear term, in the surface integral, and inside each ${\Omega }_i$ via $b$ .

Remarks

1) Every occurrence of $b$ carries a contribution from $w$ . 2) Derivatives of $w$ enter through ${\partial }_{t}b$ , $\nabla \cdot b$ , ${\nabla }^{2}b$ . 3) This PDE generalizes Equation (19) to include non-aligned perturbations $w$ .

7.6. Lambert-W Velocity Field with Singular Scaling $\delta \left(x,t\right)$

Consider the decomposition

$u\left(x,t\right)=\frac{1}{\delta \left(x,t\right)}v\left(x,t\right)+w\left(x,t\right),$ (98)

where $w$ is smooth, say $w\in C^1$ , and $v$ carries the rough part of $u$ .

Suppose the velocity $u$ is given by a Lambert-W profile:

$u\left(x,t\right)=f\left(\Xi \left(x,t\right)\right):=W\left(-e^{-\Xi \left(x,t\right)-1}\right)+1,$ (99)

with $\Xi \left(x,t\right)$ affine in space and time.

Near the branch point $\Xi =0$ , the Lambert-W function satisfies

$f\left(\Xi \right)~\sqrt{2\Xi },f^{\prime }\left(\Xi \right)~\frac{1}{\sqrt{2\Xi }}.$ (100)

7.7. Effect of Singular $\delta \left(x,t\right)$

Let

$\delta \left(x,t\right)~\text{Ξ}{\left(x,t\right)}^{-p}\text{near the singularity},$ (101)

with $0<p<1$ .

Then, the rough component is

$v\left(x,t\right)=\delta \left(x,t\right)\left(u-w\right)~{\Xi }^{-p}\cdot \sqrt{\Xi }~{\Xi }^{\frac{1}{2}-p}.$ (102)

The vorticity of $u$ is

$\omega =\nabla \times u=\nabla \left(\frac{1}{\delta }\right)\times v+\frac{1}{\delta } \nabla \times v+\nabla \times w.$ (103)

7.8. Vorticity of $v$ and BKM [12] Condition

The magnitude of vorticity of $v$ scales as

$\left|v^{\prime }\left(\Xi \right)\right|\sim \left|\frac{\text{d}}{\text{d}\Xi }{\Xi }^{\frac{1}{2}-p}\right|=\left|\frac{1}{2}-p\right| {\Xi }^{-\frac{1}{2}-p}~{\Xi }^{-\frac{1}{2}-p}.$ (104)

8. Consistent Scaling, BKM Finiteness, and Onsager $C^{1/3}$

Scaling Ansatz

Let $\Xi$ denote the relevant small spatial scale (e.g., distance to a putative singular set). Assume the velocity scales as

$v\left(\Xi \right)~{\Xi }^{\frac{1}{2}-p},$ (105)

for some real parameter $p$ .

8.1. Derivative and Vorticity Scaling

Differentiating (105) with respect to $\Xi$ gives

$v^{\prime }\left(\Xi \right)=\frac{\text{d}}{\text{d}\Xi }{\Xi }^{\frac{1}{2}-p}=\left(\frac{1}{2}-p\right){\Xi }^{-\frac{1}{2}-p}.$ (106)

Hence, up to multiplicative constants,

$\left|{\omega }_v\right|~\left|v^{\prime }\right|~{\Xi }^{-\frac{1}{2}-p}.$ (107)

This is the correct and consistent scaling: differentiation lowers the exponent by one.

8.2. BKM Integral and Finiteness

The Beale-Kato-Majda criterion [12] requires

${\displaystyle {\int }_0^T{\Vert \omega \left(\cdot ,t\right)\Vert }_{L^{\infty }}\text{d}t}<\infty .$ (108)

Assume the smallest resolved scale behaves as

$\text{Ξ}\left(t\right)~{\left(T-t\right)}^{\gamma },\gamma >0,$

which is the standard self-similar scaling hypothesis near a singular time $T$ . Using (107),

${\Vert \omega \left(\cdot ,t\right)\Vert }_{L^{\infty }}~\text{Ξ}{\left(t\right)}^{-\frac{1}{2}-p}~{\left(T-t\right)}^{-\gamma \left(\frac{1}{2}+p\right)}.$

Thus, the BKM integral behaves like

${\displaystyle {\int }_0^T{\left(T-t\right)}^{-\gamma \left(\frac{1}{2}+p\right)}\text{d}t},$ (109)

which is finite if and only if

$\gamma \left(\frac{1}{2}+p\right)<1.$ (110)

In the natural case $\gamma =1$ (parabolic scaling), this reduces to

$\boxed{p<\frac{1}{2}.}$ (111)

Hence:

• $p<\frac{1}{2}$ : BKM integral finite,

• $p=\frac{1}{2}$ : borderline (logarithmic divergence),

• $p>\frac{1}{2}$ : BKM integral diverges.

Reconciliation with Onsager $C^{1/3}$

Onsager criticality requires velocity increments to scale as

$\left|{\delta }_{\Xi }v\right|~{\Xi }^{\frac{1}{3}}.$ (112)

Comparing with (105),

$\frac{1}{2}-p=\frac{1}{3}\Rightarrow \boxed{p=\frac{1}{6}.}$

For this value,

$\left|\omega \right|~{\Xi }^{-\frac{1}{2}-\frac{1}{6}}={\Xi }^{-\frac{2}{3}},$

and hence

${\Vert \omega \left(\cdot ,t\right)\Vert }_{L^{\infty }}~{\left(T-t\right)}^{-\frac{2}{3}},$

which satisfies

${\displaystyle {\int }_0^T{\left(T-t\right)}^{-\frac{2}{3}}\text{d}t}<\infty .$

Thus, Onsager-critical scaling lies strictly within the BKM-subcritical regime.

We have that,

$\boxed{\begin{array}{l}v\left(\Xi \right)~{\Xi }^{\frac{1}{2}-p}\Rightarrow \left|\omega \right|~{\Xi }^{-\frac{1}{2}-p},\\ \text{BKM finite}\iff p<\frac{1}{2},\\ \text{Onsager}{C}^{1/3}\iff p=\frac{1}{6}<\frac{1}{2}.\end{array}}$ (113)

Therefore, Onsager-critical velocity regularity is fully compatible with bounded vorticity growth in the sense of Beale-Kato-Majda, and strong non-alignment or intermittency does not by itself imply singularity formation.

What if the spheres where $\mathcal{X}=0$ satisfy more general conditions such that the spheres overlap?

9. Generalization of the Quadratic Invariant and Deep Overlap Analysis

9.1. The Original Special Case

For general functions $f_1,f_2$ and $f_3$ , suppose that

$f_1=f_2=f_3=\frac{1}{6},$

so that $f_1+f_2+f_3=1/2$ the quadratic invariant takes the form

$y_1^2+y_2^2+y_3^2-y_1-y_2-y_3+\frac{1}{2}=0.$ (114)

Completing the square,

$y_i^2-y_i={\left(y_i-\frac{1}{2}\right)}^2-\frac{1}{4},i=1,2,3,$

so that

${\displaystyle \sum }_{i=1}^3{\left(y_i-\frac{1}{2}\right)}^2=\frac{1}{4}.$

Hence, the invariant surface is a sphere with

$C=\left(\frac{1}{2},\frac{1}{2},\frac{1}{2}\right),r=\frac{1}{2}.$

9.2. General Quadratic Form for a Sphere

A general sphere in ${ℝ}^3$ can be written as

${\displaystyle \sum }_{i=1}^3{y}_i^2-2{\displaystyle \sum }_{i=1}^3{a}_i{y}_i+\gamma =0,$ (115)

which is equivalent to

${\displaystyle \sum }_{i=1}^3{\left(y_i-a_i\right)}^2={\displaystyle \sum }_{i=1}^3{a}_i^2-\gamma .$

Thus, the center and radius are

$C=\left(a_1,a_2,a_3\right),r=\sqrt{{\displaystyle \sum }_{i=1}^3{a}_i^2-\gamma }.$

9.3. Incorporating General Parameters $f_1,f_2,f_3$

To encode general parameters $f_1,f_2,f_3$ , we consider the quadratic invariant

$y_1^2+y_2^2+y_3^2-2f_1{y}_1-2f_2{y}_2-2f_3{y}_3+\gamma =0.$ (116)

Completing the square yields

${\displaystyle \sum }_{i=1}^3{\left(y_i-f_i\right)}^2={\displaystyle \sum }_{i=1}^3{f}_i^2-\gamma .$

Therefore,

$C=\left(f_1,f_2,f_3\right),r^2={\displaystyle \sum }_{i=1}^3{f}_i^2-\gamma .$

9.4. Prescribing the Radius $r=f_1+f_2+f_3$

If we require the radius to be

$r=f_1+f_2+f_3,$

Then, the constant term in (116) must satisfy

$\gamma ={\displaystyle \sum }_{i=1}^3{f}_i^2-{\left(f_1+f_2+f_3\right)}^2=-2\left(f_1{f}_2+f_1{f}_3+f_2{f}_3\right).$

Substituting into (116), we obtain the exact generalized sphere equation

$y_1^2+y_2^2+y_3^2-2f_1{y}_1-2f_2{y}_2-2f_3{y}_3-2\left(f_1{f}_2+f_1{f}_3+f_2{f}_3\right)=0,$ (117)

or equivalently,

${\displaystyle \sum }_{i=1}^3{\left(y_i-f_i\right)}^2={\left(f_1+f_2+f_3\right)}^2.$

9.5. Deep Overlap Analysis under the Generalized Radius Choice

Consider a collection of spheres

$S_n:=\left\{y\in {ℝ}^3:\Vert y-C_n\Vert =r_n\right\},$

with

$C_n=\left(f_1^{\left(n\right)},f_2^{\left(n\right)},f_3^{\left(n\right)}\right),r_n=f_1^{\left(n\right)}+f_2^{\left(n\right)}+f_3^{\left(n\right)}.$

Assume the parameters $f_i^{\left(n\right)}$ are chosen so that $C_n\in {\left[0,1\right]}^3$ .

9.5.1. Overlap Criterion

Two spheres with identical radii $r_n=r_m=r$ overlap if

$\Vert C_n-C_m\Vert <2r.$

More generally, for distinct radii,

$\Vert C_n-C_m\Vert <r_n+r_m.$

9.5.2. Uniform Deep Overlap Condition

Assume there exists $ϵ>0$ such that for all neighboring indices $n,m$ ,

$\Vert C_n-C_m\Vert \le \left(1-ϵ\right)\left(r_n+r_m\right).$

Then, the intersection

$B_{r_n}\left(C_n\right)\cap B_{r_m}\left(C_m\right)$

has nonempty interior, and in fact contains a ball of radius at least

${\rho }_{n,m}\ge \frac{1}{2}\left(r_n+r_m-\Vert C_n-C_m\Vert \right)\ge \frac{ϵ}{2}\left(r_n+r_m\right).$

This provides a quantitative lower bound on the overlap thickness. See Figure 6, which shows the overlap thickness.

Figure 6. Deep overlap of two spheres (2D cross-section). The shaded lens represents the nondegenerate overlap region ensured by $\epsilon >0$ .

9.5.3. Relation to the Cubic Partition

If the centers $C_n$ lie on a grid of spacing $\text{Δ}$ , then

$\Vert C_n-C_m\Vert \le \sqrt{3}\text{Δ}\text{for neighboring centers}.$

Thus, a sufficient condition for uniform deep overlap is

$f_1^{\left(n\right)}+f_2^{\left(n\right)}+f_3^{\left(n\right)}\ge \frac{\sqrt{3}}{2} \text{Δ}\left(1+\epsilon \right),\epsilon >0.$

Under this condition, the generalized-radius spheres satisfy the same deep overlap property as in the fixed-radius construction. The idea here is that when there are deep overlap stronger radial components of vorticity appear as opposed to purely tangential components for which the BKM analysis was previously based on. Having both radial and tangential components may raise the question of finite time blowup and its possibility. As a conjecture, there is no finite time blowup in this case. In the sequel to this paper, BKM analysis and specific inequalities will be used to show that singularity formation is not possible under overlapping spheres assumption.

10. Conclusion

In this work, a geometric and analytic framework has been developed for the three-dimensional incompressible Navier-Stokes equations on the torus ${\mathbb{T}}^3$ , with particular emphasis on the role of nonlinear interaction geometry, periodicity, and singular structure. By reformulating the Navier-Stokes system through a componentwise and geometrically aligned decomposition, the analysis isolates a class of solutions whose evolution is governed by implicit relations naturally inverted by the Lambert $W$ function. This approach departs from purely functional analytic existence theories by resolving the internal structure of the nonlinear term and identifying a rigid normal form that remains closed under differentiation, multiplication, and nonlinear interaction. A central result of the paper is the explicit construction of Lambert $W$ -based solution profiles that capture both singular and non-blowup behavior within the same analytical framework. The appearance of singularities is shown to be intrinsically linked to the elliptic degeneracy of the inverse Weierstrass $\wp$ -function and its critical values, rather than to uncontrolled growth in the Navier-Stokes nonlinearity itself. Through careful analysis of the associated Poisson equation and the induced coefficient $\kappa \left(x,y,z,t\right)$ , it is demonstrated that these singularities occur on geometrically structured sets—specifically, spherical surfaces embedded in the torus and are governed by invariant elliptic data rather than spacetime-dependent instabilities. A key contribution of the paper is the demonstration that such singularities are movable rather than intrinsic. By employing fixed-point methods and iterated compositions of the Lambert $W$ function, the singular sets can be pushed outward in space as the underlying torus is expanded. In the limit of increasingly large tori approaching ${ℝ}^3$ , the singularities are transported to infinity, yielding globally regular behavior on any fixed bounded spatial domain. This mechanism provides a concrete realization of how finite-space regularity may coexist with formally singular solution representations, and clarifies the distinction between genuine finite-time blowup and geometric singularities tied to compactification. The construction of strictly periodic velocity fields is achieved by embedding the Lambert $W$ dynamics into a Weierstrass elliptic lattice and compensating for the quasi-periodicity of the zeta function via an explicit linear drift term. This drift is shown to be physically admissible, as it can be absorbed into a time-dependent pressure gauge without affecting the velocity field. The resulting solutions are bounded, smooth, and periodic in both space and time, providing explicit examples of nontrivial periodic attractor-type behavior for the three-dimensional Navier-Stokes equations on ${\mathbb{T}}^3$ . Beyond the specific constructions presented, the analysis supports a broader conjectural picture: that Lambert $W$ -generated profiles form invariant manifolds of finite codimension for the Navier-Stokes flow and represent a universal normal form for locally self-similar, renormalization-invariant solutions. While this work does not claim a complete resolution of the Navier-Stokes regularity problem, it establishes a concrete algebraic-geometric mechanism by which singular structures can be controlled, relocated, and neutralized through composition and scaling. The results suggest that explicit transcendental structures, rather than being pathological, may lie at the core of the Navier-Stokes dynamics and offer a promising avenue for further analytical classification of solutions. Future work will focus on strengthening the invariant-manifold interpretation, clarifying the stability properties of the periodic attractor constructed here, and extending the framework to broader classes of forcing and boundary conditions. Also, BKM analysis and inequalities will be used to show that singularity formation is not possible under overlapping spheres assumption. The interplay between elliptic geometry, implicit functional inversion, and nonlinear fluid dynamics highlighted in this paper points toward a unified geometric calculus for periodic Navier-Stokes flows and provides a new perspective on the long-standing challenges of singularity formation and global regularity.

Acknowledgements

I would like to thank the reviewers at APM journal for their helpful comments during the peer review process.

Appendix A

In [1], Equation (6) is repeated here,

$\begin{array}{l}{u}_z\frac{\partial a}{\partial t}+a\cdot \nabla a-\frac{\mu }{\rho }{u}_z^2{\nabla }^{2}b+\frac{1}{\delta \rho }{u}_z^2{\nabla }_{xy}P+b\frac{1}{\rho }{u}_z\frac{\partial P}{\partial z}\\ +{\delta }^2{u}_z^2{F}_T+{\delta }^3{u}_z{F}_{z}b+u_{z}b\nabla \cdot a-\frac{\mu }{\rho }{u}_{z}b{\nabla }^2{u}_z=0\end{array}$ (118)

The highest derivative is in $t$ in the expression above and is linked to $u_{z}b\nabla \cdot a$ , (the term that was added to give Equation (6)), hence the PDE to solve is:

$u_z\frac{\partial a}{\partial t}+u_{z}b\nabla \cdot a=0$ (119)

Since we are dealing with an incompressible fluid flow we have $\nabla \cdot b=0$ , so that,

$u_{z}b\nabla \cdot a=b^2{u}_z\nabla u_z$ (120)

where $b$ is one component of $b=\left(u_x{u}_z,u_y{u}_z,u_x{u}_y\right)$ .

Now setting $u_z=\text{LambertW}(-\text{exp}\left(x-y+2z-t-1\right)+1$ in the previous PDE, gives the PDE in terms of $b=b\left(x,y,z,t\right)$ as,

$\begin{array}{l}\left(W\left(-e^{-x+y-2z-t-1}\right)+1\right)\left(-\frac{W\left(-e^{-x+y-2z-t-1}\right)b}{W\left(-e^{-x+y-2z-t-1}\right)+1}+\left(W\left(-e^{-x+y-2z-t-1}\right)+1\right)\frac{\partial b}{\partial t}\right)\\ =-2W\left(-e^{-x+y-2z-t-1}\right)b^2\end{array}$ (121)

$b\left(x,y,z,t\right)={\left(\left(F_1\left(x,y,z\right)-1\right)W\left(-e^{-x+y-2z-t-1}\right)+F_1\left(x,y,z\right)+1\right)}^{-1}$ (122)

Thus, we have,

$u_x{u}_y{u}_z=\frac{W\left(-e^{-x+y-2z-t-1}\right)+1}{\left(F_1\left(x,y,z\right)-1\right)W\left(-e^{-x+y-2z-t-1}\right)+F_1\left(x,y,z\right)+1}$ (123)

Since $u_z=W\left(-e^{-x+y-2z-t-1}\right)+1$ . The special case of $F_1\left(x,y,z\right)=0$ gives the required result $u_x{u}_y{u}_z=-\frac{W+1}{W-1}$ .

Appendix B: Periodic Compensation for the Weierstrass Zeta Function at $g_2=4$ , $g_3=0$

1. The Elliptic Curve and Its Roots

The Weierstrass elliptic function satisfies

${\left({\wp }^{\prime }\right)}^2=4{\wp }^3-g_2\wp -g_3.$

For

$g_2=4,g_3=0,$

this reduces to

${\left({\wp }^{\prime }\right)}^2=4\left({\wp }^3-\wp \right).$

Factoring,

${\wp }^3-\wp =\wp \left(\wp -1\right)\left(\wp +1\right),$

So the roots are

$e_1=1,e_2=0,e_3=-1.$

The discriminant

$\text{Δ}=g_2^3-27g_3^2=64\ne 0$

confirms that the curve is nondegenerate and that the associated lattice is rectangular (the lemniscatic case).

2. Derivation of the Real Half-Period ${\Omega }_1$

The real half-period is defined by

${\Omega }_1={\displaystyle {\int }_{e_1}^{\infty }\frac{\text{d}t}{\sqrt{4\left(t-e_1\right)\left(t-e_2\right)\left(t-e_3\right)}}}.$

Substituting the roots,

${\Omega }_1={\displaystyle {\int }_1^{\infty }\frac{\text{d}t}{\sqrt{4\left(t-1\right)t\left(t+1\right)}}}=\frac{1}{2}{\displaystyle {\int }_1^{\infty }\frac{\text{d}t}{\sqrt{\left(t-1\right)t\left(t+1\right)}}}.$

Lemniscatic Substitution

Let

$t=\frac{1}{\sin \theta },\theta \in \left(0,\frac{\pi }{2}\right).$

Then

$\text{d}t=-\frac{\cos \theta }{{\sin }^2\theta }\text{d}\theta ,\left(t-1\right)t\left(t+1\right)=\frac{{\cos }^2\theta }{{\sin }^3\theta }.$

Hence

$\sqrt{\left(t-1\right)t\left(t+1\right)}=\frac{\cos \theta }{{\sin }^{3/2}\theta }.$

Substituting,

${\Omega }_1=\frac{1}{2}{\displaystyle {\int }_0^{\pi /2}{\sin }^{1/2}\theta \text{d}\theta }.$

Evaluation

Using the identity

${\displaystyle {\int }_0^{\pi /2}{\sin }^p\theta \text{d}\theta }=\frac{\sqrt{\pi }\Gamma \left(\frac{p+1}{2}\right)}{2\Gamma \left(\frac{p+2}{2}\right)},$

with $p=\frac{1}{2}$ , we obtain

${\Omega }_1=\frac{\sqrt{\pi }}{4}\frac{\Gamma \left(\frac{3}{4}\right)}{\Gamma \left(\frac{5}{4}\right)}.$

Numerically,

${\Omega }_1\approx \mathrm{0.927039.}$

3. Derivation of the Zeta Jump $\eta$

The zeta constant is defined by

$\eta :=\zeta \left({\Omega }_1\right).$

For a rectangular lattice, $\eta$ admits the integral representation

$\eta ={\displaystyle {\int }_{e_1}^{\infty }\left(\frac{t}{\sqrt{4\left(t-e_1\right)\left(t-e_2\right)\left(t-e_3\right)}}-\frac{1}{2t^{3/2}}\right)\text{d}t}.$

Substituting the roots,

$\eta ={\displaystyle {\int }_1^{\infty }\left(\frac{t}{\sqrt{4\left(t-1\right)t\left(t+1\right)}}-\frac{1}{2t^{3/2}}\right)\text{d}t}.$

This convergent integral evaluates to the closed form

$\eta =\frac{\sqrt{\pi }}{4}\frac{\Gamma \left(\frac{1}{4}\right)}{\Gamma \left(\frac{3}{4}\right)}.$

Numerically,

$\eta \approx \mathrm{2.6232.}$

4. Quasi-Periodicity of $\zeta \left(z;4,0\right)$

The Weierstrass zeta function satisfies the quasi-periodicity relation

$\zeta \left(z+2{\Omega }_1;4,0\right)=\zeta \left(z;4,0\right)+2\eta .$

This linear drift is universal for elliptic lattices.

5. Periodic Compensation

Define the compensated function

$F\left(z\right)=\zeta \left(z;4,0\right)-\frac{\eta }{{\Omega }_1}z.$

Then

$\begin{array}{c}F\left(z+2{\Omega }_1\right)=\zeta \left(z+2{\Omega }_1\right)-\frac{\eta }{{\Omega }_1}\left(z+2{\Omega }_1\right)=\left[\zeta \left(z\right)+2\eta \right]-\frac{\eta }{{\Omega }_1}z-2\eta \\ =\zeta \left(z\right)-\frac{\eta }{{\Omega }_1}z=F\left(z\right).\end{array}$

6. Final Result

$F\left(z+2{\Omega }_1\right)=F\left(z\right)$

with

${\Omega }_1=\frac{\sqrt{\pi }}{4}\frac{\Gamma \left(\frac{3}{4}\right)}{\Gamma \left(\frac{5}{4}\right)},\eta =\frac{\sqrt{\pi }}{4}\frac{\Gamma \left(\frac{1}{4}\right)}{\Gamma \left(\frac{3}{4}\right)}.$

Thus, although $\zeta \left(z;4,0\right)$ is quasi-periodic, the uniquely determined linear compensation removes the drift and yields a strictly periodic function. This is the canonical mechanism underlying periodic velocity potentials constructed from the Weierstrass zeta function. Here is a Maple code that accomplishes this and produces a periodic function,

7. Maple Code for Drift-Free Periodic $V_z\left(x\right)$