1. Introduction
Despite more than a century of concentrated effort, fluid turbulence remains steadfast as the oldest unsolved problem of classical physics. Much of the progress in developing our understanding of turbulence can be traced to the work of Kolmogorov and his three 1941 papers (Kolmogorov Reference Kolmogorov1941a,Reference Kolmogorovb,Reference Kolmogorovc), in which what has come to be known as the K41 theory was first described. The results in these papers are derived for an idealised form of turbulence known as homogeneous and isotropic turbulence (HIT), however, they are also remarkably applicable to real world flows under certain conditions. One of the most important predictions of the K41 theory, valid at sufficiently high Reynolds number, is the existence of a range of intermediate sized eddies in the flow referred to as the inertial range, characterised by scale invariance and a constant energy flux. This scale invariance manifests itself clearly in the power law form of the K41 energy spectrum in the inertial range
\begin{equation} E(k) = C\varepsilon^{{2}/{3}} k^{-({5}/{3})}, \end{equation}
where k is the wavenumber and 
$\varepsilon$ is the constant energy flux, which for stationary turbulence will be equal to the rate of viscous energy dissipation, and 
$C$ is a universal constant.
The scale invariance of the inertial range is reminiscent of that seen in critical phenomena close to the critical point. Following this line of argument, there have been numerous analogies comparing turbulence to critical phenomena (Nelkin Reference Nelkin1974; De Gennes Reference De Gennes1975; Bramwell, Holdsworth & Pinton Reference Bramwell, Holdsworth and Pinton1998; Aji & Goldenfeld Reference Aji and Goldenfeld2001; Yakhot Reference Yakhot2001; Giuliani, Jensen & Yakhot Reference Giuliani, Jensen and Yakhot2002; Frisch et al. Reference Frisch, Pomyalov, Procaccia and Ray2012). A salient feature of many such critical systems is the existence of an upper critical dimension, above which fluctuations are suppressed and the mean-field theory values for critical exponents become exact. These ideas have their roots in the work of Ginzburg (Reference Ginzburg1960), as well as that of Wilson and Fisher, in the application of renormalisation group methods to critical phenomena (Wilson Reference Wilson1971; Wilson & Fisher Reference Wilson and Fisher1972). For turbulence, a case can be made that the K41 theory, since it uses the mean energy flux in the form of the inertial range energy spectrum, is in fact a kind of mean-field theory itself (Siggia Reference Siggia1977; Bell & Nelkin Reference Bell and Nelkin1978).
This interpretation of K41 is interesting in light of the measurement of deviations from the exponents predicted by K41 for both the energy spectrum and the structure functions. Both intermittency and finite Reynolds number effects have been theorised as being responsible for such deviations (Kolmogorov Reference Kolmogorov1962; Frisch, Sulem & Nelkin Reference Frisch, Sulem and Nelkin1978; Benzi et al. Reference Benzi, Paladin, Parisi and Vulpiani1984; McComb Reference McComb2014), which leads to (1.1) being re-expressed in the form
\begin{equation} E(k) \propto k^{-({5}/{3}) - \mu}, \end{equation}
where 
$\mu$ is the anomalous exponent. From here an analogy can be drawn once more to critical phenomena, in which a similar anomalous exponent, which vanishes for mean-field theory, is seen when looking at two-point correlation functions. Naturally, this has led to speculation about whether an upper critical dimension for turbulence exists and, if so, what its associated value would be (Rose & Sulem Reference Rose and Sulem1978; Liao Reference Liao1990, Reference Liao1991; Nelkin Reference Nelkin2001). Results from a recent study (Berera, Ho & Clark Reference Berera, Ho and Clark2020) which performed direct numerical simulation (DNS) of four spatial dimensional HIT found, amongst other results, a suppression of energy fluctuations in going from three to four dimensions, which has raised further interesting questions relating to a critical dimension in turbulence. There have also been claims related to, and a small number of studies investigating, the possibility of simplification in infinite dimensions (Kraichnan Reference Kraichnan1974; Fournier & Frisch Reference Fournier and Frisch1978; Fournier, Frisch & Rose Reference Fournier, Frisch and Rose1978).
To date, only a handful of DNS studies of turbulence in spatial dimensions greater than three have been carried out. The work of Suzuki et al. (Reference Suzuki, Nakano, Takahashi and Gotoh2005) came first with the same group following up on this work in Gotoh et al. (Reference Gotoh, Watanabe, Shiga, Nakano and Suzuki2007). These insightfully motivated studies were focussed on the effect of the spatial dimension on the intermittent nature of turbulence. In both cases due to the computational limitations of the time only decaying turbulence was studied. Here, they found an increased efficiency of the energy transfer in four dimensions compared with three, as well as increased anomalous scaling of the longitudinal structure functions. More recently, in the DNS study of Berera et al. (Reference Berera, Ho and Clark2020), which looked at the stationary case at a higher resolution, various measurements performed pointed to an increased tendency for energy to be transferred from large to small scales in four dimensions, confirming what had been seen in Gotoh et al. (Reference Gotoh, Watanabe, Shiga, Nakano and Suzuki2007), potentially driven by an increase in vortex stretching. This interpretation was based on finding a higher velocity derivative skewness. Whether this trend of increased forward energy transfer continues into higher spatial dimensions is an interesting question, as it may indicate the possibility of the turbulent dynamics being simplified in higher dimensions. Presently, the computational cost of performing DNS of higher-dimensional turbulence is beyond even the largest of supercomputers. Indeed, we are only aware of one study which has investigated beyond four-dimensional turbulence via DNS. In the study by Yamamoto et al. (Reference Yamamoto, Shimizu, Inoshita, Nakano and Gotoh2012), a five-dimensional simulation was carried out, though this involved a relatively large lattice spacing, meaning the results are only for very low Reynolds number values where the conditions of K41 are not met. With these computational considerations in mind, we are forced to turn to closure approximations in order to conduct a feasible study.
Closure approximations have their roots in quantum field theory (QFT). Initially, work was pioneered by Kraichnan (Reference Kraichnan1959), Wyld (Reference Wyld1961) and Edwards (Reference Edwards1964) in employing QFT methods to develop a perturbation theory for the Navier–Stokes equations, and this subsequently led to various approximation schemes. In this investigation, we will make use of the eddy damped quasi-normal Markovian (EDQNM) closure, first described by Orszag (Reference Orszag1970) as a method of achieving realisability in the quasi-normal approximation (Millionshchikov Reference Millionshchikov1941). Before proceeding with the details of this method and our calculations, it is interesting to note that QFT is also a subject in which the behaviour of systems in different dimensions has been an area of sustained interest in systems including string theory, gauge theory, and anti-de Sitter/conformal field theory correspondence. Hence, aside from the computational tools, this area has helped develop for the field of fluid turbulence, there is also relevance in appreciating this conceptual point and thus in placing more focus on understanding the dimensional behaviour of fluid turbulence.
The EDQNM closure scheme has seen widespread use in both two- and three-dimensional turbulence (see Lesieur (Reference Lesieur1987) for an in-depth review), where it has produced numerous qualitative results. The EDQNM approximation allows investigation to very high Reynolds number flows at relatively low computational cost and has the added benefit that extension to any dimension incurs no additional computational expense. The EDQNM closure is compatible with the Kolmogorov energy spectrum and is well suited for the study of energy transfer in isotropic turbulence. Additionally, it was noted by Orszag (Reference Orszag1974) that the quasi-normal approximation is analogous to the random phase approximation of many-body physics, which is also closely linked to the Gaussian approximation. Therefore, if higher-dimensional turbulence shows a systematic improvement in agreement with the EDQNM approximation, it may in its own right be an indicator towards a simplification in the turbulent dynamics.
If forward energy transfer does indeed become stronger with higher spatial dimension, this may result in an increased bottleneck effect (Falkovich Reference Falkovich1994). This effect manifests itself as a pile up of energy in the near dissipative range of the flow and has been observed both experimentally (Mestayer Reference Mestayer1982; Saddoughi & Veeravalli Reference Saddoughi and Veeravalli1994) and numerically (Kerr Reference Kerr1990). It has been suggested by Herring et al. (Reference Herring, Schertzer, Lesieur, Newman, Chollet and Larcheveque1982) that this effect is a result of viscosity suppressing the nonlinear transfer of energy to the smallest scales. Hence, by varying the spatial dimension of the system it is possible to investigate these claims of viscous energy transfer suppression.
In Berera et al. (Reference Berera, Ho and Clark2020) a DNS dataset of unprecedented size for four spatial dimensions was developed. It reached a box size of 
$512^{4}$, which in terms of computational demands is similar to 
$4096^{3}$, so amongst the larger DNS datasets. Moreover, this was a forced simulation, run for a very long time to achieve good equilibration and adequate time for robust sampling of the data. The simulation focused on studying large scale properties of the four-dimensional turbulent state. In particular it examined the anomalous dissipation and total energy fluctuation, comparing behaviour between three and four spatial dimensions. The study found a significant suppression of the energy fluctuation in four dimensions, thus having some qualitative similarities to behaviour found in critical phenomenon. Additionally, the study found an increased velocity derivative skewness and asymptotic dissipation rate in four spatial dimensions compared to three. This is consistent with the interpretation on an enhanced forward energy transfer in four dimensions when compared to three dimensions. The study by Berera et al. (Reference Berera, Ho and Clark2020) also served the purpose of providing a baseline of results that will be useful for comparison in any future DNS studies in four spatial dimensions. For all these reasons it is important to have an independent check of the very new type of results seen in that study. Where possible this paper will examine these quantities in four spatial dimensions using the EDQNM approximation independently confirming several of the DNS results in Berera et al. (Reference Berera, Ho and Clark2020). We will then go further and examine similar results in higher dimensions. The structure of this paper is as follows: in § 2 we outline a number of derivations for 
$d$-dimensional turbulence in the Navier–Stokes equations. Section 3 introduces the EDQNM closure model used for this work. Section 4 presents the results of our numerical study and finally, § 5 discusses the interpretation and possible implications of these results.
2. Theory
The incompressible Navier–Stokes equations (NSE) can be expressed for spatial dimension 
$d \geq 2$ as
\begin{equation} \left.\begin{aligned} \partial_t \boldsymbol{u} + \boldsymbol{u} \cdot \boldsymbol{\nabla} \boldsymbol{u} & ={-} \boldsymbol{\nabla} P + \nu \nabla^{2} \boldsymbol{u}, \\ \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} & = 0. \end{aligned}\right\} \end{equation}
In the above, 
$\boldsymbol {u}( \boldsymbol {x},t)$ is the velocity field, 
$P( \boldsymbol {x},t)$ is the pressure field, 
$\nu$ is the kinematic viscosity and 
$\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {u} = 0$ is the incompressibility condition, which allows us to set the fluid density to unity.
Here, we discuss the relevant quantities and definitions that are utilised throughout the rest of the paper. For brevity in the following we will drop the explicit time dependence. For our purposes, we are primarily interested in the second- and third-order two-point velocity correlations, which are given by
\begin{equation} \begin{aligned} C_{\alpha \beta}( \boldsymbol{r}) & = \langle u_{\alpha}( \boldsymbol{x}) u_{\beta}( \boldsymbol{x+r})\rangle, \\ C_{\alpha \beta \gamma}( \boldsymbol{r}) & = \langle u_{\alpha}( \boldsymbol{x}) u_{\beta}( \boldsymbol{x}) u_{\gamma}( \boldsymbol{x+r})\rangle, \end{aligned} \end{equation}
with 
$\alpha , \beta , \gamma = 1,\ldots ,d$. More specifically, we will be focussed on the second- and third-order longitudinal correlations defined as
\begin{equation} \begin{aligned} C_{LL} & = \frac{r_{\alpha} r_{\beta}}{r^{2}}C_{\alpha \beta}( \boldsymbol{r}) = u^{2} f(r), \\ C_{LL,L} & = \frac{r_{\alpha} r_{\beta} r_{\gamma}}{r^{3}}C_{\alpha \beta \gamma}( \boldsymbol{r}) = u^{3} K(r), \end{aligned} \end{equation}
where 
$f(r)$ and 
$K(r)$ are scalar correlation functions and 
$u$ is the root-mean-square velocity. These functions are intimately related to the longitudinal structure functions of the same order. It will prove useful to introduce the Fourier transform of (2.1)
\begin{equation} \left(\partial_t + \nu k^{2}\right) u_{\alpha}( \boldsymbol{k}) = \frac{1}{2i}P_{\alpha \beta \gamma}( \boldsymbol{k}) \int \mathrm{d} \boldsymbol{p} \, u_{\beta}( \boldsymbol{p})u_{\gamma}( \boldsymbol{k-p}), \end{equation}
where 
$P_{\alpha \beta \gamma }( \boldsymbol {k}) = k_{\beta } P_{\alpha \gamma }( \boldsymbol {k}) + k_{\gamma }P_{\alpha \beta }( \boldsymbol {k})$ is the inertial transfer operator and 
$P_{\alpha \beta }( \boldsymbol {k}) = \delta _{\alpha \beta } - k_{\alpha } k_{\beta }/k^{2}$ is the projection operator which imposes the incompressibility condition. Homogeneity requires that the corresponding second-order velocity correlation in Fourier space takes the form
\begin{equation} C_{\alpha \beta}( \boldsymbol{k}) = \langle u_{\alpha}( \boldsymbol{k})u_{\beta}( \boldsymbol{-k}) \rangle. \end{equation}
It can be shown that this correlator is related to the energy spectrum, 
$E(k)$, through
\begin{equation} \langle u_{\alpha}( \boldsymbol{k})u_{\beta}( \boldsymbol{-k}) \rangle = \frac{2P_{\alpha \beta}( \boldsymbol{k})E(k)}{(d-1)A_{d}k^{d-1}}, \end{equation}
where 
$A_{d} = 2( ({\rm \pi} )^{d/2}/\varGamma (d/2))$ is the surface area of a 
$d$-dimensional unit sphere. We can then form an equation for 
$E(k)$
\begin{align} \left(\partial_t + 2\nu k^{2}\right)E(k) &= \frac{\textrm{i} A_{d} k^{d-1}}{2} P_{\alpha \beta \gamma}( \boldsymbol{k})\mathcal{A}_{\alpha \beta \gamma}( \boldsymbol{k}) \nonumber\\ &= T(k). \end{align}
Here, 
$T(k)$ is the energy transfer spectrum and 
$\mathcal {A}_{\alpha \beta \gamma }( \boldsymbol {k})$ is defined as
\begin{equation} \mathcal{A}_{\alpha \beta \gamma}( \boldsymbol{k}) = \int \mathrm{d} \boldsymbol{p} \, C_{\alpha \beta \gamma}( \boldsymbol{k}, \boldsymbol{p},- \boldsymbol{k}- \boldsymbol{p}), \end{equation}
where 
$C_{\alpha \beta \gamma }( \boldsymbol {k}, \boldsymbol {p},- \boldsymbol {k} - \boldsymbol {p}) = \langle u_{\alpha }( \boldsymbol {k})u_{\beta }( \boldsymbol {p}) u_{\gamma }(- \boldsymbol {k}- \boldsymbol {p})\rangle$ is the spectral third-order velocity field moment. More detailed derivations and interpretation of these quantities can be found in McComb (Reference McComb1990).
2.1. Second-order structure function
The dimensional dependence of HIT is well elucidated in the form of the velocity correlation functions, and thus in turn the structure functions. The longitudinal structure function of order 
$n$ is given by
\begin{equation} S_{n}^{(d)} = \left\langle \boldsymbol{u}( \boldsymbol{x}+ \boldsymbol{r}) - \boldsymbol{u}( \boldsymbol{x}) \right\rangle \cdot \frac{ \boldsymbol{r}}{r}. \end{equation}
Now, we derive the relationship between the second- and third-order longitudinal structure functions and the energy and transfer spectra, respectively. As a result, we are then able to evaluate these structure functions in our numerical EDQNM results. The method used here is the 
$d$-dimensional extension of that used in Bos et al. (Reference Bos, Chevillard, Scott and Rubinstein2012). We begin by considering the second-order longitudinal structure function which can be expressed as
\begin{equation} S_{2}^{(d)}(r) = 2\left(u^{2} - C_{LL}\right) = 2\left( u^{2} - \frac{r_{\alpha} r_{\beta}}{r^{2}}C_{\alpha \beta}( \boldsymbol{r}) \right). \end{equation}
These correlations are related to their spectral analogues via a 
$d$-dimensional inverse Fourier transform
\begin{equation} C_{\alpha \beta}( \boldsymbol{r}) = \int \mathrm{d} \boldsymbol{k} \, C_{\alpha \beta}( \boldsymbol{k}) \,\textrm{e}^{\textrm{i} \boldsymbol{k}\cdot \boldsymbol{r}} = \int \mathrm{d} \boldsymbol{k} \, \frac{2P_{\alpha \beta}( \boldsymbol{k})E(k)}{(d-1)A_{d}k^{d-1}} \textrm{e}^{\textrm{i} \boldsymbol{k}\cdot \boldsymbol{r}}. \end{equation}In HIT, these transforms are simplified by the fact that the correlators must be spherically symmetric, thus we have
\begin{align} C_{LL} &= \int \mathrm{d} \boldsymbol{k} \frac{2(1-\cos^{2} \theta)E(k)}{(d-1)S_{d}k^{d-1}} \textrm{e}^{\textrm{i} \boldsymbol{k}\cdot \boldsymbol{r}} \nonumber\\ &= \int_{0}^{\infty} \mathrm{d}k \, \frac{2E(k)}{(d-1)A_{d}} A_{d-1} \int_{0}^{\rm \pi}\mathrm{d}\theta \, exp({\textrm{i}kr\cos\theta})\sin^{d-2} \theta(1-\cos^{2} \theta) \nonumber\\ &= 2^{{d}/{2}}\varGamma\left(\frac{d}{2}\right)\int_{0}^{\infty} \mathrm{d}k \, E(k)\frac{J_{{d}/{2}}(kr)}{(kr)^{{d}/{2}}}. \end{align}
Where 
$\theta$ is defined as the angle between 
$\boldsymbol {k}$ and 
$\boldsymbol {r}$, such that 
$\boldsymbol {k} \cdot \boldsymbol {r} = kr\cos \theta$ and hence 
$r_\alpha r_\beta P_{\alpha \beta }( \boldsymbol {k}) = r^{2} (1-\cos ^{2} \theta )$, and 
$J_{n}(x)$ is the 
$n$th-order Bessel function of the first kind. Additionally, for 
$u^{2}$ we find
\begin{equation} u^{2} = \frac{2}{d}\int_{0}^{\infty}\mathrm{d}k \,E(k). \end{equation}Inserting these two results into (2.10) yields
\begin{equation} S_{2}^{(d)}(r) = 2\int_{0}^{\infty}\mathrm{d}k \, E(k) \left[\frac{2}{d} -2^{{d}/{2}}\varGamma\left(\frac{d}{2}\right) \frac{J_{{d}/{2}}(kr)}{(kr)^{{d}/{2}}}\right], \end{equation}
thus we have related 
$S_{2}^{(d)}(r)$ to 
$E(k)$. It can be verified that this gives the standard results for two and three dimensions as found in Davidson (Reference Davidson2015). Additionally, this result allows us to determine the integral length scale in 
$d$-dimensional HIT. The integral length scale is defined (Batchelor Reference Batchelor1953) as
\begin{equation} L_d = \int_{0}^{\infty} \mathrm{d}r\, f(r), \end{equation}
where 
$f(r)$ is the scalar correlation function defined in (2.3), such that 
$f(r) = C_{LL}/u^{2}$. Therefore, for the 
$d$-dimensional case we have
\begin{equation} L_{d} = \frac{2^{{d}/{2}}\varGamma\left(\dfrac{d}{2}\right)}{u^{2}} \int_{0}^{\infty} \mathrm{d}k \, E(k) \int_{0}^{\infty} \mathrm{d}r \, \frac{J_{{d}/{2}}(kr)}{(kr)^{{d}/{2}}} = \frac{\varGamma\left(\dfrac{d}{2}\right) \sqrt{\rm \pi}}{\varGamma\left(\dfrac{d+1}{2}\right)u^{2}}\int_{0}^{\infty}\mathrm{d}k \, E(k)k^{{-}1}. \end{equation}
We will find this expression useful when defining our integral scale Reynolds number. It is also possible to generalise the Taylor microscale (Taylor Reference Taylor1935), 
$\lambda _d$, which gives the average size of the dissipative eddies, to 
$d$-dimensions. This length scale is defined through fitting a parabola to the small 
$r$ expansion of the scalar longitudinal correlation function 
$f(r)$, i.e.
\begin{equation} f(r) = 1 - \frac{r^{2}}{2\lambda_{d}^{2}} + {O}(r^{4}). \end{equation}
From (2.3) and (2.12) we find through expansion for small 
$r$
\begin{equation} f(r) = \frac{C_{LL}}{u^{2}} = 1 -\frac{r^{2}}{d(d+2)u^{2}}\int_{0}^{\infty} \mathrm{d}k \, k^{2} E(k) + {O}(r^{4}), \end{equation}where if we recall that
\begin{equation} \varepsilon = 2 \nu\int_{0}^{\infty} \mathrm{d}k \, k^{2} E(k), \end{equation}we then have
\begin{equation} f(r) = 1 - \frac{\varepsilon r^{2}}{2d(d+2)\nu u^{2}} + {O}(r^{4}). \end{equation}
Hence, the Taylor microscale in 
$d$-dimensions is given by
\begin{equation} \lambda_{d} = \sqrt{\frac{d(d+2) \nu }{\varepsilon}}u. \end{equation}2.2. Third-order structure function
A similar analysis can be performed for the third-order structure function, 
$S_{3}(r)$, whereby it is related to the transfer spectrum. We begin with the relation
\begin{equation} S_{3}^{(d)}(r) = 6C_{LL,L} = 6\frac{r_{\alpha} r_{\beta} r_{\gamma}}{r^{3}}C_{\alpha \beta \gamma}( \boldsymbol{r}), \end{equation}and find that upon Fourier transform we have
\begin{equation} \mathcal{F}[C_{\alpha \beta \gamma}( \boldsymbol{r})] = \mathcal{A}_{\alpha \beta \gamma}( \boldsymbol{k}), \end{equation}
with 
$A_{\alpha \beta \gamma }( \boldsymbol {k})$ defined as in (2.8). From here, we can observe that
\begin{equation} S_{3}^{(d)}(r) = 6\frac{r_{\alpha} r_{\beta} r_{\gamma}}{r^{3}} \int \mathrm{d} \boldsymbol{k} \, \mathcal{A}_{\alpha \beta \gamma}( \boldsymbol{k}) \,\textrm{e}^{\textrm{i} \boldsymbol{k}\cdot \boldsymbol{r}}. \end{equation}Now all that remains is to express this in terms of the transfer function and perform the Fourier integrals. From (2.7) we have
\begin{equation} \frac{\textrm{i}A_{d} k^{d-1}}{2} P_{\alpha \beta \gamma}( \boldsymbol{k})\mathcal{A}_{\alpha \beta \gamma}( \boldsymbol{k}) = T(k), \end{equation}
and, as 
$\mathcal {A}_{\alpha \beta \gamma }( \boldsymbol {k})$ is a third rank solenoidal tensor symmetric in the indices 
$\beta$ and 
$\gamma$, we must have 
$\mathcal {A}_{\alpha \beta \gamma }( \boldsymbol {k}) = P_{\alpha \beta \gamma }( \boldsymbol {k})\mathcal {A}(k)$. Thus by evaluating the product 
$P_{\alpha \beta \gamma }( \boldsymbol {k})P_{\alpha \beta \gamma }( \boldsymbol {k}) = 2(d-1)k^{2}$ we have
\begin{equation} \mathcal{A}_{\alpha \beta \gamma}( \boldsymbol{k}) = \frac{P_{\alpha \beta \gamma}( \boldsymbol{k}) T(k)}{\textrm{i}(d-1)A_{d}k^{d+1}}. \end{equation}
Hence, for 
$S_{3}^{(d)}(r)$ we have the following
\begin{equation} S_{3}^{(d)}(r) = 6\frac{r_{\alpha} r_{\beta} r_{\gamma}}{r^{3}} \int \mathrm{d} \boldsymbol{k} \,\frac{P_{\alpha \beta \gamma}( \boldsymbol{k}) T(k)}{\textrm{i}(d-1)A_{d}k^{d+1}} \textrm{e}^{\textrm{i} \boldsymbol{k}\cdot \boldsymbol{r}}. \end{equation}
In the same way as we did for 
$S^{(d)}_2(r)$ we can evaluate this integral by taking 
$\theta$ to be the angle between 
$\boldsymbol {k}$ and 
$\boldsymbol {r}$, which upon doing so we find
\begin{equation} \frac{r_{\alpha} r_{\beta}r_{\gamma}}{r^{3}}P_{\alpha \beta \gamma}( \boldsymbol{k}) = 2k\cos \theta \left(1-\cos^{2} \theta\right). \end{equation}
We then once more evaluate all but two of the 
$d$ Fourier integrals to obtain
\begin{equation} S_{3}^{(d)}(r) = \frac{12 A_{d-1}}{\textrm{i}(d-1)A_d}\int_{0}^{\infty} \mathrm{d}k \, \frac{T(k)}{k}\int_{0}^{\rm \pi} \mathrm{d}\theta \cos \theta \left(1-\cos^{2} \theta\right)\sin^{d-2}\theta \exp({\textrm{i}kr\cos \theta}). \end{equation}The inner integral is formidable, however, it can be evaluated using computer algebra software. The result when restricted to integer dimensions is then found to be
\begin{equation} S_{3}^{(d)}(r) = 3\varGamma\left(\frac{d}{2}\right)r\int_{0}^{\infty} \mathrm{d}k \, 2^{1+{d}/{2}}T(k)\frac{J_{1+{d}/{2}}(kr)}{\left(kr\right)^{1+{d}/{2}}}. \end{equation} As a check, we compare this result with the case for 
$d=2$ derived in Cerbus & Chakraborty (Reference Cerbus and Chakraborty2017)
\begin{equation} S_{3}^{(2)}(r) = 12 r\int_{0}^{\infty} \mathrm{d}k \, T(k) \frac{J_2 (kr)}{\left(kr\right)^{2}}. \end{equation}
Now, clearly upon inserting 
$d=2$ into (2.30) we recover the result above. Furthermore using properties of Bessel functions it can also easily be shown that for 
$d=3$ the above reduces to the expected expression as seen in Bos et al. (Reference Bos, Chevillard, Scott and Rubinstein2012).
Finally, we consider a small 
$r$ expansion of the third-order structure function in 
$d$ dimensions
\begin{equation} S_{3}^{(d)} = \frac{12r}{d(d+2)}\int_{0}^{\infty}\mathrm{d}k \, T(k) - \frac{6r^{3}}{d(d+2)(d+4)} \int_{0}^{\infty}\mathrm{d}k \, k^{2} T(k) + {O}(r^{5}). \end{equation}
From this expansion and the conservation properties of 
$T(k)$, we can see 
$S_{3}^{(d)}(r) \sim r^{3}$ for small 
$r$ and 
$d\geq 3$. In two dimensions the second term also vanishes as a result of enstrophy conservation. This expansion is also of practical use for evaluation of 
$S_{3}^{(d)}(r)$ for very small 
$r$ numerically, where floating point arithmetic errors can arise
2.3. Enstrophy production and skewness
Vorticity and enstrophy play an important role in the behaviour of two- and three-dimensional turbulence. Enstrophy production is also known to be linked to the velocity derivative skewness, hereafter referred to simply as the skewness, of the flow. To generalise these concepts to arbitrary spatial dimension, we first introduce the vorticity 2-form
\begin{equation} \varOmega_{\alpha \beta}( \boldsymbol{x}) = \partial_{\alpha} u_{\beta}( \boldsymbol{x}) - \partial_{\beta} u_{\alpha}( \boldsymbol{x}). \end{equation}
We can then re-express (2.1) using 
$\varOmega _{\alpha \beta }$ as
\begin{equation} \partial_t u_{\alpha} + \varOmega_{\beta \alpha} u_{\beta} ={-}\partial_{\alpha} \left(P + \frac{u^{2}}{2}\right) + \nu \nabla^{2} u_{\alpha}, \end{equation}
which is valid in any dimension, and is equivalent to the rotational form of the NSE in three dimensions. Using (2.34) we are then able to derive an evolution equation for 
$\varOmega _{\alpha \beta }( \boldsymbol {x})$
\begin{equation} \partial_t \varOmega_{\alpha \beta} + u_{\gamma} \partial_{\gamma} \varOmega_{\alpha \beta} + \varOmega_{\alpha \gamma}S_{\gamma \beta} + \varOmega_{\gamma \beta}S_{\alpha \gamma} = \nu \nabla^{2} \varOmega_{\alpha \beta}, \end{equation}
where 
$S_{\alpha \beta }( \boldsymbol {x}) = (\partial _{\alpha } u_{\beta }( \boldsymbol {x}) + \partial _{\beta } u_{\alpha }( \boldsymbol {x}))/2$ is the strain tensor. Enstrophy in three dimensions is defined in terms of the vorticity, 
$\boldsymbol {\omega }( \boldsymbol {x})$, as
\begin{equation} Z(t) = \tfrac{1}{2}\langle\omega_{\alpha} \omega_{\alpha} \rangle = \int_{0}^{\infty} \mathrm{d}k \, k^{2} E(k), \end{equation}
where for this case we also have 
$\omega _{\alpha } = \epsilon _{\alpha \beta \gamma }\varOmega _{\beta \gamma }/2$, which suggests the correct form of enstrophy in terms of the 2-form is
\begin{equation} Z(t) = \tfrac{1}{4}\langle \varOmega_{\alpha \beta}^{2} \rangle = \int_{0}^{\infty}\mathrm{d}k \, k^{2} E(k). \end{equation}
To be confident this second equality holds, we will first form an equation for the evolution of 
$u^{2}$ using the vorticity 2-form
\begin{equation} \frac{1}{2}\partial_t u_{\alpha}u_{\alpha} + u_{\alpha}\varOmega_{\beta \alpha} u_{\beta} ={-}\partial_{\alpha}\left(P + \frac{u^{2}}{2}\right)u_{\alpha} + \nu \partial_{\beta}u_{\alpha}\varOmega_{\beta \alpha} - \frac{\nu}{2}\varOmega_{\alpha \beta}^{2}. \end{equation}Upon averaging and invoking homogeneity we find
\begin{equation} \partial_t E(t) ={-}\frac{\nu}{2}\langle \varOmega_{\alpha \beta}^{2}\rangle, \end{equation}which, when compared to the standard result, is
\begin{equation} \langle \varOmega_{\alpha \beta}^{2}\rangle = 2\langle\left(\partial_{\beta}u_{\alpha}\right)^{2} \rangle. \end{equation}From here, it can be shown that indeed the second equality in (2.37) holds as
\begin{equation} \langle\left(\partial_{\beta}u_{\alpha}\right)^{2} \rangle ={-}\langle u_{\alpha}\nabla^{2} u_{\alpha} \rangle = 2\int_{0}^{\infty}\mathrm{d}k \, k^{2} E(k), \end{equation}and the equality is proved. Thus, we are confident (2.37) is a consistent generalisation of enstrophy to all dimensions.
In order to relate the production of enstrophy to skewness, we require an equation for our generalised enstrophy, which we can obtain from the 
$d$-dimensional vorticity equation above
\begin{equation} \partial_t Z(t) ={-}\tfrac{1}{2}\langle \varOmega_{ij}\varOmega_{ik}S_{kj} + \varOmega_{ij}\varOmega_{kj}S_{ik} \rangle - \frac{\nu}{2}\langle\varOmega_{ij}\nabla^{2} \varOmega_{ij}\rangle. \end{equation}Following steps similar to those for the enstrophy, we can show that the palinstrophy generalises as
\begin{equation} P(t) =\tfrac{1}{4}\langle\varOmega_{ij}\nabla^{2} \varOmega_{ij}\rangle = \int_{0}^{\infty}\mathrm{d}k \, k^{4} E(k). \end{equation}
Now, in order to express (2.42) in terms of the skewness, we consider the von Kármán–Howarth equation (von Kármán & Howarth Reference von Kármán and Howarth1938) in 
$d$-dimensions expressed in terms of the second- and third-order two-point longitudinal correlations
\begin{equation} \partial_t C_{LL} = \frac{1}{r^{d+1}}\partial_r \left[r^{d+1}C_{LL,L} \right] + \frac{2\nu}{r^{d+1}}\partial_r\left[r^{d+1}\partial_r C_{LL}\right]. \end{equation}From the preceding discussion we now recognise the integrals in (2.18) as being the total energy, enstrophy and palinstrophy, hence we have
\begin{equation} C_{LL} = \frac{2}{d}E(t) -\frac{r^{2}}{d(d+2)}Z(t) + \frac{r^{4}}{4d(d+2)(d+4)}P(t) + {O}(r^{5}). \end{equation}
Using this expansion in (2.44), and the fact that 
$C_{LL,L} \sim r^{3}$ for small 
$r$, we find to zeroth order in 
$r$
\begin{equation} \partial_t E(t) ={-} 2\nu Z(t). \end{equation}
This is entirely equivalent to (2.39) and represents the decay of energy in turbulent flows without external forcing. Continuing now to 
${O}(r^{2})$
\begin{equation} \partial_t Z(t) ={-}\left. d(d+2)(d+4)\frac{C_{LL,L}}{r} \right|_{r\rightarrow 0} -2 \nu P(t). \end{equation}
Both these expressions are consistent with what is derived directly from the NSE (see Davidson (Reference Davidson2015) for the two- and three-dimensional cases). Also, we note here that, since this derivation required 
$C_{LL,L}\sim r^{3}$ for small 
$r$, in two dimensions the first term on the right hand side vanishes. Recalling that 
$S_{3}^{(d)}(r) = 6C_{LL,L}$ and the skewness, 
$S_0$, can be expressed as
\begin{equation} S_0 = \left.\frac{S_{3}^{(d)}(r)}{\left[S_{2}^{(d)}(r)\right]^{{3}/{2}}}\right|_{r\rightarrow 0}, \end{equation}we then write the enstrophy equation as
\begin{equation} \partial_t Z(t) ={-}S_{0}\varLambda(d)Z(t)^{{3}/{2}} -2 \nu P(t), \end{equation}
where the 
${O}(r^{2})$ term of the small 
$r$ expansion of 
$S_{2}^{(d)}(r)$ has been used and the dimensional factor is
\begin{equation} \varLambda(d) = \frac{(d+4)}{3}\sqrt{\frac{2}{d(d+2)}}. \end{equation}
The dimensional factor, 
$\varLambda (d)$, in this equation is a decreasing function of 
$d$ with an asymptotic limit of 
$\sqrt {2}/3$. As such, with increasing dimension a larger skewness is required to generate the same amount of enstrophy.
In the above, we have demonstrated that, in any dimension, enstrophy production is governed by the action of the strain field on the generalised vorticity. This action can be thought of as the stretching and folding of structures analogous to vortices in all dimensions. From another viewpoint, this stretching is seen to be caused by a non-zero skewness. However, a larger skewness value is required with increasing dimension to produce the same level of vortex stretching.
We can also consider the skewness as a function of dimension using (2.48) and the small 
$r$ expansions for the structure functions. This gives
\begin{equation} S_{0}(d) ={-}\frac{1}{\varLambda(d)}\int_{0}^{\infty}\mathrm{d}k \, k^{2} T(k) \left[\int_{0}^{\infty}\mathrm{d}k \, k^{2} E(k)\right]^{-({3}/{2})}. \end{equation}
This result will be used in § 4 in order to measure the dimensional dependence of 
$S_0$.
3. Methods
As was highlighted earlier in this work, the cost of DNS in dimensions higher than four is prohibitive. Therefore, in order to study the effects of spatial dimension on turbulent fluid flows, we utilise the 
$d$-dimensional EDQNM closure approximation. Under this approximation we find the equation for the time evolution of the energy spectrum is
\begin{align} \left(\partial_t + 2\nu k^{2}\right) E(k) &= 8K_d \iint_{{\rm \Delta} k} \mathrm{d}p \, \mathrm{d}q \, \frac{k}{pq}b^{(d)}_{kpq}\theta_{kpq}(t) \nonumber\\ &\quad \times \bigg[\sin^{d-3}(\alpha) k^{2} E(p)E(q) - \sin^{d-3}(\beta) p^{2} E(q)E(k)\bigg] + f(k). \end{align}
It should be noted that the 
$kpq$ subscripts are simply labels and not indices, so no summation is implied. For a derivation of (3.1), see Lesieur (Reference Lesieur1987) or Sagaut & Cambon (Reference Sagaut and Cambon2008) for the three-dimensional case, and Rose & Sulem (Reference Rose and Sulem1978) for the extension to arbitrary dimension. In (3.1) the integration is performed over wave-vector triads which can form triangles, i.e. 
$\lbrace \boldsymbol {k}, \boldsymbol {p}, \boldsymbol {q} \rbrace$ satisfying 
$\boldsymbol {k}+ \boldsymbol {p}+ \boldsymbol {q}= \boldsymbol {0}$, and we have
\begin{equation} b^{(d)}_{kpq} = \frac{p}{2k}\left((d-3)Z + (d-1)XY + 2Z^{3}\right), \end{equation}
where 
$X, Y$ and 
$Z$ are the cosines of the angles 
$\alpha , \beta$ and 
$\gamma$ which lie opposite the sides 
$k,p$ and 
$q$ respectively. As such, 
$b^{(d)}_{kpq}$, coupled with the two sine terms, contains all the information regarding the geometry of the triadic interactions in the EDQNM closure. We also have the triad relaxation time
\begin{equation} \theta_{kpq}(t) = \frac{1 - \exp({-(\mu_{k} +\mu_{p} + \mu_{q} + \nu(k^{2} + p^{2} + q^{2}))t})}{\mu_{k} + \mu_{p} + \mu_{q} + \nu(k^{2} + p^{2} + q^{2})}, \end{equation}in which the eddy-damping rate is given by
\begin{equation} \mu_{k} = \lambda_1 \sqrt{\int_{0}^{k} \mathrm{d}s \, s^{2} E(s,t)}. \end{equation}
Finally the dimensional factor 
$K_d$ is given by
\begin{equation} K_d = \frac{A_{d-1}}{(d-1)^{2} A_d}. \end{equation}
In these equations there exists a free parameter 
$\lambda _1$, which can be shown to set the Kolmogorov constant (André & Lesieur Reference André and Lesieur1977), see appendix A for further details. For certain applications setting the value of the Kolmogorov constant will be important, for example in investigating the dimensionless dissipation rate, whilst for others we will be more interested in the triadic interactions which are independent of the choice of 
$\lambda _1$. The forcing term 
$f(k)$ is chosen such that injected energy is distributed evenly across the forcing range of wavenumbers. This forcing allows for the energy dissipation, 
$\varepsilon$, to be set exactly at stationary state. As most of this work will be focused on quantities that are determined by the small scales of the flow, the exact form of forcing at the large scales is not an important choice. Moreover, we have made use of different forms of forcing and verified that the small-scale behaviour is unaltered.
A major advantage of this closure is that, since the energy spectrum is a smooth function of 
$k$, we can make use of a logarithmic discretisation of the wavenumber space when solving numerically (Leith Reference Leith1971). This is the key feature that allows very high Reynolds numbers to be achieved. Throughout this work we discretise our numerical simulations with
\begin{equation} k_i = k_{{min}} 2^{i/F} \end{equation}
and we choose 
$k_{{min}} = 1$ and vary 
$F$ to maintain resolution. In particular, we choose 
$F$ such that the choice of 
$\lambda _1$ accurately determines the value of the Kolmogorov constant in the inertial range, see appendices A and B. Due to the presence of oscillatory terms raised to the power of the spatial dimension in the EDQNM equation, a finer mesh, and thus higher value for 
$F$, is required. Much consideration was given to this point by Leith & Kraichnan (Reference Leith and Kraichnan1972) using the test field model. In this case they performed an averaging procedure on the oscillatory terms, using a geometric approach as opposed to the elegant analytic approach used by Bowman (Reference Bowman1996). In order to simplify our approach we choose simply to use a finer mesh with increasing dimension rather than invoke an averaging procedure. If the simulation is under-resolved the value of the Kolmogorov constant will not be correct, as such, measurement of this constant provides a useful check for the accuracy of a simulation, we thus ensure that the value of the Kolmogorov constant measured in our simulations is in agreement with the results of appendix B.
It is important to note that it has been shown (Lesieur & Schertzer Reference Lesieur and Schertzer1978) that, when using such a logarithmic discretisation, interactions between certain non-local triads will not be accounted for in the EDQNM closure. It is possible to handle these missing triads analytically, however, for this study the extension of such methods to 
$d$-dimensions has not been carried out. Consequently, we have verified our main results are not influenced by the choice of discretisation through cross-checks with linearly discretised simulations where feasible. In our numerical work, we employ a parallel EDQNM code using second order predictor corrector method for time stepping. Details of this code can be found in Clark (Reference Clark2019).
In order to determine an appropriate value for 
$\lambda _1$ in any dimension, we turn to results from a parameterless closure approximation. The Lagrangian renormalised approximation (LRA) (Kaneda Reference Kaneda1981) was used by Gotoh et al. (Reference Gotoh, Watanabe, Shiga, Nakano and Suzuki2007) to obtain values for the Kolmogorov constant in arbitrary dimension. Using these values we are then able to set the free parameter in our EDQNM simulations, such that we have an appropriate value for the Kolmogorov constant, see appendix B for further details. In our numerical work above four dimensions 
$\lambda _1$ has been set such that the Kolmogorov constant in each dimension is equal to the value predicted by the LRA. Unfortunately, this ad hoc method is all that is available as there are no DNS results beyond four dimensions from which a Kolmogorov constant can be approximated.
In all simulations we start from a zero initial energy spectrum. Simulations are then evolved from this spectrum until a statistically stationary state is reached. The exact form of the initial condition is not important for this work: the statistically steady state is independent of such choices. In fact, as we begin from a zero state the steady state spectra is guaranteed to be generated entirely by the triadic interactions of (3.1). In appendix B we give further details on the simulations performed included the largest Reynolds numbers reached.
4. Results
4.1. Energy and transfer spectra
In the DNS study of Berera et al. (Reference Berera, Ho and Clark2020), the observed scaling of the energy spectrum was consistent with that predicted by K41 in both the three- and four-dimensional cases insofar as within a given dimension the energy spectra were found to collapse upon scaling by Kolmogorov variables. However, the spatial dimension was found to have an influence on the shape of the energy spectrum in the near dissipative region. It was observed that, when comparing four to three dimensions, dissipative effects did not become dominant until smaller scales, evidenced by the presence of a seemingly extended inertial range. Before we investigate this behaviour in higher dimensions using the EDQNM closure, we need to understand to what extent the model is capable of reproducing the effects seen in three- and four-dimensional DNS. To this end, in figure 1 we show energy spectra from both DNS and EDQNM in three and four dimensions scaled by the Kolmogorov constant. In both dimensions we have the same 
$\nu$ and 
$\varepsilon$ across DNS and EDQNM. A good collapse of the data can be seen in the inertial range and, approaching the dissipative region, the EDQNM closure captures the extended inertial range reasonably well. However, once the dissipative region is reached, the EDQNM results begin to diverge from those of the DNS in both dimensions, and it does not appear that this divergence is worse in one dimension over the other. The discrepancies at the large scales, small 
$k\eta$, are due to the forcing differences between DNS and EDQNM simulations.
Figure 1. (a) Comparison between DNS (solid lines) and EDQNM (dashed lines) of energy spectra scaled by the Kolmogorov microscale, 
$\eta$, for three (blue) and four (black) dimensions. (b) Compensated energy spectra for the same data.
Having verified that the EDQNM approximation can satisfactorily reproduce properties of the energy spectra seen in three- and four-dimensional DNS, we now turn to even higher spatial dimensions. In figure 2(a) we plot the energy spectra from EDQNM simulations for three, six, seven, ten and twenty dimensions, scaled by the appropriate Kolmogorov constant for each dimension, once more we keep 
$\nu$ and 
$\varepsilon$ constant across dimensions. Here, it can be seen that increasing spatial dimension is accompanied by a growing accumulation of energy on the edge of the inertial range. Such behaviour may be indicative of an enhanced forward transfer of energy as the spatial dimension increases. This view is consistent with theoretical arguments which conjecture that as the spatial dimension tends to infinity the nature of the triadic interactions leads to all energy being transferred in the forward direction to the small scales (Fournier et al. Reference Fournier, Frisch and Rose1978). Further arguments have suggested that the appearance of such an energy bottleneck is the result of triad interactions being damped by viscosity at the smallest scales (Herring et al. Reference Herring, Schertzer, Lesieur, Newman, Chollet and Larcheveque1982). Hence, if the forward energy transfer is enhanced and, as appears to be the case, the aforementioned viscous damping of certain triads is not influenced by the dimension, then this pile up should be expected to increase with dimension. In figure 2(b), we show the compensated spectra which give an even clearer demonstration that the bottleneck effect becomes more pronounced with increasing dimension. In DNS results (Berera et al. Reference Berera, Ho and Clark2020), this view is further supported by an increased skewness in four dimensions. We will return to skewness in § 4.2 where it is evaluated for EDQNM results. It is also clear that in all cases we observe an inertial range with a 
$k^{-5/3}$ power law scaling which persists over a number of decades in wavenumber space. As we go to higher dimensions we find this scaling region appears to become progressively shortened by the increased bottleneck effect. Although not shown, in all dimensions we find a collapse within a given dimension of energy spectra across a range of 
${{\textit {Re}}}$ values when rescaled by 
$\nu$ and 
$\varepsilon$. That is, the energy spectrum is found to take on a universal shape in each dimension.
Figure 2. (a) EDQNM energy spectra scaled by the Kolmogorov microscale, 
$\eta$, for three, six, seven, ten and twenty dimensions, the darker the shade of the line the higher the dimension with twenty being black. Dashed line shows 
$k^{-5/3}$ scaling. All dimensions have the same viscosity and energy dissipation. (b) Compensated spectra for the same data, colours the same.
The energy transfer spectrum, 
$T(k)$, provides a natural measure of the exchange of energy between different scales in turbulent flows. However, in closing the infinite moment hierarchy using the quasi-normal hypothesis, the transfer spectrum is directly affected by the closure. Moreover, 
$T(k)$ is also directly influenced by the eddy-damping assumption. As such, if the quasi-normal and/or the eddy-damping assumptions are not sound we should expect the transfer spectra produced in the EDQNM simulations to show differences when compared to DNS transfer spectra. Indeed, in figure 3(a) we see far larger discrepancies between EDQNM and DNS for the nonlinear energy transfer than we did in the corresponding energy spectra. However, the qualitative behaviour going from three to four dimensions is the same in both DNS and EDQNM results. The peak nonlinear transfer is greater and found at smaller scales in four dimensions compared to three dimensions. However, in contrast with results for the energy spectra, the agreement between the nonlinear transfer in DNS and EDQNM appears to be better in four dimensions than in three dimensions, insofar as the peak transfer occurs at similar scales in both DNS and EDQNM. Without higher-dimensional DNS results we cannot say whether this better agreement between DNS and EDQNM is purely coincidental or in fact evidence that four-dimensional turbulence is in effect more mean-field like. In figure 3(b), the energy flux is displayed. For both DNS and EDQNM the energy flux remains roughly constant until smaller scales in four dimensions relative to three dimensions. Once more, due to the forcing differences at larger scales there is a greater disagreement between DNS and EDQNM results for both the nonlinear transfer and the energy flux. In light of these comparisons, we should be more cautious in our interpretation of results derived from the transfer spectrum in EDQNM.
Figure 3. (a) Comparison between DNS (solid lines) and EDQNM (dashed lines) of nonlinear transfer in three (blue) and four (black) dimensions. (b) Normalised energy flux for DNS (solid lines) and EDQNM (dashed lines) in three (blue) and four (black) dimensions.
Turning once more to purely EDQNM results, we look at the dependence on the spatial dimension of the nonlinear energy transfer. In figure 4(a) we plot the nonlinear energy transfer for a number of dimensions. Here, the trend of the peak nonlinear energy transfer moving to smaller scales as the dimension increases is observed to continue to around six dimensions, at which point it begins to move to larger scales again. This potential cross-over at six dimensions is interesting given the work of Liao (Reference Liao1990, Reference Liao1991), where a possible upper critical dimension for turbulence at six dimensions is conjectured. Of course, since our results are obtained via a closure approximation they should not be over-interpreted. In figure 4(b) we show the spectral energy flux for a range of spatial dimensions. In line with what was found for the compensated energy spectra, we observe a scaling range of several decades in all dimensions. In this figure it is clear that above three dimensions there is an increased energy transfer to smaller scales as evidenced by the flux dropping slower as we enter the dissipative region.
Figure 4. (a) EDQNM nonlinear energy transfer for three, six, seven, ten and twenty dimensions, the darker the shade of the line the higher the dimension with twenty being black. (b) Normalised energy flux for same data coloured in the same manner.
4.2. Skewness
In § 2.3 we derived (2.49), relating the production of the generalised enstrophy in 
$d$-dimensions to the skewness, 
$S_0$. Higher values of 
$S_0$ are then associated with greater vortex stretching, which would then provide the mechanism for the increased forward energy transfer in higher dimensions. Using (2.51) we can measure the effect of spatial dimension on 
$S_0$ in our EDQNM simulations. These results are presented in figure 5(a). It is observed that in the EDQNM equations the skewness reaches a maximum value of approximately 
$-$0.72 at seven dimensions, before remaining roughly constant until ten dimensions, beyond which 
$S_{0}$ decreases. If the trend seen in figure 5(a) continues then the skewness may vanish for infinite spatial dimension. When compared to the results for skewness in the DNS of Berera et al. (Reference Berera, Ho and Clark2020) both three- and four-dimensional EDQNM results exhibit a lower value of 
$S_0$. This is not a surprising result and has been observed in EDQNM simulations in three dimensions (Bos et al. Reference Bos, Chevillard, Scott and Rubinstein2012) and is likely a result of the assumptions made in the EDQNM model e.g. the quasi-normal or eddy-damping assumptions. The EDQNM approximation is also known to exhibit a constant asymptotic value for 
$S_0$ at sufficiently high Reynolds number and the results presented in figures 5(a) and 5(b) are of this asymptotic value in all cases. This existence of this asymptotic value is predicated on the exponent in the inertial range being 
$-5/3$ hence given our energy spectrum results in all dimensions we can be confident in these asymptotic values.
Figure 5. (a) EDQNM velocity derivative skewness vs dimension. (b) Enstrophy production term 
$-S_{0} \varLambda (d)$ vs dimension.
In (2.49), due to the dimensional pre-factor of the skewness term, at larger spatial dimensions higher skewness values are required to generate the same level of enstrophy production. In figure 5(b) we show the pre-factor, 
$-S_{0}\varLambda (d)$, of the enstrophy production term in (2.49) for a range of dimensions. These results suggest that until five dimensions there is an increase in enstrophy production, then starting at six dimensions there is a reduction in enstrophy production and thus in vortex stretching. If we consider that the action of vortex stretching produces smaller scales in the flow, then this is consistent with the accumulation of energy at the end of the inertial range. Once more, if this trend continues, then for infinite dimension there will be no enstrophy production and hence no vortex stretching. The vanishing of vortex stretching, and thus of velocity derivative skewness, at infinite dimension is quite an extreme scenario, and it may be that a non-zero but finite asymptotic value is reached instead. This is consistent with what was seen in Gotoh et al. (Reference Gotoh, Watanabe, Shiga, Nakano and Suzuki2007) where a finite but non-zero asymptotic skewness value at infinite dimension was predicted using the LRA. This finite value is reached after a consistent increase with dimension. This is at odds with what is seen in our simulations, where a maximum skewness is reached at a finite dimension. However, if the skewness does vanish at infinite dimension a more drastic statement may then be that in infinite dimensions there is no energy flux in either direction and thus no turbulence. We stress here that since these are closure results it is not possible to make any definitive claims.
The significance of this maximum skewness dimension is unknown, and may be a result of the closure assumptions. As such, without higher-dimensional DNS results, we cannot make a definitive statement. Furthermore, to ensure that the skewness reaching a maximum at finite dimension is not a result of the values of free parameter, 
$\lambda _1$, used in this work we have carried out an additional set of simulations detailed in appendix C. As a result of these simulations we are confident the appearance of a maximum skewness dimension is a real effect in the EDQNM equations, with the caveat that the Kolmogorov constant does not increase with dimension. Given the observed Kolmogorov constants in the four-dimensional DNS of both Gotoh et al. (Reference Gotoh, Watanabe, Shiga, Nakano and Suzuki2007) and Berera et al. (Reference Berera, Ho and Clark2020) and the predictions of the LRA there is evidence that this caveat holds.
An important further point to consider for the infinite-dimensional problem is that of the energy spectrum normalisation in this limit. In Fournier et al. (Reference Fournier, Frisch and Rose1978) it was found that, to each order in perturbation theory, the energy spectrum had a finite limit if a rescaled time variable is used. However, this time rescaling can be shown to be equivalent to a rescaling of the energy spectrum. The latter rescaling is equivalent to a finite energy per velocity component as opposed to a finite total energy. If the total energy is to remain finite then as we tend to infinite dimension each individual velocity component will tend to zero, suggesting zero skewness. However, if the components remain non-zero then the skewness will be finite but non-zero.
Once more, these results are interesting considering the work of Liao (Reference Liao1990, Reference Liao1991) suggesting the possibility of a critical dimension of 6 for turbulence. Studying this infinite-dimensional limit within numerical EDQNM studies becomes difficult with increasing dimensions. For example, at 
$d=20$ with a numerical resolution of 
$F=50$ we find an error of approximately 2 % in the expected value for the Kolmogorov constant using the free parameter in appendix B. Compared with 
$d=3$ where a resolution of 
$F=16$ gives an error much less than a per cent it is clear that the resolution costs of higher dimensions quickly becomes an issue.
4.3. Third-order structure functions
The longitudinal structure functions have frequently been measured in experimental (Anselmet et al. Reference Anselmet, Gagne, Hopfinger and Antonia1984) and numerical studies of turbulence (Gotoh, Fukayama & Nakano Reference Gotoh, Fukayama and Nakano2002; Ishihara, Gotoh & Kaneda Reference Ishihara, Gotoh and Kaneda2009). In such studies, it is found that, particularly at higher orders, these structure functions show deviations from the scaling predicted by the K41 theory. Such deviations are typically attributed to intermittency corrections, although there are also arguments suggesting these are simply finite Reynolds number corrections due to K41 being an asymptotic theory (Kolmogorov Reference Kolmogorov1962; Benzi et al. Reference Benzi, Paladin, Parisi and Vulpiani1984; Qian Reference Qian1997, Reference Qian1999; Antonia & Burattini Reference Antonia and Burattini2006; Tchoufag, Sagaut & Cambon Reference Tchoufag, Sagaut and Cambon2012; McComb Reference McComb2014; Tang et al. Reference Tang, Antonia, Djenidi, Danaila and Zhou2017; Antonia et al. Reference Antonia, Tang, Djenidi and Zhou2019).
The effect of the spatial dimension on these corrections is an interesting question, and comparisons to critical phenomena, in particular anomalous exponents, have been made (Nelkin Reference Nelkin1974; Rose & Sulem Reference Rose and Sulem1978). Of course, the EDQNM model does not account for intermittency and therefore measurements of structure functions in such simulations cannot answer questions regarding these deviations. Furthermore, we are not aware of a method by which to calculate beyond third-order structure functions from spectral quantities. However, if such a method were available it is likely it could, potentially with some difficulty, be applied to the EDQNM model. Indeed, it was shown by Kraichnan (Reference Kraichnan1959) that in the direct interaction approximation (DIA) it is possible to study quantities beyond third order. This was then demonstrated by Chen et al. (Reference Chen, Herring, Kerr and Kraichnan1989) for the DIA and by Bos & Rubinstein (Reference Bos and Rubinstein2013) for the EDQNM model. As such, even if intermittency was present in the EDQNM model we would be unlikely to find substantial corrections at these low orders, even in three dimensions.
Putting questions regarding intermittency and anomalous exponents aside, we wish to test (2.30) using our DNS and EDQNM results. In figures 6(a) and 6(b), the third-order structure functions computed using (2.30) in EDQNM are compared to those in DNS. Good agreement is seen in the inertial range with both DNS and EDQNM exhibiting the expected scaling. Looking to the dissipative region, we find that all our data follow 
$r^{3}$ scaling, however, this scaling begins at a different point in DNS compared to EDQNM. Given our energy and transfer spectrum results, these deviations in small-scale behaviour are expected. The most interesting feature of these figures can be found in figure 6b) where the agreement between DNS and EDQNM appears to be better in four dimensions when compared with three. This is consistent with what was seen in the nonlinear transfer, perhaps not surprisingly given that the transfer spectrum is used in determining 
$S^{(d)}_{3}(r)$. Without higher-dimensional DNS data it is impossible to know if this better agreement is due to the EDQNM approximation becoming more accurate with increasing dimension or simply a coincidence.
Figure 6. Comparison of third-order structure functions in DNS (solid lines) and EDQNM (short dashed lines) for three (blue) and four (black) dimensions: (a) scaled by appropriate Kolmogorov quantities; (b) scaled by energy dissipation and 
$r$. In (a) the dotted lines correspond to the expected 
$12r/d(d+2)$ inertial range scaling behaviours whilst in (b) they are for the values 
$12/d(d+2)$.
As already stated, the EDQNM model does not exhibit intermittency. However, it does capture well the finite Reynolds number effect. Indeed, in Bos et al. (Reference Bos, Chevillard, Scott and Rubinstein2012) the EDQNM model was compared with the multi-fractal model for three-dimensional turbulence giving comparable results for low 
${{\textit {Re}}}$ suggesting it is difficult to distinguish between intermittency effects and finite Reynolds effects in this region. We have not studied the multi-fractal model in four dimensions to make this comparison, however, given the better agreement between EDQNM and DNS in four dimensions for 
$S^{(d)}_{3}(r)$ it may be the case that the finite Reynolds number effect becomes dominant over intermittency in higher dimensions.
Looking now at figure 7(a), we see the third-order longitudinal structure function for our EDQNM data in four, five and six dimensions scaled by appropriate Kolmogorov quantities. It can be seen that in all dimensions we have 
$r^{3}$ scaling in the small 
$r$ limit, as is seen in three-dimensional turbulence and is predicted from the small 
$r$ expansion of 
$S^{(d)}_{3}(r)$. Turning to the 
$d$-dimensional von Kármán–Howarth equation, it can be shown that in the inertial range for the third-order structure function, we should find
\begin{equation} S_{3}^{(d)}(r) \simeq{-}\frac{12}{d(d+2)}\varepsilon r, \end{equation}
which reduces to the standard four-fifths law of three-dimensional turbulence. Indeed, in figure 7(a) we can see that each dimension follows its own 
$-12/d(d+2)$ law in the inertial range. Once more in all dimensions we observe a long scaling region. For a clearer comparison in figure 7(b) we normalise each dimension by the expect inertial range value such that all cases show an scaling range at 1. In doing so we find differences across dimensions, in particular by 
$d=20$ the scaling region begins at higher 
$r/\eta$ than in lower dimensions.
Figure 7. EDQNM third-order structure functions scaled by appropriate Kolmogorov quantities for three, six, seven, ten and twenty dimensions, the darker the shade of the line the higher the dimension with twenty being black. Dashed lines represent appropriate power law scaling for each dimension, i.e. 
$12r/d(d+2)$ for the inertial range. (b) Third-order structure function scaled by energy dissipation and 
$r$ for the same data with the same colouring.
4.4. Dissipative anomaly
In both experimental (Sreenivasan Reference Sreenivasan1984; Burattini, Lavoie & Antonia Reference Burattini, Lavoie and Antonia2005) and numerical studies (Wang et al. Reference Wang, Chen, Brasseur and Wyngaard1996; Kaneda et al. Reference Kaneda, Ishihara, Yokokawa, Itakura and Uno2003; Ishihara et al. Reference Ishihara, Morishita, Yokokawa, Uno and Kaneda2016) of three-dimensional turbulence, including in EDQNM (Bos, Shao & Bertoglio Reference Bos, Shao and Bertoglio2007), there is a large body of evidence which indicates the existence of a non-zero energy dissipation rate, even in the limit of zero viscosity. This is known as the dissipative anomaly. In Berera et al. (Reference Berera, Ho and Clark2020) an increased value for this asymptotic dissipation rate was observed in four dimensions when compared with three. This result once again suggests an enhancement of the forward energy cascade in four dimensions compared with three. We should then expect that beyond four dimensions this asymptotic dissipation rate should increase further given our spectral and skewness results.
The dimensionless dissipation rate is defined as
\begin{equation} C_{\varepsilon} = \frac{\varepsilon L}{u^{3}}, \end{equation}and its Reynolds number dependence can be shown to be approximately described by the relationship (Doering & Foias Reference Doering and Foias2002; McComb et al. Reference McComb, Berera, Yoffe and Linkmann2015)
\begin{equation} C_{\varepsilon}({{\textit{Re}}}) = C_{\varepsilon,\infty} + \frac{C}{{{\textit{Re}}}}, \end{equation}
where 
${{\textit {Re}}}=uL/\nu$ is the integral-scale Reynolds number in which 
$L$ is defined by (2.16).
In figure 8(a) we show the dimensionless dissipation rate against 
${{\textit {Re}}}$ for a wide range of 
${{\textit {Re}}}$ values for three, four and five dimensions. In all cases we observe 
$C_{\varepsilon }$ tending to a constant asymptotic value. This asymptotic value is seen to grow with dimension as would be expected from increased forward energy transfer. We find that in three dimensions this asymptotic value in EDQNM is 
$C^{3d}_{\varepsilon ,\infty } = 0.38$ lower than what is seen in DNS where 
$C^{3d}_{\varepsilon ,\infty } \approx 0.5$ in the forced case. A higher value was found in the EDQNM of Bos et al. (Reference Bos, Shao and Bertoglio2007), however, they take a different choice of 
$\lambda _1$ which will directly influence the value found. In four dimensions we find 
$C^{4d}_{\varepsilon ,\infty } = 0.96$ lower than the value of 1.26 found in DNS (Berera et al. Reference Berera, Ho and Clark2020). Since 
$C_{\varepsilon }$ is defined in terms of large-scale quantities, where the forcing is active, these discrepancies between DNS and EDQNM are not surprising given the differences at small 
$k$ seen in figure 3. In figure 8(b) we show the asymptotic dimensionless dissipation rate against the spatial dimension. Here we find the asymptotic value grows with dimension until at least 10 dimensions. Beyond here the resolution requirements render further calculation increasingly difficult. As 
$C_{\varepsilon , \infty }$ is related to the Kolmogorov constant such measurements would become exceptionally sensitive to resolution errors.
Figure 8. (a) Plot of 
$C_{\varepsilon }$ vs 
${{\textit {Re}}}$ for three, four and five dimensions, darker colour indicates higher dimension with five dimensions in black. Dashed lines indicate the value for 
$C_{\varepsilon , \infty }$; (b) 
$C_{\varepsilon , \infty }$ against spatial dimension, 
$d$.
5. Conclusion
Motivated by the four-dimensional DNS results presented in Berera et al. (Reference Berera, Ho and Clark2020), as well as analogies to critical phenomena, in this paper we have performed a thorough investigation into the effect of the spatial dimension in the EDQNM model of turbulence. While this is only an approximation to true fluid turbulence, we find that it is able to satisfactorily reproduce many of the dimensional effects seen in DNS. To facilitate this study, a number of standard results from three-dimensional turbulence have been extended to 
$d$-dimensional turbulence. Some of these quantities have been discussed in the literature with theoretical ideas as to how they will behave in higher dimensions. This paper has presented for the first time both numerical results using EDQNM in a range of dimensions above three, and expressions for the second- and third-order structure functions in terms of spectral quantities for any dimension 
$d$. Furthermore, an equation relating the production of enstrophy to skewness in 
$d$-dimensional turbulence was derived from the von Kármán–Howarth equation.
The energy and transfer spectra were measured across a wide range of spatial dimensions. For the three- and four-dimensional cases comparisons were made with DNS results where it was observed that the EDQNM model accurately captures the dimensional differences seen in DNS. In these measurements along with those of the skewness as function of spatial dimension, we have found a consistent picture suggesting the forward energy cascade becomes enhanced as spatial dimension increases. In terms of spectra, this can be seen as an increase in the bottleneck effect in the near dissipation region of the energy spectra, which grows with dimension. This suggests that the viscous damping of triad interactions is either not effected by dimension, or is enhanced at a lesser rate than the energy transfer. Further to this, in the transfer spectra, we see a larger peak in the nonlinear transfer. We also observe that the position of this peak first moves to smaller scales before then reversing and moving back to larger scales, which we posit is a result of suppression of transfer to the dissipative modes as suggested by Herring et al. (Reference Herring, Schertzer, Lesieur, Newman, Chollet and Larcheveque1982). We find that enstrophy production reaches a maximum in five dimensions, possibly going to zero at very high dimension. In light of the enstrophy–skewness equation we have derived, this corresponds to a reduction in small-scale vortex stretching, which is again consistent with the enhanced forward transfer bottleneck effect. The possibility of a zero enstrophy production in the limit of infinite dimensions poses interesting questions for the fate of turbulence in this limit.
Additionally, we have measured the third-order structure functions in higher dimensions using the spectral relations we have derived. We find here that each dimension has its own analogue to the four-fifths law of three-dimensional turbulence. Interestingly, when looking at the third order structure function, the comparison between EDQNM and DNS appears better in four dimensions compared with three suggesting possible changes in the turbulent dynamics. Finally, we studied the effect of the spatial dimension on the asymptotic dissipation rate. Due to being defined in terms of large-scale quantities, this is more difficult to accurately measure in the EDQNM model. However, we do find an increase in this asymptotic dissipation rate with dimension, which is a continuation of the trend seen in the DNS performed in Berera et al. (Reference Berera, Ho and Clark2020). This is also consistent with the existence of an enhanced forward transfer of energy in higher dimensions.
These results are interesting for a number of reasons. Importantly, they confirm many of the results found in four spatial dimensions from DNS (Berera et al. Reference Berera, Ho and Clark2020). The fact the EDQNM results show consistency with DNS in three and four dimensions then suggests the trends found by this method at even higher dimensions, which at present are computationally too demanding for DNS, should have some reliability. Thus, this paper has helped to examine some of the theoretical ideas that have been in the literature for decades on the behaviour of turbulence in dimensions greater than three. Furthermore, the appearance of an increased bottleneck effect in higher dimensions may help to shed light on the standard three-dimensional bottleneck effect. Perhaps counter-intuitively, these results are also of interest for how little turbulence changes from three to twenty dimensions. Indeed, compared with the dramatic changes observed going from three to two dimensions the differences between three and twenty are minimal, which is certainly intriguing.
In Berera et al. (Reference Berera, Ho and Clark2020) the scaling behaviour of fluctuations with Reynolds number was measured and found to decrease in four dimensions compared with three. In this work, temporal fluctuations in the total energy were measured and found to scale slower with 
${{\textit {Re}}}$ in four dimensions. It would be interesting to understand how such fluctuations would scale in even higher dimensions, especially in light of the results in this paper concerning the bottleneck effect and velocity derivative skewness. However, measurement of these fluctuations is unfortunately outside the scope of EDQNM calculations and would require future DNS study.
Funding
This work has used resources from ARCHER (http://www.archer.ac.uk) via the Director's Time budget. D.C. is supported by the University of Edinburgh, R.D.J.G.H is supported by the U.K. Engineering and Physical Sciences Research Council (EP/M506515/1). A.B. acknowledges funding from the U.K. Science and Technology Facilities Council.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Setting the Kolmogorov constant
In order to make use of the EDQNM model, it is necessary to specify the value of the free parameter 
$\lambda _1$ seen in (3.1). The value of this constant can be shown to fix the value of the Kolmogorov constant, 
$C_d$, which is important for a number of numerical measurements. Here, we extend the method used by McComb (Reference McComb1990) to derive of the relationship between 
$\lambda _1$ and 
$C_d$ to the 
$d$-dimensional case. An alternative derivation of the relationship for the three-dimensional case can also be found in André & Lesieur (Reference André and Lesieur1977).
To begin, we consider the eddy-damping rate defined in (3.4) in the limit 
$\nu \rightarrow 0$, such that the energy dissipation rate remains constant and, taking the energy spectrum to be a Kolmogorov spectrum to infinity and find
\begin{equation} \mu_{k} = \frac{\sqrt{3}}{2}\lambda_1 \sqrt{C_d} \varepsilon^{{1}/{2}}k^{{2}/{3}}. \end{equation}
On dimensional grounds, in the inertial range we can take 
$\mu _{k}$ to have the form
\begin{equation} \mu_{k} = \beta_d \varepsilon^{{1}/{2}}k^{{2}/{3}}, \end{equation}and thus we have
\begin{equation} C_{d} = \frac{4\beta^{2}_{d}}{3\lambda^{2}_{1}}. \end{equation}
To make further progress we now need to relate 
$C_{d}$ and 
$\beta$ in the EDQNM model. We begin by considering the forced Lin equation at stationary state
\begin{equation} -T(k) ={-}2\nu k^{2} E(k) + F(k), \end{equation}
where 
$F(k)$ is the forcing spectrum. Now, in taking the limit of zero viscosity, the dissipation range will move to the smallest possible scales, hence we have
\begin{equation} 2\nu k^{2} E(k) \rightarrow\varepsilon\delta(k-\infty), \end{equation}and energy conservation then implies that
\begin{equation} F(k) \rightarrow \varepsilon\delta(k). \end{equation}
Upon integrating both sides of the Lin equation up to a value 
$\kappa$, which by the symmetry properties of the integral is arbitrary so long as it is neither 0 nor 
$\infty$, we find
\begin{equation} -\int_{0}^{\kappa} \textrm{d} k \, T(k) = \varepsilon. \end{equation}
In the EDQNM model from (3.1), we have a closed expression for 
$T(k)$, using which we can perform this integration and relate 
$C_d$ and 
$\lambda _1$. Hence, we have
\begin{align} & 8K_d \int_{0}^{\kappa} \textrm{d} k \iint_{{\rm \Delta} k} \mathrm{d}p \, \mathrm{d}q \, \frac{k}{pq}b^{(d)}_{kpq}\theta_{kpq} \nonumber\\ &\quad \times \bigg[\sin^{d-3}(\beta) p^{2} E(q)E(k)- \sin^{d-3}(\alpha) k^{2} E(p)E(q)\bigg] = \varepsilon. \end{align}We can re-express the above integral as
\begin{align} & -8K_d C^{2}_{d}\int_{0}^{\kappa} \textrm{d} k \int_{\kappa}^{\infty} \textrm{d} p \int_{|k-p|}^{k+p} \textrm{d} q \frac{k}{pq}b^{(d)}_{kpq} \nonumber\\ &\quad \times \frac{\sin^{d-3}(\beta) p^{2} q^{-({5}/{3})}k^{-({5}/{3})} - \sin^{d-3}(\alpha) k^{2} p^{-({5}/{3})}q^{-({5}/{3})} }{\beta_d\left(k^{{2}/{3}} + p^{{2}/{3}} + q^{{2}/{3}}\right)}= 1. \end{align}
The two sine terms can also be expressed in terms of 
$k,p$ and 
$q$ as
\begin{equation} \sin (\alpha) = \sqrt{1 - \left(\frac{p^{2} + q^{2} - k^{2}}{2pq}\right)^{2}} \quad \mathrm{and} \quad \sin (\beta) = \sqrt{1 - \left(\frac{k^{2} + q^{2} - p^{2}}{2kq}\right)^{2}}. \end{equation}
The resulting integral must be evaluated numerically and, if we denote it by 
$I_d$ then we have
\begin{align} & I_d = \int_{0}^{\kappa} \textrm{d} k \, \int_{\kappa}^{\infty} \textrm{d} p \, \int_{|k-p|}^{k+p} \textrm{d} q \, \frac{k}{pq}b^{(d)}_{kpq} \nonumber\\ &\quad \times \frac{\sin^{d-3}(\beta) p^{2} q^{-({5}/{3})}k^{-({5}/{3})} - \sin^{d-3}(\alpha) k^{2} p^{-({5}/{3})}q^{-({5}/{3})} }{\left(k^{{2}/{3}} + p^{{2}/{3}} + q^{{2}/{3}}\right)}, \end{align}and
\begin{equation} \frac{C^{2}_{d} I_d}{\beta_d} = \frac{1}{8K_d}. \end{equation}
Then using (A 3) we can eliminate 
$\beta _d$, after which we find
\begin{equation} C_d = \left(\frac{\sqrt{3}}{16 K_d I_d}\right)^{{2}/{3}}\lambda^{{2}/{3}}_{1}. \end{equation}
This is our desired result and after numerical evaluation of 
$I_d$ this can be used to fix 
$C_d$ in simulation.
To make this more concrete we present the results of this procedure for the cases 
$d=3$ and 
$4$ where we take 
$\kappa =1$. For 
$d=3$ the dimensional factor 
$K_d$, which results from performing spherical integration in 
$d$-dimensions, is 
$K_3 = 1/8$, hence, we look to solve
\begin{equation} C_3 = \left(\frac{\sqrt{3}}{2 I_d}\right)^{{2}/{3}}\lambda^{{2}/{3}}_{1}. \end{equation}
It is found that 
$I_d \approx 0.19038$ and therefore for 
$d=3$ we have
\begin{equation} C_3 \approx 2.75 \lambda^{{2}/{3}}_{1}. \end{equation}
This pre-factor differs very slightly compared with André & Lesieur (Reference André and Lesieur1977), however, this is likely a result of different methods of numerical integration. Following a similar procedure for 
$d=4$ it is found that
\begin{equation} C_4 \approx 2.6 \lambda^{{2}/{3}}_{1}. \end{equation}Appendix B. Simulation parameters
For completeness, here we present the values used for the Kolmogorov constant in each dimension alongside the corresponding value of the free parameter. For the Kolmogorov constant the values were obtained from DNS results for three and four dimensions (Berera et al. Reference Berera, Ho and Clark2020) whilst in higher dimensions results obtained in Gotoh et al. (Reference Gotoh, Watanabe, Shiga, Nakano and Suzuki2007) by using the LRA (Kaneda Reference Kaneda1981) are used. The appropriate value for the free parameter in each dimension is obtained using the method in appendix A. These values are displayed in table 1. Additionally we give both the integral length scale and Taylor length scale Reynolds numbers for the largest cases in each dimension.
Table 1. Simulation parameters for all dimensions. In all cases we have 
$\nu = 1 \times 10^{-9}$ and 
$\varepsilon = 0.1$. Here, 
$Re_L$ is the integral-scale Reynolds number and 
$Re_\lambda$ is the Taylor Reynolds number.
Appendix C. Influence of the free parameter on the velocity derivative skewness
When using the 
$d$-dimensional skewness equation given in (2.51) to measure the velocity derivative skewness in our EDQNM simulations we should be careful in our interpretation of the results. The major concern we have is that since this equation depends on 
$E(k)$ then the adjustable parameter 
$\lambda _1$, through its influence on the Kolmogorov constant, may be the cause of any trend we observe. This is an important point, as the values for the Kolmogorov constant in higher dimensions cannot be known. Indeed, in this work the values used are obtained from the self-consistent LRA closure. Under this closure approximation the Kolmogorov constant decreases from three dimensions before achieving an asymptotic constant value at high dimension.
To understand the influence of 
$\lambda _!$ on our results we have performed another set of simulations where 
$\lambda _1$ is chosen such that 
$C_d = 1$ in all cases. The results of these simulations are presented in figure 9. In figure 9(a) it is clear that even when the influence of 
$\lambda _1$ on the free parameter is removed the skewness still looks to decrease at high dimension after a peak somewhere below 10 dimensions. Without data at far higher dimension we cannot speculate on the asymptotic behaviour of the velocity derivative skewness. What we can say is that, at least for the EDQNM closure, a peak value is observed which is solely determined by the triadic interactions of the system. The significance, if there is any, of this peak skewness dimension is unknown.
Figure 9. (a) EDQNM velocity derivative skewness vs dimension for 
$C_d =1$ runs. (b) Enstrophy production term 
$-S_{0} \varLambda (d)$ vs dimension for 
$C_d = 1$ runs.
Open access
