3.1. Local governing equations

For the purposes of emphasising a probabilistic perspective that focuses on control volumes, we regard the spatial domain $\varOmega$ as a sample space, $\boldsymbol{X}:\varOmega \rightarrow \mathbb{R}^{d}$ as coordinate functions and $\boldsymbol{Y}_{\!t}:\varOmega \rightarrow \mathbb{R}^{n}$ as random variables, as illustrated in figure 3. Let the Eulerian evolution of $\boldsymbol{Y}_{\!t}$ at a given point in the domain be determined by the stochastic differential equation

(3.1) \begin{equation} \,\textrm{d} \boldsymbol{Y}_{\!t}=(\boldsymbol{I}_{\!t}-\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{J}_{\!t})\,\textrm{d} t+\boldsymbol{\sigma }\,\textrm{d} \boldsymbol{W}_{\!t}, \end{equation}

where $\boldsymbol{I}_{\!t}\in \mathbb{R}^{n}$ represents internal sources/sinks, $\boldsymbol{J}_{\!t}\in \mathbb{R}^{d\times n}$ is a matrix of fluxes, $\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{J}_{\!t}\in \mathbb{R}^{n}$ is the divergence of the flux of each component of $\boldsymbol{Y}_{\!t}$ and $\boldsymbol{W}_{\!t}\in \mathbb{R}^{m}$ is a vector of independent Wiener processes. In general, the matrix coefficient $\boldsymbol{\sigma }$ could depend on both $\boldsymbol{X}$ and $\boldsymbol{Y}_{\!t}$ , such that $\boldsymbol{\sigma }(\boldsymbol{X},\boldsymbol{Y}_{\!t})\in \mathbb{R}^{n\times m}$ determines $m$ modes of forcing from $\boldsymbol{W}_{\!t}$ for each of the $n$ components of $\boldsymbol{Y}_{\!t}$ .

Figure 3. Sample space $\varOmega$ as the domain of random variables corresponding to coordinates $\boldsymbol{X}$ and field variables $\boldsymbol{Y}_{\!t}=\varphi _{\!t}(\boldsymbol{X})$ . (a) Two fields $Y_{\!t}^{1}$ (dark blue isolines) and $Y_{\!t}^{2}$ (red/pink filled isoregions) over the region in space defined by the density $f_{\boldsymbol{X}}$ . In panel (a), the arrows labelled ‘in’ and ‘out’ denote fluxes of $\boldsymbol{Y}_{\!t}$ at a fixed value of $Y_{\!t}^{1}=c_{1}$ over a boundary region defined by $\boldsymbol{\nabla }f_{\boldsymbol{X}}$ (see § 3.2). (b) Section through the current and future probability density $f_{\boldsymbol{Y}}$ with thick and dashed blue lines, respectively. As demonstrated in § 3 and illustrated by the ‘in‘ and ‘out‘ arrows in panel (b), the boundary fluxes lead to drift in the equation governing the joint distribution $f_{\boldsymbol{Y}}$ (for $y^{1}=c_{1}$ , the drift is in the positive $y^{2}$ direction for those values of $Y^{2}$ for which $\boldsymbol{J}\boldsymbol{\cdot }\boldsymbol{\nabla }f_{\boldsymbol{X}}\gt 0$ and in the negative $y^{2}$ direction for those values of $Y^{2}$ for which $\boldsymbol{J}\boldsymbol{\cdot }\boldsymbol{\nabla }f_{\boldsymbol{X}}\lt 0$ ). At the same time, $f_{\boldsymbol{Y}}$ is made narrower due to irreversible mixing, which homogenises the fields.

In what follows, it should be remembered that with $\varOmega$ regarded as a sample space, the precise relationship $\boldsymbol{Y}_{\!t}=\varphi _{\!t}(\boldsymbol{X})$ is not assumed to be known. In this context, the gradient $\boldsymbol{\nabla }\boldsymbol{Y}_{\!t}$ should be regarded as an unknown random variable, rather than a known gradient of $\varphi _{\!t}(\boldsymbol{X})$ . It is unknown in the sense that (3.1) does not provide $\boldsymbol{\nabla }\boldsymbol{Y}_{\!t}$ with its own governing equation. Noting that $\boldsymbol{J}_{\!t}$ typically depends on $\boldsymbol{\nabla }\boldsymbol{Y}_{\!t}$ , the system (3.1) is therefore unclosed because it consists of more unknowns than equations. If spatial derivatives were included in $\boldsymbol{Y}_{\!t}$ , their evolution equations would include higher derivatives and one would be faced with the infinite hierarchy of equations that constitutes the closure problem of turbulence (see McComb Reference McComb1990, and § 8 for a discussion of closures in the context of the present work).

To determine the evolution of the probability distribution of $\boldsymbol{Y}_{\!t}$ , consider the infinitesimal generator $\mathscr{L}$ (Pavliotis Reference Pavliotis2014, § 2.3) acting on an observable $g:\mathbb{R}^{n}\rightarrow \mathbb{R}$ :

(3.2) \begin{equation} \mathscr{L}g(\boldsymbol{y},s):=\lim \limits _{t\rightarrow s}\frac {\mathbb{E}[g(\boldsymbol{Y}_{\!t})|\boldsymbol{Y}_{\!s}=\boldsymbol{y}]-g(\boldsymbol{y})}{t-s}, \end{equation}

where the conditional expectation accounts for the behaviour of all points within the domain where $\boldsymbol{Y}_{\!s}=\boldsymbol{y}$ . In a deterministic context, without the expectation around $g(\boldsymbol{Y}_{\!t})$ , (3.2) corresponds to the infinitesimal generator of the Koopman operator (Brunton et al. Reference Brunton, Budišić, Kaiser and Kutz2022).

Using Itô’s formula to determine the time rate of change of the observable $g$ , the limit in (3.2) can, noting that the time integral of the resulting Wiener process has an expected value of zero, be expressed as (see § 3.4 and Lemma 3.2, Pavliotis Reference Pavliotis2014):

(3.3) \begin{equation} \mathscr{L}g(\boldsymbol{y},s)=\mathbb{E}\left [\left(I_{s}^{i}-\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{J}_{s}^{i}\right)\dfrac {\partial g(\boldsymbol{Y}_{\!s})}{\partial Y_{s}^{i}}\bigg |\boldsymbol{Y}_{\!s}=\boldsymbol{y}\right ]+\frac {1}{2}\varSigma ^{\textit{ij}}_{s}\frac {\partial ^{2}g(\boldsymbol{y})}{\partial y^{i}\partial y^{j}}, \end{equation}

where

(3.4) \begin{equation} \boldsymbol{\varSigma }_{s}:=\mathbb{E}\big[\boldsymbol{\sigma }\boldsymbol{\sigma }^{\top }|\boldsymbol{Y}_{\!s}=\boldsymbol{y}\big], \end{equation}

and Einstein summation notation has been used over components of $\boldsymbol{Y}_{\!s}$ .

We will follow standard arguments to construct backward and forward Kolmogorov equations in terms of the generator $\mathscr{L}$ (see, e.g. Gardiner Reference Gardiner1990; Pavliotis Reference Pavliotis2014). Along the way, integration by parts will be applied to produce boundary and irreversible interior mixing terms that will appear in the associated Kolmogorov equations as forcing and anti-diffusion terms, respectively.

3.2. Transport and mixing

To understand the physics behind the evolution of the probability density of $\boldsymbol{Y}_{\!t}$ with respect to a prescribed spatial sample distribution $f_{\boldsymbol{X}}$ , it is useful to decompose the transport terms in (3.3) into those associated with internal irreversible mixing and transport across the sample distribution’s ‘boundary’, which will now be defined.

If the sample distribution $f_{\boldsymbol{X}}$ is uniform over a subset of $\varOmega$ , the boundary associated with $f_{\boldsymbol{X}}$ is simply the boundary of the associated subset. However, more generally, $f_{\boldsymbol{X}}$ could be non-uniform and would therefore not correspond to a classical set with a sharp bounding surface. As explained later and illustrated in figure 4, a natural definition of the bounding region in this regard proves to be the unnormalised density $|\boldsymbol{\nabla }f_{\boldsymbol{X}}|$ , which defines a sharp surface as $f_{\boldsymbol{X}}$ becomes proportional to an indicator function. The approach therefore generalises the traditional notion of a bounding surface in a way that can be expressed in terms of probability distributions. The manipulations described later are relatively simple because of (rather than in spite of) the generality of this approach.

To establish a connection between expectations obtained via conditioning on $\boldsymbol{X}$ compared with conditioning on $\boldsymbol{Y}$ , recall the elementary property

(3.5) \begin{equation} \mathbb{E}[\square ]=\mathbb{E}[\mathbb{E}[\square |\boldsymbol{Y}]] =\mathbb{E}[\mathbb{E}[\square |\boldsymbol{X}]]=\int \limits _{\mathbb{R}^{d}}\mathbb{E}[\square |\boldsymbol{X}=\boldsymbol{x}]f_{\boldsymbol{X}}(\boldsymbol{x})\,\textrm{d} \boldsymbol{x}, \end{equation}

which means that expectations can be computed from those conditioned on $\boldsymbol{Y}$ or those conditioned on $\boldsymbol{X}$ . The latter are useful because they correspond to variables in physical space to which it is possible to apply integration by parts to distinguish boundary fluxes from interior mixing. In this regard, note that use of the product and chain rule implies

(3.6) \begin{equation} -(\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{J})^{i}\dfrac {\partial g(\boldsymbol{Y})}{\partial Y^{i}} = -\boldsymbol{\nabla }\boldsymbol{\cdot }\left (\boldsymbol{J}^{i}\dfrac {\partial g(\boldsymbol{Y})}{\partial Y^{i}}\right ) +\boldsymbol{J}^{i}\boldsymbol{\cdot }\boldsymbol{\nabla }Y^{j}\frac {\partial ^{2}g(\boldsymbol{Y})}{\partial Y^{i}\partial Y^{j}}. \end{equation}

Figure 4. (a) the spatial sample distribution $f_{\boldsymbol{X}}(\boldsymbol{x}):=N_{k}/(1+\exp (k p(\boldsymbol{x})))$ , where $p(\boldsymbol{x}):=x_{1}^{4}+x_{2}^{4}+x_{1}x_{2}-1$ and $N_{k}$ is a normalisation constant. (b) the bounding region of $f_{\boldsymbol{X}}$ defined by the unnormalised density $|\boldsymbol{\nabla }f_{\boldsymbol{X}}|$ (red) and vector field $\boldsymbol{\nabla }f_{\boldsymbol{X}}$ (blue). As $k\rightarrow \infty$ , the distribution describes well-defined subsets of $\mathbb{R}^{2}$ with a sharp boundary.

It will be assumed that $f_{\boldsymbol{X}}(\boldsymbol{x})$ either has compact support on $\mathbb{R}^{d}$ or tends to zero sufficiently fast as $|\boldsymbol{x}|\rightarrow \infty$ . Then, making use of (3.5), (3.6) and applying integration by parts to the first term on the right-hand side of (3.6) yields

(3.7) \begin{align} -\mathbb{E}\left [(\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{J})^{i}\dfrac {\partial g(\boldsymbol{Y})}{\partial Y^{i}}\right ] &= \int \limits _{\mathbb{R}^{d}}\mathbb{E}\left [\boldsymbol{J}^{i}\dfrac {\partial g(\boldsymbol{Y})}{\partial Y^{i}}\bigg |\boldsymbol{X}=\boldsymbol{x}\right ]\boldsymbol{\cdot }\boldsymbol{\nabla }f_{\boldsymbol{X}}(\boldsymbol{x})\,\textrm{d} \boldsymbol{x}\nonumber\\ &\quad +\int \limits _{\mathbb{R}^{d}}\mathbb{E}\left [\boldsymbol{J}^{i}\boldsymbol{\cdot }\boldsymbol{\nabla }Y^{j}\frac {\partial ^{2}g(\boldsymbol{Y})}{\partial Y^{i}\partial Y^{j}}\bigg |\boldsymbol{X}=\boldsymbol{x}\right ] f_{\boldsymbol{X}}(\boldsymbol{x})\,\textrm{d} \boldsymbol{x}, \end{align}

where the first and second term on the right-hand side correspond to boundary transport and internal mixing, respectively. The first can be expressed in terms of the original density by including the weight $\boldsymbol{\nabla }f_{\boldsymbol{X}}/f_{\boldsymbol{X}}=\boldsymbol{\nabla }\log f_{\boldsymbol{X}}$ inside the expectation:

(3.8) \begin{equation} \int \limits _{\mathbb{R}^{d}}\mathbb{E}[\square |\boldsymbol{X}=\boldsymbol{x}]\boldsymbol{\cdot }\boldsymbol{\nabla }f_{\boldsymbol{X}}(\boldsymbol{x})\,\textrm{d} \boldsymbol{x} = \int \limits _{\mathbb{R}^{d}}\mathbb{E}\left [\square \boldsymbol{\cdot }\boldsymbol{\nabla }\log f_{\boldsymbol{X}}\big |\boldsymbol{X}=\boldsymbol{x}\right ]f_{\boldsymbol{X}}(\boldsymbol{x})\,\textrm{d} \boldsymbol{x}; \end{equation}

hence, using (3.5),

(3.9) \begin{equation} -\mathbb{E}\left [(\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{J})^{i}\dfrac {\partial g(\boldsymbol{Y})}{\partial Y^{i}}\bigg |\boldsymbol{Y}\right ] = \underbrace {\mathbb{E}\left [\boldsymbol{J}^{i}\dfrac {\partial g(\boldsymbol{Y})}{\partial Y^{i}}\boldsymbol{\cdot }\boldsymbol{\nabla }\log f_{\boldsymbol{X}}\big |\boldsymbol{Y}\right ]}_{\textrm{boundary transport}}+ \underbrace {\mathbb{E}\left [\boldsymbol{J}^{i}\boldsymbol{\cdot }\boldsymbol{\nabla }Y^{j}\frac {\partial ^{2}g(\boldsymbol{Y})}{\partial Y^{i}\partial Y^{j}}\bigg |\boldsymbol{Y}\right ]}_{\textrm{internal mixing}}. \end{equation}

Physically, the factor $\boldsymbol{\nabla }\log f_{\boldsymbol{X}}$ conditions the flux against the normal direction of the boundary region. To see why, note that $-\boldsymbol{\nabla }f_{\boldsymbol{X}}/|\boldsymbol{\nabla }f_{\boldsymbol{X}}|$ defines an outward unit vector that is perpendicular to isosurfaces of $f_{\boldsymbol{X}}$ and that $|\boldsymbol{\nabla }f_{\boldsymbol{X}}|$ is the (unnormalised) density associated with the bounding region. Therefore, introducing the normalisation factor $N$ and using $\tilde {\mathbb{E}}$ to denote the expectation of a random variable $\square$ over the boundary region:

(3.10) \begin{equation} \tilde {\mathbb{E}}[\square ]:=\frac {1}{N}\int \limits _{\mathbb{R}^{d}}\square |\boldsymbol{\nabla }f_{X}|\,\textrm{d} \boldsymbol{x}, \end{equation}

the formula for conditional expectations under a change in measure, removing a factor of $N^{-1}$ from both sides, yields (see, for example Shreve Reference Shreve2004, § 5.2.1)

(3.11) \begin{equation} \mathbb{E}\left [\square \boldsymbol{\cdot }\boldsymbol{\nabla }\textrm{log}f_{\boldsymbol{X}}|\boldsymbol{Y}\right ]= \tilde {\mathbb{E}}\left [\square \boldsymbol{\cdot }\frac {\boldsymbol{\nabla }f_{\boldsymbol{X}}}{|\boldsymbol{\nabla }f_{\boldsymbol{X}}|}\big |\boldsymbol{Y}\right ] \mathbb{E}\left [\frac {|\boldsymbol{\nabla }f_{\boldsymbol{X}}|}{f_{\boldsymbol{X}}}\big |\boldsymbol{Y}\right ]. \end{equation}

The first term on the right-hand side of (3.11) is the expected flux $\square$ into the region defined by $f_{\boldsymbol{X}}$ . The second term is the expected value of the normalised density ratio $N^{-1}|\boldsymbol{\nabla }f_{\boldsymbol{X}}|/f_{\boldsymbol{X}}$ multiplied by $N$ and conditioned on $\boldsymbol{Y}$ . Now consider the probability that $\boldsymbol{Y}\in \mathcal{Y}$ (where $\mathcal{Y}\subseteq \mathbb{R}^{n}$ is an arbitrary subset) over the boundary region using the indicator function $\mathbb{I}_{\mathcal{Y}}$ :

(3.12) \begin{equation} \tilde {\mathbb{P}}\{\boldsymbol{Y}\in \mathcal{Y}\}:=\tilde {\mathbb{E}}[\mathbb{I}_{\mathcal{Y}}]=\mathbb{E}\left [\mathbb{I}_{\mathcal{Y}}\frac {|\boldsymbol{\nabla }f_{\boldsymbol{X}}|}{Nf_{\boldsymbol{X}}}\right ]=\int \limits _{\mathcal{Y}}\underbrace {\mathbb{E}\left [\frac {|\boldsymbol{\nabla }f_{\boldsymbol{X}}|}{f_{\boldsymbol{X}}}\big |\boldsymbol{Y}=\boldsymbol{y}\right ]\frac {f_{\boldsymbol{Y}}(\boldsymbol{y})}{N}}_{=:\tilde {f}_{\boldsymbol{Y}}(\boldsymbol{y})}\,\textrm{d}\boldsymbol{y}, \end{equation}

which defines the probability density $\tilde {f}_{\boldsymbol{Y}}$ of $\boldsymbol{Y}$ when $\boldsymbol{X}$ is sampled from the boundary region and enables the final term in (3.11) to be expressed as $N\tilde {f}_{\boldsymbol{Y}}/f_{\boldsymbol{Y}}$ . The normalisation constant $N$ can be interpreted as the relative size of the boundary region. In particular, as $f_{\boldsymbol{X}}$ tends to the indicator function $\mathbb{I}_{\mathcal{X}}$ , selecting a subset $\mathcal{X}\subseteq \mathbb{R}^{d}$ , $N$ tends to the ratio of the surface area of $\mathcal{X}$ to its volume. If, for example, $\mathcal{X}$ is a unit interval of $\mathbb{R}$ , the unnormalised distribution $|\boldsymbol{\nabla }f_{\boldsymbol{X}}|$ would tend to a delta measure at either end of the interval and the required normalisation constant $N$ would, therefore, tend to $2$ .

3.3. Backward Kolmogorov equation

For random variables $\square$ and $\boldsymbol{Y}$ , and a function $g:\mathbb{R}^{n}\rightarrow \mathbb{R}$ , the expected value of $g(\boldsymbol{Y})$ , conditional on $\boldsymbol{Y}=\boldsymbol{y}$ , is $g(\boldsymbol{y})$ . Therefore, $\mathbb{E}[\square g(\boldsymbol{Y})|\boldsymbol{Y}=\boldsymbol{y}]=\mathbb{E}[\square |\boldsymbol{Y}=\boldsymbol{y}]g(\boldsymbol{y})$ and (3.3) becomes

(3.13) \begin{equation} \mathscr{L}g(\boldsymbol{y})=\unicode{x1D63F}_{1}^{\,i}\dfrac {\partial g}{\partial y^{i}}+\unicode{x1D63F}_{2}^{\textit{,ij}}\frac {\partial ^{2}g}{\partial y^{i}\partial y^{j}}. \end{equation}

Using the results from § 3.2, the so-called drift velocity is

(3.14) \begin{equation} \unicode{x1D63F}_{1}:=\mathbb{E}\left [\boldsymbol{I}_{s}+\boldsymbol{J}_{s}\boldsymbol{\cdot }\boldsymbol{\nabla }\log f_{\boldsymbol{X}}\big |\boldsymbol{Y}_{\!s}=\boldsymbol{y}\right ] \end{equation}

and the symmetric diffusion coefficient is

(3.15) \begin{equation} \unicode{x1D63F}_{2}:=\frac {1}{2}\mathbb{E}\left [\boldsymbol{J}^{\top }_{s}\boldsymbol{\nabla }\boldsymbol{Y}_{\!s}+\boldsymbol{\nabla }\boldsymbol{Y}_{\!s}^{\top }\boldsymbol{J}_{s}|\boldsymbol{Y}_{\!s}=\boldsymbol{y}\right ]+\frac {1}{2}\boldsymbol{\varSigma }. \end{equation}

If $u(\boldsymbol{y},s):=\mathbb{E}[g(\boldsymbol{Y}_{\!t})|\boldsymbol{Y}_{\!s}=\boldsymbol{y}]$ , then

(3.16) \begin{equation} -\partial _{s}u(\boldsymbol{y},s)=\mathscr{L}u(\boldsymbol{y},s) \end{equation}

is solved backwards in time from the end condition $u(\boldsymbol{y},t)=g(\boldsymbol{y})$ , so that $u(\boldsymbol{y},s)$ is the expected value of $g(\boldsymbol{Y}_{\!t})$ given $\boldsymbol{Y}_{\!s}=\boldsymbol{y}$ for $s\leqslant t$ . The generator $\mathscr{L}$ therefore ‘pulls back’ the observation $g$ along the dynamics specified by (3.1).

A local (in space) diffusivity with what looks like the same form as (3.15) is found in classical p.d.f. methods (Pope Reference Pope1981). Here, however, the spatial conditioning is distributed over an arbitrary control volume, as defined by the sample distribution $f_{\boldsymbol{X}}$ (see item (i) at the end of § 1), which leads to the boundary fluxes in (3.14) and has been derived from the dual perspective of observables $g$ (see item (ii) at the end of § 1). Consequently, (3.15) readily accounts for distributed stochastic forcing through $\boldsymbol{\varSigma }$ . The following section explains how $\unicode{x1D63F}_{1}$ and $\unicode{x1D63F}_{2}$ determine the evolution of probability density using duality.

3.4. Forward Kolmogorov equation

The ‘forward’ equation corresponding to (3.16) is derived by expressing the time rate of change of the conditional expectation $\mathbb{E}[g(\boldsymbol{Y}_{\!t})|\boldsymbol{Y}_{\!s}=\boldsymbol{y}]$ in terms of a conditional probability density $f_{\boldsymbol{Y}_{\!t}|\boldsymbol{Y}_{\!s}}$ , as described by Pavliotis (Reference Pavliotis2014, pp. 45–50). Here, we present an abridged version to highlight the main idea without introducing new notation.

If the density $f_{\boldsymbol{Y}}(\boldsymbol{y},t)$ is understood as being conditional on a specified density $f_{\boldsymbol{Y}}(\boldsymbol{y},0)$ at $t=0$ , the expected value of $g$ resulting from the initial density $f_{\boldsymbol{Y}}(\boldsymbol{y},0)$ satisfies

(3.17) \begin{equation} \partial _{t}\mathbb{E}[g(\boldsymbol{Y}_{\!t})]=\partial _{t}\int \limits _{\mathbb{R}^{n}}g(\boldsymbol{y})f_{\boldsymbol{Y}}(\boldsymbol{y},t)\,\textrm{d}\boldsymbol{y}=\int \limits _{\mathbb{R}^{n}}g(\boldsymbol{y})\partial _{t}f_{\boldsymbol{Y}}(\boldsymbol{y},t)\,\textrm{d}\boldsymbol{y}. \end{equation}

Alternatively, using $\mathscr{L}$ , integrating by parts and assuming that either $f_{\boldsymbol{Y}}$ or $g$ decay sufficiently fast as $|\boldsymbol{y}|\rightarrow \infty$ ,

(3.18) \begin{equation} \partial _{t}\mathbb{E}[g(\boldsymbol{Y}_{\!t})]=\int \limits _{\mathbb{R}^{n}}\mathscr{L}g(\boldsymbol{y})f_{\boldsymbol{Y}}(\boldsymbol{y},t)\,\textrm{d} \boldsymbol{y}=\int \limits _{\mathbb{R}^{n}}g(\boldsymbol{y})\mathscr{L}^{\dagger }f_{\boldsymbol{Y}}(\boldsymbol{y},t)\,\textrm{d} \boldsymbol{y}. \end{equation}

Therefore, equating (3.17) and (3.18), and noting that they are valid for all suitable observables $g$ ,

(3.19) \begin{equation} \partial _{t}f_{\boldsymbol{Y}}(\boldsymbol{y},t)=\mathscr{L}^{\dagger }f_{\boldsymbol{Y}}(\boldsymbol{y},t), \end{equation}

where

(3.20) \begin{equation} \mathscr{L}^{\dagger }f_{\boldsymbol{Y}}(\boldsymbol{y},t):=-\dfrac {\partial }{\partial y^{i}}(\unicode{x1D63F}_{1}^{\,i}f_{\boldsymbol{Y}}(\boldsymbol{y},t))+\frac {\partial ^{2}}{\partial y^{i}\partial y^{j}}(\unicode{x1D63F}^{\,\textit{ij}}_{2}f_{\boldsymbol{Y}}(\boldsymbol{y},t)). \end{equation}

The forward equation evolves $f_{\boldsymbol{Y}}$ forwards in time from the initial density $f_{\boldsymbol{Y}}(\boldsymbol{y},0)$ .

3.5. Physical interpretation

Sources $\boldsymbol{I}$ and boundary fluxes $\boldsymbol{J}\boldsymbol{\cdot }\boldsymbol{\nabla }f_{\boldsymbol{X}}/|\boldsymbol{\nabla }f_{\boldsymbol{X}}|$ into/out of the control volume contribute to positive/negative drift $\unicode{x1D63F}_{1}$ in the evolution of the density $f_{\boldsymbol{Y}}$ . As discussed in § 3.2 beneath (3.12), the boundary fluxes are scaled to account for the relative size of the boundary. Both processes contribute to drift because they determine a rate of gradual change in $\boldsymbol{Y}$ , rather than an immediate addition or removal of fluid with a given value of $\boldsymbol{Y}$ . An exception to this interpretation occurs when advection due to a velocity field is included in $\boldsymbol{J}$ , which will be discussed in § 4.

Expectations of the cross-gradient mixing terms $\boldsymbol{J}^{\top }\boldsymbol{\nabla }\boldsymbol{Y}$ determine the diffusion coefficient $\unicode{x1D63F}_{2}$ in (3.15). The sign of $\unicode{x1D63F}_{2}$ will be discussed in detail in § 5. In many practical applications, and according to Fick’s law, each column of $\boldsymbol{J}$ points in the direction of each column of $-\boldsymbol{\nabla }\boldsymbol{Y}$ , making it reasonable to expect $\unicode{x1D63F}_{2}$ to be negative semidefinite. This property corresponds to the fact that down-gradient molecular transport homogenises $\boldsymbol{Y}$ , which leads to greater certainty in the value of $\boldsymbol{Y}_{\!t}$ and lower entropy of the distribution $f_{\boldsymbol{Y}}$ .

The following sections (§§ 4 and 5) consider fluxes that comprise advection and diffusion:

(3.21) \begin{equation} \boldsymbol{J}=\boldsymbol{U}\otimes \boldsymbol{Y}-\boldsymbol{\nabla }\boldsymbol{Y}\boldsymbol{\alpha }, \end{equation}

where $\boldsymbol{U}$ is a solenoidal velocity field (i.e. $\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{U}\equiv 0$ ) and $\boldsymbol{\alpha }\in \mathbb{R}^{n\times n}$ is a (typically diagonal) multicomponent diffusion matrix. More generally, $\boldsymbol{\alpha }$ could also account for anisotropic diffusion, which would make it a four-dimensional tensor.

7.1. Advection and diffusion by the ABC flow

An Arnold–Beltrami–Childress (ABC) flow with coefficients equal to unity is the three-dimensional divergence-free velocity field

(7.1) \begin{equation} \boldsymbol{U} := (\sin (X^3) + \cos (X^2),\ \sin (X^1) + \cos (X^3),\ \sin (X^2) + \cos (X^1))^{\top }. \end{equation}

It is an exact solution of the Euler equations in a periodic domain and is known to exhibit chaotic streamlines (Dombre et al. Reference Dombre, Frisch, Greene, Hénon, Mehr and Soward1986), which makes it an ideal candidate to study mixing. Couched in the format of (3.1), the evolution of a passive scalar concentration $Y_t$ due to the combined effects of advection by $\boldsymbol{U}$ and diffusion is

(7.2) \begin{equation} \dfrac {\textrm{d}{Y_t}}{\textrm{d}{t}} = -\boldsymbol{U} \boldsymbol{\cdot }\boldsymbol{\nabla }Y_t + \alpha \Delta Y_t, \end{equation}

where $\alpha$ denotes the constant scalar diffusivity. The corresponding forward Kolmogorov equation governing the evolution of $f_{Y}(y,t)$ is given by (3.19):

(7.3) \begin{equation} \frac {\partial f_Y}{ \partial t} = -\frac {\partial ^2 }{ \partial y^2} \big ( \underbrace {\alpha \mathbb{E}[|\boldsymbol{\nabla }Y_t|^2|Y_t=y]}_{=:-\unicode{x1D63F}_{2}} f_{Y} \big ). \end{equation}

The coefficients $V$ and $\unicode{x1D63F}_{1}$ from § 3.3 are both zero because the periodic domain does not have boundaries. Although (7.3) cannot be solved, as $\unicode{x1D63F}_{2}$ is unknown, we can examine estimates of $\unicode{x1D63F}_{2}$ by solving (7.2) numerically. Space is treated using a Fourier pseudo-spectral method and time using a Crank–Nicolson discretisation with a resolution of $64^3$ modes and time step of $\Delta t = 5 \times 10^{-3}$ , respectively. Choosing an initial condition consisting of a front separating two regions of different concentration $Y_{0} = \tanh (10\boldsymbol{X}^{\top }\boldsymbol{1})$ and setting $\alpha =1/10$ , we numerically integrate with respect to $t\in [0,4]$ . Then, by approximating the density $f_{Y}(y,t)$ and $\mathbb{E}[|\boldsymbol{\nabla }Y_t|^2|Y_t=y]$ with a histogram, constructed from a single snapshot of field data, we estimate the terms in (7.3).

Figure 8. (a) Time evolution of a cross-section of the concentration field $Y_{\!t}$ and (b) its corresponding probability density $f_{Y}(y,t)$ (shaded) and (negative) diffusion coefficient $\alpha \mathbb{E}[|\boldsymbol{\nabla }Y_t|^2|Y_{\!t}=y]$ . https://www.cambridge.org/S0022112025106319/JFM-Notebooks/files/fig08/fig08.ipynb.

A time evolution of a cross-section of the scalar field $Y_{t}$ alongside its corresponding density $f_{Y}(y,t)$ and the conditional expectation $\mathbb{E}[|\boldsymbol{\nabla }Y_t|^2|Y_t=y]$ are shown in figure 8. At $t=0.5$ , the scalar field, consisting of two regions of approximately uniform concentration separated by a relatively sharp interface, is represented in probability space by spikes at $y = \pm 1$ . As the concentration field is mixed irreversibly, the amplitude of these spikes reduces and the density associated with values of $y$ in the vicinity of zero increases. Mixing subsequently homogenises the scalar field to the extent that no evidence of the initial distribution remains by $t=2$ . As there is no injection of concentration to balance the homogenising effects of diffusion, $Y_t$ will ultimately tend to a constant uniformly over the domain and the density $f_{Y}(y,t)$ in the limit $t\rightarrow \infty$ would tend to a Dirac distribution.

The expectation $\mathbb{E}[|\boldsymbol{\nabla }Y_t|^2|Y_{\!t}=y]$ , which quantifies the rate of mixing and, therefore, negative diffusion in (7.3) is maximised by $y=0$ . For $t=0.5$ and $t=1.0$ , scalar concentrations that are less probable are associated with the most mixing. This observation is not surprising because, as discussed in the examples in § 2 and illustrated in figure 1, large gradients are associated with small probability densities, unless they occur over sufficiently many branches of the function’s inverse. In this regard, it should also be noted that $\mathbb{E}[|\boldsymbol{\nabla }Y_t|^2|Y_t=y]$ is zero for the minimum ( $y=-1$ ) and maximum ( $y=1$ ) values of $Y_{\!t}$ , at which $\boldsymbol{\nabla }Y_t=\boldsymbol{0}$ . In the absence of diffusion at these values of $y$ , the density $f_{Y}$ consequently maintains compact support (i.e. mixing interpolates and cannot produce values that lie outside the range of values that were there in the first place).

7.2. Stochastic boundary conditions

To sustain a finite variance of a diffusive scalar field, it is necessary to apply a forcing. Given that a wide variety of physical systems are sustained by a boundary forcing, often characterised by substantial fluctuations and uncertainty, it is natural to consider how these forces manifest in the evolution equation governing the system’s probability density. A simple example of such a problem that is sufficient to highlight several of the terms discussed in § 3 is the diffusion of heat through a one-dimensional rod of unit length forced by Ornstein–Uhlenbeck processes (i.e. normally distributed thermal fluctuations) at the boundary. When cast in the form of (3.1), for unit thermal diffusivity ( $\alpha =1$ ), the temperature $Y_{\!t}$ evolves according to

(7.4) \begin{equation} {\dfrac {\textrm{d}{Y_t}}{\textrm{d}{t}}} = \partial ^{2}_{X}Y_t\quad \textrm{for}\quad X\in [0,1], \end{equation}

where values of $Y_{\!t}$ for points on the boundary $\partial \varOmega$ (corresponding to $X\in \{0,1\}$ ) evolve according to

(7.5) \begin{equation} \,\textrm{d} Y_t(\omega ) = -aY_t(\omega )\,\textrm{d} t + \sigma \,\textrm{d} W_t(\omega )\quad \textrm{for}\quad \omega \in \partial \varOmega , \end{equation}

where $a,\sigma$ are real constants and $W_t(\omega )$ is a Wiener process. In the context of (3.1), $J_{\!t}=-\partial _{X}Y_{\!t}$ is the (diffusive) scalar flux. Since we are interested in the entire domain, the sampling distribution is $f_{X}\equiv 1$ and $\partial _{X} \log f_{X}$ in the calculation of $\unicode{x1D63F}_{1}$ should be regarded as changing the distribution under which the expectation is calculated from $f_{X}$ to point evaluation at the (two) boundaries (see discussion following (3.9) in § 3.2):

(7.6) \begin{equation} \mathbb{E}[J_{\!t}\partial _{X}\log f_{X}|Y_{\!t}]f_{Y}=2\mathbb{E}[\partial _{n}Y_{\!t}|Y_{\!t}]\tilde {f}_{Y}, \end{equation}

where $\partial _{n}Y_{\!t}$ is the boundary normal derivative of $Y_{\!t}$ , $\tilde {f}_{Y}$ is the probability density of $Y_{\!t}$ , conditional on evaluation at a boundary, and the factor $2$ accounts for there being two boundaries separated by a distance of one unit.

Figure 9. (a) Space–time plot of $Y_t$ for the one-dimensional diffusion (7.4) forced by an Ornstein–Ulhenbeck process (7.5) with $a=5,\sigma =1$ , and (b) joint density $\tilde {f}_{Y, \partial _{n}Y}$ conditioned on the boundaries. The positive covariance of this plot reflects the tendency of heat to flow down gradients at the boundaries. https://www. cambridge.org/S0022112025106319/JFM-Notebooks/files/fig09-10/fig09-10.ipynb.

Figure 9(a) shows a space–time plot of the system, for which the temperature at the boundaries fluctuates on the integral time scale ( $1/a$ ) of the Ornstein–Uhlenbeck processes (7.5). The solution is obtained numerically using second-order finite differences in space with $32$ grid points and an Euler–Maruyama discretisation in time with a time step of $2.5 \times 10^{-3}$ . Towards the centre of the domain, the effects of diffusion reduce the amplitude of the fluctuations.

Using (7.6), the forward Kolmogorov equation, which describes how the probability density $f_{Y}(y,t)$ evolves in time, is given by (3.20):

(7.7) \begin{equation} \frac {\partial f_Y}{ \partial t} = -\frac {\partial }{ \partial y} \Big ( \underbrace {2\mathbb{E}[ \partial _{n}Y_t|Y_{\!t}=y] \frac {\tilde {f}_{Y}}{f_{Y}}}_{=:\unicode{x1D63F}_{1}} f_{Y} \Big ) - \frac {\partial ^2 }{ \partial y^2} \big ( \underbrace {\mathbb{E}[|\partial _{X} Y_t|^2|Y_{\!t}=y]}_{=:-\unicode{x1D63F}_{2}} f_{Y} \big ), \end{equation}

According to (7.5), $Y_{\!t}$ sampled at the boundary is an Ornstein–Uhlenbeck process and therefore has a stationary density given by $\tilde {f}_{Y}(y)=\sqrt {a/(\pi \sigma ^{2})}\,\textrm{exp}(-ay^{2}/\sigma ^{2})$ (Gardiner Reference Gardiner1990), which, using the approach described following (3.12), is used to express $\unicode{x1D63F}_{1}$ .

As there is no advection in (7.4), the coefficient $V$ from § 4 is zero. However, due to conduction and, in turn, boundary fluxes of heat, (7.7) has non-zero drift and diffusion terms. Although (7.7) cannot be solved prognostically without a model for the unknown coefficients, the coefficients can be approximated numerically by constructing an ensemble satisfying (7.4) and (7.5).

Figure 10 shows the drift $\unicode{x1D63F}_{1}$ and diffusion $\unicode{x1D63F}_{2}$ coefficients, as well as the stationary density $f_{Y}(y)$ using time and ensemble averaging. An ensemble of 500 paths each containing 1950 snapshots in time was constructed numerically by simulating (7.4), (7.5) for $t\in [0,5]$ and discarding the initial 50 snapshots. In § 7.1, molecular diffusion acted to homogenise the scalar field $Y_t$ in the ABC flow and drive its corresponding density towards a Dirac measure. Were the value of $Y_t$ at the boundaries in this example to be fixed at zero, the distribution would, due to the negative diffusion $\unicode{x1D63F}_{2}$ (figure 10 d), also tend towards a Dirac measure at zero, irrespective of the precise initial conditions. Instead, random forcing at the system’s boundary creates variance that balances the destruction of variance described by $\unicode{x1D63F}_{2}$ . To illustrate how, figure 9(b) shows the conditional joint density $\tilde {f}_{Y, \partial _{n} Y_{\!t}}$ of $Y_t$ and its boundary-normal gradient $\partial _{n}Y_{\!t}$ . The conditional density reveals the expected positive correlation between $Y_t$ and $\partial _{n}Y_t$ , which means that heat flux into/out of the domain is typically accompanied by positive/negative temperatures at the boundary. Therefore, $\unicode{x1D63F}_{1}\lessgtr 0$ for $y\lessgtr 0$ , as shown in figure 10(a), which implies that $\unicode{x1D63F}_{1}$ corresponds to divergent transport of probability density away from the origin.

Figure 10. (a,b,d) Drift coefficient $\unicode{x1D63F}_{1}$ , diffusion coefficient $\unicode{x1D63F}_{2}$ and probability density $f_Y(y)$ , respectively, of the forward Kolmogorov equation corresponding to a stationary state of the system (7.4) with $a=5,\sigma =1$ . (c) A numerical verification of the balance between boundary forcing and mixing. https://www.cambridge.org/S0022112025106319/JFM-Notebooks/files/fig09-10/fig09-10.ipynb.

Computing the balance of the terms in the right-hand side of (7.7), as shown in figure 10(c), verifies (7.7). It is worth noting that while incorporating uncertainty into (7.4) required a Monte Carlo approach to be employed in order to estimate $f_Y(y)$ , such boundary conditions pose no particular additional difficulty in (7.7), which accommodates stochastic forcing naturally.

7.3. Boussinesq equations

The following example demonstrates that coordinates, regarded here as functions $\boldsymbol{X}$ over the sample space $\varOmega$ , can be included in the state vector $\boldsymbol{Y}_{\!t}$ . For instance, if one wishes to understand probability distributions over horizontal slices of a domain, the vertical coordinate can be included in $\boldsymbol{Y}_{\!t}$ . More generally, other functions, such as geopotential height, could be incorporated into $\boldsymbol{Y}_{\!t}$ . The equations developed in § 3 are therefore very general in being able to generate equations that are specific to a particular problem and circumvent the need for case-specific derivations.

In a Boussinesq context, the pointwise deterministic equations governing the behaviour of the velocity field $\boldsymbol{U}_{\!t}$ and the buoyancy $B_{t}$ are

(7.8) \begin{align} {\dfrac {\textrm{d}{\boldsymbol{U}_{\!t}}}{\textrm{d}{t}}}&=B_{t}\boldsymbol{e}_{3}-\boldsymbol{\nabla }P_{\!t} -\boldsymbol{U}_{\!t}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{U}_{\!t}+\alpha _{1}\Delta \boldsymbol{U}_{\!t},\nonumber\\ \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{U}_{\!t}&=0,\nonumber \\ {\dfrac {\textrm{d}{B_{t}}}{\textrm{d}{t}}}&=-\boldsymbol{U}_{\!t}\boldsymbol{\cdot }\boldsymbol{\nabla }{B}_{\!t} +\alpha _{2}\Delta B_{t}, \end{align}

where $\boldsymbol{e}_{3}$ is the unit vector in the vertical direction. Since the vertical direction plays a distinguished role in corresponding to the direction of gravity, let $\boldsymbol{Y}_{\!t}:=(W_{\!t},B_{t},Z)^{\top }$ , where $W_{\!t}:=\boldsymbol{U}_{\!t}\boldsymbol{\cdot }\boldsymbol{e}_{3}$ is the vertical velocity (not to be confused with the vector of Wiener processes $\boldsymbol{W}_{\!t}$ used in § 3) and $Z=X^{3}$ is the (time independent) vertical coordinate. Writing $\,\textrm{d} Z/\,\textrm{d} t=W_{\!t}-\boldsymbol{U}_{\!t}\boldsymbol{\cdot }\boldsymbol{\nabla }Z=0$ , and noting that $\alpha _{11}=\alpha _{1}$ , $\alpha _{22}=\alpha _{2}$ , $\alpha _{33}=0$ and $\alpha _{\textit{ij}}=0$ for $i\neq j$ implies that $V$ is given by (4.3),

(7.9) \begin{equation} \unicode{x1D63F}_{1}= \begin{pmatrix} b\\ 0\\ w \end{pmatrix} - \mathbb{E} \left [\begin{array}{@{}c|c@{}} \partial _{Z} P_{\!t}+\alpha _{1}\boldsymbol{\nabla }W_{\!t}\boldsymbol{\cdot }\boldsymbol{\nabla }\log f_{\boldsymbol{X}} & \\ \alpha _{2}\boldsymbol{\nabla }B_{t}\boldsymbol{\cdot }\boldsymbol{\nabla }\log f_{\boldsymbol{X}} & \boldsymbol{Y}_{\!t}=\boldsymbol{y} \\ 0 & \end{array}\right ], \end{equation}

and

(7.10) \begin{equation} \unicode{x1D63F}_{2}=-\frac {1}{2} \mathbb{E} \left [\begin{array}{@{}ccc|c@{}} 2\alpha _{1} |\boldsymbol{\nabla }W_{\!t}|^{2} & (\alpha _{1}+\alpha _{2})\boldsymbol{\nabla }W_{\!t}\boldsymbol{\cdot }\boldsymbol{\nabla }B_{t} & \alpha _{1} \partial _{Z}W_{\!t} & \\ (\alpha _{1}+\alpha _{2})\boldsymbol{\nabla }W_{\!t}\boldsymbol{\cdot }\boldsymbol{\nabla }B_{t} & 2\alpha _{2}|\boldsymbol{\nabla }B_{t}|^{2} & \alpha _{2} \partial _{Z}B_{t} & \boldsymbol{Y}_{\!t}=\boldsymbol{y}\\ \alpha _{1} \partial _{Z}W_{\!t} & \alpha _{2} \partial _{Z}B_{t} & 0 & \end{array}\right ], \end{equation}

where $\boldsymbol{y}:=(w, b, z)^{\top }$ . In $\unicode{x1D63F}_{1}$ , $b$ is responsible for drift in the $w$ direction, because buoyancy increases $W_{\!t}$ and, in turn, $w$ is responsible for drift in the $z$ direction, because vertical velocity increases $Z$ . When interpreting the latter, it should be borne in mind that for the closed domain in this example $\mathbb{E}[W_{\!t}|Z]=0$ , which, therefore, does not affect the marginal distribution of $Z$ (i.e. the domain does not change its shape over time). The unknown conditionally averaged vertical pressure gradient also affects the evolution of the joint density through $\unicode{x1D63F}_{1}$ , whose remaining terms account for boundary fluxes of $W_{\!t}$ and $B_{t}$ , scaled by the relative size of the boundary.

While the drift $\unicode{x1D63F}_{1}$ is responsible for moving and stretching the joint density, the diffusion coefficient $\unicode{x1D63F}_{2}$ accounts for the effects of irreversible mixing. As discussed in § 5, $\unicode{x1D63F}_{2}$ is expected to be negative semidefinite ( $\unicode{x1D63F}_{2}\preceq 0$ ) in most applications, which means that it typically represents anti-diffusion. In particular, if $\alpha _{1}=\alpha _{2}=1$ , then $\unicode{x1D63F}_{2}\preceq 0$ can be guaranteed. For other combinations of $\alpha _{1}$ and $\alpha _{2}$ , $\unicode{x1D63F}_{2}\preceq 0$ depends on the correlations among $\boldsymbol{\nabla }W_{\!t}$ , $\boldsymbol{\nabla }B_{t}$ and $\boldsymbol{\nabla }Z$ , as discussed in § 5. The gradients, $\boldsymbol{\nabla }W_{\!t}$ and $\boldsymbol{\nabla }B_{t}$ are a priori unknown and therefore require closure (see § 8), which would involve postulating their dependence on the joint density $f_{\boldsymbol{Y}}$ and/or the independent variables $\boldsymbol{y}:=(w, b, z)^{\top }$ .

7.4. Rayleigh–Bénard convection with a stochastic boundary condition

This section describes an example of the domain partition discussed in § 6 applied to Rayleigh–Bénard convection. Here, the zonal sample distributions are defined using the diffusive vertical buoyancy flux and therefore separate the bulk central part of the flow from the top and bottom boundary layers (see figure 11, which displays instantaneous fields alongside weighted sample distributions $f_{\boldsymbol{X}|Z}$ for each zone). To illustrate how the partition yields useful information, we focus on the contributions to the equations governing $f_{\boldsymbol{Y}|Z}$ from $V$ and $\unicode{x1D63F}_{1}$ arising from transport between the zones.

Figure 11. (a) Snapshot of the vertical velocity $Y^{1}_{\!t}=W_{\!t}$ (dashed/solid white contours for negative/positive vertical velocity, respectively) and buoyancy field $Y^{2}_{\!t}=B_{t}$ (red is positive; blue is negative) following a negative fluctuation in the buoyancy at the bottom boundary. (b) Partition of the vertical domain into three zones $Z\in \{z_{1},z_{2},z_{3}\}$ defined by $f_{\boldsymbol{X}|Z}$ , where $f_{\boldsymbol{X}|Z}(\boldsymbol{x}|z_{2})$ (the central zone) is proportional to the horizontally and time-averaged vertical buoyancy flux $W_{\!t}B_{t}$ . In the boundary-layer zones above and below the bottom and top boundary, $f_{\boldsymbol{X}|Z}$ is proportional to the horizontally and time-averaged diffusive buoyancy flux $-\alpha _{2}\partial _{X^{3}}B_{t}$ . The constant of proportionality $\mathbb{P}\{Z=z\}$ for each zone is chosen to ensure that the sum of $f_{\boldsymbol{X}|Z}\mathbb{P}$ over the zones is the uniform distribution (black line) (cf. (6.1) in § 6). https://www.cambridge.org/S0022112025106319/JFM-Notebooks/files/fig11-12/fig11-12.ipynb.

The data for this example were obtained from direct numerical simulation of the Boussinesq equations (7.8) on a two-dimensional and horizontally periodic domain of unit height, as shown in figure 11(a). A spatially uniform buoyancy equal to $-1/2$ is imposed on the top boundary. On the bottom boundary, the imposed buoyancy is spatially uniform but follows an Ornstein–Uhlenbeck process in time, as described in § 7.2. The Ornstein–Uhlenbeck process has a mean of $1/2$ , unit variance and a time scale (corresponding to $1/a$ in § 7.2) of $2$ units. The governing equations can therefore be regarded as non-dimensionalised on the domain height and mean buoyancy difference between the horizontal boundaries. In terms of these scales, the Rayleigh number is $10^{7}$ , which corresponds to a (non-dimensionalised) viscosity and diffusivity of $\alpha _{1}=\sqrt {10^{7}}$ and $\alpha _{2}=\sqrt {10^{7}}$ , respectively, for a Prandtl number of unity.

The governing equations are approximated numerically using the pseudo-spectral method provided by Dedalus (Burns et al. Reference Burns, Vasil, Oishi, Lecoanet and Brown2020). Fourier and Chebyshev bases are used for the horizontal ( $X^{1}$ ) and vertical ( $X^{3}$ ) directions, with $512$ and $128$ grid points, respectively. The grid points are spaced uniformly in the horizontal direction and non-uniformly in the vertical direction, where they correspond to Gauss–Chebyshev quadrature nodes. A single simulation is run for $4000$ time units, where one time unit corresponds to a typical turnover time. The probability densities discussed later and presented in figure 12 were constructed from field data sampled at regular time intervals of $2$ units following a discarded transitory period of $100$ time units.

Following §§ 3.4 and 6, the joint probability density for each zone evolves according to $\partial _{t}f_{\boldsymbol{Y}|Z}=(\mathscr{L}^{\dagger }-V)f_{\boldsymbol{Y}|Z}$ , where $\mathscr{L}$ and $V$ are specific to each zone. For example, the source/sink term is

(7.11) \begin{equation} V(\boldsymbol{Y}_{\!t}|Z,t):=-\mathbb{E}[\boldsymbol{U}_{\!t}\boldsymbol{\cdot }\boldsymbol{\nabla }\log f_{\boldsymbol{X}|Z}|\boldsymbol{Y}_{\!t},Z] \end{equation}

and the drift coefficient is

(7.12) \begin{equation} \unicode{x1D63F}_{1}(\boldsymbol{Y}_{\!t}|Z,t):= \begin{pmatrix} b\\ 0 \end{pmatrix} - \mathbb{E} \left [\begin{array}{@{}c|c@{}} \partial _{X^{3}} P_{\!t} & \\ 0 & {\boldsymbol{Y}_{\!t},Z} \end{array}\right ] - \underbrace {\mathbb{E} \left [\begin{array}{@{}c|c@{}} \alpha _{1}\boldsymbol{\nabla }W_{\!t}\boldsymbol{\cdot }\boldsymbol{\nabla }\log f_{\boldsymbol{X}|Z} & \\ \alpha _{2}\boldsymbol{\nabla }B_{t}\boldsymbol{\cdot }\boldsymbol{\nabla }\log f_{\boldsymbol{X}|Z} & {\boldsymbol{Y}_{\!t},Z} \end{array}\right ]}_{\unicode{x1D63F}_{1}^{\textit{wall}}+\unicode{x1D63F}_{1}^{\textit{zone}}}, \end{equation}

where $\unicode{x1D63F}_{1}^{\textit{wall}}$ and $\unicode{x1D63F}_{1}^{\textit{zone}}$ distinguish between transport across the bottom and top boundaries ( $X^{3}\in \{0,1\}$ ), which was discussed in the example in § 7.2, from transport between zones, labelled $\unicode{x1D63F}_{1}^{\textit{zone}}$ . In each of these expressions, it is important to appreciate that the expectation is conditioned, and therefore particular, to a zone $Z$ . Integration with respect to $w$ of the governing equation for $f_{\boldsymbol{Y}|Z}$ recovers an equation for $f_{B|Z}$ , the marginal probability density for the buoyancy in each zone:

(7.13) \begin{align} \partial _{t}f_{B|Z}&=\int _{\mathbb{R}}(\mathscr{L}^{\dagger }-V)f_{\boldsymbol{Y}|Z} \,\textrm{d} w=(\overline {\mathscr{L}}^{\dagger }-\overline {V})f_{B|Z}\nonumber\\ &=-\overline {V}f_{B|Z}-\partial _{b}((\overline {\unicode{x1D63F}}^{\textit{wall}}_{1}+\overline {\unicode{x1D63F}}^{\textit{zone}}_{1})f_{B|Z})+\partial _{b}^{2}(\overline {\unicode{x1D63F}}_{2}f_{B|Z}) \end{align}

for which $\overline {\square }$ denotes the projection of the operator $\square$ defined by the integral in (7.13).

Figure 11(a) displays a snapshot of the vertical velocity $Y^{1}_{\!t}:=W_{\!t}$ and the buoyancy $Y^{2}_{\!t}:=B_{t}$ , which correspond to the fields displayed in the example derived from the Lorenz equations presented in § 2, a short time after a negative fluctuation in the buoyancy of the bottom boundary. Figure 11(b) displays the three sample distributions used to define the control volume zones $z_{1}$ (bottom boundary layer), $z_{2}$ (bulk central region) and $z_{3}$ (top boundary layer). For $z_{1}$ , $f_{\boldsymbol{X}|Z}$ is proportional to the horizontally and time-averaged diffusive buoyancy flux $-\alpha _{2}\partial _{X^{3}}\mathbb{E}[B_{t}|X^{3}]$ in the lower half of the domain and zero above. Conversely, for $z_{3}$ , $f_{\boldsymbol{X}|Z}$ is proportional to $-\alpha _{2}\partial _{X^{3}}\mathbb{E}[B_{t}|X^{3}]$ in the upper half of the domain and zero below. In the central region ( $Z=z_{2}$ ), $f_{\boldsymbol{X}|Z}$ is proportional to the mean convective buoyancy flux $\mathbb{E}[WB|X^{3}]$ . For a given zone $z$ , the constant of proportionality for each $f_{\boldsymbol{X}|Z}$ can be expressed as the total (diffusive plus convective) volume-integrated vertical buoyancy flux divided by $\mathbb{P}\{Z=z\}$ , which implies that $\mathbb{P}\{Z=z\}$ corresponds to the fraction of the volume-integrated vertical buoyancy flux over a given zone, leading to the partition of unity by $f_{\boldsymbol{X}|Z}\mathbb{P}\{Z=z\}$ illustrated in figure 11(b).

Figure 12. (a), (b) and (c) Information relating to the joint density $f_{\boldsymbol{Y}|Z}$ conditioned on zone $z_{1}$ (bottom boundary layer), $z_{2}$ (bulk central zone) and $z_{3}$ (top boundary layer), respectively. Panel (b) displays the joint density $f_{\boldsymbol{Y}|Z}$ (red contours) along with the flux $\unicode{x1D63F}^{\,\textit{zone}}_{1}f_{\boldsymbol{Y}|Z}$ (blue arrows) and isoregions of $-Vf_{\boldsymbol{Y}|Z}$ corresponding to $(-\infty ,-0.1]$ (light grey and labelled $\underline {\downarrow }, \overline {\uparrow }$ ) and $[0.1,\infty )$ (dark grey and labelled $\underline {\uparrow }, \overline {\downarrow }$ ). Panels (a) and (c) display the marginal distribution for buoyancy $f_{B|Z}$ (red bars) along with the source/sink term $-\overline {V}f_{B|Z}$ (light grey) and zonal flux term $-\partial _{b}(\overline {\unicode{x1D63F}}_{1}^{\textit{zone}}f_{B|Z})$ (blue line) in (7.13). The horizontal lines in panels (a) and (c) demarcate the range of the vertical axis in panel (b). https://www.cambridge.org/S0022112025106319/JFM-Notebooks/files/fig11-12/fig11-12.ipynb.

Figure 12 displays stationary joint and marginal distributions of $\boldsymbol{Y}_{t}$ , along with contributions to $-Vf_{\boldsymbol{Y}|Z}$ and $-\partial _{\boldsymbol{y}}\unicode{x1D63F}_{1}^{\textit{zone}}f_{\boldsymbol{Y}|Z}$ arising from transport between the control volumes. The range of buoyancy in zone $z_{1}$ (figure 12 a) is larger than it is in zone $z_{3}$ due to the stochastic boundary condition. Indeed, the marginal distribution $f_{B|Z}$ in zone $z_{3}$ is noticeably skewed; buoyancy in the top zone does not appear to fall below the fixed value of $-1/2$ imposed on the top boundary, but reaches values up to $1/2$ due to exchange with the bulk zone below. The projected convective transport $-\overline {V}f_{B|Z}$ and the diffusive transport $-\partial _{b}(\overline {\unicode{x1D63F}}^{\textit{zone}}_{1}f_{B|Z})$ for both boundary layer zones are the same order of magnitude and exhibit similar behaviour. In zone $z_{1}$ , fluid with relatively large ( $b\approx 0.5$ ) buoyancy is transported out of the zone, whereas fluid with buoyancy that is close to the volume average ( $b\approx 0$ ) is transported into the zone. In this regard, it is interesting that fluid with a buoyancy ( $b\lt -0.5$ ), which must have been created in the vicinity of a large negative fluctuation in the buoyancy imposed on the bottom boundary (see figure 11 a), is transported out of the zone. Similar remarks apply to zone $z_{3}$ , except that $-Vf_{B|Z}$ , which accounts for convective transport, tends to zero as $b$ tends to the value $-0.5$ of buoyancy imposed on the top boundary. In contrast, transport from zone $z_{3}$ into zone $z_{2}$ (the central region) is negative and large in magnitude in the vicinity of the buoyancy $b=-0.5$ of the top boundary.

Zones $z_{1}$ and $z_{3}$ interact with the central bulk zone $z_{2}$ , whose joint probability density for vertical velocity $W_{t}$ and buoyancy $B_{t}$ is displayed in figure 12(b). The vector field depicted with blue arrows represents the flux of probability density $\unicode{x1D63F}^{\textit{zone}}_{1}f_{\boldsymbol{Y}|Z}$ due to diffusion to the boundary layers, whose divergence corresponds to a sink in the budget for $f_{\boldsymbol{Y}|Z}$ . For this particular problem, the vector field resembles a saddle node, converging with respect to vertical velocity ( $w$ ) and diverging with respect to buoyancy ( $b$ ). This structure is due to the fact that the boundary layers are sinks of velocity/momentum and sources of buoyancy.

The light and dark shaded regions in figure 12(b) correspond to negative (velocity directed out of the control volume) and positive (velocity directed into the control volume) contributions to $\partial _{t} f_{\boldsymbol{Y}|Z}$ from $-Vf_{\boldsymbol{Y}|Z}$ , respectively. Positively buoyant fluid tends to be transported by positive vertical velocity from the bottom boundary layer ( $\underline {\uparrow }$ ), while negatively buoyant fluid tends to be transported by negative vertical velocity from the top boundary layer ( $\overline {\downarrow }$ ). Both of these processes correspond to positive forcing of $f_{\boldsymbol{Y}|Z}$ (i.e. $-V f_{\boldsymbol{Y}|Z}\gt 0$ ). The remaining negative regions identified in figure 12(b) typically correspond to mixed fluid, with buoyancy close to zero, transported from the bulk into the bottom or top boundary layer (denoted $\underline {\downarrow }$ and $\overline {\uparrow }$ , respectively).

7.5. Available potential energy

Available potential energy is the part of potential energy of a body of fluid that is theoretically ‘available’ for conversion into kinetic energy. The other part, known as background potential energy, cannot be converted into kinetic energy and is associated with a stable equilibrium (Lorenz Reference Lorenz1955). For example, kinetic energy cannot typically be extracted from a stably stratified environment (in which a fluid’s density decreases with height) and therefore possesses zero available potential energy.

When divided by the volume of the domain, the potential energy of the fluid modelled in § 7.3 corresponds to the expectation of the product $-B_{t}Z$ (see, for example Winters et al. Reference Winters, Lombard, Riley and D’Asaro1995):

(7.14) \begin{equation} -\mathbb{E}\left [B_{t}Z\right ]=-\int \limits _{\mathbb{R}^{n}}bz f_{\boldsymbol{Y}}(\boldsymbol{y},t)\,\textrm{d} \boldsymbol{y}, \end{equation}

where $f_{\boldsymbol{Y}}(\boldsymbol{y},t)$ is the joint density for the variables $\boldsymbol{Y}_{\!t}:=(W_{\!t},B_{t},Z)^{\top }$ from the example in § 7.3 evaluated at $\boldsymbol{y}:=(w,b,z)^{\top }$ and the minus sign accounts for relatively warm parcels of fluid having greater potential energy when they are moved downwards. The so-called ‘reference’ buoyancy profile $b^{*}:\mathbb{R}\rightarrow \mathbb{R}$ in available potential energy theory is a volume-preserving function of height that minimises potential energy (see e.g. Winters et al. Reference Winters, Lombard, Riley and D’Asaro1995). ‘Volume-preserving’ in this context means that the mapping $b^{*}$ does not change the marginal distribution of buoyancy, which, in the present context, means that $\mathbb{E}[g\circ b^{*}(Z)]=\mathbb{E}[g(B_{t})]$ for any observable $g$ .

The potential energy (7.14) is minimised by placing positively buoyant parcels at the top of the domain and negatively buoyant parcels at the bottom of the domain. The required mapping $b^{*}$ is monotonic and corresponds to

(7.15) \begin{equation} b^{*}:=F_{B}^{-1}\circ F_{Z}, \end{equation}

where $F_{B}$ and $F_{Z}$ are the marginal cumulative distribution functions for $B_{t}$ and $Z$ , respectively. In the field of optimal transport, $b^{*}$ minimises the expected squared distance between two variables and has been known for a long time (Knott & Smith Reference Knott and Smith1984; McCann Reference McCann1995). The reference state allows the available potential energy to be computed as the difference between the actual potential energy and the (minimum) potential energy associated with the reference state:

(7.16) \begin{equation} \mathbb{E}\left [b^{*}(Z)Z\right ]-\mathbb{E}\left [B_{t}Z\right ]=\int \limits _{\mathbb{R}^{n}}(b^{*}(z)-b)zf_{\boldsymbol{Y}}(\boldsymbol{y},t)\,\textrm{d} \boldsymbol{y}\geqslant 0. \end{equation}

As described by Craske (Reference Craske2021), available potential energy can be decomposed into contributions from subvolumes of a domain. Each contribution consists of ‘inner’ and ‘outer’ parts, corresponding to the potential energy that can be released within a subvolume and by combining the subvolume with its parent domain, respectively. The approach generalises to subvolumes that are specified by the conditional sample distributions $f_{\boldsymbol{X}|Z}$ described in § 6. The resulting constructions formally resemble Bayes’ theorem due to their connection with probability theory.

A key related quantity in stratified turbulence is the horizontally averaged vertical buoyancy flux $\mathbb{E}[W_{\!t}B_{t}|Z=z]$ , which is responsible for the reversible conversion of available potential energy to and from kinetic energy (see, for example, Caulfield Reference Caulfield2020). In the averaged Boussinesq equations, the horizontally averaged vertical buoyancy flux requires closure, but from the perspective of the joint density for $W_{\!t}$ and $B_{t}$ , conditional on $Z$ , it is known exactly in terms of independent and dependent variables:

(7.17) \begin{equation} \mathbb{E}[W_{\!t}B_{t}|Z=z]=\iint \limits _{\mathbb{R}^{2}}wbf_{BW|Z}(b,w|z,t)\,\textrm{d} w\,\textrm{d} b,\quad \quad \forall z: f_{Z}(z)\neq 0, \end{equation}

which suggests that the evolution equations for $f_{\boldsymbol{Y}}$ might be a useful perspective to adopt for modelling stratified turbulence and making use of available energetics in a prognostic capacity.