1 Introduction
Turbulence has a stochastic nature, so that the randomness and complexity of turbulence make theoretical analysis difficult. However, such randomness may permit the statistical description of turbulence, and it is believed that there is a certain kind of universality at sufficiently small scales in fully developed turbulence away from boundaries (Kolmogorov Reference Kolmogorov1941a,Reference Kolmogorovb). A striking feature of such universality is the existence of exact statistical laws for homogeneous isotropic turbulence at sufficiently high Reynolds numbers. Starting with Kraichnan (Reference Kraichnan1959), much effort has been dedicated to deriving the statistical theory of turbulence using various field-theoretic approaches from first principles such as the Navier–Stokes equation (Leslie Reference Leslie1973; Monin & Yaglom Reference Monin and Yaglom1975; McComb Reference McComb1989, Reference McComb2014). The goal is to accurately describe and predict statistical quantities of turbulence while maintaining a strong connection to the underlying dynamics of the Navier–Stokes equation. However, it has been long known that the equations for any statistical quantities are never closed because of the nonlinearity of the Navier–Stokes equation. The so-called closure problem has remained as a fundamental unsolved problem of fluid mechanics despite much effort by many researchers (Leslie Reference Leslie1973; Monin & Yaglom Reference Monin and Yaglom1975; Davidson Reference Davidson2004; McComb Reference McComb1989, Reference McComb2014). The closure related to the present study belongs to the so-called spectral closure theory that postulates a dynamic equation for the energy transfer function and its spirit is to address in theory the characteristics of different scales that constitute a turbulent field. A major concern of the spectral closure theory of turbulence is a derivation of the second-order moments such as spectrum from the Navier–Stokes equation without any adjustable parameters.
As one of the great achievements in the statistical theory of turbulence, the direct interaction approximation (DIA) will be highlighted. In DIA, the two-point two-time correlation and the response function are introduced, and the nonlinear interactions are involved in the response function. This yields the non-negativity of the energy spectrum. Although DIA fails in the derivation of Kolmogorov’s 
$-5/3$ law, DIA is attractive in the sense that there are several ways to derive the DIA equations, e.g. the weak dependence principle (Kraichnan Reference Kraichnan1959), self-consistent method (Herring Reference Herring1965) and renormalized perturbation method (Kraichnan Reference Kraichnan1977). As a reason for the failure of DIA, the Eulerian coordinate system is unsuited to distinguishing between two possible effects such as uniform convection and distortion. To maintain Galilean invariance, the Lagrangian coordinate system is suited to distinguishing between the above effects, and Lagrangian two-point two-time spectral closures have succeeded in the derivation of Kolmogorov’s 
$-5/3$ law (Kraichnan Reference Kraichnan1965; Kraichnan & Herring Reference Kraichnan and Herring1978; Kaneda Reference Kaneda1981). Kolmogorov’s 
$-5/3$ law has been observed in various turbulent motions, the interstellar medium (Armstrong, Rickett & Spangler Reference Armstrong, Rickett and Spangler1995), solar wind (Matthaeus & Goldstein Reference Matthaeus and Goldstein1982) and atmosphere (Nastrom & Gage Reference Nastrom and Gage1985; Tsuji Reference Tsuji2004). Among various closures, abridged Lagrangian history DIA (ALHDIA), strain-based ALHDIA (SBALHDIA) and Lagrangian renormalized approximation (LRA) are representative closures which can provide Kolmogorov’s 
$-5/3$ law through systematic ways without any adjustable parameters. As an exception to closures in the sense of the Eulerian viewpoint, local energy transfer is highlighted (McComb & Shanmugasundaram Reference McComb and Shanmugasundaram1984). Local energy transfer is derived from a purely Eulerian viewpoint using a mapping function, whereas the result is compatible with Kolmogorov’s 
$-5/3$ law. The Kolmogorov constant 
$K_{0}$ for an infinite Reynolds number is 1.77 and 1.72 for ALHDIA and LRA, respectively. On the other hand, 
$K_{0}$ for a finite Reynolds number is 1.78, 2.0, 1.67, 1.69 and 2.3 for ALHDIA, SBALHDIA, LRA, Markovianized LRA and local energy transfer, respectively (Herring & Kraichnan Reference Herring and Kraichnan1979; Gotoh, Kaneda & Bekki Reference Gotoh, Kaneda and Bekki1988; McComb & Shanmugasundaram Reference McComb and Shanmugasundaram1984). Taking into account the facts that 
$K_{0}=1.62\pm 0.17$ (Sreenivasan Reference Sreenivasan1995) in experiments and 
$K_{0}=1.8\pm 0.1$ for state-of-the-art high-resolution direct numerical simulation (DNS) (Ishihara et al. Reference Ishihara, Morishita, Yokokawa, Uno and Kaneda2016), it is difficult to say which closure is superior. On the other hand, because of facts such as (i) the so-called DIA families have been mainly limited to the exploration of universal constants of extremely large-Reynolds-number turbulent flow since the closed equation is a complex integro-differential equation having memory effect and (ii) there has recently been a growing tendency to study the effect of finite Reynolds number on the Kolmogorov theory since the clear Kolmogorov’s 
$-5/3$ law is not observed even for state-of-the-art high-resolution DNS (Ishihara et al. Reference Ishihara, Morishita, Yokokawa, Uno and Kaneda2016), practical closure, e.g. eddy-damped quasi-normal Markovianization (EDQNM), has been used in a wide range of turbulence studies (Sagaut & Cambon Reference Sagaut and Cambon2008). One of the advantages of EDQNM is to deepen our understanding of the phenomenology of turbulence including Reynolds number dependency (Bos, Clark & Rubinstein Reference Bos, Clark and Rubinstein2007a; Bos, Shao & Bertoglio Reference Bos, Shao and Bertoglio2007b), decaying law regardless of Reynolds and Schmidt number (Lesieur & Ossia Reference Lesieur and Ossia2000; Briard et al. Reference Briard, Gomez, Sagaut and Memari2015) and the role of the nonlinear term in comparison with linear theory (Sagaut & Cambon Reference Sagaut and Cambon2008). On the other hand, the disadvantage of EDQNM is that the results quantitatively depend on adjustable parameters (Bos & Fang Reference Bos and Fang2015). Therefore, it is of great importance for the further development of turbulence theory to develop closure theory derived from the Navier–Stokes equation without the introduction of any adjustable parameters or taking phenomenology into account for closure itself. In this paper, we extend Kaneda’s LRA, which has been recognized as the most sophisticated closure (Sagaut & Cambon Reference Sagaut and Cambon2008), to two-point single-time spectral closure under a few assumptions and legitimate mathematical procedures.
The present paper is organized as follows. In § 2.1, LRA is briefly summarized. The mathematical procedure for the extension to single-time spectral closure is described in § 2.2. In § 3, the numerical methods for closure and DNS are described and the numerical results for single- and two-point statistics are shown in § 4. Discussion and conclusions are presented in § 5.
2 Spectral closure theory
2.1 Lagrangian renormalized approximation
Here, we briefly summarize the LRA (more details can be found in Kaneda (Reference Kaneda1981) or Kida & Goto (Reference Kida and Goto1997)). The equations of incompressible fluid are
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}+(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{u}=-\unicode[STIX]{x1D735}p+\unicode[STIX]{x1D708}\unicode[STIX]{x0394}\boldsymbol{u}+\boldsymbol{f}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{f}=0, & \displaystyle\end{eqnarray}$$
$\boldsymbol{u}$ is velocity, 
$p$ is pressure normalized by density, 
$\boldsymbol{f}$ is stirring force and 
$\unicode[STIX]{x1D708}$ is kinematic viscosity. Kaneda (Reference Kaneda1981) introduced the Lagrangian position function defined by 
$$\begin{eqnarray}\unicode[STIX]{x1D713}(\boldsymbol{y},t|\boldsymbol{x},s)=\unicode[STIX]{x1D6FF}(\boldsymbol{y}-\boldsymbol{z}(\boldsymbol{x},s|t)),\end{eqnarray}$$ where 
$\unicode[STIX]{x1D6FF}$ is the three-dimensional Dirac delta function and 
$\boldsymbol{z}(\boldsymbol{x},s|t)$ is the position of a fluid element at time 
$t$ which passes 
$\boldsymbol{x}$ at time 
$s$. The Lagrangian position function obeys
$$\begin{eqnarray}\displaystyle & \displaystyle \left(\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}+\boldsymbol{u}(\boldsymbol{y},t)\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{y}\right)\unicode[STIX]{x1D713}(\boldsymbol{y},t|\boldsymbol{x},s)=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D713}(\boldsymbol{y},t|\boldsymbol{x},t)=\unicode[STIX]{x1D6FF}(\boldsymbol{x}-\boldsymbol{y}). & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\boldsymbol{v}(\boldsymbol{x},s|t)=\int \,\text{d}\boldsymbol{y}\unicode[STIX]{x1D713}(\boldsymbol{y},t|\boldsymbol{x},s)\boldsymbol{u}(\boldsymbol{y},t),\end{eqnarray}$$and from (2.1a), (2.1b) and (2.3a), it obeys
$$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{v}(\boldsymbol{x},s|t)}{\unicode[STIX]{x2202}t} & = & \displaystyle \int \,\text{d}\boldsymbol{y}\unicode[STIX]{x1D713}(\boldsymbol{y},t|\boldsymbol{x},s)\boldsymbol{f}(\boldsymbol{y},t)+\unicode[STIX]{x1D708}\int \,\text{d}\boldsymbol{y}\unicode[STIX]{x1D713}(\boldsymbol{y},t|\boldsymbol{x},s)\unicode[STIX]{x1D6E5}_{y}\boldsymbol{u}(\boldsymbol{y},t)\nonumber\\ \displaystyle & & \displaystyle -\,\int \,\text{d}\boldsymbol{y}\unicode[STIX]{x1D713}(\boldsymbol{y},t|\boldsymbol{x},s)\unicode[STIX]{x1D735}_{y}p(\boldsymbol{y},t).\end{eqnarray}$$ The closure problem considered here is to determine the subsequent statistical second-order moments such as 
$\langle \boldsymbol{u}(\boldsymbol{x},t)\boldsymbol{u}(\boldsymbol{y},t)\rangle$ and 
$\langle \boldsymbol{v}(\boldsymbol{x},s|t)\boldsymbol{u}(\boldsymbol{y},s)\rangle$, where 
$\langle \cdots \rangle$ represents an ensemble average. In LRA, the two-point two-time Lagrangian velocity correlation 
$\boldsymbol{Q}$ and Lagrangian velocity response function 
$\boldsymbol{G}$ are chosen as representatives. Here 
$\boldsymbol{Q}$ and 
$\boldsymbol{G}$ are defined as
$$\begin{eqnarray}\displaystyle & \displaystyle Q_{ij}(\boldsymbol{x},t|\boldsymbol{x}^{\prime },s)\equiv \langle P_{i\unicode[STIX]{x1D6FC}}(\unicode[STIX]{x1D735}_{x})v_{\unicode[STIX]{x1D6FC}}(\boldsymbol{x},s|t)u_{j}(\boldsymbol{x}^{\prime },s)\rangle ,\quad t\geqslant s, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle G_{ij}(\boldsymbol{x},t|\boldsymbol{x}^{\prime },s)\equiv \left\langle \frac{P_{i\unicode[STIX]{x1D6FC}}(\unicode[STIX]{x1D735}_{x})\unicode[STIX]{x1D6FF}v_{\unicode[STIX]{x1D6FC}}(\boldsymbol{x},s|t)}{\unicode[STIX]{x1D6FF}f_{j}(\boldsymbol{x}^{\prime },s)}\right\rangle ,\quad t\geqslant s, & \displaystyle\end{eqnarray}$$
$\unicode[STIX]{x1D6FF}$ represents the functional derivative. The operator 
$\boldsymbol{P}$ plays a role of projecting the arbitrary vector field with respect to 
$\boldsymbol{x}$ onto the solenoidal field. This operator plays an important role in LRA to eliminate the irrelevant terms in the expansion. In homogeneous turbulence, the closure theory is formulated in terms of the Fourier representation of (2.1a)–(2.5). After lengthy algebra (the details are shown in appendix A), we have 
$$\begin{eqnarray}\displaystyle & \displaystyle \left(\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}+2\unicode[STIX]{x1D708}k^{2}\right)Q(k,t)=S(k,t)+Q_{f}(k,t), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \left(\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D708}k^{2}+\unicode[STIX]{x1D702}(k,t,s)\right)Q(k,t,s)=0,\quad t>s, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \left(\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D708}k^{2}+\unicode[STIX]{x1D702}(k,t,s)\right)G(k,t,s)=0,\quad t>s, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle G(k,t,t)=1, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}S(k,t)=2\unicode[STIX]{x03C0}\int _{t_{0}}^{t}\,\text{d}s\iint _{\unicode[STIX]{x1D6E5}}\,\text{d}p\,\text{d}qkpqb_{kpq}Q(q,t,s)[G(k,t,s)Q(p,t,s)-Q(k,t,s)G(p,t,s)],\end{eqnarray}$$
$$\begin{eqnarray}\unicode[STIX]{x1D702}(k,t,s)=\unicode[STIX]{x03C0}\iint _{\unicode[STIX]{x1D6E5}}\,\text{d}p\,\text{d}qkpq(1-y^{2})(1-z^{2})\int _{s}^{t}\,\text{d}s^{\prime }Q(q,t,s^{\prime }),\end{eqnarray}$$ where 
$t_{0}$ is the initial time; 
$Q(k,t)=Q(k,t,t)=E(k,t)/(2\unicode[STIX]{x03C0}k^{2})$; 
$Q_{F}(k,t)$ is the forcing spectrum; 
$E(k)$ is the energy spectrum; 
$\iint _{\unicode[STIX]{x1D6E5}}$ denotes the integration over all regions of the 
$p$–
$q$ plane such that the three wavenumber vectors 
$\boldsymbol{k}$, 
$\boldsymbol{p}$ and 
$\boldsymbol{q}$ form a triangle; 
$x$, 
$y$ and 
$z$ are the cosines of the angles opposite to 
$k$, 
$p$ and 
$q$ in this triangle; and 
$b_{kpq}$ is the geometrical factor given by 
$b_{kpq}=p/k(xy+z^{3})$. The LRA satisfies the fluctuation–dissipation relation 
$Q(k,t,s)=G(k,t,s)Q(k,s)$ for 
$t\geqslant s$. In LRA, an eddy damping term 
$\unicode[STIX]{x1D702}(k,t,s)$ comes from the Lagrangian acceleration of the pressure. In (2.7b) and (2.7c), contributions from the random force are neglected because they are not important for 
$t>s$ in the same way as Gotoh, Nagaki & Kaneda (Reference Gotoh, Nagaki and Kaneda2000) for passive scalar.
2.2 Extension to the single-time Markovianized LRA
When turbulence is quasi-stationary in the sense that the decay with respect to 
$t$ of 
$Q(k,t,t)$ is much slower than that of 
$G(k,t,s)$ for 
$s$, then Markovianization 
$Q(k,t,s)=Q(k,t)G(k,t,s)$ is applicable (Gotoh et al. Reference Gotoh, Kaneda and Bekki1988; Gotoh & Kaneda Reference Gotoh and Kaneda2000). Applying Markovianization in (2.7a), (2.7c) and (2.7d), we have the Markovianized LRA (MLRA) equations:
$$\begin{eqnarray}\left(\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}+2\unicode[STIX]{x1D708}k^{2}\right)Q(k,t)=2\unicode[STIX]{x03C0}\iint _{\unicode[STIX]{x1D6E5}}\,\text{d}p\,\text{d}qkpqb_{kpq}\unicode[STIX]{x1D703}_{kpq}(t)Q(q,t)[Q(p,t)-Q(k,t)]+Q_{F}(k,t),\end{eqnarray}$$
$$\begin{eqnarray}\left(\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}}+\unicode[STIX]{x1D708}k^{2}+\unicode[STIX]{x1D702}(k,t,\unicode[STIX]{x1D70F})\right)G(k,t,\unicode[STIX]{x1D70F})=0,\quad \unicode[STIX]{x1D70F}>0,\end{eqnarray}$$
$$\begin{eqnarray}G(k,t,0)=1,\end{eqnarray}$$in which
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D703}_{kpq}(t)=\int _{0}^{t-t_{0}}\,\text{d}sG(k,t,s)G(p,t,s)G(q,t,s), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D702}(k,t,\unicode[STIX]{x1D70F})=\unicode[STIX]{x03C0}\iint _{\unicode[STIX]{x1D6E5}}\,\text{d}p\,\text{d}qkpq(1-y^{2})(1-z^{2})Q(q,t)\int _{0}^{\unicode[STIX]{x1D70F}}\,\text{d}s^{\prime }G(q,t,s^{\prime }), & \displaystyle\end{eqnarray}$$ where note that the evolution of 
$G$ is with respect to time difference 
$\unicode[STIX]{x1D70F}=t-s$. Using (2.10b), (2.10c) and (2.12), when 
$\unicode[STIX]{x1D70F}\rightarrow 0$, the Lagrangian velocity response function is expanded as
$$\begin{eqnarray}G(k,t,\unicode[STIX]{x1D70F})=c_{0}(k,t)+c_{1}(k,t)\unicode[STIX]{x1D70F}+\frac{c_{2}(k,t)}{2}\unicode[STIX]{x1D70F}^{2}+\cdots \,,\end{eqnarray}$$where
$$\begin{eqnarray}\displaystyle & \displaystyle c_{0}(k,t)=G(k,t,0)=1, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle c_{1}(k,t)=\left.\frac{\unicode[STIX]{x2202}G(k,t,\unicode[STIX]{x1D70F})}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}}\right|_{\unicode[STIX]{x1D70F}=0}=-\unicode[STIX]{x1D708}k^{2}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle c_{2}(k,t)=\left.\frac{\unicode[STIX]{x2202}^{2}G(k,t,\unicode[STIX]{x1D70F})}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}^{2}}\right|_{\unicode[STIX]{x1D70F}=0}=(\unicode[STIX]{x1D708}k^{2})^{2}-\unicode[STIX]{x1D707}(k,t), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\unicode[STIX]{x1D707}(k,t)=\iint _{\unicode[STIX]{x1D6E5}}\,\text{d}p\,\text{d}q\unicode[STIX]{x03C0}kpq(1-y^{2})(1-z^{2})Q(q,t)=2\unicode[STIX]{x03C0}\int _{0}^{\infty }\,\text{d}qkq^{3}Q(q,t)J\left(\frac{q}{k}\right),\end{eqnarray}$$ where the latter expression is derived from the exact integral on 
$p$, and 
$J(x)$ is given by
$$\begin{eqnarray}J(x)=\left[\frac{(1-x)^{4}(1+x)^{4}}{32x^{4}}\log \frac{1+x}{|1-x|}-\frac{1+x^{2}}{48x^{3}}(3x^{4}-14x^{2}+3)\right].\end{eqnarray}$$ It is noted that 
$J(x)=J(1/x)$ and 
$J(x\rightarrow 0)=8/15x+O(x^{3})$. As mentioned by Kaneda (Reference Kaneda1993) and Gotoh & Kaneda (Reference Gotoh and Kaneda2000), the short-time behaviours of 
$Q$ and 
$G$ have important information about the statistical quantities. For the long-time behaviour of 
$G$, we assume
$$\begin{eqnarray}G(k,t,\unicode[STIX]{x1D70F})\simeq \exp \left[-C_{1}(k,t)\unicode[STIX]{x1D70F}-\frac{C_{2}(k,t)}{2}\unicode[STIX]{x1D70F}^{2}\right]\equiv {\mathcal{G}}(k,t,\unicode[STIX]{x1D70F}),\end{eqnarray}$$ where unknown coefficients 
$C_{1}$ and 
$C_{2}$ are chosen to match asymptotic behaviours of (2.17) at 
$\unicode[STIX]{x1D70F}\rightarrow 0$ with the exact short-time behaviours (2.14b) and (2.14c). We then obtain
$$\begin{eqnarray}C_{1}(k,t)=\unicode[STIX]{x1D708}k^{2},\quad C_{2}(k,t)=\unicode[STIX]{x1D707}(k,t).\end{eqnarray}$$ This choice with respect to 
${\mathcal{G}}$ is equivalent to solving the following differential equation:
$$\begin{eqnarray}\left(\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}}+\unicode[STIX]{x1D708}k^{2}+\unicode[STIX]{x1D707}(k,t)\unicode[STIX]{x1D70F}\right){\mathcal{G}}(k,t,\unicode[STIX]{x1D70F})=0,\quad \unicode[STIX]{x1D70F}>0.\end{eqnarray}$$ Using (2.17) and (2.18), 
$\unicode[STIX]{x1D703}_{kpq}(t)$ in the nonlinear term can be calculated in a straightforward way as follows:
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D703}_{kpq}(t) & \simeq & \displaystyle \int _{0}^{t-t_{0}}\,\text{d}s{\mathcal{G}}(k,t,s){\mathcal{G}}(p,t,s){\mathcal{G}}(q,t,s)\nonumber\\ \displaystyle & = & \displaystyle \sqrt{\frac{\unicode[STIX]{x03C0}}{2d_{kpq}(t)}}\exp \left(\frac{c_{kpq}^{2}}{2d_{kpq}(t)}\right)\nonumber\\ \displaystyle & & \displaystyle \times \,\left[\text{erf}\left(\frac{c_{kpq}}{\sqrt{2d_{kpq}(t)}}+\sqrt{\frac{d_{kpq}(t)}{2}}(t-t_{0})\right)-\text{erf}\left(\frac{c_{kpq}}{\sqrt{2d_{kpq}(t)}}\right)\right]\!,\qquad\end{eqnarray}$$with
$$\begin{eqnarray}\displaystyle & c_{kpq}=\unicode[STIX]{x1D708}(k^{2}+p^{2}+q^{2}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & d_{kpq}=\unicode[STIX]{x1D707}(k,t)+\unicode[STIX]{x1D707}(p,t)+\unicode[STIX]{x1D707}(q,t). & \displaystyle\end{eqnarray}$$
$t\gg t_{0}$ or 
$t_{0}=-\infty$, 
$\unicode[STIX]{x1D703}_{kpq}(t)$ can be simplified as 
$$\begin{eqnarray}\unicode[STIX]{x1D703}_{kpq}(t)=\sqrt{\frac{\unicode[STIX]{x03C0}}{2d_{kpq}(t)}}\exp \left(\frac{c_{kpq}^{2}}{2d_{kpq}(t)}\right)\text{erfc}\left(\frac{c_{kpq}}{\sqrt{2d_{kpq}(t)}}\right).\end{eqnarray}$$ In the numerical calculation, there are various rational approximations for the error function (Abramowitz & Stegun Reference Abramowitz and Stegun1964). We, however, calculated it by means of a built-in function. For 
$x=c_{kpq}/\sqrt{(2d_{kpq}(t))}\rightarrow \infty$, the integrand is calculated using the asymptotic property of the error function (Gradshteyn & Ryzhik Reference Gradshteyn and Ryzhik2014) as follows:
$$\begin{eqnarray}\unicode[STIX]{x1D703}_{kpq}(t)=\sqrt{\frac{1}{2d_{kpq}}}\frac{1}{x}\mathop{\sum }_{k=0}^{\infty }(-1)^{k}\frac{(2k-1)!!}{(2x^{2})^{k}}.\end{eqnarray}$$Thus, equations (2.20), (2.21a) and (2.21b) are used to calculate the nonlinear term in the single-time MLRA. The computational cost of the single-time MLRA by this approximation is of the same order as that of the EDQNM, test field model (Kraichnan Reference Kraichnan1971) and Lagrangian Markovianized field approximation (Bos & Bertoglio Reference Bos and Bertoglio2013). The present method is also promising for anisotropic turbulence. In EDQNM, an eddy damping term is derived phenomenologically on the dimensional ground and anisotropic effects are not incorporated in an eddy damping term. Thus, the single-time MLRA is more sophisticated compared to the other spectral closures since it can calculate the statistical quantities by theory without introduction of any adjustable parameters.
3 Numerical method
Hereinafter, the single-time MLRA is termed closure for the sake of simplicity. In closure, we solved (2.10a) with (2.22). The time derivative was estimated with a first-order forward time scheme. The viscous, energy transfer and forcing terms are estimated explicitly. We used 
$Q_{F}\sim \exp (-k^{2}/k_{f}^{2})$ in this study and we have confirmed that the form of forcing spectrum is negligible for statistical quantities, where 
$k_{f}=2.5$. The wavenumbers were discretized logarithmically, namely 
$k_{i}=k_{min}\times 2^{i-1/F}$, where 
$k_{min}=1$, 
$F=8$ and 
$i\in [1,\ldots ,N]$. Here 
$N$ was chosen to satisfy 
$k_{max}\unicode[STIX]{x1D702}\geqslant 2$ at the statistical steady state. Statistical quantities were evaluated when a criterion 
$\sqrt{(K(t)-K(t-\text{d}t))^{2}}/K(t)\leqslant 1\times 10^{-6}$ is satisfied, where 
$\text{d}t$ is the time step and 
$K(t)$ is the turbulence kinetic energy (TKE) at time 
$t$.
In DNS, the spectral method was used with periodic boundary conditions of periods of 
$2\unicode[STIX]{x03C0}$ in each of the three Cartesian coordinate directions. For time integration, we used the fourth-order Runge–Kutta scheme and the so-called phase shift method was used for de-aliasing, in which the maximum wavenumber of the retained Fourier modes is about 
$\sqrt{2}/3N$, where 
$N^{3}$ is the number of grid points. In this study, we force the turbulence as
$$\begin{eqnarray}\boldsymbol{f}(\boldsymbol{k},t)=\frac{\unicode[STIX]{x1D716}(t)-(K(t)-K_{\infty })/\unicode[STIX]{x1D70F}_{f}}{2K_{f}(t)}\hat{\boldsymbol{u}}(\boldsymbol{k},t)[H(k)-H(k-k_{f})],\end{eqnarray}$$ where 
$K_{\infty }$ is the ideal TKE, 
$\unicode[STIX]{x1D70F}_{f}$ is the time scale to control forcing, 
$K_{f}$ is the TKE composed of wavenumbers lower than forced wavenumber 
$k_{f}=2.5$ and 
$H$ is the Heaviside step function. Here, the parameters were set as 
$K_{\infty }=1/2$ and 
$\unicode[STIX]{x1D70F}_{f}=5dt=O(\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}})$, where 
$\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}$ is the Kolmogorov time scale. Under this condition, TKE is well controlled, e.g. the value is 
$1/2$ with standard deviation 
$O(10^{-7})$ for lowest Reynolds number DNS and with 
$O(10^{-5})$ for highest Reynolds number DNS. Here, 
$E(k,0)\sim k^{4}\exp (-k^{2}/k_{L}^{2})$ with 
$k_{L}=2.5$ was used for the initial energy spectrum in all DNSs. We started all runs with the same initial energy spectrum instead of using the final state of the run with the smaller 
$N$ as the initial state of the new runs as done by Kaneda et al. (Reference Kaneda, Ishihara, Yokokawa, Itakura and Uno2003). After 5 eddy-turnover times based on the initial energy spectrum, we carried out DNSs a further 10 eddy-turnover times for calculating the single-point statistics. For two-point statistics, the statistics were calculated every eddy-turnover time, i.e. 10 data gatherings in total. In the present code, hybrid parallelization was implemented using the message-passing interface library and OpenMP.
Normalized TKE dissipation rate 
$C_{\unicode[STIX]{x1D716}}$ versus turbulence Reynolds number 
$Re_{\unicode[STIX]{x1D706}}$ with 
$\pm 3\unicode[STIX]{x1D70E}$ uncertainty error bars for present DNS, where 
$\unicode[STIX]{x1D70E}$ is the standard deviation. Also included are DNS data from Jiménez et al. (Reference Jiménez, Wray, Saffman and Rogallo1993), Cao, Chen & Doolen (Reference Cao, Chen and Doolen1999), Gotoh, Fukayama & Nakano (Reference Gotoh, Fukayama and Nakano2002) and Kaneda et al. (Reference Kaneda, Ishihara, Yokokawa, Itakura and Uno2003).
Turbulence characteristics in the numerical simulations.
4 Numerical results
Turbulence characteristics are summarized in table 1. Here, 
$K$, 
$\unicode[STIX]{x1D716}$, 
$L$, 
$\unicode[STIX]{x1D706}$ and 
$\unicode[STIX]{x1D702}$ are the TKE, TKE dissipation rate, integral length scale, Taylor microscale and Kolmogorov scale, respectively. These variables are, in spectral space, defined as 
$K=\int _{0}^{\infty }\,\text{d}kE(k)$, 
$\unicode[STIX]{x1D716}=2\unicode[STIX]{x1D708}\int _{0}^{\infty }\,\text{d}kk^{2}E(k)$, 
$L=3\unicode[STIX]{x03C0}/(4K)\int _{0}^{\infty }\,\text{d}kE(k)/k$, 
$\unicode[STIX]{x1D706}^{2}=10\unicode[STIX]{x1D708}K/\unicode[STIX]{x1D716}$ and 
$\unicode[STIX]{x1D702}=(\unicode[STIX]{x1D708}^{3}/\unicode[STIX]{x1D716})^{1/4}$, respectively (Davidson Reference Davidson2004). Here 
$Re_{\unicode[STIX]{x1D706}}$ is the turbulence Reynolds number based on the Taylor microscale 
$Re_{\unicode[STIX]{x1D706}}\equiv \sqrt{2K/3}\unicode[STIX]{x1D706}/\unicode[STIX]{x1D708}$ and 
$C_{\unicode[STIX]{x1D716}}$ is the normalized TKE dissipation rate 
$C_{\unicode[STIX]{x1D716}}\equiv \unicode[STIX]{x1D716}L/(2K/3)^{3/2}$. Figure 1 shows the relation between 
$C_{\unicode[STIX]{x1D716}}$ and 
$Re_{\unicode[STIX]{x1D706}}$ obtained by closure and DNS together with DNS data of other researchers. The functional form (Doering & Foias Reference Doering and Foias2002) 
$C_{\unicode[STIX]{x1D716}}=a[1+\sqrt{1+(b/Re_{\unicode[STIX]{x1D706}})^{2}}]$ with 
$a=0.217$ and 
$b=88.65$ agrees with the present results obtained by closure and DNS. The parameters were evaluated by means of the Levenberg–Marquardt least-squares method for closure results. The fitting curve is within the acceptable range of 
$\pm 3$ standard deviations for DNS data.
Normalized energy spectra for (
$a$) DNS and (
$b$) closure and compensated energy spectra for (
$c$) DNS and (
$d$) closure. Insets in (
$c$) and (
$d$) show the log-linear plots. The dashed grey line in (
$a$) is a closure result for the highest Reynolds number. In (
$c$), DNS data only for the top four highest Reynolds numbers are shown to emphasize the spectral bump.
Figure 2 shows the normalized energy spectra and compensated energy spectra. As shown in figure 2, the spectra collapse well in the inertial subrange and dissipation range for both DNS and closure. The differences between closure and DNS are mainly observed in the energy-containing range due to the different forcing function and the degree of spectral bump. However, the energy spectra obtained in closure qualitatively agree with those of DNS. As shown in figure 2(
$b$), a slight bump is observed in the closure; however, its degree is suppressed compared to EDQNM (Lesieur & Ossia Reference Lesieur and Ossia2000). The estimated Kolmogorov constant is 
$K_{0}=1.51$ for the closure and is slightly smaller than values in DNSs and theoretical estimation of LRA (
$K_{0}=1.72$). Figure 3 shows the energy flux 
$\unicode[STIX]{x1D6F1}(k)$ normalized by 
$\unicode[STIX]{x1D716}$, where energy flux is defined as 
$\unicode[STIX]{x1D6F1}(k)=\int _{k}^{\infty }\,\text{d}kT(k)=-\int _{0}^{k}\,\text{d}kT(k)$, where 
$T(k)$ is the energy transfer function. For high-Reynolds-number turbulent flows, 
$\unicode[STIX]{x1D6F1}(k)=\unicode[STIX]{x1D716}$ is satisfied in the inertial subrange according to Kolmogorov (Reference Kolmogorov1941a). As shown in figure 3, 
$\unicode[STIX]{x1D6F1}(k)=\unicode[STIX]{x1D716}$ in the inertial subrange is satisfied for high-Reynolds-number cases in both closure and DNS. The difference between closure and DNS is mainly observed around the energy-containing and bump ranges.
Normalized energy fluxes for (
$a$) DNS and (
$b$) closure. For a description, refer to figure 2.
Normalized second-order structure functions for (
$a$) DNS and (
$b$) closure and compensated second-order structure functions for (
$c$) DNS and (
$d$) closure. For a description, refer to figure 2. Here, 
$v_{\unicode[STIX]{x1D702}}=(\unicode[STIX]{x1D708}\unicode[STIX]{x1D716})^{1/4}$ is the Kolmogorov velocity.
Normalized third-order structure functions for (
$a$) DNS and (
$b$) closure and skewness functions for (
$c$) DNS and (
$d$) closure. For a description, refer to figure 2.
The spatial second- and third-order structure functions are given by
$$\begin{eqnarray}S_{2}(r)=4\int _{0}^{\infty }\,\text{d}kE(k)f(kr),\quad f(x)=1-\frac{\sin x-x\cos x}{x^{3}},\end{eqnarray}$$
$$\begin{eqnarray}S_{3}(r)=12r\int _{0}^{\infty }\,\text{d}kT(k)g(kr),\quad g(x)=\frac{3\sin x-3x\cos x-x^{2}\sin x}{x^{5}}.\end{eqnarray}$$ Here 
$S_{2}(r)\sim (\unicode[STIX]{x1D716}r)^{2/3}$ and 
$S_{3}(r)=-4/5\unicode[STIX]{x1D716}r$ are believed to be satisfied in the inertial subrange for high-Reynolds-number turbulent flows. These expressions in physical space are equivalent to the 
$-5/3$ law of energy spectrum and 
$\unicode[STIX]{x1D6F1}(k)=\unicode[STIX]{x1D716}$ in wavenumber space. Using Taylor expansions, it is found that 
$S_{2}(r\rightarrow 0)=\unicode[STIX]{x1D716}r^{2}/15\unicode[STIX]{x1D708}+O(r^{4})$ and 
$S_{3}(r\rightarrow 0)=-{\textstyle \frac{2}{35}}r^{3}\int _{0}^{\infty }\,\text{d}kk^{2}T(k)+O(r^{5})$ for dissipation range. Figure 4 shows the second- and third-order structure functions and skewness function. The Kolmogorov constant in the second-order structure function is given as 
$C_{0}=27/55\unicode[STIX]{x1D6E4}(1/3)K_{0}$. The estimated 
$C_{0}$ in the present closure is 2.0, and its value is slightly lower than the DNS value. The Reynolds number dependency of 
$S_{3}(r\rightarrow 0)/S_{2}(r\rightarrow 0)^{3/2}=\langle (\unicode[STIX]{x2202}u/\unicode[STIX]{x2202}x)^{3}\rangle /\langle (\unicode[STIX]{x2202}u/\unicode[STIX]{x2202}x)^{2}\rangle ^{3/2}$ is observed for both closure and DNS as in Sreenivasan & Antonia (Reference Sreenivasan and Antonia1997). In the present closure, 
$S_{3}(r\rightarrow 0)/S_{2}(r\rightarrow 0)^{3/2}$ approaches 
$-0.6$ with an increase of 
$Re_{\unicode[STIX]{x1D706}}$ and is slightly smaller than for LRA, 
$S_{3}(r\rightarrow 0)/S_{2}(r\rightarrow 0)^{3/2}=-0.66$, for an infinite Reynolds number (Kaneda Reference Kaneda1993). Compared to the EDQNM (Bos et al. Reference Bos, Clark and Rubinstein2007a), the present closure is closer to DNS at moderate Reynolds numbers (Gotoh et al. Reference Gotoh, Fukayama and Nakano2002). It is noted that the present closure does not contradict the boundedness of the velocity derivative skewness for high-Reynolds-number turbulent flows (Antonia et al. Reference Antonia, Tang, Djenidi and Danaila2015).
As shown in figures 1–5, the present closure yields quantitative agreement with DNS and its extension to anisotropic turbulence will be performed in the future.
Eddy damping for (
$a$) present closure and (
$b$) EDQNM.
Before closing this paper, we discuss the eddy damping and briefly comment on the capability of the intermittency effect. In EDQNM, the response function is regarded as having the form 
$G(k,t,\unicode[STIX]{x1D70F})=\exp [-(\unicode[STIX]{x1D708}k^{2}+\unicode[STIX]{x1D707}^{EDQNM}(k,t))\unicode[STIX]{x1D70F}]$, where 
$\unicode[STIX]{x1D707}^{EDQNM}(k,t)=\unicode[STIX]{x1D6FC}\sqrt{\int _{0}^{k}\,\text{d}pp^{2}E(p,t)}$ is the inverse of the rotation time or that of strain time. In EDQNM, one can adjust 
$K_{0}$ through the parameter 
$\unicode[STIX]{x1D6FC}$ and its relation is given as 
$K_{0}\approx 2.76\unicode[STIX]{x1D6FC}^{2/3}$ (Lesieur & Ossia Reference Lesieur and Ossia2000). Figure 6 shows the eddy damping for the present closure and EDQNM. In the inertial subrange using Kolmogorov’s 
$-5/3$ law, 
$\unicode[STIX]{x1D707}(k)$ in the present closure takes the form 
$\unicode[STIX]{x1D707}(k)=1.06K_{0}\unicode[STIX]{x1D716}^{2/3}k^{4/3}$ (Kaneda Reference Kaneda1993), whereas 
$[\unicode[STIX]{x1D707}^{EDQNM}(k)]^{2}=3/4\unicode[STIX]{x1D6FC}^{2}K_{0}\unicode[STIX]{x1D716}^{2/3}k^{4/3}$ in EDQNM. In the present closure, 
$\unicode[STIX]{x1D707}_{k}$ is regarded as the inverse of the squared Lagrangian time scale, where the Lagrangian time scale in wavenumber space is given as 
$\unicode[STIX]{x1D716}^{-1/3}k^{-2/3}$ (Sagaut & Cambon Reference Sagaut and Cambon2008). The behaviours of energy-containing range and dissipation range can be understood by analysing the non-local interaction of eddy damping. For 
$k\ll k_{c}\leqslant p\sim q$, we assume that the contribution of integration 
$\iint _{\unicode[STIX]{x1D6E5}}\,\text{d}p\,\text{d}q$ mainly comes from 
$p\geqslant k_{c}$ and 
$q\geqslant k_{c}$, where 
$k_{c}$ is the convenient cut-off wavenumber. We then have
$$\begin{eqnarray}\unicode[STIX]{x1D707}(k\ll k_{c})=2\int _{k_{c}}^{\infty }\,\text{d}p\int _{p-k}^{p}\,\text{d}q\frac{kp}{q}(x+yz)^{2}E(q)\approx \frac{8}{15}k^{2}\int _{k_{c}}^{\infty }\,\text{d}pE(p),\end{eqnarray}$$ where the change of variable using 
$z=(k^{2}+p^{2}-q^{2})/(2kp)$ and Taylor expansions at 
$k/p\ll 1$ were used in the derivation. It is found from (4.2) that 
$\unicode[STIX]{x1D707}(k)\sim Kk^{2}$ at low wavenumber for the present closure, whereas 
$[\unicode[STIX]{x1D707}^{EDQNM}(k)]^{2}\sim k^{3}E(k)$ at low wavenumber for EDQNM. In this study, the forcing spectrum is given as 
$k^{2}$ at low wavenumber, so that 
$k^{5}$ behaviours can be seen for the eddy damping in EDQNM as shown in figure 6. In other words, eddy damping of the present closure is independent of the energy spectrum at low wavenumber. On the other hand, for 
$q\leqslant k_{c}\ll k\sim p$ (or 
$p\leqslant k_{c}\ll k\sim q$), we assume that the contribution of integration 
$\iint _{\unicode[STIX]{x1D6E5}}\,\text{d}p\,\text{d}q$ mainly comes from 
$q\leqslant k_{c}$ (or 
$p\leqslant k_{c}$). We then have
$$\begin{eqnarray}\unicode[STIX]{x1D707}(k\gg k_{c})=2\int _{0}^{k_{c}}\,\text{d}q\int _{k-q}^{k+q}\,\text{d}p\frac{kp}{q}(x+yz)^{2}E(q)\approx \frac{8}{15}\int _{0}^{k_{c}}\,\text{d}qq^{2}E(q),\end{eqnarray}$$ where the change of variable using 
$y=(k^{2}+q^{2}-p^{2})/(2kq)$ and Taylor expansions at 
$q/k\ll 1$ were used in the derivation. It is found from (4.3) that 
$\unicode[STIX]{x1D707}(k)$ approaches 
$4/15\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}^{2}$ for high wavenumber for the present closure, in which 
$\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}=\sqrt{\unicode[STIX]{x1D708}/\unicode[STIX]{x1D716}}$ is the Kolmogorov time scale. For EDQNM, 
$[\unicode[STIX]{x1D707}^{EDQNM}(k)]^{2}=\unicode[STIX]{x1D6FC}^{2}\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}^{2}/2$ for the high-wavenumber limit. Thus, the eddy damping in the present closure qualitatively agrees with squared eddy damping in EDQNM for high wavenumber, whereas there is a qualitative difference for low wavenumber.
It is known that the drawback of spectral closures is neglect of the intermittency effects. However, as discussed by Yoshida, Ishihara & Kaneda (Reference Yoshida, Ishihara and Kaneda2003), it is non-trivial whether or not the closure theories do not yield anomalous scaling since the present closure is similar in a sense to the exact closure equation for the Kraichnan model (Kraichnan Reference Kraichnan1994) which has a solution exhibiting anomalous scaling. Thus, the problem related to intermittency still remains unsolved and further studies are required.
5 Discussion and conclusion
In the present study, we derived the two-point single-time spectral closure starting from LRA and the use of assumptions (i) Markovianization and (ii) the Lagrangian velocity response function being expressed by 
$G(k,\unicode[STIX]{x1D70F})=\exp (-C_{1}(k)\unicode[STIX]{x1D70F}-C_{2}(k)\unicode[STIX]{x1D70F}^{2}/2)$. The unknown functions 
$C_{1}(k)$ and 
$C_{2}(k)$ are theoretically derived from being consistent with the exact short-time behaviour of 
$G(k,\unicode[STIX]{x1D70F})$ and the asymptotic short-time behaviour of assumed exponential form of 
$G(k,\unicode[STIX]{x1D70F})$. It is found that 
$C_{1}(k)$ is the inverse of the viscous time scale and 
$C_{2}(k)$ is the inverse of the squared Lagrangian time scale in the present closure. The main difference between the present closure and other single-time spectral closures (EDQNM, test field model and Lagrangian Markovianized field approximation) is in the form of the response function, where 
$G(k,\unicode[STIX]{x1D70F})=\exp (-C(k)\unicode[STIX]{x1D70F})$ is assumed for the other single-time spectral closures. In the inertial subrange, it is found that 
$C_{2}(k)$ and 
$C(k)^{2}$ are proportional to 
$\unicode[STIX]{x1D716}^{2/3}k^{4/3}$ regardless of closure. On the other hand, a difference can be seen for eddy damping in the energy-containing range.
More generally speaking, most closures proposed so far take the form of (2.8) for the energy transfer function in homogeneous isotropic turbulence and the differences in closures originate from the Lagrangian velocity response function and Lagrangian velocity correlation. From the perspective of the renormalization theory, these differences come from the different representatives, and a more suitable choice of representatives may improve the level of approximation. However, taking into account attractive points in LRA, e.g. (i) satisfactory nature of the fluctuation–dissipation theorem and (ii) the simple closed equations compared to ALHDIA and SBALHDIA, the present closure based on LRA is not only promising for an extension to closure in anisotropic turbulence but also reduces the computational cost as well as the other two-point single-time closures.
The numerical results show that the single-time MLRA has good agreement with DNS for single- and two-point statistics including the finite Reynolds number dependency. Thus, the single-time MLRA is more sophisticated compared to other two-point single-time spectral closures in the sense that (i) it can calculate the statistical quantities by theory without the introduction of any adjustable parameters, (ii) it involves simple closed equations and (iii) it shows quantitatively good agreements with DNS. An extension of the present closure to anisotropic turbulence will be performed in future work.
Acknowledgements
The author would like to thank colleagues Professor D. Sakaguchi (Nagasaki University), Professor H. Ueki (Nagasaki University) and emeritus Professor M. Ishida (Nagasaki University) for discussions. He would like to acknowledge Professor T. Ishihara (Okayama University) for informing him of the useful paper on closure and to thank Ms H. Kitamura (Nagasaki Nihon University JHS) for checking and correcting the grammatical errors. He is also grateful to the three anonymous referees for interesting discussions and helpful comments on this topic. The DNS was carried out using the computer resource offered by Research Institute for Information Technology, Kyushu University. This work was supported by the Japanese Ministry of Education, Science and Culture through grants-in-aid (no. 19K14891) and Harada Memorial Foundation.
Declaration of interests
The author reports no conflict of interest.
Appendix A
The Fourier coefficient of Lagrangian velocity is defined as
$$\begin{eqnarray}\hat{v}_{i}(\boldsymbol{k},s|t)=\frac{1}{(2\unicode[STIX]{x03C0})^{3}}\int v_{i}(\boldsymbol{x},s|t)\exp (-\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{x})\,\text{d}\boldsymbol{x}\end{eqnarray}$$and that of the Lagrangian position function is defined as
$$\begin{eqnarray}\hat{\unicode[STIX]{x1D713}}(\boldsymbol{k},t|\boldsymbol{k}^{\prime },s)=\frac{1}{(2\unicode[STIX]{x03C0})^{6}}\iint \unicode[STIX]{x1D713}(\boldsymbol{y},t|\boldsymbol{x},s)\exp (-\text{i}(\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{y}+\boldsymbol{k}^{\prime }\boldsymbol{\cdot }\boldsymbol{x}))\,\text{d}\boldsymbol{y}\,\text{d}\boldsymbol{x},\end{eqnarray}$$ with an initial condition 
$\hat{\unicode[STIX]{x1D713}}(\boldsymbol{k},s|\boldsymbol{k}^{\prime },s)=\unicode[STIX]{x1D6FF}(\boldsymbol{k}+\boldsymbol{k}^{\prime })/(2\unicode[STIX]{x03C0})^{3}$ as understood from (2.3b). Here 
$\text{i}=\sqrt{-1}$. From (2.4), (A 1) and (A 2), we have
$$\begin{eqnarray}\hat{v}_{i}(\boldsymbol{k},s|t)=(2\unicode[STIX]{x03C0})^{3}\int \,\text{d}\boldsymbol{k}^{\prime }\hat{u} _{i}(\boldsymbol{k}^{\prime },t)\hat{\unicode[STIX]{x1D713}}(-\boldsymbol{k}^{\prime },t|\boldsymbol{k},s),\end{eqnarray}$$ with an initial condition 
$\hat{v}_{i}(\boldsymbol{k},s|s)=\hat{u} _{i}(\boldsymbol{k},s)$. Hereinafter, 
$\hat{f}(\boldsymbol{k},t)$ represents the Fourier representation of 
$f(\boldsymbol{x},t)$. From (2.5) and (A 1), the equation for 
$\hat{\boldsymbol{v}}(\boldsymbol{k},s|t)$ is given as
$$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}\hat{v}_{i}(\boldsymbol{k},s|t) & = & \displaystyle (2\unicode[STIX]{x03C0})^{3}\!\int \,\text{d}\boldsymbol{k}^{\prime }\hat{f}_{i}(\boldsymbol{k}^{\prime },t)\hat{\unicode[STIX]{x1D713}}(-\boldsymbol{k}^{\prime },t|\boldsymbol{k},s)-\unicode[STIX]{x1D708}(2\unicode[STIX]{x03C0})^{3}\!\int \,\text{d}\boldsymbol{k}^{\prime }k^{\prime 2}\hat{u} _{i}(\boldsymbol{k}^{\prime },t)\hat{\unicode[STIX]{x1D713}}(-\boldsymbol{k}^{\prime },t|\boldsymbol{k},s)\nonumber\\ \displaystyle & & \displaystyle +\,\unicode[STIX]{x1D6FC}\text{i}(2\unicode[STIX]{x03C0})^{3}\int \,\text{d}\boldsymbol{k}^{\prime }\hat{\unicode[STIX]{x1D713}}(-\boldsymbol{k}^{\prime },t|\boldsymbol{k},s)\frac{k_{i}^{\prime }k_{j}^{\prime }k_{k}^{\prime }}{k^{\prime 2}}\int _{\boldsymbol{k}^{\prime }=\boldsymbol{p}+\boldsymbol{q}}\hat{u} _{j}(\boldsymbol{p},t)\hat{u} _{k}(\boldsymbol{q},t)\,\text{d}\boldsymbol{p},\end{eqnarray}$$ where 
$\unicode[STIX]{x1D6FC}=1$ is the bookkeeping parameter for the convenience of perturbation expansion and is introduced for the convolution terms. From (A 4), the equations of the Lagrangian velocity correlation 
$Q_{ij}(\boldsymbol{k},t,s)=P_{i\unicode[STIX]{x1D6FC}}(\boldsymbol{k})\hat{v}_{\unicode[STIX]{x1D6FC}}(\boldsymbol{k},s|t)\hat{v}_{j}(-\boldsymbol{k},s|s)$ and the Lagrangian velocity response function 
$G_{ij}=P_{i\unicode[STIX]{x1D6FC}}(\boldsymbol{k})\unicode[STIX]{x1D6FF}\hat{v}_{\unicode[STIX]{x1D6FC}}(\boldsymbol{k},s|t)/\unicode[STIX]{x1D6FF}\hat{f}_{j}(\boldsymbol{k}^{\prime },s)$ are given as
$$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}Q_{ij}(\boldsymbol{k},t,s) & = & \displaystyle \underbrace{(2\unicode[STIX]{x03C0})^{3}P_{i\unicode[STIX]{x1D6FC}}(\boldsymbol{k})\int \,\text{d}\boldsymbol{k}^{\prime }\hat{f}_{\unicode[STIX]{x1D6FC}}(\boldsymbol{k}^{\prime },t)\hat{u} _{j}(-\boldsymbol{k},s)\hat{\unicode[STIX]{x1D713}}(-\boldsymbol{k}^{\prime },t|\boldsymbol{k},s)}_{A}\nonumber\\ \displaystyle & & \displaystyle \underbrace{-\unicode[STIX]{x1D708}(2\unicode[STIX]{x03C0})^{3}P_{i\unicode[STIX]{x1D6FC}}(\boldsymbol{k})\int \,\text{d}\boldsymbol{k}^{\prime }k^{\prime 2}\hat{u} _{\unicode[STIX]{x1D6FC}}(\boldsymbol{k}^{\prime },t)\hat{u} _{j}(-\boldsymbol{k},s)\hat{\unicode[STIX]{x1D713}}(-\boldsymbol{k}^{\prime },t|\boldsymbol{k},s)}_{B}\nonumber\\ \displaystyle & & \displaystyle \underbrace{+\unicode[STIX]{x1D6FC}\text{i}(2\unicode[STIX]{x03C0})^{3}P_{i\unicode[STIX]{x1D6FC}}(\boldsymbol{k})\int \,\text{d}\boldsymbol{k}^{\prime }\hat{\unicode[STIX]{x1D713}}(-\boldsymbol{k}^{\prime },t|\boldsymbol{k},s)\frac{k_{\unicode[STIX]{x1D6FC}}^{\prime }k_{\unicode[STIX]{x1D6FD}}^{\prime }k_{\unicode[STIX]{x1D6FE}}^{\prime }}{k^{\prime 2}}\int _{\boldsymbol{k}^{\prime }=\boldsymbol{p}+\boldsymbol{q}}\hat{u} _{j}(-\boldsymbol{k},s)\hat{u} _{\unicode[STIX]{x1D6FD}}(\boldsymbol{p},t)\hat{u} _{\unicode[STIX]{x1D6FE}}(\boldsymbol{q},t)\,\text{d}\boldsymbol{p}}_{C},\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}G_{ij}(\boldsymbol{k},t|\boldsymbol{k}^{\prime },s) & = & \displaystyle \underbrace{(2\unicode[STIX]{x03C0})^{3}P_{i\unicode[STIX]{x1D6FC}}(\boldsymbol{k})\int \,\text{d}\boldsymbol{k}^{\prime \prime }\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D6FC}j}\unicode[STIX]{x1D6FF}(\boldsymbol{k}^{\prime \prime }+\boldsymbol{k}^{\prime })\unicode[STIX]{x1D6FF}(t-s)\unicode[STIX]{x1D713}(-\boldsymbol{k}^{\prime \prime },t|\boldsymbol{k},s)}_{D}\nonumber\\ \displaystyle & & \displaystyle \underbrace{-\unicode[STIX]{x1D708}(2\unicode[STIX]{x03C0})^{3}P_{i\unicode[STIX]{x1D6FC}}(\boldsymbol{k})\int \,\text{d}\boldsymbol{k}^{\prime \prime }k^{\prime \prime 2}G_{\unicode[STIX]{x1D6FC}j}^{E}(\boldsymbol{k}^{\prime \prime },t|\boldsymbol{k}^{\prime },s)\unicode[STIX]{x1D713}(-\boldsymbol{k}^{\prime \prime },t|\boldsymbol{k},s)}_{E}\nonumber\\ \displaystyle & & \displaystyle \underbrace{+\unicode[STIX]{x1D6FC}2\text{i}(2\unicode[STIX]{x03C0})^{3}P_{i\unicode[STIX]{x1D6FC}}(\boldsymbol{k})\int \,\text{d}\boldsymbol{k}^{\prime \prime }\hat{\unicode[STIX]{x1D713}}(-\boldsymbol{k}^{\prime \prime },t|\boldsymbol{k},s)\frac{k_{\unicode[STIX]{x1D6FC}}^{\prime \prime }k_{\unicode[STIX]{x1D6FD}}^{\prime \prime }k_{\unicode[STIX]{x1D6FE}}^{\prime \prime }}{k^{\prime \prime 2}}\int _{\boldsymbol{k}^{\prime \prime }=\boldsymbol{p}+\boldsymbol{q}}\,\text{d}\boldsymbol{p}\hat{u} _{\unicode[STIX]{x1D6FD}}(\boldsymbol{p},t)G_{\unicode[STIX]{x1D6FE}j}^{E}(\boldsymbol{q},t|\boldsymbol{k}^{\prime },s)}_{F}\nonumber\\ \displaystyle & & \displaystyle \underbrace{+(2\unicode[STIX]{x03C0})^{3}P_{i\unicode[STIX]{x1D6FC}}(\boldsymbol{k})\int \,\text{d}\boldsymbol{k}^{\prime }\hat{f}_{\unicode[STIX]{x1D6FC}}(\boldsymbol{k}^{\prime \prime },t)\unicode[STIX]{x1D6F9}_{j}(-\boldsymbol{k}^{\prime \prime },t|\boldsymbol{k},\boldsymbol{k}^{\prime },s)}_{G}\nonumber\\ \displaystyle & & \displaystyle \underbrace{-\unicode[STIX]{x1D708}(2\unicode[STIX]{x03C0})^{3}P_{i\unicode[STIX]{x1D6FC}}(\boldsymbol{k})\int \,\text{d}\boldsymbol{k}^{\prime \prime }\hat{u} _{\unicode[STIX]{x1D6FC}}(\boldsymbol{k}^{\prime \prime },t)\unicode[STIX]{x1D6F9}_{j}(-\boldsymbol{k}^{\prime \prime },t|\boldsymbol{k},\boldsymbol{k}^{\prime },s)}_{H}\nonumber\\ \displaystyle & & \displaystyle \underbrace{+\unicode[STIX]{x1D6FC}\text{i}(2\unicode[STIX]{x03C0})^{3}P_{i\unicode[STIX]{x1D6FC}}(\boldsymbol{k})\int \,\text{d}\boldsymbol{k}^{\prime \prime }\hat{\unicode[STIX]{x1D6F9}}_{j}(-\boldsymbol{k}^{\prime \prime },t|\boldsymbol{k},\boldsymbol{k}^{\prime },s)\frac{k_{\unicode[STIX]{x1D6FC}}^{\prime \prime }k_{\unicode[STIX]{x1D6FD}}^{\prime \prime }k_{\unicode[STIX]{x1D6FE}}^{\prime \prime }}{k^{\prime \prime 2}}\int _{\boldsymbol{k}^{\prime \prime }=\boldsymbol{p}+\boldsymbol{q}}d\boldsymbol{p}\hat{u} _{\unicode[STIX]{x1D6FD}}(\boldsymbol{p},t)\hat{u} _{\unicode[STIX]{x1D6FE}}(\boldsymbol{q},t)}_{I},\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$ where 
$G_{\unicode[STIX]{x1D6FC}j}^{E}(\boldsymbol{k}^{\prime \prime },t|\boldsymbol{k}^{\prime },s)$ is the Eulerian response function defined as 
$G_{\unicode[STIX]{x1D6FC}j}^{E}(\boldsymbol{k}^{\prime \prime },t|\boldsymbol{k}^{\prime },s)=\unicode[STIX]{x1D6FF}\hat{u} _{\unicode[STIX]{x1D6FC}}(\boldsymbol{k}^{\prime \prime },t)/\unicode[STIX]{x1D6FF}\hat{f}_{j}(\boldsymbol{k}^{\prime },s)$ (Kraichnan Reference Kraichnan1959). Function 
$\unicode[STIX]{x1D6F9}_{j}(\boldsymbol{k}^{\prime \prime },t|\boldsymbol{k},\boldsymbol{k}^{\prime },s)$ is the response function of 
$\hat{\unicode[STIX]{x1D713}}(\boldsymbol{k}^{\prime \prime },t|\boldsymbol{k},s)$ for 
$\hat{f}_{j}(\boldsymbol{k}^{\prime },s)$ and is defined as 
$\unicode[STIX]{x1D6F9}_{j}(\boldsymbol{k}^{\prime \prime },t|\boldsymbol{k},\boldsymbol{k}^{\prime },s)=\unicode[STIX]{x1D6FF}\hat{\unicode[STIX]{x1D713}}(\boldsymbol{k}^{\prime \prime },t|\boldsymbol{k},s)/\unicode[STIX]{x1D6FF}\hat{f}_{j}(\boldsymbol{k}^{\prime },s)$ for 
$t\geqslant s$ and zero for 
$t<s$.
The role of 
$\boldsymbol{P}$ is to enforce the dilatational component of Lagrangian velocity to be zero since 
$\hat{\boldsymbol{v}}(\boldsymbol{k},s|t)$ satisfies the solenoidal condition at 
$t=s$, but does not necessarily satisfy it at 
$t\neq s$. As shown later, this operator eliminates the irrelevant terms in the expansion and simplifies the closed equations.
Assuming that variables can be expanded in terms of 
$\unicode[STIX]{x1D6FC}$, e.g. 
$\hat{\boldsymbol{u}}(\boldsymbol{k},t)=\sum _{n=0}^{\infty }\unicode[STIX]{x1D6FC}^{n}\hat{\boldsymbol{u}}^{(n)}(\boldsymbol{k},t)$, 
$\hat{\unicode[STIX]{x1D713}}(\boldsymbol{k},t|\boldsymbol{k}^{\prime },s)=\sum _{n=0}^{\infty }\unicode[STIX]{x1D6FC}^{n}\hat{\unicode[STIX]{x1D713}}^{(n)}(\boldsymbol{k},t|\boldsymbol{k},s)$, 
$G_{ij}^{E}(\boldsymbol{k},t|\boldsymbol{k}^{\prime },s)=\sum _{n=0}^{\infty }\unicode[STIX]{x1D6FC}^{n}G_{ij}^{E(n)}(\boldsymbol{k},t|\boldsymbol{k}^{\prime },s)$ and so on, and substituting the perturbation expansions into each equation and equating powers of 
$\unicode[STIX]{x1D6FC}$, we find that, for example, 
$\hat{\unicode[STIX]{x1D713}}^{(0)}(\boldsymbol{k},t|\boldsymbol{k}^{\prime },s)=\unicode[STIX]{x1D6FF}(\boldsymbol{k}+\boldsymbol{k}^{\prime })/(2\unicode[STIX]{x03C0})^{3}$, 
$G_{ij}^{E(0)}(\boldsymbol{k},t|\boldsymbol{k}^{\prime },s)=G_{ij}^{(0)}(\boldsymbol{k},t,s)\unicode[STIX]{x1D6FF}(\boldsymbol{k}+\boldsymbol{k}^{\prime })/(2\unicode[STIX]{x03C0})^{3}$ with 
$G_{ij}^{(0)}(\boldsymbol{k},t,s)=P_{ij}(\boldsymbol{k})\exp (-\unicode[STIX]{x1D708}k^{2}(t-s))$ and so on for zeroth-order solutions. The first-order solutions are also obtained using the zeroth-order solutions, e.g. 
$\hat{u} _{i}^{(1)}(\boldsymbol{k},t)=\int _{t_{0}}^{t}\,\text{d}sG_{i\unicode[STIX]{x1D6FC}}^{(0)}(\boldsymbol{k},t,s)M_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FE}}(\boldsymbol{k})\int _{\boldsymbol{k}=\boldsymbol{p}+\boldsymbol{q}}\,\text{d}\boldsymbol{p}\hat{u} _{\unicode[STIX]{x1D6FD}}^{(0)}(\boldsymbol{p},s)\hat{u} _{\unicode[STIX]{x1D6FE}}^{(0)}(\boldsymbol{q},s)$ and so on.
In the derivation of closed equations of LRA, the following assumptions are made:
(a)
$\hat{\boldsymbol{u}}^{(0)}$, $G^{(0)}$, $\hat{\unicode[STIX]{x1D713}}^{(0)}$ and $\hat{\unicode[STIX]{x1D6F9}}^{(0)}$ are statistically independent of each other, and(b) the probability distribution of
$\hat{\boldsymbol{u}}^{(0)}$ is assumed to be Gaussian, i.e. the odd moments are zero and the even moments can be expressed by the second moments.
It is noted that the basic ideas for the derivations are different between LRA and Lagrangian DIA, although the closed equations of LRA and Lagrangian DIA have the same forms. Here is shown the systematic derivation of LRA equations.
A.1 Two-point single-time Eulerian velocity correlation
After lengthy algebra using the perturbation expansions and assumptions, we have the following equation for the two-point single-time Eulerian velocity correlation 
$Q_{ij}(\boldsymbol{k},t,t)$:
$$\begin{eqnarray}\left(\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}+2\unicode[STIX]{x1D708}k^{2}\right)Q_{ij}(\boldsymbol{k},t,t)=\unicode[STIX]{x1D6FC}^{2}D_{ij}(\boldsymbol{k},t,t)+Q_{ij}^{F}(\boldsymbol{k},t,t)+Q_{ji}^{F\ast }(\boldsymbol{k},t,t),\end{eqnarray}$$where
$$\begin{eqnarray}D_{ij}(\boldsymbol{k},t)=H_{ij}(\boldsymbol{k},t)+H_{ji}(-\boldsymbol{k},t)+O(\unicode[STIX]{x1D6FC}),\end{eqnarray}$$with
$$\begin{eqnarray}\displaystyle H_{ij}(\boldsymbol{k},t) & = & \displaystyle M_{i\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}(\boldsymbol{k})\int _{\boldsymbol{k}=\boldsymbol{p}+\boldsymbol{q}}\,\text{d}\boldsymbol{p}\nonumber\\ \displaystyle & & \displaystyle \times \,\left[2\int _{t_{0}}^{t}\,\text{d}s^{\prime }G_{j\unicode[STIX]{x1D6FE}}^{(0)}(-\boldsymbol{k},t,s^{\prime })M_{\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D716}}(-\boldsymbol{k})Q_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FF}}^{(0)}(\boldsymbol{p},t,s^{\prime })Q_{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D716}}^{(0)}(\boldsymbol{q},t,s^{\prime })\right.\nonumber\\ \displaystyle & & \displaystyle +\left.4\int _{t_{0}}^{t}\,\text{d}s^{\prime }G_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FE}}^{(0)}(\boldsymbol{p},t,s^{\prime })M_{\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D716}}(\boldsymbol{p})Q_{j\unicode[STIX]{x1D6FF}}^{(0)}(-\boldsymbol{k},t,s^{\prime })Q_{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D716}}^{(0)}(\boldsymbol{q},t,s^{\prime })\right]+O(\unicode[STIX]{x1D6FC}),\qquad\end{eqnarray}$$ where 
$Q_{ij}^{(0)}(\boldsymbol{k},t,s)=\langle \hat{u} _{i}^{(0)}(\boldsymbol{k},t)\hat{u} _{j}^{(0)}(-\boldsymbol{k},s)\rangle$, 
$Q_{ij}^{F}(\boldsymbol{k},t)=\langle \,\hat{f}_{i}(\boldsymbol{k},t)\hat{u} _{j}^{\ast }(\boldsymbol{k},t)\rangle$, the superscript 
$\ast$ denotes the complex conjugate and 
$M_{i\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}(\boldsymbol{k})=-\text{i}/2[k_{\unicode[STIX]{x1D6FC}}P_{i\unicode[STIX]{x1D6FD}}(\boldsymbol{k})+k_{\unicode[STIX]{x1D6FD}}P_{i\unicode[STIX]{x1D6FC}}(\boldsymbol{k})]$. See Leslie (Reference Leslie1973) and McComb (Reference McComb1989, Reference McComb2014) for the derivation of equation of 
$Q_{ij}(\boldsymbol{k},t,t)$.
A.2 Two-point two-time Lagrangian velocity correlation
We consider the ensemble average of (A 5). Substituting the perturbation expansions into (A 5) and taking the ensemble average, 
$\langle A\rangle$ and 
$\langle B\rangle$ are given as
$$\begin{eqnarray}\displaystyle & \langle A\rangle =Q_{ij}^{F(0)}(\boldsymbol{k},t,s)+O(\unicode[STIX]{x1D6FC}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \langle B\rangle =-\unicode[STIX]{x1D708}k^{2}Q_{ij}^{(0)}(\boldsymbol{k},t,s)+O(\unicode[STIX]{x1D6FC}). & \displaystyle\end{eqnarray}$$ After lengthy algebra using similar procedures, 
$\left\langle C\right\rangle$ is given as
$$\begin{eqnarray}\displaystyle \langle C\rangle & = & \displaystyle \unicode[STIX]{x1D6FC}^{2}4\text{i}P_{i\unicode[STIX]{x1D6FC}}(\boldsymbol{k})\frac{k_{\unicode[STIX]{x1D6FC}}k_{\unicode[STIX]{x1D6FD}}k_{\unicode[STIX]{x1D6FE}}}{k^{2}}\int _{\boldsymbol{k}=\boldsymbol{p}+\boldsymbol{q}}\,\text{d}\boldsymbol{p}\int _{t_{0}}^{t}\,\text{d}s^{\prime }G_{\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D6FF}}^{(0)}(\boldsymbol{q},t,s^{\prime })M_{\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D716}\unicode[STIX]{x1D701}}(\boldsymbol{q})Q_{j\unicode[STIX]{x1D716}}^{(0)}(\boldsymbol{k},s^{\prime },s)Q_{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D701}}^{(0)}(\boldsymbol{p},t,s^{\prime })\nonumber\\ \displaystyle & & \displaystyle +\,\unicode[STIX]{x1D6FC}^{2}2\text{i}P_{i\unicode[STIX]{x1D6FC}}(\boldsymbol{k})\frac{k_{\unicode[STIX]{x1D6FC}}k_{\unicode[STIX]{x1D6FD}}k_{\unicode[STIX]{x1D6FE}}}{k^{2}}\int _{\boldsymbol{k}=\boldsymbol{p}+\boldsymbol{q}}\,\text{d}\boldsymbol{p}\int _{t_{0}}^{s}\,\text{d}s^{\prime }G_{j\unicode[STIX]{x1D6FF}}^{(0)}(-\boldsymbol{k},s,s^{\prime })\nonumber\\ \displaystyle & & \displaystyle \times \,M_{\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D716}\unicode[STIX]{x1D701}}(-\boldsymbol{k})Q_{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D716}}^{(0)}(\boldsymbol{p},t,s^{\prime })Q_{\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D701}}^{(0)}(\boldsymbol{q},t,s^{\prime })\nonumber\\ \displaystyle & & \displaystyle -\,\unicode[STIX]{x1D6FC}^{2}2P_{i\unicode[STIX]{x1D6FC}}(\boldsymbol{k})\int _{\boldsymbol{k}=\boldsymbol{p}+\boldsymbol{q}}\,\text{d}\boldsymbol{p}\int _{s}^{t}\,\text{d}s^{\prime }\frac{p_{\unicode[STIX]{x1D6FC}}p_{\unicode[STIX]{x1D6FD}}p_{\unicode[STIX]{x1D6FE}}p_{\unicode[STIX]{x1D6FF}}}{p^{2}}Q_{\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D6FF}}^{(0)}(-\boldsymbol{q},t,s^{\prime })Q_{\unicode[STIX]{x1D6FD}j}^{(0)}(\boldsymbol{k},t,s)+O(\unicode[STIX]{x1D6FC}^{3}),\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$ where the first and second terms disappear because of 
$P_{i\unicode[STIX]{x1D6FC}}(\boldsymbol{k})k_{\unicode[STIX]{x1D6FC}}k_{\unicode[STIX]{x1D6FD}}k_{\unicode[STIX]{x1D6FE}}/k^{2}=0$. In the derivation of (A 12), symmetries of the integral with respect to 
$\boldsymbol{p}$ and 
$\boldsymbol{q}$ and 
$M_{i\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}(\boldsymbol{k})\unicode[STIX]{x1D6FF}(\boldsymbol{k})=0$ were used. Summarizing the above equations, the equation for two-point two-time Lagrangian velocity correlation is given as
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}Q_{ij}^{L}(\boldsymbol{k},t,s)=-\unicode[STIX]{x1D708}X_{ij}(\boldsymbol{k},t,s)-\unicode[STIX]{x1D6FC}^{2}I_{ij}(\boldsymbol{k},t,s)+Z_{ij}(\boldsymbol{k},t,s),\end{eqnarray}$$where
$$\begin{eqnarray}\displaystyle & \displaystyle X_{ij}(\boldsymbol{k},t,s)=k^{2}Q_{ij}^{(0)}(\boldsymbol{k},t,s)+O(\unicode[STIX]{x1D6FC}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle I_{ij}(\boldsymbol{k},t,s)=2P_{i\unicode[STIX]{x1D6FC}}(\boldsymbol{k})\int _{s}^{t}\,\text{d}s^{\prime }\int _{\boldsymbol{k}=\boldsymbol{p}+\boldsymbol{q}}\,\text{d}\boldsymbol{p}\frac{p_{\unicode[STIX]{x1D6FC}}p_{\unicode[STIX]{x1D6FD}}p_{\unicode[STIX]{x1D6FE}}p_{\unicode[STIX]{x1D6FF}}}{p^{2}}Q_{\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D6FF}}^{(0)}(-\boldsymbol{q},t,s^{\prime })Q_{\unicode[STIX]{x1D6FD}j}^{(0)}(\boldsymbol{k},t,s)+O(\unicode[STIX]{x1D6FC}), & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle Z_{ij}(\boldsymbol{k},t,s)=Q_{ij}^{F(0)}(\boldsymbol{k},t,s)+O(\unicode[STIX]{x1D6FC}). & \displaystyle\end{eqnarray}$$A.3 Lagrangian velocity response function
As in two-point two-time Lagrangian velocity correlation, we consider the ensemble-averaged equation of (A 6) for 
$\boldsymbol{k}^{\prime }=-\boldsymbol{k}$, i.e. 
$G_{ij}(\boldsymbol{k},t|-\boldsymbol{k},s)=G_{ij}(\boldsymbol{k},t,s)$. Substituting the perturbation expansions into (A 6), replacing 
$\boldsymbol{k}^{\prime }=-\boldsymbol{k}$ and taking the ensemble average, 
$\langle D\rangle$ with 
$\boldsymbol{k}^{\prime }=-\boldsymbol{k}$ is given as
$$\begin{eqnarray}\langle D\rangle =P_{ij}(\boldsymbol{k})\unicode[STIX]{x1D6FF}(t-s)+O(\unicode[STIX]{x1D6FC}).\end{eqnarray}$$ Similarly, 
$\langle E\rangle$ for 
$\boldsymbol{k}^{\prime }=-\boldsymbol{k}$ is given as
$$\begin{eqnarray}\langle F\rangle =-\unicode[STIX]{x1D708}k^{2}G_{ij}^{(0)}(\boldsymbol{k},t,s)+O(\unicode[STIX]{x1D6FC}).\end{eqnarray}$$ After lengthy algebra using similar procedures, 
$\left\langle F\right\rangle$ for 
$\boldsymbol{k}^{\prime }=-\boldsymbol{k}$ in (A 6) is given as
$$\begin{eqnarray}\displaystyle \langle F\rangle & = & \displaystyle -\unicode[STIX]{x1D6FC}^{2}2P_{i\unicode[STIX]{x1D6FC}}(\boldsymbol{k})\int _{\boldsymbol{k}=\boldsymbol{p}+\boldsymbol{q}}\,\text{d}\boldsymbol{p}\frac{p_{\unicode[STIX]{x1D6FC}}p_{\unicode[STIX]{x1D6FD}}p_{\unicode[STIX]{x1D6FE}}p_{\unicode[STIX]{x1D6FF}}}{p^{2}}\int _{s}^{t}\,\text{d}s^{\prime }Q_{\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D6FF}}^{(0)}(-\boldsymbol{q},t,s^{\prime })G_{\unicode[STIX]{x1D6FD}j}^{(0)}(\boldsymbol{k},t,s)\nonumber\\ \displaystyle & & \displaystyle +\,\unicode[STIX]{x1D6FC}^{2}4\text{i}P_{i\unicode[STIX]{x1D6FC}}(\boldsymbol{k})\int _{\boldsymbol{k}=\boldsymbol{p}+\boldsymbol{q}}\,\text{d}\boldsymbol{p}\frac{k_{\unicode[STIX]{x1D6FC}}k_{\unicode[STIX]{x1D6FD}}k_{\unicode[STIX]{x1D6FE}}}{k^{2}}\int _{s}^{t}\,\text{d}s^{\prime }G_{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FF}}^{(0)}(\boldsymbol{p},t,s^{\prime })\nonumber\\ \displaystyle & & \displaystyle \times \,M_{\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D716}\unicode[STIX]{x1D701}}(\boldsymbol{p})Q_{\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D716}}^{(0)}(\boldsymbol{q},t,s^{\prime })G_{\unicode[STIX]{x1D701}j}^{(0)}(\boldsymbol{k},s^{\prime },s)+O(\unicode[STIX]{x1D6FC}^{3}),\end{eqnarray}$$ where the second term disappears because of 
$P_{ij}(\boldsymbol{k})k_{\unicode[STIX]{x1D6FC}}k_{\unicode[STIX]{x1D6FD}}k_{\unicode[STIX]{x1D6FE}}/k^{2}=0$. Substituting perturbation expansions into 
$G$ and 
$H$ in (A 6), replacing 
$\boldsymbol{k}^{\prime }=-\boldsymbol{k}$ and taking the ensemble average, it is found that the zeroth-order terms of 
$\langle G\rangle$ and 
$\langle H\rangle$ are zeros, i.e. 
$\langle G\rangle$ and 
$\langle H\rangle$ start with 
$O(\unicode[STIX]{x1D6FC})$. Furthermore, the second-order term of 
$\langle I\rangle$ is zero, i.e. 
$\langle I\rangle$ in (A 6) starts with 
$O(\unicode[STIX]{x1D6FC}^{3})$.
Summarizing the above equations, the equation for the Lagrangian velocity response function is given as
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}G_{ij}(\boldsymbol{k},t,s)=-\unicode[STIX]{x1D708}Y_{ij}(\boldsymbol{k},t,s)-\unicode[STIX]{x1D6FC}^{2}J_{ij}(\boldsymbol{k},t,s)+W_{ij}(\boldsymbol{k},t,s),\end{eqnarray}$$where
$$\begin{eqnarray}\displaystyle & \displaystyle Y_{ij}(\boldsymbol{k},t,s)=k^{2}G_{ij}^{(0)}(\boldsymbol{k},t,s)+O(\unicode[STIX]{x1D6FC}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle J_{ij}(\boldsymbol{k},t,s)=2P_{i\unicode[STIX]{x1D6FC}}(\boldsymbol{k})\int _{\boldsymbol{k}=\boldsymbol{p}+\boldsymbol{q}}\,\text{d}\boldsymbol{p}\frac{p_{\unicode[STIX]{x1D6FC}}p_{\unicode[STIX]{x1D6FD}}p_{\unicode[STIX]{x1D6FE}}p_{\unicode[STIX]{x1D6FF}}}{p^{2}}\int _{s}^{t}\,\text{d}s^{\prime }Q_{\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D6FF}}^{(0)}(-\boldsymbol{q},t,s^{\prime })G_{\unicode[STIX]{x1D6FD}j}^{(0)}(\boldsymbol{k},t,s)+O(\unicode[STIX]{x1D6FC}), & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle W_{ij}(\boldsymbol{k},t,s)=P_{ij}(\boldsymbol{k})\unicode[STIX]{x1D6FF}(t-s)+O(\unicode[STIX]{x1D6FC}). & \displaystyle\end{eqnarray}$$A.4 Renormalized expansion
In order to close equations by the representatives, the renormalized expansion is used as in Kraichnan (Reference Kraichnan1977) and Kaneda (Reference Kaneda1981). The procedure is as follows:
(a) Begin with the primitive perturbation expansion for
$\boldsymbol{Q}(\boldsymbol{k},t,t)$, $\boldsymbol{Q}(\boldsymbol{k},t,s)$ and $\boldsymbol{G}(\boldsymbol{k},t,s)$ in terms of $\boldsymbol{Q}^{(0)}$ and $\boldsymbol{G}^{(0)}$(A 24)$$\begin{eqnarray}\displaystyle & Q_{ij}(\boldsymbol{k},t,t)=Q_{ij}^{(0)}(\boldsymbol{k},t,t)+\unicode[STIX]{x1D6FC}{\mathcal{U}}(\boldsymbol{Q}^{(0)},\boldsymbol{G}^{(0)})+\cdots \,, & \displaystyle\end{eqnarray}$$(A 25)$$\begin{eqnarray}\displaystyle & Q_{ij}(\boldsymbol{k},t,s)=Q_{ij}^{(0)}(\boldsymbol{k},t,s)+\unicode[STIX]{x1D6FC}{\mathcal{V}}(\boldsymbol{Q}^{(0)},\boldsymbol{G}^{(0)})+\cdots \,, & \displaystyle\end{eqnarray}$$(A 26)$$\begin{eqnarray}\displaystyle & G_{ij}(\boldsymbol{k},t,s)=G_{ij}^{(0)}(\boldsymbol{k},t,s)+\unicode[STIX]{x1D6FC}{\mathcal{W}}(\boldsymbol{Q}^{(0)},\boldsymbol{G}^{(0)})+\cdots \,. & \displaystyle\end{eqnarray}$$(b) Expand
$\boldsymbol{D}$, $\boldsymbol{I}$ and $\boldsymbol{J}$ in terms of $\boldsymbol{Q}^{(0)}$ and $\boldsymbol{G}^{(0)}$(A 27)$$\begin{eqnarray}\displaystyle & D_{ij}(\boldsymbol{k},t,s)=\unicode[STIX]{x1D6FC}^{2}{\mathcal{D}}(\boldsymbol{Q}^{(0)},\boldsymbol{G}^{(0)})+O(\unicode[STIX]{x1D6FC}^{3})+\cdots \,, & \displaystyle\end{eqnarray}$$(A 28)$$\begin{eqnarray}\displaystyle & I_{ij}(\boldsymbol{k},t,s)=\unicode[STIX]{x1D6FC}^{2}{\mathcal{I}}(\boldsymbol{Q}^{(0)},\boldsymbol{G}^{(0)})+O(\unicode[STIX]{x1D6FC}^{3})+\cdots \,, & \displaystyle\end{eqnarray}$$(A 29)$$\begin{eqnarray}\displaystyle & J_{ij}(\boldsymbol{k},t,s)=\unicode[STIX]{x1D6FC}^{2}{\mathcal{J}}(\boldsymbol{Q}^{(0)},\boldsymbol{G}^{(0)})+O(\unicode[STIX]{x1D6FC}^{3})+\cdots \,. & \displaystyle\end{eqnarray}$$(c) Revert these primitive expansions to obtain
$\boldsymbol{Q}^{(0)}$ and $\boldsymbol{G}^{(0)}$ as functional power expansions in $\boldsymbol{Q}$ and $\boldsymbol{G}$(A 30)$$\begin{eqnarray}\displaystyle & Q_{ij}^{(0)}(\boldsymbol{k},t,s)=Q_{ij}(\boldsymbol{k},t,s)+\unicode[STIX]{x1D6FC}{\mathcal{X}}(\boldsymbol{Q},\boldsymbol{G})+\cdots \,, & \displaystyle\end{eqnarray}$$(A 31)$$\begin{eqnarray}\displaystyle & Q_{ij}^{(0)}(\boldsymbol{k},t,s)=Q_{ij}(\boldsymbol{k},t,s)+\unicode[STIX]{x1D6FC}{\mathcal{Y}}(\boldsymbol{Q},\boldsymbol{G})+\cdots \,, & \displaystyle\end{eqnarray}$$(A 32)$$\begin{eqnarray}\displaystyle & G_{ij}^{(0)}(\boldsymbol{k},t,s)=G_{ij}(\boldsymbol{k},t,s)+\unicode[STIX]{x1D6FC}{\mathcal{Z}}(\boldsymbol{Q},\boldsymbol{G})+\cdots \,. & \displaystyle\end{eqnarray}$$(d) Substituting (A 30)–(A 32) into (A 27)–(A 29),
$\boldsymbol{D}$, $\boldsymbol{I}$ and $\boldsymbol{J}$ can be expressed in terms of representatives $\boldsymbol{Q}$ and $\boldsymbol{G}$.(e) Truncate the renormalized expansion at the lowest order and put
$\unicode[STIX]{x1D6FC}$ equal to unity.
The following equations are the final equations of LRA:
$$\begin{eqnarray}\displaystyle & \displaystyle \left(\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}+2\unicode[STIX]{x1D708}k^{2}\right)Q_{ij}(\boldsymbol{k},t)=H_{ij}(\boldsymbol{k},t)+H_{ji}^{\ast }(\boldsymbol{k},t)+Q_{ij}^{F}(\boldsymbol{k},t)+Q_{ji}^{F\ast }(\boldsymbol{k},t), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \left(\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D708}k^{2}\right)Q_{ij}(\boldsymbol{k},t,s)=I_{ij}(\boldsymbol{k},t,s)\quad t>s, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \left(\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D708}k^{2}\right)G_{ij}(\boldsymbol{k},t,s)=J_{ij}(\boldsymbol{k},t,s)\quad t>s, & \displaystyle\end{eqnarray}$$where
$$\begin{eqnarray}\displaystyle H_{ij}(\boldsymbol{k},t) & = & \displaystyle M_{i\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}(\boldsymbol{k})\int _{\boldsymbol{k}=\boldsymbol{p}+\boldsymbol{q}}d\boldsymbol{p}\nonumber\\ \displaystyle & & \displaystyle \times \,\left[2\int _{0}^{t}\,\text{d}s^{\prime }G_{j\unicode[STIX]{x1D6FE}}(-\boldsymbol{k},t,s^{\prime })M_{\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D716}}(-\boldsymbol{k})Q_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FF}}(\boldsymbol{p},t,s^{\prime })Q_{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D716}}(\boldsymbol{q},t,s^{\prime })\right.\nonumber\\ \displaystyle & & \displaystyle +\left.4\int _{0}^{t}\,\text{d}s^{\prime }G_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FE}}(\boldsymbol{p},t,s^{\prime })M_{\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D716}}(\boldsymbol{p})Q_{j\unicode[STIX]{x1D6FF}}(-\boldsymbol{k},s,s^{\prime })Q_{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D716}}(\boldsymbol{q},t,s^{\prime })\right],\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle I_{ij}(\boldsymbol{k},t,s)=-2P_{i\unicode[STIX]{x1D6FC}}(\boldsymbol{k})\int _{\boldsymbol{k}=\boldsymbol{p}+\boldsymbol{q}}\,\text{d}\boldsymbol{p}\frac{p_{\unicode[STIX]{x1D6FC}}p_{\unicode[STIX]{x1D6FD}}p_{\unicode[STIX]{x1D6FE}}p_{\unicode[STIX]{x1D6FF}}}{p^{2}}\int _{s}^{t}\,\text{d}s^{\prime }Q_{\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D6FF}}(-\boldsymbol{q},t,s^{\prime })Q_{\unicode[STIX]{x1D6FD}j}(\boldsymbol{k},t,s), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle J_{ij}(\boldsymbol{k},t,s)=-2P_{i\unicode[STIX]{x1D6FC}}(\boldsymbol{k})\int _{\boldsymbol{k}=\boldsymbol{p}+\boldsymbol{q}}\,\text{d}\boldsymbol{p}\frac{p_{\unicode[STIX]{x1D6FC}}p_{\unicode[STIX]{x1D6FD}}p_{\unicode[STIX]{x1D6FE}}p_{\unicode[STIX]{x1D6FF}}}{p^{2}}\int _{s}^{t}\,\text{d}s^{\prime }Q_{\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D6FF}}(-\boldsymbol{q},s^{\prime },s)G_{\unicode[STIX]{x1D6FD}j}(\boldsymbol{k},t,s), & \displaystyle\end{eqnarray}$$ where the contributions from the random force in (A 35) and (A 35) are neglected because they are not important for 
$t>s$. In homogeneous isotropic turbulence, the Lagrangian velocity correlation and Lagrangian velocity response function are given as
$$\begin{eqnarray}\displaystyle & Q_{ij}(\boldsymbol{k},t,s)={\textstyle \frac{1}{2}}Q(k,t,s)P_{ij}(\boldsymbol{k}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & G_{ij}(\boldsymbol{k},t,s)=G(k,t,s)P_{ij}(\boldsymbol{k}). & \displaystyle\end{eqnarray}$$The closed equations for homogeneous isotropic turbulence can be obtained by substituting these expressions into (A 33)–(A 38). For details of the derivation of the closed equation for isotropic turbulence, the reader is referred to Leslie (Reference Leslie1973) and McComb (Reference McComb1989, Reference McComb2014). For anisotropic turbulence, the reader is referred to Sagaut & Cambon (Reference Sagaut and Cambon2008).
