1. Introduction
Direct numerical simulation (DNS) (Orszag Reference Orszag1970; Rogallo Reference Rogallo1981; She et al. Reference She, Jackson and Orszag1990; Moin & Mahesh Reference Moin and Mahesh1998; Coleman & Sandberg Reference Coleman and Sandberg2010; Alexakis et al. Reference Alexakis, Marino, Mininni, van Kan, Foldes and Feraco2024) has been widely applied to solve Navier–Stokes (NS) turbulence, i.e. turbulence governed by NS equations (Navier Reference Navier1822; Stokes Reference Stokes1845). However, given that exact solutions of NS turbulence are unknown, there exists a reliability crisis in turbulence simulation: it is unknown whether numerical results are true solutions in mathematics or just artefacts caused by numerical noise.
Currently, many researchers (Deissler Reference Deissler1986; Boffetta & Musacchio Reference Boffetta and Musacchio2001, Reference Boffetta and Musacchio2017; Berera & Ho Reference Berera and Ho2018; Ge et al. Reference Ge, Rolland and Vassilicos2023; Okamoto et al. Reference Okamoto, Ishihara, Yokokawa and Kaneda2025) have pointed out that NS turbulence is chaotic in that the spatiotemporal trajectories of NS turbulence are very sensitive to small disturbances. Numerical algorithms had to be used because the exact solutions of chaotic systems are unknown. However, all numerical algorithms contain numerical noise, in the form of truncation and round-off errors. Lorenz (Reference Lorenz2006) further discovered that numerical noise might exert distinct influences on the macroscopic properties of chaotic systems. Invariably, the traditional algorithms of DNS for NS turbulence introduce numerical error, which may increase exponentially to a macro-level in a chaotic flow. This suggests that the DNS result of NS turbulence should rapidly be badly polluted by numerical noise. To confirm this logical deduction, Liao (Reference Liao2009) proposed the ‘clean numerical simulation’ (CNS). Unlike DNS whose algorithms are solved using double precision (i.e. with a fixed level of round-off error), CNS (Liao Reference Liao2013; Liao & Wang Reference Liao and Wang2014; Li & Liao Reference Li and Liao2017; Li et al. Reference Li, Jing and Liao2018; Hu & Liao Reference Hu and Liao2020; Qin & Liao Reference Qin and Liao2020; Liao et al. Reference Liao, Li and Yang2022; Liao & Qin Reference Liao and Qin2022; Liao Reference Liao2023) reduced both the truncation error (by using high-order Taylor expansions in time and a Fourier spectral method in space) and the round-off error (by using multiple precision arithmetic) to such a predetermined minimum level that numerical noise remained negligible (compared with the exact solution) throughout a finite time period that was sufficiently long to accurately compute the flow statistics. Using CNS as a new simulation tool to gain clean turbulence, Qin & Liao (Reference Qin and Liao2022) found that numerical noise, as a stochastic disturbance, could indeed have a significant impact on sustained turbulence, not only in terms of its spatiotemporal trajectory, but also in its statistics and macro-level properties (such as flow type). Notably, for NS turbulence subject to a symmetric initial condition, all CNS results maintain the same spatial symmetry as for the initial conditions in a finite but long enough interval of time. However, all corresponding DNS results (without enforcing the same spatial symmetry for
$t \gt 0$
) maintain the same spatial symmetry as for the initial conditions only over a short initial period after which symmetry is completely lost (Qin et al. Reference Qin, Yang, Huang, Mei, Wang and Liao2024; Liao & Qin Reference Liao and Qin2025; Qin & Liao Reference Qin and Liao2025; Qin et al. Reference Qin, Xu and Liao2026), indicating that numerical noise, which is random in space, indeed increases to a macro-level and then destroys the spatial symmetry of DNS results.
Kida (Reference Kida1985) used numerical simulation to investigate a case of two-dimensional turbulence. According to the covariation of NS equations, Kida assumed that the turbulent flow should maintain the same spatial symmetry as its initial condition. To increase computational efficiency, Kida (Reference Kida1985) used spatial symmetry to reduce the size of the computational domain by a factor of 64 so as to achieve higher Reynolds number flows with limited resources. This approach treated spatial symmetry as an observed utility rather than a proven law (Kida Reference Kida1985). However, in violation of this assumption, almost all current state-of-the-art DNS (without using the property of spatial symmetry to increase computational efficiency) of NS turbulence subject to a symmetric initial condition maintains the same spatial symmetry as the initial condition only over a short period of time after the simulation commences, after which the flow then rapidly loses its spatial symmetry (Qin et al. Reference Qin, Yang, Huang, Mei, Wang and Liao2024; Liao & Qin Reference Liao and Qin2025; Qin & Liao Reference Qin and Liao2025; Qin et al. Reference Qin, Xu and Liao2026). Note that DNS results have been widely used as benchmark solutions by the turbulence community for more than half a century (Orszag Reference Orszag1970; Rogallo Reference Rogallo1981; She et al. Reference She, Jackson and Orszag1990; Moin & Mahesh Reference Moin and Mahesh1998; Coleman & Sandberg Reference Coleman and Sandberg2010; Alexakis et al. Reference Alexakis, Marino, Mininni, van Kan, Foldes and Feraco2024). Note that all current state-of-the-art DNS results (obtained without using spatial symmetry to increase computational efficiency) of NS turbulence subject to an initial condition with spatial symmetry rapidly violate Kida’s assumption of fully symmetric flow (Qin et al. Reference Qin, Yang, Huang, Mei, Wang and Liao2024; Liao & Qin Reference Liao and Qin2025; Qin & Liao Reference Qin and Liao2025; Qin et al. Reference Qin, Xu and Liao2026). Should we trust the DNS results in such cases? In addition, it has been established that CNS results maintain the same spatial symmetry as the initial condition over a prescribed long, finite period (Qin et al. Reference Qin, Yang, Huang, Mei, Wang and Liao2024; Liao & Qin Reference Liao and Qin2025; Qin & Liao Reference Qin and Liao2025; Qin et al. Reference Qin, Xu and Liao2026). It is therefore unknown whether or not the exact solution of NS turbulence subject to a symmetric initial condition should maintain the same spatial symmetry forever. Note that, despite the widespread use of spatial symmetry, a global mathematical proof of symmetry persistence is a missing ‘cornerstone’ of fluid mechanics. This implies the existence of a ‘reliability crisis’ in turbulence simulation, as mentioned earlier.
Physical insights concerning turbulent flows can be obtained by physical experiments involving actual turbulence, numerical simulation and mathematical investigation of a mathematical model of turbulence. Although numerical simulation is routinely used nowadays, it is very beneficial to have access to theoretical results based on purely mathematical proofs of NS turbulence without recourse to numerical simulations. In § 2, a mathematical proof is given of the spatial symmetry invariance theorem (SSIT), whereby the exact solution of the NS equations under periodic boundary conditions maintains the same spatial symmetry as its symmetric initial flow, forever. Numerical confirmation and applications of the SSIT are described in § 3, which illustrates that certain important properties of NS turbulence, including ‘ultra-chaos’ and the ‘noise-expansion cascade’, can be straightforwardly deduced from the SSIT. It is demonstrated that the spatiotemporal trajectory of traditional DNS is quickly badly polluted by numerical noise, as the DNS results rapidly violate the SSIT. It is found that the loss of symmetry in traditional DNS is not a ‘physical instability’, but instead a ‘numerical artefact’ caused by numerical nose and the chaotic nature of NS turbulence. The foregoing enables deeper physical insight into turbulence. Section 4 presents a detailed discussion of the findings. Section 5 draws the main conclusions.
2. Spatial symmetry invariance of NS turbulence
Consider dimensionless three-dimensional (3-D) incompressible NS equations (Batchelor Reference Batchelor2000):
\begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol {u} &= 0,\\[-10pt]\nonumber \end{align}
\begin{align} \frac {\partial \boldsymbol {u}}{\partial t} +\left ( \boldsymbol {u}\boldsymbol{\cdot }\boldsymbol{\nabla }\right ) \boldsymbol {u} &=-\boldsymbol{\nabla }p + \frac {1}{Re} \Delta \boldsymbol {u} + \boldsymbol {f}, \end{align}
subject to periodic boundary conditions in a square domain
$\boldsymbol {x}=(x_1,x_2,x_3)\in [0,2\pi ]^3$
, where
$t\in [0,+\infty )$
denotes time,
$\boldsymbol {u}=(u_1,u_2,u_3)$
is the fluid flow velocity,
$p$
is pressure,
$\boldsymbol {f}(\boldsymbol {x})$
is a steady-state external force,
$Re$
is the Reynolds number based on a characteristic length,
$\boldsymbol{\nabla}$
is a Hamilton operator and
$\varDelta =\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\nabla}$
is a Laplace operator. Let
$O-x_{1}x_{2}x_{3}$
denote its corresponding coordinate system.
Define
\begin{align} \boldsymbol {u}^{(m)}(\boldsymbol {x},t) =\frac {1}{m!} \frac {\partial \boldsymbol {u}(\boldsymbol {x},t)}{\partial t^{m}},\;\;\; p^{(m)}(\boldsymbol {x},t) =\frac {1}{m!} \frac {\partial p(\boldsymbol {x},t)}{\partial t^{m}}, \;\;\; m\geq 0, \end{align}
with
\begin{align} \boldsymbol {u}^{(0)}(\boldsymbol {x},t) = \boldsymbol {u}(\boldsymbol {x},t), \;\;\; p^{(0)}(\boldsymbol {x},t) = p(\boldsymbol {x},t), \end{align}
for integer
$m\geq 0$
.
From (2.2), it holds that
\begin{align} \frac {\partial (\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol {u})}{\partial t} +\boldsymbol{\nabla }\boldsymbol{\cdot }\left [\left ( \boldsymbol {u}\boldsymbol{\cdot }\boldsymbol{\nabla }\right ) \boldsymbol {u}\right ] & = -\Delta p + \frac {1}{Re} \varDelta \left (\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol {u}\right ) + \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol {f}, \end{align}
which, using the continuation equation (2.1), gives that
\begin{align} \Delta p =\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol {f} -\boldsymbol{\nabla }\boldsymbol{\cdot }\left [\left ( \boldsymbol {u}\boldsymbol{\cdot }\boldsymbol{\nabla }\right ) \boldsymbol {u}\right ]. \end{align}
Differentiating (2.2) and (2.6)
$m$
times with respect to
$t$
and then dividing by
$m!$
, we have
\begin{align} (m+1)\boldsymbol {u}^{(m+1)} +\sum _{i=0}^{m}( \boldsymbol {u}^{(i)}\boldsymbol{\cdot }\boldsymbol{\nabla }) \boldsymbol {u}^{(m-i)} =-\boldsymbol{\nabla }p^{(m)} + \frac {1}{Re} \Delta \boldsymbol {u}^{(m)} , \;\; \end{align}
\begin{align} \Delta p^{(m)} =-\sum _{i=0}^{m}\boldsymbol{\nabla }\boldsymbol{\cdot }[( \boldsymbol {u}^{(i)}\boldsymbol{\cdot }\boldsymbol{\nabla }) \boldsymbol {u}^{(m-i)}], \end{align}
for
$m\geq 1$
and
$t\geq 0$
.
Let
$\boldsymbol {u}'(\boldsymbol {x}',t)$
,
$p'(\boldsymbol {x}',t)$
,
$\boldsymbol {f}'(\boldsymbol {x}')$
denote the velocity and pressure of fluid, and the external force in the coordinate system
$O'-x'_{1}x'_{2}x'_{3}$
, where
$\boldsymbol {x}'\in [0,2\pi ]^3$
. Given that they also satisfy (2.2) and (2.6)–(2.8) for
$t\geq 0$
, then
\begin{align} \boldsymbol {u}'^{(1)} +\left ( \boldsymbol {u}'\boldsymbol{\cdot }{\nabla }'\right ) \boldsymbol {u}' & = -{\nabla }' p' + \frac {1}{Re} \Delta ' \boldsymbol {u}' + \boldsymbol {f}',\\[-10pt]\nonumber \end{align}
\begin{align} (m+1)\boldsymbol {u}'^{(m+1)} +\sum _{i=0}^{m} ( \boldsymbol{u}'^{(i)}\boldsymbol{\cdot }{\nabla }') \boldsymbol{u'^{(m-i)}} & = -{\nabla }' p'^{(m)} + \frac {1}{Re} \Delta ' \boldsymbol {u}'^{(m)} , \;\; m\geq 1, \end{align}
and
\begin{align} \Delta ' p' & = {\nabla }'\boldsymbol{\cdot }\boldsymbol {f}' -{\nabla }'\boldsymbol{\cdot }\left [\left ( \boldsymbol {u'}\boldsymbol{\cdot }{\nabla }'\right ) \boldsymbol {u'}\right ],\\[-5pt]\nonumber \end{align}
\begin{align} \Delta ' p'^{(m)} & = -\sum _{i=0}^{m}{\nabla }'\boldsymbol{\cdot }\big[( \boldsymbol {u}'^{(i)}\boldsymbol{\cdot }{\nabla }') \boldsymbol {u}'^{(m-i)}\big], \;\; m\geq 1. \end{align}
In general,
$\boldsymbol {u}'(\boldsymbol {x}',t) \neq \boldsymbol {u}(\boldsymbol {x},t)$
. However, if there exists
$\boldsymbol {u}'(\boldsymbol {x}',t) = \boldsymbol {u}(\boldsymbol {x},t)$
, then there also exists spatial symmetry under the coordinate translation
$x'_{i}=\alpha _i+x_{i}$
and/or the coordinate rotation
$x'_{i}=\beta _i-x_{i}$
with constant
$\alpha _i$
and
$\beta _i$
$(i=1,2,3)$
.
Lemma A (Operator Covariance). Under the coordinate translation
$x'_{i}=x_{i}+\alpha _i$
and/or the coordinate rotation
$x'_{i}=\beta _i-x_{i}$
, where
$\alpha _i$
and
$\beta _i$
$(i=1,2,3)$
are constant, it holds that
\begin{align} {\nabla }' = \boldsymbol{\nabla }, \;\; \varDelta ' = \varDelta . \end{align}
Proof. Let
$\boldsymbol {e}_{i}, \boldsymbol {e}'_{i}$
denote axis unit vectors for
$x_{i}, x'_{i}$
, where
$i=1,2,3$
. For the coordinate translation
$x'_{i}=x_{i}+\alpha _i$
with constant
$\alpha _i$
$(i=1,2,3)$
, it holds that
\begin{align} \boldsymbol {e}_{i} = \boldsymbol {e}'_{i}, \;\;\; \frac {\partial }{\partial x_{i}} = \frac {\partial }{\partial x'_{i}} \frac {\partial x'_{i}}{\partial x_{i}} = \frac {\partial }{\partial x'_{i}}, \;\;\; i=1,2,3, \end{align}
and so
\begin{align} {\nabla }' =\sum _{i=1}^{3} \boldsymbol {e}'_{i} \; \frac {\partial }{\partial x'_{i}} =\sum _{i=1}^{3} \boldsymbol {e}_{i} \; \frac {\partial }{\partial x_{i}} = \boldsymbol{\nabla }. \end{align}
Similarly, for the coordinate rotation
$x'_{i}=\beta _i-x_{i}$
, with constant
$\beta _i$
$(i=1,2,3)$
, it holds that
\begin{align} \boldsymbol {e}_{i} = -\boldsymbol {e}'_{i}, \;\;\; \frac {\partial }{\partial x_{i}} = \frac {\partial }{\partial x'_{i}} \frac {\partial x'_{i}}{\partial x_{i}} = -\frac {\partial }{\partial x'_{i}}, \end{align}
and so
\begin{align} {\nabla }' =\sum _{i=1}^{3} \boldsymbol {e}'_{i} \; \frac {\partial }{\partial x'_{i}} =\sum _{i=1}^{3} \left (-\boldsymbol {e}_{i}\right ) \; \left (- \frac {\partial }{\partial x_{i}}\right ) =\sum _{i=1}^{3} \boldsymbol {e}_{i} \; \frac {\partial }{\partial x_{i}} = \boldsymbol{\nabla }. \end{align}
So, for
$x'_{i}=x_{i}+\alpha _i$
and/or
$x'_{i}=\beta _i-x_{i}$
$(i=1,2,3)$
, we have
\begin{align} \varDelta ' = {\nabla }' \boldsymbol{\cdot }{\nabla }' = \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\nabla }=\varDelta . \end{align}
This ends the proof of Lemma A whereby the Hamilton and Laplacian operators are covariant under specific spatial transformations.
Lemma B (Recursive Symmetry). For the incompressible Navier–Stokes equations (2.1) and (2.2) subject to periodic boundary conditions, if there exist the spatial symmetry
$\boldsymbol {u}'(\boldsymbol {x}',t_{0}) = \boldsymbol {u}(\boldsymbol {x},t_{0})$
at
$t=t_{0}$
and
$\boldsymbol {f}'(\boldsymbol {x}')= \boldsymbol {f}(\boldsymbol {x})$
, where
$x'_i = \alpha _i+x_i$
and/or
$x'_i = \beta _i-x_i$
with constant
$\alpha _i$
and
$\beta _i$
$(i=1,2,3)$
, and in addition there exists a temporal Taylor series of
$\boldsymbol { u}(\boldsymbol {x},t)$
at
$t=t_{0}$
, then there exist spatial symmetry
$\boldsymbol {u}'^{(m+1)}(\boldsymbol {x}',t_{0}) = \boldsymbol {u}^{(m+1)}(\boldsymbol {x},t_{0})$
and
$p'^{(m)}(\boldsymbol {x}',t_{0}) = p^{(m)}(\boldsymbol {x},t_{0})$
at
$t=t_{0}$
for all integer
$m\geq 0$
.
Proof. Lemma B can be proved by means of the recursive method.
(A) When
$m=0$
, (2.11) at
$t=t_{0}$
becomes
\begin{align} \Delta p'=\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol {f} -\boldsymbol{\nabla }\boldsymbol{\cdot }\left [\left ( \boldsymbol {u}\boldsymbol{\cdot }\boldsymbol{\nabla }\right ) \boldsymbol {u}\right ], \end{align}
given that
$\varDelta ' =\varDelta , {\nabla }' = \boldsymbol{\nabla}$
according to Lemma A, and, in addition, there exist the spatial symmetry
$\boldsymbol {u}'(\boldsymbol {x}',t_0) = \boldsymbol {u}(\boldsymbol {x},t_{0})$
and
$\boldsymbol {f}'(\boldsymbol {x}') = \boldsymbol {f}(\boldsymbol {x})$
. Comparing (2.19) with (2.6), we have
$\Delta p =\Delta p'$
such that
$p'(\boldsymbol {x}',t_{0}) = p(\boldsymbol {x},t_{0})$
, i.e.
$p'^{(0)}(\boldsymbol {x}',t_{0}) = p^{(0)}(\boldsymbol {x},t_{0})$
.
Similarly, (2.9) at
$t=t_{0}$
becomes
\begin{align} \boldsymbol {u}'^{(1)} +\left ( \boldsymbol {u}\boldsymbol{\cdot }\boldsymbol{\nabla }\right ) \boldsymbol {u} =-\boldsymbol{\nabla }p + \frac {1}{Re} \Delta \boldsymbol {u} + \boldsymbol {f}, \end{align}
because
$\varDelta ' =\varDelta ,{\nabla }' = \boldsymbol{\nabla}$
according to Lemma A and there exist the spatial symmetries
$\boldsymbol {u}'(\boldsymbol {x}',t_0) = \boldsymbol {u}(\boldsymbol { x},t_{0})$
,
$\boldsymbol {f}'(\boldsymbol {x}') = \boldsymbol {f}(\boldsymbol {x})$
and
$p'(\boldsymbol {x}',t_{0}) = p(\boldsymbol {x},t_{0})$
as mentioned earlier. Comparing (2.20) with (2.2) gives
$\boldsymbol {u}'^{(1)}(\boldsymbol {x}',t_0) = \boldsymbol {u}^{(1)}(\boldsymbol {x},t_{0})$
.
Thus, Lemma B holds when
$m=0$
.
(B) Assume that Lemma B holds when
$m = n \geq 0$
, and there exists the spatial symmetry
\begin{align} \boldsymbol {u}'^{(i+1)}(\boldsymbol {x}',t_{0}) = \boldsymbol {u}^{(i+1)}(\boldsymbol {x},t_{0}), \;\; p'^{(i)}(\boldsymbol {x}',t_{0}) = p^{(i)}(\boldsymbol {x},t_{0}), \;\; 0\leq i\leq n. \end{align}
Then, (2.12 ) at
$t=t_{0}$
becomes
\begin{align} \Delta p'^{(n+1)} =-\sum _{i=0}^{n+1}\boldsymbol{\nabla }\boldsymbol{\cdot }\big[( \boldsymbol {u}^{(i)}\boldsymbol{\cdot }\boldsymbol{\nabla }) \boldsymbol {u}^{(n+1-i)}\big], \end{align}
because
$\varDelta ' =\varDelta , {\nabla }' = \boldsymbol{\nabla}$
according to Lemma A and there exist the spatial symmetries given by (2.21). Comparing (2.22) with (2.8) gives
$p'^{(n+1)}(\boldsymbol {x}',t_{0}) = p^{(n+1)}(\boldsymbol {x},t_{0})$
.
Then, according to the spatial symmetry (2.21), (2.10) (when
$m=n+1$
) at
$t=t_{0}$
becomes
\begin{align} (n+2)\boldsymbol {u}'^{(n+2)} +\sum _{i=0}^{n+1}( \boldsymbol {u}^{(i)}\boldsymbol{\cdot }\boldsymbol{\nabla }) \boldsymbol {u}^{(n+1-i)} = -\boldsymbol{\nabla }p^{(n+1)} + \frac {1}{Re} \Delta \boldsymbol {u}^{(n+1)}, \end{align}
because
$\varDelta ' =\varDelta , {\nabla }' = \boldsymbol{\nabla}$
according to Lemma A and
$p'^{(n+1)}(\boldsymbol {x}',t_{0}) = p^{(n+1)}(\boldsymbol {x},t_{0})$
as proved earlier. Setting
$m=n+1$
in (2.7) at
$t=t_{0}$
gives
\begin{align} (n+2)\boldsymbol {u}^{(n+2)} +\sum _{i=0}^{n+1}( \boldsymbol {u}^{(i)}\boldsymbol{\cdot }\boldsymbol{\nabla }) \boldsymbol {u}^{(n+1-i)} =-\boldsymbol{\nabla }p^{(n+1)} + \frac {1}{Re} \Delta \boldsymbol {u}^{(n+1)}. \;\; \end{align}
Comparing (2.23) with (2.24) gives
$\boldsymbol {u}'^{(n+2)}(\boldsymbol {x}',t_{0}) =\boldsymbol {u}^{(n+2)}(\boldsymbol {x},t_{0})$
.
Thus, if Lemma B holds for
$m=n\geq 0$
, it holds for
$m=n+1$
.
(C) According to parts (A) and (B), Lemma B holds for arbitrary integer
$m\geq 0$
.
This ends the proof of Lemma B. It has used mathematical induction to show that if the initial condition is symmetric, all higher-order temporal derivatives must inherit exactly the same symmetry.
Lemma C (Symmetry Continuation). For the incompressible continuity and Navier–Stokes equations (2.1) and (2.2) subject to periodic boundary conditions, if there exist the spatial symmetries
$\boldsymbol {u}'(\boldsymbol {x}',t_{0}) = \boldsymbol {u}(\boldsymbol {x},t_{0})$
at
$t=t_{0}$
and
$\boldsymbol {f}'(\boldsymbol {x}')= \boldsymbol {f}(\boldsymbol {x})$
, where
$x'_i = \alpha _i+x_i$
and/or
$x'_i = \beta _i-x_i$
with constant
$\alpha _i$
and
$\beta _i$
$(i=1,2,3)$
, and in addition there exists a temporal Taylor series for
$\boldsymbol {u}(\boldsymbol {x},t)$
at
$t=t_{0}$
, with a non-zero radius
$\rho$
of convergence, then there exists the spatial symmetry
$\boldsymbol {u}'(\boldsymbol {x}',t) = \boldsymbol {u}(\boldsymbol {x},t)$
in
$t\in [t_{0},\rho )$
.
Proof. According to Lemma B, the temporal Taylor series of
$\boldsymbol {u}(\boldsymbol {x},t)$
at
$t=t_{0}$
reads
\begin{align} \boldsymbol {u}(\boldsymbol {x},t) &=& \sum _{m=0}^{+\infty } \boldsymbol {u}^{(m)}(\boldsymbol {x},t_{0}) (t-t_{0})^{m}= \sum _{m=0}^{+\infty } \boldsymbol {u}'^{(m)}(\boldsymbol {x}',t_{0}) (t-t_{0})^{m}=\boldsymbol {u}'(\boldsymbol {x}',t) \end{align}
in
$|t-t_{0}|\lt \rho _{0}$
, which includes
$t\in [t_{0},t_{0}+\rho _{0})$
. This ends the proof of Lemma C.
Finally, let us give the mathematical proof of the following.
Spatial Symmetry Invariance Theorem (SSIT). The exact solution
$\boldsymbol {u}$
of incompressible three-dimensional continuity and Navier–Stokes equations (2.1) and (2.2) subject to periodic boundary conditions have the same spatial symmetry
$\boldsymbol {u}'(\boldsymbol {x}',t)$
$ = \boldsymbol {u}(\boldsymbol {x},t)$
for all
$t\gt 0$
as its smooth initial condition, where
$x'_i = \alpha _{i}+x_i$
and/or
$x'_i = \beta _{i}-x_i$
with constants
$\alpha _{i}$
and
$\beta _{i}$
$(i=1,2,3)$
,provided there exist spatial symmetries
$\boldsymbol {u}'(\boldsymbol {x}',0)= \boldsymbol {u}(\boldsymbol {x},0)$
for initial condition and
$\boldsymbol {f}'(\boldsymbol {x}')= \boldsymbol {f}(\boldsymbol {x})$
for the external force, and in addition a temporal Taylor series for
$\boldsymbol {u}(\boldsymbol {x},t)$
exists with non-zero radius of convergence for
$t\geq 0$
.
Proof. Since there exist the spatial symmetries
$\boldsymbol {u}'(\boldsymbol {x}',0) = \boldsymbol {u}(\boldsymbol {x},0)$
and
$\boldsymbol {f}'(\boldsymbol {x}')= \boldsymbol { f}(\boldsymbol {x})$
, where
$x'_i = \alpha _i + x_i$
and/or
$x'_i = \beta _i - x_i$
with constant
$\alpha _i$
and
$\beta _i$
$(i=1,2,3$
), we have according to Lemma C (setting
$t_{0}=0$
) the spatial symmetry
\begin{align} \boldsymbol {u}'(\boldsymbol {x}',t) = \boldsymbol {u}(\boldsymbol {x},t), \hspace {1.0cm} t\in [0,\delta _{0}], \end{align}
where
$\delta _{0} \lt \rho _{0}$
, and
$\rho _{0}\gt 0$
is the non-zero radius of convergence of its Taylor series at
$t=0$
.
Similarly, since there exists the spatial symmetries
$\boldsymbol {u}'(\boldsymbol {x}',\delta _{0}) = \boldsymbol {u}(\boldsymbol {x},\delta _{0})$
at
$t=\delta _{0}$
and
$\boldsymbol { f}'(\boldsymbol {x}')= \boldsymbol {f}(\boldsymbol {x})$
, we have according to Lemma C (setting
$t_{0}=\delta _{0}$
) the spatial symmetry
\begin{align} \boldsymbol {u}'(\boldsymbol {x}',t) = \boldsymbol {u}(\boldsymbol {x},t), \hspace {1.0cm} t\in [\delta _{0},\delta _{1}], \end{align}
where
$\delta _{1} \lt \rho _{1}$
, and
$\rho _{1}\gt 0$
is the non-zero radius of convergence of its Taylor series at
$t=\delta _{0}$
.
Note that
\begin{align} [0,\delta _{0}] \cup [\delta _{0},\delta _{1}] = [0,\delta _{1}]. \end{align}
Hence, there exists the spatial symmetry
\begin{align} \boldsymbol {u}'(\boldsymbol {x}',t) = \boldsymbol {u}(\boldsymbol {x},t), \hspace {1.0cm} t\in [0,\delta _{1}], \end{align}
over a larger interval of time. The same process can continue and the time interval becomes larger and larger because the temporal Taylor series of
$\boldsymbol { u}(\boldsymbol {x},t)$
invariably has a non-zero radius of convergence. Thus, it holds the spatial symmetry
\begin{align} \boldsymbol {u}'(\boldsymbol {x}',t) = \boldsymbol {u}(\boldsymbol {x},t), \hspace {1.0cm} t\in [0,+\infty ). \end{align}
This ends the proof of the spatial symmetry invariance theorem (SSIT).
Note that, for the two-dimensional Kolmogorov flow in the form of stream function, Liao (Reference Liao2026a , Reference Liaob ) proved the similar theorem. The mathematical proof described earlier for the three-dimensional deterministic NS equations indicates that the SSIT has the general meanings.
3. Numerical confirmation and application of SSIT
It should be emphasised that the spatial symmetry invariance theorem (SSIT) provides an ‘absolute scale’ for computational fluid dynamics (CFD). In other words, if a DNS fails to maintain the same symmetry of a symmetric initial state, the numerical simulation is mathematically invalidated as a true solution of NS equations.
Assume that exact solution
${\mathcal{F}}(\boldsymbol {x},t)$
of the continuity and NS equations has a spatial symmetry
${\mathcal{F}}(\boldsymbol {x},t)={\mathcal{F}}(\boldsymbol {x}',t)$
in
$\boldsymbol {x}\in [a,b]^{3}$
for
$t\in [0,+\infty ]$
. In any numerical simulation, there should invariably exist some difference between
${\mathcal{F}}(\boldsymbol {x},t)$
and
${\mathcal{F}}(\boldsymbol {x}',t)$
. So, we define spatial symmetry loss (denoted by
$\delta$
) of a numerical simulation by
\begin{align} \delta (t) =\left \{ \frac {1}{(b-a)^{3}}\int _{a}^{b}\int _{a}^{b}\int _{a}^{b} \left [ {\mathcal{F}}(\boldsymbol {x},t) - {\mathcal{F}}(\boldsymbol {x}',t) \right ]^{2}\,{\rm d}x_{1}\, {\rm d} x_{2}\, {\rm d} x_{3}\right \}^{1/2}. \end{align}
Obviously, spatial symmetry is maintained when
$\delta \ll 1$
, but is destroyed when
$\delta \sim O(1)$
.
3.1. Two-dimensional Kolmogorov flow
Qin et al. (Reference Qin, Yang, Huang, Mei, Wang and Liao2024) applied CNS to solve a two-dimensional (2-D) Kolmogorov flow in the form of a stream-function
$\psi$
:
\begin{align} \frac {\partial }{\partial t} ( {\nabla} ^{2}\psi )+\frac {\partial (\psi ,{\nabla} ^{2}\psi )}{\partial (x,y)}-\frac {1}{Re}{\nabla} ^{4}\psi +n_K\cos (n_Ky)=0, \end{align}
where
$x,y\in [0,2\pi ]$
are horizontal and vertical coordinates,
$t$
denotes time,
$Re$
the Reynolds number,
${\nabla} ^{2}$
a Laplace operator,
${\nabla} ^{4}={\nabla} ^{2}{\nabla} ^{2}$
, and
\begin{align} \frac {\partial (a,b)}{\partial (x,y)}=\frac {\partial a}{\partial x}\frac {\partial b}{\partial y}-\frac {\partial b}{\partial x}\frac {\partial a}{\partial y} \end{align}
is the Jacobi operator, subject to the periodic boundary condition
\begin{align} \psi (x,y,t) = \psi (x\pm 2\pi ,y\pm 2\pi ,t). \end{align}
Following Chandler & Kerswell (Reference Chandler and Kerswell2013), Qin et al. (Reference Qin, Yang, Huang, Mei, Wang and Liao2024) prescribed the physical parameters
$n_K$
= 4,
$Re$
= 40 and the symmetric initial condition
\begin{align} & \psi (x,y,0)=-\frac {1}{2}[ \cos (x+y)+\cos (x-y)], \end{align}
corresponding to a sustained turbulent Kolmogorov flow.
The initial condition (3.5) has the spatial symmetries in translation (
$s_{1}$
) and rotation (
$s_{2}$
) given by
\begin{align} \left \{\!\! \begin{array} {lll} \mbox{translation} &:& s_{1} \circ \psi (x,y,t) \mapsto \psi (x+\pi ,y+\pi ,t),\\[3pt] \mbox{rotation} &:& s_{2} \circ \psi (x,y,t) \mapsto \psi (-x,-y,t), \end{array} \right . \end{align}
respectively. For
$n_k=4$
, (3.2) is covariant under the above-mentioned spatial symmetries in translation and rotation, so that, according to the spatial symmetry invariance theorem (SSIT), its exact solution must maintain the same spatial symmetries (3.6) forever, i.e. for
$t\in [0,+\infty )$
.
Qin et al. (Reference Qin, Yang, Huang, Mei, Wang and Liao2024) used both of CNS and DNS to solve this 2-D Kolmogorov flow. Unlike Kida (Reference Kida1985), they did not use the spatial symmetry (3.6) to increase computational efficiency of DNS and CNS. In other words, their DNS and CNS algorithms are valid for a general initial condition, no matter whether it is symmetric or not. As shown in figure 1, the spatial symmetry losses of the CNS result are much less than 1 over the whole simulation duration
$t\in [0,1500]$
, indicating that the CNS results maintain spatial symmetry throughout the whole duration of simulation. However, as shown in figure 1, the spatial symmetry losses of the DNS result increase to a macro-level
$O(1)$
at
$t\approx 110$
, indicating that the DNS result maintains the spatial symmetries only in a short duration
$t\in [0,110]$
, after which time, the results rapidly lose spatial symmetry. Hence, the DNS result for
$t\gt 110$
violates the SSIT. In other words, based on the SSIT, we are quite certain that the DNS result after
$t\gt 110$
departs far away from the exact solution for turbulent Kolmogorov flow. This test also clearly indicates that loss of symmetry in traditional DNS is not due to a ‘physical instability’, but is in fact a ‘numerical artefact’ caused by the chaotic nature of NS turbulence. Moreover, Qin et al. (Reference Qin, Yang, Huang, Mei, Wang and Liao2024) found that the statistics given by the CNS are distinctly different from those given by the DNS. This is quite reasonable from the viewpoint of SSIT because the CNS results maintain spatial symmetry (3.6), but the DNS results rapidly lose symmetry due to stochastic numerical noise (that has no spatial symmetry at all) and the chaotic property of NS turbulence.
Time histories of spatial symmetry loss given by CNS (red) and DNS (blue) for 2-D Kolmogorov flow (Qin et al. Reference Qin, Yang, Huang, Mei, Wang and Liao2024). (a) Spatial symmetry of translation (
$s_{1}$
); (b) spatial symmetry of rotation (
$s_{2}$
). Spatial symmetry is defined by (3.6).
3.2. Three-dimensional Kolmogorov flow
Qin & Liao (Reference Qin and Liao2025) considered 3-D Kolmogorov flow (Wu et al. Reference Wu, Schmitt, Calzavarini and Wang2021), governed by the 3-D incompressible Navier–Stokes equation
\begin{align} \frac {\partial \boldsymbol {u}}{\partial t}+(\boldsymbol {u}\boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol {u} = -\boldsymbol{\nabla }p + \frac {1}{Re}\Delta \boldsymbol {u} + \boldsymbol {f}, \end{align}
and the continuity equation
$\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol {u}=0$
, where
$t$
denotes time,
$\boldsymbol {x}=(x, y, z)\in [0,1]^3$
is the Cartesian coordinate,
$\boldsymbol {u}= u\; \boldsymbol {e}_x +v\; \boldsymbol {e}_y + w\; \boldsymbol {e}_z$
is the velocity field,
$Re=1211.5$
is the Reynolds number, and the external force
$\boldsymbol {f}$
is of the form
$\boldsymbol {f} = \sin (8\pi z)\,\boldsymbol {e}_x$
, subject to a periodic boundary condition and the symmetric initial condition
\begin{align} \left \{\!\!\begin{array}{l} u(x,y,z,0)=\cos (2\pi x)\cos (2\pi y)\sin (2\pi z)/2, \\[3pt] v(x,y,z,0)=\sin (2\pi x)\sin (2\pi y) \sin (2\pi z), \\[3pt] w(x,y,z,0)=\sin (2\pi x)\cos (2\pi y)\cos (2\pi z)/2, \\ \end{array} \right . \end{align}
which has the spatial symmetry of rotation (
$s_{3}$
)
\begin{align} s_{3}\circ \boldsymbol {u}(x,y,z,t)\mapsto -\boldsymbol {u}(1-x,1-y,1-z,t). \end{align}
Note that the external force
$\boldsymbol {f} = \sin (8\pi z)\,\boldsymbol {e}_x$
has the same spatial symmetry as (3.9), too. Thus, according to the SSIT, the exact solution must maintain the same spatial symmetry (3.9) forever, i.e. for
$t\in [0,+\infty )$
.
Time histories of spatial symmetry loss of velocity given by CNS (red) and DNS (blue) for the 3-D Kolmogorov flow (Qin & Liao Reference Qin and Liao2025). (a)
$u(x,y,z,t)$
; (b)
$v(x,y,z,t)$
; (c)
$w(x,y,z,t)$
. The spatial symmetry (
$s_{3}$
) is defined by (3.9).
Qin & Liao (Reference Qin and Liao2025) applied both of CNS and DNS to solve the above-mentioned case of 3-D Kolmogorov flow. As shown in figure 2, the spatial symmetry loss of the CNS result is much less than 1 in the whole simulation of duration
$t\in [0,500]$
, indicating that the CNS result maintains the same spatial symmetry (3.9) as the initial condition (3.8) throughout the whole duration of simulation. However, as shown in figure 2, the spatial symmetry loss of the DNS result quickly increases to macro-level
$O(1)$
at
$t\approx 100$
, indicating that the corresponding DNS result maintains the spatial symmetry only over a brief initial period and then completely loses the symmetry when
$t\gt 100$
. In other words, the CNS result complies with the SSIT throughout the whole simulation, but the DNS result violates the SSIT once
$t \gt 100$
. Therefore, based on the SSIT, we find that the spatiotemporal trajectories given by DNS are quickly badly polluted by numerical noise, which is random in space and destroys the spatial symmetry when reaching the macro-level. This again clearly illustrates that loss of symmetry in traditional DNS is not created by a ‘physical instability’, but is rather a ‘numerical artefact’ caused by numerical noise and the chaotic nature of NS turbulence.
3.3. Revealing mathematical reality of NS turbulence
3.3.1. Ultra-chaotic property
Let us consider such a set of different initial flows of NS turbulence
\begin{align} \boldsymbol {u}_{i} = \boldsymbol {u}_{0} + \epsilon \; \boldsymbol {u}'_{i}, \hspace {1.0cm} i=1,2,\ldots , \end{align}
where
$\boldsymbol {u}_{i}$
has a different spatial symmetry from each other when the constant
$\epsilon$
is at a macroscopic level, i.e.
$\epsilon \sim O(1)$
. According to the SSIT, the exact solution of the NS turbulence given by the initial flow
$\boldsymbol {u}_{i}$
should have different spatial symmetries that correspond to different statistics. Then, let us consider the case when
$\epsilon$
is very small, for example,
$\epsilon =10^{-10}$
so that the above-mentioned initial flows are almost the same and
$\epsilon \; \boldsymbol {u}'_{i}$
can be regarded as a small disturbance of the initial flow. However, since NS turbulence has the chaotic property, each small disturbance
$\epsilon \; \boldsymbol {u}'_{i}$
certainly increases to macro-level, and then, according to the SSIT, their exact solutions of the corresponding NS turbulence should have different spatial symmetries (and also different statistics). Recalling the fact that these initial flows are almost the same when
$\epsilon =10^{-10}$
, their statistics are certainly rather sensitive to the initial flows. In other words, NS turbulence should have ultra-chaotic property, a novel concept proposed by Liao & Qin (Reference Liao and Qin2022).
The above-mentioned logical deduction has been verified by CNS. Qin et al. (Reference Qin, Xu and Liao2026) applied both of CNS and DNS to solve the 2-D Kolmogorov equation (3.2) subject to the periodic boundary condition (3.4), but using the physical parameters
$n_K=16$
,
$Re=2000$
and the following three different initial conditions:
\begin{align} \varPhi _{1,0}(x,y)=\psi (x,y,0) & = -\frac {1}{2}\big [\cos (4x+4y) + \cos (4x-4y) + \sin (4x+4y)\big ] \nonumber \\[4pt]& \quad + 10^{-10}\big [\sin (x+y) + \sin (x-y)\big ], \end{align}
\begin{align} \varPhi _{2,0}(x,y)= \psi (x,y,0) & = -\frac {1}{2}\big[\cos (4x+4y) + \cos (4x-4y) + \sin (4x+4y)\big] \nonumber \\[4pt] & \quad + 10^{-10}\big[\sin (2x+y) + \sin (2x-y)\big ], \end{align}
\begin{align} \varPhi _{3,0}(x,y)= \psi (x,y,0) & = -\frac {1}{2}\big[\cos (4x+4y) + \cos (4x-4y) + \sin (4x+4y)\big] \nonumber \\[4pt]& \quad + 10^{-10}\big[\sin (4x+y) + \sin (4x-y)\big], \end{align}
corresponding to three turbulent Kolmogorov flows with different spatial symmetries (
$s_{4}, s_{5}, s_{6}$
) of translation
\begin{align} \varPhi _{1,0}(x,y) &: s_{4}\circ \psi (x,y,t) \mapsto \psi (x+\pi ,y+\pi ,t),\\[-10pt]\nonumber \end{align}
\begin{align} \varPhi _{2,0}(x,y) &: s_{5}\circ \psi (x,y,t) \mapsto \psi (x+\pi /2,y+\pi ,t),\\[-10pt]\nonumber \end{align}
\begin{align} \varPhi _{3,0}(x,y) &: s_{6}\circ \psi (x, y,t) \mapsto \psi (x+\pi /2, y,t). \end{align}
It should be emphasised that the above-mentioned three spatial symmetries are quite different, although the deviations between the three initial conditions (3.11), (3.12) and (3.13) are extremely small, say,
\begin{align} \left | \varPhi _{i,0}(x,y) - \varPhi _{j,0}(x,y) \right | \leq 4 \times 10^{-10}, \hspace {0.5cm} 1\leq i,j\leq 3, i\neq j, \end{align}
so that they are almost the same.
Time histories of spatial symmetry loss given by CNS (red) and DNS (blue) for the 2-D Kolmogorov flow (Qin et al. Reference Qin, Xu and Liao2026). (a) Comparison between Flow CNS-1 (red) and Flow DNS-1 (blue) for the spatial symmetry (3.14); (b) comparison between Flow CNS-2 (red) and Flow DNS-2 (blue) for the spatial symmetry (3.15); (c) comparison between Flow CNS-3 (red) and Flow DNS-3 (blue) for the spatial symmetry (3.16).
As reported by Qin et al. (Reference Qin, Xu and Liao2026), the CNS results maintain the same spatial symmetries (3.14), (3.15) and (3.16) as the corresponding initial conditions (3.11), (3.12) and (3.13), respectively, and thus comply with the SSIT throughout the whole simulation duration
$t\in [0,300]$
. Figure 3 confirms that this is indeed true in that the spatial symmetry loss
$\delta (t)$
of each set of CNS results is much less than 1, indicating that Flow CNS-1 corresponding to initial condition (3.11), Flow CNS-2 corresponding to initial condition (3.12) and Flow CNS-3 corresponding to initial condition (3.13) always maintain their corresponding spatial symmetries (3.14), (3.15) and (3.16). However, as shown in figure 3, the spatial symmetry loss of each set of DNS results quickly escalates to the macro-level at
$t \approx 100$
, clearly indicating that the DNS results only maintain the spatial symmetry over a brief interval of time
$t\in [0,100]$
, but thereafter quickly lose spatial symmetry so that the DNS results for
$t\gt 100$
violate the SSIT! In other words, according to the SSIT, it is absolutely certain that the DNS results for
$t\gt 100$
depart far away from the exact solutions of the turbulent Kolmogorov flow.
Qin et al. (Reference Qin, Xu and Liao2026) found that, although the deviations between the three initial conditions (3.11), (3.12) and (3.13) are very small, the statistics of the corresponding CNS results are entirely different. In other words, the statistical results of the NS turbulence exhibit sensitivity dependence on the initial conditions. Thus, the corresponding Kolmogorov flow is ultra-chaotic (Liao & Qin Reference Liao and Qin2022). Note that ultra-chaos, whose statistics are sensitive to small disturbances, has a higher level of disorder than normal-chaos, whose statistics are not sensitive to small disturbances. A typical example of normal-chaos is given by the famous Lorenz equation, for which it has been established that although a chaotic numerical trajectory given by a traditional algorithm (using double precision arithmetic) deviates significantly at
$t\approx 35$
from the corresponding CNS result (Liao & Wang Reference Liao and Wang2014) that is convergent and reliable in
$t\in [0,10000]$
; its statistics agree well with those given by the CNS result. Here, the concept ‘ultra-chaos’ (Liao & Qin Reference Liao and Qin2022) complies with the spirit of SSIT because the three initial conditions (3.11), (3.12) and (3.13), although almost the same, have different spatial symmetries which should always persist when the small initial differences evolve to the macro-level due to the chaotic property (i.e. butterfly-effect) of NS turbulence. Thus, the statistics of the three sets of CNS results certainly should be obviously different from each other. This reveals the following non-uniqueness conjecture for NS turbulence: for almost the same initial conditions, distinct globally smooth solutions of NS turbulent flows may be obtained with different spatial symmetries and statistics. In other words, the ultra-chaotic characteristic (Qin et al. Reference Qin, Xu and Liao2026) of NS turbulence is quite reasonable and straight-forward regarding both of the spatial symmetry invariance theorem (SSIT) and the chaotic property of NS turbulence!
3.3.2. Noise-expansion cascade
Let us consider such an initial flow of NS turbulence,
\begin{align} \boldsymbol {u}_{ 0}= \boldsymbol {u}'_{0} + \sum _{i=1}^{2}\epsilon _{i} \; \boldsymbol {u}'_{i}, \hspace {1.0cm} \end{align}
where
$\boldsymbol {u}'_{0}$
,
$\boldsymbol {u}'_{0}+\epsilon _{1} \boldsymbol {u}'_{1}$
and
$\boldsymbol {u}'_{0}+\epsilon _{1} \boldsymbol {u}'_{1}+\epsilon _{2} \boldsymbol {u}'_{2}$
have different spatial symmetries, namely
$\alpha , \beta$
and
$\gamma$
, respectively, when the constants
$\epsilon _{1}$
and
$\epsilon _{2}$
are at macroscopic level, i.e.
$\epsilon _{1} \sim O(1)$
and
$\epsilon _{2} \sim O(1)$
. According to the SSIT, the exact solutions of NS turbulence subject to each of these three symmetric initial flows should have different spatial symmetries (and different statistics). Now, let us consider the initial flow (3.18) in the case of
$\epsilon _{2}\ll \epsilon _{1}\ll 1$
, for example,
$\epsilon _{1}=10^{-20}$
and
$\epsilon _{2}=10^{-40}$
. According to the chaotic property of NS turbulence, the first small disturbance
$\epsilon _{1} \boldsymbol {u}'_{1}$
increases to macro-level after a period of time
$T_{1}$
, and then the second small disturbance
$\epsilon _{2} \boldsymbol {u}'_{2}$
increases to macro-level after a period of time
$T_{2}$
. Thus, according to the SSIT, the exact solution of the NS turbulence subject to (3.18) should maintain the spatial symmetry
$\alpha$
from the beginning until
$t \approx T_{1}$
when the first small disturbance
$\epsilon _{1} \boldsymbol {u}'_{1}$
increases to macro-level to replace the spatial symmetry
$\alpha$
by
$\beta$
, and the symmetry
$\beta$
should maintain until
$t\approx T_{1} + T_{2}$
when the second small disturbance
$\epsilon _{2} \boldsymbol {u}'_{2}$
increases to macro-level to replace the spatial symmetry
$\beta$
by
$\gamma$
. In other words, NS turbulence should have the so-called ‘noise-expansion cascade’, which can be directly deduced in logic from the SSIT.
The so-called ‘noise-expansion cascade’ has been first discovered by Liao & Qin (Reference Liao and Qin2025) who applied CNS to solve the 2-D Kolmogorov equation (3.2) subject to the periodic boundary condition (3.4) for
$n_K=16$
and
$Re=2000$
, but subject to the initial condition
\begin{align} \psi (x,y,0) = -\frac {1}{2}\left [ \cos (x+y)+\cos (x-y) \right ] + \delta _{1} \sin (x+y) +\delta _{2} \sin (2x + y), \end{align}
where
$\delta _{1}$
and
$\delta _{2}$
are constants. This involves the two kinds of spatial symmetries
\begin{align} \left \{\!\! \begin{array}{lll} \mbox{translation} &:& s_{1}\circ \psi (x,y,t) \mapsto \psi (x+\pi ,y+\pi ,t), \\[3pt] \mbox{rotation} &:& s_{2}\circ \psi (x,y,t) \mapsto \psi (-x,-y,t), \end{array} \right . \hspace {0.3cm} \end{align}
when
$\delta _{1}=0, \delta _{2}=0$
, and only the spatial symmetry of rotation
\begin{align} s_{2}\circ \psi (x,y,t) \mapsto \psi (-x,-y,t), \hspace {0.5cm} \end{align}
when
$\delta _{1}\neq 0, \delta _{2}=0$
, but no spatial symmetry when
$\delta _{1}\neq 0, \delta _{2} \neq 0$
. In the special case of
$\delta _{1}=10^{-20}$
and
$\delta _{2}=10^{-40}$
for the initial condition (3.19), Liao & Qin (Reference Liao and Qin2025) found that the CNS results maintain both kinds of spatial symmetry (3.20) until
$t\approx 35$
when the first small disturbance
$\delta _{1}\sin (x+y)$
of the initial flow increases to a macro-level and destroys the rotation symmetry (3.20), thereafter maintains the translational symmetry (3.21) until
$t\approx 93$
when the second small disturbance
$\delta _{2}\sin (2 x+y)$
of the initial flow increases to macro-level and destroys the spatial symmetry completely. This is indeed true: as shown in figure 4, the rotation spatial symmetry loss increases to macro-level
$O(1)$
at
$t\approx 35$
with the translation spatial symmetry loss remaining well below. At
$t\approx 93$
, the translation spatial symmetry loss reaches macro-level
$O(1)$
and the flow loses both spatial symmetries thereafter. This behaviour is logically reasonable and readily understand from the viewpoint of SSIT.
The results reveal the following mathematical truth underpinning NS turbulence called the ‘noise-expansion cascade’ (Liao & Qin Reference Liao and Qin2025): each micro-level noise and/or disturbance of a particular order of magnitude in the initial conditions of NS turbulence grows progressively, following a kind of inverse cascade until reaching the macro level. Note that, when
$\delta _{2} \ll \delta _{1}$
, the second disturbance to the initial condition (3.19) could be neglected from the traditional DNS viewpoint. However, this is incorrect, because the DNS results violate the SSIT, as illustrated by Liao & Qin (Reference Liao and Qin2025). In other words, the SSIT inevitably leads to the ‘noise-expansion cascade’, first discovered by Liao & Qin (Reference Liao and Qin2025) using CNS.
Time histories of spatial symmetry loss of CNS result for the 2-D Kolmogorov flow (Liao & Qin Reference Liao and Qin2025) under the initial condition (3.19) when
$\delta _{1}=10^{-20}$
and
$\delta _{2}=10^{-40}$
. Blue line, spatial symmetry of translation
$s_{1}\circ \psi (x,y,t) \mapsto \psi (x+\pi ,y+\pi ,t)$
; Red line, spatial symmetry of rotation
$s_{2}\circ \psi (x,y,t) \mapsto \psi (-x,-y,t)$
.
3.3.3. Permanent memory of symmetric initial flow condition
According to the SSIT, the exact solution of NS turbulence should remember the spatial symmetry of initial flow, forever. In other words, the memory about spatial symmetry of initial flow will never been forgotten by the exact solution of NS turbulence. Thus, if initial flow has spatial symmetry, the exact solution of the corresponding NS turbulence should have no mixing. This is a very surprising conclusion, which however is a logical deduction from the SSIT.
Normalised temporal autocorrelation coefficient
$\rho _p(\tau )$
of the pressure field
$p(x,y,z,\tau )$
with respect to the initial pressure field
$p(x,y,z,0)$
of the 3-D Kolmogorov flow (Qin & Liao Reference Qin and Liao2025), given by the CNS benchmark solution (red line) and the DNS result (blue line), where the black dashed line denotes
$\rho =0$
.
The above-mentioned logical deduction was confirmed by Qin & Liao (Reference Liao and Qin2025) who used both of CNS and DNS to solve a 3-D Kolmogorov flow (Wu et al. Reference Wu, Schmitt, Calzavarini and Wang2021), governed by the 3-D incompressible Navier–Stokes equation. As proved in § 2, the exact solution of its pressure field
$p(x,y,z,t)$
should remain the spatial symmetry in the whole interval of time
$t\in [0,+\infty )$
. Consider the temporal autocorrelation coefficient
$\rho (\tau )$
of the pressure field
$p(x,y,z,\tau )$
with respect to the initial pressure
$p(x, y, z, 0)$
normalised by
$p(x, y, z,0)$
, i.e.
\begin{align} \rho (\tau ) = \frac {\big \lt p(x, y, z,0) \; p(x, y, z,\tau ) \big \gt }{\big \lt p^{2}(x, y, z,0) \big \gt }, \end{align}
where
$\langle\rangle$
is an operator of spatial average. Note that
$\rho (\tau )$
might be larger than 1 since it is normalised by
$p(x, y, z,0)$
. As shown in figure 5, the autocorrelation coefficient
$\rho (\tau )$
of the CNS solution of the 3-D Kolmogorov flow (Qin & Liao Reference Qin and Liao2025) persistently oscillates with a large amplitude in the whole interval
$t\in [0,500]$
of simulation, say, almost always departs far away from zero, clearly indicating that the CNS solution almost always remembers its initial condition. In other words, the CNS results of the 3-D Kolmogorov flow (Qin & Liao Reference Qin and Liao2025) have no mixing. So, the CNS result (Qin & Liao Reference Qin and Liao2025) complies with the logical deduction from the SSIT mentioned earlier, say, the exact solution of NS turbulence has permanent memory of initial flow. However, the autocorrelation coefficient
$\rho (\tau )$
of the DNS solution of the 3-D Kolmogorov flow (Qin & Liao Reference Qin and Liao2025) quickly decreases to be close to zero, indicating that the DNS trajectory, which is quickly polluted by random numerical noise, forgets its initial flow quickly and thus has ‘mixing’. Unfortunately, this kind of ‘mixing’ of DNS violates the logical deduction from the SSIT mentioned earlier. The SSIT and our CNS results (Qin & Liao Reference Qin and Liao2025) highly suggest that spatial asymmetry, and/or more importantly, random disturbances, should be necessary for ‘mixing’ of NS turbulence (and even physical turbulence).
3.3.4. Topological restriction on the evolution
It should be emphasised that the SSIT imposes a topological straitjacket on the evolution of NS turbulence with symmetric initial flow conditions. This is for two reasons.
-
(i) Mandatory synchronisation: in a symmetric flow, vortices cannot merge randomly; they must merge in synchronised pairs or clusters to preserve the global topological structure.
-
(ii) Path locking: spatial symmetry dictates ‘allowed’ and ‘forbidden’ paths for energy transfer, effectively governing the rate and direction of the energy cascade in NS turbulence.
Note that mandatory synchronisation and path locking of exact solution of NS turbulence have been confirmed by CNS (Qin et al. Reference Qin, Yang, Huang, Mei, Wang and Liao2024; Liao & Qin Reference Liao and Qin2025; Qin & Liao Reference Qin and Liao2025; Qin et al. Reference Qin, Xu and Liao2026).
It should be emphasised that the above-mentioned interesting features of the exact solution of NS turbulence, such as ultra-chaos, the so-called noise-expansion cascade, permanent memory of symmetric initial flow and topological restriction on the evolution, can be directly deduced in logic from the SSIT. The fact that the CNS results (Qin et al. Reference Qin, Yang, Huang, Mei, Wang and Liao2024; Liao & Qin Reference Liao and Qin2025; Qin & Liao Reference Qin and Liao2025; Qin et al. Reference Qin, Xu and Liao2026) comply with the SSIT and all of these interesting features highly suggest that CNS could indeed reveal the mathematical truth concerning the exact solution of NS turbulence. Unfortunately, DNS results violate the SSIT and thus cannot reveal these mathematical truths, clearly indicating that DNS results would rapidly become badly polluted by numerical noise and thus depart far from the exact solution of NS turbulence for given initial conditions.
4. Discussion
First, let us consider the viewpoints of Corless (Reference Corless1994), winner of the Paul R. Halmos–Lester R. Ford Prize of the Mathematical Association of America, that ‘…there are four levels of abstraction used:
-
(i) the physical reality of the problem under study;
-
(ii) the continuous mathematical model of that physical reality;
-
(iii) the numerical discretisation of that mathematical model;
-
(iv) and the floating-point simulation of that discretisation.
Results from one level may or may not transfer easily to another level, and in particular, even qualitative features may not be preserved in that transfer.’ The present paper focuses on two aspects: the continuous mathematical model of physical reality (i.e. the NS turbulence governed by continuity equation and deterministic NS equations), and the numerical discretisation of that mathematical model and the floating-point simulation of that discretisation (i.e. CNS and traditional DNS of NS turbulence). As pointed out by Corless (Reference Corless1994), we should note that there may be certain numerical simulations that may not preserve even qualitative features of NS turbulence.
4.1. Scientific value of SSIT
Do numerical simulations of NS turbulence preserve features of NS turbulence? To answer this important question, we should consider key features of NS turbulence. For example, the most widely used feature of NS turbulence relates to the famous Kolmogorov −5/3 power law, which has a statistical form. The present paper reaches beyond this by providing a mathematical proof of the spatial symmetry invariance theorem (SSIT), whereby the spatial symmetry of the exact spatiotemporal solution of 3-D incompressible NS turbulence subject to periodic boundary conditions and symmetric external force maintains the same symmetry as its initial flow, forever. The SSIT has important, scientific value for two reasons.
-
(i) SSIT provides an ‘absolute scale’ for CFD: if a DNS fails to maintain spatial symmetry of a symmetric initial state, the result is mathematically invalidated as a true solution of NS turbulence.
-
(ii) SSIT provides the theoretical foundation (and necessity) for CNS, demonstrating that ultra-high precision is not just a preference but a requirement to capture the ‘mathematically clean’ trajectory of NS turbulence.
Our published CNS results confirm the SSIT by means of several benchmark tests. The detailed comparisons between CNS and DNS results strongly suggest that:
-
(i) the spatiotemporal trajectory of NS turbulence given by DNS is rapidly badly polluted by numerical noise, given that the DNS result quickly violates the SSIT;
-
(ii) loss of spatial symmetry in traditional DNS of NS turbulence is not a ‘physical instability’ but instead a ‘numerical artefact’ caused by unavoidable numerical noise and the chaotic nature of the turbulent flow;
-
(iii) the CNS results comply with the SSIT over a prescribed finite period of sufficiently long duration to reveal the mathematical ground-truth of NS turbulence, such as its ultra-chaotic characteristics, the noise-expansion cascade, the permanent memory of symmetric initial flow conditions and so on.
Notably, the NS equations lie at the heart of one of the famous Clay Millennium Problems. The SSIT could offer strategic insights into the Navier–Stokes Existence and Smoothness problem (Millennium Challenge) (Clay Mathematics Institute of Cambridge, Massachusetts 2000) for the following reasons.
-
(i) SSIT provides a constraint on numerical singularities and could restrict the search for potential singularities to within the symmetric function space.
-
(ii) The SSIT is a valuable diagnostic tool. If spatial symmetry is lost before a singularity occurs in a simulation, then the associated numerical breakdown is most likely a numerical artefact rather than a mathematical property of the NS equations.
-
(iii) The SSIT leads to the following non-uniqueness conjecture: due to the ultra-chaotic characteristic of NS turbulence, almost the same initial conditions admit distinct globally smooth solutions of NS equations with different spatial symmetries and statistics.
Note that, unlike Kida (Reference Kida1985) who used ‘high-symmetry’ as a computational tool, the SSIT represents a fundamental shift in paradigm for the following reasons:
-
(i) Kida focused on computational efficiency and used symmetry to reduce the size of the simulation domain by a factor of 64 to achieve higher Reynolds numbers with limited resources. Kida treated spatial symmetry as an observed utility rather than a proven law;
-
(ii) the SSIT establishes spatial symmetry as a theoretical benchmark. The focus is on the mathematical certainty that initial spatial symmetry must persist forever, say, in
$t\in [0,+\infty )$
. Any loss of spatial symmetry in a simulation is redefined not as physical instability, but as numerical pollution – providing mathematical evidence that the simulation has deviated from the true solution of the NS equations.
Expansion of numerical noise is essentially due to the chaotic property (trajectory instability caused by the uncontrolled growth of small disturbances) (Deissler Reference Deissler1986; Boffetta & Musacchio Reference Boffetta and Musacchio2001, Reference Boffetta and Musacchio2017; Berera & Ho Reference Berera and Ho2018; Ge et al. Reference Ge, Rolland and Vassilicos2023; Okamoto et al. Reference Okamoto, Ishihara, Yokokawa and Kaneda2025) of NS turbulence, but has nothing to do with spatial symmetry itself. This is the reason why DNS of laminar flow can maintain spatial symmetry of initial flow, given that laminar flow is insensitive to numerical noise. In other words, the numerical noise of NS turbulence certainly increases to macro-level no matter whether or not the initial flow is symmetric. It is due to the SSIT that we can use spatial symmetry to reveal properly the exponential expansion of numerical noise in DNS. Similarly, the ‘noise-expansion cascade’ is also a fundamental property of NS turbulence for general flows. Thus, the ‘noise-expansion cascade’ is a fundamental property of NS turbulence for general initial flows, say, no matter whether or not the initial flow is symmetric. Note that it is the SSIT that provides a mathematical cornerstone for the ultra-chaotic characteristics, the so-called noise-expansion cascade, and permanent memory of initial flow condition of NS turbulence.
It is well known that a numerical solution of a uniformly hyperbolic system may depart from the exact trajectory while still shadowing a true trajectory of the same dynamical system (the shadowing lemma) (Anosov Reference Anosov1967; Bowen Reference Bowen1975). However, hardly any chaotic system is uniformly hyperbolic and the non-existence of such shadowing trajectories may be caused by the finite-time Lyapunov exponent fluctuating about zero (Dawson et al. Reference Dawson, Grebogi, Sauer and Yorke1994). In addition, no true trajectory of reasonable length can be found to exist near any computer-generated trajectories for non-hyperbolic chaotic systems with unstable dimension variability (Do & Lai Reference Do and Lai2004). It seems that the good promise suggested by the Shadow Lemma (Anosov Reference Anosov1967; Bowen Reference Bowen1975) for a uniformly hyperbolic system does not work for many non-hyperbolic systems (Anishchenko et al. Reference Anishchenko, Kopeikin, Kurths, Vadivasova and Strelkova2000). Most investigations on shadowing are based on low dimensional dynamic systems. It has been found (Sauer Reference Sauer2002) that, even for a low-dimensional dynamic system with only two degrees of freedom, extremely small levels of numerical noise might generate macroscopic errors in simulation average statistics that are several orders of magnitude larger than the noise level. To the best of the author’s knowledge, practical shadowing algorithms have not been previously proposed for spatiotemporal chaotic systems such as turbulent flows governed by nonlinear partial differential equations (PDEs). In short, it is certain that the spatiotemporal trajectory of DNS of NS turbulence subject to a given initial flow should be quickly badly polluted by numerical noise.
4.2. Some open questions
Note that the mathematical proof of the SSIT is based on the assumption that the solution is analytic with a non-zero Taylor series in
$t\in [0,+\infty )$
. For the time scales of interest in most physical processes, this assumption holds and nowadays is the foundation of CFD; thus, the SSIT provides the only rigorous bound available. However, it is an open question whether or not this assumption is absolutely necessary for the proof of the spatial symmetry invariance of NS turbulence.
Ergodicity requires that time averages are the same for almost all initial conditions; in other words, the flow should rapidly forget its initial flow conditions; however, this inherently fails if different initial conditions lead to different attractors (corresponding to different statistics). Unfortunately, NS turbulence has an ultra-chaotic property, which can be deduced directly using logic from the SSIT, as mentioned in § 3.3.1, meaning that different initial conditions might lead to different attractors (and statistics). In addition, it can be logically deduced from the SSIT that NS turbulence has permanent memory of symmetric initial flow condition, as mentioned in § 3.3.3. Thus, according to the SSIT, the exact solution of NS turbulence subject to spatially symmetric initial flow should have no ergodicity. This is a very surprising mathematical conclusion based on the SSIT.
The mathematical proofs of ergodicity given by Kuksin & Shirikyan (Reference Kuksin and Shirikyan2001, Reference Kuksin and Shirikyan2002, Reference Kuksin and Shirikyan2012) for the 2-D NS equations rely on random forcing to provide the ‘mixing’ that destroys memory of initial conditions. Random white noise acts as a powerful mixing agent: it continuously injects energy across the spatial modes on which it acts and forces the system to explore all possible states. A deterministic force, even if chaotic, does not provide this uniform stochastic exploration. Without such randomness, one cannot guarantee that correlations decay to zero. Thus, as pointed out by Kuksin & Shirikyan (Reference Kuksin and Shirikyan2001, Reference Kuksin and Shirikyan2002, Reference Kuksin and Shirikyan2012), the random forcing is a physically essential ingredient that breaks the system’s symmetry, destroys multiple attractors, and enforces a unique, mixing stationary state. Without it, the system is generally not ergodic, and its long-term statistics depend critically on the initial condition and the exact form of the forcing. The permanent memory of NS turbulence concerning symmetric initial flow conditions, which is revealed herein by the SSIT, complies with the mathematical proofs of ergodicity given by Kuksin & Shirikyan (Reference Kuksin and Shirikyan2001, Reference Kuksin and Shirikyan2002, Reference Kuksin and Shirikyan2012) for the 2-D NS equations, since it clearly indicates the importance and necessity of randomness and/or spatial asymmetry in initial flow to ergodicity of NS turbulence, which is a cornerstone of turbulence theory. However, it is still an open question whether or not random disturbance in initial flow and/or external force should be absolutely necessary for three-dimensional NS turbulence to have ergodicity.
Corless (Reference Corless1994) pointed out that ‘…there are four levels of abstraction used: the physical reality of the problem under study, the continuous mathematical model of that physical reality, the numerical discretisation of that mathematical model, and the floating-point simulation of that discretisation. Results from one level may or may not transfer easily to another level, and in particular, even qualitative features may not be preserved in that transfer.’ Note that all CNS results comply with the SSIT over a prescribed simulation period that is finite but long enough to reveal mathematically exact solutions of NS turbulence, and thus may preserve qualitative features of NS turbulence such as the spatial symmetry invariance, ultra-chaotic property, noise-expansion cascade and permanent memory of the symmetric initial flow condition of NS turbulence. Unfortunately, the SSIT is violated by all current state-of-the-art DNS of NS turbulence (Qin et al. Reference Qin, Yang, Huang, Mei, Wang and Liao2024; Qin & Liao Reference Qin and Liao2025; Qin et al. Reference Qin, Xu and Liao2026), indicating that DNS may not preserve qualitative features of NS turbulence. Thus, the SSIT as a theoretical benchmark reveals the non-preservation and violation between numerical simulation given by traditional DNS and these qualitative features of exact solution of NS turbulence. How to eliminate such kind of violation and non-preservation of traditional DNS for NS turbulence? This is an open question for future research.
In addition, as mentioned by Corless (Reference Corless1994), the real turbulence (including physical experiments of real turbulent flow) and its mathematical model (i.e. the NS turbulence here) are two levels of abstraction. Similarly, the SSIT as a theoretical benchmark of NS turbulence may be used to verify whether the features of real turbulent flows may be preserved in NS turbulence or not. In other words, it is an open question whether or not NS turbulence as a mathematical model of turbulence may preserve all features of real turbulence. Note that this paper focused on NS turbulence governed by the continuity and deterministic NS momentum equations, which neglect all small stochastic disturbances, including artificial numerical noise and physical disturbance. It is also an open question whether these small stochastic disturbances may exert significant influences on real turbulent flows so that it would be better to use stochastic NS equations rather than deterministic ones.
5. Conclusion
A full mathematical proof was given of a spatial symmetry invariance theorem (SSIT) for the incompressible fluid flow equations. The theorem states that the exact solution of the continuity and deterministic NS momentum equations with periodic boundary conditions maintains the same spatial symmetry as that of the symmetric initial flow, forever. The SSIT has significant scientific value in the fluid mechanics of turbulent flows for the following reasons.
-
(i) Using the SSIT, it is straightforward to deduce several key properties of NS turbulence, including ultra-chaos, the so-called noise-expansion cascade, permanent memory of symmetric initial flow conditions and so on. All of them are confirmed and verified by CNS.
-
(ii) SSIT provides an absolute scale by which to check the validity and accuracy of numerical algorithms used to simulate NS turbulence.
-
(iii) It was found that CNS complies with SSIT and the interesting properties deduced from it mentioned in (i). Thus, CNS indeed could obtain true spatiotemporal trajectories of NS turbulence in a finite but long enough period of time. CNS provides a reliable numerical tool to investigate NS turbulence.
-
(iv) It was found that DNS quickly violates SSIT, revealing that the spatiotemporal trajectory of NS turbulence predicted by traditional DNS soon becomes badly polluted by numerical noise. Importantly, it is found that, under certain circumstances, even the statistical results of DNS might be unstable, i.e. exhibit sensitivity dependence on the initial flow conditions. In other words, the loss of spatial symmetry in traditional DNS of NS turbulence is not a physical instability but instead a numerical artefact caused by the chaotic nature of the flow. This implies that, under certain circumstances, modellers should be very cautious when taking traditional DNS results as benchmark solutions for NS turbulence.
-
(v) SSIT and the interesting properties of NS turbulence mentioned in (i) may be used as a theoretical benchmark to check and verify (by comparison with physical experiments) whether or not deterministic NS equations as a mathematical model of turbulence may indeed preserve all features of real turbulence.
Thus, SSIT provides a new mathematical cornerstone of NS turbulence. SSIT also gives us a novel way of reconsidering the foundation of the deterministic NS equations in turbulence modelling. SSIT is therefore capable of greatly deepening our physical insight concerning turbulence. In the future, SSIT could provide a completely new theoretical cornerstone and a completely novel perspective regarding turbulence, considered by Feynman to be the most important unsolved problem of classical physics. Some open questions are proposed for future investigations.
Acknowledgements
Sincere thanks to the anonymous reviewers for their valuable comments.
Funding
This work is partly supported by National Natural Science Foundation of China (No. 12521002).
Declaration of interests
The authors declare no competing interests.
Data availability statement
All data are available by sending requirement to the corresponding author.
s1
s2
u(x,y,z,t)
v(x,y,z,t)
w(x,y,z,t)
s3
δ1=10−20
δ2=10−40
s1∘ψ(x,y,t)↦ψ(x+π,y+π,t)
s2∘ψ(x,y,t)↦ψ(−x,−y,t)
ρp(τ)
p(x,y,z,τ)
p(x,y,z,0)
ρ=0
Open access
