Box Stokes: Theory and Applications

• Box Stokes is a multifaceted topic that unifies methods from singular set analysis in Navier–Stokes equations to finite-volume and formal cubical approaches.
• It employs box-based geometric scaffolds—such as Voronoi dual meshes and axis-aligned cubes—to derive regularity bounds, stabilize discretizations, and formalize Stokes’ theorem.
• Applications span improved singularity dimension estimates, first-order accurate numerical schemes, and scalable periodic boundary-integral solvers for complex fluid dynamics.

“Box Stokes” names several technically distinct constructions in the contemporary literature. In the cited arXiv corpus, the phrase connects box-counting and parabolic covering arguments for singular sets of the 3D Navier–Stokes equations, Box-Method discretizations of the Stokes problem on Voronoi dual meshes, the box Stokes theorem for axis-aligned cubes and smooth singular cubical chains, oscillatory flows in a fully oscillating box where Stokes-layer effects govern particle-pair dynamics, and periodic-box boundary-integral formulations for Stokes flow (Wang et al., 2016, Wang et al., 2018, Negrini et al., 2023, Hulak et al., 1 May 2026, Overveld et al., 2022, Li et al., 29 May 2026).

Context	Meaning of “box”	Stokes/Navier–Stokes role
Regularity theory	box-counting cover by parabolic cylinders	singular-set dimension
Finite-volume/Petrov–Galerkin discretization	Voronoi dual “Box mesh”	Stokes problem
Differential topology and formalization	axis-aligned cube or box	Stokes theorem
Oscillatory viscous flow	closed oscillating box	Stokes boundary layer
Boundary-integral periodization	periodic cell box	periodic Stokes flow

1. Box-counting singularity theory for suitable weak solutions

For a bounded set $X\subset\mathbb{R}^3\times\mathbb{R}$, let $N(X,\varepsilon)$ denote the minimal number of closed parabolic cylinders of radius $\varepsilon$ needed to cover $X$. The upper box-counting dimension is

$$\overline{\dim}_B(X) = \limsup_{\varepsilon\to0} \frac{\log N(X,\varepsilon)}{-\log\varepsilon}.$$

In Wang–Wu’s analysis of suitable weak solutions $(u,\Pi)$ to the 3D Navier–Stokes system, the possible singular set $\mathcal S$ consists of space-time points $(x_0,t_0)$ at which $(u,\Pi)$ can fail to be locally bounded, and the main theorem is

$$\overline{\dim}_B(\mathcal S)\le \frac{135}{104}\approx 1.30.$$

This improves the earlier bounds $N(X,\varepsilon)$0, $N(X,\varepsilon)$1, and $N(X,\varepsilon)$2 (Wang et al., 2016).

The underlying solution class is the standard suitable weak-solution framework: $N(X,\varepsilon)$3 with distributional satisfaction of Navier–Stokes and the local energy inequality

$N(X,\varepsilon)$4

The proof has two parts. The first is an $N(X,\varepsilon)$5-regularity criterion with improved scaling: there exists $N(X,\varepsilon)$6 such that if, on a parabolic cylinder $N(X,\varepsilon)$7,

$N(X,\varepsilon)$8

then $N(X,\varepsilon)$9 is regular. The decisive ingredient is the pressure treatment: $\varepsilon$0 with $\varepsilon$1 in $\varepsilon$2, combined with scaling-invariant estimates and explicit use of $\varepsilon$3. This is coupled to an interpolation–Gagliardo–Nirenberg inequality,

$\varepsilon$4

and to the pressure decay estimate

$\varepsilon$5

The second part is the covering argument. If $\varepsilon$6 with $\varepsilon$7, then for scales $\varepsilon$8 there exist

$\varepsilon$9

disjoint cylinders $X$0, each centered at a singular point. On each cylinder the scaled norm must exceed $X$1, so

$X$2

Since $X$3, the left-hand side diverges, a contradiction. The result is a quantitative restriction on the fractal size of the possible singular set.

2. One-scale $X$4-regularity and the bounds $X$5 and $X$6

A later refinement considers both interior and boundary suitable weak solutions of the 3D Navier–Stokes equations. For interior singular points $X$7 and boundary singular points $X$8, the main bounds are

$X$9

The relevant parabolic sets are

$$\overline{\dim}_B(X) = \limsup_{\varepsilon\to0} \frac{\log N(X,\varepsilon)}{-\log\varepsilon}.$$0

with $$\overline{\dim}_B(X) = \limsup_{\varepsilon\to0} \frac{\log N(X,\varepsilon)}{-\log\varepsilon}.$$1 in the boundary case (Wang et al., 2018).

The proof is organized around one-scale $$\overline{\dim}_B(X) = \limsup_{\varepsilon\to0} \frac{\log N(X,\varepsilon)}{-\log\varepsilon}.$$2-regularity criteria. In the interior case, there exist $$\overline{\dim}_B(X) = \limsup_{\varepsilon\to0} \frac{\log N(X,\varepsilon)}{-\log\varepsilon}.$$3 and any $$\overline{\dim}_B(X) = \limsup_{\varepsilon\to0} \frac{\log N(X,\varepsilon)}{-\log\varepsilon}.$$4 such that, if for some $$\overline{\dim}_B(X) = \limsup_{\varepsilon\to0} \frac{\log N(X,\varepsilon)}{-\log\varepsilon}.$$5,

$$\overline{\dim}_B(X) = \limsup_{\varepsilon\to0} \frac{\log N(X,\varepsilon)}{-\log\varepsilon}.$$6

then $$\overline{\dim}_B(X) = \limsup_{\varepsilon\to0} \frac{\log N(X,\varepsilon)}{-\log\varepsilon}.$$7 is regular. In the boundary case there is an analogous criterion with right-hand side $$\overline{\dim}_B(X) = \limsup_{\varepsilon\to0} \frac{\log N(X,\varepsilon)}{-\log\varepsilon}.$$8. These criteria are paired with higher-scale regularity tests of Guevara–Phuc type, including

$$\overline{\dim}_B(X) = \limsup_{\varepsilon\to0} \frac{\log N(X,\varepsilon)}{-\log\varepsilon}.$$9

for the interior problem, and a boundary version using mixed $(u,\Pi)$0 bounds for $(u,\Pi)$1 and $(u,\Pi)$2.

Two scale-transfer estimates are central. The interpolation lemma states that for $(u,\Pi)$3 and $(u,\Pi)$4,

$(u,\Pi)$5

while the pressure-decay lemma gives, for $(u,\Pi)$6,

$(u,\Pi)$7

Analogous boundary lemmas supply the same mechanism near $(u,\Pi)$8.

The covering argument is standard in form but sharper in exponent balance. If $(u,\Pi)$9 exceeded $\mathcal S$0, with $\mathcal S$1 in the interior case or $\mathcal S$2 at the boundary, then at arbitrarily small scales $\mathcal S$3 one could find $\mathcal S$4 disjoint parabolic balls centered at singular points. On each such ball the one-scale criterion fails, so

$\mathcal S$5

with $\mathcal S$6 or $\mathcal S$7. Summing over the disjoint family gives $\mathcal S$8, hence $\mathcal S$9. Choosing $(x_0,t_0)$0 yields the interior exponent $(x_0,t_0)$1.

In the comparison recorded with the result, the interior bound $(x_0,t_0)$2 strictly improves earlier values $(x_0,t_0)$3, $(x_0,t_0)$4, $(x_0,t_0)$5, $(x_0,t_0)$6, while the boundary bound $(x_0,t_0)$7 improves the previous $(x_0,t_0)$8. The reduction of the box dimension is a quantitative statement that any blow-up set, if it exists, occupies a thinner subset of space-time.

3. Hyperdissipative Navier–Stokes and box dimension at blow-up times

For the hyperdissipative system

$(x_0,t_0)$9

with $(u,\Pi)$0, the relevant weak solutions are Leray–Hopf solutions

$(u,\Pi)$1

satisfying the weak formulation and the global energy inequality

$(u,\Pi)$2

The singular set is now spatial: if $(u,\Pi)$3 is smooth on disjoint time intervals $(u,\Pi)$4, then $(u,\Pi)$5 consists of points singular at the blow-up time $(u,\Pi)$6, and $(u,\Pi)$7. The resulting estimates are

$(u,\Pi)$8

These bounds are derived without the suitable-solution local energy inequality machinery (Ożański, 2020).

The analysis is frequency-localized. For dyadic frequency $(u,\Pi)$9, spatial cube $$\overline{\dim}_B(\mathcal S)\le \frac{135}{104}\approx 1.30.$$0 of side $$\overline{\dim}_B(\mathcal S)\le \frac{135}{104}\approx 1.30.$$1, and cutoff $$\overline{\dim}_B(\mathcal S)\le \frac{135}{104}\approx 1.30.$$2, define

$$\overline{\dim}_B(\mathcal S)\le \frac{135}{104}\approx 1.30.$$3

A differential inequality for $$\overline{\dim}_B(\mathcal S)\le \frac{135}{104}\approx 1.30.$$4 isolates the dissipative term $$\overline{\dim}_B(\mathcal S)\le \frac{135}{104}\approx 1.30.$$5 and decomposes the nonlinear remainder into low-to-high, local, and high-to-high interactions. This is the analogue of a local energy inequality in dyadic form.

The geometric input is a good/bad cube decomposition. At scale $$\overline{\dim}_B(\mathcal S)\le \frac{135}{104}\approx 1.30.$$6, $$\overline{\dim}_B(\mathcal S)\le \frac{135}{104}\approx 1.30.$$7 is partitioned into cubes of side $$\overline{\dim}_B(\mathcal S)\le \frac{135}{104}\approx 1.30.$$8. A $$\overline{\dim}_B(\mathcal S)\le \frac{135}{104}\approx 1.30.$$9-cube $N(X,\varepsilon)$00 is $N(X,\varepsilon)$01-good if

$N(X,\varepsilon)$02

otherwise it is $N(X,\varepsilon)$03-bad. If a $N(X,\varepsilon)$04-cube and all its dyadic ancestors back to scale $N(X,\varepsilon)$05 are good, then

$N(X,\varepsilon)$06

uniformly up to blow-up. Any point outside infinitely many bad cubes therefore enjoys slightly supercritical decay.

A Vitali-type argument constructs coverings $N(X,\varepsilon)$07 of the bad $N(X,\varepsilon)$08-cubes with

$N(X,\varepsilon)$09

This yields $N(X,\varepsilon)$10. For box dimension, a more refined finite-generation covering is required. A first estimate uses scales $N(X,\varepsilon)$11, but the refined argument shows that it is enough to use $N(X,\varepsilon)$12, giving

$N(X,\varepsilon)$13

and therefore

$N(X,\varepsilon)$14

When $N(X,\varepsilon)$15, the formula gives $N(X,\varepsilon)$16, recovering the classical Scheffer–CKN box-counting bound. As $N(X,\varepsilon)$17, the bound decreases to zero, matching the regularizing effect of the critical hyperdissipation threshold. The method situates box dimension inside a frequency-space partial regularity program based on Littlewood–Paley packets, Bony paraproducts, and commutator estimates.

4. The Rhie–Chow stabilized Box Method for the Stokes problem

In numerical analysis, “Box Method” refers to a piecewise linear Petrov–Galerkin formulation on the Voronoi dual mesh of a Delaunay triangulation. Let $N(X,\varepsilon)$18, $N(X,\varepsilon)$19, be polygonal or polyhedral, and let $N(X,\varepsilon)$20 be a conforming Delaunay triangulation. Its Voronoi dual $N(X,\varepsilon)$21 is the Box mesh, with one box attached to each interior vertex $N(X,\varepsilon)$22. The trial spaces are

$N(X,\varepsilon)$23

while the pressure test space is the boxwise-constant space $N(X,\varepsilon)$24, accessed through the lumping map $N(X,\varepsilon)$25, $N(X,\varepsilon)$26 for $N(X,\varepsilon)$27. The resulting scheme is the Rhie–Chow stabilized Box Method (RCBM) for the Stokes problem (Negrini et al., 2023).

The unstabilized auxiliary formulation uses the bilinear forms

$N(X,\varepsilon)$28

$N(X,\varepsilon)$29

The associated finite-volume-like momentum balance on each box $N(X,\varepsilon)$30 is

$N(X,\varepsilon)$31

and continuity is imposed through

$N(X,\varepsilon)$32

To suppress spurious pressure modes on the co-located grid, continuity is augmented by a Rhie–Chow stabilization

$N(X,\varepsilon)$33

Algebraically,

$N(X,\varepsilon)$34

is a consistent pressure Laplacian, while

$N(X,\varepsilon)$35

is the off-diagonal term induced by $N(X,\varepsilon)$36 acting on the pressure; their difference recovers a stable Schur complement. This is the discrete analogue of checkerboard-pressure suppression.

The analysis introduces the norm

$N(X,\varepsilon)$37

with

$N(X,\varepsilon)$38

Under standard mesh regularity and the properties

$N(X,\varepsilon)$39

$N(X,\varepsilon)$40

and the generalized inf–sup condition

$N(X,\varepsilon)$41

a Babuška–Brezzi argument yields discrete inf–sup stability, coercivity, existence, and uniqueness.

The a priori estimate proved under $N(X,\varepsilon)$42, $N(X,\varepsilon)$43 is

$N(X,\varepsilon)$44

Manufactured-solution tests in 2D and 3D show approximately first-order convergence in both the $N(X,\varepsilon)$45 velocity norm and the $N(X,\varepsilon)$46 pressure norm. The method is described as locally conservative, valid on arbitrary polygonal or polyhedral meshes, and easily implementable in industrial solvers such as OpenFOAM. The stated limitations are equally specific: only first-order convergence in the present lowest-order setting, the need for upwind variants for faster advection-dominated flows, and theoretical coercivity and inf–sup proofs that remain based on numerically validated conjectures.

5. Box Stokes theorem, cubical chains, and formal verification in Lean 4

In geometric and formalized mathematics, a box is an axis-aligned product of intervals

$N(X,\varepsilon)$47

with unit cube

$N(X,\varepsilon)$48

Each box can be regarded as a singular $N(X,\varepsilon)$49-cube through the affine inclusion

$N(X,\varepsilon)$50

Smooth singular chains in $N(X,\varepsilon)$51 are finitely supported $N(X,\varepsilon)$52-linear combinations of such $N(X,\varepsilon)$53 parametrizations, implemented as

$N(X,\varepsilon)$54

This is the setting of the Lean 4 formalization of Stokes’ theorem for smooth singular cubes (Hulak et al., 1 May 2026).

The cubical boundary operator on a single smooth singular cube $N(X,\varepsilon)$55 is

$N(X,\varepsilon)$56

where $N(X,\varepsilon)$57 inserts the constant $N(X,\varepsilon)$58 in the $N(X,\varepsilon)$59-th coordinate. Extending $N(X,\varepsilon)$60-linearly gives the chain-complex differential. The identity $N(X,\varepsilon)$61 is verified by a sign-reversing involution on double-face indices,

$N(X,\varepsilon)$62

which pairs the contributions in the double boundary.

The box Stokes theorem itself is stated on $N(X,\varepsilon)$63 for an $N(X,\varepsilon)$64-form

$N(X,\varepsilon)$65

as

$N(X,\varepsilon)$66

The boundary integral is the alternating sum

$N(X,\varepsilon)$67

and the interior integral is the Bochner integral of the top-degree coefficient of $N(X,\varepsilon)$68. In the formalization, pullback is defined using the Fréchet derivative: $N(X,\varepsilon)$69 Applying box Stokes to $N(X,\varepsilon)$70 on $N(X,\varepsilon)$71 yields chain-level Stokes for singular cubes.

The development also records the compatibility of $N(X,\varepsilon)$72 with the nilpotency of the exterior derivative,

$N(X,\varepsilon)$73

through the chain-level identity

$N(X,\varepsilon)$74

This identifies integration as a cochain map from the cubical chain complex to the de Rham complex. In dimension $N(X,\varepsilon)$75, the formalized theorem specializes to Green’s theorem on the unit square $N(X,\varepsilon)$76, with the oriented boundary decomposition

$N(X,\varepsilon)$77

6. Oscillating boxes, Stokes layers, and particle-pair dynamics

In viscous oscillatory-flow mechanics, the relevant box is a closed container whose top and bottom walls oscillate horizontally. The dimensional incompressible Navier–Stokes equations are

$N(X,\varepsilon)$78

with wall conditions

$N(X,\varepsilon)$79

Using

$N(X,\varepsilon)$80

one obtains the nondimensional form with viscous length scale

$N(X,\varepsilon)$81

The derivation of $N(X,\varepsilon)$82 comes from Stokes’s second problem: for $N(X,\varepsilon)$83, substitution into $N(X,\varepsilon)$84 gives

$N(X,\varepsilon)$85

hence

$N(X,\varepsilon)$86

This is the Stokes boundary-layer thickness (Overveld et al., 2022).

The nondimensional parameters introduced in the study are

$N(X,\varepsilon)$87

where $N(X,\varepsilon)$88, $N(X,\varepsilon)$89, and $N(X,\varepsilon)$90 or sometimes $N(X,\varepsilon)$91. The central equivalence statement is that the oscillating box and the oscillating channel become equivalent only when

$N(X,\varepsilon)$92

equivalently $N(X,\varepsilon)$93 and $N(X,\varepsilon)$94. The reported numerical criterion for excellent collapse is

$N(X,\varepsilon)$95

The direct numerical simulations use a uniform Cartesian grid with $N(X,\varepsilon)$96, domain size $N(X,\varepsilon)$97, periodic boundary conditions in $N(X,\varepsilon)$98, no-slip at $N(X,\varepsilon)$99 or symmetry at $\varepsilon$00, about $\varepsilon$01 Lagrangian markers per sphere, explicit three-stage Runge–Kutta time advancement in a projection-method pressure-corrector, and a soft-sphere collision model with lubrication correction.

For two identical spheres placed side-by-side in the oscillating box, the nonlinear convective term generates a nonzero time-averaged steady-streaming flow of order $\varepsilon$02. For $\varepsilon$03, this streaming is organized into two half-vortex rings on the upstream and downstream sides of each sphere. In the spanwise gap, overlap of the two streaming patterns produces a net inward mean flow. Equilibrium occurs when the viscous repulsion from thin vorticity layers of thickness $\varepsilon$04 balances the advective attraction from streaming rings of size $\varepsilon$05.

The resulting “Box Stokes” scaling laws are

$\varepsilon$06

$\varepsilon$07

for the nondimensional mean gap $\varepsilon$08, and

$\varepsilon$09

for the magnitude of the steady streaming. The study further states that in a fully oscillating box no Stokes boundary layer modifies the streaming flow, so the only relevant parameters are $\varepsilon$10 and the relative excursion $\varepsilon$11. Outside the regime $\varepsilon$12 and $\varepsilon$13, box and channel dynamics diverge because the oscillating channel has a bottom-generated Stokes layer that distorts the half-vortex rings.

7. Periodic-box boundary-integral formulations for Stokes flow

A further usage of “box” appears in scalable boundary-integral solvers for periodic Stokes flow. The governing free-space representation uses the Stokeslet, stresslet, and rotlet kernels,

$\varepsilon$14

$\varepsilon$15

with $\varepsilon$16, $\varepsilon$17, and associated layer potentials $\varepsilon$18, $\varepsilon$19, and $\varepsilon$20. The periodic geometry is a box of side $\varepsilon$21, and the key algorithmic object is a one-time precomputed operator that maps outgoing proxy strengths to incoming proxy strengths for all far periodic images (Li et al., 29 May 2026).

The auxiliary basis comes from kernel-independent FMM equivalent surfaces. Proxy sources are placed on an outgoing equivalent surface $\varepsilon$22 surrounding the central box, and strengths $\varepsilon$23 are determined by matching the field at check points. Triply periodicity is imposed by splitting the infinite lattice into the central box plus its $\varepsilon$24 neighbors and a far field comprising all remaining images. The level-0 multipole-to-local operator is

$\varepsilon$25

and higher levels are generated by the scaling

$\varepsilon$26

Absolute convergence requires the zero net-force compatibility condition

$\varepsilon$27

implemented through a rank-one projector that removes any nonzero monopole component. The final dense matrix $\varepsilon$28 depends only on the periodic-box geometry, the periodicity dimension, the proxy and check surfaces, and the KIFMM multipole order $\varepsilon$29; it does not depend on the embedded surfaces $\varepsilon$30 or on the physical density $\varepsilon$31. The same $\varepsilon$32 is reused verbatim across the Stokeslet, stresslet, and rotlet.

Each GMRES application of the periodized operator consists of four stages: an FMM upward pass to compute outgoing proxy strengths, a dense multiplication $\varepsilon$33, an FMM downward pass with modified near-field lists including the $\varepsilon$34 images, and local Nyström corrections. The complexity is

$\varepsilon$35

because $\varepsilon$36 is fixed by the requested accuracy. Only one near-field layer of image boxes is needed explicitly.

The reported numerical performance includes spectral convergence in the azimuthal Fourier modes $\varepsilon$37, high-order convergence in the number of panels $\varepsilon$38, and relative errors $\varepsilon$39–$\varepsilon$40 with $\varepsilon$41–$\varepsilon$42 panels and $\varepsilon$43–$\varepsilon$44 Fourier modes. The hierarchical far-field sum converges exponentially in both $\varepsilon$45 and the number of levels: $\varepsilon$46 levels suffice to reduce the error below $\varepsilon$47 for $\varepsilon$48, while $\varepsilon$49 levels were used to saturate machine precision. On a single $\varepsilon$50-core node, adding triply periodicity increased setup time by $\varepsilon$51–$\varepsilon$52 and per-apply evaluation time by $\varepsilon$53–$\varepsilon$54. On up to $\varepsilon$55 cores, systems with up to $\varepsilon$56 unknowns exhibited $\varepsilon$57 weak-scaling efficiency, while strong scaling from $\varepsilon$58 to $\varepsilon$59 cores yielded $\varepsilon$60 efficiency at $\varepsilon$61 cores and per-iteration wall times below $\varepsilon$62 for a single boundary-integral evaluation on roughly $\varepsilon$63 million unknowns.

Taken together, these usages show that “Box Stokes” is not a single theorem or method. In analysis it denotes box-counting restrictions on singular sets; in discretization it denotes a Voronoi-dual Petrov–Galerkin or finite-volume formulation; in formalized geometry it is Stokes’ theorem on axis-aligned cubes; in oscillatory-flow mechanics it is the clean box realization of steady streaming; and in periodic boundary-integral computation it is a box-geometry framework for $\varepsilon$64 Stokes solvers. The common structural motif is the replacement of a continuous domain or singular set by a box-based geometric scaffold that makes regularity, conservation, topology, or long-range hydrodynamic interaction quantitatively tractable.