From Navier-Stokes millennium-prize problem to soft matter computing

Navier-Stokes equations, are Newton’s laws of motion to model the motion of a fluid, in two or three spatial dimensions ($\mathcal{R}^{2}\textrm{or}\,\mathcal{R}^{3}$), and single time $t$ dimension. The velocity vector field $\vec{u}(x,t)=(u_{i}(x,t))_{1\leq i\leq n}\in\mathcal{R}^{n}$ and scalar pressure field $p(x,t)\in\mathcal{R}$ are the unknowns in the equations. Given an assumption of incompressibility of the fluid, these equations are given by

where $\nu$ is the kinematic viscosity (units: $m^{2}s^{-1}$) and $\vec{f}$ is the external body force. The initial conditions of the velocity vector field are given as

where (1) and (2) are defined for $t\geq 0$ and $x\in\mathcal{R}^{n}$. If $\nu=0$, then the fluid theoretically has no viscosity (inviscid fluid), and the equations are termed Euler equations. The solutions to either of these equations are considered physically reasonable if

for all $t\geq 0$. Here, $C^{\infty}$ refers to the functions $p,u$ being smooth and infinitely differentiable, and the kinematic energy of the fluid is bounded above by a constant $C$.

Attempts to find solutions to the Navier-Stokes equations (NSE) goes back to early 20th century, long before the problem, stated above, was formally formulated. In 1911, Oseen (1911) resolved a paradox posed by George G. Stokes known as Stokes’ paradox. Stokes’ paradox stated that there are no bounded (finite) solutions to the Stokes equations (Reynolds number is negligibly small) for the case of a infinitely long cylinder or a disk in 2D falling through a fluid. Oseen resolved the paradox by taking into account the inertial term ($u_{i}\partial_{i}\vec{u}$) and a singular point force (called stokeslet). Stokeslet allows the velocity to grow as the logarithm of the radial distance, and thus incorporated as the inner boundary condition to yield bounded solutions $\vec{u}(x,t)$ outside the cylinder. In 1934, Leray (Leray, 1934) found a set of solutions (called weak solutions) $u_{0}$ belonging to $L^{2}$ norm¹¹1An $L^{2}$ norm or the Euclidean norm of a function gives the length of a vector $\mathbf{x}=(x_{1},x_{2},...,x_{n})$ as: $||\mathbf{x}||_{2}=\sqrt{x_{1}^{2}+x_{2}^{2}+..+x_{n}^{2}}$, for a manifold $\mathcal{R}^{n}$ such that the weak solution in addition to energy inequality shall be unique if a strong solution exists (or the weak solution is regular enough to be a strong solution). However, this dream has not been achieved till date because weak solutions have not been proved yet to be unique, thus regular, in nature. The physical meaning of such a partial regularity is that the solutions $\vec{u}(x,t)$ match analytically well with the smooth solutions only in the large regions of spacetime, which excludes the regions of null sets (measurable sets of measure zero) ²²2Cantor sets are one kind of null sets. corresponding to time $t$ and tiny length scale $x$. Such null sets are finite sets of arbitrary small length, and are usually disconnected, or in other words, discrete. Hence, it is possible to have Leray weak solutions which would be different at different times $t$ (only because of unknown regularity in the null sets), leading one to have non-unique general class of solutions.

In 1993, Constantin, Fefferman Constantin and Fefferman (1993) proved that, for given constants $\Omega$ and $\rho$, if the angle between unit vorticity vectors at positions $x$ and $y$ is $\varphi(x,y,t)$, then the magnitude of the vorticity $\omega$ at both locations is above $\Omega$ if

such that the vortex tubes are formed in an antiparallel fashion osculating (observed physically and numerically in experiments) each other. The vorticity vector is regular only in $L^{2}$ norm and not until $L^{\infty}$ norm. As a result, this well-behaved solution of vorticity vector (in the sense of parallel or antiparallel alignment) exists only because of the viscous mechanism, leading to partial regularity within the nonlinear stretching mechanism.
One type of scaling of NSE is given as

and there are various partial regularity results known for this. One of the most prominent one is called Ladyzhenskaya-Prodi-Serrin conditions ((Ladyzhenskaya and Kiselev, 1957), (Prodi, 1959), (Serrin, 1961)), which states that if the Leray weak solution lies in $L^{p}_{t}L^{q}_{x}$ with $2/p+3/q\leq 1$, then the solution is unique and smooth in positive time. The endpoint in the above equality for $p=\infty$ and $q=3$, given for $L^{\infty}_{t}L^{3}_{x}$, was proved by Escauriaza-Seregin-S̆verák Escauriaza et al. (2003). It should be noted that while the natural scaling of weak solutions has a regularity for the condition, $2/p+3/q=3/2$, the energy equality holds in NSE with an additional regularity for the condition, $2/p+3/q=5/4$ (Shinbrot, 1974). This means that there is a gap open between the natural scaling of the equations and the kinetic energy, which can lead to non-uniqueness of weak solutions. In fact, this nonuniqueness was hinted by Jia-S̆verák (Jia and Šverák, 2014), where the authors proved that in the class $L_{t}^{\infty}L_{x}^{3,\infty}$, Leray solutions are nonunique if a certain spectral assumption holds for a linearised Navier-Stokes operator. In 2019, Buckmaster and Vicol (Buckmaster and Vicol, 2019) proved a stronger result that the weak solutions to the 3D NSE, $v\in C_{t}^{0}H_{x}^{\beta}$ (where $H^{\beta}$ is the $L^{2}$ space with the regularity index $\beta$), are nonunique. The nonuniqueness refers to ill-posedness and existence of infinitely many solutions so that a clear statement about the nature of the solution can not be made. The idea used in the proof of this work is the convex integration method (De Lellis and Székelyhidi, 2015) and the authors expect that these tools might in the future establish nonuniqueness of Leray weak solutions.

While the results shown above are obtained by functional analysts for past one century, the NSE question has also percolated the community of applied mathematicians. Kimura and Moffatt (Moffatt and Kimura, 2019) showed that two vortex rings, of finite radii and counter-rotating circulation $\Gamma$ (Fig. 2(a)), when approaching each other can exhibit a finite-time singularity, which is physical (not mathematical) in nature, when the key length scale of the system becomes smaller than the intermolecular scale, making the continuum hypothesis ³³3A fundamental physical assumption in the derivation of NSE, continuum hypothesis states that the the effect of statistical and thermal fluctuations because of the motion of molecules in a fluid are not necessary to be taken into account while deriving continuum laws of motion, given that the length scale of system is at least 10 times greater than the mean free path of molecules in the continua. See (Batchelor, 2000, see p. 6) becomes invalid. In fact, such finite-time physical singularities are fairly ubiquitous in household phenomena such as: the pinch-off of pendant drop (Fig. 2(b).i); fingering patterns in a Hele-Shaw flow cell consisting of parallel plates separated by a distance (Fig. 2(b).ii); eddies formed at an edge of a wedge-shaped container driven by the rotation of a cylinder on the top of the container (Fig. 2(b).iii); slender-like filaments formed by a receding drop of biological fluid secreted by pitcher plants (Fig. 2(b).iv); finite-time singularity formed by an unstable thin-film (Fig. 2(b).v). It is to be noted that the case of vortex reconnection has also been considered by functional analysts to check if the NSE solutions are smooth and regular in nature or not. It turns out that the vortex reconnection has two stages. The first stage is when the strong coherent vortices come together and stretch to form intense narrow vortex tubes which approach at various angles by violating the smoothness assumption on the gradient of direction $\nabla\zeta$. The second stage is when the viscosity plays an important role by having an upper bound

on the vorticity. As a result, the vortex tubes become antiparallel, reconnect, and finally decay. While this picture has been numerically shown in several studies (Meiron et al., 1989), there is no rigorous proof on the regularity of solutions in this problem.

In Section 2, numerous efforts to either prove uniqueness or hint towards nonuniqueness of the solutions to NSE were outlined. This section, specifically, highlights a new scheme to solve the problem by using ideas from theoretical computer science.

where the authors constructed solutions for the steady Euler flow (without viscosity) on a Riemannian 3-sphere $\mathcal{S}^{3}$ that are Turing-complete. Here, Turing completeness means that for given points on the sphere, the problem of bounding the dynamical trajectory of those points, as the equation evolves, is an undecidable problem. Following this work, the authors extended the analysis to a standard Euclidean 3-dimensional space. The tools used to prove such a result are contact topology, symplectic geometry, h-principle, generalised shift maps, and Cantor sets. The key step in the construction of their proof is translating the problem in Fluid Dynamics as a problem in Geometry by using a key master point on establishing a correspondence between stationary Euler flows (Beltrami fields) and Reeb vector fields - a former result by Etnyre, Ghrist and also envisaged formerly by Dennis Sullivan (2022 Abel Prize winner (Papadopoulos, 2022)). Such a mirror is used in the proof to construct a Reeb vector field in dimension three. This construction of the Reeb vector field reminds Poincare’s section and promotes a 2D to 3D construction proving that a certain map in the plan can be seen as time-1-map of the flow of a Reeb vector field. The 2D construction builds on former work of Cristopher Moore in 1991 (Moore, 1991) on whether hydrodynamics is able to perform computations, so the construction in this work answers Moore’s question. An intuitive sense of how the proof is derived can be given as follows. For a given string of symbols of the alphabet $(t^{*}_{-k},....,t^{*}_{k})$ and a Turing machine $T$ with an input tape $t$, the authors showed that it is possible to construct explicitly a point $p$ and an open set $U$ in $\mathcal{S}^{3}$ such that when the velocity field $X$, which is the solution of inviscid Navier-Stokes equations, passes through $p$, it also intersects $U$ if and only if $T$ halts with the output tape at positions corresponding to symbols $(t^{*}_{-k},....,t^{*}_{k})$ (Fig. 3(d)). This work marks as an important step in towards treating the general solution of different subclasses of NSE as embeddable in a Turing machine. The authors found that if the initial vector field $\vec{X}$ for Euler equations is given as an input in the NSE, the resulting velocity field can only simulate a finite number of steps, and thus, the Turing completeness of Navier-Stokes equations would still remain unproven.

The results presented above are only for manifolds which are of higher dimensions, and thus, are not scalar invariant. Proving finite-time blowup is an issue for such manifolds, as at finer scales, the metric structure becomes flat and the equivalence between Euler flow and arbitrary vector fields does not hold true. The key next step here is to find solutions to Euler equations on $\mathcal{R}^{3}$ which is scalar invariant, as the original Clay problem ((1)) states. Tao discussed in MSRI seminar that certain type of hybrid manifolds which exhibit scalar invariance in one region and universality in other region might potentially be a way to move forward to prove Euler flows Turing-complete in $\mathcal{R}^{3}$ be importantly noted that Turing-completeness of Euler equations in $\mathcal{R}^{3}$ is not a straightforward path as it seems, as once it is established that the equations are Turing-complete, the proof of finite-time blowup will require lots of care in the proof.

As discussed in Tao (Tao, 2016), proving the Turing-completeness of the Euler or Navier-Stokes equations is the first step towards building, designing, and quantifying the design of logic gates that can be assembled to simulate the classical laws of conservation of energy and helicity that the fluid equations will obey. After making sure that the physics is right and design of logic gates is rigorous enough, such a “fluid program” can be shown to blowup for a specific set of initial data, and the debate of regularity versus blowup can be given more concrete answer, for a given choice of initial data.

As discussed by George Batchelor (Batchelor, 2000), the mathematical basis for the continuum treatment of liquids in motion is incomplete and thus, it is difficult to deduce the properties of a hypothetical continuous medium that moves the same way as a real fluid. Very recently, there is a growing evidence that ‘fluids’ do not much show the ‘regular’ or symmetric behaviour, at certain regimes where the flow is turbulent or the length scale is close to the molecular scale. We recently argued that there lacks a fundamental framework which connects the physics of fluids at two different scales Sharma et al. (2022). If the analysis of fluids is restricted only to one continuum scale without incorporating the underlying molecular contributions (Bandak et al. (2022)), then the Navier-Stokes equations yield singular solutions in an unpredictable fashion. To solve this problem, the topic of “emergence” within complexity sciences field comes to our rescue as it celebrates the idea that physics at multiple scales can be correlated or conceptually connected to each other. One framework built off complexity view of fluid mechanics suggested by us (Sharma et al. (2022)) is the use of cantor sets (Lindenmayer systems) as a computational framework to derive continuum level properties from molecular level. These systems consists of certain (recursive) rules being applied at one scale to yield the physics at another scale.

If one accepts the use of L-systems or related computational frameworks in analysing fluids at multiple scales, then another corollary from this result is to start accepting that fluids, in fact, ((Wheeler, 2018)) consist of information (bits) at their core. This changes the way we see fluids; instead of fluids being modelled mathematically, they can be computationally programmed. Such a program can be both experimental (Adamatzky (2019); Brun (2022)) and virtual (Sharma et al. (2022); Cardona et al. (2021)) in nature.

Given that the fundamental nature of fluid dynamical behaviour is computational in nature, the question that immediately comes up is: is it digital or analog? To answer this question, it is important to resource two important results. First is that the the fluid dynamical equations are chaotic in nature, as discussed extensively by Doering & Gibbons (Doering and Gibbon, 1995, see Ch. 5). Second is that the chaotic behaviour of a physical system can be encoded in Siegelmann (1995) in analog shift maps which are more powerful than Turing machine, and computes similar to how a neural network and an analog machine would compute. These two results can be combined to argue that a fluid can be programmed to behave like an analog computer, and perform computations similar to a neural network. This result leads one to explore new avenues of research, particularly in neuromorphic and unconventional computing communities, as the tools from these fields can be borrowed to program the analog behaviour of fluids, soft matter, or active matter.

In the increasing community of scientists working on analog and unconventional computing, neuromorphic computing (NMC) is the branch that has been gaining the most attention since the beginning of the fourth industrial revolution Its main objective is to overcome the challenge of engineering machines based on brain-like algorithms, namely, capable of solving problems following a holistic approach, a peculiarity exclusive to the human brain.

The way NMC addresses this challenge involves both cutting-edge software architectures and novel hardware paradigms, in a never-ending interchange. Indeed, while it is possible to trace the path that led the scientific community from the need for efficient big-data analysis to the development of artificial neural networks (ANNs) (Fernando and Sojakka, 2003), nowadays separating the development of software and hardware in NMC is becoming an impossible task. For this reason, modeling fluid dynamics using NSE and analyzing the computing power of a liquid state machine (LSM) can be seen as two sides of the same coin, a concept that justifies the presence of this section in this paper, and on which the authors are actively working and will show original results in forthcoming publications.

From an algorithm perspective, the two most important architectures of NMC are spiking neural networks (SNNs) (Fernando and Sojakka, 2003), and reservoir computing (RC) (Fernando and Sojakka, 2003). SNNs are ANNs based on time-dependent synapses. At variance with an ensemble of perceptrons (basic units of ANNs with output defined as $y=f(\sum_{i}w_{i}x_{i}+b)$, where $f$ is the activation function, $x_{i}$ are the nodes, $w_{i}$ the corresponding trainable weights, and $b$ is the bias term), the neurons of a SNN are allowed to transmit information only when the activation potential crosses its threshold.

RC represents a subclass of recurrent neural networks (RNNs) with training performed only at the readout stage. In other words, RC can be considered as a RNN where the total ensemble of hidden layers is flattened into a sparse highly-dimensional layer with fixed dynamics (i.e., a black box). It benefits from the RNN’s Turing completeness in most of its architectures, excluded the limit ones, like the extreme learning machine (ELM) which is a RC feedforward version. The most widespread RC paradigms are echo state networks (ESNs), LSM, ELM, and all RC quantum versions, outside the purpose of this paper. In terms of equations, we can generalize their models to

with $\mathbf{y}(t_{j})$ the network readout at time $t_{j}$ satisfying the condition

where $\mathbf{y}_{T}$ is the target output and $\epsilon\lesssim 0$. In Eqs. (8), $\mathbf{W}^{out}$ is the matrix of trainable weights obtained from Eq. (9) as

where $\dagger$ stands for the Moore-Penrose pseudo-inverse matrix; $\mathbf{x}(t_{j})$ and $\mathbf{u}(t_{j})$ are the vectors of nodes in the output and input layers at time $t_{j}$, respectively; $\mathbf{b}$ is the bias vector; $f$ is the activation function; $\mathbf{W}^{s}$ and $\mathbf{W}$ are the matrices of the non-trainable weights of the reservoir.

The ELM is obtained from Eqs. (8) when $\mathbf{W}=\mathbf{0}$, the limit case characterized by absence of recurrence. It has already been shown that the nonlinear Schrödinger equation can model an ELM (Marcucci et al., 2020), a fact that suggests that also a reservoir governed by the NSE is able to perform interpolation, classification, and logic gates. Indeed, the nonlinear Schrödinger equation can be derived from the NSE in deep-water approximation (Peregrine, 1983). More precisely, in Marcucci et al. (2020) (Marcucci et al., 2020), the authors demonstrated how a continuous-wave laser beam propagating in a self-focusing crystal can realize a neuromorphic computer by encoding information in a low-intensity structured beam and exciting the crystal nonlinearity through a bias beam with rectangular spatial intensity distribution. Except for the presence of the structured beam, this model describes also the hydrodynamics of the dam break problem (El et al., 2016), tracing the path for further investigations in the field.

When Eqs. (8) are modified in a way that $\mathbf{x}$ represents spike neurons, we obtain an LSM (Maass, 2011), a NMC architecture proven to be a universal approximator. A hardware implementation of the LSM was realized at the beginning of this millennium (Fernando and Sojakka, 2003). This system uses mechanical excitation of linear waves in a bucket filled with water and, by imaging the wave interference at the center of the bucket, it accomplishes complex tasks such as vowel recognition or logic operations. Even if the understanding of this machine could clearly benefit from an analysis of its reservoir in terms of NSE, no research has been done in this direction because of the complexity of the system. However, we believe that such an investigation would be an important milestone to achieve in order to develop better models of the NMC’s overall performance.

An answer to the NSE regularity problem will probably be the result of influx of ideas from other disciplines such as, theory of computation, theory of neural networks, and computational complexity. Fluid mechanics has been a topic of interest for engineers, physicists, mathematicians, biologists, and computer. The turn of 20th century saw advancements in the mathematical side when questions of existence and uniqueness of PDEs were asked, particularly by Ladyzhenskaya, Leray, and Oseen. There was a shift in thinking, in the mathematical physics community, from proving/finding a classical solution of a problem towards accepting the notion of a ‘weak solution’ of the PDE. On the other hand, the experimental side of fluid mechanics, over the last few decades, has gained immense traction because of rapid developments in high resolution microscopy, high-speed imaging, experimentation in adverse conditions like space station, alongside the influences of molecular chemistry, polymer physics, quantum physics, and many other disciplines. The vast amounts of data available from experimental observations has invited an array of diverse neural network architectures to model the phenomenon. However, there still remains a lot to be done on the fundamental side of neural networks: there is a need for the theoretical-computational theory of fluid mechanics that can support ongoing remarkable feats of neural networks in modeling the experimental observations. In fact, a recent paper, Wang, Lai, Gomez-Serrano, Buckmaster (Wang et al., 2022) showed that a physics-informed neural networks (PINN) can be employed to find smooth self-similar solutions for the Boussinesq equations. The authors argue that their numerical approach, introduced in their paper can guide a computer-assisted proof of finite-time singularity formation from finite-energy initial data. In fact, very recently, 2D Boussinesq and 3D Euler equations were proven to blow-up in self-similar fashion (Chen and Hou (2022)). The technique involved decomposing the solution into semi-analytical construction that captures the long-range behaviour and the numerical computation that yields the behaviour of the next time step, such that that combination of the two is useful to prove self-similar blowup of the solutions. This result present first of its own kind of computer-assisted proof of a blowup scenario. It is the hope that in the future, fluid equations without any boundary or support might well be proved to exhibit blowup for smooth initial data, however, that probably would involve highly computer-assisted programming of fluid parcels.

Given that the Navier-Stokes regularity problem lies at the heart of current scientific advances of the 21st century, it is highly likely that the theoretical and numerical advances in physics-informed machine learning techniques and its connections to the computational complexity of ‘computable’ fluid parcels are going to metamorphosise the physics of fluids, soft matter, and active matter in coming few decades and potentially lead to emergence of a new discipline along the lines of “soft matter complexity, computing, and learning”.

We thank Prof. Eva Miranda for additional references and points to add to this article; David Kaloper Meršinjak for long hours of discussion on various related problems that led to establishing the narrative of this article; Dr Shamit Shrivastava for the deeply insightful scientific conversations and for sharing with us his enthusiasm in understanding the complexity of the world.