The instanton method and its numerical implementation in fluid mechanics

The functional integral description of stochastic systems started in the seventies with the MSRJD formalism [41–43]. The method transforms the stochastic system into a functional integral representation. Although there is little hope to solve the path integral analytically, it is the starting point for systematic perturbation theory in stochastic dynamical systems (see especially [44] for applications in critical dynamics). It is also the basis for the instanton calculus for stochastic dynamical systems. To introduce the formalism, consider a SPDE of the form

$$\begin{eqnarray}&&\dot{u}=b[u]+\eta (x,t)\end{eqnarray} \tag{ 1 }$$

in d space-dimensions, u having n components: $u(x,t):{{\mathbb{R}}}^{d}\times [-T,0]\to {{\mathbb{R}}}^{n}$ and η being a Gaussian forcing with correlation

$$\begin{eqnarray}&&\langle \eta (x,t)\eta (x+r,t+s)\rangle =\chi (r)\delta (s),\end{eqnarray} \tag{ 2 }$$

i.e. white in time and some correlation function $\chi (r)$ in space.

In fluid dynamics applications, the exact form of the spatial correlation is often irrelevant and can be characterized solely by its amplitude $\chi (0)$ and correlation length L, where $\chi (r)$ does not change significantly for $x\lt L$ and decays to zero quickly for $x\gt L$. Dimensionally, for example $L=\sqrt{\chi (0)/{\chi }^{\prime\prime }(0)}$ (assuming the forcing is isotropic and χ depends only on the scalar distance r). The forcing is then completely characterized by $\chi (0)$ and ${\chi }^{\prime\prime }(0)$. Furthermore, we will only be considering additive noise. The drift b in general is a nonlinear, non-gradient operator. A (quick) derivation (similar to [45]) of the MSRJD formalism follows by considering the expectation of an observable $O[u]$ and writing this expectation as a path integral taken over all noise realizations (using the fact that η is Gaussian and δ-correlated in time). For a more rigorous discussion of path integrals, in particular in the presented context, we refer the reader to the works [46, 47]. Here, we keep in mind that the field $u[\eta ]$ has a functional dependence on the forcing η implicitly given by equation (1):

$$\begin{eqnarray}&&\langle O[u]\rangle =\int {\mathcal{D}}\eta \;O[u[\eta ]]{{\rm{e}}}^{-\int \langle \eta ,{\chi }^{-1}\eta \rangle /2{\rm{d}}t},\end{eqnarray} \tag{ 3 }$$

where $\langle \cdot ,\cdot \rangle$ is an appropriate inner product, e.g. L². At this stage, we can consider the transformation of the noise paths to the paths of the field u given by $\eta \to u$ given by (1), hence $\eta =\dot{u}-b[u]$. This coordinate transform leads to a Jacobian in ${\mathcal{D}}\eta =J[u]{\mathcal{D}}u$ with

$$\begin{eqnarray}&&J[u]=\mathrm{det}\left(\displaystyle \frac{\delta \eta }{\delta u}\right)=\mathrm{det}\left({\partial }_{t}-\displaystyle \frac{\delta b}{\delta u}\right).\end{eqnarray} \tag{ 4 }$$

Performing this coordinate change results in the Onsager–Machlup functional [48, 49]

$$\begin{eqnarray}&&\langle O[u]\rangle =\int {\mathcal{D}}u\;O[u]J[u]{{\rm{e}}}^{-\int \langle \dot{u}-b[u],{\chi }^{-1}(\dot{u}-b[u])\rangle /2{\rm{d}}t}.\end{eqnarray} \tag{ 5 }$$

This formulation is the starting point for directly minimizing the Lagrangian action

$$\begin{eqnarray}&&{S}_{{\mathcal{L}}}[u,\dot{u}]=\displaystyle \frac{1}{2}\int \langle \dot{u}-b[u],{\chi }^{-1}(\dot{u}-b[u])\rangle \;{\rm{d}}t.\end{eqnarray} \tag{ 6 }$$

Since it is often more convenient to work with the original correlation function χ instead of working with its inverse, we delay this coordinate change and introduce an auxiliary field μ and obtain (by virtue of the Fourier transform, completion of the square) from (5):

$$\begin{eqnarray}&&\langle O[u]\rangle =\int {\mathcal{D}}\eta \;{\mathcal{D}}\mu \;O[u[\eta ]]\;{{\rm{e}}}^{-\int [\langle \mu ,\chi \mu \rangle /2-{\rm{i}}\langle \mu ,\eta \rangle ]\;{\rm{d}}t}.\end{eqnarray} \tag{ 7 }$$

Now again considering the coordinate change $\eta \to u$ we arrive at the path integral representation of $O[u]$

$$\begin{eqnarray}&&\langle O[u]\rangle =\int {\mathcal{D}}\eta \;{\mathcal{D}}\mu \;O[u]J[u]\;{{\rm{e}}}^{-S[u,\mu ]},\end{eqnarray} \tag{ 8 }$$

with the action function $S[u,\mu ]$ given by

$$\begin{eqnarray}&&S[u,\mu ]=\int \left[-{\rm{i}}\langle \mu ,\dot{u}-b[u]\rangle +\displaystyle \frac{1}{2}\langle \mu ,\chi \mu \rangle \right]\;{\rm{d}}t,\end{eqnarray} \tag{ 9 }$$

also termed the MSRJD response functional. In many cases it can be shown explicitly that the Jacobian J(u) is not relevant to the discussion as it reduces to a constant, the value of which depending on the choice of either Itō or Stratonovich discretization (see e.g. [45]).

Two remarks should be added: (i) the change introducing an auxiliary field μ to linearize the action with respect to the forcing is also known as Hubbard–Stratonovich transformation [50, 51] and (ii) the MSRJD action S is closely related to classical limit of the Keldysh action [52–54].

The main idea of the instanton method is to compute the dominating contribution to this path integral via a saddle point approximation, i.e. by finding extremal configurations of the action functional (9).

Let us be more specific and consider an observable

$$\begin{eqnarray}&&O[u]=\delta (F[u(x,t=0)]-a)\end{eqnarray} \tag{ 10 }$$

for a scalar functional $F[u]$, which is defined at the end point t = 0 of a path that started at $-T\lt 0$ (keeping in mind that we might be interested in the limit $T\to \infty$). The idea is that we observe an extreme event at t = 0 that has been created by the noise (and we give the noise infinite time to create the extreme event). In a turbulent flow, for example, an extreme event of interest could be a large negative gradient of the velocity profile (i.e. $F[u]={u}_{x}\delta (x)$), an event of high vorticity, an event of extreme local energy dissipation, etc.

Seeking a path integral representation of the probability distribution $P(a)$ for the events that fulfill $F[u]=a$ at t = 0, we find

$$\begin{eqnarray*}P(a) & = & \langle \delta (F[u]\delta (t)-a)\rangle \\ & = & \displaystyle \int {\mathcal{D}}\eta \;{\mathcal{D}}\mu \;\displaystyle \frac{1}{2\pi }{\displaystyle \int }_{-\infty }^{\infty }{\rm{d}}\lambda J[u]\;{{\rm{e}}}^{-S[u,\mu ]}\;{{\rm{e}}}^{-{\rm{i}}\lambda (F[u]\delta (t)-a)}\end{eqnarray*}$$

By now computing functional derivatives, we find

$$\begin{eqnarray*}\displaystyle \frac{\delta S}{\delta \mu } & = & -{\rm{i}}\left(\dot{u}-b[u]\right)+\chi \mu \;\\ \displaystyle \frac{\delta S}{\delta u} & = & {\rm{i}}\dot{\mu }+{\rm{i}}{({\nabla }_{u}b[u])}^{T}\mu \end{eqnarray*}$$

and it is easy to see that the saddle point equations (the instanton equations) are given by

$$\begin{eqnarray}&&\dot{u}=b[u]-{\rm{i}}\chi \mu \end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&\dot{\mu }=-{({\nabla }_{u}b[u])}^{T}\mu -{\rm{i}}\lambda {\nabla }_{u}F[u]\delta (t),\end{eqnarray} \tag{ 12 }$$

where we have incorporated the Lagrange multiplier λ at the right-hand side of the equation for μ as a final condition equivalent to

$$\begin{eqnarray}&&\dot{\mu }=-{({\nabla }_{u}b[u])}^{T}\mu ,\qquad \mu (0)=-{\rm{i}}\lambda {\nabla }_{u}F[u(t=0,x)].\end{eqnarray} \tag{ 13 }$$

A solution $(\tilde{u},\tilde{\mu })$ of this set of equations is termed the instanton configuration or instanton. It represents an extremal point of the action functional (9).

We want to make a few remarks on the instanton configuration and the structure of the chosen path integral representation:

• In the limit of applicability of the saddle point approximation, the instanton configuration corresponds to the most probable trajectory connecting the initial conditions to a final configuration which complies with the additional constraint defined by the observable $O[u]$. In the language of quantum mechanics, it corresponds to the classical trajectory obtained for the limit ${\hslash }\to 0$. For stochastic differential equations (SDEs) of the form (1), we similarly need a smallness parameter to justify the saddle point approximation. This might either be the limit of vanishing forcing, $\chi (0)\to 0$, or the limit of extreme events, $| a| \to \infty$ (which corresponds to the limit $| \lambda | \to \infty$, see discussion in sections 2.2 and 5.3).
• It is instructive to apply a change of variables $(u,p)\equiv (u,-{\rm{i}}\mu )$. The response functional (9) can then be written in terms of a Hamiltonian $H[x,p]$,
  $$\begin{eqnarray}&&S[u,p]=\int \left(\langle p,\dot{u}\rangle -H[u,p]\right){\rm{d}}t,\end{eqnarray} \tag{ 14 }$$
  $$\begin{eqnarray}&&H[u,p]=\langle p,b[u]\rangle -\displaystyle \frac{1}{2}\langle p,\chi p\rangle .\end{eqnarray} \tag{ 15 }$$
  The field variable u and the auxiliary field p then play the role of generalized coordinate and momentum for the Hamiltonian system defined by $H[u,p]$. The instanton equations (11), (12) are the corresponding Hamiltonian equations of motion. We remark in particular that the Hamiltonian $H[u,p]$ is a conserved quantity even if the original dynamical system (1) is dissipative.
• In the above setup, note that the choice p = 0, $\dot{u}=b[u]$ is a solution of the equations of motion with vanishing action, corresponding to a deterministic motion of the original dynamical system without any perturbation by noise. Since the action in general is always positive, $S[u,p]\geqslant 0$, this implies that a deterministic trajectory connecting the initial and final point (if it exists) will always be the global minimizer of the action functional.
• Another special solution with H = 0 is the choice $p=-2{\chi }^{-1}b[u]$, which implies $\dot{u}=-b$, i.e. the minimizer follows reversed deterministic trajectories. This choice is only consistent with the auxiliary equation of motion if ${\nabla }_{u}b[u]={({\nabla }_{u}b[u])}^{T}$, as is easily verified by direct comparison between $\dot{p}$ and equation (12). This restriction is only fulfilled, if the drift term $b[u]$ is gradient, $b[u]={\nabla }_{u}V[u]$, which forms an important sub-class of problems. Here, minimum action paths become minimum energy paths.

An alternative approach is given by Wentzell and Freidlin [55, 56]. It has its roots in the application of LDT [57–59] to dynamical systems under random perturbations. In general, we say that a family of random processes ${u}^{\epsilon }(t)$ defined on $t\in [-T,0]$ fulfills a large deviation principle, if

$$\begin{eqnarray}&&P[{u}^{\epsilon }(0)\in \Omega ]\asymp \mathrm{exp}\left(-\displaystyle \frac{1}{\epsilon }\underset{\psi }{\mathrm{inf}}\;{I}_{T}(\psi )\right),\end{eqnarray} \tag{ 16 }$$

where '$\asymp$' is to be understood as the ratio of the logarithms of both sides tending to unity as $\epsilon \to 0$. The left-hand side describes the probability of the random process to 'end up' in some set Ω. On the right-hand side, the infimum is taken over all realizations of the path that fulfill the boundary conditions. The functional I_T is called the rate function, which plays the same role as the action functional in the previous section. The minimizer ${\psi }_{*}(t)$ of the rate function represents the maximum likelihood realization of the random process fulfilling the boundary conditions, and corresponds to the instanton configuration of the MSRJD-formalism. In the context of the Freidlin–Wentzell theory, the rate function can be written down explicitly as

$$\begin{eqnarray}{I}_{T}(\psi )=\left\{\begin{array}{ll}{\displaystyle \int }_{-T}^{0}L(\psi ,\dot{\psi }){\rm{d}}t & {\rm{if}}\;{\rm{the}}\;{\rm{integral}}\;{\rm{exists}},\\ \infty & {\rm{otherwise}},\end{array}\right.\end{eqnarray} \tag{ 17 }$$

where

$$\begin{eqnarray}&&L(x,y)=\underset{p}{\mathrm{sup}}(\langle y,p\rangle -H(x,p))\end{eqnarray} \tag{ 18 }$$

is the Lagrangian of the system with Hamiltonian $H(x,p)$, and $\langle \cdot ,\cdot \rangle$ is the inner product of the considered space.

We are interested in systems of the form (1),

$$\begin{eqnarray}&&{\rm{d}}{u}^{\epsilon }(t)=b[{u}^{\epsilon }(t)]\;{\rm{d}}t+\sigma \sqrt{\epsilon }{\rm{d}}W,\end{eqnarray} \tag{ 19 }$$

with a d-dimensional Brownian increment dW and σ such that $\chi ={\sigma }^{T}\sigma$, where we can write down the Hamiltonian

$$\begin{eqnarray}&&H(u,p)=\langle b[u],p\rangle +\displaystyle \frac{1}{2}\langle p,\chi p\rangle \end{eqnarray} \tag{ 20 }$$

and consequently the Lagrangian

$$\begin{eqnarray}&&L(u,\dot{u})=\displaystyle \frac{1}{2}\langle \dot{u}-b[u],{\chi }^{-1}(\dot{u}-b[u])\rangle .\end{eqnarray} \tag{ 21 }$$

The rate function thus looks like

$$\begin{eqnarray}&&{I}_{T}(u)=\displaystyle \frac{1}{2}{\displaystyle \int }_{-T}^{0}\langle \dot{u}-b[u],{\chi }^{-1}(\dot{u}-b[u])\rangle \;{\rm{d}}t.\end{eqnarray} \tag{ 22 }$$

In a similar manner as before we consider observables at the end point t = 0 for a scalar functional $F[{u}^{\epsilon }(t=0)]$ and want to compute the probability P(a) of the random process evolving such that at the final time $F[{u}^{\epsilon }(t=0)]={a}^{*}$, with ${a}^{*}\in {\mathbb{R}}$. As we will see, this corresponds to considering the large deviations of the moment-generating functions of the random process, $\langle \mathrm{exp}(\frac{\lambda }{\epsilon }F[{u}^{\epsilon }(0)]\rangle$, with $\lambda \in {\mathbb{R}}$. Starting from the large deviation principle (16) we define the set ${\Omega }_{a}=\{\phi | F(\phi )={a}^{*}\}$ of final conditions that comply with the requirement of our observable. Now, defining ${S}^{*}(\lambda )=\epsilon \mathrm{log}\langle \mathrm{exp}(\frac{\lambda }{\epsilon }F[{u}^{\epsilon }(0)]\rangle$, for large λ

$$\begin{eqnarray}&&{S}^{*}(\lambda )=\mathrm{log}{\displaystyle \int }_{{\mathbb{R}}}\mathrm{exp}(\lambda a/\epsilon )P(a){\rm{d}}a\sim \underset{a\in {\mathbb{R}}}{\mathrm{max}}(\lambda a/\epsilon -S(a)/\epsilon ),\end{eqnarray} \tag{ 23 }$$

where $S(a):= -\epsilon \mathrm{log}P(a)=-\epsilon \mathrm{log}P({u}^{\epsilon }(0)\in {\Omega }_{a})$. The last step of (23) makes use of Laplace's method to show that ${S}^{*}(\lambda )$ is the Legendre transform of S(a), and we conclude that $\epsilon {S}^{*}(\lambda )=\lambda {a}^{*}-S({a}^{*})$. On the other hand

$$\begin{eqnarray}&&{I}_{T}({\psi }_{*}^{a})=\underset{\psi }{\mathrm{inf}}\;{I}_{T}(\psi )\sim \displaystyle \frac{\lambda }{\epsilon }{a}^{*}-\displaystyle \frac{1}{\epsilon }{S}^{*}(\lambda )\sim S({a}^{*}),\end{eqnarray} \tag{ 24 }$$

where the infimum is taken over all realizations that fulfill the boundary conditions, ${\psi }_{*}^{a}$ is its minimizer, and of course $F[{\psi }_{*}^{a}]={a}^{*}$. In short, the rate function corresponds to the Legendre transform of the logarithm of moment-generating functions. We obtain the PDF of our observable from the Freidlin–Wentzell action.

To obtain the minimizer ${\psi }_{*}^{a}$, find the solution to the equations of motion (instanton equations) corresponding to (22),

$$\begin{eqnarray}\dot{u} & = & b[u]+\chi p\\ \dot{p} & = & -{({\nabla }_{u}b[u])}^{T}p+\lambda {\nabla }_{u}F[u]\delta (t),\end{eqnarray} \tag{ 25 }$$

where the Lagrange multiplier λ assures that the constraint $F[u]=a$ at t = 0 is satisfied. We immediately see that these equations are equivalent to (11), (12) by setting $\mu =-{\rm{i}}p$.

In the context of theoretical chemistry and the computation of reaction rates, a common scenario is a simplification of system (1), where the drift-term is a gradient, $b[u]=-{\nabla }_{u}V(u)$. Here, minimum action paths can be translated into minimum energy paths, which in turn allows for a number of simplifications of the numerical method, as a lot more is known about the properties of the minimizer: between stable fixed points of the deterministic dynamics, which correspond to the minima of the potential landscape V(u), the transition trajectory has to cross through the saddle point u_* of the potential with the lowest potential value $V({u}_{*})$, when moving between neighboring basins of attraction. This simplified structure of the minimizer allows for a direct estimate of the associated transition probability in terms of the height of the potential barrier and the Hessian at the saddle point, ultimately leading to Kramers' law [11], which is a refinement of Arrhenius' law, known from the problem of chemical reaction rates [128]. Even though in the context of fluid dynamics the systems involved are usually not gradient, some of the mentioned methods are used as a basis for more generalized schemes. Notable methods for finding minimizing trajectories in this setup include the nudged elastic band method [129] and the string method [130, 131].

Returning to the problem of an arbitrary action functional, more general methods have been developed, most notably the MAM [96], which can be seen as a generalization of [132] and variants, for example adaptive MAM [68, 133], parallel MAM [134], and gentlest ascent dynamics [135]. The main idea for an efficient computation is to exploit the structure of the underlying Euler–Lagrange (instanton) equation (25) for efficient computation of the solution. This structure might depend on the particular problem and thus, the computational method often needs to be adapted to the concrete problem. An important step in the development of a method that can be applied in a very general context is the geometric minimum action method (gMAM). The gMAM can be viewed as a generalization of the string method for non-gradient fields. Its starting point is the construction of the quasi-potential that exists even if the drift b does not possess a potential (for details see chapter 5 of [56]). If the drift b does not depend explicitly on time, time can be viewed as a parametrization of the instanton equations. This parametrization can, in principle, be changed to one that is more favorable from a computational point of view. Based on this idea, the gMAM [120, 121, 136] was developed which offers an efficient way to compute minimizers for transitions between an initial and a final state and works even in cases where b does not possess a potential. In particular, one arrives at a modified action functional

$$\begin{eqnarray}&&\hat{S}[u,\dot{u}]={\displaystyle \int }_{0}^{1}\left(\parallel \dot{u}{\parallel }_{\chi }\parallel b[u]{\parallel }_{\chi }-{\langle \dot{u},b[u]\rangle }_{\chi }\right){\rm{d}}s,\end{eqnarray} \tag{ 32 }$$

where the search space is restricted to arclength parametrized trajectories. The metric is given by the correlation matrix χ of the problem via

$$\begin{eqnarray}&&{\langle u,v\rangle }_{\chi }\equiv {\langle u,{\chi }^{-1}v)\rangle }_{{L}^{2}},\end{eqnarray} \tag{ 33 }$$

and the norm $\parallel \cdot {\parallel }_{\chi }$ is induced by this scalar product. The action functional (32) is also termed the geometric action functional, and its minimizers are identical to the infinite time minimizers of the original action, if reparametrized to physical time.

An important class of problems that can be solved efficiently by the gMAM are noise-driven transitions between stable fixed points in the context of meta-stability. These problems are difficult to solve computationally in the original parametrization using time, since it can be shown that it takes infinite time for the minimizer to approach the stable fixed point. To illustrate this fact, in the simplest case, consider a one-dimensional Ornstein–Uhlenbeck process ${\rm{d}}u=-\gamma u{\rm{d}}t+\sqrt{\epsilon }{\rm{d}}W$ for the exit from a stable fix point u = 0, hence the transition from $u({T}_{\mathrm{min}})=0$ at ${T}_{\mathrm{min}}\lt 0$ to $u(0)=a\gt 0$. In this case, the explicit solution of the instanton equations yields

$$\begin{eqnarray}&&u(t)=a{{\rm{e}}}^{\gamma t}\left(\displaystyle \frac{1-{{\rm{e}}}^{2\gamma ({T}_{\mathrm{min}}-t)}}{1-{{\rm{e}}}^{2\gamma {T}_{\mathrm{min}}}}\right),\end{eqnarray} \tag{ 34 }$$

such that the action S depends on ${T}_{\mathrm{min}}$, hence

$$\begin{eqnarray}&&S({T}_{\mathrm{min}})=\displaystyle \frac{\gamma {a}^{2}}{\epsilon (1-{{\rm{e}}}^{2\gamma {T}_{\mathrm{min}}})},\qquad \underset{{T}_{\mathrm{min}}\in (-\infty ,0)}{\mathrm{inf}}S({T}_{\mathrm{min}})=\displaystyle \frac{\gamma {a}^{2}}{\epsilon }.\end{eqnarray} \tag{ 35 }$$

The geometric reparametrization used in both the string method and the gMAM avoids the computational difficulties that arise from the fact that the minimizer requires ${T}_{\mathrm{min}}\to \infty$, which makes them in particular efficient for the problems involving the transition between fix points.

As well known from classical mechanics, the alternative to minimizing the action functional is the solution of the corresponding equations of motion, which are equivalent to the instanton equation (25). This dualism transfers into the algorithmic context: Instead of employing global numerical minimization of the action functional, we can choose to solve a pair of coupled deterministic PDEs instead. In general, it is not obvious which choice is advantageous. In the context of fluid problems, though, we will show in the following that many problems are more easily treated by solving the equations of motion.

At the center of this lies the realization that commonly we are interested in observables of the form of (10), measuring a single degree of freedom in the final field configuration. Most importantly, we do not want to specify the full final condition of our field, or have no access to it a priori. For example when finding dissipative structures in Burgers turbulence, we do not know the exact form of the shock structure that will evolve. While the creation of a steep gradient can still be viewed as an exit problem (interpreting the initial state of the system ${u}_{1}=u({T}_{\mathrm{min}},x)=0$ as the stable fix point for ${T}_{\mathrm{min}}\to -\infty$), we would choose $O[u]=\delta (F[u(x,t=0)]-a)$ as observable with $F[u(x,t=0)]={\partial }_{x}u(x=0,t=0)$, notably without specifying the final state u₂. The choice of the observable restricts u₂ to have a gradient equal to a at $x=0,t=0$, but makes no other assumption about its form. In other words, while it is intuitively plausible that the system will develop a diffusion-limited shock-like structure in order to produce a large negative gradient at t = 0, we cannot compute u₂ without actually obtaining the full solution of the instanton equation (25). This is a quite powerful realization: the only input being a restriction on a single degree of freedom of the final condition, the instanton formalism (if applicable) provides us with the most probable final state which fulfills the constraint, the most probable evolution in time from a given initial condition into this state, and the force that was necessary to achieve this evolution.

As a direct consequence, problems of this type are difficult to frame in the context of minimization of the action functional, as the final condition is not fixed, but only subject to a constraint. In contrast to this they are readily solved by considering the instanton equation (25)

$$\begin{eqnarray*}\dot{u} & = & b[u]+\chi p\\ \dot{p} & = & -{({\nabla }_{u}b[u])}^{T}p+\lambda {\nabla }_{u}F[u]\delta (t)\end{eqnarray*}$$

as the term $\lambda {\nabla }_{u}F[u]\delta (t)$ defines an initial condition $p(x,t=0)=-\lambda {\nabla }_{u}F[u]$ for the auxiliary equation, propagated backwards in time. The problem is therefore reduced to solving two mutually dependent deterministic PDEs with known initial conditions. In what follows, we explain in detail how to implement this formalism in the case of Burgers shocks, noting that the scheme is similarly applicable to different problems in fluid dynamics of the same functional form.

In 2001, Chernykh and Stepanov published a seminal paper describing an algorithm suited for such boundary conditions [104]. Their idea was to solve (25) iteratively, demonstrating that the iterative scheme converges to the instanton given by the solution of (25). The original motivation of Chernykh and Stepanov was to confirm analytical predictions regarding the scaling exponent of $3/2$ of the left tail of the velocity gradient, when diffusion counteracts the steepening of the shock profile. For the case of extreme negative velocity gradients, they chose $O[u]=\delta (F[u(x,t=0)]-a)$ with $F[u(x,t=0)]={\partial }_{x}u(x=0,t=0)$ and $a\ll 0$. This choice leads to the boundary condition $p(x,t=0)=-\lambda {\delta }^{\prime }(x)$. In order to iteratively solve the equation, the basic scheme suggested by Chernykh and Stepanov was to view the instanton equations (25) from a dynamical point of view on the time interval $(-\infty ,0]$. For a specified velocity field, for example the field ${u}^{(k)}(x,t)$ of the kth iteration, we can solve the auxiliary equation backwards in time up to a time ${T}_{\mathrm{min}}\ll 0$. From our previous considerations of the simple Ornstein–Uhlenbeck case, it is clear that one would like to choose ${T}_{\mathrm{min}}$ as small as possible in order to approximate the minimizer (which is obtained in the limit ${T}_{\mathrm{min}}\to -\infty$). Once the next iteration ${p}^{(k+1)}(x,t)$ is found on $[{T}_{\mathrm{min}},0]$, this solution is used in the equation for u, to find the next iterate ${u}^{(k+1)}(x,t)$. Note that the convolution $\chi p$ represents the forcing of the system, while the rest of the u-equation resembles the deterministic part of the original Burgers equation (26). Therefore, having found the fixed point of our iteration, we can recover the optimal forcing by evaluating $\chi p$. To start off the iteration, an initial guess ${u}^{0}(x,t)$ is required, e.g. ${u}^{0}(x,t)=0$. The iteration is continued until reaching a fixed point. The following figure illustrates the iteration scheme:

Note that one needs to specify the value of λ in this algorithm as it is present in the term $p(t=0,x)=-\lambda {\delta }^{\prime }(x)$. Usually, it is not known a priori how λ relates to the parameter a in the observable $O[u]=\delta (F[u(x,t=0)]-a)$. One can view λ as the Lagrange multiplier that ensures that at time t = 0 we are observing indeed $F[u(x,t=0)]=a$. In the above scheme, however, one can measure a after each iteration by computing $F[{u}^{(k)}(x,t=0)]$. Once the differences in the values of a are below a chosen threshold, the desired approximation of the instanton has been computed. An example implementation is given in the appendix in listing 1 for the case of a nonlinear Ornstein–Uhlenbeck process.

Although this scheme appeals by its simplicity and generality, there are several computational challenges associated with it. For the Burgers case, this was already noted by Chernykh and Stepanov. Some of them are so severe that they make a direct application of the scheme almost impossible. In the following, we discuss two particular challenges: (a) the numerical instability of the fixed point and (b) convergence issues related to the fact that one needs to take ${T}_{\mathrm{min}}\to -\infty$ in order to compute the minimizer.

Regarding (a), we confirmed in our simulations results from Chernykh and Stepanov: for the Burgers case, for large λ (corresponding to large negative gradients a), taking directly the updated ${u}^{(k+1)}$ as the new advection field in the equation for p can lead to instabilities. These instabilities can be suppressed by a softer update in the spirit of a Gauss–Seidel iteration: Assume that, in the above scheme, we have computed the new velocity field u on $({T}_{\mathrm{min}},0]$ and let us call this field ${\tilde{u}}^{(k+1)}$. Then we can take as an update a weighted average of the newly computed field ${\tilde{u}}^{(k+1)}$ and the old field ${\tilde{u}}^{(k)}$ using a weight parameter $\alpha \in [0,1)$

$$\begin{eqnarray}&&{u}^{(k+1)}(x,t)=\alpha {u}^{(k)}(x,t)+(1-\alpha ){\tilde{u}}^{(k)}(x,t).\end{eqnarray} \tag{ 36 }$$

During the iterative procedure, we monitor the numerical value of the action at each iteration step to ensure that the update yields a minimizer. Note that one could have also used a gradient based method to descent into the minimizer. However the choice above enables us to re-use existing implementations of the deterministic dynamics.

Regarding (b), the difficulties associated with ${T}_{\mathrm{min}}\to -\infty$ require a more extensive modification of the algorithm. One approach is to borrow an idea from geometric methods, in particular the string method and the minimum action method (gMAM), in order to find a suitable reparametrization of the evolution variable [122]. Hamilton's equations, which describe the trajectory of the minimizer of the Freidlin–Wentzell functional, are given by

$$\begin{eqnarray}&&\dot{u}=\displaystyle \frac{\partial H}{\partial p},\qquad \dot{p}=-\displaystyle \frac{\partial H}{\partial u}.\end{eqnarray} \tag{ 37 }$$

If H does not depend explicitly on time, we can reparametrize these equations using $t=t(s)$ (assuming $t\prime (s)\gt 0$) to obtain

$$\begin{eqnarray}&&u\prime =t\prime (s)\displaystyle \frac{\partial H}{\partial p},\qquad p\prime =-t\prime (s)\displaystyle \frac{\partial H}{\partial u}.\end{eqnarray} \tag{ 38 }$$

We know from the gMAM that a suitable parametrization to circumvent the computational difficulties associated with ${T}_{\mathrm{min}}\to -\infty$ is to choose t(s) in such a way that $\parallel u\prime \parallel =C$. The norm $\parallel \cdot \parallel$ is simply the Euclidean length for δ-correlated noise. For a correlation function χ, the induced norm $\parallel \cdot {\parallel }_{\chi }$ is given by $\parallel f{\parallel }_{\chi }^{2}=\langle f,{\chi }^{-1}f\rangle$ and we choose t(s) such that $\parallel u\prime {\parallel }_{\chi }=C$. In analogy, we call this arclength-parametrization of Hamilton's equations for the minimizer. For the Burgers case, working directly with the correlation function given in (27) can lead to instabilities for high frequencies. Therefore, one can truncate the correlation function in Fourier space, e.g. work with

$$\begin{eqnarray}&&\hat{\chi }(\omega )={\omega }^{2}{{\rm{e}}}^{-{\omega }^{2}/2}{\mathcal{H}}({\omega }_{c}-| \omega | ),\end{eqnarray} \tag{ 39 }$$

where ${\mathcal{H}}$ denotes the Heaviside function. Then, the correlation function χ can be inverted on its support in Fourier domain and we can implement the induced norm $\parallel \cdot {\parallel }_{\chi }$ via

where $\hat{f}$ is the Fourier transform of f and ${{\mathcal{F}}}^{-1}$ denotes the inverse Fourier transform. Then, for Burgers equation, the arclength-parametrized Hamilton equation (38) are written as

$$\begin{eqnarray}&&{u}_{s}=\displaystyle \frac{\parallel {u}_{s}{\parallel }_{\chi }}{\parallel b[u]{\parallel }_{\chi }}(b[u]+\chi p),\qquad {p}_{s}=-\displaystyle \frac{\parallel {u}_{s}{\parallel }_{\chi }}{\parallel b[u]{\parallel }_{\chi }}({{up}}_{x}+\nu {p}_{{xx}}).\end{eqnarray} \tag{ 41 }$$

Note that it follows from the above Hamilton's equations that we have at the minimizer

$$\begin{eqnarray}&&\parallel b[u]{\parallel }_{\chi }=\parallel b[u]+\chi p{\parallel }_{\chi }.\end{eqnarray} \tag{ 42 }$$

This equation can be used in order to regularize the term $\parallel b[u]{\parallel }_{\chi }$ in numerical simulations if necessary. The arclength parametrized equations of motion (41) can similarly be obtained by considering the variation of the geometric action (32).

It is noteworthy that Chernykh and Stepanov were aware of the difficulties related to taking ${T}_{\mathrm{min}}\to -\infty$. They addressed this difficulty in their paper in two ways. First, they replaced the initial condition for u at ${T}_{\mathrm{min}}\to -\infty$ with a modified initial condition found from the linearization of the instanton equations around the fixed point. Second, they used a self-similar transform related to the kernel of the heat equation which resulted in an exponential rescaling of the time. While this approach proved to lead to efficient numerics for the Burgers case, the above chosen geometric reformulation leads to similar results but furthermore generalizes to other equations in a straightforward way.

We remark that for known initial and final conditions of the field variable, the method outlined above is not easily applied. In this case there is no restriction on the auxiliary variable, while the field variable has to hit a particular point u₂ starting from u₁. To solve the instanton equations, one can then only rely on shooting methods [119], which is nearly hopeless in such high-dimensional systems. Therefore, the direct minimization techniques discussed in section 5.2 are in general preferred.

In order to show the simplicity of the numerical filtering, let us consider a simple one-dimensional Ornstein–Uhlenbeck process, hence

$$\begin{eqnarray}&&{\rm{d}}u=b(u){\rm{d}}t+\sigma {\rm{d}}W,\qquad b(u)=-\gamma u,\qquad \gamma \gt 0.\end{eqnarray} \tag{ 43 }$$

We consider the exit from a stable fixed point at ${\mathrm{lim}}_{t\to -\infty }\;u(t)=0$ such that $u(0)=a$. A simple analytic solution of the instanton equations (11), (12) shows that the minimizer of the Freidlin–Wentzell action is given by

$$\begin{eqnarray}&&u(t)=a{{\rm{e}}}^{\gamma t}.\end{eqnarray} \tag{ 44 }$$

How can we obtain a numerical approximation of this minimizer via DNS? One way is to generate an ensemble of realizations and filter the paths that end in a small neighborhood of the target value $u(0)=a$, for example given by the interval $a-\epsilon \lt u(0)\lt a+\epsilon$. Since the probability distribution will have its maximum around the minimizer, and the fluctuations around the minimizer are (at the leading order) Gaussian, we can expect the mean of the filtered sample paths to approximate the minimizer. Figure 2 shows the result of 10⁷ realizations. The figure on the left presents the evolution of the mean of the filtered paths. The figure on the right shows the result of the simulations in a slightly different way: here, a rectangle $[-T,0]\times [-1,5]$ was divided up in cells and the number of hits per cell was computed. From the Freidlin–Wentzell theory, we know that the paths will cluster around the minimizer in the weak-noise limit of $\sigma \to 0$. An example implementation of this algorithm is given in listings 1 and 2 in the appendix.

Zoom In Zoom Out Reset image size

In this section we will comment briefly on the extension of the filtering to jump-processes, in particular related to SDEs driven by symmetric alpha-stable processes. In the context of fluid mechanics, one might think of a system that is driven by an already non-Gaussian noise, such as the ocean driven by wind. In general, such SDEs have attracted increasing interest over the last decade, not only in the context of turbulence, but also in biology and particularly in finance [137]. An area of active research regarding fluids is the question whether Levy walks can be used to model pair dispersion in turbulent flows related to the Richardson law [138, 139].

In the present work, we restrict ourselves to the simplest basic approach to review how one can use path integrals in the context of Levy-driven SDEs. Levy processes are also called fractional noise processes as their characteristic function (see below) involves a non-integer exponent α. It is well known that, using the MSRJD approach, for such Levy processes, a generalization of the Wiener-path integral can be formulated [140]. The starting point of our discussion is the stochastic differential equation (19) where the Brownian increment dW is now replaced by the Levy-increment dL,

$$\begin{eqnarray}&&{\rm{d}}u=b[u]\;{\rm{d}}t+\sigma {\rm{d}}L.\end{eqnarray} \tag{ 45 }$$

For simplicity, let us discuss the above equation in a one-dimensional context. For symmetric alpha-stable processes, the characteristic function of L_t is given by

$$\begin{eqnarray}&&{\mathbb{E}}({{\rm{e}}}^{{\rm{i}}{{sL}}_{t}})={{\rm{e}}}^{-t| s{| }^{\alpha }}.\end{eqnarray} \tag{ 46 }$$

In the following we assume $1\lt \alpha \lt 2$. The characteristic function lacks smoothness around the origin, which gives rise to fat tails. Indeed, it can be shown [141] that (set t = 1) we have

$$\begin{eqnarray}&&P(u)=\displaystyle \frac{1}{2\pi }{\displaystyle \int }_{{\mathbb{R}}}{{\rm{e}}}^{-{\rm{i}}{su}}\varphi (s)\;{\rm{d}}s=\displaystyle \frac{1}{2\pi }{\displaystyle \int }_{{\mathbb{R}}}{{\rm{e}}}^{-{\rm{i}}{su}-| s{| }^{\alpha }}\;{\rm{d}}s\sim \displaystyle \frac{1}{{| u| }^{\alpha +1}},\end{eqnarray} \tag{ 47 }$$

as $| u| \to \infty$. While the theory of SDEs driven by Levy processes is much more involved compared to SDEs driven by Brownian motion, it is relatively simple to adapt the filtering technique to such equations. Equation (45) can be discretized in the following form:

$$\begin{eqnarray}&&{u}_{n+1}={u}_{n}+b({u}_{n}){\rm{d}}t+{r}_{n}{\rm{d}}{t}^{1/\alpha },\end{eqnarray} \tag{ 48 }$$

where the random number r_n is generated, for example, by a Box–Muller-like algorithm: draw v_n from a uniform distribution on $(-\pi /2,\pi /2)$ and draw w_n from an exponential distribution with mean 1. Then compute

$$\begin{eqnarray}&&{r}_{n}=\displaystyle \frac{\mathrm{sin}(\alpha {v}_{n})}{{(\mathrm{cos}({v}_{n}))}^{1/\alpha }}{\left(\displaystyle \frac{\mathrm{cos}({v}_{n}-\alpha {v}_{n})}{{w}_{n}}\right)}^{(1-\alpha )/\alpha }.\end{eqnarray} \tag{ 49 }$$

This simulation technique can be extended to multi-dimensional processes in a straight forward way, using a Cholesky decomposition in order to account for a specific correlation function. Even in simple cases, however, the behavior of the solution differs fundamentally from the diffusive case. In order to illustrate these differences, as a simple example, we will discuss a Levy-driven Ornstein–Uhlenbeck process. Note that, for such processes, the corresponding path-integrals can be explicitly evaluated [142]. For our numerical simulations, we have the Ornstein–Uhlenbeck process u follow the SDE given by

$$\begin{eqnarray}&&{\rm{d}}u=b(u){\rm{d}}t+\sigma {\rm{d}}L,\qquad b(u)=-\gamma u.\end{eqnarray} \tag{ 50 }$$

Here, again, we discuss the exit from a stable fixed point at negative infinity such that $u(t=0)=a$. For the Levy case, consider as an example $\alpha =1.8$. Applying the filtering technique, we see that the filtered mean $\langle u\rangle$ of the process for the Levy case follows a similar form given by

$$\begin{eqnarray}&&\langle u\rangle \approx a\;{{\rm{e}}}^{\gamma (\alpha -1)t}.\end{eqnarray} \tag{ 51 }$$

figure 3 shows again the result of 10⁷ simulations. As for the diffusive case, the figure on the left presents the evolution of the mean of the filtered paths and the figure on the right shows the histogram for a rectangle $[-T,0]\times [-1,5]$. From here we see that, although the mean $\langle u\rangle$ appears to yield a smooth curve, the typical exit path exhibits an entirely different evolution: It deviates only very little from the stable fixed point at u = 0 until a (random) time t_r, only to then perform a jump of height $a\;{{\rm{e}}}^{\gamma {t}_{r}}$, such that the deterministic drift b will take it to the point a at the time t = 0. In the diffusive case, as we have seen above, the filtered paths were distributed around the mean paths which coincided with the minimizer of the Freidlin–Wentzell action. An example implementation for the case of Levy noise given in listings 2 and 4 in the appendix.

Zoom In Zoom Out Reset image size

We now describe how to apply filtering in the context of the stochastic Burgers equation. Here, the situation is more expensive than for the Ornstein–Uhlenbeck process since we are dealing with many more dimensions: for each realization we need to solve equation (26) with the appropriate right-hand side. The basic idea of the filtering remains the same: we generate an enormous number of stochastic simulations and filter the realizations that produce the rare event that our observable measures. In the case under consideration here, this observable is the slope of the velocity field and we are in particular interested in the probability of large negative gradients.

There are a variety of numerical methods available to generate the noise term on the right-hand side of equation (26) in DNS, often the preferred methods are Fourier-based [110, 143]. The idea is that the main results, for instance in our case the scaling of the probability distribution of the velocity gradient, should show universal behavior and not depend too much on the particular choice of the noise term as long as the noise is sufficiently weak and does not dominate the system. In the following, however, we discuss the results of a particular numerical experiment aimed at showing the relevance of instantons [111]. Here, we are not only interested in the scaling of the probability distributions but in the detailed structure of the instantons, in particular their time history. For this purpose, it is desirable to generate noise that imitates the fluctuations assumed in the instanton analysis as closely as possible. This can be done as follows:

• (i)  
  Draw a vector r of appropriately scaled normally distributed random numbers. The size of the vector corresponds to the discretization in x.
• (ii)  
  Multiply this vector by a matrix A resulting from the Cholesky decomposition of the (discretized) correlation matrix C. Note that naive discretization of χ is likely to lead to a $\tilde{C}$ that, due to finite machine-precision, is not positive-semidefinite. This is due to the fact that, for a large number of discretization points, the rows of $\tilde{C}$ are almost linearly dependent. In this case it is possible to use the algorithm introduced in [144] in order to obtain a matrix C that is positive-semidefinite and sufficiently close to $\tilde{C}$.

When implementing the DNS, the available computational hardware needs to be used in an efficient way. In our work from [111], we used a combination of accelerator graphics cards (CUDA) and multiple processors connected via MPI: for the required resolution, the size of one simulation was small enough to fit on a graphics card such that a single realization could be performed directly on the card. In addition, other numerical procedures necessary to extract the instanton (see following section) could be performed as well on the graphics card. MPI was used in order to average over the different CUDA realizations. This made the actual implementation of the algorithm fairly simple but efficient: Since the realizations are stochastically independent, we obtain a linear scaling with the number of graphics cards.

This approach, in principle, is applicable to any generic SPDE with one spatial dimension. Already for two-dimensions (see section 7), the memory requirements grow significantly, so that in these cases more sophisticated numerical algorithms have to be employed.

In order to obtain the instanton from an ensemble of DNS, the following filtering can be applied: prescribe a small interval around the desired value of the observable and keep the relations for which the value of the observable falls into this interval. This is similar to the procedure previously discussed and illustrated for the Ornstein–Uhlenbeck process. In the concrete case for Burgers equation, we prescribe the interval around a specified velocity gradient at x = 0 at t = 0. We assume that the realizations start at ${T}_{\mathrm{min}}$ from a zero initial condition and are integrated to the final time t = 0. For computational efficiency, our filtering algorithm does not only look at the point x = 0 but makes use of spatial translation invariance. This requires shifting of both the velocity and the forcing field. In the simple case of the Ornstein–Uhlenbeck process discussed above, such shifting was unnecessary since we were working in only one dimension. In the case of a stochastic partial differential equation like the stochastic Burgers equation, however, the rare event (here the large negative gradient) can occur at any place in the profile. The above described procedure for searching the maximum gradient and shifting the fields allows for detecting many more rare events.

The averaging procedure now consists of taking the average of all those shifted fields ${u}_{\mathrm{shifted}}(t,x)$. Note that, from the same simulations, we obtain also a corresponding force field ${\eta }_{\mathrm{shifted}}(t,x)$ and, for both ensembles, we can compute the ensemble average $\langle {u}_{\mathrm{shifted}}(t,x)\rangle$ and $\langle {\eta }_{\mathrm{shifted}}(t,x)\rangle$ in space and time. Due to the δ-correlation of the driving noise (in time), we expect convergence only after many realizations. This filtering approach can be interpreted as the numerical counterpart to the MSRJD formulation of the path integral taking into account the chosen observable $O(u)=\langle \delta ({u}_{x}(0,0)=a)\rangle$ as the conditional mean of the observable (where the conditioning is happening in our case due to the constraint of the velocity gradient at the end point of the evolution). For sufficiently strong gradients, corresponding to sufficiently rare events, we therefore expect that the ensemble averages $\langle {u}_{\mathrm{shifted}}(t,x)\rangle$ and $\langle {\eta }_{\mathrm{shifted}}(t,x)\rangle$ will be close to the instanton solution of (11), (12). The corresponding optimal force $\langle {\eta }_{\mathrm{shifted}}(t,x)\rangle$ is related to the instanton solution of (11) via $\langle {\eta }_{\mathrm{shifted}}(t,x)\rangle =-{\rm{i}}\chi \mu$.

The first important choice in the direct numerical simulation is setting the initial time ${T}_{\mathrm{min}}$. Clearly, from an analytical point of view, one requires ${T}_{\mathrm{min}}\to -\infty$ and, therefore $-{T}_{\mathrm{min}}$ should be as large as possible in order to let the instanton develop. The maximum $\Delta t$ for the resolution in time (usually the choice of $\Delta t$ is enforced by a stability condition of the numerical scheme and related to the spatial discretization) sets the number N_t of discretization points in time. We found that, in fact, the results are very sensitive to the choice of ${T}_{\mathrm{min}}$ if $-{T}_{\mathrm{min}}$ is too small. Therefore, it is recommended to perform several runs with a variety of values for the parameter ${T}_{\mathrm{min}}$ in order to be sure that $-{T}_{\mathrm{min}}$ is sufficiently large. In our simulations, we found that a decent value of ${T}_{\mathrm{min}}$ is obtained by comparing ${T}_{\mathrm{min}}$ to the integral time T_L and we chose $| {T}_{\mathrm{min}}| \gt 10{T}_{L}$. The averages were obtained from about 10⁷ stochastically independent realizations.

Using this data set of simulations at different forcing strengths (or equivalently, different Reynolds numbers), we can compare the applicability of the instanton approximation for events of extreme negative gradients in turbulent Burgers flows. In figures 4 (left) and 5 (left) this comparison is demonstrated for low and high values of the forcing prefactor, respectively. The instanton prediction for the left velocity gradient PDF tail, $\mathrm{exp}(-{S}_{\mathrm{inst}})$, agrees over the whole range of gradients in the low forcing limit, where the instanton approximation becomes asymptotically exact. As expected, for a higher value of forcing, the instanton prediction captures only the scaling of the tail for events of extreme negative gradients. Figures 4 (right) and 5 (right) on the other hand demonstrate the agreement of the final configuration ${u}_{\mathrm{inst}}(x,t=0)$ of the instanton against the filtered field $\langle {u}_{\mathrm{shifted}}(t,x)\rangle$ for low and high values of forcing, respectively. While in the limit of low forcing the average shock event resembles the instanton prediction almost exactly over the whole range of gradients, for higher values of the forcing the agreement is only visible in the far left tail of the velocity gradient PDF.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As expected, the filtering approach does not only extract the final configuration at time t = 0 but also the entire time history of the evolution of the instanton and of the optimal force which is compared in figure 6 to the solution of the instanton equations. The agreement is remarkable.

Zoom In Zoom Out Reset image size

The instanton equations for the field variable (11) and the auxiliary variable (12) are a mixed initial/final value problem. While the field variable is integrated forward in time, starting from an initial condition u₁ at $t=-T$, the auxiliary variable is integrated backwards in time from the final condition p₂ at t = 0. Both equations mutually depend on each other. A similar form of a mixed initial/final value problem is encountered in a different area of fluid dynamics, namely specific questions regarding a passively transported scalar in a flow: analyzing where the measured density at a given point originates from. A method for solving the resulting coupled equations—one for the evolving fluid, forward in time, and one for the passive scalar, integrated from its given final density backwards in time—was described in [145]. To elucidate this method, consider a simpler system of equations on the time interval $t\in [-T,0]$ of the form

$$\begin{eqnarray}&&{u}_{t}=f(u,t),\qquad u(-T)={u}_{1}\end{eqnarray} \tag{ 52 }$$

$$\begin{eqnarray}&&{p}_{t}=g(u,p,t),\qquad p(0)={p}_{2},\end{eqnarray} \tag{ 53 }$$

as described in [145] and which is closely related to the situation in control theory [151]. There are two opposing ways of solving this system numerically: (i) solving (52) forward in time starting from u₁ at $t=-T$ and saving the complete evolution of u(t) along the way, then subsequently solving (53) backwards in time, using the stored solution u(t) to evaluate $g(u,p,t)$, and (ii) solve (53) backwards in time, starting at t = 0 from p₂, while computing u(t) at each timestep by integrating equation (52) forward in time from u₁ at $t=-T$. Method (i) is ${\mathcal{O}}({N}_{t})$ in both memory and computing time, which is the worst case in memory, while method (ii) scales only ${\mathcal{O}}(1)$ in memory (which is optimal), while being ${\mathcal{O}}({N}_{t}^{2})$ in computing time. Note that variant (ii) is only possible due to the fact that equation (52) is independent of the 'auxiliary' field p. In the case of mutual dependence, such as the system of instanton equations, this variant does not apply.

These two building blocks can now be employed to obtain an efficient compromise—also known as checkpointing in control theory [152,153]—between the two scalings and balancing memory efficiency and computational cost. First, split the interval $[-T,0]$ into k sub-intervals. After integrating (52) once and storing the result at the beginning of each sub-interval, we reduce the original problem to a similar problem on a shorter domain of size ${N}_{t}/k$. Recursively, for each of these sub-intervals, we can break down the problem by splitting the domain, until we reach an interval length of only a single integration step. At this point we can carry out the integration. Inspired by similar principles in multi-grid algorithms or the fast Fourier-transform, a natural choice is k = 2. In this case, the memory requirement scales as ${\mathcal{O}}(\mathrm{log}{N}_{t})$ and computing time as ${\mathcal{O}}({N}_{t}\;\mathrm{log}{N}_{t})$. A schematic depiction of the algorithm for k = 2 and N_t = 16 is shown in figure 7. For the initial solution of the field equation, the field is stored at the intermediate timesteps $i\in \{8,12,14,15,16\}$, of which the least three can be used immediately for the backwards propagating auxiliary equation. Whenever a timestep is encountered for which the field configuration is not stored, such as i = 13, it is propagated forward from the last known position (in this case i = 12), while storing intermediate values recursively in the same fashion.

Zoom In Zoom Out Reset image size

The algorithm described above cannot be used to solve the instanton equations (11), (12) without modification, because in contrast to the model problem (52), (53), both equations are mutually dependent and neither the u- nor the p-equation can be solved without referring to the other field. As a consequence, at least one of the fields has to be stored for each time-step. Note though that in equation (11) the auxiliary field only acts on u through a convolution with the forcing correlation. Since in fluid dynamical applications the forcing is usually restricted to only a few active modes, $\chi p$ is very compact to store. In other words, the optimal forcing $\chi p$ is the only field that needs to be kept for every timestep to be able to integrate the complete velocity field and auxiliary field of the instanton configuration. We remark that this approach cannot recover the ${\mathcal{O}}(\mathrm{log}{N}_{t})$ scaling of the original algorithm. On the other hand, since the number of active modes of the forcing are constant, the memory cost of the auxiliary field becomes independent of the number of degrees of freedom in space N_x. In contrast to this, the expensive storage of the field variable scales logarithmically in time. This interplay of savings allows for the drastic savings in memory necessary to compute higher-dimensional instantons.

As an example application for the numerical computation of higher-dimensional instantons in fluid applications, we will consider the Burgers equation in 2D [146]

$$\begin{eqnarray}&&{\partial }_{t}u+u\cdot \nabla u-\nu \Delta u=f,\end{eqnarray} \tag{ 54 }$$

for a forcing $f(x,t)$ correlated white in time and with a finite correlation length L in space, active only on large-scale modes (i.e. truncating the correlation function in Fourier space above a mode ${\omega }_{c}$, as in (39)). In contrast to the 2D incompressible Navier–Stokes equation, for which was shown in [27] that the naive instanton will be trivial, the 2D Burgers equation exhibits a direct cascade. Furthermore, equation (54) preserves irrotationality of the flow under irrotational forcing, which is the scenario we will restrict ourselves to in this section.

As generalization of the observable taken in 1D, i.e. $F[u]={\partial }_{x}u(x=0,t=0)$, several choices are possible: taking $F[u]=\nabla \cdot u$ conditions on events exhibiting extreme velocity divergence at a single point, which corresponds to 'explosion' and 'implosion' events. Instead of this, we will focus on the case

$$\begin{eqnarray}&&F[u]={\partial }_{x}{u}_{x}(0,0),\end{eqnarray} \tag{ 55 }$$

which selects events of high velocity gradient in a single direction. This observable is closely related to the energy dissipation, $\nu {| \nabla u| }^{2}$, but additionally breaks rotational symmetry. We will identify the instanton solution with shock structures known from 2D Burgers turbulence. The instanton equations corresponding to equation (54) read

$$\begin{eqnarray}&&{\partial }_{t}u+u\cdot \nabla u-\nu \Delta u=\chi p\end{eqnarray} \tag{ 56 }$$

$$\begin{eqnarray}&&{\partial }_{t}p+u\cdot \nabla p-(p\times \nabla ){u}^{\perp }+\nu \Delta p=0,\end{eqnarray} \tag{ 57 }$$

with ${u}^{\perp }=(-{u}_{y},{u}_{x})$ and the convolution ${(\chi p)}_{i}={\sum }_{j}{\chi }_{{ij}}{p}_{j}$. Figure 8 (left) depicts the solution of these equations with the numerical scheme laid out above. Shown is a surface-plot of the x-component of the velocity field for the instanton around the origin at t = 0. As expected, the size of the shock structure across the shock is determined by the viscosity $\nu \ll 1$, while its length along the shock is given by the forcing correlation length L = 1.

Zoom In Zoom Out Reset image size

This example demonstrates the drastic savings in memory that are possible by employing the recursive technique: The highest resolution achieved for the 2D Burgers problem is ${N}_{x}=1024\times 1024$, N_t = 8192, with a total memory usage of 577MB. This fits completely into the memory of a single graphics card on which the computation was undertaken. In contrast, the extrapolated memory usage of the un-optimized algorithm would be about 300 times higher. Note also that the computational overhead of the optimization, due to logarithmic scaling, is slightly lower than a factor 3 in the total computation time. Figure 8 (right) shows a comparison of the memory cost for a lower resolution setup of ${N}_{x}=256\times 256$, N_t = 2048 in order to fit the unoptimized case on the machine: Both the recursive optimization and the projection method amount to a saving of roughly one half, the combination of both leads to memory savings of a factor 20.

As last example of higher-dimensional fluid systems to be treated with the instanton formalism we will consider the 3D incompressible Navier–Stokes equations

$$\begin{eqnarray}&&{\partial }_{t}u+u\cdot \nabla u-\nabla p-\nu \Delta u=f,\quad \nabla \cdot u=0.\end{eqnarray} \tag{ 58 }$$

This equation of course lies at the very center of turbulence research. In spirit of the program laid out above there is hope to link the dissipative structures of turbulent flows to the configuration obtained by the instanton formalism. This in turn might elucidate not only the form and emergence of such structures, but serve to answer questions about the intermittent nature of turbulent flows. Yet, it is far from obvious that the naive approach of formulating the corresponding MSRJD-action (9) for the Navier–Stokes equation (58) and then computing its minimizing trajectory (which in itself is a considerable amount of work) leads to any meaningful results: it is not clear which observable to consider, a choice which can change the structure of the instanton trajectory completely. One might be tempted to look at events of high gradients (like in Burgers) or high energy dissipation

$$\begin{eqnarray}&&F[u]=\nu | \nabla u(x=0,t=0){| }^{2}\end{eqnarray} \tag{ 59 }$$

to obtain the dissipative structures of turbulence. Yet, these events might be well into the dissipative range and thus not capture the multi-fractal nature of the turbulent cascade in the inertial range.

On the other hand, one can obviously commence on the path laid out in sections 5 and 6 for the Burgers equation: comparing the events dominating the extreme tails of the PDF of a chosen observable to the structures obtained by solving the corresponding instanton equations. This way, one gets access to the tail scaling of PDFs of turbulent quantities that are much harder to obtain by classical methods. In order to demonstrate that such an approach is feasible, we show preliminary results of a numerical computation of a 3D Navier–Stokes instanton in figure 9. As observable

$$\begin{eqnarray}&&F[u]=\nabla \times u(x=0,t=0)\end{eqnarray} \tag{ 60 }$$

was chosen, which amounts to events of extreme vorticity in the origin at t = 0. The resulting structure in quality resembles the well-known vorticity filaments observed in developed 3D Navier–Stokes turbulence. For the filtered fields, the sample size is $\approx {10}^{3}$.

Zoom In Zoom Out Reset image size

Structures at the final time t = 0 are very similar to those first obtained by Novikov [147, 148] and analyzed in detail by Wilczek [149]. They consider the conditionally averaged vorticity field

$$\begin{eqnarray}&&{\omega }_{\mathrm{filter}}=\langle \omega (x,t)| \omega (x=0,t=0)=a\rangle ,\end{eqnarray} \tag{ 61 }$$

i.e. the vorticity field in space and time, conditioned on the fact that it will attain the vorticity $\omega (0,0)=a$ in the course of its evolution. This is of course the same as the 'filtering'-approach laid out above for Burgers shock structures. It therefore makes sense to ask whether the instanton resembles the filtered structure for events of extreme vorticity. Preliminary results of this are shown in figure 10. Due to the symmetry of the observable, the instanton solution observes rotational symmetry around the axis prescribed by the vorticity in the origin. The filtered velocity field is obtained by averaging fields of the same vorticity $\omega (x,t)=a$, after translating and rotating them into the origin. Depicted is a direct comparison of all three components of the vorticity field in cylindrical coordinates $(r,\theta ,z)$ between the instanton configuration and the filtered vorticity field of a turbulent DNS of the 3D Navier–Stokes equation, both with ${N}_{x}={128}^{3}$ and N_t = 128 for the instanton. The instanton field reproduces the conditionally averaged one remarkably well. In particular, features like the spreading of the vortex and the appearance of a swirl for larger radii (also visible in figure 9) is well-reproduced. These results serve as evidence that the techniques laid out above are powerful enough to compute instantons even for the 3D setup.

Zoom In Zoom Out Reset image size

In how far this Navier–Stokes vorticity instanton can be analytically captured by the approximation derived in [150] remains a future task.