3.1. Dissipation-range turbulence with molecular fluctuations

We first probe dissipation-range intermittency by analysing the probability density function (PDF) of local vorticity obtained from direct numerical simulations (see § I of the supplementary material is available at https://doi.org/10.1017/jfm.2025.10796 for details on the numerical computation of local vorticity) averaged over at least $8\tau _\lambda$ , where $\tau _\lambda$ is the eddy turnover time. Intermittency in turbulent flows results in extreme bursts of local vorticity that are spatially interspersed within regions of relatively quiescent flow; as a result, the statistics of vorticity become highly non-Gaussian (Frisch Reference Frisch1995). This is confirmed in figure 1(a), which shows non-Gaussian tails in the PDF of the vector components of local vorticity $\omega$ normalised by the ensemble standard deviation $\sigma _{\omega }$ . Remarkably, when molecular fluctuations are included (labelled FHD), a more Gaussian-like PDF is obtained that indicates the homogenising effect of molecular fluctuations at dissipation scales that are approximately three times larger than the molecular mean free path. In FHD simulations at thermodynamic equilibrium in the absence of external turbulent forcing, the PDF is completely Gaussian. For this case, the ensemble standard deviation of local vorticity $\sigma _{\omega }^{{eq}}$ matches well with theoretical predictions of equilibrium thermodynamics (Landau & Lifshitz Reference Landau and Lifshitz1980), to within less than $1\,\%$ . The homogenising effect of molecular fluctuations is readily observed in the visualisation of local vorticity magnitude $|\omega |$ normalised by $\sigma _{\omega }^{{eq}}$ in figures 1(b) and 1(c). Whereas in deterministic simulations, regions of high vorticity are highly localised around large regions of quiescence, FHD simulations exhibit a more diffuse distribution of vorticity. Here, localised regions of high vorticity are overlaid on homogeneously distributed fluctuating velocity (and vorticity) as a result of thermal equipartition from molecular fluctuations. In FHD simulations at thermodynamic equilibrium, the local vorticity is a completely Gaussian random field (see figure S1 of the supplementary material).

Figure 1. (a) The PDFs of local vorticity $\omega$ normalised by their ensemble standard deviation $\sigma _{\omega }$ averaged over at least $8\tau _{\lambda }$ , where $\tau _{\lambda }$ is the eddy turnover time for deterministic and FHD simulations. The PDF from an FHD simulation at thermodynamic equilibrium without turbulent forcing, FHD (eq.), is also plotted. Three-dimensional visualisations of local vorticity magnitude $|\omega |$ in (b) deterministic and (c) FHD simulations. Here, $|\omega |$ is normalised by the standard deviation of vorticity fluctuations at thermodynamic equilibrium $\sigma _{\omega }^{{eq}}\approx 5\times 10^6\,\text{s}^{-1}$ ; the standard deviations of vorticity fluctuations are $\sigma _{\omega }\approx 7.3\times 10^6\,\text{s}^{-1}$ and $\sigma _{\omega }\approx 6.3\times 10^6\,\text{s}^{-1}$ for deterministic and FHD simulations, respectively.

Figure 2. (a) The PDF of local divergence $\mathcal{D}$ normalised by its ensemble standard deviation $\sigma _{\mathcal{D}}$ for deterministic and FHD simulations. The inset shows the PDF of local Mach number $\textit {Ma}$ in FHD (orange) and deterministic (blue) simulations. Three-dimensional visualisations of local divergence in (b) deterministic and (c) FHD simulations. Here, $\mathcal{D}$ is normalised by the standard deviation of divergence fluctuations that are $\sigma _{\mathcal{D}}\approx 3.1\times 10^5\,\text{s}^{-1}$ and $\sigma _{\mathcal{D}}\approx 8.7\times 10^6\,\text{s}^{-1}$ for deterministic and FHD simulations, respectively.

Compressible turbulence exhibits strong hydrodynamic shocks (Federrath et al. Reference Federrath, Klessen, Iapichino and Beattie2021); however, even weakly compressible subsonic compressible turbulent flows can exhibit large density gradients without shocks (Benzi et al. Reference Benzi, Biferale, Fisher, Kadanoff, Lamb and Toschi2008). Here, we restrict ourselves to nonlinear subsonic flows without any strong shock effects (Sagaut & Cambon Reference Sagaut and Cambon2018), but where the local Mach numbers can be as high as $0.5$ (see inset of figure 2 a) such that compressibility effects are not negligible, and we observe regions of large density variations (see figure S2 of the supplementary material for three-dimensional visualisations of local density fields). The dilatational behaviour of turbulence is analysed by the PDF of local divergence $\mathcal{D}=\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}$ normalised by the ensemble standard deviation $\sigma _{\mathcal{D}}$ in figure 2(a) (see § I of the supplementary material for details on the numerical computation of local divergence). The PDF is nearly Gaussian for FHD simulations, and is coincident with the fully Gaussian PDF for FHD simulations without turbulent forcing. Deterministic simulations exhibit modest non-Gaussian tails for both positive and negative divergence. Furthermore, the instantaneous PDFs exhibit significant temporal variability in deterministic simulations, whereas the variability is very small for FHD simulations (see figure S3 of the supplementary material for PDFs of local divergence). On average, however, divergence in deterministic simulations is negatively skewed, with skewness $\mathcal{S}\approx -0.12\pm 0.19$ , whereas $\mathcal{S}\approx 0$ for FHD simulations. More spatial volume is associated with expansion than compression in deterministic simulations (Sagaut & Cambon Reference Sagaut and Cambon2018), whereas FHD simulations exhibit nearly equal volumes of expansion and compression.

The strength of dilatation is much stronger in FHD simulations ( $\sigma _{\mathcal{D}}\approx 8.7\times 10^6\,\text{s}^{-1}$ ) compared to deterministic simulations ( $\sigma _{\mathcal{D}}\approx 3.1\times 10^5\,\text{s}^{-1}$ ). Molecular fluctuations in FHD simluations excite both vortical and dilatational modes of fluid motion via equipartition, whereas dilatational modes are indirectly excited through nonlinear coupling with the fluid vorticity in deterministic simulations (Sagaut & Cambon Reference Sagaut and Cambon2018), which is a much weaker effect for pure solenoidally forced turbulent flows considered here. In FHD simulations with no turbulent forcing, $\sigma _{\mathcal{D}}^{{eq}}\approx 8.6\times 10^6\,\text{s}^{-1}$ , which is nearly equal to its value in FHD simulations with turbulent forcing, thus demonstrating that molecular fluctuations completely dominate the dilatational dynamics. The differences are apparent in figures 2(b) and 2(c), which visualise local $\mathcal{D}/\sigma _{\mathcal{D}}$ fields. While deterministic simulations exhibit extended regions of both positive and negative divergence separated by contact discontinuities, the local divergence field is spatially nearly Gaussian in FHD simulations.

Here, we remark that in order to derive various hydrodynamic quantities, such as vorticity and divergence discussed above, we computed the numerical derivatives of the velocity field on the finite-volume grid. We emphasise that the discrete numerical operators that are used to derive these quantities are the same operators that were used to evaluate derivatives in the numerical solution algorithm for the FHD equations, thus making them consistent with the underlying numerical algorithm. As with the numerical solution of the FHD equations, the derived hydrodynamic quantities also depend on the mesh resolution; however, this resolution dependence is physically correct since the variance of thermal fluctuations depends on the scale of measurement.

3.2. Turbulence and thermal dissipation: separation of scales

In order to provide an objective comparison between deterministic Navier–Stokes and FHD simulations, we compute various microscale and dissipation (Kolmogorov) scale turbulence quantities from the simulations. Unlike deterministic Navier–Stokes equations, the computation of velocity gradients in FHD is highly scale-dependent, and as such, they do not represent an objective physical quantity. Therefore, any derived microscale and dissipation scale turbulence quantities from local velocity gradients will depend on the low-pass filter cut-off for the hydrodynamic and thermodynamic fields. In order to define an objective and meaningful turbulent energy dissipation rate, we compute the mean low-pass filtered enstrophy $\langle \varOmega ^{\lt }(k) \rangle$ from the kinetic energy spectrum $\langle E(k)\rangle =({1}/{2})\langle \hat {\boldsymbol{u}}(k) \boldsymbol{\cdot } \hat {\boldsymbol{u}}(k)^{*}\rangle$ as

(3.1) \begin{equation} \langle \varOmega ^{\lt }(k) \rangle = \int _{0}^{k} q^{2}\langle E(q) \rangle\, {\rm d}q, \end{equation}

where $\hat {\boldsymbol{u}}(k)$ is the total velocity in the Fourier space. Subsequently, a mean low-pass filtered dissipation rate is computed as $\langle \epsilon ^{\lt }(k)\rangle = ( {2\langle \eta \rangle }/{\langle \rho \rangle })\, \varOmega ^{\lt }(k)$ . Figure 3 shows $\langle \epsilon ^{\lt }(k)\rangle$ as a function of the filtering wavenumber for deterministic Navier–Stokes and FHD simulations. At large cut-off wavenumbers, $\langle \epsilon ^{\lt }(k)\rangle$ plateaus to a constant value owing to very small velocities at small scales, whereas in FHD simulations, $\langle \epsilon ^{\lt }(k)\rangle$ plateaus to a nearly similar constant value before rapidly increasing at even higher wavenumbers. This increase is attributed to dissipation primarily occurring from molecular fluctuations at small scales, which is an effect distinct from turbulent eddy fluctuations (Eyink & Jafari Reference Eyink and Jafari2022). As such, the plateau value of $\langle \epsilon ^{\lt }(k)\rangle$ , hereby denoted $\langle \epsilon ^{\lt }\rangle$ , provides a physically meaningful and objective definition of turbulent energy dissipation rate in deterministic Navier–Stokes and FHD simulations. Furthermore, the current experimental techniques for turbulence measure coarse-grained fluid velocities and dissipation rates at scales much larger than the Kolmogorov scale, and as such, are consistent with the low-pass filtered definition of these quantities. Future experiments that can measure sub-Kolmogorov-scale velocities can potentially disentangle dissipation due to molecular fluctuations from turbulence dissipation (Bandak et al. Reference Bandak, Goldenfeld, Mailybaev and Eyink2022).

Figure 3. Mean low-pass filtered dissipation rate $\langle \epsilon ^{\lt }(k)\rangle$ as a function of the wavenumber $k$ computed from the mean mean low-pass filtered enstrophy in (3.1) for deterministic Navier–Stokes and FHD simulations of compressible turbulence.

Using the prescription for low-pass filtered turbulent energy dissipation rate discussed above, we derived various microscale and dissipation-scale quantities from the simulations. In particular, we computed the following microscale quantities: (i) turbulent Mach number $\textit {Ma}_t=u'/\langle c \rangle$ , where $c$ is the local speed of sound and $u'$ is the root mean square velocity that is computed from the kinetic energy spectrum as

(3.2) \begin{equation} u'^2 = \frac {2}{3}\int _{0}^{\infty }\langle E(k) \rangle\, {\rm d}k; \end{equation}

(ii) microscale Reynolds number ${\textit{Re}}_{\lambda } = \langle \rho \rangle u' l_{\lambda } /\langle \eta \rangle$ corresponding to the Taylor microscale length (Pope Reference Pope2001)

(3.3) \begin{equation} l_{\lambda } = \sqrt {\frac {2u^{^\prime 2}}{\left \langle \left (\frac {\partial u_1}{\partial x_1}\right )^2\right \rangle }}, \end{equation}

where $\langle ({\partial u_1}/{\partial x_1})^2\rangle = ({2}/{9}) \langle \varOmega ^{\lt }\rangle$ assuming isotropy of flow. Per the discussion above, we use the plateau value of $\langle \varOmega ^{\lt }\rangle$ to estimate the velocity gradients. A microscale eddy turnover time is also computed as $\tau _{\lambda }=l_{\lambda }/u'$ . To compute dissipation-scale quantities, we use the plateau value of mean low-pass filtered dissipation rate $\langle \epsilon ^{\lt }\rangle$ as described above. The dissipation (Kolmogorov) length scale is calculated as $l_{\eta } = (\langle \eta \rangle ^3 / \langle \rho \rangle ^3 \langle \epsilon ^{\lt }\rangle )^{1/4}$ , and the corresponding Kolmogorov (small eddy turnover) time scale is calculated as $\tau _{\eta } = (\langle \eta \rangle / \langle \rho \rangle \langle \epsilon ^{\lt }\rangle )^{1/2}$ .

Table 1 lists the microscale and dissipation-scale turbulence statistics for deterministic Navier–Stokes and FHD simluations. By using a low-pass filter for velocity gradients and dissipation rates in the Fourier space as described above, we obtain a meaningful comparison between the two simulations. We note that in the case of deterministic simulation, the turbulent Mach number ${\textit{Re}}_{\lambda }$ computed above matches reasonably well with its value ${\textit{Re}}_{\lambda }=41.8$ computed directly from the velocities in the real space on the finite-volume grid. The small discrepancy between the two values is possibly attributed to complex enstrophy budgeting among the nonlinearly coupled dilatational and solenoidal components of the turbulence velocity.

Table 1. Mean turbulence statistics obtained from the simulations. Here, D-NS denotes deterministic Navier–Stokes, $\textit {Ma}_t$ is the turbulent Mach number, ${\textit{Re}}_{\lambda }$ is the microscale Reynolds number, $l_{\lambda }$ is the Taylor microscale length, $\tau _{\lambda }$ is the eddy turnover time, $l_{\eta }$ is the Kolmogorov length corresponding to the total dissipation rate, and $\tau _{\eta }$ is the Kolmogorov time scale.

3.3. Thermal energy crossover scale in the energy spectrum

We now discuss the length scales at which molecular fluctuations have an appreciable influence on compressible turbulence beyond the dissipation scale. The total energy spectra $E(k)=({1}/{2})\langle \hat {\boldsymbol{u}}(k) \boldsymbol{\cdot } \hat {\boldsymbol{u}}(k)^{*}\rangle$ of a turbulent flow can be approximately divided into the following three regimes (see figure 4 a). (i) The far-dissipation range (FDR) represents the smallest length scales, specifically wavenumbers larger than the Kolmogorov wavenumber $k_{\eta } = \nu ^{-3/4}\langle \epsilon \rangle ^{1/4}$ , where $\nu$ is the kinematic viscosity, and $\langle \epsilon \rangle$ is the total mean dissipation rate. This regime is dominated by viscous dissipation and strong intermittency (Kraichnan Reference Kraichnan1967), and molecular fluctuations strongly dominate turbulence at these length scales, as shown above. (ii) The inertial sub-range (ISR) represents length scales where energy cascades from larger eddies to smaller eddies in a scale-invariant manner, and energy spectra have the form $E(k) \propto \langle \epsilon \rangle ^{2/3}k^{-5/3}$ (Frisch Reference Frisch1995). (iii) The near-dissipation range (NDR) (Frisch & Vergassola Reference Frisch and Vergassola1991; Buaria & Sreenivasan Reference Buaria and Sreenivasan2020) that extends approximately from $k_{\eta }/30$ to $k_{\eta }$ represents the transition between ISR and FDR where the viscous effects start to become important and intermittency starts growing rapidly (Chevillard et al. Reference Chevillard, Castaing and Lévêque2005). Here, the turbulent spectra drop exponentially as $E(k) = u_{\eta }^{2} l_{\eta } \exp(-\beta k l_{\eta })$ , where $u_{\eta } = (\langle \epsilon \rangle \nu )^{1/4}$ is the Kolmogorov velocity scale, $l_{\eta } = (\nu ^{3}/\langle \epsilon \rangle )^{1/4}$ is the Kolmogorov length, and $\beta$ is the rate of exponential decay of the spectrum that typically ranges from 3 to 7 (we have fixed $\beta =5$ in our analysis) (Khurshid, Donzis & Sreenivasan Reference Khurshid, Donzis and Sreenivasan2018).

Figure 4. (a) Comparison of the total kinetic energy spectrum $\langle E(k) \rangle$ in FHD versus deterministic simulations. Three approximate ranges of length scales are highlighted: inertial sub-range (ISR, in blue), near-dissipation range (NDR, in pink) and far-dissipation range (FDR, in green). In FHD simulations the thermal spectrum $E_{\textit{th}}(k) = ({3k_B \langle T \rangle }/{2\langle \rho \rangle }) 4\unicode{x03C0} k^{2}$ (red dash-dotted line) dominates for wavenumbers larger than the thermal crossover scale $k_{\textit{th}}$ , where $k_B$ is the Boltzmann constant. (b) Standard deviation in total kinetic energy spectrum $\delta E(k) = \langle (E(k) - \langle E(k)\rangle )^{2} \rangle ^{1/2}$ normalised by $\langle E(k) \rangle$ .

Molecular fluctuations introduce another length scale in the turbulence spectrum (Bandak et al. Reference Bandak, Eyink, Mailybaev and Goldenfeld2021). From equilibrium thermodynamics, the contribution of molecular fluctuations to the energy spectrum (assuming no net flow, i.e. $\langle \boldsymbol{u}\rangle = 0$ ) is

(3.4) \begin{equation} E_{\textit{th}}(k) = \frac {3k_B \langle T \rangle }{2\langle \rho \rangle } 4\unicode{x03C0} k^{2}, \end{equation}

which is ‘equipartitioned’ white noise with variance ${3k_B \langle T \rangle }/{2\langle \rho \rangle }$ at all scales. The wavenumber $k_{\textit{th}}$ at which molecular fluctuations are approximately equal in magnitude to the turbulent spectrum is (Bandak et al. Reference Bandak, Eyink, Mailybaev and Goldenfeld2021)

(3.5) \begin{equation} u_{\eta }^{2} l_{\eta } \exp\left (-\beta k_{\textit{th}} l_{\eta }\right ) \approx \frac {k_B \langle T \rangle }{\langle \rho \rangle } k_{\textit{th}}^{2}. \end{equation}

Indeed, in figure 4(a), we observe that for FHD simulations, the total energy spectrum crosses over from an exponential decay in the NDR to being dominated by the thermal spectrum $E_{\textit{th}}(k)$ at high wavenumbers. The agreement with $E_{\textit{th}}(k)$ is remarkable without any fitting parameters. The thermal crossover wavenumber $k_{\textit{th}}$ is approximately three times smaller than the Kolmogorov wavenumber $k_{\eta }$ , and its predicted value from (3.5) (shown by the dashed vertical black line) matches well with the observed crossover to $E_{\textit{th}}(k)$ (shown by the dash-dotted red line). While the ratio $k_{\textit{th}}/k_{\eta }$ depends on turbulence conditions, such as density, viscosity, temperature and mean dissipation rate (Bandak et al. Reference Bandak, Goldenfeld, Mailybaev and Eyink2022; Bell et al. Reference Bell, Nonaka, Garcia and Eyink2022), the relationship between $k_{\textit{th}}$ and $k_{\eta }$ is fairly robust, and varies only very marginally across a wide range of turbulence conditions (Bell et al. Reference Bell, Nonaka, Garcia and Eyink2022).

The crossover into the thermal regime is also observed for the dilatational part of the energy spectrum $E_{d}(k)=({1}/{2})\langle \hat {\boldsymbol{u}}_{d}(k) \boldsymbol{\cdot } \hat {\boldsymbol{u}}_{d}(k)^{*}\rangle$ , as shown in figure 5(a), where $\hat {\boldsymbol{u}}_{d}$ is the dilatational (curl-free) part of the total velocity $\hat {\boldsymbol{u}}$ . At low wavenumbers, the total kinetic energy is dominated by solenoidal modes since the external turbulence forcing is solenoidal (see figure S4 of the supplementary material for $\langle E_d(k)/E(k) \rangle$ ). However, following a rapid decay in the NDR, $E_{d}(k)$ crosses over to $E_{d,{\textit{th}}}(k)=(1/3)\,E_{\textit{th}}(k)$ at the wavenumber $k_{\textit{th}}$ in FHD simulations. The factor $1/3$ appears because one-third of the thermal energy of molecular fluctuations is ‘equipartitioned’ into the dilatational part, and two-thirds into the solenoidal part of the total kinetic energy.

The picture that emerges from these observations is that the impact of molecular fluctuations on turbulence is not limited to dissipation scales in the FDR, but appears at larger thermal crossover scales in the NDR. While the simulations in this study have been conducted at low Reynolds numbers due to computational constraints, we can estimate the scales at which molecular fluctuations will be significant in several practical scenarios. For example, following Garratt (Reference Garratt1994) and Bandak et al. (Reference Bandak, Eyink, Mailybaev and Goldenfeld2021), in an atmospheric boundary layer assumed to be composed entirely of nitrogen at $T=300\,\text{K}$ , the energy dissipation rate is $\epsilon =400\,\text{cm}^{2}\,\,\text{s}^{- 3}$ , kinematic viscosity of nitrogen is $\nu =0.16\,\text{cm}^{2}\,\text{s}^{- 1}$ , and density is $\rho =1.1\times 10^{-3}\text{g}\,\text{cm}^{- 3}$ . The mean free path is $l_{\textit{mfp}}\approx 70\,\text{nm}$ , while the Kolmogorov length scale is $l_{\eta }=0.57\,\text{mm}$ . From (3.2), the thermal crossover length scale at which molecular fluctuations will dominate is $l_{\textit{th}}\approx 1.3\,\text{mm}$ , which is over four orders of magnitude larger than the mean free path.

Figure 5. (a) Comparison of dilatational kinetic energy $\langle E_d(k) \rangle$ in FHD versus deterministic simulations. The FHD simulations transition over to the thermal energy spectrum is $E_{d,{\textit{th}}}(k)=(1/3)\,E_{\textit{th}}(k)$ (red dash-dotted line) at $k_{\textit{th}}$ . (b) Standard deviation in the dilatational kinetic energy spectrum $\delta E_d(k) = \langle (E_d(k) - \langle E_d(k)\rangle )^{2} \rangle ^{1/2}$ normalised by $\langle E_d(k) \rangle$ .

3.4. Molecular fluctuations impact turbulence statistics across the near-dissipation range

It is apparent that mean turbulence properties are significantly modified in the NDR at all length scales smaller than $1/k_{\textit{th}}$ . However, it is well known that intermittency in turbulence starts building up in the ISR, and rapidly increases in the NDR, where viscous effects start to intensify (Frisch & Vergassola Reference Frisch and Vergassola1991; Chevillard et al. Reference Chevillard, Castaing and Lévêque2005). Therefore, even though molecular fluctuations do not affect the ensemble-averaged turbulence properties such as the energy spectrum $\langle E(k) \rangle$ for $k\lt k_{\textit{th}}$ , we can expect them to modify the statistical properties of turbulence.

Indeed, a remarkable picture emerges where the large temporal statistical variability of turbulence in the NDR is significantly reduced due to molecular fluctuations. Figures 4(b) and 5(b) respectively show the standard deviation of the total energy $\delta E(k)$ and dilatational energy spectra $\delta E_d(k)$ normalised by the mean value averaged over at least $8\tau _{\lambda }$ . The growth of $\delta E(k)$ and $\delta E_d(k)$ is much slower in FHD than in deterministic simulations for $k\lt k_{\textit{th}}$ , thus implying increased statistical stability of the dynamical turbulent system with molecular fluctuations. For $k\gt k_{\textit{th}}$ , the statistical variability plummets by two orders of magnitude in FHD simulations, whereas it keeps increasing with $k$ for deterministic simulations up to the beginning of the FDR. The eventual drop-off in $\delta E(k)$ and $\delta E_d(k)$ at very high $k$ results from limitations in numerical precision.

Next, we quantify scale-dependent spatial intermittency of turbulence through high-pass filtered skewness $\mathcal{S}^{\gt }(k)$ and kurtosis (flatness) $\mathcal{K}^{\gt }(k)$ of the velocity gradient $\partial _{x} \boldsymbol{u}^{\gt }$ that are computed as

(3.6) \begin{equation} \mathcal{S}^{\gt }(k_i) = \frac {\overline {\left (\partial _{x} \boldsymbol{u}^{\gt }\right )^{3}}}{\left [\overline {\left (\partial _{x} \boldsymbol{u}^{\gt }\right )^{2}}\right ]^{3/2}}, \quad \mathcal{K}^{\gt }(k_i) = \frac {\overline {\left (\partial _{x} \boldsymbol{u}^{\gt }\right )^{4}}}{\left [\overline {\left (\partial _{x} \boldsymbol{u}^{\gt }\right )^{2}}\right ]^{2}}, \end{equation}

where

(3.7) \begin{equation} \overline {\left (\partial _{x} \boldsymbol{u}^{\gt }\right )^{n}} = \frac {1}{V}\int \text{d}\boldsymbol{r}\, \left (\partial _{x} \boldsymbol{u}^{\gt }(\boldsymbol{r})\right )^{n}, \end{equation}

and $\boldsymbol{u}^{\gt }$ is the high-pass filtered velocity. Numerically, $\boldsymbol{u}^{\gt }$ is obtained by first computing the discrete Fourier transform of the velocity field over the finite-volume grid and zeroing out the Fourier modes for wavenumbers smaller than $k$ , followed by a discrete inverse Fourier transform to obtain the high-pass filtered velocity on the same finite-volume grid. Once $\boldsymbol{u}^{\gt }$ is obtained, $\mathcal{S}^{\gt }(k)$ and $\mathcal{K}^{\gt }(k)$ are calculated by numerically computing the derivative $\partial _{x} \boldsymbol{u}^{\gt }$ using the same gradient operators as employed in the numerical simluation of the FHD equations.

In an intermittent dynamical system, $\mathcal{K}^{\gt }(k)$ is expected to grow unboundedly with $k$ in the NDR and into the FDR as regions of intense turbulent activity become increasingly localised in smaller fractions of the system volume (Frisch Reference Frisch1995). A negative skewness for a turbulent system implies energy cascade from large to small scales (Frisch Reference Frisch1995), and its magnitude ranges from $\mathcal{S}\approx -0.5$ to $\mathcal{S}\approx -0.3$ . In a fully Gaussian distribution, $\mathcal{S}=0$ and $\mathcal{K}=3$ .

Figure 6. (a) Filtered kurtosis (flatness) $\mathcal{K}^{\gt }(k)$ and (b) filtered skewness $\mathcal{S}^{\gt }(k)$ of the velocity gradient $\partial _{x} \boldsymbol{u}^{\gt }$ , where $\boldsymbol{u}^{\gt }$ is the high-pass filtered velocity obtained by zeroing out all the Fourier modes for wavenumbers lesser than $k$ in the velocity field. The horizontal dashed line corresponds to the kurtosis and skewness of a Gaussian random field with $\mathcal{K}^{\gt }=3$ and $\mathcal{S}^{\gt }=0$ for all wavenumbers. The error bars denote the ensemble standard deviation.

In the present simulations, rapidly increasing intermittency from its build up in the ISR and propagation through the NDR and into the FDR is observed in the deterministic case, as seen by the variation of $\mathcal{K}^{\gt }$ in figure 6(a). In a remarkable contrast, $\mathcal{K}^{\gt }(k)\approx 3$ at all wavenumbers in FHD simulations, thus demonstrating that the intermittent dynamics is completely inhibited not just in the FDR, but well into the NDR. Furthermore, large variations in $\mathcal{K}^{\gt }(k)$ in deterministic simulations at high $k$ , which are indicative of highly intermittent behaviour, are not observed in FHD simulations. On the other hand, the skewness of velocity gradient $\mathcal{S}^{\gt }(k)$ in figure 6(b) saturates to its Gaussian value, as expected, for both FHD and deterministic simulations at high $k$ . However at low $k$ , deterministic simulations exhibit a negative skewness with large variability, whereas it is of a much smaller magnitude and variability in FHD simulations. We note that in a recent study on the role of molecular fluctuations in incompressible turbulence (Bell et al. Reference Bell, Nonaka, Garcia and Eyink2022), the skewness and kurtosis of the velocity gradient were reported to be unaffected by molecular fluctuations. Furthermore, through the analysis of structure functions in recent studies on stochastic shell modelling of incompressible turbulence (Bandak et al. Reference Bandak, Goldenfeld, Mailybaev and Eyink2022) and molecular gas dynamics simulations of compressible turbulence (McMullen, Torczynski & Gallis Reference McMullen, Torczynski and Gallis2023), it was observed that while the far-dissipation range intermittency is replaced by Gaussian fluctuations, the intermittency in the intermediate range persists. While our results are consistent with the studies in the far-dissipation range, our observations of drastically reduced intermittency in the near-dissipation range can potentially be attributed to low- ${\textit{Re}}$ flows simulated here and/or compressibility effects.

Figure 7. Cross-sectional visualisations of the local vorticity magnitude $|\omega |$ (normalised by the ensemble mean $\langle |\omega |\rangle$ ) only for wavenumbers $k\lt k_{\textit{th}}$ in (a) deterministic and (b) FHD simulations. (c,d) Same as (a,b), respectively, but only for wavenumbers $k\gt k_{\textit{th}}$ . Cross-sectional visualisation of the local divergence $\mathcal{D}$ (normalised by the ensemble standard deviation $\sigma _{\mathcal{D}}$ ) only for wavenumbers $k\lt k_{\textit{th}}$ in (e) deterministic and ( f) FHD simulations. (g,h) Same as (e,f), respectively, but only for wavenumbers $k\gt k_{\textit{th}}$ .

A visual analysis of the filtered invariants of velocity gradient (i.e. vorticity magnitude $|\omega |$ and divergence $\mathcal{D}$ ) highlights our observations. Figures 7(a) and 7(b) show two-dimensional slices of vorticity magnitude $|\omega |$ , and figures 7(e) and 7( f) show two-dimensional slices of divergence $\mathcal{D}$ filtered for wavenumbers $k\lt k_{\textit{th}}$ . Similarly, figures 7(c) and 7(d), and figures 7(g) and 7(h), show the same data but filtered for wavenumbers $k\gt k_{\textit{th}}$ . While these fields ‘appear’ similar at large wavelengths, $k\lt k_{\textit{th}}$ , in FHD and deterministic simulations, the visual differences are significant wavenumbers $k\gt k_{\textit{th}}$ . Here, FHD simulations exhibit a nearly homogeneous spatial distribution of vorticity and divergence with no signs of intermittency, whereas deterministic simulations exhibit classic signs of dissipation-range intermittency with localised bursts of high vorticity and divergence in a ‘sea’ of quiescent fluid.