1. Introduction
One of the fascinating aspects of turbulence is the intricate interplay between vorticity and velocity. Indeed, the vorticity, defined as the curl of the velocity, is both advected and stretched by a turbulent velocity field. The importance of vortex stretching was recognized early on in turbulence research (see e.g. Taylor Reference Taylor1938; Betchov Reference Betchov1956; Ashurst et al. Reference Ashurst, Kerstein, Kerr and Gibson1987), and authors such as Tennekes & Lumley (Reference Tennekes and Lumley1972) and Tsinober (Reference Tsinober2009) have highlighted its important role in the nonlinear dynamics of turbulence. However, not all aspects of vortex stretching are understood. In recent investigations, the precise role of vortex stretching was investigated (Johnson & Meneveau Reference Johnson and Meneveau2016; Buaria, Bodenschatz & Pumir Reference Buaria, Bodenschatz and Pumir2020; Carbone & Bragg Reference Carbone and Bragg2020) and models were proposed to obtain a better understanding of the effects of vorticity stretching and velocity-gradient dynamics in general (Chertkov, Pumir & Shraiman Reference Chertkov, Pumir and Shraiman1999; Chevillard & Meneveau Reference Chevillard and Meneveau2006).
Different approaches can be used to obtain a better understanding of a particular feature of turbulence. One method is to attempt to disentangle, in a simulation or experiment, the influence of a particular term or structure from other features. Such an attempt is not without difficulty, because in an instantaneous flow-field it is non-trivial to recognize which flow features are caused by vortex stretching. Indeed, the whole Lagrangian history must be taken into account (Guala et al. Reference Guala, Lüthi, Liberzon, Tsinober and Kinzelbach2005). Furthermore, separating vortex stretching from other effects often requires identifying regions of flow dominated by either stretching or advection, both of which depend on thresholds and arbitrary definitions.
Another approach which can be used to identify the influence of a certain effect is to modify the physical system in order to isolate the influence of a particular feature. This is achieved by deliberately removing the feature from the dynamics, and comparing the resulting system with the original. Yet another approach is to alter the spatial dimension of the system. These methods have been used to study turbulence for decades. The resulting system does not, in general, correspond to a physical system, but shows, by comparison with the original set-up, what the influence of the additional or missing feature is.
Typical examples of such modifications are the removal of the pressure terms from the governing equations as first investigated by Burgers (Reference Burgers1950), and later by Polyakov (Reference Polyakov1995) and Boldyrev (Reference Boldyrev1999), the removal of a certain class of modes upon which the flow-field is projected to focus on certain triadic interactions (Biferale, Musacchio & Toschi Reference Biferale, Musacchio and Toschi2012; Alexakis Reference Alexakis2017; Briard, Biferale & Gomez Reference Briard, Biferale and Gomez2017; Qu, Naso & Bos Reference Qu, Naso and Bos2018), or the decimation of Fourier space to change the fractal dimension of space (Frisch, Lesieur & Sulem Reference Frisch, Lesieur and Sulem1976; Frisch et al. Reference Frisch, Pomyalov, Procaccia and Ray2012; Lanotte et al. Reference Lanotte, Benzi, Malapaka, Toschi and Biferale2015). The change of the dimension of space can also be directly investigated by reformulating turbulence in more than three dimensions (Gotoh et al. Reference Gotoh, Watanabe, Shiga, Nakano and Suzuki2007; Yamamoto et al. Reference Yamamoto, Shimizu, Inoshita, Nakano and Gotoh2012; Berera, Ho & Clark Reference Berera, Ho and Clark2020), or by considering intermediate systems such as axisymmetric turbulence, with properties of both two- and three-dimensional systems (Leprovost, Dubrulle & Chavanis Reference Leprovost, Dubrulle and Chavanis2006; Naso et al. Reference Naso, Monchaux, Chavanis and Dubrulle2010; Qu, Bos & Naso Reference Qu, Bos and Naso2017; Qin et al. Reference Qin, Faller, Dubrulle, Naso and Bos2020), or thin-layer turbulence (Celani, Musacchio & Vincenzi Reference Celani, Musacchio and Vincenzi2010; Benavides & Alexakis Reference Benavides and Alexakis2017; Favier, Guervilly & Knobloch Reference Favier, Guervilly and Knobloch2019).
The present work follows the approach of altering the Navier–Stokes equations. The modification is drastic because it involves the removal of vortex stretching from the system, which in turn changes the nonlinearity of the Navier–Stokes equations. One way to do this is to consider two-dimensional turbulence, because in a two-dimensional velocity field, the vorticity is perpendicular to the velocity (and its gradient) so that the vortex stretching term drops out of the vorticity equation. However, this method also changes the dimension of the system. In this investigation vortex stretching is removed from the three-dimensional Navier–Stokes equations without changing the space dimension, and the resulting statistical properties of the system are analysed. This approach can also be seen as extending the case of pure advection of vorticity, as in two-dimensional turbulence, to a system with a higher dimension.
The approach by which we investigate the resulting system uses closure theory, which provides a method of analysing non-physical systems. Furthermore, for an exploratory investigation such as this, the method allows us to explore, at low-computational cost, assumptions of asymptotic scaling. The assumptions underlying closure theory do not violate the detailed conservation properties of invariants of the governing equations and do, in general, correctly capture asymptotic scaling of second- and third-order velocity correlations (Sagaut & Cambon Reference Sagaut and Cambon2008). Higher-order moments, reflecting the intermittency properties of the flow can in principle be addressed by closure (Chen et al. Reference Chen, Herring, Kerr and Kraichnan1989; Bos & Rubinstein Reference Bos and Rubinstein2013), but are not always reproduced correctly. The instantaneous structure of the flow cannot be reproduced by a purely statistical approach. It is therefore recommended that the present results be validated using direct numerical simulations, but this will be left for future work.
In the following section (§ 2) we present the theoretical framework and modified Navier–Stokes system. Then, in § 3 we derive a closed expression for the evolution equation of the kinetic energy using an eddy-damped quasi-normal Markovian (EDQNM) approach. The equations are integrated numerically and the results are presented in § 4. Section 5 contains the conclusions.
2. Navier–Stokes equations without vortex stretching
We consider three-dimensional, incompressible, statistically homogeneous turbulence, maintained by an external force term 
$\boldsymbol {f}(\boldsymbol {x},t)$. The velocity 
$\boldsymbol {u}(\boldsymbol {x},t)$ of this unmodified flow is governed by the Navier–Stokes equations
\begin{equation} \frac{\partial \boldsymbol{u}}{\partial t}-\nu \Delta \boldsymbol{u}={-}\boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}-\boldsymbol{\nabla} \mathcal{P} +\boldsymbol{f}, \end{equation}
with 
$\mathcal {P}(\boldsymbol {x},t)$ the pressure, 
$\nu$ the kinematic viscosity and 
$\boldsymbol {\nabla }\boldsymbol {\cdot } \boldsymbol {u}=0$. The time dependence of the velocity, force and pressure is omitted here.
In order to remove the vortex stretching from the dynamics, we consider the curl of (2.1), yielding the vorticity equation
\begin{equation} \frac{\partial \boldsymbol{\omega}}{\partial t}+\underbrace{\boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{\omega}}_{{Advection}}-\nu \Delta \boldsymbol{\omega}= \underbrace{\boldsymbol{\omega}\boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}}_{{Stretching}}+ \boldsymbol{\nabla} \times \boldsymbol{f}. \end{equation}If the vorticity is only advected and not stretched, we remove the stretching term, leading to
\begin{equation} \frac{\partial \boldsymbol{\omega}}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{\omega}-\nu \Delta \boldsymbol{\omega}=\boldsymbol{\nabla} \times \boldsymbol{f}. \end{equation}This equation is similar to the two-dimensional Navier–Stokes equations, but with the important difference that vorticity is now a three-component vector. The corresponding Navier–Stokes equation is as follows
\begin{equation} \frac{\partial \boldsymbol{u}}{\partial t}-\nu \Delta \boldsymbol{u}={-}\boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}-\boldsymbol{\nabla} \mathcal{P} - \boldsymbol{\phi}+\boldsymbol{f}, \end{equation}
with 
$\boldsymbol {\phi }$ a force, or damping, applied to the velocity field defined such that
\begin{equation} \boldsymbol{\nabla} \times \boldsymbol{\phi}=\boldsymbol{\omega}\boldsymbol{\cdot} \boldsymbol{\nabla}\boldsymbol{u}. \end{equation}Such an artificial forcing term, applied to all scales is analogous to the helical forcing used in Plunian et al. (Reference Plunian, Teimurazov, Stepanov and Verma2020). It is important to note that the present investigation does not address the question of whether it is vortex stretching or strain self-amplification that is more significant in the process of energy transfer (Johnson & Meneveau Reference Johnson and Meneveau2016; Carbone & Bragg Reference Carbone and Bragg2020), because removing the vortex-stretching term from the vorticity equation will simultaneously suppress the strain self-amplification from the dynamics. It is, therefore, implied in the following discussion that references to vortex stretching also include strain self-amplification.
From (2.3) for a periodic or statistically homogeneous system
\begin{equation} \frac{\textrm{d}Z}{\textrm{d}t}=\beta_{in}-\beta, \end{equation}
with the enstrophy 
$Z$ given by
\begin{equation} Z=\tfrac{1}{2}\langle \|\boldsymbol{\omega}\|^2\rangle, \end{equation}the enstrophy injection by external forcing
\begin{equation} \beta_{in}=\langle \boldsymbol{\omega}\boldsymbol{\cdot} \boldsymbol{\nabla} \times \boldsymbol{f}\rangle, \end{equation}and the enstrophy dissipation
\begin{equation} \beta={-}\nu\langle \boldsymbol{\omega}\boldsymbol{\cdot} \Delta \boldsymbol{\omega}\rangle , \end{equation}
so that the enstrophy of the unforced inviscid system (
$\nu =0, \boldsymbol {f}=0$) obeys
\begin{equation} \frac{\textrm{d}Z}{\textrm{d}t}=0, \end{equation}and is thus conserved by the nonlinear interactions of the system.
The kinetic energy balance is
\begin{equation} \frac{\textrm{d}K}{\textrm{d}t}=\epsilon_{in}-\epsilon-\varPsi, \end{equation}
where kinetic energy 
$K$, energy input and viscous dissipation are defined, respectively, by
\begin{gather} K=\tfrac{1}{2}\langle \|\boldsymbol{u}\|^2\rangle , \end{gather}
\begin{gather}\epsilon_{in}=\langle \boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{f}\rangle , \end{gather}
\begin{gather}\epsilon={-}\nu\langle \boldsymbol{u}\boldsymbol{\cdot} \Delta \boldsymbol{u}\rangle. \end{gather}The energy input or destruction due to the absence of the vortex stretching term is
\begin{equation} \varPsi=\langle \boldsymbol{\phi} \boldsymbol{\cdot} \boldsymbol{u} \rangle. \end{equation}In the unforced, inviscid case, the energy balance reads
\begin{equation} \frac{\textrm{d}K}{\textrm{d}t}={-}\varPsi, \end{equation}
so that the inviscid system ((2.4) with 
$\nu =0$) does not necessarily conserve energy.
From these considerations it follows that 
$Z$ (or other moments of the vorticity) are conserved by the nonlinearity of the system. We have not identified other invariants. Here we consider only the mirror-symmetric case. Whether an invariant, such as the volume averaged helicity, is conserved for the non-mirror-symmetry case will be left for future research.
3. Fourier representation and closure of the system
In Fourier space the Biot–Savart operator becomes an algebraic operation, which allows us to rewrite the Navier–Stokes equations without vortex stretching in a convenient form. The explicit form of the evolution of the Fourier modes is derived in Appendix A and can be written in the form
\begin{equation} \frac{\partial u_i(\boldsymbol{k})}{\partial t}={-}\textrm{i} \int \delta(\boldsymbol{k}-\boldsymbol{p}-\boldsymbol{q}) \varGamma_{ijm}(\boldsymbol{k},\boldsymbol{p},\boldsymbol{q})u_j(\boldsymbol{p})u_m(\boldsymbol{q}) \, \textrm{d}\boldsymbol{p}\,\textrm{d}\boldsymbol{q}-\nu k^2 u_i(\boldsymbol{k})+\,f_i(\boldsymbol{k}), \end{equation}where
\begin{equation} \varGamma_{ijm}(\boldsymbol{k},\boldsymbol{p},\boldsymbol{q})=\left[ \left(\lambda \frac{\boldsymbol{p}\boldsymbol{\cdot} \boldsymbol{k}}{k^2}+ \frac{\boldsymbol{q}\boldsymbol{\cdot} \boldsymbol{k}}{k^2}\right)k_jP_{im}(\boldsymbol{k})\right] \end{equation}
and 
$P_{ij}(\boldsymbol {k})=\delta _{ij}-k_ik_j k^{-2}$. Note that for 
$\lambda =1$ we have 
$\varGamma _{ijm}({\boldsymbol {k},\boldsymbol {p},\boldsymbol {q}})=k_jP_{im}(\boldsymbol {k})$ and we retrieve the unmodified Navier–Stokes equations. The case 
$\lambda =0$ corresponds to the dynamics without vortex stretching. Here and in the following we will distinguish the Fourier coefficients 
$\boldsymbol {u}(\boldsymbol {k})$ from the associated velocity field 
$\boldsymbol {u}$ by their explicit dependence on the wave vector.
The energy spectrum, defined as the spherically averaged energy density in Fourier space, is then governed by the Lin equation
\begin{equation} \frac{\partial E(k)}{\partial t}=T(k)-2\nu k^2 E(k)+P(k), \end{equation}
where 
$P(k)$ represents the energy input,
\begin{equation} \int P(k) \,\textrm{d}k=\langle \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{f} \rangle\equiv \epsilon_{{in}}. \end{equation}The nonlinear transfer term is given by
\begin{equation} T(k)={-} 4\textrm{i}{\rm \pi} k^2 \int \delta(\boldsymbol{k}-\boldsymbol{p}-\boldsymbol{q}) \varGamma_{ijm}({\boldsymbol{k},\boldsymbol{p},\boldsymbol{q}})\langle u_i( -\boldsymbol{k})u_j(\boldsymbol{p})u_m(\boldsymbol{q})\rangle \,\textrm{d}\boldsymbol{p} \,\textrm{d}\boldsymbol{q}, \end{equation}and its integral is
\begin{equation} \int T(k) \,\textrm{d}k={-}\varPsi, \end{equation}which is strictly zero in non-modified turbulence. The dissipation spectrum is related to the viscous dissipation by
\begin{equation} \int 2\nu k^2 E(k) \,\textrm{d}k=\epsilon. \end{equation}The nonlinear transfer contains a triple velocity correlation. For this correlation we derive a closed expression within the framework of the EDQNM theory (Orszag Reference Orszag1970) and the details of the derivation can be found in Appendix B. For the case without vortex stretching we have
\begin{align} T(k)&=\frac{1}{2}\int_\varDelta \frac{\textrm{d}p}{p}\frac{\textrm{d}q}{q}\varTheta_{k,p,q}\left\{\left[f^{\lambda}_{(1)} X_{(1)}+f^{\lambda}_{(2)} X_{(2)}\right]k^3 E(p)E(q) \right.\nonumber\\ &\quad \left. -\left[f^{\lambda}_{(3)} X_{(3)}+ f^{\lambda}_{(4)} X_{(4)}\right]p^3 E(k)E(q) -\left[f^{\lambda}_{(5)} X_{(5)}+ f^{\lambda}_{(6)} X_{(6)}\right]q^3 E(p)E(k)\right\}, \end{align}
where the 
$\varDelta$ denotes the integration domain in the 
$p$–
$q$ plane, where 
$k,p,q$ can form the sides of a triangle. The 
$X$ terms are
\begin{align} X_{(1)}&=(1-z^2)(1+y^2),\quad X_{(2)}={-}xyz-y^2z^2 \end{align}
\begin{align}X_{(3)}&=xy(1-z^2),\quad\quad\quad\ \,X_{(4)}=z({-}y^2-xyz) \end{align}
\begin{align}X_{(5)}&=y({-}x^2-xyz),\quad\quad\ X_{(6)}=(y+zx)(1+y^2) \end{align}
and the 
$f^{\lambda }$ terms read for 
$\lambda =0$,
\begin{align} f^{\lambda}_{(1)}&=y^2\left(\frac{q}{k}\right)^2,\quad f^{\lambda}_{(2)}=yz\frac{pq}{k^2} \end{align}
\begin{align}f^{\lambda}_{(3)}&=xy\frac{q^2}{pk},\quad\quad\ \,f^{\lambda}_{(4)}=yz\frac{q}{p} \end{align}
\begin{align}f^{\lambda}_{(5)}&=xy\frac{p}{k},\quad\quad\ \ \ f^{\lambda}_{(6)}=y^2. \end{align}The triad interaction time is defined as
\begin{equation} \varTheta(k,p,q)=\frac{1-\exp[-(\theta_k+\theta_p+\theta_q)t]}{\theta_k+\theta_p+\theta_q} \end{equation}with
\begin{equation} \theta_k=\alpha\sqrt{\int_0^k s^2 E(s)\,\textrm{d}s} +\nu k^2, \end{equation}
and 
$\alpha =0.5$. For the classical EDQNM closure we have, using 
$\lambda =1$ in the derivation, that all 
$f^{\lambda }$ terms are equal to unity, yielding
\begin{align} T(k)&=\frac{1}{2}\int_\varDelta \frac{\textrm{d}p}{p}\frac{\textrm{d}q}{q}\varTheta_{k,p,q} \left\{\left[X_{(1)}+X_{(2)}\right]k^3 E(p)E(q) \right. \nonumber\\ &\quad\left. -\left[ X_{(3)}+ X_{(4)}\right]p^3 E(k)E(q) -\left[X_{(5)}+ X_{(6)}\right]q^3 E(p)E(k)\right\}, \end{align}which can be symmetrized to obtain the familiar expression found in the literature.
4. Results
In this section we integrate the closure equations. We first highlight the differences with respect to unmodified three-dimensional turbulence. Subsequently we investigate, by dimensional analysis and variation of the Reynolds number, the inertial range of the system. Finally, we consider the inviscid system to show the equilibrium properties of the system.
4.1. Numerical set-up and parameters
To integrate the equations we use the same in-house EDQNM code that has been used over the last few decades in our laboratory (Sagaut & Cambon Reference Sagaut and Cambon2008), containing a routine written by Leith (Reference Leith1971). For all simulations a logarithmic discretization was used with a minimum resolution of 20 points per decade. Simulations were performed on a domain 
$k\in [0.05 k_f,5k_\eta ]$, where 
$k_f$ is the forcing frequency (
$k_f=1$) and 
$k_\eta =\nu ^{-3/4}\epsilon _{in}^{1/4}$. The forcing term in (3.3) is defined by (Briard et al. Reference Briard, Biferale and Gomez2017),
\begin{equation} P(k)=A\exp\left({-}100 [\ln(k/k_f)]^2\right) \end{equation}
with 
$A$ determined such that the integral of 
$P(k)$ is unity. This ensures an energy input which is fixed at 
$\epsilon _{{in}}=1$ and an enstrophy input close to unity (
$\beta _{{in}}=1.02$). With the exception of the inviscid relaxation simulations in § 4.4 all results are reported when a statistically steady state is reached.
4.2. Comparison with classical turbulence and scaling ranges
In order to highlight the differences between unmodified high-Reynolds-number turbulence and the system without vortex stretching, we first present the results of the integration of the two systems, using the same parameters, defined in the foregoing paragraph, with viscosity 
$\nu =1\times 10^{-4}$. This corresponds, for the unmodified turbulence, to a Taylor-scale Reynolds number 
$R_\lambda \approx 10^3$.
The resulting spectra are shown in figure 1. The Navier–Stokes system yields a spectrum with for 
$k>k_f$ a Kolmogorov 
$k^{-5/3}$ scaling and for 
$k<k_f$ a 
$k^2$ dependence, reminiscent of energy-equipartition. For large 
$k$, a dissipation range is observed, where the energy spectrum falls off more rapidly than a power law.
Figure 1. Steady-state kinetic energy spectra of forced Navier–Stokes turbulence compared with the energy spectra of the system without vortex stretching for the same parameters. VS, vortex stretching.
The system without vortex stretching shows a peak, representing the forcing scale, similar to that of the Navier–Stokes system. However, for larger and smaller wave numbers, the spectrum is steeper or shallower, respectively. An inertial range is observed with a wave-number dependence close to 
$k^{-3}$. A dissipation range is observed which starts at a wave number about 10 times smaller than for the reference case. Furthermore, for small 
$k$, an approximately flat spectrum is observed. This range is associated with the statistical equilibrium properties of the new system and we return to this in § 4.4.
In figure 2, transfer spectra are presented. Multiplication by 
$k$ allows to assess in this semi-logarithmic representation the conservation of energy, by comparing the positive and negative areas (or lobes) delimited by the spectra (figure 2a). In the Navier–Stokes case the energy is transferred towards the small scales. This is illustrated by the negative dip around the forcing frequency, where the nonlinear interaction absorbs the energy, and the positive lobe near the dissipation scale where it is expelled. In between these two scales the energy is conserved. The energy transfer of the turbulent flow without vortex stretching shows that virtually no energy is transferred, and all energy is locally destroyed by nonlinear interaction. Indeed, no positive lobe in the transfer is observed, illustrating that the energy is absorbed and destroyed by the nonlinearity. The amount of energy destroyed by nonlinear effects corresponds to the term 
$\varPsi$ in the energy balance equation (2.11).
Figure 2. Spectra of (a) energy transfer and (b) enstrophy transfer for both systems. The energy transfer is normalized by the energy injection rate 
$\epsilon _{in}$, whereas the enstrophy spectra are normalized by the enstrophy dissipation rate 
$\beta$. VS, vortex stretching.
By multiplying the transfer spectrum by 
$k^2$, the enstrophy-transfer spectrum is obtained. It is observed (figure 2b) that the system without vortex stretching conserves enstrophy, transferring it from the injection scale to the enstrophy dissipation scale. The Navier–Stokes enstrophy balance shows that enstrophy is strongly enhanced throughout the cascade. Normalizing the transfer spectra by the enstrophy dissipation 
$\beta$ shows that the amount of enstrophy at the viscous end of the cascade is so much larger than its injected value, that the latter is negligible. Indeed, considering the enstrophy transfer near the injection scales would appear to indicate that no enstrophy is injected. However, both systems are subject to equal amounts of injected enstrophy, and owing to the very strong production of enstrophy, caused by vortex stretching, the normalized spectra mask the injected enstrophy in this representation.
We quantify this by the ratios of injected to dissipated energy and injected to dissipated enstrophy. We obtain, for the case of Navier–Stokes turbulence,
\begin{equation} \frac{\epsilon}{\epsilon_{in}}=1,\quad \frac{\beta}{\beta_{in}}=1.3\times 10^4 \end{equation}and without vortex stretching,
\begin{equation} \frac{\epsilon}{\epsilon_{in}}=2.5\times 10^{{-}3},\quad \frac{\beta}{\beta_{in}}=1. \end{equation}Comparing these values illustrates the enormous amount of enstrophy which is generated by vortex stretching at high-Reynolds-number turbulence. It also shows the importance of vortex stretching to the process of energy conservation.
In order to complete the picture, we show in figure 3 the fluxes associated with the conserved quantities. In figure 3(a,b) we show the energy and enstrophy flux towards small scales, defined, respectively, as
\begin{equation} \varPi_E(k)={-}\int_0^k T(k) \,\textrm{d}k,\quad \varPi_Z(k)={-}\int_0^k k^2 T(k) \,\textrm{d}k. \end{equation}
A clear inertial range is observed for both quantities, where the flux is approximately constant. In both cases, fluxes are in the direction of the small scales. Indeed, a steady state is observed where scales 
$k < k_f$ are in statistical equilibrium. This equipartition state is associated with zero net transfer. Extending the wave-number domain to smaller 
$k$ extends this equipartition state, where 
$E(k)\sim k^2$ for Navier–Stokes turbulence and 
$E(k)\sim k^0$ for the system without vortex stretching. It is, therefore, not necessary to add large-scale friction to the system in order to achieve steady state, unlike the case of two-dimensional turbulence, where energy accumulates in the forced system in the absence of large-scale damping terms. Further analysis of the equipartition range is postponed to § 4.4, where the truncated inviscid system is considered.
Figure 3. Fluxes of (a) energy in Navier–Stokes turbulence and (b) enstrophy in the system without vortex stretching (VS). The fluxes are normalized by the energy injection rate 
$\epsilon _{in}$ and enstrophy injection rate 
$\beta _{in}$, respectively.
4.3. Dimensional analysis and scaling
The foregoing analysis shows that, without vortex stretching, enstrophy is conserved by the nonlinear interactions and is transferred to small scales, where it is dissipated. No inverse (or direct) cascade of energy is observed. Kolmogorov-type arguments, assuming scale locality, will lead to a scaling that is dependent on the enstrophy flux and the local length scale (or wave number). In a steady state the enstrophy flux is equal to the enstrophy dissipation 
$\beta$, so that we obtain from dimensional arguments that
\begin{equation} E(k)\sim \beta^{2/3} k^{{-}3}. \end{equation}The equivalent of the viscous Kolmogorov scale will now become
\begin{equation} \zeta \equiv \frac{\nu^{1/2}}{\beta^{1/6}}. \end{equation}
The energy spectra should then collapse in the high-wave-number range, for large Reynolds numbers using the length scale 
$\zeta$ and the viscosity
\begin{equation} E(k)=\nu^{3/2}\beta^{1/6}f(k\zeta), \end{equation}
where 
$f$ is a unique function. This is assessed in figure 4 where we show plots of the energy spectrum associated with turbulence in the absence of vortex stretching, for 
$\nu =0.01,0.001,0.0001$. In figure 4(a) we show the non-normalized spectra, which coincide at the large scales. Normalizing using the above scaling arguments allows data collapse for all three cases in the dissipation and inertial ranges (figure 4b).
Figure 4. Reynolds number scaling of the energy spectra of turbulence without vortex stretching. (a) Non-normalized spectra for three distinct values of the viscosity. (b) The same spectra normalized by enstrophy dissipation and viscosity.
These scaling arguments can also be used to explain why in figure 1 the viscous cut-off of both systems is different by an order of magnitude. The Kolmogorov scale in turbulence is of order 
$\eta =\nu ^{3/4}\epsilon _{{in}}^{-1/4}$ and the Kolmogorov-like scale in the modified system is given by expression (4.6). The ratio is then
\begin{equation} \frac{\eta}{\zeta}=\frac{\beta_{{in}}^{1/6}}{\epsilon_{{in}}^{1/4}}\nu^{1/4}. \end{equation}
As both 
$\epsilon _{{in}}$ and 
$\beta _{{in}}$ are order unity, and 
$\nu =10^{-4}$, this ratio is around 10 for the spectra shown in figure 1.
4.4. Inviscid equilibrium
Previously it was observed that for scales larger than the forcing scale, i.e. for wave numbers 
$k<k_f$, the energy spectrum is flat in turbulence without vortex stretching. In the current section, we show that this scaling is associated with the inviscid equilibrium state of the system.
The inviscid equilibrium properties of turbulence have received much interest motivated by the establishment of a direct connection between thermodynamic and fluid mechanics properties of such flows. Early investigations showed that a Galerkin truncation of the Navier–Stokes system, in the absence of viscosity, allows an equilibrium solution where all Fourier modes contain, statistically, the same amount of energy (Lee Reference Lee1952; Kraichnan Reference Kraichnan1973). The resulting system shows then an energy spectrum proportional to 
$k^2$. Kraichnan (Reference Kraichnan1967) extended these ideas to two-dimensional turbulence. In the present case, where only one particular invariant is present in the system, it is plausible that the equilibrium distribution corresponds to an equidistribution of enstrophy between modes. As the enstrophy spectrum is related to the energy spectrum by
\begin{equation} E_Z(k)\sim k^2 E(k), \end{equation}we can expect an equilibrium spectrum,
\begin{equation} E(k)\sim k^0. \end{equation}We check this by integrating the inviscid system, starting from a concentrated energy (and enstrophy) distribution,
\begin{equation} E(k,0)=H(2-k) \end{equation}
with 
$H$ the Heaviside function. The domain is 
$k\in [1,100]$. In figure 5(a) we show the short-time evolution of the system. At very short times we observe the staircase scaling recently discussed in Fang, Wu & Bos (Reference Fang, Wu and Bos2020). Then, at intermediate times, as shown in figure 5(b) an enstrophy cascade coexists with a thermalized part as also observed in the three-dimensional case for the energy cascade (Cichowlas et al. Reference Cichowlas, Bonaïti, Debbasch and Brachet2005; Bos & Bertoglio Reference Bos and Bertoglio2006). The relaxation towards an 
$E(k)\sim k^0$ spectrum shows the equipartition of enstrophy. It also explains the wave-number dependence of the energy spectra for scales larger than the forcing scale observed in figure 1. Indeed, for such large scales in forced three-dimensional turbulence the modes are shown to be in thermal equilibrium, showing a 
$k^2$ equipartition energy spectrum, as observed in figure 1 (see, for instance, Alexakis & Brachet (Reference Alexakis and Brachet2019) for a discussion). Transposing this to the enstrophy-conserving dynamics observed in the present investigation suggests the observed 
$k^0$ scaling observed in figure 1 for turbulence without vortex stretching.
Figure 5. Inviscid relaxation to an equilibrium state where enstrophy is equipartitioned over the different modes. (a) Short time behaviour, showing the kinetic energy spectra at 
$t=0$ and 
$t=0.04;0.08;0.16;0.32;0.64;1.28;2.56;5.12$. (b) Long time evolution, showing the spectra at 
$t=5.12;10;20;40;80;160;320;640;1280$ and 
$t=10^4$.
5. Conclusion
Three-dimensional turbulence without vortex stretching is different from two-dimensional turbulence. Both systems conserve enstrophy and cascade that quantity to the small scales. However, whereas in two-dimensional turbulence energy is transferred towards large scales, displaying thereby a double cascade, the modified three-dimensional system does not conserve energy and a simple cascade is observed, associated with a 
$k^{-3}$ inertial range.
This absence of vortex stretching also alters the absolute equilibrium states of the truncated inviscid system. In the two-dimensional case the statistical equilibrium is a function of both invariants, whereas in three-dimensional turbulence without vortex stretching the equilibrium distribution is a simple equipartition of enstrophy, corresponding to a flat 
$k^0$ energy spectrum. This behaviour is also observed in the forced system for scales larger than the forcing scale.
What we can therefore safely state is that vortex stretching is inseparable from the energy cascade mechanism. Indeed, in its absence, energy is not conserved and the dynamics of the flow is radically changed. In classical turbulence at high Reynolds numbers, vortex stretching amplifies enstrophy by several orders of magnitude between the injection scale and the dissipation scale. In the absence of vortex stretching, the enstrophy becomes scale independent in the inertial range. However, in the same range for the system without stretching, the energy is destroyed and only a very small fraction survives the cascade towards the dissipation scale. We repeat here that the present investigation does not address the dynamical differences between vortex stretching and strain self-amplification (Johnson & Meneveau Reference Johnson and Meneveau2016; Carbone & Bragg Reference Carbone and Bragg2020), because both effects are suppressed by removing the vortex-stretching term from the vorticity equation. It does not seem easy to remove only one of the two effects from the Navier–Stokes equations without altering the other.
The present investigation opens up several possibilities for future work. In particular, direct numerical simulations will provide enough data to assess the fine properties of the flow, such as intermittency, and allow researchers to investigate the physical space structure of this new type of turbulence.
Recently, turbulence was investigated using local surgery of the velocity field, where strongly vortical regions were locally damped (Buzzicotti, Biferale & Toschi Reference Buzzicotti, Biferale and Toschi2020). A variant of the present approach, where vortex stretching is suppressed locally, might be of interest. A parametric study, varying 
$\lambda$ in (3.2) between zero and unity would allow an assessment of how exactly the statistics depend on the strength of the vortex stretching. Indeed, the procedure of fractal decimation (Lanotte et al. Reference Lanotte, Benzi, Malapaka, Toschi and Biferale2015) shows that statistical properties such as intermittency can be extremely sensitive to the variation of a control parameter around a critical value or dimension.
A final possibility is the mathematical investigation of turbulence without vortex stretching. Indeed, because the steep spectral energy distribution in the inertial range suggests that the flow is statistically smooth, considerations about the existence and uniqueness of solutions of the present system might suggest new methods of assessing the mathematical properties of Navier–Stokes turbulence.
Acknowledgement
I thank T. Wu for useful comments.
Declaration of interest
The author reports no conflict of interest.
Appendix A. Fourier representation of the Navier–Stokes equation without vortex stretching
The vorticity equation without vortex stretching reads
\begin{equation} \frac{\partial \omega_i(\boldsymbol{k})}{\partial t}={-}\textrm{i} k_j \int_\delta u_j(\boldsymbol{p})\omega_i(\boldsymbol{q}), \end{equation}where we use the short-hand notation
\begin{equation} \int_\delta\equiv \int \textrm{d}\boldsymbol{p} \,\textrm{d}\boldsymbol{q} \delta(\boldsymbol{k}-\boldsymbol{p}-\boldsymbol{q}). \end{equation}Here, and in the following, we omit time arguments and Fourier modes are indicated by their dependence on the wave vector. Forcing and viscous terms can be added afterwards.
As for a solenoidal field we have 
$\boldsymbol {\nabla }\times \boldsymbol {\nabla } \times \boldsymbol {u}=-\Delta \boldsymbol {u}$, and in Fourier space the Laplacian becomes an algebraic operator, the vorticity equation can be easily uncurled, yielding
\begin{equation} \frac{\partial u_i(\boldsymbol{k})}{\partial t}=\frac{\textrm{i}}{k^2}\epsilon_{iab}\epsilon_{bcd} \int_\delta k_a k_j q_c u_j(\boldsymbol{p})u_d(\boldsymbol{q}). \end{equation}Developing the permutation tensor gives
\begin{equation} \frac{\partial u_i(\boldsymbol{k})}{\partial t} =\frac{\textrm{i}}{k^2}\int_\delta ( k_d k_j q_i u_j(\boldsymbol{p})u_d(\boldsymbol{q}) - k_aq_a k_j u_j(\boldsymbol{p})u_i(\boldsymbol{q})). \end{equation}The first term in parentheses can be symmetrized
\begin{equation} \frac{\partial u_i(\boldsymbol{k})}{\partial t} =\frac{\textrm{i}}{k^2}\int_\delta ( k_d k_j q_i u_j(\boldsymbol{p})u_d(\boldsymbol{q})/2+k_d k_j p_i u_j(\boldsymbol{p})u_d(\boldsymbol{q})/2 - k_aq_a k_j u_j(\boldsymbol{p})u_i(\boldsymbol{q})), \end{equation}
and using 
$\boldsymbol {p}+\boldsymbol {q}=\boldsymbol {k}$ gives
\begin{equation} \frac{\partial u_i(\boldsymbol{k})}{\partial t} =\frac{\textrm{i}}{k^2}\int_\delta ( k_d k_j k_i u_j(\boldsymbol{p})u_d(\boldsymbol{q})/2 - k_aq_a k_j u_j(\boldsymbol{p})u_i(\boldsymbol{q})). \end{equation}
Removing the potential part by multiplying both sides with 
$P_{in}(\boldsymbol {k})$ and relabeling gives
\begin{equation} \frac{\partial u_i(\boldsymbol{k})}{\partial t} ={-}\textrm{i}k_jP_{im}(\boldsymbol{k})\int_\delta\frac{(k_aq_a)}{k^2} u_j(\boldsymbol{p})u_m(\boldsymbol{q})). \end{equation}Comparison with the Navier–Stokes equations,
\begin{equation} \frac{\partial u_i(\boldsymbol{k})}{\partial t} ={-}\textrm{i}k_jP_{im}(\boldsymbol{k})\int_\delta u_j(\boldsymbol{p})u_m(\boldsymbol{q})), \end{equation}allows us to write the general form
\begin{equation} \frac{\partial u_i(\boldsymbol{k})}{\partial t} ={-}\textrm{i}\int_\delta \left(\lambda\frac{(k_ap_a)}{k^2}+\frac{(k_aq_a)}{k^2}\right)k_jP_{im}(\boldsymbol{k}) u_j(\boldsymbol{p})u_m(\boldsymbol{q})), \end{equation}
which gives the case of turbulence without vortex stretching for 
$\lambda =0$, and which, for 
$\lambda =1$, reduces to the Navier–Stokes equations, because the term in parentheses yields unity.
Appendix B. EDQNM closure of the nonlinear transfer
We follow the procedure outlined in Bos & Bertoglio (Reference Bos and Bertoglio2013), leading to equations of the EDQNM family. Alternative procedures (Orszag Reference Orszag1970; Sagaut & Cambon Reference Sagaut and Cambon2008) should yield the same closure.
We start with
\begin{equation} \frac{\partial u_i(\boldsymbol{k})}{\partial t}= \int_\delta\varGamma_{ijm}(\boldsymbol{k},-\boldsymbol{p},-\boldsymbol{q})\left[u_j(\boldsymbol{p})u_m(\boldsymbol{q})\right]-\nu k^2 u_i(\boldsymbol{k}). \end{equation}The velocity correlations obey
\begin{equation} \left(\frac{1}{2}\frac{\partial }{\partial t}+\nu k^2\right) \left\langle u_i(\boldsymbol{k})u_i(-\boldsymbol{k})\right\rangle= \int_\delta\varGamma_{ijm}(\boldsymbol{k},\boldsymbol{p},\boldsymbol{q})\langle u_j(\boldsymbol{p})u_m(\boldsymbol{q})u_i(-\boldsymbol{k})\rangle, \end{equation}
where the triple correlations need to be determined. We can formally invert equation (B1) to obtain expressions for evolution of the three modes constituting the triple correlations, 
$u=u^{(0)}+u^{(1)}$, with 
$u^{(0)}$ the Gaussian or independent velocity estimate, and 
$u^{(1)}$ the perturbation by nonlinear direct triad interaction within the mode 
$\boldsymbol {k}=\boldsymbol {p}+\boldsymbol {q}$,
\begin{align} \left.\begin{gathered} u^{(1)}_i(-\boldsymbol{k})= \int_0^t \,\textrm{d}s G'(k)\left(\varGamma_{iab}(-\boldsymbol{k},\boldsymbol{p},\boldsymbol{q}) \left[u'_a(-\boldsymbol{p})u'_b(-\boldsymbol{q})\right] +\varGamma_{iab}(-\boldsymbol{k}, \boldsymbol{q},\boldsymbol{p})\left[u'_a(-\boldsymbol{q})u'_b(-\boldsymbol{p})\right]\right),\\ u^{(1)}_j(\boldsymbol{p})=\int_0^t \,\textrm{d}s G'(p)\left(\varGamma_{jab}(\boldsymbol{p},-\boldsymbol{k},\boldsymbol{q})\left[u'_a(\boldsymbol{k})u'_b (-\boldsymbol{q})\right] +\varGamma_{jab}(\boldsymbol{p},\boldsymbol{q},-\boldsymbol{k})\left[u'_a(-\boldsymbol{q})u'_b(\boldsymbol{k})\right]\right),\\ u^{(1)}_m(\boldsymbol{q})= \int_0^t \,\textrm{d}s G'(q)\left(\varGamma_{mab}(\boldsymbol{q},-\boldsymbol{k},\boldsymbol{p}) \left[u'_a(\boldsymbol{k})u'_b(-\boldsymbol{p})\right] +\varGamma_{mab}(\boldsymbol{q},\boldsymbol{p},-\boldsymbol{k}) \left[u'_a(-\boldsymbol{p})u'_b(\boldsymbol{k})\right]\right). \end{gathered}\right\} \end{align}
The primed quantities depend on the time 
$t=s$ and 
$G'$ is Green's function. There is a difference in the treatment of Eulerian and Lagrangian theory here (Kraichnan Reference Kraichnan1965), but the final expressions for the single-time closure which we derive here are insensitive to this difference (see Bos & Bertoglio Reference Bos and Bertoglio2013). Substituting the triple correlation, the velocity modes 
$u=u^{(0)}+u^{(1)}$ and retaining the first-order terms, we obtain
\begin{equation} \frac{1}{2}\frac{\partial }{\partial t} \left\langle u_i(\boldsymbol{k})u_i(-\boldsymbol{k})\right\rangle=\sum_{i=1..6}\int_\delta \int_0^t \,\textrm{d}s F_{(i)}(k), \end{equation}with
\begin{gather} F_{(1)}=G'(k)\varGamma_{ijm}(\boldsymbol{k},-\boldsymbol{p},-\boldsymbol{q})\varGamma_{iab}(-\boldsymbol{k},\boldsymbol{p},\boldsymbol{q})\left[u'_a(-\boldsymbol{p})u'_b(-\boldsymbol{q})u_j(\boldsymbol{p})u_m(\boldsymbol{q})\right], \end{gather}
\begin{gather}F_{(2)}=G'(k)\varGamma_{ijm}(\boldsymbol{k},-\boldsymbol{p},-\boldsymbol{q})\varGamma_{iab}(-\boldsymbol{k},\boldsymbol{q},\boldsymbol{p})\left[u'_a(-\boldsymbol{q})u'_b(-\boldsymbol{p})u_j(\boldsymbol{p})u_m(\boldsymbol{q})\right], \end{gather}
\begin{gather}F_{(3)}=G'(p)\varGamma_{ijm}(\boldsymbol{k},-\boldsymbol{p},-\boldsymbol{q})\varGamma_{jab}(\boldsymbol{p},-\boldsymbol{k},\boldsymbol{q})\left[u'_a(\boldsymbol{k})u'_b(-\boldsymbol{q})u_m(\boldsymbol{q})u_i(-\boldsymbol{k})\right], \end{gather}
\begin{gather}F_{(4)}=G'(p)\varGamma_{ijm}(\boldsymbol{k},-\boldsymbol{p},-\boldsymbol{q})\varGamma_{jab}(\boldsymbol{p},\boldsymbol{q},-\boldsymbol{k})\left[u'_a(-\boldsymbol{q})u'_b(\boldsymbol{k})u_m(\boldsymbol{q})u_i(-\boldsymbol{k})\right], \end{gather}
\begin{gather}F_{(5)}=G'(q)\varGamma_{ijm}(\boldsymbol{k},-\boldsymbol{p},-\boldsymbol{q})\varGamma_{mab}(\boldsymbol{q},-\boldsymbol{k},\boldsymbol{p})\left[u'_a(\boldsymbol{k})u'_b(-\boldsymbol{p})u_j(\boldsymbol{p})u_i(-\boldsymbol{k})\right], \end{gather}
\begin{gather}F_{(6)}=G'(q)\varGamma_{ijm}(\boldsymbol{k},-\boldsymbol{p},-\boldsymbol{q})\varGamma_{mab}(\boldsymbol{q},\boldsymbol{p},-\boldsymbol{k})\left[u'_a(-\boldsymbol{p})u'_b(\boldsymbol{k})u_j(\boldsymbol{p})u_i(-\boldsymbol{k})\right], \end{gather}and using the definition
\begin{equation} \left\langle u'_i(\boldsymbol{k})u_j(-\boldsymbol{k})\right\rangle=P_{ij}(\boldsymbol{k})U(k)R'(k) \end{equation}
with 
$R'(k)$ a time-correlation function
\begin{equation} R'(k)=\frac{\left\langle u'_i(\boldsymbol{k})u_i(-\boldsymbol{k})\right\rangle}{\left\langle u_i(\boldsymbol{k})u_i(-\boldsymbol{k})\right\rangle}, \end{equation}
and 
$U(k)=E(k)/4{\rm \pi} k^2$. This can be written
\begin{gather} \frac{\partial E(k)}{\partial t} =4{\rm \pi} k^2\sum_{i=1..6}\int_\delta \int_0^t \,\textrm{d}s F_{(i)}(k), \end{gather}
\begin{gather}F_{(1)}=G'(k)R'(p)R'(q)\varGamma_{ijm}(\boldsymbol{k},-\boldsymbol{p},-\boldsymbol{q})\varGamma_{iab}(-\boldsymbol{k},\boldsymbol{p},\boldsymbol{q})P_{aj}(\boldsymbol{p})P_{bm}(\boldsymbol{q})U(p)U(q), \end{gather}
\begin{gather}F_{(2)}=G'(k)R'(p)R'(q)\varGamma_{ijm}(\boldsymbol{k},-\boldsymbol{p},-\boldsymbol{q})\varGamma_{iab}(-\boldsymbol{k},\boldsymbol{q},\boldsymbol{p})P_{bj}(\boldsymbol{p})P_{am}(\boldsymbol{q})U(p)U(q), \end{gather}
\begin{gather}F_{(3)}=G'(p)R'(k)R'(q)\varGamma_{ijm}(\boldsymbol{k},-\boldsymbol{p},-\boldsymbol{q})\varGamma_{jab}(\boldsymbol{p},-\boldsymbol{k},\boldsymbol{q})P_{ia}(\boldsymbol{k})P_{bm}(\boldsymbol{q})U(k)U(q), \end{gather}
\begin{gather}F_{(4)}=G'(p)R'(k)R'(q)\varGamma_{ijm}(\boldsymbol{k},-\boldsymbol{p},-\boldsymbol{q})\varGamma_{jab}(\boldsymbol{p},\boldsymbol{q},-\boldsymbol{k})P_{ib}(\boldsymbol{k})P_{am}(\boldsymbol{q})U(k)U(q), \end{gather}
\begin{gather}F_{(5)}=G'(q)R'(k)R'(p)\varGamma_{ijm}(\boldsymbol{k},-\boldsymbol{p},-\boldsymbol{q})\varGamma_{mab}(\boldsymbol{q},-\boldsymbol{k},\boldsymbol{p})P_{ia}(\boldsymbol{k})P_{jb}(\boldsymbol{p})U(k)U(p), \end{gather}
\begin{gather}F_{(6)}=G'(q)R'(k)R'(p)\varGamma_{ijm}(\boldsymbol{k},-\boldsymbol{p},-\boldsymbol{q})\varGamma_{mab}(\boldsymbol{q},\boldsymbol{p},-\boldsymbol{k})P_{ib}(\boldsymbol{k})P_{ja}(\boldsymbol{p})U(k)U(p). \end{gather}
All the two-time dependence is contained in the quantities 
$R', G'$. Markovianization consists here in assuming exponential time dependence for these quantities. Furthermore, assuming 
$G'=R'$ for 
$s<t$, allows us to write
\begin{equation} \int_0^t G'(k)R'(p)R'(q) \,\textrm{d}s=\varTheta(k,p,q), \end{equation}resulting in a Markovian closure.
To advance we need to contract and substitute the 
$\varGamma$s. These are defined as (see Appendix A)
\begin{equation} \varGamma_{ijm}(\boldsymbol{k},-\boldsymbol{p},-\boldsymbol{q})={-}\textrm{i} \left(\lambda \frac{\boldsymbol{p}\boldsymbol{\cdot} \boldsymbol{k}}{k^2}+ \frac{\boldsymbol{q}\boldsymbol{\cdot} \boldsymbol{k}}{k^2}\right)k_jP_{im}(\boldsymbol{k}) \end{equation}so that
\begin{gather} \varGamma_{iab}(-\boldsymbol{k},\boldsymbol{p},\boldsymbol{q})=\textrm{i} \left(\lambda \frac{\boldsymbol{p}\boldsymbol{\cdot} \boldsymbol{k}}{k^2}+ \frac{\boldsymbol{q}\boldsymbol{\cdot} \boldsymbol{k}}{k^2}\right)k_aP_{ib}(\boldsymbol{k}), \end{gather}
\begin{gather}\varGamma_{iab}(-\boldsymbol{k},\boldsymbol{q},\boldsymbol{p})=\textrm{i} \left(\lambda \frac{\boldsymbol{q}\boldsymbol{\cdot} \boldsymbol{k}}{k^2}+ \frac{\boldsymbol{p}\boldsymbol{\cdot} \boldsymbol{k}}{k^2}\right)k_aP_{ib}(\boldsymbol{k}), \end{gather}
\begin{gather}\ \varGamma_{jab}(\boldsymbol{p},-\boldsymbol{k},\boldsymbol{q})={-}\textrm{i} \left(\lambda \frac{\boldsymbol{k}\boldsymbol{\cdot} \boldsymbol{p}}{p^2}+ \frac{-\boldsymbol{q}\boldsymbol{\cdot} \boldsymbol{p}}{p^2}\right)p_aP_{jb}(\boldsymbol{p}), \end{gather}
\begin{gather}\varGamma_{jab}(\boldsymbol{p},\boldsymbol{q},-\boldsymbol{k})={-}\textrm{i} \left(\lambda \frac{-\boldsymbol{q}\boldsymbol{\cdot} \boldsymbol{p}}{p^2}+ \frac{\boldsymbol{k}\boldsymbol{\cdot} \boldsymbol{p}}{p^2}\right)p_aP_{jb}(\boldsymbol{p}), \end{gather}
\begin{gather}\varGamma_{mab}(\boldsymbol{q},-\boldsymbol{k},\boldsymbol{p})={-}\textrm{i} \left(\lambda \frac{\boldsymbol{k}\boldsymbol{\cdot} \boldsymbol{q}}{q^2}+ \frac{-\boldsymbol{p}\boldsymbol{\cdot} \boldsymbol{q}}{q^2}\right)q_aP_{mb}(\boldsymbol{q}), \end{gather}
\begin{gather}\varGamma_{mab}(\boldsymbol{q},\boldsymbol{p},-\boldsymbol{k})={-}\textrm{i} \left(\lambda \frac{-\boldsymbol{p}\boldsymbol{\cdot} \boldsymbol{q}}{q^2}+ \frac{\boldsymbol{k}\boldsymbol{\cdot} \boldsymbol{q}}{q^2}\right)q_aP_{mb}(\boldsymbol{q}) \end{gather}
and the product of the 
$\varGamma$s is
\begin{equation} \left.\begin{gathered} \varGamma_{ijm}(\boldsymbol{k},-\boldsymbol{p},-\boldsymbol{q})\varGamma_{iab}(-\boldsymbol{k},\boldsymbol{p}, \boldsymbol{q})= \left(\lambda \frac{\boldsymbol{p}\boldsymbol{\cdot} \boldsymbol{k}}{k^2}+ \frac{\boldsymbol{q}\boldsymbol{\cdot} \boldsymbol{k}}{k^2}\right) \left(\lambda \frac{\boldsymbol{p} \boldsymbol{\cdot} \boldsymbol{k}}{k^2}+ \frac{\boldsymbol{q} \boldsymbol{\cdot} \boldsymbol{k}}{k^2}\right)k_aP_{ib}(\boldsymbol{k}) k_jP_{im}(\boldsymbol{k})\\ \varGamma_{ijm}(\boldsymbol{k},-\boldsymbol{p},-\boldsymbol{q})\varGamma_{iab}(-\boldsymbol{k},\boldsymbol{q}, \boldsymbol{p})= \left(\lambda \frac{\boldsymbol{q}\boldsymbol{\cdot} \boldsymbol{k}}{k^2}+ \frac{\boldsymbol{p}\boldsymbol{\cdot} \boldsymbol{k}}{k^2}\right) \left(\lambda \frac{\boldsymbol{p}\boldsymbol{\cdot }\boldsymbol{k}}{k^2}+ \frac{\boldsymbol{q}\boldsymbol{\cdot} \boldsymbol{k}}{k^2}\right)k_jk_aP_{ib}(\boldsymbol{k}) P_{im}(\boldsymbol{k})\\ \varGamma_{ijm}(\boldsymbol{k},-\boldsymbol{p},-\boldsymbol{q})\varGamma_{jab}(\boldsymbol{p},-\boldsymbol{k},\boldsymbol{q})={-} \left(\lambda \frac{\boldsymbol{k}\boldsymbol{\cdot} \boldsymbol{p}}{p^2}+ \frac{-\boldsymbol{q}\boldsymbol{\cdot} \boldsymbol{p}}{p^2}\right) \left(\lambda \frac{\boldsymbol{p}\boldsymbol{\cdot} \boldsymbol{k}}{k^2}+ \frac{\boldsymbol{q}\boldsymbol{\cdot} \boldsymbol{k}}{k^2}\right)k_jp_aP_{jb}(\boldsymbol{p})P_{im}( \boldsymbol{k})\\ \varGamma_{ijm}(\boldsymbol{k},-\boldsymbol{p},-\boldsymbol{q})\varGamma_{jab}(\boldsymbol{p}, \boldsymbol{q},-\boldsymbol{k})={-} \left(\lambda \frac{-\boldsymbol{q}\boldsymbol{\cdot} \boldsymbol{p}}{p^2}+ \frac{\boldsymbol{k}\boldsymbol{\cdot} \boldsymbol{p}}{p^2}\right) \left(\lambda \frac{\boldsymbol{p}\boldsymbol{\cdot} \boldsymbol{k}}{k^2}+ \frac{\boldsymbol{q}\boldsymbol{\cdot} \boldsymbol{k}}{k^2}\right)k_jp_aP_{jb}(\boldsymbol{p}) P_{im}(\boldsymbol{k})\\ \varGamma_{ijm}(\boldsymbol{k},-\boldsymbol{p},-\boldsymbol{q})\varGamma_{mab}(\boldsymbol{q},-\boldsymbol{k}, \boldsymbol{p})={-} \left(\lambda \frac{\boldsymbol{k}\boldsymbol{\cdot} \boldsymbol{q}}{q^2}+ \frac{-\boldsymbol{p}\boldsymbol{\cdot} \boldsymbol{q}}{q^2}\right) \left(\lambda \frac{\boldsymbol{p} \boldsymbol{\cdot} \boldsymbol{k}}{k^2}+ \frac{\boldsymbol{q}\boldsymbol{\cdot} \boldsymbol{k}}{k^2}\right)k_jq_aP_{mb}(\boldsymbol{q})P_{im}(\boldsymbol{k})\\ \varGamma_{ijm}(\boldsymbol{k},-\boldsymbol{p},-\boldsymbol{q})\varGamma_{mab}(\boldsymbol{q},\boldsymbol{p},-\boldsymbol{k})={-} \left(\lambda \frac{-\boldsymbol{p}\boldsymbol{\cdot} \boldsymbol{q}}{q^2}+ \frac{\boldsymbol{k}\boldsymbol{\cdot} \boldsymbol{q}}{q^2}\right) \left(\lambda \frac{\boldsymbol{p}\boldsymbol{\cdot} \boldsymbol{k}}{k^2}+ \frac{\boldsymbol{q}\boldsymbol{\cdot} \boldsymbol{k}}{k^2}\right)k_jq_aP_{mb}(\boldsymbol{q})P_{im}(\boldsymbol{k}) \end{gathered}\right\} \end{equation}
The terms in parentheses all yield the value 1 for 
$\lambda =1$. For 
$\lambda =0$ we have, using the definitions 
$\boldsymbol {q}\boldsymbol {\cdot } \boldsymbol {k}=kqy$, 
$\boldsymbol {p}\boldsymbol {\cdot } \boldsymbol {k}=pkz$, 
$\boldsymbol {q}\boldsymbol {\cdot } \boldsymbol {p}=-pqx$,
\begin{gather} \varGamma_{ijm}(\boldsymbol{k},-\boldsymbol{p},-\boldsymbol{q})\times\varGamma_{iab}(-\boldsymbol{k}, \boldsymbol{p},\boldsymbol{q})=y^2\left(\frac{q}{k}\right)^2 k_aP_{ib}(\boldsymbol{k}) k_jP_{im}(\boldsymbol{k}), \end{gather}
\begin{gather}\varGamma_{ijm}(\boldsymbol{k},-\boldsymbol{p},-\boldsymbol{q})\times\varGamma_{iab}(-\boldsymbol{k}, \boldsymbol{q},\boldsymbol{p})= yz\frac{pq}{k^2} k_aP_{ib}(\boldsymbol{k}) k_jP_{im}(\boldsymbol{k}), \end{gather}
\begin{gather}\varGamma_{ijm}(\boldsymbol{k},-\boldsymbol{p},-\boldsymbol{q})\times\varGamma_{jab}(\boldsymbol{p},-\boldsymbol{k},\boldsymbol{q})={-}xy\frac{q^2}{pk}p_aP_{jb}(\boldsymbol{p}) k_jP_{im}(\boldsymbol{k}), \end{gather}
\begin{gather}\varGamma_{ijm}(\boldsymbol{k},-\boldsymbol{p},-\boldsymbol{q})\times\varGamma_{jab}(\boldsymbol{p},\boldsymbol{q},-\boldsymbol{k})={-}yz\frac{q}{p}p_aP_{jb}(\boldsymbol{p}) k_jP_{im}(\boldsymbol{k}), \end{gather}
\begin{gather}\varGamma_{ijm}(\boldsymbol{k},-\boldsymbol{p},-\boldsymbol{q})\times\varGamma_{mab}(\boldsymbol{q},-\boldsymbol{k},\boldsymbol{p})={-}xy\frac{p}{k} q_aP_{mb}(\boldsymbol{q}) k_jP_{im}(\boldsymbol{k}), \end{gather}
\begin{gather}\varGamma_{ijm}(\boldsymbol{k},-\boldsymbol{p},-\boldsymbol{q})\times\varGamma_{mab}(\boldsymbol{q},\boldsymbol{p},-\boldsymbol{k})={-}y^2 q_aP_{mb}(\boldsymbol{q}) k_jP_{im}(\boldsymbol{k}), \end{gather}thus
\begin{gather} F_{(1)}=y^2\left(\frac{q}{k}\right)^2 k_aP_{ib}(\boldsymbol{k}) k_jP_{im}(\boldsymbol{k})P_{aj}(\boldsymbol{p})P_{bm}(\boldsymbol{q})\varTheta(k,p,q)U(p)U(q), \end{gather}
\begin{gather}F_{(2)}= yz\frac{pq}{k^2} k_aP_{ib}(\boldsymbol{k}) k_jP_{im}(\boldsymbol{k})P_{bj}(\boldsymbol{p})P_{am}(\boldsymbol{q})\varTheta(k,p,q)U(p)U(q), \end{gather}
\begin{gather}F_{(3)}={-}xy\frac{q^2}{pk}p_aP_{jb}(\boldsymbol{p}) k_jP_{im}(\boldsymbol{k})P_{ia}(\boldsymbol{k})P_{bm}(\boldsymbol{q})\varTheta(k,p,q)U(k)U(q), \end{gather}
\begin{gather}F_{(4)}={-}yz\frac{q}{p}p_aP_{jb}(\boldsymbol{p}) k_jP_{im}(\boldsymbol{k})P_{ib}(\boldsymbol{k})P_{am}(\boldsymbol{q})\varTheta(k,p,q)U(k)U(q), \end{gather}
\begin{gather}F_{(5)}={-}xy\frac{p}{k} q_aP_{mb}(\boldsymbol{q}) k_jP_{im}(\boldsymbol{k})P_{ia}(\boldsymbol{k})P_{jb}(\boldsymbol{p})\varTheta(k,p,q)U(k)U(p), \end{gather}
\begin{gather}F_{(6)}={-}y^2 q_aP_{mb}(\boldsymbol{q}) k_jP_{im}(\boldsymbol{k})P_{ib}(\boldsymbol{k})P_{ja}(\boldsymbol{p})\varTheta(k,p,q)U(k)U(p), \end{gather}which is written as
\begin{gather} F_{(1)}=y^2\left(\frac{q}{k}\right)^2 [k_1k_2P_{34}^kP_{12}^pP_{34}^q] \varTheta(k,p,q)U(p)U(q), \end{gather}
\begin{gather}F_{(2)}= yz\frac{pq}{k^2} [k_1k_2P_{34}^kP_{13}^pP_{24}^q] \varTheta(k,p,q)U(p)U(q), \end{gather}
\begin{gather}F_{(3)}={-}xy\frac{q^2}{pk}[k_1p_2P_{23}^kP_{14}^pP_{34}^q]\varTheta(k,p,q)U(k)U(q), \end{gather}
\begin{gather}F_{(4)}={-}yz\frac{q}{p}[k_1p_2P_{34}^kP_{13}^pP_{24}^q]\varTheta(k,p,q)U(k)U(q), \end{gather}
\begin{gather}F_{(5)}={-}xy\frac{p}{k} [k_1q_2P_{23}^kP_{14}^pP_{34}^q]\varTheta(k,p,q)U(k)U(p), \end{gather}
\begin{gather}F_{(6)}={-}y^2 [k_1q_2P_{34}^kP_{12}^pP_{34}^q]\varTheta(k,p,q)U(k)U(p). \end{gather}
Note that we prefer for purely technical reasons to replace wave vector dependence by superscripts and the indices by numbers, which allows us to more easily organize and order the different terms. Summation over repeated indices is still assumed. A procedure, explained in Leslie (Reference Leslie1973), allows us to rewrite the integral over wave vectors as a scalar integral over 
$p,q$ space. This allows us to rewrite the expression as
\begin{align} \frac{\partial E(k)}{\partial t}&=\frac{1}{2}\int_\varDelta \frac{\textrm{d}p}{p}\frac{\textrm{d}q}{q}\varTheta(k,p,q)\nonumber\\ &\quad \times \left[f^{\lambda}_{(1)} X_{(1)}+f^{\lambda}_{(2)} X_{(2)}\right]k^3 E(p)E(q) \end{align}
\begin{align} &\quad -\left[f^{\lambda}_{(3)} X_{(3)}+ f^{\lambda}_{(4)} X_{(4)}\right]p^3 E(k)E(q) \end{align}
\begin{align} &\quad -\left[f^{\lambda}_{(5)} X_{(5)}+ f^{\lambda}_{(6)} X_{(6)}\right]q^3 E(p)E(k), \end{align}
where the 
$X$ are
\begin{gather} X_{(1)}=k^{{-}2}k_1k_2P_{34}^kP_{12}^pP_{34}^q=(1-z^2)(1+y^2), \end{gather}
\begin{gather}X_{(2)}=k^{{-}2}k_1k_2P_{34}^kP_{13}^pP_{24}^q={-}xyz-y^2z^2, \end{gather}
\begin{gather}X_{(3)}=(kp)^{{-}1}k_1p_2P_{23}^kP_{14}^pP_{34}^q=xy(1-z^2), \end{gather}
\begin{gather}X_{(4)}=(kp)^{{-}1}k_1p_2P_{34}^kP_{13}^pP_{24}^q=z({-}y^2-xyz), \end{gather}
\begin{gather}X_{(5)}=(kq)^{{-}1}k_1q_2P_{23}^kP_{14}^pP_{34}^q=y({-}x^2-xyz), \end{gather}
\begin{gather}X_{(6)}=(kq)^{{-}1}k_1q_2P_{34}^kP_{12}^pP_{34}^q=(y+zx)(1+y^2). \end{gather}
As a consistency check we can assess the case 
$\lambda =1$ leading to the classical EDQNM closure.
