4.1. Basic statistics

We start with a few results that characterise the basic properties of the flows and the synchronisation process. For reference, the average energy spectra are documented in figure 1, which is compared with the $k^{-5/3}$ scaling law. To monitor the synchronisation between ${\boldsymbol {u}}^{(M)}$ and ${\boldsymbol {u}}^{(S)}$ , we use the synchronisation error

(4.1) \begin{equation} {\mathcal E}_{MS}(t) = \Vert {\boldsymbol {u}}^{(M)}-{\boldsymbol {u}}^{(S)}\Vert, \end{equation}

which will decay exponentially when the two flows synchronise, and the rate of decay of ${\mathcal E}_{MS}$ is related to the conditional LLE $\lambda (k_m)$ , as having been shown by Henshaw et al. (Reference Henshaw, Kreiss and Ystróm2003), Yoshida et al. (Reference Yoshida, Yamaguchi and Kaneda2005), Inubushi et al. (Reference Inubushi, Saiki, Kobayashi and Goto2023) and Li et al. (Reference Li, Tian, Li, Si and Mohammed2024b ). Figure 2(a) compares the results for ${\mathcal E}_{MS}(t)$ (symbols) with $\exp (\Lambda t/\tau )$ (lines), where

(4.2) \begin{equation} \Lambda (k_m) \equiv \lambda (k_m) \tau \end{equation}

is the conditional LLE non-dimensionalised with $\tau$ . Some small discrepancies can be seen between the two quantities, which we attribute to statistical uncertainty. This result is consistent with previous findings (Nikolaidis & Ioannou Reference Nikolaidis and Ioannou2022; Li et al. Reference Li, Tian, Li, Si and Mohammed2024b ), which shows that the non-dimensional decay rate of the synchronisation error can be given by $\Lambda$ .

Figure 1.

Mean energy spectra $\langle E \rangle (k)$ of the base flow. Dash-dotted line is the $k^{-5/3}$ power law.

Figure 2.

(a) Comparison between the decay rates of the synchronisation error ${\mathcal E}_{MS}(t)$ and the conditional LLEs. (b) Normalised conditional LLEs $\Lambda \equiv \lambda \tau$ , with the short-dashed line showing $-(k_m \eta )^2$ and the vertical line marking $k_m\eta =0.2$ . Only four data points are shown for the case F1. The empty circles are the decay rate data from Yoshida et al. (Reference Yoshida, Yamaguchi and Kaneda2005) (see the main text for more detail).

Figure 3.

Production $P$ , dissipation $D$ and $P {-} D$ as functions of $k_m\eta$ . $P$ , lines with both solid and empty squares; $D$ , lines with both solid and empty triangles; $P{-}D$ , lines with both solid and empty diamonds. Solid lines, R1; long dashed lines, R2; dash-dotted lines, R3; short-dashed lines, R4.

Figure 2(b) shows $\Lambda$ as functions of $k_m\eta$ . The conditional LLEs are known to depend on the Reynolds number weakly (see, e.g., Inubushi et al. Reference Inubushi, Saiki, Kobayashi and Goto2023). Here, the curves for different $Re_\lambda$ show some differences at small $k_m\eta$ , but they collapse on each other for larger $k_m \eta$ (for, e.g. $k_m\eta \geqslant 0.15$ ) because $Re_\lambda$ varies only mildly. The conditional LLEs decrease as $k_m \eta$ increases, but the variation appears to be small for small $k_m\eta$ . The threshold wavenumber $k_c$ , for which $\Lambda = 0$ , is found to be $k_c \eta \approx 0.2$ , which is also the value obtained by Yoshida et al. (Reference Yoshida, Yamaguchi and Kaneda2005) for isotropic turbulence. The dashed line without symbols is the function $-(k_m\eta )^2$ , which is an estimate of the non-dimensional conditional LLE when the evolution of the velocity perturbation is determined solely by viscous diffusion (see, e.g., Inubushi et al. Reference Inubushi, Saiki, Kobayashi and Goto2023). We will discuss this estimate below together with figure 3.

Figure 2(b) also includes data from two other sources to cross-check the results from groups R1–R4. The conditional LLEs from case F1 at four different $k_m \eta$ are plotted with empty squares. They display only negligible differences with those from R1, which shows that the results are essentially grid independent. The empty circles are calculated from the decay rates $\tilde {\alpha }$ obtained by Yoshida et al. (Reference Yoshida, Yamaguchi and Kaneda2005), plotted in figure 3 therein. More specifically, the empty circles are the values of $-\tilde {\alpha }/2$ . Per Yoshida et al. (Reference Yoshida, Yamaguchi and Kaneda2005), $\tilde {\alpha }$ is defined by $\Vert {\boldsymbol {u}}^{(M)} - {\boldsymbol {u}}^{(S)}\Vert ^2 \sim \exp (- \tilde {\alpha } t/\tau )$ . Since $\Vert {\boldsymbol {u}}^{(M)} - {\boldsymbol {u}}^{(S)}\Vert \sim \exp (\Lambda t/\tau )$ , as is shown in figure 2(a), we expect $\Lambda = - \tilde {\alpha }/2$ . This relation is verified by the data, as the empty circles fall closely on the curves for $\Lambda$ . Note that the empty circles correspond to a wide range of Reynolds numbers. Also, $\tilde {\alpha }$ is obtained by measuring the decay rate of $\Vert {\boldsymbol {u}}^{(M)} - {\boldsymbol {u}}^{(S)}\Vert ^2$ , which is a procedure that is very different from how $\Lambda$ is calculated (Yoshida et al. Reference Yoshida, Yamaguchi and Kaneda2005). Thus, the agreement between $\Lambda$ and $-\tilde {\alpha }/2$ is a strong validation of our results.

Further insights on $\Lambda$ can be explored according to (2.17). We use

(4.3) \begin{equation} P \equiv \tau \left \langle {\mathcal P}_e \right \rangle, \qquad D \equiv \tau \left \langle {\mathcal D}_e \right \rangle \end{equation}

to denote the non-dimensional mean production and mean dissipation, respectively. It follows from (2.17) that

(4.4) \begin{equation} \Lambda = P - D. \end{equation}

Note that a pre-requisite for (4.4) is that the external forcing has no impact on the evolution of the velocity perturbation (cf. (2.6) and the comments following (2.7)). Otherwise, the equation would contain a forcing term.

The values of $P$ , $D$ and their difference $P-D$ are shown in figure 3 for different coupling wavenumber $k_m$ . As dictated by (4.4), the curves for $P-D$ agree precisely with those of the conditional LLEs $\Lambda$ shown in figure 2(b). The main observation about figure 3 is that $P$ decreases, whilst $D$ increases, with increasing $k_m$ . The two curves intersect at the threshold wavenumber $k_c \eta$ . Both seem to contribute roughly equally to the change in $\Lambda$ as $k_m$ varies. The results are slightly different for different $Re_\lambda$ .

We now return to the discussion of the estimate $-(k_m\eta )^2$ for $\Lambda (k_m)$ shown in figure 2. It was found by Inubushi et al. (Reference Inubushi, Saiki, Kobayashi and Goto2023) that $\Lambda (k_m)$ is always larger than $-(k_m\eta )^2$ , and approaches $-(k_m \eta )^2$ from above (see Figure 3 therein). These trends suggest that $P$ becomes negligible for large $k_m\eta$ , and $D$ approaches the pure viscous estimate as $k_m\eta$ increases. Our results, however, appear to display different behaviours, as we can see from figures 2 and 3. First, figure 2 shows that $-(k_m\eta )^2$ can be larger than $\Lambda (k_m)$ in our simulations. Second, figure 3 shows that, though $P$ decreases as $k_m\eta$ increases, it is still significantly larger than zero for the largest $k_m \eta$ , even when $\Lambda (k_m )$ is already smaller than $-(k_m \eta )^2$ . In sum, the dissipation contribution in our simulations is significantly higher than what is implied by the estimate $-(k_m\eta )^2$ . We will discuss this further in § 4.2 when we look into the spectral dynamics of the LLVs.

Finding the threshold wavenumber $k_c$ via the conditional LLEs requires considerable computational cost, as it entails calculating $\lambda (k_m)$ for many $k_m$ . On this issue, Li et al. (Reference Li, Tian, Li, Si and Mohammed2024b ) made an interesting observation in the context of rotating turbulence. They found that the average energy spectrum of the unconditional LLV peaks at a wavenumber which appears to be close to, or the same as, $k_c$ . The observation is reproduced by Li et al. (Reference Li, Si, Li and Xu2024a ) for the synchronisation of LES of periodic turbulence. Elementary results for the spectra of the LLV in the present study, i.e. $E_v(k,t)$ , are shown in figure 4. In figure 4(a), the early evolution of $E_v(k,t)$ is plotted. We initialise the Fourier modes of the perturbation with independent random numbers with identical probability distribution functions (p.d.f.s). As a consequence, $E_v(k,t=0) \sim k^2$ for large $k$ , as the number of modes in the spherical shell with radius $k$ is proportional to the area $4\pi k^2$ of the shell, as shown by the green dashed line. The spectrum at time $t$ , $E_v(k,t)$ , exhibits a period of transient evolution, as depicted by the think red lines, which show that the peak of the spectrum moves towards lower wavenumbers. At $t \approx 20 \tau$ , the spectrum converges towards a distribution shown with the thick black line, which then fluctuates over time. The long time average $\langle E_v \rangle$ is shown in figure 2(b). Clearly, the peaks of the spectra are found at $k\eta \approx 0.2$ , reproducing previous findings. The flow here is very different from the rotating turbulence investigated by Li et al. (Reference Li, Tian, Li, Si and Mohammed2024b ). For example, the base flow energy spectrum $\langle E \rangle (k)$ follows the $k^{-2}$ or $k^{-3}$ power laws for Li et al. (Reference Li, Tian, Li, Si and Mohammed2024b ), whereas here, it follows the canonical $k^{-5/3}$ scaling. Therefore, it is non-trivial for the same relationship to hold in both cases. As a step towards understanding the origin of this relationship, we look into the spectral dynamics of the LLVs and the conditional LLVs in what follows.

Figure 4.

(a) Short time evolution of $E_v(k,t)$ for $t$ between $0$ and approximately $ 20 \tau$ , with the dashed line showing $E_v(k,t=0)$ and the black thick line showing $E_v(k,t)$ for $t \approx 20\tau$ . The arrow indicates the increasing direction of time. (b) Average energy spectra $\langle E_v\rangle$ for the unconditional LLV, with the vertical dashed line marking the value $k\eta =0.2$ .

4.2. Production and dissipation: spectral analyses

To understand the dependence of $P$ and $D$ (hence $\Lambda$ ) on $k_m$ better, we look into the production spectrum ${\mathcal P}_v(k,t)$ and the dissipation spectrum ${\mathcal D}_v(k,t)$ . We first consider their non-dimensional ensemble averages

(4.5) \begin{equation} P_v(k) \equiv \tau \langle {\mathcal P}_v (k,t)\rangle, \qquad D_v(k) \equiv \tau \langle {\mathcal D}_v(k,t) \rangle . \end{equation}

The expression for $D_v(k)$ can also be written as

(4.6) \begin{equation} D_v(k) = 2 \bar {\nu } k^2 \tau \langle E_v \rangle (k) = 2 (k\eta )^2 \langle E_v \rangle (k). \end{equation}

The behaviours of $P_v$ and $D_v$ are related by (2.27). Taking the average of (2.27) and noting $\partial _t \langle E_v \rangle = 0$ , we find

(4.7) \begin{equation} P_v(k) - D_v(k) = 2\tau \langle \gamma E_v\rangle . \end{equation}

Our data show that $\gamma$ is essentially uncorrelated to $E_v$ (figure omitted). Therefore, $\langle \gamma E_v\rangle \approx \langle \gamma \rangle \langle E_v \rangle = \lambda \langle E_v \rangle$ and we obtain

(4.8) \begin{equation} P_v (k) = D_v(k) + 2 \Lambda \langle E_v \rangle = 2 \left [ (k\eta )^2 + \Lambda (k_m) \right ] \langle E_v \rangle . \end{equation}

The equation delineates the spectral balance of the energetics of the velocity perturbation. Integrating (4.8) over $k$ , we obtain

(4.9) \begin{equation} \Lambda (k_m) = \int _0^\infty P_v(k) {\rm d}k - \int _0^\infty 2 (k\eta )^2 \langle E_v \rangle (k) {\rm d}k, \end{equation}

which makes clear that (4.9) is the spectral version of (4.4).

Figure 5.

(a) Production spectrum $P_v$ , dissipation spectrum $D_v$ and $2\Lambda \langle E_v \rangle$ as functions of $k\eta$ for $k_m=0$ with the dashed line showing the residual $P_v - D_v - 2\Lambda \langle E_v \rangle$ ; (b) $P_v$ ; (c) $ D_v$ ; (d) $\langle E_v \rangle$ . The non-dimensionalised coupling wavenumbers $k_m\eta$ are $0$ , $0.12$ , $0.17$ , $0.22$ , $0.26$ and $0.31$ . For group R3.

Figure 6.

Same as figure 5 but for group R4. The non-dimensionalised coupling wavenumbers $k_m\eta$ are $0$ , $0.08$ , $0.14$ , $0.21$ , $0.27$ and $0.34$ .

Figure 5 shows the spectra $P_v$ , $D_v$ and $\langle E_v \rangle$ with different $k_m$ for the cases in group R3. The same results for group R4 are shown in figure 6 as corroboration. The two figures depict the same behaviours. Therefore, we will only discuss figure 5 in detail. Figure 5(a) compares the three distributions for $k_m=0$ , i.e. for the cases where no coupling is imposed. The dashed line shows the residual $P_v - D_v - 2\Lambda \langle E_v \rangle$ , which, according to (4.8), should be essentially zero. Though it is not exactly zero, the dashed line shows that it is negligible compared with the dominant terms at all wavenumbers. Not surprisingly, the dissipation $D_v$ peaks at a higher wavenumber compared with $P_v$ and $2\Lambda \langle E_v\rangle$ , as dissipation dominantly comes from small scales. On the other hand, $2\Lambda \langle E_v \rangle$ peaks at $k\eta \approx 0.2$ (same as that of $\langle E_v \rangle$ ). The peak of $P_v$ is found at a wavenumber in between, which shows that the strongest production is found at wavenumbers well inside the dissipation range of the base flow. Also, $P_v(k)$ is positive definite, implying that the velocity perturbation is amplified by the production mechanism at all scales.

The production spectra for several different $k_m$ are shown in figure 5(b). Because $\hat {\boldsymbol {v}}(\boldsymbol {k},t) = 0$ for $\vert \boldsymbol {k}\vert \lt k_m$ due to the coupling between the synchronised flows, $P_v(k) = 0$ for $k \lt k_m$ . What is noteworthy is that $P_v(k)$ is also reduced by the coupling for $k\gt k_m$ and it is increasingly smaller for larger $k_m$ . The attenuating effect of the coupling is localised in the wavenumbers around $k_m$ , with $P_v(k)$ at large $k$ little affected. To interpret this feature, we refer back to (2.29). Since $\hat {\boldsymbol {v}}(\boldsymbol {q},t) = 0$ when $\vert \boldsymbol {q}\vert \lt k_m$ , the summation in (2.29) does not include the Fourier modes with $\vert \boldsymbol {q} \vert \lt k_m$ . Thus, increasing $k_m$ means excluding more Fourier modes from the summation, which thus likely leads to smaller $P_v(k)$ , because the panel (a) shows that the contributions from these excluded modes are likely to be positive.

Another observation is that, though $P_v(k)$ is positive definite for $k \gt k_m$ in most cases, it assumes negative values for some wavenumbers when $k_m$ is large enough (e.g. for case R $3$ K $13$ ). This observation can be understood from (4.8). Obviously, $P_v(k)$ would be positive definite for $k\gt k_m$ if $(k\eta )^2 + \Lambda (k_m ) \gt 0$ for all $k\gt k_m$ . This is clearly satisfied if $\Lambda (k_m) \gt 0$ , which is the case when $k_m=0$ or $k_m$ is small. If $\Lambda (k_m) \lt 0$ , however, then $(k\eta )^2 + \Lambda \lt 0$ when $k\eta \lt (-\Lambda )^{1/2}$ . As a result, $P_v(k)$ would be negative for wavenumber $k$ if $k_m\eta \lt k\eta \lt (-\Lambda )^{1/2}$ . The inequality can be satisfied by some wavenumbers if $ - (k_m \eta )^2 \gt \Lambda (k_m)$ , which is satisfied when $k_m\eta$ is large enough as shown in figure 2(b). Therefore, when $k_m$ is large enough, the velocity perturbations at the scales just below the coupling scale would be suppressed by the production term.

The dissipation spectrum $D_v$ is shown in figure 5(c), while the energy spectrum $\langle E_v \rangle$ is shown in figure 5(d). As $D_v = 2 (k\eta )^2 \langle E_v \rangle$ , the two parameters are similar in many ways. The most conspicuous feature for both is that their distributions are elevated as $k_m$ increases. This behaviour could be a simple consequence of the normalisation condition $\langle e \rangle _v = 1/2$ , which fixes the total integral of $\langle E_v \rangle$ . As the support of $\langle E_v \rangle$ is reduced when $k_m$ increases, its values have to increase. The values of $D_v(k)$ , as a consequence, have to increase too. Nevertheless, there is another mechanism by which $D_v(k)$ is enhanced. Equation (4.8) implies that

(4.10) \begin{equation} D_v(k) = P_v(k) \frac {(k\eta )^2}{(k\eta )^2 + \Lambda (k_m)}. \end{equation}

Equation (4.10) shows that reduced $P_v(k)$ tends to reduce $D_v(k)$ . However, reduced $P_v(k)$ also tends to reduce $\Lambda (k_m)$ , which in turns enhances the factor $ (k\eta )^2/[\Lambda (k_m) + (k\eta )^2]$ , thus potentially increases $D_v(k)$ . That is, there is a mechanism by which $D_v(k)$ increases as a consequence of reduced $P_v(k)$ . This effect is stronger at lower wavenumbers as the factor $ (k\eta )^2/[\Lambda (k_m) + (k\eta )^2]$ is more sensitive to the change in $\Lambda (k_m)$ when $k\eta$ is smaller. Unfortunately, it is unclear which of the above two mechanisms contributes more to the enhancement of $D_v(k)$ .

We now turn to a brief discussion on the viscous estimate $-(k_m\eta )^2$ for $\Lambda (k_m)$ . Note that $\langle E_v \rangle (k) = 0$ for $k\lt k_m$ . Therefore,

(4.11) \begin{align} D & = \int _{k_m}^\infty 2 (k\eta )^2 \langle E_v \rangle {\rm d}k = 2(k_m\eta )^2 \int _{k_m}^\infty (k/k_m)^2 \langle E_v \rangle {\rm d}k \notag \\ &\geqslant 2 (k_m\eta )^2 \int _{k_m}^\infty \langle E_v \rangle {\rm d}k = (k_m\eta )^2, \end{align}

which implies that $(k_m\eta )^2$ tends to underestimate the dissipation $D$ . Therefore, it might not be surprising that in our simulations, $\Lambda (k_m)$ becomes smaller than $-(k_m\eta )^2$ for some $k_m$ . The estimate $(k_m\eta )^2$ would, however, become exact if $\langle E_v\rangle$ was proportional to the Dirac delta function concentrated at $k_m$ (with the strength being $1/2$ ). Our results in figure 5(d) do show a tendency for $\langle E_v \rangle$ to concentrate around $k_m$ as $k_m$ increases. The argument can be recast in more physical terms. Let $k_d$ be the peak wavenumber for $D_v(k)$ . Figures 5 and 6 show that $k_d \gt k_m$ in our simulations (cf. panels c therein). In this case, one may estimate $D$ as follows:

(4.12) \begin{equation} D \sim 2 (k_m \eta )^2 \int _{k_m}^\infty \frac {k^2_d}{k^2_m}\langle E_v \rangle (k) {\rm d}k = (k_d \eta )^2, \end{equation}

assuming the dissipation dominantly comes from the wavenumbers around $k_d$ . If the production term is negligible, it then follows from (4.4) that $ \Lambda \sim - (k_d \eta )^2 \lt -(k_m\eta )^2$ for some $k_m$ . This argument is consistent with our results, although we note that we find $\Lambda \lt -(k_m\eta )^2$ for large $k_m$ without neglecting the production term, as shown in figures 2(b) and 3. On the other hand, it is possible that $D_v(k)$ peaks at $k_m$ for cases with sufficiently large $k_m$ . In this case, the dissipation would be concentrated around $k_m$ , and an argument in the spirit of (4.12) would lead to $D \sim (k_m\eta )^2$ . If the production term is negligible, one then obtains $\Lambda (k_m) \sim -(k_m \eta )^2$ , as proposed by Inubushi et al. (Reference Inubushi, Saiki, Kobayashi and Goto2023).

We have been able to present some semi-analytical discussions of the results in figure 5 based on (4.8). Deriving the relationship between the coupling wavenumber $k_c\eta$ and the peak wavenumber of $\langle E_v \rangle$ analytically requires, as a very first step, finding an analytical expression for the peak wavenumber for $\langle E_v \rangle$ at $k_m=0$ . Attempts at such analyses, however, quickly run into the classical closure problem, as we do not have an expression for $P_v(k)$ in terms of $\langle E_v \rangle$ . Nevertheless, some elementary results can be obtained. Since the peak wavenumber is given by ${\rm d} \langle E_v \rangle /{\rm d} (k\eta ) = 0$ , we can find from $\langle E_v \rangle = P_v(k)/[2(k\eta )^2 + 2 \Lambda (0)]$ that the peak wavenumber $k\eta$ satisfies

(4.13) \begin{equation} \frac { 2 k\eta }{ (k \eta )^2 + \Lambda (0)} = \frac {1}{P_v} \frac { {\rm d} P_v}{{\rm d} (k\eta )} = \frac { {\rm d} \ln P_v}{{\rm d} (k\eta )}, \end{equation}

which can be solved for $k\eta$ if we have an expression for $P_v(k)$ . Purely as a demonstration, we let $P_v(k) \sim a (k\eta )^b$ , assuming the expression provides a good approximation for the slope of $\ln P_v(k)$ around the peak wavenumber for some constants $a$ and $b$ . The peak wavenumber is then given by

(4.14) \begin{equation} k\eta = \left (\frac {b \Lambda (0)}{2-b}\right )^{1/2}. \end{equation}

If we let $\Lambda (0) = 0.12$ (hence ignoring $\Lambda (0)$ ’s dependence on $Re_\lambda$ ), then $b=1/2$ would give $k\eta = 0.2$ . Progress may be made by developing an EDQNM-type model for $P_v(k)$ , but this is beyond the scope of this investigation.

Figure 7.

Normalised spectrum $ E_v(k,t)/L_\delta$ as a function of $k L_\delta$ at different times. For cases with $k_m=0$ only.

We now explore some aspects of the time evolution of the spectrum $E_v(k,t)$ . An interesting observation is made by Ge et al. (Reference Ge, Rolland and Vassilicos2023) and Yoshimatsu & Ariki (Reference Yoshimatsu and Ariki2019), which shows that $E_\delta (k,t)$ evolves in a self-similar manner over a period of time. To examine this phenomenon in our simulations, we follow Ge et al. (Reference Ge, Rolland and Vassilicos2023) and define an integral length scale for the velocity perturbation by

(4.15) \begin{equation} L_\delta (t) = \frac {3\pi }{4K_\delta (t)} \int _0^\infty k^{-1} E_\delta (k,t) {\rm d}k = \frac {3\pi }{2} \int _0^\infty k^{-1} E_v(k,t) {\rm d}k. \end{equation}

Self-similar evolution takes place if $E_\delta (k,t)/K_\delta L_\delta$ is a function of $k L_\delta$ alone, independent of time. In terms of $E_v(k,t)$ , it implies that we have

(4.16) \begin{equation} E_v(k,t) = L_\delta (t) g(kL_\delta ) \end{equation}

for some function $g(\cdot )$ . Equation (4.16) implies that $E_v / L_\delta$ , when plotted against $k L_\delta$ , should collapse on a single curve. Figure 7 plots the results obtained from our data over a period of time spanning over $100 \tau$ . The immediate observation is that the curves mostly fall on each other. The agreement is the best for wavenumbers somewhat larger than the wavenumber where the curves peak. This feature is also observed by Ge et al. (Reference Ge, Rolland and Vassilicos2023). The discrepancies are larger at the two ends of the spectra. The peak of the normalised spectra is found approximately at $kL_\delta = 2$ , the same as by Ge et al. (Reference Ge, Rolland and Vassilicos2023). There are attempts to deduce analytically the slope of the spectra as $kL_\delta \to 0$ (Yoshimatsu & Ariki Reference Yoshimatsu and Ariki2019), but various values have been observed in DNS. For example, $k^4$ is found by Yoshimatsu & Ariki (Reference Yoshimatsu and Ariki2019) and slopes closer to $k^{3.3}$ are reported by Ge et al. (Reference Ge, Rolland and Vassilicos2023). For the cases in group R4, which have the largest Reynolds number in our simulations, the slope appears to scale with $k^{2.8}$ , as shown by the dash-double-dotted line.

Our observation broadly agrees with those of Ge et al. (Reference Ge, Rolland and Vassilicos2023) and Yoshimatsu & Ariki (Reference Yoshimatsu and Ariki2019). Note that the spectra shown in figure 7 are calculated from the long time limit of $\boldsymbol {v}$ . In contrast, Ge et al. (Reference Ge, Rolland and Vassilicos2023) and Yoshimatsu & Ariki (Reference Yoshimatsu and Ariki2019) observe self-similarity in an intermediate stage of the evolution of the velocity perturbation where no rescaling is applied to keep the perturbation small. The agreement between the results shows that the intermediate stage of evolution observed in the latter appears to be the same as the long-time asymptotic state of an infinitesimal perturbation.

Figure 8.

Normalised integral length scale $L_\delta /\eta$ for the perturbation velocity.

Figure 9.

Normalised instantaneous energy spectrum $ E_v (k,t)/ L_\delta$ as a function of $k L_\delta$ for different $k_m$ at different times. For cases in group R4.

It is natural to explore the relationship between the peak location of the self-similar spectrum and the peak location of $\langle E_v \rangle$ shown in figure 4(b). This can be inferred from figure 8, which shows that the ratio $L_\delta /\eta$ fluctuates around $10$ . Therefore, $kL_\delta \approx 2$ is equivalent to $k\eta \approx 0.2$ , which is the peak location of $\langle E_v \rangle (k)$ . This relationship lends further support to the conjecture that the peak wavenumber of $\langle E_v \rangle$ is a physically significant parameter for turbulence synchronisation. Another finding is that self-similarity is also observed for the spectra of the conditional LLVs, as shown in figure 9. For clarity, the curves for larger $k_m$ are shifted upwards by a factor of $10$ successively. Evidently, there is a very good agreement between the curves at different times, as required by self-similarity.

Figure 10.

Normalised instantaneous production spectrum $ P_v(k,t)/(\gamma L_\delta )$ as a function of $k L_\delta$ . For cases with $k_m = 0$ only.

The self-similar evolution can be examined quantitatively via the equation for the spectrum $E_v (k,t)$ , i.e. (2.27). Substituting (4.16) into (2.27), we obtain

(4.17) \begin{equation} g(\xi ) + \xi \frac { {\rm d} g}{{\rm d}\xi } = \frac {P_v(k,t)}{\dot {L}_\delta } - \frac {2\bar {\nu }}{L_\delta \dot {L}_\delta } \xi ^2 g(\xi ) - \frac {\gamma (t)L_\delta (t)}{\dot {L}_\delta } g(\xi ), \end{equation}

where $\xi \equiv k L_\delta$ and $\dot {L}_\delta$ is the time derivative of $L_\delta$ . Therefore, a fully self-similar solution (over all wavenumbers) is possible only if

(4.18) \begin{equation} \frac {P_v(k,t)}{\dot {L}_\delta } = h(\xi ), \qquad \frac {2\bar {\nu }}{L_\delta \dot {L}_\delta } = \alpha, \qquad \frac {\gamma (t) L_\delta (t)}{\dot {L}_\delta } = \beta, \end{equation}

where $h(\xi )$ is some function to be determined, and $\alpha$ and $\beta$ are constants. Figure 7 shows that $E_v(k,t)$ is mainly self-similar for intermediate wavenumbers. For these wavenumbers, we may drop the viscous effect, i.e. the second term on the right-hand side of (4.17). With that, only the first and the third equations in (4.18) are required for there to be a self-similar solution. The third equation establishes a relation between $\gamma (t)$ and $L_\delta (t)$ . It shows that $L_\delta$ grows exponentially if $\gamma (t)$ is a constant. The growth rate of $L_\delta$ is given by $\gamma /\beta$ and $L_\delta \sim \exp (\gamma t/\beta )$ . This regime appears to be that observed by Ge et al. (Reference Ge, Rolland and Vassilicos2023). However, figure 8 shows that $L_\delta (t)$ does not grow exponentially for our data (and $\gamma (t)$ is generally not a constant). Therefore, the self-similarity in our data belongs to a different regime, characterised more generally by the third equation in (4.18).

For the self-similar solution to exist, $P_v(k,t)$ must also have a self-similar form, as shown by the first equation in (4.18). Together with the third equation in (4.18), we obtain

(4.19) \begin{equation} P_v(k,t) = \dot {L}_\delta h(\xi ) = \beta ^{-1} \gamma (t) L_\delta (t) h(\xi ) \sim \gamma (t) L_\delta (t) h(\xi ). \end{equation}

We plot $P_v(k,t)/(\gamma (t) L_\delta (t))$ against $k L_\delta$ in figure 10. The agreement between the curves at different times is less satisfactory compared with that shown in figure 7, but the curves still largely fall on each other. The deviation from a clear self-similarity in $P_v$ could be due to the contamination from the two ends of $E_v(k,t)$ . Note that $E_v(k,t)$ is self-similar mainly in the mid-wavenumber range, which means $g(\xi )$ is well-defined only for a finite range of values for $\xi$ . Since $h(\xi ) \sim {\rm d} g/{\rm d}\xi$ according to (4.17), it is plausible that $h(\xi )$ is well defined over a narrower range of $\xi$ . This argument suggests that simulations covering a wider range of wavenumbers are required to ascertain whether strict self-similarity in $P_v(k,t)$ exists or not.

4.3. Production and dissipation: physical space analyses

Additional understanding of the production and dissipation of the LLV can be obtained with complementary analyses in the physical space. It has been known for a while that the spatial structures of the velocity perturbations (Nikitin Reference Nikitin2008, Reference Nikitin2018; Ge et al. Reference Ge, Rolland and Vassilicos2023) are non-trivial. Physical space analyses are well suited if one is interested in the impacts of these spatial structures.

In physical space analyses, it is more meaningful to express the production term in the intrinsic coordinates formed by the eigenvectors of the strain rate tensor. Let $s^+_{ij} \equiv \tau s_{ij}$ be the non-dimensional strain rate tensor, so that we can write $P = -\langle v_i v_j s^+_{ij} \rangle$ . We let $\lambda _\alpha \geqslant \lambda _\beta \geqslant \lambda _\gamma$ be the eigenvalues of $s^+_{ij}$ , with corresponding eigenvectors ${\boldsymbol {e}}_i$ $(i = \alpha, \beta, \gamma )$ . Due to incompressibility, we have $\lambda _\alpha + \lambda _\beta + \lambda _\gamma = 0$ . Thus, $\lambda _\alpha$ is always non-negative whereas $\lambda _\gamma$ is always non-positive. Letting $\theta _i$ be the angle between ${\boldsymbol {e}}_i$ and $\boldsymbol {v}$ , we may write

(4.20) \begin{equation} P = P_\alpha + P_\beta + P_\gamma, \end{equation}

with

(4.21) \begin{equation} P_\alpha = - 2\big\langle e\lambda _\alpha \cos ^2\theta _\alpha \big\rangle, \quad P_\beta = - 2\big\langle e \lambda _\beta \cos ^2\theta _\beta \big\rangle, \quad P_\gamma = - 2\big\langle e \lambda _\gamma \cos ^2\theta _\gamma \big\rangle . \end{equation}

Equation (4.21) captures the fact that the production hinges on the correlation among the eigenvalues of $s^+_{ij}$ , the alignment between $\boldsymbol {v}$ and the eigenvectors, and the magnitude of $\boldsymbol {v}$ .

Figure 11.

The $P_\alpha$ , $P_\beta$ and $P_\gamma$ components of the mean production $P$ as functions of $k_m\eta$ .

The production components $P_\alpha$ , $P_\beta$ and $P_\gamma$ as functions of $k_m \eta$ are shown in figure 11. One can observe that $P_\alpha$ , $P_\beta$ and $P_\gamma$ exhibit negligible dependence on $Re_\lambda$ . As expected, $P_\alpha \lt 0$ and $P_\gamma \gt 0$ . We also observe $P_\beta \lt 0$ with much smaller magnitudes, and $P_\gamma$ is the dominant one among the three terms. Of particular note is that $P_\alpha$ and $P_\beta$ are nearly independent of $k_m \eta$ , whereas $P_\gamma$ decreases as $k_m \eta$ increases. Thus, the change in $P$ with respect to $k_m \eta$ (as shown in figure 3) is predominantly due to the contribution from $P_\gamma$ . As a result, we focus on $P_\gamma$ only in what follows.

Figure 12 presents the p.d.f. of $\cos \theta _\gamma$ . As can be observed in figure 12(a), results for different Reynolds numbers are essentially the same. The p.d.f.s increase with $\cos \theta _\gamma$ , suggesting a strong tendency for the eigenvector ${\boldsymbol {e}}_\gamma$ to align with the vector $\boldsymbol {v}$ . There are noticeable differences between the p.d.f.s for different $k_m$ , as shown in figure 12(b). As $k_m$ increases, the preferential alignment between ${\boldsymbol {e}}_\gamma$ and $\boldsymbol {v}$ is weakened, manifested in the lower peaks. This behaviour clearly is one of the reasons why $P$ decreases with $k_m$ , as shown in figure 3. The p.d.f.s at other Reynolds numbers exhibit similar trends, thus not shown to avoid redundancy. The behaviours shown in figure 12(b) can be qualitatively understood from the characteristic length scales or wavenumbers of $s^+_{ij}$ and $v_i$ . The characteristic wavenumber of $s^+_{ij}$ can be estimated by the wavenumber where the dissipation spectrum of the base flow peaks, which is found to be approximately $0.15\eta ^{-1}$ . The characteristic wavenumber for $v_i$ can be estimated by $\max (0.2 \eta ^{-1}, k_m)$ , with $0.2 \eta ^{-1}$ being the peak wavenumber for $\langle E_v \rangle$ when it is bigger than $k_m$ . Thus, as $k_m$ increases, the mismatch between the two characteristic wavenumbers tends to increases, which tends to weaken the correlation between $s^+_{ij}$ and $v_i$ , hence the alignment in figure 12(b).

Figure 12.

The p.d.f.s of $\cos \theta _\gamma$ . (a) Cases from all groups with $k_m = 0$ . (b) Cases in group R4 with different $k_m$ .

Incidentally, the preferential alignment discussed above is reminiscent of the behaviours of the gradient of a passive scalar in isotropic turbulence, which also tends to align with ${\boldsymbol {e}}_\gamma$ of the strain rate tensor (Ashurst et al. Reference Ashurst, Chen and Rogers1987). However, the statistics of $\boldsymbol {v}$ are different from those of a passive scalar on many aspects, as we can see from the statistics of $B_{ij}$ which will be discussed later.

The impacts of the correlation between the perturbation and the strain rate tensor can be explored with suitable conditional statistics. Note that

(4.22) \begin{equation} P_\gamma = -2 \int \langle e \cos ^2 \theta _\gamma \vert \lambda _\gamma \rangle \lambda _\gamma f_\gamma (\lambda _\gamma ) {\rm d} \lambda _\gamma, \end{equation}

where $f_\gamma (\lambda _\gamma )$ is the p.d.f. of the eigenvalue $\lambda _\gamma$ . Therefore, how the correlation between $\lambda _\gamma$ , the alignment and $e$ contributes to $P_\gamma$ can be understood from the conditional average $\langle e \cos ^2 \theta _\gamma \vert \lambda _\gamma \rangle$ . Our tests show that $\langle e \cos ^2 \theta _\gamma \vert \lambda _\gamma \rangle$ tends to be smaller than $\langle e \vert \lambda _\gamma \rangle$ due to the factor $\cos ^2 \theta _\gamma$ , but the two distributions display similar shapes. To keep the discussion succinct, we consider only $\langle e \vert \lambda _\gamma \rangle$ , which is shown in figure 13. Figure 13(a) shows that $ \langle e \vert \lambda _\gamma \rangle$ changes with the Reynolds number quite significantly, in contrast to the alignment trend shown in figure 12. The impact of the Reynolds numbers is especially strong for large $\vert \lambda _\gamma \vert$ , where the conditional average generally is larger for larger Reynolds numbers. This behaviour is likely due to the fact that the probability for strong strain rate increases with the Reynolds number due to the intermittency effects. Figure 13( b) plots $\langle e \vert \lambda _\gamma \rangle$ for different $k_m$ with a fixed Reynolds number. The conditional average increases with $\vert \lambda _\gamma \vert$ for all $k_m$ , but it is generally smaller for larger $k_m$ . This behaviour is another factor by which $P$ decreases as $k_m$ increases, in addition to the weakened alignment shown in figure 12(b). The reduction in $\langle e \vert \lambda _\gamma \rangle$ as $k_m$ increases may also be attributed to the mismatch between the characteristic length scales of $s^+_{ij}$ and $v_i$ . Overall, figure 13 shows that the perturbation tends to be stronger at regions with stronger strain rate (shown with larger $\vert \lambda _\gamma \vert$ ), though it is reduced as $k_m$ increases.

Figure 13.

Conditional average $ \langle e \vert \lambda _\gamma \rangle$ as a function of $\lambda _\gamma$ . (a) Cases from all groups with $k_m = 0$ . (b) Cases in group R4 with different $k_m$ .

Figure 14.

Conditional average $ \langle e \vert {\omega }_i^+{\omega }_i^+ \rangle$ . (a) Cases from all groups with $k_m = 0$ . (b) Cases in group R4 with different $k_m$ .

The descriptions of the velocity perturbation can be further elaborated by considering the correlation between $e$ and the base flow vorticity. Let ${\boldsymbol {\omega }} = \nabla \times {\boldsymbol {u}}$ be the vorticity of the base flow and $\boldsymbol {\omega }^+ = \tau {\boldsymbol {\omega }}$ be the non-dimensionalised version of $\boldsymbol {\omega }$ . Figure 14 shows the conditional average $\langle e \vert {\omega }_i^+{\omega }_i^+ \rangle$ , which characterises the correlation between the magnitude of the perturbation and the vorticity of the base flow. When $k_m = 0$ , figure 14 shows that the conditional average increases with the magnitude of base flow vorticity and it depends only weakly on the Reynolds number. When $k_m$ is increased, the dependence of $ \langle e \vert {\omega }_i^+ {\omega }_i^+\rangle$ on ${\omega }_i^+ {\omega }_i^+$ is weakened. As a result, the conditional average increases with $k_m$ for smaller ${\omega }_i^+ {\omega }_i^+$ and decreases with $k_m$ for larger ${\omega }_i^+ {\omega }_i^+$ . Therefore, perturbations associated with regions of strong vorticity in the base flow are stronger on average, and this trend tends to be weakened as the coupling wavenumber $k_m$ increases. The conditional average $\langle e\vert {\omega }_i^+{\omega }_i^+\rangle$ depends on the Reynolds number and $k_m$ in the same way as $\langle e \vert \lambda _\gamma \rangle$ , thus its behaviours can be explained in a similar way. The correlation between $e$ and strong strain rate and strong vorticity is also observed in channel flows to some extent (Nikitin Reference Nikitin2018).

Figure 15.

(a) P.d.f.s $f_e(e)$ from all groups with $k_m=0$ with the inset showing the variation of the p.d.f.s with the Reynolds number near $e=0$ . (b) Conditional average $\langle \lambda _\gamma \vert e \rangle$ for group R4 with different $k_m$ .

To understand how the fluctuations in the perturbation velocity contribute to the mean production term, we may write the mean production as the weighted integral of the average conditioned on given $e$ , i.e.

(4.23) \begin{equation} P_\gamma = -2 \int \big\langle \lambda _\gamma \cos ^2 \theta _\gamma \vert e \big\rangle e f_e (e) {\rm d}e, \end{equation}

where $f_e(e)$ is the p.d.f. of $e$ . Figure 15(a) plots $f_e(e)$ for different Reynolds numbers with $k_m=0$ only. The p.d.f.s display very elongated tails, showing high probabilities for large fluctuations in $e$ . The tail is very slightly fatter for higher Reynolds numbers. The peaks of the p.d.f.s are found at small $e$ values and they are slightly sharper for higher Reynolds numbers. The distributions indicate that the spatial distribution of the perturbation velocity is highly intermittent, with small fluctuations covering large part of the spatial domain and strong fluctuations observed in localised spots. The results for $\langle \lambda _\gamma \vert e\rangle$ are shown in figure 15(b). The magnitude of the conditional average increases with $e$ and $k_m$ . These behaviours are consistent with the results for $\langle e \vert \lambda _\gamma \rangle$ . Given the highly intermittent nature of the distribution of $e$ , one might ask how important are the large fluctuations to the mean production. Though its figure is omitted, one can readily see that the product $-\langle \lambda _\gamma \cos ^2\theta _\gamma \vert e \rangle e f_e(e ) \approx - \langle \lambda _\gamma \vert e \rangle e f_e(e)$ , as a function of $e$ , peaks at an intermediate value of $e$ . Therefore, the main contribution to the mean production does not come from the very large fluctuations.

We now explore the behaviours of the dissipation term ${\mathcal D}_e$ in the physical space. Equation (2.22) shows that ${\mathcal N}_e$ is one of the main mechanisms that determines the dissipation rate ${\mathcal D}_e$ , and it will be the focus below. We let

(4.24) \begin{equation} N_{es} = \frac {\tau \langle {\mathcal N}_{es}\rangle }{\langle {\mathcal D}_e \rangle }, \qquad N_{eo} = \frac {\tau \langle {\mathcal N}_{eo}\rangle }{\langle {\mathcal D}_e \rangle }, \end{equation}

which are both dimensionless (cf. (2.22)). Equivalently, we may introduce a length scale $(\bar {\nu }/\langle {\mathcal D}_e \rangle )^{1/2}$ , and then define non-dimensional perturbation strain rate and perturbation vorticity, denoted by $s_{ij}^-$ and ${\omega }^-_{i}$ , respectively, with

(4.25) \begin{equation} s^-_{ij} =\frac {\bar {\nu }^{1/2} s^v_{ij}}{\langle {\mathcal D}_e \rangle ^{1/2}}, \qquad {\omega }^-_{i} = \frac {\bar {\nu }^{1/2} {\omega }^v_{i}}{\langle {\mathcal D}_e \rangle ^{1/2}}. \end{equation}

Recalling that $\boldsymbol {v}$ is dimensionless, therefore, the dimension of $s^v_{ij}$ and ${\omega }^v_i$ is that of the reciprocal of length, so that $s^-_{ij}$ and ${\omega }^-_i$ are dimensionless. As a consequence, we obtain

(4.26) \begin{equation} N_{es} = -4 \langle s^+_{ij} s^-_{jk} s^-_{ki} \rangle, \quad N_{eo} = \langle s^+_{ij} {\omega }^-_i {\omega }^-_j\rangle . \end{equation}

In terms of $N_{es}$ and $N_{eo}$ , we may re-write (2.24). Assuming the correlation between $\gamma$ and ${\mathcal D}_e$ is negligible, we obtain $\langle \gamma {\mathcal D}_e \rangle \approx \langle \gamma \rangle \langle {\mathcal D}_e \rangle = \lambda \langle {\mathcal D}_e \rangle$ . Therefore, (2.24) becomes

(4.27) \begin{equation} 2 \Lambda \approx N_{es} + N_{eo} + r_e, \end{equation}

where

(4.28) \begin{equation} r_e \equiv - \frac {2 \bar {\nu } \tau \langle s^v_{ij} \partial ^2_{ij} p_e \rangle }{\langle {\mathcal D}_e \rangle } - \frac {2 \tau \langle \bar {\nu } \nabla B_{ij} \cdot \bar {\nu } \nabla B_{ij} \rangle }{ \langle {\mathcal D}_e \rangle } + \frac {2 \tau \bar {\nu } \langle B_{ij} \nabla \cdot ( \boldsymbol {v} A_{ij})\rangle }{ \langle {\mathcal D}_e \rangle } \end{equation}

is considered a ‘residual’ term. Therefore, the values of $N_{eo}$ and $N_{es}$ can be compared with $2\Lambda$ to gauge their contributions.

The expression of $N_{eo}$ is similar to that of $P$ in form with $v_i$ replaced by ${\omega }^-_i$ . It is also similar in form to the vortex stretching term $s_{ij}^+ {\omega }^+_i {\omega }^+_j$ for the enstrophy of the base flow. Using the eigen-frame defined by the eigenvectors of $s^+_{ij}$ introduced previously, we may write

(4.29) \begin{equation} N_{eo} = N_{eo\alpha } + N_{eo\beta } + N_{eo\gamma }, \end{equation}

with

(4.30) \begin{equation} N_{eo\alpha } = \langle \lambda _\alpha \vert {\boldsymbol {\omega }}^{-}\vert ^2 \cos ^2\theta ^o_\alpha \rangle, \, N_{eo\beta } = \langle \lambda _\beta \vert {\boldsymbol {\omega }}^{-}\vert ^2 \cos ^2\theta ^o_\beta \rangle, \, N_{eo\gamma } = \langle \lambda _\gamma \vert {\boldsymbol {\omega }}^{-} \vert ^2 \cos ^2 \theta ^o_\gamma \rangle, \end{equation}

where $\theta ^o_i$ ( $i=\alpha, \beta, \gamma$ ) denotes the angle between ${\boldsymbol {\omega }}^-$ and ${\boldsymbol {e}}_i$ . Equations (4.29) and (4.30) are similar to those for $P$ . Similarly, $N_{es}$ has an expression in terms of the eigenvalues and eigenvectors of $s^+_{ij}$ and $s^-_{ij}$ . Let ${\boldsymbol {e}}_\ell ^-$ ( $\ell =\alpha, \beta, \gamma$ ) be the eigenvectors of $s^-_{ij}$ , with corresponding eigenvalues $\lambda ^-_\ell$ . We follow the same tradition where $\lambda ^-_\alpha \geqslant \lambda ^-_\beta \geqslant \lambda ^-_\gamma$ . Letting $\theta ^s_{ij}$ be the angle between ${\boldsymbol {e}}^-_i$ and ${\boldsymbol {e}}_j$ , we obtain

(4.31) \begin{align} N_{es} &= -\,4 \Big[\Big\langle (\lambda _\alpha ^-)^2 \lambda _\alpha \cos ^2 \theta ^s_{\alpha \alpha }\Big \rangle + \Big\langle (\lambda _\alpha ^-)^2 \lambda _\beta \cos ^2 \theta ^s_{\alpha \beta }\Big\rangle + \Big\langle (\lambda _\alpha ^-)^2 \lambda _\gamma \cos ^2 \theta ^s_{\alpha \gamma } \Big\rangle\notag \\ &\quad+ \Big \langle (\lambda _\beta ^-)^2 \lambda _\alpha \cos ^2 \theta ^s_{\beta \alpha }\Big\rangle + \Big\langle (\lambda _\beta ^-)^2 \lambda _\beta \cos ^2 \theta ^s_{\beta \beta }\Big\rangle + \Big\langle (\lambda _\beta ^-)^2 \lambda _\gamma \cos ^2 \theta ^s_{\beta \gamma }\Big\rangle \notag \\ &\quad+ \Big\langle (\lambda _\gamma ^-)^2 \lambda _\alpha \cos ^2 \theta ^s_{\gamma \alpha }\Big\rangle + \Big\langle (\lambda _\gamma ^-)^2 \lambda _\beta \cos ^2 \theta ^s_{\gamma \beta }\Big\rangle + \Big\langle (\lambda _\gamma ^-)^2 \lambda _\gamma \cos ^2\theta ^s_{\gamma \gamma }\Big\rangle \Big] \end{align}

The equation provides a decomposition of $N_{es}$ into contributions associated with different eigenvalues of the two tensors and makes explicit how the relative orientation of the eigenvectors affects $N_{es}$ . We will use $N_{es\alpha \alpha }$ , $N_{es\alpha \beta }$ , …, $N_{es\gamma \gamma }$ to denote the nine components on the right-hand side of the equation. We also let $N_{es\alpha } = N_{es\alpha \alpha } + N_{es \alpha \beta } + N_{es\alpha \gamma }$ represent the sum of the contributions involving the eigenvalue $\lambda ^-_\alpha$ , $N_{es\beta } = N_{es\beta \alpha } + N_{es\beta \beta } + N_{es\beta \gamma }$ represent the sum of the contributions involving $\lambda ^-_\beta$ and $N_{es\gamma } = N_{es\gamma \alpha } + N_{es\gamma \beta } + N_{es \gamma \gamma }$ represent the sum of the contributions involving $\lambda ^-_\gamma$ . Obviously, we have $N_{es} = N_{es\alpha } + N_{es\beta } + N_{es\gamma }$ .

Figure 16.

(a) Plots of $N_{eo}$ and its components as functions of $k_m\eta$ . (b) Plots of $N_{es}$ and the contributions $N_{es\alpha }$ , $N_{es\beta }$ and $N_{es\gamma }$ . (c) Nine components of $N_{es}$ . For the cases in group R4.

The data for $N_{eo\alpha }$ , $N_{eo \beta }$ and $N_{eo\gamma }$ are plotted in figure 16(a). The magnitudes of these values decrease only slightly as $k_m\eta$ increases. Recalling the normalisation shown in (4.24) and the fact that $\langle {\mathcal D}_e \rangle$ increases as $k_m\eta$ increases, the conclusion one may draw is thus that the magnitudes of the non-normalised versions of $N_{eo\alpha }$ , $N_{eo\beta }$ and $N_{eo\gamma }$ all increase with $k_m \eta$ , but at rates that are slightly smaller than that of $\langle {\mathcal D}_e \rangle$ , so that the magnitudes of $N_{eo\alpha }$ , $N_{eo\beta }$ and $N_{eo\gamma }$ decrease slightly as $k_m\eta$ increases. The figure shows that $N_{eo\beta }$ and $N_{eo\gamma }$ have opposite signs with similar magnitudes. Consequently, $N_{eo}$ is only slightly different from $N_{eo\alpha }$ .

Figure 16(b) shows $N_{es\alpha }$ , $N_{es\beta }$ , $N_{es\gamma }$ , together with $N_{es}$ . In terms of the dependence on $k_m\eta$ , these parameters all behaviour similarly to $N_{eo}$ , i.e. their magnitudes all decrease with $k_m$ , but only weakly. Among the three components, $N_{es\beta }$ is the smallest and essentially negligible. Both $N_{es\alpha }$ and $N_{es\gamma }$ are much larger, and they appear to be almost the same as each other. The breakdown into the nine components is given in figure 16(c). We can see that $N_{es\beta \alpha }$ , $N_{es\beta \beta }$ , $N_{es\beta \gamma }$ , $N_{es\gamma \beta }$ and $N_{es\alpha \beta }$ are all quite small. The largest contributions come from $N_{es\alpha \gamma }$ and $N_{es\gamma \gamma }$ , which are positive by definition, and appear to be identical. The contributions from $N_{es\alpha \alpha }$ and $N_{es\gamma \alpha }$ are also significant though somewhat smaller than those from $N_{es\alpha \gamma }$ and $N_{es\gamma \gamma }$ . They also appear to be identical.

The results given in the panel (c) shows that the close agreement between $N_{es\alpha }$ and $N_{es\gamma }$ (panelb) is a consequence of the close agreement between $N_{es\alpha \ell }$ and $N_{es\gamma \ell }$ ( $\ell = \alpha, \beta, \gamma$ ). The agreement between the latter two is, in fact, a mathematical consequence of the linearity of (2.19). The linearity of (2.19) dictates that $s^-_{ij}$ and $-s^-_{ij}$ must have the same statistics, which implies that the largest eigenvalues of the two, $\lambda ^-_\alpha$ and $-\lambda ^-_\gamma$ , respectively, should have the same statistics too. As a result, $N_{es\alpha \ell } = N_{es\gamma \ell }$ exactly for all $k_m$ , which is reflected in the figure.

We will not discuss the residual term $r_e$ in detail to keep the scope of this investigation manageable. Nevertheless, we may use (4.27) to obtain an estimate of its impact by comparing the values of $N_{eo}$ and $N_{es}$ with $ 2\Lambda$ . Recall that, according to figure 2(b ), $2 \Lambda \approx 0.24$ when $k_m = 0$ . However, figure 16 shows that $N_e = N_{eo} + N_{es} \approx 0.7$ for $k_m=0$ . Therefore, the residual term $r_e$ has a significant contribution and appears to be acting to counter the effects of $N_e$ . Furthermore, $\Lambda$ decreases from $0.12$ to $-0.20$ as $k_m \eta$ increases according to figure 2 (for cases in group R4). Though $N_e$ also decreases with $k_m$ , the change is not large enough to account for the change in $\Lambda$ , which shows that the dependence of $r_e$ on $k_m$ also plays a role.

Equations (4.30) and (4.31) also suggest that the relative orientations between ${\boldsymbol {\omega }}^-$ and $s^+_{ij}$ or between $s^-_{ij}$ and $s^+_{ij}$ might impact the values of $N_{eo}$ and $N_{es}$ . We thus look into relevant results, for $N_{eo}$ to begin with. The alignment between ${\boldsymbol {\omega }}^-$ and $s^+_{ij}$ can be characterised by the angles $\theta ^o_\ell$ ( $\ell = \alpha, \beta, \gamma$ ) introduced previously. However, since this problem has not been investigated before, we opt for a more complete description based on the polar angle $\theta ^o_P$ and the azimuthal angle $\psi ^o_P$ that the vector ${\boldsymbol {\omega }}^-$ make in the eigen-frame formed by the eigenvectors of $s^+_{ij}$ . The definitions of the two angles are illustrated in figure 17. Specifically, $\theta ^o_P$ is the angle between ${\boldsymbol {\omega }}^-$ and the polar direction ${\boldsymbol {e}}_\alpha$ , and $\psi ^o_P$ is the angle between ${\boldsymbol {e}}_\beta$ and the projection of ${\boldsymbol {\omega }}^-$ on the equatorial plane. The relations between $(\theta ^o_P, \psi ^o_P)$ and $\theta ^o_i$ can be derived readily. The joint p.d.f. of $\cos \theta ^o_P$ and $\psi ^o_P$ is shown in figure 18. The joint p.d.f. has a very sharp peak at the origin, i.e. at $\theta ^o_P = 90^\circ$ and $\psi ^o_P = 0$ . Thus, ${\boldsymbol {\omega }}^-$ tends to very strongly align with the intermediate eigenvector ${\boldsymbol {e}}_\beta$ of $s^+_{ij}$ (cf. figure 17). This geometrical feature appears to have not be reported before.

Figure 17.

Definitions of $\theta ^o_P$ and $\psi ^o_P$ .

Figure 18.

Alignment between ${\boldsymbol {\omega }}^-$ and the eigenvectors of $s^+_{ij}$ , for case R $4$ K $0$ .

The alignment between $s^+_{ij}$ and $s^-_{ij}$ can be described by the polar angles that each individual eigenvector of $s^-_{ij}$ makes in the eigen-frame of $s^+_{ij}$ . The polar angles are defined in the same way shown in figure 17, with ${\boldsymbol {\omega }}^-$ replaced by one of the eigenvectors, such as ${\boldsymbol {e}}^-_\alpha$ . We use $\theta ^s_P$ and $\psi ^s_P$ to denote the angles. Figure 19 plots the joint p.d.f.s of $\cos \theta ^s_P$ and $\psi ^s_P$ for the three eigenvectors, ${\boldsymbol {e}}^-_\alpha$ , ${\boldsymbol {e}}^-_\beta$ and ${\boldsymbol {e}}^-_\gamma$ , in panels (a), (b) and (c), respectively. Distinct peaks can be identified for all three distributions, though the peaks are not as sharp as in, e.g. figure 18. Panel (a) shows that ${\boldsymbol {e}}^-_\alpha$ displays a bi-modal behaviour, with the alignment switching between $(\theta ^s_P, \psi ^s_P) = (90^\circ, 90^\circ )$ and $\theta ^s_P = 0^\circ$ (the value of $\psi ^s_P$ is not defined when $\theta ^s_P = 0$ ). In the first configuration, ${\boldsymbol {e}}^-_\alpha$ aligns with ${\boldsymbol {e}}_\gamma$ , whereas in the second configuration, ${\boldsymbol {e}}^-_\alpha$ aligns with ${\boldsymbol {e}}_\alpha$ . The eigenvector ${\boldsymbol {e}}^-_\beta$ , as shown in panel (b), tends to align with ${\boldsymbol {e}}_\beta$ , since the p.d.f. peaks at $(\theta ^s_P, \psi ^s_P) = (90^\circ, 0)$ . Panel (c) shows the joint p.d.f. for ${\boldsymbol {e}}^-_\gamma$ . Due to the linearity of the equation for $\boldsymbol {v}$ , it should be exactly the same as that for ${\boldsymbol {e}}^-_\alpha$ shown in panel (a). Due to statistical fluctuations, the two joint p.d.f.s are not exactly the same, but they are very close, as expected. For example, they display exactly same peak locations.

Figure 19.

Joint p.d.f.s of $\cos \theta ^s_P$ and $\psi ^s_P$ for (a) ${\boldsymbol {e}}^-_\alpha$ , (b) ${\boldsymbol {e}}^-_\beta$ and (c) ${\boldsymbol {e}}^-_\gamma$ . For case R $4$ K $0$ .