2.1 Evolution of passive scalar fluctuations

Let us assume that the mean scalar field vanishes. In order to derive the evolution equation for the scalar two-point correlation function, we combine (2.10) and (2.11) to obtain,

(2.12) $$\begin{eqnarray}\langle \hat{\unicode[STIX]{x1D703}}(\boldsymbol{k},t)\hat{\unicode[STIX]{x1D703}}^{\ast }(\boldsymbol{p},t)\rangle =\text{e}^{-\unicode[STIX]{x1D705}\unicode[STIX]{x1D70F}(\boldsymbol{k}^{2}+\boldsymbol{p}^{2})}\int \langle \text{e}^{-\text{i}(\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{x}-\boldsymbol{p}\boldsymbol{\cdot }\boldsymbol{y})}\rangle M(\boldsymbol{r}_{0},t_{0})\,\text{d}^{3}\boldsymbol{x}\,\text{d}^{3}\boldsymbol{y},\end{eqnarray}$$

where we have defined $\boldsymbol{r}_{0}=\boldsymbol{x}_{0}-\boldsymbol{y}_{0}$ . The statistical homogeneity of the field allows us to write the two-point scalar correlator in Fourier space as,

(2.13) $$\begin{eqnarray}\langle \hat{\unicode[STIX]{x1D703}}(\boldsymbol{k},t)\hat{\unicode[STIX]{x1D703}}^{\ast }(\boldsymbol{p},t)\rangle =(2\unicode[STIX]{x03C0})^{3}\unicode[STIX]{x1D6FF}^{3}(\boldsymbol{k}-\boldsymbol{p})\hat{M}(\boldsymbol{k},t).\end{eqnarray}$$

We now follow a method similar to that used in BS15. We can transform $(\boldsymbol{x},\boldsymbol{y})$ in (2.12) to $(\boldsymbol{x}_{0},\boldsymbol{y}_{0})$ using (2.9). The Jacobian of this transformation is unity due to incompressibility of the flow. Then we write $\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{x}_{0}-\boldsymbol{p}\boldsymbol{\cdot }\boldsymbol{y}_{0}=\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{r}_{0}+\boldsymbol{y}_{0}\boldsymbol{\cdot }(\boldsymbol{k}-\boldsymbol{p})$ in (2.12), transform from $(\boldsymbol{x}_{0},\boldsymbol{y}_{0})$ to a new set of variables $(\boldsymbol{r}_{0},\boldsymbol{y}_{0}^{\prime }=\boldsymbol{y}_{0})$ and integrate over $\boldsymbol{y}_{0}^{\prime }$ . This leads to a delta function in $(\boldsymbol{k}-\boldsymbol{p})$ and (2.12) becomes,

(2.14) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\hat{M}(\boldsymbol{k},t)=\text{e}^{-2\unicode[STIX]{x1D705}\unicode[STIX]{x1D70F}\boldsymbol{k}^{2}}\displaystyle \int \langle R\rangle M(\boldsymbol{r}_{0},t_{0})\text{e}^{-\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{r}_{0}}\,\text{d}^{3}\boldsymbol{r}_{0},\\ \langle R\rangle =\langle \text{e}^{-\text{i}\unicode[STIX]{x1D70F}(\boldsymbol{a}\boldsymbol{\cdot }\boldsymbol{k})(\sin \bar{A}-\sin \bar{B})}\rangle ,\end{array}\right\} & & \displaystyle\end{eqnarray}$$

where $\bar{A}=(\boldsymbol{x}_{0}\boldsymbol{\cdot }\boldsymbol{q}+\unicode[STIX]{x1D713})$ and $\bar{B}=(\boldsymbol{y}_{0}\boldsymbol{\cdot }\boldsymbol{q}+\unicode[STIX]{x1D713})$ . Note that $\langle R\rangle$ is expected to be only a function of $\boldsymbol{r}_{0}$ , due to statistical homogeneity of the renewing flow, as we also see explicitly in § 3.

We are interested in the limit of large Péclet number, $R_{\unicode[STIX]{x1D705}}\gg 1$ . In this limit, we can expand the exponential in the diffusive Green’s function and consider only leading-order term in $\unicode[STIX]{x1D705}$ . Then we have

(2.15) $$\begin{eqnarray}\displaystyle \hat{M}(\boldsymbol{k},t) & = & \displaystyle (1-2\unicode[STIX]{x1D705}\unicode[STIX]{x1D70F}\boldsymbol{k}^{2})\int \langle R\rangle M(\boldsymbol{r}_{0},t_{0})\text{e}^{-\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{r}_{0}}\,\text{d}^{3}\boldsymbol{r}_{0}\nonumber\\ \displaystyle & = & \displaystyle \int \langle R\rangle \hat{M}(\boldsymbol{p},t_{0})\text{e}^{\text{i}(\boldsymbol{p}-\boldsymbol{k})\boldsymbol{\cdot }\boldsymbol{r}_{0}}\,\text{d}^{3}\boldsymbol{r}_{0}\frac{\text{d}^{3}\boldsymbol{p}}{(2\unicode[STIX]{x03C0})^{3}}-2\unicode[STIX]{x1D705}\unicode[STIX]{x1D70F}\boldsymbol{k}^{2}\hat{M}(\boldsymbol{k},t_{0}).\end{eqnarray}$$

In the diffusive term, we consider only the $\unicode[STIX]{x1D70F}$ independent term of $\langle R\rangle$ to multiply with $2\unicode[STIX]{x1D705}\unicode[STIX]{x1D70F}\boldsymbol{k}^{2}$ , as all the other terms will be of the order of $\unicode[STIX]{x1D705}\unicode[STIX]{x1D70F}^{2}$ or higher. (Specifically, in the independent $\unicode[STIX]{x1D705}\rightarrow 0$ limit, we neglect terms proportional to $\unicode[STIX]{x1D705}\unicode[STIX]{x1D70F}^{2}$ compared to those $\propto \unicode[STIX]{x1D70F}^{3}$ and $\unicode[STIX]{x1D70F}^{4}$ .) We will use (2.15) when we generalize the Kraichnan evolution equation for the passive scalar to incorporate finite $\unicode[STIX]{x1D70F}$ effects.

We can also study the evolution of scalar fluctuations in terms of the real space correlation function. For this we take the inverse Fourier transform of $\hat{M}(\boldsymbol{k},t)$ and get

(2.16) $$\begin{eqnarray}\displaystyle M(\boldsymbol{r},t) & = & \displaystyle \int (1-2\unicode[STIX]{x1D705}\unicode[STIX]{x1D70F}\boldsymbol{k}^{2})\langle R\rangle M(\boldsymbol{r}_{0},t_{0})\text{e}^{\text{i}\boldsymbol{k}\boldsymbol{\cdot }(\boldsymbol{r}-\boldsymbol{r}_{0})}\,\text{d}^{3}\boldsymbol{r}_{0}\frac{\text{d}^{3}\boldsymbol{k}}{(2\unicode[STIX]{x03C0})^{3}}\nonumber\\ \displaystyle & = & \displaystyle \int \left\langle R\right\rangle M(\boldsymbol{r}_{0},t_{0})\text{e}^{\text{i}\boldsymbol{k}\boldsymbol{\cdot }(\boldsymbol{r}-\boldsymbol{r}_{0})}\,\text{d}^{3}\boldsymbol{r}_{0}\frac{\text{d}^{3}\boldsymbol{k}}{(2\unicode[STIX]{x03C0})^{3}}+2\unicode[STIX]{x1D705}\unicode[STIX]{x1D70F}\unicode[STIX]{x1D6FB}^{2}M(\boldsymbol{r},t).\end{eqnarray}$$

Here we again expanded the exponential in the diffusive Green’s function relevant in the independent small $\unicode[STIX]{x1D705}$ (or $R_{\unicode[STIX]{x1D705}}\gg 1$ ) limit. And in the diffusive term of (2.16) we consider only the $\unicode[STIX]{x1D70F}$ independent term in $\langle R\rangle$ as above. Then we write $(-\boldsymbol{k}^{2})(\text{e}^{\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{r}})$ as $\unicode[STIX]{x1D6FB}^{2}\text{e}^{\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{r}}$ and can take $\unicode[STIX]{x1D6FB}^{2}$ out of the integral.

Both the Fourier space and real space approaches are useful in different contexts. We now turn to the evaluation of $\langle R\rangle$ .

3.1 Kraichnan passive scalar equation from terms up to order $\unicode[STIX]{x1D70F}^{2}$

We first consider terms in (3.1) up to the order $\unicode[STIX]{x1D70F}^{2}$ and average over $\unicode[STIX]{x1D713}$ , $\hat{\boldsymbol{a}}$ and $\hat{\boldsymbol{q}}$ . To begin with we note that $\langle \unicode[STIX]{x1D70E}\rangle =0$ as it is proportional to $\cos (\unicode[STIX]{x1D713}+\cdots \,)$ and the integral over $\unicode[STIX]{x1D713}$ vanishes. In the third term of (3.1), we have

(3.3) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70E}^{2} & = & \displaystyle 4(\boldsymbol{a}\boldsymbol{\cdot }\boldsymbol{k})^{2}\sin ^{2}\left(\frac{\boldsymbol{q}\boldsymbol{\cdot }\boldsymbol{r}_{0}}{2}\right)\cos ^{2}(\unicode[STIX]{x1D713}+\boldsymbol{q}\boldsymbol{\cdot }\boldsymbol{R}_{0})\nonumber\\ \displaystyle & = & \displaystyle (\boldsymbol{a}\boldsymbol{\cdot }\boldsymbol{k})^{2}(1-\cos \boldsymbol{q}\boldsymbol{\cdot }\boldsymbol{r}_{0})(1+\cos (2\unicode[STIX]{x1D713}+2\boldsymbol{q}\boldsymbol{\cdot }\boldsymbol{R}_{0}))\end{eqnarray}$$

and on averaging over $\unicode[STIX]{x1D713}$ , the term proportional to $\cos (2\unicode[STIX]{x1D713}+\cdots \,)$ vanishes. Therefore we have for this term

(3.4) $$\begin{eqnarray}\frac{\unicode[STIX]{x1D70F}^{2}\langle \unicode[STIX]{x1D70E}^{2}\rangle }{2}=\frac{\unicode[STIX]{x1D70F}^{2}}{2}k_{i}k_{j}\langle a_{i}a_{j}(1-\cos \boldsymbol{q}\boldsymbol{\cdot }\boldsymbol{r}_{0})\rangle =2\unicode[STIX]{x1D70F}k_{i}k_{j}[T_{ij}(0)-T_{ij}(r_{0})],\end{eqnarray}$$

where we have defined the turbulent diffusion tensor, in terms of the two-point velocity correlator, as

(3.5) $$\begin{eqnarray}T_{ij}(\boldsymbol{r}_{0})=\frac{\unicode[STIX]{x1D70F}}{2}\langle u_{i}(\boldsymbol{x}_{0})u_{j}(\boldsymbol{y}_{0})\rangle =\frac{\unicode[STIX]{x1D70F}}{2}\langle a_{i}a_{j}\sin (\bar{A})\sin (\bar{B})\rangle =\frac{\unicode[STIX]{x1D70F}}{4}\langle a_{i}a_{j}\cos (\boldsymbol{q}\boldsymbol{\cdot }\boldsymbol{r}_{0})\rangle .\end{eqnarray}$$

We have included a factor of $\unicode[STIX]{x1D70F}$ in the above definition of $T_{ij}$ , since in the limit of $\unicode[STIX]{x1D70F}\rightarrow 0$ (as for delta correlated in time flow) one still wants the turbulent diffusion to have a finite limit. Using (3.4), we then have up to order $\unicode[STIX]{x1D70F}^{2}$ ,

(3.6) $$\begin{eqnarray}\langle R\rangle =1-2\unicode[STIX]{x1D70F}k_{i}k_{j}[T_{ij}(0)-T_{ij}(r_{0})].\end{eqnarray}$$

We substitute (3.6) in (2.16) for the passive scalar correlation function. For the term multiplying unity, the integration in (2.16) over $\boldsymbol{k}$ gives a delta function $\unicode[STIX]{x1D6FF}^{3}(\boldsymbol{r}-\boldsymbol{r}_{0})$ , and thus the integration over $\boldsymbol{r}_{0}$ simply replaces it with $\boldsymbol{r}$ in $M(\boldsymbol{r}_{0},t_{0})$ . For the second term containing $k_{i}$ , we can first write it as a derivative with respect to $r_{i}$ , pull it out of the integral and then again integrate over $\boldsymbol{k}$ and $\boldsymbol{r}_{0}$ as above. We then get from (3.6) and (2.16),

(3.7) $$\begin{eqnarray}\displaystyle M(r,t) & = & \displaystyle \int (1-2\unicode[STIX]{x1D70F}k_{i}k_{j}[T_{ij}(0)-T_{ij}(r_{0})])M(\boldsymbol{r}_{0},t_{0})\text{e}^{\text{i}\boldsymbol{k}\boldsymbol{\cdot }(\boldsymbol{r}-\boldsymbol{r}_{0})}\,\text{d}^{3}\boldsymbol{r}_{0}\frac{\text{d}^{3}\boldsymbol{k}}{(2\unicode[STIX]{x03C0})^{3}}+2\unicode[STIX]{x1D705}\unicode[STIX]{x1D70F}\unicode[STIX]{x1D6FB}^{2}M(\boldsymbol{r},t)\nonumber\\ \displaystyle & = & \displaystyle M(r,t_{0})+2\unicode[STIX]{x1D70F}\unicode[STIX]{x2202}_{i}\unicode[STIX]{x2202}_{j}\int [T_{ij}(0)-T_{ij}(\boldsymbol{r}_{0})]M(\boldsymbol{r}_{0},t_{0})\text{e}^{\text{i}\boldsymbol{k}\boldsymbol{\cdot }(\boldsymbol{r}-\boldsymbol{r}_{0})}\,\text{d}^{3}\boldsymbol{r}_{0}\frac{\text{d}^{3}\boldsymbol{k}}{(2\unicode[STIX]{x03C0})^{3}}\nonumber\\ \displaystyle & & \displaystyle +\,2\unicode[STIX]{x1D705}\unicode[STIX]{x1D70F}\unicode[STIX]{x1D6FB}^{2}M(\boldsymbol{r},t)\nonumber\\ \displaystyle & = & \displaystyle M(r,t_{0})+2\unicode[STIX]{x1D70F}\unicode[STIX]{x2202}_{i}\unicode[STIX]{x2202}_{j}[(T_{ij}(0)-T_{ij}(\boldsymbol{r}))M(r,t_{0})]+2\unicode[STIX]{x1D705}\unicode[STIX]{x1D70F}\unicode[STIX]{x1D6FB}^{2}M(r,t).\end{eqnarray}$$

Here $\unicode[STIX]{x2202}_{i}=\unicode[STIX]{x2202}/\unicode[STIX]{x2202}r_{i}$ and we have used the fact that for statistically isotropic scalar correlations $M(\boldsymbol{r},t)=M(r,t)$ . Note that for a divergence free velocity, $\unicode[STIX]{x2202}_{i}(T_{ij})=0$ and so the spatial derivatives can be taken through the $T_{ij}$ terms (3.7) to directly act on $M(r,t)$ .

A differential equation for the time evolution of $M(r,t)$ can be derived from (3.7), by dividing throughout by $\unicode[STIX]{x1D70F}$ , and noting that in the limit $\unicode[STIX]{x1D70F}\rightarrow 0$ , $(M(r,t)-M(r,t_{0}))/\unicode[STIX]{x1D70F}=\unicode[STIX]{x2202}M/\unicode[STIX]{x2202}t$ . We obtain

(3.8) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}M(r,t)}{\unicode[STIX]{x2202}t}=[T_{ij}(0)-T_{ij}(\boldsymbol{r})]\frac{\unicode[STIX]{x2202}^{2}M}{\unicode[STIX]{x2202}r^{i}\unicode[STIX]{x2202}r^{j}}+2\unicode[STIX]{x1D705}\unicode[STIX]{x1D6FB}^{2}M(r,t).\end{eqnarray}$$

This can be further simplified by first noting that

(3.9) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}^{2}M}{\unicode[STIX]{x2202}r^{i}\unicode[STIX]{x2202}r^{j}}=(\unicode[STIX]{x1D6FF}_{ij}-\hat{r}_{i}\hat{r}_{j})\frac{1}{r}\frac{\unicode[STIX]{x2202}M}{\unicode[STIX]{x2202}r}+\hat{r}_{i}\hat{r}_{j}\frac{\unicode[STIX]{x2202}^{2}M}{\unicode[STIX]{x2202}r^{2}}.\end{eqnarray}$$

Also for a statistically homogeneous, isotropic and non-helical velocity field, the correlation function

(3.10) $$\begin{eqnarray}T_{ij}=(\unicode[STIX]{x1D6FF}_{ij}-\hat{r}_{i}\hat{r}_{j})T_{N}(r,t)+\hat{r}_{i}\hat{r}_{j}T_{L}(r,t).\end{eqnarray}$$

Here $\hat{r_{i}}=r_{i}/r$ , and $T_{ij}(0)=\unicode[STIX]{x1D6FF}_{ij}T_{L}(0)$ , with

(3.11a,b ) $$\begin{eqnarray}T_{L}(r,t)=\hat{r}_{i}\hat{r}_{j}T_{ij}\quad \text{and}\quad T_{N}(r,t)=\frac{1}{2r}\frac{\unicode[STIX]{x2202}(r^{2}T_{L})}{\unicode[STIX]{x2202}r}\end{eqnarray}$$

respectively being, the longitudinal and transversal correlation functions of the velocity field. Using (3.9) and (3.10) in (3.8), we then obtain

(3.12) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}M(r,t)}{\unicode[STIX]{x2202}t}=\frac{2}{r^{2}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}\left[r^{2}(\unicode[STIX]{x1D705}+T_{L}(0)-T_{L}(r))\frac{\unicode[STIX]{x2202}M}{\unicode[STIX]{x2202}r}\right].\end{eqnarray}$$

This is the Kraichnan equation for the evolution of the passive scalar fluctuations in real space. It has a very simple interpretation; at the lowest order in $\unicode[STIX]{x1D70F}$ the effect of the random velocity is to add scale-dependent turbulent diffusivity $[T_{L}(0)-T_{L}(r)]$ to the microscopic diffusivity $\unicode[STIX]{x1D705}$ . Although we derived (3.12) in the context of renewing flows, the only reflection of the velocity field is in $T_{ij}$ . In fact the same equation results when one uses a delta correlated in time velocity field with the spatial part of the correlation being given by $T_{ij}$ (see also Zeldovich et al. Reference Zeldovich, Molchanov, Ruzmaikin and Sokoloff1988).

3.2 Kraichnan equation in Fourier space from up to $\unicode[STIX]{x1D70F}^{2}$ terms

In order to derive the passive scalar spectrum (both for the steady state and the decaying case), in a transparent manner, it is more useful to directly work in Fourier space. In Fourier space, equation (3.12) becomes an integro–differential equation. However, the Batchelor regime is supposed to occur only for wavenumbers $k$ , much larger than the viscous cutoff wavenumber (in the case of single-scale renewing flow this corresponds to $q$ ) and much smaller than the wavenumber where scalar diffusion is important; or in the range $q\ll k\ll k_{\unicode[STIX]{x1D705}}$ . For $k\gg q$ , one has $qr\ll 1$ and so we can make the small $qr$ approximation in evaluating $T_{ij}(r)$ .

For this we expand $\cos (\boldsymbol{q}\boldsymbol{\cdot }\boldsymbol{r})$ in (3.5) for small $qr$ , use (2.3a,b ) to write $a_{i}$ in terms of $A_{i}$ . We also use $\langle A_{i}A_{j}\rangle =A^{2}\unicode[STIX]{x1D6FF}_{ij}/3$ , $\langle q_{i}q_{j}\rangle =q^{2}\unicode[STIX]{x1D6FF}_{ij}/3$ and $\langle q_{i}q_{j}q_{l}q_{m}\rangle =(4/15)[\unicode[STIX]{x1D6FF}_{ij}\unicode[STIX]{x1D6FF}_{lm}+\unicode[STIX]{x1D6FF}_{il}\unicode[STIX]{x1D6FF}_{jm}+\unicode[STIX]{x1D6FF}_{im}\unicode[STIX]{x1D6FF}_{lj}]$ to get $T_{ij}(0)=\unicode[STIX]{x1D6FF}_{ij}T_{L}(0)$ , with $T_{L}(0)=\unicode[STIX]{x1D705}_{t}=A^{2}\unicode[STIX]{x1D70F}/18$ and

(3.13) $$\begin{eqnarray}[T_{ij}(0)-T_{ij}(\boldsymbol{r})]=\frac{\unicode[STIX]{x1D70F}}{4}\langle a_{i}a_{j}[1-\cos (\boldsymbol{q}\boldsymbol{\cdot }\boldsymbol{r})]\rangle ={\mathcal{A}}[2r^{2}\unicode[STIX]{x1D6FF}_{ij}-r_{i}r_{j}]\end{eqnarray}$$

with ${\mathcal{A}}=A^{2}q^{2}\unicode[STIX]{x1D70F}/90=\unicode[STIX]{x1D705}_{t}q^{2}/5$ . From (3.13) we then have for small $qr$ , and up to order $\unicode[STIX]{x1D70F}^{2}$ ,

(3.14) $$\begin{eqnarray}\langle R\rangle =1-\unicode[STIX]{x1D70F}k_{i}k_{j}{\mathcal{A}}(2r_{0}^{2}\unicode[STIX]{x1D6FF}_{ij}-r_{0i}r_{0j})\equiv 1-R_{2}.\end{eqnarray}$$

We use this in (2.15) to derive the Fourier space Kraichnan equation and passive scalar spectrum in the viscous–convective regime.

We substitute (3.14) in (2.15) for the passive scalar spectrum, in the limit $k\gg q$ , and up to order $\unicode[STIX]{x1D70F}^{2}$ ,

(3.15) $$\begin{eqnarray}\hat{M}(\boldsymbol{k},t)=\int (1-R_{2})\hat{M}(\boldsymbol{p},t_{0})\text{e}^{\text{i}(\boldsymbol{p}-\boldsymbol{k})\boldsymbol{\cdot }\boldsymbol{r}_{0}}\,\text{d}^{3}\boldsymbol{r}_{0}\frac{\text{d}^{3}\boldsymbol{p}}{(2\unicode[STIX]{x03C0})^{3}}-2\unicode[STIX]{x1D705}\unicode[STIX]{x1D70F}\boldsymbol{k}^{2}\hat{M}(\boldsymbol{k},t_{0}).\end{eqnarray}$$

In the first term proportional to unity one can straightaway do the integration over $\boldsymbol{r}_{0}$ to get a delta function $\unicode[STIX]{x1D6FF}^{3}(\boldsymbol{p}-\boldsymbol{k})$ and then the $\boldsymbol{p}$ integral becomes trivial and gives $\hat{M}(\boldsymbol{k},t_{0})$ . In the second term, proportional to $R_{2}$ , we convert each factor or $r_{0i}$ to a derivative with respect to $k_{i}$ and integrate by parts. The integral over $\boldsymbol{r}_{0}$ again gives $\unicode[STIX]{x1D6FF}^{3}(\boldsymbol{p}-\boldsymbol{k})$ and again the integral over $\boldsymbol{p}$ can be done trivially. Then

(3.16) $$\begin{eqnarray}\displaystyle -\int R_{2}\hat{M}(\boldsymbol{p},t_{0})\text{e}^{\text{i}(\boldsymbol{p}-\boldsymbol{k})\boldsymbol{\cdot }\boldsymbol{r}_{0}}\,\text{d}^{3}\boldsymbol{r}_{0}\frac{\text{d}^{3}\boldsymbol{p}}{(2\unicode[STIX]{x03C0})^{3}} & = & \displaystyle {\mathcal{A}}\unicode[STIX]{x1D70F}k_{i}k_{j}\left[2\unicode[STIX]{x1D6FF}_{ij}\frac{\unicode[STIX]{x2202}^{2}\hat{M}}{\unicode[STIX]{x2202}k_{l}\unicode[STIX]{x2202}k_{l}}-\frac{\unicode[STIX]{x2202}^{2}\hat{M}}{\unicode[STIX]{x2202}k_{i}\unicode[STIX]{x2202}k_{j}}\right]\nonumber\\ \displaystyle & = & \displaystyle {\mathcal{A}}\unicode[STIX]{x1D70F}\left[k^{2}\frac{\unicode[STIX]{x2202}^{2}\hat{M}}{\unicode[STIX]{x2202}k^{2}}+4k\frac{\unicode[STIX]{x2202}\hat{M}}{\unicode[STIX]{x2202}k}\right].\end{eqnarray}$$

Substituting (3.16) into (3.15), dividing throughout by $\unicode[STIX]{x1D70F}$ and taking the $\unicode[STIX]{x1D70F}\rightarrow 0$ limit we then obtain,

(3.17) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\hat{M}(k,t)}{\unicode[STIX]{x2202}t}=\frac{\unicode[STIX]{x1D705}_{t}q^{2}}{5}\left[k^{2}\frac{\unicode[STIX]{x2202}^{2}\hat{M}}{\unicode[STIX]{x2202}k^{2}}+4k\frac{\unicode[STIX]{x2202}\hat{M}}{\unicode[STIX]{x2202}k}\right]-2\unicode[STIX]{x1D705}k^{2}\hat{M}(k,t).\end{eqnarray}$$

This is the Kraichnan (Reference Kraichnan1968) evolution equation for the passive scalar in Fourier space applicable to the viscous–convective regime, with the first term on the right-hand side of (3.17) the same as $T(k,t)/4\unicode[STIX]{x03C0}k^{2}$ given in equation (3.6) of K68. We can also rewrite (3.17) in terms of an equation for the 1-D scalar spectrum $E_{\unicode[STIX]{x1D703}}(k)=4\unicode[STIX]{x03C0}k^{2}\tilde{M}(k)$ . We find,

(3.18) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}E_{\unicode[STIX]{x1D703}}(k,t)}{\unicode[STIX]{x2202}t}=\frac{\unicode[STIX]{x1D705}_{t}q^{2}}{5}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}k}\left[k^{4}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}k}\left(\frac{E_{\unicode[STIX]{x1D703}}}{k^{2}}\right)\right]-2\unicode[STIX]{x1D705}k^{2}E_{\unicode[STIX]{x1D703}}(k,t).\end{eqnarray}$$

We look at its consequences for the passive scalar spectrum after also working out the finite $\unicode[STIX]{x1D70F}$ corrections to (3.17).

3.3 Finite $\unicode[STIX]{x1D70F}$ correction to Kraichnan equation from up to $\unicode[STIX]{x1D70F}^{4}$ terms

Now let us consider the higher-order terms up to $\unicode[STIX]{x1D70F}^{4}$ . Firstly, the $\unicode[STIX]{x1D70F}^{3}$ term in (3.1) is proportional to $\langle \cos ^{3}(\unicode[STIX]{x1D713}+\boldsymbol{q}\boldsymbol{\cdot }\boldsymbol{R}_{0})\rangle$ and always contains a $n\unicode[STIX]{x1D713}$ in the argument of a cosine. Thus it goes to zero on averaging over $\unicode[STIX]{x1D713}$ . Now consider the term in (3.1) of order $\unicode[STIX]{x1D70F}^{4}$ . This is given by

(3.19) $$\begin{eqnarray}\displaystyle R_{4}=\frac{\langle \unicode[STIX]{x1D70F}^{4}\unicode[STIX]{x1D70E}^{4}\rangle }{24} & = & \displaystyle \frac{\unicode[STIX]{x1D70F}^{4}}{24}\left\langle (\boldsymbol{a}\boldsymbol{\cdot }\boldsymbol{k})^{4}\sin ^{4}\left(\frac{\boldsymbol{q}\boldsymbol{\cdot }\boldsymbol{r}_{0}}{2}\right)\cos ^{4}\left(\unicode[STIX]{x1D713}+\boldsymbol{q}\boldsymbol{\cdot }\boldsymbol{R}_{0}\right)\right\rangle \nonumber\\ \displaystyle & = & \displaystyle \frac{\unicode[STIX]{x1D70F}^{4}}{16}k_{m}k_{n}k_{r}k_{s}\left\langle a_{m}a_{n}a_{r}a_{s}\left(\frac{3}{2}-2\cos (\boldsymbol{q}\boldsymbol{\cdot }\boldsymbol{r}_{0})+\frac{\cos (2\boldsymbol{q}\boldsymbol{\cdot }\boldsymbol{r}_{0})}{2}\right)\right\rangle \nonumber\\ \displaystyle & = & \displaystyle \unicode[STIX]{x1D70F}^{2}k_{m}k_{n}k_{r}k_{s}\left[\frac{T_{mnrs}^{x^{2}y^{2}}}{4}-\frac{T_{mnrs}^{x^{3}y}}{3}+\frac{T_{mnrs}^{x^{4}}}{12}\right],\end{eqnarray}$$

where we have carried out the averaging over $\unicode[STIX]{x1D713}$ . Further, in the last line of (3.19), we have expressed $R_{4}$ in terms of the two-point fourth-order velocity correlators as defined by BS15. Note that three such velocity correlators can be defined as,

(3.20) $$\begin{eqnarray}\displaystyle & \displaystyle T_{mnrs}^{x^{2}y^{2}}=\unicode[STIX]{x1D70F}^{2}\langle u_{m}(\boldsymbol{x})u_{n}(\boldsymbol{y})u_{r}(\boldsymbol{x})u_{s}(\boldsymbol{y})\rangle =\frac{\unicode[STIX]{x1D70F}^{2}}{4}\left\langle a_{m}a_{n}a_{r}a_{s}\left(1+\frac{\cos (2\boldsymbol{q}\boldsymbol{\cdot }\boldsymbol{r}_{0})}{2}\right)\right\rangle , & \displaystyle\end{eqnarray}$$

(3.21) $$\begin{eqnarray}\displaystyle & \displaystyle T_{mnrs}^{x^{3}y}=\unicode[STIX]{x1D70F}^{2}\langle u_{m}(\boldsymbol{x})u_{n}(\boldsymbol{x})u_{r}(\boldsymbol{x})u_{s}(\boldsymbol{y})\rangle =\frac{3\unicode[STIX]{x1D70F}^{2}}{8}\left\langle a_{m}a_{n}a_{r}a_{s}\cos (\boldsymbol{q}\boldsymbol{\cdot }\boldsymbol{r}_{0})\right\rangle , & \displaystyle\end{eqnarray}$$

(3.22) $$\begin{eqnarray}\displaystyle & \displaystyle T_{mnrs}^{x^{4}}=\unicode[STIX]{x1D70F}^{2}\langle u_{m}(\boldsymbol{x})u_{n}(\boldsymbol{x})u_{r}(\boldsymbol{x})u_{s}(\boldsymbol{x})\rangle =\frac{3\unicode[STIX]{x1D70F}^{2}}{8}\left\langle a_{m}a_{n}a_{r}a_{s}\right\rangle , & \displaystyle\end{eqnarray}$$

where in the last parts of each equation the $\unicode[STIX]{x1D713}$ averaged expressions are given. Further, we have also multiplied the fourth-order velocity correlators by $\unicode[STIX]{x1D70F}^{2}$ , as we envisage that $T_{ijkl}$ will be finite even in the $\unicode[STIX]{x1D70F}\rightarrow 0$ limit. Essentially, they behave like products of turbulent diffusion. Interestingly, the renewing flow is not Gaussian random, and therefore higher-order correlators of $\boldsymbol{u}$ are not just the product of two-point correlators.

We can now add the $R_{4}$ contribution to $\langle R\rangle$ , substitute in to (2.16) for the scalar correlator. Again each $k_{i}$ can be converted to a derivative over $r_{i}$ and pulled out of the integral in (2.16). The resulting integrals over $\boldsymbol{k}$ and $\boldsymbol{r}_{0}$ , including the $R_{4}$ term, can then be carried out as earlier in deriving (3.7). We then divide all contributions by $\unicode[STIX]{x1D70F}$ and take the limit as $\unicode[STIX]{x1D70F}\rightarrow 0$ . For the $\unicode[STIX]{x1D70F}^{4}$ term on division by $\unicode[STIX]{x1D70F}$ , a factor of $\unicode[STIX]{x1D70F}$ remains as a small effective finite time parameter. We then get for the extended Kraichnan evolution equation,

(3.23) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}M(r,t)}{\unicode[STIX]{x2202}t}=[T_{ij}(0)-T_{ij}(\boldsymbol{r})]\frac{\unicode[STIX]{x2202}^{2}M}{\unicode[STIX]{x2202}r^{i}\unicode[STIX]{x2202}r^{j}}+2\unicode[STIX]{x1D705}\unicode[STIX]{x1D70F}\unicode[STIX]{x1D6FB}^{2}M(r,t)+\unicode[STIX]{x1D70F}\left[\frac{T_{mnrs}^{x^{2}y^{2}}}{4}-\frac{T_{mnrs}^{x^{3}y}}{3}+\frac{T_{mnrs}^{x^{4}}}{12}\right]M_{,mnrs}. & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

The first line in (3.23) has the terms which give the standard Kraichnan equation (3.8) in real space, whereas the second line contains the finite correlation time corrections. This generalized Kraichnan equation can also be further simplified by using an explicit form for the fourth-order correlators and the fourth derivative of $M$ , following very similar algebra as in BS15. It can then be used to examine the finite correlation time modification to the passive scalar spectrum both in the steady state and decaying cases. However this is again derived in a more transparent manner in Fourier space, to which we now turn.

3.4 Generalized Kraichnan equation in Fourier space

For calculating the finite $\unicode[STIX]{x1D70F}$ corrections to the Fourier space K68 equation (3.17) and the passive scalar spectrum in the Batchelor regime, we need only the small $qr_{0}$ expansion to the $R_{4}$ term. This is because, as we discussed earlier, an extended Batchelor regime occurs for large Schmidt numbers, in the range of wavenumbers when $k_{\unicode[STIX]{x1D705}}\gg k\gg q$ , from (3.19), in the $qr_{0}\ll 1$ limit we have

(3.24) $$\begin{eqnarray}R_{4}=\frac{\unicode[STIX]{x1D70F}^{4}}{64}\langle (\boldsymbol{a}\boldsymbol{\cdot }\boldsymbol{k})^{4}(\boldsymbol{q}\boldsymbol{\cdot }\boldsymbol{r}_{0})^{4}\rangle .\end{eqnarray}$$

Substituting (3.24) into (2.15), the additional $\unicode[STIX]{x1D70F}^{4}$ contribution of $R_{4}$ to $M(\boldsymbol{k},t)$ is given by

(3.25) $$\begin{eqnarray}\displaystyle \hat{M}_{4}(k,t) & = & \displaystyle \frac{\unicode[STIX]{x1D70F}^{4}}{64}\langle (\boldsymbol{a}\boldsymbol{\cdot }\boldsymbol{k})^{4}(q_{a}q_{b}q_{c}q_{d})\rangle \int r_{0a}r_{0b}r_{0c}r_{0d}\hat{M}(\boldsymbol{p},t_{0})\text{e}^{\text{i}(\boldsymbol{p}-\boldsymbol{k})\boldsymbol{\cdot }\boldsymbol{r}_{0}}\,\text{d}^{3}\boldsymbol{r}_{0}\frac{\text{d}^{3}\boldsymbol{p}}{(2\unicode[STIX]{x03C0})^{3}}\nonumber\\ \displaystyle & = & \displaystyle \frac{\unicode[STIX]{x1D70F}^{4}}{64}\langle (\boldsymbol{a}\boldsymbol{\cdot }\boldsymbol{k})^{4}(q_{a}q_{b}q_{c}q_{d})\rangle \frac{\unicode[STIX]{x2202}^{4}\hat{M}(k,t_{0})}{\unicode[STIX]{x2202}k_{a}\unicode[STIX]{x2202}k_{b}\unicode[STIX]{x2202}k_{c}\unicode[STIX]{x2202}k_{d}}.\end{eqnarray}$$

In the second step of (3.25), we have converted each $r_{0a}$ factor to a derivative of the exponential factor with respect to $p_{a}$ and integrated by parts to transfer these to derivatives of $\hat{M}$ . The integral over $\boldsymbol{r}_{0}$ then gives a delta function $\unicode[STIX]{x1D6FF}^{3}(\boldsymbol{k}-\boldsymbol{p})$ and the integral over $\boldsymbol{p}$ then simply replaces $\boldsymbol{p}$ everywhere by $\boldsymbol{k}$ .

We can again write for each $\boldsymbol{a}\boldsymbol{\cdot }\boldsymbol{k}=a_{i}k_{i}=[\unicode[STIX]{x1D6FF}_{ij}-\hat{q}_{i}\hat{q}_{j}]A_{j}k_{i}$ , and average over the independent $A_{i}$ . This leaves averaging over $\boldsymbol{q}$ . If one does this $q_{i}$ averaging naively, we would need to deal with averaging an eight-point $q_{i}$ correlator, which involves 105 terms. To avoid this we first evaluate the fourth derivative of $\hat{M}$ explicitly, which gives

(3.26) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}^{4}\hat{M}}{\unicode[STIX]{x2202}k_{a}\unicode[STIX]{x2202}k_{b}\unicode[STIX]{x2202}k_{c}\unicode[STIX]{x2202}k_{d}} & = & \displaystyle \frac{1}{k}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}k}\left(\frac{1}{k}\frac{\unicode[STIX]{x2202}M}{\unicode[STIX]{x2202}k}\right)\unicode[STIX]{x1D6FF}_{(ab}\unicode[STIX]{x1D6FF}_{cd )}+\frac{1}{k}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}k}\left[\frac{1}{k}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}k}\left(\frac{1}{k}\frac{\unicode[STIX]{x2202}M}{\unicode[STIX]{x2202}k}\right)\right]\unicode[STIX]{x1D6FF}_{(ab}k_{c}k_{d )}\nonumber\\ \displaystyle & & \displaystyle +\,\frac{1}{k}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}k}\left(\frac{1}{k}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}k}\left[\frac{1}{k}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}k}\left(\frac{1}{k}\frac{\unicode[STIX]{x2202}M}{\unicode[STIX]{x2202}k}\right)\right]\right)k_{a}k_{b}k_{c}k_{d}.\end{eqnarray}$$

We have used the brackets $()$ in the subscripts of two second-rank tensors to denote addition of all terms from different permutations of the four indices considered in pairs. This brings out components of $k_{a}$ , which can combine with corresponding values of $q_{a}$ , and one only gets up to eighth power of the dot product $(\boldsymbol{q}\boldsymbol{\cdot }\boldsymbol{k})=qk\unicode[STIX]{x1D707}$ , where $\unicode[STIX]{x1D707}$ is uniformly random over the interval $(-1,1)$ . Thus one only has to do integrations of the form $\int \unicode[STIX]{x1D707}^{n}\,\text{d}\unicode[STIX]{x1D707}$ for the $q_{i}$ averaging, which become trivial. Carrying out the above steps, a lengthy but straightforward calculation gives,

(3.27) $$\begin{eqnarray}\frac{\hat{M}_{4}(k,t)}{\unicode[STIX]{x1D70F}}=\frac{9}{350}\unicode[STIX]{x1D70F}(\unicode[STIX]{x1D705}_{t}q^{2})^{2}\left[k^{4}\frac{\unicode[STIX]{x2202}^{4}\hat{M}}{\unicode[STIX]{x2202}k^{4}}+12k^{3}\frac{\unicode[STIX]{x2202}^{3}\hat{M}}{\unicode[STIX]{x2202}k^{3}}+24k^{2}\frac{\unicode[STIX]{x2202}^{2}\hat{M}}{\unicode[STIX]{x2202}k^{2}}-24k\frac{\unicode[STIX]{x2202}\hat{M}}{\unicode[STIX]{x2202}k}\right].\end{eqnarray}$$

We define dimensionless time coordinates $\bar{t}=t\unicode[STIX]{x1D705}_{t}q^{2}$ and $\bar{\unicode[STIX]{x1D70F}}=\unicode[STIX]{x1D70F}\unicode[STIX]{x1D705}_{t}q^{2}$ , dividing by a turbulent diffusion time scale $(\unicode[STIX]{x1D705}_{t}q^{2})^{-1}$ . Again taking the limit $\unicode[STIX]{x1D70F}\rightarrow 0$ , keeping $\unicode[STIX]{x1D705}_{t}$ fixed and adding in the contribution due to $\hat{M}_{4}$ , we get for the generalized Kraichnan equation

(3.28) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\hat{M}}{\unicode[STIX]{x2202}\bar{t}} & = & \displaystyle \frac{1}{5}\left[k^{2}\frac{\unicode[STIX]{x2202}^{2}\hat{M}}{\unicode[STIX]{x2202}k^{2}}+4k\frac{\unicode[STIX]{x2202}\hat{M}}{\unicode[STIX]{x2202}k}\right]-2\frac{\unicode[STIX]{x1D705}}{\unicode[STIX]{x1D705}_{t}}\frac{k^{2}}{q^{2}}\hat{M}\nonumber\\ \displaystyle & & \displaystyle +\,\bar{\unicode[STIX]{x1D70F}}\frac{9}{350}\left[k^{4}\frac{\unicode[STIX]{x2202}^{4}\hat{M}}{\unicode[STIX]{x2202}k^{4}}+12k^{3}\frac{\unicode[STIX]{x2202}^{3}\hat{M}}{\unicode[STIX]{x2202}k^{3}}+24k^{2}\frac{\unicode[STIX]{x2202}^{2}\hat{M}}{\unicode[STIX]{x2202}k^{2}}-24k\frac{\unicode[STIX]{x2202}\hat{M}}{\unicode[STIX]{x2202}k}\right].\end{eqnarray}$$

Again the first line of (3.28) is the standard K68 equation in Fourier space while the second line gives the leading-order finite $\unicode[STIX]{x1D70F}$ correction to the Kraichnan equation. The generalized Kraichnan equation allows eigensolutions of the form $\hat{M}(k,t)=\tilde{M}(k)\text{e}^{-\unicode[STIX]{x1D6FE}\bar{t}}$ . We consider such solutions in the next section. Note that to get a steady state solution with $\unicode[STIX]{x1D6FE}=0$ , one would also need a source term at small $k$ as discussed earlier; if there is no source the scalar fluctuations would decay. We consider both the steady state and decaying cases below. The decaying case, which is more subtle, needs a different treatment, obtained by solving for the evolution of $\hat{M}$ as an initial value problem in terms of the Green function solution of (3.28). The eigensolutions are useful in providing basis eigenmodes and eigenvalues for deriving the Green’s function, and thus we keep $\unicode[STIX]{x1D6FE}$ to be non-zero below.

4.1 The scalar spectrum from a scaling solution

Consider first the scalar spectrum in the viscous–convective regime, where $q\ll k\ll k_{\unicode[STIX]{x1D705}}$ . Thus (3.28) or (4.4) is applicable and the diffusive term can also be neglected. In this regime we see that (4.4), and in fact the original (3.28), are scale free, as scaling $k\rightarrow ck$ leaves them invariant. Therefore they admit power-law solutions, of the form $\tilde{M}(k)=M_{0}k^{-\unicode[STIX]{x1D706}}$ , where from (4.4), $\unicode[STIX]{x1D706}$ is determined by

(4.5a-c ) $$\begin{eqnarray}\unicode[STIX]{x1D706}^{2}-3\unicode[STIX]{x1D706}+\frac{5[\unicode[STIX]{x1D6FE}+(9/7)\unicode[STIX]{x1D6FE}_{0}\bar{\unicode[STIX]{x1D70F}}]}{1-(9/14)\unicode[STIX]{x1D6FE}_{0}\bar{\unicode[STIX]{x1D70F}}}=0;\quad \unicode[STIX]{x1D706}=\frac{3}{2}\pm D;\quad D=\frac{3}{2}\left[1-\frac{(20/9)[\unicode[STIX]{x1D6FE}+(9/7)\unicode[STIX]{x1D6FE}_{0}\bar{\unicode[STIX]{x1D70F}}]}{(1-(9/14)\unicode[STIX]{x1D6FE}_{0}\bar{\unicode[STIX]{x1D70F}})}\right]^{1/2}.\end{eqnarray}$$

One can consider the following two cases.

4.1.2 The case of finite eigenvalue $\unicode[STIX]{x1D6FE}$

Now turn to the case when there is no source of scalar fluctuations, and so they necessarily decay with possibly a superposition of eigenfunctions with finite $\unicode[STIX]{x1D6FE}\neq 0$ . An approximate value for the leading (smallest) eigenvalue $\unicode[STIX]{x1D6FE}$ , for $R_{\unicode[STIX]{x1D705}}\gg 1$ , can be got following an argument from Gruzinov, Cowley & Sudan (Reference Gruzinov, Cowley and Sudan1996), who looked at a similar power-law solution of the small-scale dynamo problem. In order to satisfy the boundary condition that the $k^{-\unicode[STIX]{x1D706}}$ eigenmode is regular at both large and small $k$ , they argued that the limiting eigenvalue $\unicode[STIX]{x1D6FE}$ (or the growth rate in the dynamo problem) is given by $\unicode[STIX]{x1D6FE}(\unicode[STIX]{x1D706}_{m})$ , with $\unicode[STIX]{x1D706}_{m}$ the value of $\unicode[STIX]{x1D706}$ where $\text{d}\unicode[STIX]{x1D6FE}/\text{d}\unicode[STIX]{x1D706}=0$ . This leads to

(4.6a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D6FE}_{0}\approx {\textstyle \frac{9}{20}}=\bar{\unicode[STIX]{x1D6FE}}_{0},\quad \text{and}\quad \unicode[STIX]{x1D6FE}\approx {\textstyle \frac{9}{20}}\left(1-{\textstyle \frac{63}{40}}\bar{\unicode[STIX]{x1D70F}}\right)=\bar{\unicode[STIX]{x1D6FE}}.\end{eqnarray}$$

Note that (4.6) also implies that $D\approx 0$ . We will see below that, a more exact treatment picks out a small imaginary value of $D$ , so as to also satisfy the regularity condition on $\tilde{M}(k)$ at small $k$ . The fact that $D$ is imaginary and small implies that the two roots of $\unicode[STIX]{x1D706}$ in (4.5) are complex conjugate, with a small imaginary value of $\unicode[STIX]{x1D706}$ . More importantly, the real part is $\unicode[STIX]{x1D706}_{R}=3/2$ independent of $\unicode[STIX]{x1D70F}$ , and independent of the value of the eigenvalue $\unicode[STIX]{x1D6FE}$ . The eigenfunctions (from the power-law solution, $k^{-3/2\pm \text{i}|D|}=k^{-3/2}\exp (\pm \text{i}|D|\ln k)$ ) are then of the form

(4.7) $$\begin{eqnarray}\tilde{M}(k)=k^{-3/2}\sin (|D|\ln (k)+c_{1}),\end{eqnarray}$$

where the only $\unicode[STIX]{x1D6FE}$ and $\unicode[STIX]{x1D70F}$ dependencies are in $D$ and the constant $c_{1}$ which are to be fixed from the boundary conditions at large and small $k$ . However, due to the expected small value of $D$ and the fact that the phase of the sine function is only logarithmically dependent on $k$ , we expect the dominant behaviour of the eigenfunctions to be independent of $\unicode[STIX]{x1D70F}$ . This implies that the spectrum, which would be superposition of such eigenfunctions (with different $\unicode[STIX]{x1D6FE}$ or $D$ ), will also be independent of $\unicode[STIX]{x1D70F}$ . This is an important new result, that the passive scalar spectrum in the decaying case is also independent of $\unicode[STIX]{x1D70F}$ to leading order in $\unicode[STIX]{x1D70F}$ . Further, one may also expect that the spectrum itself dominantly falls off as $\tilde{M}(k)\propto k^{-3/2}$ or $E_{\unicode[STIX]{x1D703}}\propto k^{1/2}$ , independent of $\unicode[STIX]{x1D70F}$ . We will check this in more detail below, when solving the initial value problem for the decay. (We note in passing that for the corresponding two-dimensional problem for passive scalar decay in the $\unicode[STIX]{x1D70F}\rightarrow 0$ limit, Nazarenko & Laval (Reference Nazarenko and Laval2000) obtained a flat $k^{0}$ spectrum in the Batchelor regime. Both the two- and three-dimensional results are consistent with the general expectation that, for the Kraichnan model in the Batchelor regime, the scalar eigenfunction during decay goes as $E_{\unicode[STIX]{x1D703}}\propto k^{d/2-1}$ in $d$ -dimensions (Chertkov & Lebedev Reference Chertkov and Lebedev2003; Schekochihin, Haynes & Cowley Reference Schekochihin, Haynes and Cowley2004).) Further, from (4.6), we see that the limiting decay rate $\unicode[STIX]{x1D6FE}$ is of order of the turbulent diffusion rate at the forcing scale, since we normalized $t$ by $(\unicode[STIX]{x1D705}_{t}q^{2})^{-1}$ . We also see that a finite $\unicode[STIX]{x1D70F}$ decreases this rate at which the passive scalar decays, and thus the Kraichnan model overestimates the decay rate. In the scaling solution, we neglected the diffusive term. We now consider the effect of including this term.

4.2 The scalar spectrum including diffusive effects

In fact one can get an exact solution of (4.4) which includes the diffusive term. For this let us substitute $\tilde{M}(k)=k^{-3/2}K(k)$ into (4.4). This leads to a modified Bessel equation for $K(k)$ given by

(4.8) $$\begin{eqnarray}x^{2}\frac{\text{d}^{2}K}{\text{d}x^{2}}+x\frac{\text{d}K}{\text{d}x}-[x^{2}+\unicode[STIX]{x1D708}^{2}]K=0,\end{eqnarray}$$

where we have defined $\unicode[STIX]{x1D708}^{2}\equiv D^{2}$ , and

(4.9a-d ) $$\begin{eqnarray}x=\frac{k}{k_{d}},\quad k_{d}=q\left(\frac{\unicode[STIX]{x1D705}_{t}}{10\unicode[STIX]{x1D705}}\right)^{1/2},\quad \unicode[STIX]{x1D6FC}=1-\frac{9}{14}\unicode[STIX]{x1D6FE}_{0}\bar{\unicode[STIX]{x1D70F}},\quad \unicode[STIX]{x1D708}^{2}=\frac{9}{4\unicode[STIX]{x1D6FC}}\left[1-\frac{20}{9}\unicode[STIX]{x1D6FE}-\frac{7}{2}\unicode[STIX]{x1D6FE}_{0}\bar{\unicode[STIX]{x1D70F}}\right].\end{eqnarray}$$

We note that the Péclet number $R_{\unicode[STIX]{x1D705}}\sim \unicode[STIX]{x1D705}_{t}/\unicode[STIX]{x1D705}$ , and so $k_{d}\sim qR_{\unicode[STIX]{x1D705}}^{1/2}=k_{\unicode[STIX]{x1D705}}$ , the diffusive wavenumber and $x\sim k/k_{d}$ is basically $k$ normalized by the diffusive wavenumber. Thus in the viscous–convective regime $x\ll 1$ , while in the diffusive regime, $x>1$ . Further, consistent with neglecting the effect of $\bar{\unicode[STIX]{x1D70F}}$ on the diffusive term in the extended Kraichnan equation, for $R_{\unicode[STIX]{x1D705}}\gg 1$ , we have also neglected its effect on $k_{d}$ above. (We note in passing that even at $k=q$ , $x\sim 1/R_{\unicode[STIX]{x1D705}}^{1/2}\ll 1$ , for large Péclet number $R_{\unicode[STIX]{x1D705}}\gg 1$ , although the small $qr$ approximation will become inaccurate).

4.2.1 Steady state solution

First consider the steady state solution with $\unicode[STIX]{x1D6FE}=0=\unicode[STIX]{x1D6FE}_{0}$ . As we remarked earlier, this assumes that there is a source of scalar fluctuations which contributes at low wavenumbers, outside the viscous convective range, so that (4.4) or (4.8) is still valid. For the steady state case, we have $\unicode[STIX]{x1D708}^{2}=9/4$ . Then the two independent solutions of (4.8) are the modified Bessel functions of the first kind, $I_{3/2}(x)$ and of the second kind $K_{3/2}(x)$ . We also have to satisfy the boundary condition that as $x\rightarrow \infty$ , $K(x)\rightarrow 0$ . This rules out the solution proportional to $I_{3/2}$ , and so we have

(4.10) $$\begin{eqnarray}\tilde{M}(k)=k^{-3/2}K_{3/2}(k/k_{d}).\end{eqnarray}$$

Note that as $x\rightarrow \infty$ , $K_{3/2}(x)\rightarrow x^{-1/2}\text{e}^{-x}$ and so $\tilde{M}(k)\rightarrow k^{-2}\exp (-k/k_{d})$ . This then implies that $E_{\unicode[STIX]{x1D703}}(k)=4\unicode[STIX]{x03C0}k^{2}\tilde{M}(k)\rightarrow \exp (-k/k_{d})$ , and so dies exponentially beyond the diffusive wavenumber. Such an exponential fall off for the steady state case has been shown earlier by K68 and Kraichnan (Reference Kraichnan1974). And in the viscous–convective regime when $x\ll 1$ , we have $K_{3/2}(x)\rightarrow \unicode[STIX]{x1D6E4}(3/2)2^{1/2}x^{-3/2}$ , so $\tilde{M}(k)\propto k^{-3}$ which implies $E_{\unicode[STIX]{x1D703}}(k)\propto 1/k$ , or the Batchelor spectrum. Thus as in the scaling solution above, the more complete solution in (4.8) also shows that the Batchelor spectrum occurs in the viscous–convective range $q\ll k\ll k_{\unicode[STIX]{x1D705}}$ , even for a finite $\unicode[STIX]{x1D70F}$ , to the leading order in $\unicode[STIX]{x1D70F}$ . In addition (4.10) extrapolates the Batchelor spectrum smoothly to an exponential fall off in the diffusive regime, in agreement with K68 and Kraichnan (Reference Kraichnan1974).

4.3 The scalar spectrum during decay

Now consider the case when there is no source for scalar fluctuations. We then expect that the passive scalar fluctuations will decay. We are interested in how the decay rate and the scalar power spectrum in this case are altered due to the effect of a finite $\unicode[STIX]{x1D70F}$ . For the Kraichnan problem itself, where one considers passive scalar decay with $\bar{\unicode[STIX]{x1D70F}}\rightarrow 0$ , a number of different approaches have been adopted, either by looking at the decaying eigenfunction (Kulsrud & Anderson Reference Kulsrud and Anderson1992; Schekochihin et al. Reference Schekochihin, Haynes and Cowley2004) or in terms of solving (3.28) as an initial value problem (Balkovsky & Fouxon Reference Balkovsky and Fouxon1999; Chertkov & Lebedev Reference Chertkov and Lebedev2003). Here we solve for the evolution of $\hat{M}$ as an initial value problem in terms of the Green’s function solution of (3.28) in the Landau–Lifshitz approximation.

For finding this Green’s function, we seek eigenfunction solutions to (4.8) which vanish at $x\rightarrow \infty$ and are regular as $x\rightarrow 0$ . Such eigensolutions need to have $\unicode[STIX]{x1D708}^{2}=-(\bar{\unicode[STIX]{x1D708}})^{2}<0$ , and are given by the modified Bessel function of imaginary order $\bar{\unicode[STIX]{x1D708}}>0$ . That is

(4.11) $$\begin{eqnarray}\tilde{M}(k)=\tilde{M}_{\bar{\unicode[STIX]{x1D708}}}(k)=k^{-3/2}K_{i\bar{\unicode[STIX]{x1D708}}}(k/k_{d}),\end{eqnarray}$$

where from (4.9), $\bar{\unicode[STIX]{x1D708}}$ is related to the eigenvalue $\unicode[STIX]{x1D6FE}$ ,

(4.12) $$\begin{eqnarray}\bar{\unicode[STIX]{x1D708}}^{2}=\frac{9}{4\unicode[STIX]{x1D6FC}}\left[\frac{20}{9}\unicode[STIX]{x1D6FE}+\frac{7}{2}\unicode[STIX]{x1D6FE}_{0}\bar{\unicode[STIX]{x1D70F}}-1\right]=\frac{9}{4}\frac{[20\unicode[STIX]{x1D6FE}/9+7\unicode[STIX]{x1D6FE}_{0}\bar{\unicode[STIX]{x1D70F}}/2-1]}{(1-9\unicode[STIX]{x1D6FE}_{0}\bar{\unicode[STIX]{x1D70F}}/14)}.\end{eqnarray}$$

In the limit $x\rightarrow \infty$ , $\tilde{K}_{i\bar{\unicode[STIX]{x1D708}}}(x)\rightarrow x^{-1/2}\text{e}^{-x}$ still and the eigenfunction and hence the spectrum vanishes for $k\gg k_{d}$ , again falling off exponentially with $k$ . In the opposite limit of $x\rightarrow 0$ , we have (Abramowitz & Stegun Reference Abramowitz and Stegun1972) (see also chap. 10 by Olver and Maximon in (DLMF Reference Olver, Olde Daalhuis, Lozier, Schneider, Boisvert, Clark, Miller and Saunders2016))

(4.13) $$\begin{eqnarray}\tilde{M}_{\bar{\unicode[STIX]{x1D708}}}(k)\rightarrow -\left(\frac{\unicode[STIX]{x03C0}}{\bar{\unicode[STIX]{x1D708}}\sinh (\unicode[STIX]{x03C0}\bar{\unicode[STIX]{x1D708}})}\right)^{1/2}k^{-3/2}\sin (\bar{\unicode[STIX]{x1D708}}\ln (k/2k_{d})+c_{\bar{\unicode[STIX]{x1D708}}}).\end{eqnarray}$$

Here $c_{\bar{\unicode[STIX]{x1D708}}}$ is a $\bar{\unicode[STIX]{x1D708}}$ -dependent real constant defined by the relation,

(4.14) $$\begin{eqnarray}\unicode[STIX]{x1D6E4}(1+\text{i}\bar{\unicode[STIX]{x1D708}})=\left(\frac{\unicode[STIX]{x03C0}\bar{\unicode[STIX]{x1D708}}}{\sinh (\unicode[STIX]{x03C0}\bar{\unicode[STIX]{x1D708}})}\right)^{1/2}\exp (-\text{i}c_{\bar{\unicode[STIX]{x1D708}}}),\end{eqnarray}$$

and so $c_{\bar{\unicode[STIX]{x1D708}}}\rightarrow \bar{c}\bar{\unicode[STIX]{x1D708}}$ for $\bar{\unicode[STIX]{x1D708}}\ll 1$ , where $\bar{c}$ is the Euler constant. Expectedly the solution $\tilde{M}$ in small $k/k_{d}$ limit, given in (4.13) is the same as the scaling solution given by (4.7), with now the constants in that solution fixed by matching also to the diffusive range.

It is interesting to note that the differential operator, say ${\mathcal{L}}$ , defined on the right-hand side of the generalized Kraichnan equation (3.28) is self-adjoint with respect to the inner product

(4.15) $$\begin{eqnarray}\langle \tilde{M}_{1}(k),\tilde{M}_{2}(k)\rangle =\int _{0}^{\infty }k^{2}\,\text{d}k\tilde{M}_{1}^{\ast }(k)\tilde{M}_{2}(k).\end{eqnarray}$$

(We thank a referee for bringing this to our attention.) The eigenfunctions of ${\mathcal{L}}$ are given by (4.11) (in the Landau–Lifshitz approximation), with a real eigenvalue $-\unicode[STIX]{x1D6FE}$ , where $\unicode[STIX]{x1D6FE}$ is related to $\bar{\unicode[STIX]{x1D708}}^{2}=-\unicode[STIX]{x1D708}^{2}$ as in (4.12). They can then be labelled by the value of $\bar{\unicode[STIX]{x1D708}}$ , and under the inner product (4.15), also satisfy the orthogonality condition,

(4.16) $$\begin{eqnarray}\langle \tilde{M}_{\bar{\unicode[STIX]{x1D708}}},\tilde{M}_{\bar{\unicode[STIX]{x1D708}}^{\prime }}\rangle =\int _{0}^{\infty }\frac{\text{d}k}{k}K_{i\bar{\unicode[STIX]{x1D708}}}(k/k_{d})K_{i\bar{\unicode[STIX]{x1D708}}^{\prime }}(k/k_{d})=\unicode[STIX]{x1D6FF}(\bar{\unicode[STIX]{x1D708}}-\bar{\unicode[STIX]{x1D708}}^{\prime })\left[\frac{\unicode[STIX]{x03C0}^{2}}{2\bar{\unicode[STIX]{x1D708}}\sinh (\unicode[STIX]{x03C0}\bar{\unicode[STIX]{x1D708}})}\right].\end{eqnarray}$$

The latter orthogonality relation of the modified Bessel functions of imaginary order is discussed for example in Szmytkowski & Bielski (Reference Szmytkowski and Bielski2010). Thus the eigenfunctions $\tilde{M}_{\bar{\unicode[STIX]{x1D708}}}(k)$ given in (4.11) are orthogonal under (4.15) and also delta function normalizable, like the free particle eigenfunctions in quantum mechanics. This allows us to expand the scalar spectrum $\hat{M}(k,t)$ in terms of the orthogonal eigenfunctions of ${\mathcal{L}}$ as,

(4.17) $$\begin{eqnarray}\hat{M}(k,\bar{t})=\int _{0}^{\infty }\,\text{d}\unicode[STIX]{x1D707}g(\unicode[STIX]{x1D707},\bar{t})\tilde{M}_{\unicode[STIX]{x1D707}}=\int _{0}^{\infty }\,\text{d}\unicode[STIX]{x1D707}g(\unicode[STIX]{x1D707},\bar{t})k^{-3/2}K_{i\unicode[STIX]{x1D707}}(k/k_{d}).\end{eqnarray}$$

(The above transform without the $k^{-3/2}$ factor is called the Kontorovich–Lebedev transform (DLMF Reference Olver, Olde Daalhuis, Lozier, Schneider, Boisvert, Clark, Miller and Saunders2016).) This transform can be inverted by multiplying (4.17) by $\tilde{M}_{\unicode[STIX]{x1D707}^{\prime }}(k)$ , integrating over $k^{2}\,\text{d}k$ and using the orthogonality condition (4.16), to obtain

(4.18) $$\begin{eqnarray}g(\unicode[STIX]{x1D707},\bar{t})=\frac{2\unicode[STIX]{x1D707}\sinh (\unicode[STIX]{x03C0}\unicode[STIX]{x1D707})}{\unicode[STIX]{x03C0}^{2}}\int _{0}^{\infty }\,\text{d}kk^{2}\tilde{M}_{\unicode[STIX]{x1D707}}(k)\hat{M}(k,\bar{t}).\end{eqnarray}$$

We use the evolution equation (3.28) in (4.17), and note that the eigenfunctions of ${\mathcal{L}}$ are given by $\tilde{M}_{\unicode[STIX]{x1D707}}$ (in the Landau–Lifshitz approximation), with a real eigenvalue $-\unicode[STIX]{x1D6FE}$ . Then the function $g$ is given by

(4.19) $$\begin{eqnarray}g(\unicode[STIX]{x1D707},\bar{t})=g(\unicode[STIX]{x1D707},0)\text{e}^{-\unicode[STIX]{x1D6FE}(\unicode[STIX]{x1D707})\bar{t}}=\text{e}^{-\unicode[STIX]{x1D6FE}(\unicode[STIX]{x1D707})\bar{t}}\frac{2\unicode[STIX]{x1D707}\sinh (\unicode[STIX]{x03C0}\unicode[STIX]{x1D707})}{\unicode[STIX]{x03C0}^{2}}\int _{0}^{\infty }\frac{\text{d}k^{\prime }}{k^{\prime }}K_{i\unicode[STIX]{x1D707}}(k^{\prime }/k_{d}){k^{\prime }}^{3/2}\hat{M}(k^{\prime },0).\end{eqnarray}$$

Here $g(\unicode[STIX]{x1D707},0)$ is the initial value of $g$ , and has itself been obtained by evaluating (4.18) at $\bar{t}=0$ , and $\hat{M}(k,0)$ is the initial scalar spectrum. Substituting $g(\unicode[STIX]{x1D707},\bar{t})$ from (4.19) into (4.17), allows one to write down the Green’s function solution for the evolution of $\hat{M}(k,t)$ as

(4.20) $$\begin{eqnarray}\hat{M}(k,t)=k^{-3/2}\int _{0}^{\infty }\,\text{d}\unicode[STIX]{x1D707}\text{e}^{-\unicode[STIX]{x1D6FE}(\unicode[STIX]{x1D707})\bar{t}}\frac{2\unicode[STIX]{x1D707}\sinh \unicode[STIX]{x03C0}\unicode[STIX]{x1D707}}{\unicode[STIX]{x03C0}^{2}}\int _{0}^{\infty }\frac{\text{d}k^{\prime }}{k^{\prime }}K_{i\unicode[STIX]{x1D707}}(k/k_{d})K_{i\unicode[STIX]{x1D707}}(k^{\prime }/k_{d}){k^{\prime }}^{3/2}\hat{M}(k^{\prime },0).\end{eqnarray}$$

We note that this solution is identical to the solution got by Balkovsky & Fouxon (Reference Balkovsky and Fouxon1999) in the $\unicode[STIX]{x1D70F}\rightarrow 0$ limit, but now generalized to case of finite $\unicode[STIX]{x1D70F}$ .

To study this solution in greater detail, we need the explicit form of $\unicode[STIX]{x1D6FE}(\unicode[STIX]{x1D707})$ . To obtain this relation, it is useful to replace $\unicode[STIX]{x1D6FE}_{0}\unicode[STIX]{x1D70F}$ by $\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D70F}$ to leading order in $\unicode[STIX]{x1D70F}$ before inverting (4.12). We find

(4.21) $$\begin{eqnarray}\unicode[STIX]{x1D6FE}(\unicode[STIX]{x1D707})=\bar{\unicode[STIX]{x1D6FE}}+\frac{4}{9}\unicode[STIX]{x1D707}^{2}\bar{\unicode[STIX]{x1D6FE}}~\frac{(1-9\bar{\unicode[STIX]{x1D6FE}}\bar{\unicode[STIX]{x1D70F}}/14)}{(1+2\unicode[STIX]{x1D707}^{2}\bar{\unicode[STIX]{x1D6FE}}\bar{\unicode[STIX]{x1D70F}}/7)}.\end{eqnarray}$$

Note that $\unicode[STIX]{x1D6FE}(\unicode[STIX]{x1D707})$ increases monotonically with $\unicode[STIX]{x1D707}$ from a value $\unicode[STIX]{x1D6FE}=\bar{\unicode[STIX]{x1D6FE}}$ , when $\unicode[STIX]{x1D707}=0$ to $\unicode[STIX]{x1D6FE}\rightarrow \unicode[STIX]{x1D6FE}_{u}=14/(9\bar{\unicode[STIX]{x1D70F}})$ as $\unicode[STIX]{x1D707}\rightarrow \infty$ . Thus $\unicode[STIX]{x1D6FC}$ never becomes negative as required for the consistency of the solution. We also have $\unicode[STIX]{x1D6FE}_{u}/\bar{\unicode[STIX]{x1D6FE}}=280/81\bar{\unicode[STIX]{x1D70F}}\gg 1$ for $\bar{\unicode[STIX]{x1D70F}}\ll 1$ . Substituting this $\unicode[STIX]{x1D6FE}(\unicode[STIX]{x1D707})$ into (4.20) we then have

(4.22) $$\begin{eqnarray}\displaystyle \hat{M}(k,t) & = & \displaystyle k^{-3/2}\text{e}^{-\bar{\unicode[STIX]{x1D6FE}}\bar{t}}\int _{0}^{\infty }\,\text{d}\unicode[STIX]{x1D707}\exp \left(-\frac{4}{9}\unicode[STIX]{x1D707}^{2}\bar{\unicode[STIX]{x1D6FE}}\bar{t}\frac{(1-9\bar{\unicode[STIX]{x1D6FE}}\bar{\unicode[STIX]{x1D70F}}/14)}{(1+2\unicode[STIX]{x1D707}^{2}\bar{\unicode[STIX]{x1D6FE}}\bar{\unicode[STIX]{x1D70F}}/7)}\right)\frac{2\unicode[STIX]{x1D707}\sinh \unicode[STIX]{x03C0}\unicode[STIX]{x1D707}}{\unicode[STIX]{x03C0}^{2}}\nonumber\\ \displaystyle & & \displaystyle \times \,\int _{0}^{\infty }\frac{\text{d}k^{\prime }}{k^{\prime }}K_{i\unicode[STIX]{x1D707}}(k/k_{d})K_{i\unicode[STIX]{x1D707}}(k^{\prime }/k_{d}){k^{\prime }}^{3/2}\hat{M}(k^{\prime },0).\end{eqnarray}$$

We see that there is an overall $k^{-3/2}$ dependence of the spectrum independent of $\unicode[STIX]{x1D70F}$ , and also an overall exponential decay of the spectrum with a decay rate $\bar{\unicode[STIX]{x1D6FE}}$ (as in the scaling solution). The remaining $k$ and $\bar{t}$ dependence will come from evaluating the integral in (4.22) over $\unicode[STIX]{x1D707}$ .

Note that the $\exp (-\unicode[STIX]{x1D6FE}(\unicode[STIX]{x1D707})\bar{t})$ factor in this integral leads to a damping of the integral, which increases with $\unicode[STIX]{x1D707}$ . For example, with $\bar{\unicode[STIX]{x1D70F}}=0.1$ , the extra damping factor over and above $\text{e}^{-\bar{\unicode[STIX]{x1D6FE}}\bar{t}}$ , as $\unicode[STIX]{x1D707}\rightarrow \infty$ goes as $\text{e}^{-15.6\bar{t}}$ , while there is no extra damping at $\unicode[STIX]{x1D707}=0$ . Meanwhile, from a power series expansion of $K_{i\unicode[STIX]{x1D707}}(k/k_{d})$ about $k=0$ (Dunster Reference Dunster1990), the combination $\sqrt{\unicode[STIX]{x1D707}\sinh (\unicode[STIX]{x03C0}\unicode[STIX]{x1D707})}K_{i\unicode[STIX]{x1D707}}$ is bounded as $\unicode[STIX]{x1D707}$ increases. Then the effect of $\unicode[STIX]{x1D6FE}(\unicode[STIX]{x1D707})$ damping results in the integral over $\unicode[STIX]{x1D707}$ , being very rapidly dominated by contributions from the lower range of $\unicode[STIX]{x1D707}$ . In fact as time increases, the range of $\unicode[STIX]{x1D707}$ which contributes significantly to the integral decreasing with time, and is limited to $\unicode[STIX]{x1D707}<\unicode[STIX]{x1D707}_{0}\approx 1/\sqrt{\bar{\unicode[STIX]{x1D6FE}}\bar{t}}$ , which at late times can become much smaller than unity. Further we have $K_{i\unicode[STIX]{x1D707}}(k/k_{d})\propto \sin (\unicode[STIX]{x1D707}\ln (k/2k_{d})+c_{\unicode[STIX]{x1D707}})$ , for $k/k_{d}\ll 1$ . Thus, due to the weak logarithmic dependence on $k$ compounded by the fact that $\unicode[STIX]{x1D707}<\unicode[STIX]{x1D707}_{0}\ll 1$ at late times, the integral has a weak $k$ dependence below the dissipative scales at late times. As a result, the $k$ dependence of the scalar spectrum coming from the integral can become subdominant to the overall $k^{-3/2}$ dependence of the spectrum coming from outside the integral. This occurs for a generic initial spectrum and is also independent of $\bar{\unicode[STIX]{x1D70F}}$ . On general grounds, we see from (4.22) the scalar spectrum during decay in the Batchelor regime has a dominant form $\tilde{M}\propto k^{-3/2}$ or $E_{\unicode[STIX]{x1D703}}(k)\propto k^{1/2}$ at late times, independent of $\bar{\unicode[STIX]{x1D70F}}$ . We will see these features more explicitly in the worked example below.

4.3.1 Numerical integration for scalar spectral evolution

One can integrate (4.22) numerically, to determine the scalar spectral evolution for a fixed $\bar{\unicode[STIX]{x1D70F}}$ and any initial spectrum. In order to obtain some physical insight, we focus on a simple example of the initial spectrum, which nevertheless illustrates all the features of interest. We choose a case where we assume that the initial spectrum is very strongly peaked at some small $k_{0}/k_{d}\ll 1$ ; specifically that $\hat{M}(k,0)=\unicode[STIX]{x1D6FF}(k-k_{0})/(4\unicode[STIX]{x03C0}k^{2})$ such that the total initial scalar power $\int 4\unicode[STIX]{x03C0}k^{2}\hat{M}(k,0)=1$ . We ask how this evolves in time. Substituting the above form of the initial spectrum allows one to do the $k^{\prime }$ integral in (4.22) trivially. In the small $k/k_{d}\ll 1$ limit, one can also do the $\unicode[STIX]{x1D707}$ integral analytically in an approximate manner, and we discuss that in appendix B. Here we do not make any approximations and simply do the $\unicode[STIX]{x1D707}$ integration in (4.22) numerically using Mathematica.

In figure 1 we show the evolution of scalar spectrum multiplied by a $k^{3/2}$ factor, to compensate the overall $k^{-3/2}$ factor in (4.22). Figures 1(a) and 1(b) correspond to the cases of $\bar{\unicode[STIX]{x1D70F}}=0.1$ and $\bar{\unicode[STIX]{x1D70F}}=0$ respectively. The decay starts from the initial power peaked at $k_{0}/k_{d}=0.1$ . A flat region of the curves in this figure would correspond to where the 3-D spectrum behaves as $\tilde{M}(k)\propto k^{-3/2}$ or a 1-D spectrum $E_{\unicode[STIX]{x1D703}}(k)\propto k^{1/2}$ . The topmost curve shows the compensated spectrum at $\bar{t}=1$ or one turbulent diffusion time after the decay began. Subsequent curves, from top to bottom, are shown for times $\bar{t}=2,4,6,8,\ldots ,20$ . We see from figure 1, that the compensated scalar power spreads with time as it decays, to wavenumbers both smaller and larger than $k_{0}$ . The initial delta function is already broadened significantly due to turbulent diffusion by $\bar{t}=1$ . The spectrum cannot spread to large $k>k_{d}$ due to diffusion damping. However, it can spread to arbitrarily small $k/k_{d}\ll 1$ in the Batchelor regime, as time increases. We see that as this happens, the compensated power spectrum becomes flat for a range of wavenumbers at the low $k$ end, and this range keeps increasing secularly with time. At much smaller $k$ (beyond what has been plotted) the spectrum is expected to be cutoff as the scalar power has not yet spread to those wavenumbers. And at large $k$ there is a faster cutoff due to diffusion. It is for a set of intermediate range of wavenumbers, in the Batchelor regime, that one would see $E_{\unicode[STIX]{x1D703}}(k)\propto k^{1/2}$ behaviour of the scalar spectrum during decay. More generally, we note by comparing the left and right panels of figure 1 that the qualitative form of the scalar spectrum is independent of $\bar{\unicode[STIX]{x1D70F}}$ . However, figure 1 also shows that scalar power spectrum for the $\bar{\unicode[STIX]{x1D70F}}=0.1$ case decays slower than $\bar{\unicode[STIX]{x1D70F}}=0$ , as expected.

Note that although we started with a specific initial spectrum, the qualitative features of the scalar spectrum during decay for $k/k_{d}\ll 1$ are expected to be the same, for any peaked initial spectrum.

Figure 1. (a,b) Show the decaying scalar spectrum, multiplied by $k^{3/2}$ factor, for the cases of $\bar{\unicode[STIX]{x1D70F}}=0.1$ and $\bar{\unicode[STIX]{x1D70F}}=0$ respectively. The wavenumber $k$ here is measured in units of $k_{d}$ . Initially power is peaked at $k_{0}/k_{d}=0.1$ . The topmost curve in each case shows the spectrum at $\bar{t}=1$ or one turbulent diffusion time after the decay began. Subsequent curves, from top to bottom, are shown for times $\bar{t}=2,4,6,8,\ldots ,20$ . A flat curve corresponds to the 3-D spectrum behaving as $\tilde{M}(k)\propto k^{-3/2}$ or a 1-D spectrum $E_{\unicode[STIX]{x1D703}}(k)\propto k^{1/2}$ .

Once the scalar power spectrum has spread to scales larger than the forcing scale (to $k<q$ ), the slower decay of the larger-scale modes by turbulent mixing could influence the decay of scalar power even for $k>q$ (see § 5.2). It could lead a more complicated behaviour of the decay phase (see Schekochihin et al. (Reference Schekochihin, Haynes and Cowley2004) for a model problem) than what occurs in the purely Batchelor regime which is the focus of the present work. Further theoretical study, also taking into account the large-scale (small $k$ ) cutoff of the velocity correlations, is important and will be useful to improve the understanding of the scalar decay, in both the Kraichnan limit and its finite $\unicode[STIX]{x1D70F}$ extensions. We now consider to what extent our theoretical predictions are borne out in direct numerical simulations of passive scalar mixing and decay.

5.1 Steady state case

We have used a resolution of $N=1024$ for two simulations choosing $k_{f}=1.5$ (run A) and $k_{f}=4$ (run B) respectively. The basic parameters of the simulations are summarized in table 1. The initial velocity is zero and the form of initial passive scalar is Gaussian random noise, with an arbitrary amplitude of 0.5. In order to obtain a steady state for the passive scalar evolution, we impose a constant gradient of the passive scalar in an arbitrarily chosen direction which is along the diagonal of the box. This corresponds to a force $f_{\unicode[STIX]{x1D703}}=\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D71E}$ with a constant $\unicode[STIX]{x1D71E}$ , and is a standard technique for achieving a steady state; it retains homogeneity but breaks isotropy at the largest scale.

Figure 2. (a) Shows the growing scalar spectrum as solid lines for run A from times $t=5{-}80$ every 5 time units. Time increases from the bottom to the top curve and the topmost curve shows the steady state scalar spectrum. The steady velocity spectrum is shown as a dashed line. Within $t=20$ the root mean square (r.m.s.) velocity rises to $u_{rms}\sim 0.15$ and varies between 0.13–0.17 of the sound speed. (b) Shows the corresponding spectra for run B. The scalar spectra are shown from $t=5{-}45$ every 5 time units. Within $t=10$ the r.m.s. velocity settles to $u_{rms}\sim 0.11$ in units of the sound speed.

Table 1. Summary of runs discussed in this paper.

In figure 2 we show as solid lines the evolution of the passive scalar power spectrum, $E_{\unicode[STIX]{x1D703}}$ , in equal time intervals. The steady state velocity spectrum is also shown, by a dashed line. Figure 2(a) is for run A with $k_{f}=1.5$ and figure 2(b) is for run B with $k_{f}=4$ . The root mean square (r.m.s.) velocities in units of the sound speed are $u_{rms}\sim 0.15$ and $u_{rms}\sim 0.11$ in runs A and B respectively, such that the incompressibility condition is well satisfied, and $\unicode[STIX]{x1D70C}$ is nearly constant. The viscosity has been set high enough that the flow is almost single scale, and the diffusivity low enough that the Péclet and Schmidt numbers are high. We have $Sc=400$ , $R_{\unicode[STIX]{x1D705}}\sim 4000$ and $Sc=100$ , $R_{\unicode[STIX]{x1D705}}=250$ for the runs A and B respectively. Thus one would expect indeed to obtain a significant viscous–convective range of wavenumbers. The simulations are then suited to test if one obtains the Batchelor $E_{\unicode[STIX]{x1D703}}\propto 1/k$ spectrum, in steady state. The simulation were run for a sufficiently long time to obtain steady state in evolution of the passive scalar, as can be seen from the overlap of the spectrum at final times. As can be seen in figure 2(a,b), $E_{\unicode[STIX]{x1D703}}$ reaches a steady state and exhibits a scale dependence of $k^{-1}$ to a reasonable degree in the viscous–convective range larger than $k_{f}$ . While Kraichnan’s idealistic model of delta correlation predicted $k^{-1}$ spectrum, the solution to our extended Kraichnan equation showed that to the leading order, the finite time correlation does not change this scale dependence. This is also predicted by the Lagrangian analysis (Falkovich et al. Reference Falkovich, Gawȩdzki and Vergassola2001). We find that our DNS which explores Schmidt numbers up to 400 bears this out. We note that earlier DNS, with higher $Re$ than obtained here, but with $Sc$ up to 64, have also found evidence for the Batchelor spectrum in the viscous–convective range (see Donzis et al. (Reference Donzis, Sreenivasan and Yeung2010), Gotoh & Yeung (Reference Gotoh, Yeung, Davidson, Kaneda and Sreenivasan2013) and references therein). During the course of this work, we also came across some recent very high resolution hybrid simulations of passive scalar mixing with $Sc\sim 200{-}1000$ comparable to ours, which also report reasonable agreement with the steady state Batchelor spectrum (Gotoh, Watanabe & Miura Reference Gotoh, Watanabe and Miura2014). There have been interesting predictions for the behaviour of higher-order correlators in the steady state and the influence of the diffusive scale even on larger scales (Balkovsky et al. Reference Balkovsky, Chertkov, Kolokolov and Lebedev1995; Chertkov et al. Reference Chertkov, Kolokolov and Lebedev2007). It would be of interest to check for such effects in future work.

Figure 3. (a) Shows the decaying scalar spectrum (solid lines) from $t=80{-}140$ in steps of 10 time units for run A. The corresponding decay phase for run B is shown in (b), from $t=40{-}105$ in steps of 5 time units. Time now increases from the top to the bottom curves. The topmost curve again shows the steady state spectrum from which the decay began, once the scalar forcing is turned off. The dashed lines again show the steady state velocity spectrum.

Figure 4. The figure shows the decaying scalar spectrum (solid lines) when one starts from initial power peaked at $k_{0}=30$ . The topmost curve shows the spectrum at $t=5$ or after approximately 1 eddy turnover time after the decay began, with the scalar forcing is turned off. Subsequent curves, from top to bottom, are shown for times $t=5,10,20,\ldots ,150$ . The blue solid straight line shown at the bottom is $E(k)\propto k^{1/2}$ . The dashed line again shows the (almost) steady state velocity spectrum. Within $t=5$ the r.m.s. velocity rises to $u_{rms}\sim 0.1$ and varies between 0.11–0.15, of the sound speed.

5.2 DNS of passive scalar decay

We next examine the decay of the passive scalar through our DNS. We carry out two types of DNS to explore this phase. In the first case, to obtain such a decay, we turn off the imposed constant gradient in the passive scalar equation. The passive scalar spectrum then decays from its initial steady state, as is shown in figure 3. We find that the spectrum flattens as it decays, but never becomes close to the $k^{1/2}$ form predicted by the analytic argument. This can be seen to occur in both run A (figure 3 a) and in run B (b). It could be that the box-scale mode decays slower due to turbulent diffusion, at a rate ${\sim}k_{box}^{2}\unicode[STIX]{x1D705}_{t}$ , with $k_{box}=1$ in our DNS, compared to the small-scale modes. Our analytical considerations predict that they decay at the rate ${\sim}q^{2}\unicode[STIX]{x1D705}_{t}$ , with $q\sim k_{f}$ . Some evidence for such an idea is seen in run B, where there is a clear separation of scales. Then the box-scale modes could act as a source for the small-scale fluctuations, leading to a spectrum flatter than $k^{1/2}$ but not as steep as $k^{-1}$ . A similar idea has been explored in Schekochihin et al. (Reference Schekochihin, Haynes and Cowley2004).

In view of the above, we carry out a second type of DNS where we make sure that there is no large-scale power in the scalar fluctuations at the initial time. In particular, we start the decay with scalar power peaked initially at some wavenumber $k_{0}\gg k_{f}$ and let it decay due to the random motions driven at the forcing wavenumber $k_{f}$ . This is analogous to the case of the solved example, presented in § 4.3.1 and appendix B. In figure 4, we show the results from such a high resolution run ( $N=1024$ , run C), where the initial scalar spectrum is a delta function at $k_{0}=30$ and $k_{f}=1.5$ . The subsequent decay of the scalar spectrum is shown as a sequence of solid curves in figure 4. The curves from top to bottom correspond to times $t=5,10,20,30,\ldots ,150$ . The blue straight line is the 1-D scalar spectrum which is predicted at the low $k$ end in the Batchelor regime, $E_{\unicode[STIX]{x1D703}}(k)\propto k^{1/2}$ .

The decay now is more in line with expectations. By $t=5$ , which is roughly one eddy turnover time, the spectrum has spread to have almost a Gaussian shape, although still peaked at around $k\sim k_{0}=30$ . As the decay proceeds, the scalar power spreads, with its half-width increasing in time. The slope at low $k$ is always positive now, but steeper than the $k^{1/2}$ form because of the influence of the expected rise of power towards $k_{0}$ , the Gaussian part in (B 3). Power reaches the diffusive cutoff scale $k_{d}\sim 90$ , by $t=20$ , after which it cannot spread further to larger $k$ , but continues to spread to smaller and smaller $k<k_{0}$ . The peak of the spectrum thus shifts to smaller and smaller $k<k_{0}$ as the scalar power decays. The power at $k=1$ increases until approximately $t=60$ , stays roughly constant and later after $t=90$ starts decaying. In fact after approximately $t=90$ until $t=140$ , the spectra shift down almost like an eigenfunction, preserving their shape. Towards the end of the simulation, the low- $k$ slope of the scalar spectrum approaches the $k^{1/2}$ form (see figure 4), with the peak of the spectrum at $k\sim 10$ . Eventually, the box-scale modes are expected to decay slower than the modes with $k>k_{f}$ and one could see a further flattening of the spectrum. All in all this DNS shows that the theory of scalar decay developed for the Kraichnan model and its generalization, in § 4.3, is borne out to a reasonable extent.