The conditional Lyapunov exponents and synchronisation of rotating turbulent flows

We present some results in this subsection to illustrate the basic features of the flow fields. The energy spectra normalised by Kolmogorov parameters are shown in Fig. 1. For the flows driven by Kolmogorov forcing shown in the left panel, the normalised spectra collapse onto a single curve except for the few lowest wavenumbers. At the lowest wavenumbers, the spectra increase with the rotation rate, which shows increased energetics for the large scales, consistent with our understanding of rotating turbulence.

The Reynolds number for the flows is relatively small so no clear inertial range can be identified. Nevertheless, the spectra appear to be consistent with the $k^{-2}$ scaling law which has been reported in previous research (Yeung & Zhou, 1998; Dallas & Tobias, 2016).

For the flows driven by constant power forcing, similar behaviours are observed for lower rotation rates, as shown in the right panel. However, for $\Omega=1$, the spectra have steeper slopes in the mid-wavenumber range, and they appear to be more consistent with the $k^{-3}$ power law. The spectra in the dissipation range also appear to drop off at a faster rate. The contrast between the left and right panels shows that the forcing terms can lead to significant quantitative differences in the flows.

In both flows, energy pile-up is observed at the lowest wavenumber end of the spectra, and the pile-up increases slightly with the rotation rate. The pile-up is an indication of the emergence of large scale columnar vortices, which is a common feature of rotating turbulence. Columnar vortices are indeed visually observable in our simulations with the larger rotation rates, which are illustrated in Fig. 2 for two simulations with $\Omega=1$. The figure shows a snapshot of the distribution of $|{\bm{\omega}}|$ on three horizontal cross sections of the flow domain, where ${\bm{\omega}}\equiv\nabla\times{\bm{u}}$ is the vorticity. A columnar vortex is visible at the left corner in both flows shown in the two panels. The left panel shows a simulation with a smaller Reynolds number. In this case, the diameter of the columnar vortex is roughly half of the size of the domain. For the flow with a larger Reynolds number (right panel), the background vorticity is stronger, and the columnar vortex appears to be slightly smaller in size but it is still clearly visible. We will not show the results for other rotation rates, but we can confirm that columnar vortices are also quite prevalent for $\Omega=0.5$, while they are rare for $\Omega=0.1$.

The probability density function (PDF) of the vorticity component along the rotation axis is also of interest because it is well known that the PDF displays a positive skewness (Bartello et al., 1994; Morize et al., 2005) in rotating turbulence, due to the prevalence of cyclonic vortices over the anti-cyclonic ones. The skewness emerges as rotation is introduced, peaks at an intermediate rotation rate, and then decreases when the rotation rate further increases as the flow is two-dimensionalized under strong rotation. The PDFs for our simulations are plotted in Fig. 3. The PDFs are indeed skewed towards the positive values with the corresponding skewness given in the parentheses. For flows driven by constant power forcing with $N=128$, the skewness for $\Omega=1$ is slightly smaller than that for $\Omega=0.5$. In other cases, the skewness increases with the rotation rate. These PDFs show, from another angle, that the effects of rotation are clearly significant.

Table 1 shows that, compared with the flows driven by constant power forcing, those driven by Kolmogorov forcing tend to have larger micro-scale Rossby numbers $Ro_{k}$ for a given rotation rate $\Omega$. In order to obtain even smaller $Ro_{k}$ for the latter flows, we computed a few test cases with $\Omega=5$, and found that the flows are strongly two-dimensionalised at this rotation rate. Let $E_{2D}(t)$ be the kinetic energy in the two dimensional Fourier modes with $k_{z}=0$, and $E(t)$ be the total kinetic energy, i.e.,

The results for $E_{2D}(t)$ and $E(t)$ for the flows with $\Omega=5$ (i.e., cases F1N128$\Omega$5 and F1N192$\Omega$5) are shown in the left panel of Fig. 4. As a comparison, the results for $\Omega=1$ are shown in the middle panel. It can be observed that, for $\Omega=5$, both $E(t)$ and $E_{2D}(t)$ are an order of magnitude higher than for $\Omega=1$, and almost all energy is contained in the two dimensional modes as $E_{2D}(t)$ deviates from $E(t)$ only slightly. There are regular periods of time in which $E_{2D}(t)$ is indistinguishable from $E(t)$. These behaviours suggest that, at $\Omega=5$, the flows are quasi-two-dimensionalised with large-scale, 2D columnar vortices, where instability sets in periodically which leads to temporary small deviation between $E_{2D}(t)$ and $E(t)$. A detailed discussion of this process can be found in Alexakis (2015). The energy spectra for the flow with $\Omega=5$ at various times are shown in the right panel of Fig. 4 in green lines, together with those for $\Omega=1$ with both forcing terms (shown in black or red). The high wavenumber ends of the spectra swing violently over time, in a range spanning five orders of magnitude. Though oscillations are also seen in the spectra for the flows with $\Omega=1$, the amplitude is much smaller.

To summarize, the results in this subsection show that for $\Omega=0.1$, $0.5$ and $1$, the flows are still predominantly turbulent while displaying strong effects of rotation. The flows where $\Omega=5$, on the other hand, appear to be mostly two-dimensionalised and only display weakly turbulent behaviours. We will limit our interests in the synchronisation of flows where turbulence dominates. Therefore we will focus on the first three rotation rates, and the cases with $\Omega=5$ will not be discussed further in what follows.

We now look into the synchronisation of the flows. To obtain smoother results, the data shown in this subsection are the averages of five realisations.

Fig. 5 shows the decay of the synchronisation error $\Delta(t)/\Delta(0)$ for different $k_{m}$ and $\Omega$ with Kolmogorov forcing. The top-left, top-right, and bottom-left panels correspond to three different Reynolds numbers. There are three common trends across all cases included in these three panels. Firstly, the error decays exponentially when $k_{m}$ is sufficiently large. Secondly, the decay rate increases with $k_{m}$. Thirdly, the error decays only when $k_{m}$ is greater than some threshold value, referred to as $k_{c}$, and $k_{c}$ clearly is different in different cases. For $k_{m}$ close to but still greater than $k_{c}$, the error still decays over time, but the rate of decay fluctuates, so exponential functions do not always provide a good fit.

Comparison across the above three panels in Fig. 5 shows that the decay rate of the error displays the known dependence on the Reynolds number, namely, everything else being equal, the decay rate decreases as the Reynolds number increases. This trend is illustrated in the bottom-right panel with selected cases with $\Omega=0.5$ and $k_{m}=9$. As this effect has been reported multiple times in previous research, we will not delve too much into it. For the same reason, we consider only the cases with $N=128$ and $N=192$ for flows with constant power forcing.

More pertinent to our objectives is the observation that rotation has a strong effect on the decay rate. Fig. 5 shows that, for the same $k_{m}$, the decay rate decreases with $\Omega$. The same trend is observed for different Reynolds numbers, as is shown in the first three panels of the figure.

The results corresponding to constant power forcing are plotted in Fig. 6. Not surprisingly, $\Delta(t)$ decays exponentially for sufficiently large $k_{m}$. Moreover, the dependence of the decay rate on $k_{m}$ and $Re_{\lambda}$ is qualitatively similar to that which is observed in Fig. 5. However, interestingly, the dependence on rotation is significantly different. The black lines in the left panel of Fig. 6 illustrate the difference clearly. The three black lines correspond to same $k_{m}$ but three different rotation rates. While the decay rates for $\Omega=0.1$ and $0.5$ show no clear differences, the decay rate for $\Omega=1$ is clearly larger. That is, in this case, it appears that the decay rate for $\Delta(t)$ increases with rotation. The same trend is seen in the right panel of the figure, which is for flows with a larger Reynolds number. This observation is opposite to the trend we observe in the cases with Kolmogorov forcing (c.f. Fig. 5), where the decay rate for the same $k_{m}$ is found to decrease with rotation. The difference in the results for the two forcing terms has not been reported before.

The synchronisability of the slaved flow is related to the conditional Lyapunov exponents. We calculate the CLEs $\lambda(k_{m})$ as well as the local CLEs $\gamma(k_{m},t)$ using the algorithm outlined in Section 2. The results are presented in terms of the non-dimensionalised CLEs $\Lambda$ and the non-dimensionalised local CLEs $\Gamma$, which are defined as

$\Gamma$ is time dependent and fluctuates over time. Without showing the time sequences, we note that, after a period of transience, $\Gamma$ stabilizes and fluctuates around a constant value. The magnitude of the fluctuations appears to increase with rotation, but decreases as $k_{m}$ increases. We will quantify some of these behaviours in what follows, starting with $\Lambda$, which is the average of $\Gamma$ in the stationary stage.

Fig. 7 is shown first to establish the relationship between the decay rate of $\Delta(t)$ and the CLE $\Lambda$. Shown with symbols in the figure are $\Delta(t)/\Delta(0)$ for a number of cases already discussed in Figs. 5 and 6. The lines without symbols represent functions $\exp(\Lambda t/\tau_{k})$, where $\Lambda$ is the CLE for the corresponding flow. Some small discrepancies are seen between the two, which we attribute to statistical uncertainty in the data. We note that the discrepancies are in-line with those found in previous research (e.g., Nikolaidis & Ioannou (2022)). The overall agreement between the two shows that in most cases the error decays exponentially and the decay rate $\alpha$ equals $\Lambda$. For the case shown with the black line and solid triangles, $\Delta(t)$ does not decay exponentially. However, it undulates mildly around the exponential function in such a way that $\Lambda$ appears to capture the long time mean decay rate. Overall, we may conclude that the decay rate of $\Delta(t)$ is equal to $\Lambda$, and the synchronisation between two flows can be fully characterised by $\Lambda$.

As a side note, we note that larger discrepancy is observed for case F1N128$\Omega$01K7 than for case F2N192$\Omega$1K9. This observation appears counter intuitive at first sight, since $\Omega$ is larger in the latter case which should lead to larger fluctuation in $\Gamma$ hence larger statistical error in $\Lambda$ (or the corresponding decay rate $\alpha$). However, there is another difference between these two cases, which is that case F2N192$\Omega$1K9 is computed with a larger $k_{m}$. As the fluctuation in $\Gamma$ is smaller for larger $k_{m}$, it is possible that the statistical discrepancy in case F2N192$\Omega$1K9 is smaller despite the fact that it is computed with a larger $\Omega$.

We now focus on the results for $\Lambda$. The dependence of $\Lambda$ on the rotation rate $\Omega$ and the coupling wavenumber $k_{m}$ is shown in Fig. 8, including cases with $k_{m}=0$ where $\Lambda$ represents the unconditional Lyapunov exponent. The left panel presents the cases with Kolmogorov forcing. In these cases, $\Lambda$ always increases with $\Omega$, and $\Lambda$ increases with $\Omega$ quicker for larger $k_{m}$. The blue lines, which correspond to $k_{m}=5$ for $N=128$ and $k_{m}=7$ for $N=192$, are particularly instructive. In these cases $\Lambda$ increases from a negative value to a positive one as $\Omega$ increases from $0.1$ to $1$. Therefore, the two flows synchronise when $\Omega=0.1$, but they do not when $\Omega=1$, which emphatically shows that rotation makes the flows more difficult to synchronise when the flow is driven by Kolmogorov forcing. However, the observation is different for the flows maintained by constant power forcing, which is shown in the right panel of Fig. 8. In fact, the trend is reversed in this case: here $\Lambda$ decreases as $\Omega$ increases, so the flow is easier to synchronise as rotation is increased. Also, the unconditional Lyapunov exponent appears more sensitive to the rotation rate.

Another observation we can make from Fig. 8 is that $\Lambda$ decreases with $k_{m}$, which can be seen by comparing different curves in the same panel. This trend is further investigated by plotting $\Lambda$ as a function of $k_{m}\eta$, which is given in Fig. 9. We first note the values of $\Lambda$ at $k_{m}=0$ for $\Omega=0.1$. As $\Omega$ is relatively small, one expects $\Lambda$ to be close to the value found in non-rotating turbulence. Fig. 9 shows that $\Lambda$ in this case is around $0.1$, though it weakly depends on the Reynolds number as well as the forcing term. This value is indeed close to those found previously for non-rotating turbulence (Boffetta & Musacchio, 2017).

Fig. 9 shows that, for cases with Kolmogorov forcing, $\Lambda$ decreases as $k_{m}\eta$ increases. More interestingly, the curves corresponding to different cases collapse on each other approximately. The one for $N=192$ and $\Omega=1.0$ is slightly larger than the rest. Nevertheless, overall, as a function of $k_{m}\eta$, $\Lambda$ depends on rotation only weakly. Note that this observation does not contradict with the results in Fig. 8, as the values of $k_{m}$ in the latter are not non-dimensionalised by $\eta$ and $\eta$ is different for different $\Omega$.

For the cases with constant power forcing, the right panel in Fig. 9 shows that $\Lambda$ decreases with $k_{m}\eta$ in a similar manner. However, the curves corresponding to different $\Omega$ do not collapse well. In fact, $\Lambda(k_{m}\eta)$ tends to decrease with $\Omega$, in particular for stronger rotations.

The threshold wavenumber $k_{c}$ where $\Lambda$ is zero is of particular interests, as it is the value of $k_{m}$ for which synchronisation fails. The values of $k_{c}$ can be found from Fig. 9, as they are the values of $k_{m}$ where the curves cross the horizontal axis $\Lambda=0$, which can be read from the figures directly. The values are plotted as functions of the micro-scale Rossby number $Ro_{k}$ in Fig. 10.

Interestingly, the figure shows that $k_{c}\eta$ essentially does not depend on rotation when the flows are driven by Kolmogorov forcing, within the range of rotation rates we have considered. To two decimal places, $k_{c}\eta=0.20$ or $0.19$ for all rotation rates, which is the value obtained in Yoshida et al. (2005) for non-rotating turbulence.

This result seems to be contradictory to the observation that the decay rate for a given $k_{m}$ decreases with rotation. However, it can be explained as follows. The decay rates of the synchronisation error are reduced when rotation is introduced, which leads to increased $k_{c}$. However, as Table 1 shows, $\eta$ is decreased by rotation in this type of flows. The end result is that $k_{c}\eta$ remains roughly a constant.

For the flows driven by constant power forcing, it appears from Fig. 10 that there is a consistent trend where $k_{c}\eta$ decreases as $Ro_{k}$ decreases (i.e., as rotation rate increases). For the smallest $Ro_{k}$, $k_{c}\eta$ is reduced to below $0.15$. Therefore, for the flows driven by constant power forcing, rotation does increase the synchronizability of the flow, and this is reflected in both an increased decay rate for the synchronisation error, and a reduced $k_{c}\eta$.

Leoni et al. (2020) reported that rotating turbulence was easier to synchronize and they attributed it to the coherent vortices induced by rotation. Our results for constant power forcing are consistent with their finding, which thus might be explained qualitatively in a similar way. However, the results for Kolmogorov forcing show that the physical picture can depend on the forcing scheme.

It is instructive to cross-check the results for $k_{c}$ with the energy spectra of the flows. Note that in the majority cases the spectra are consistent with a $k^{-2}$ power law (c.f. Fig. 1). Therefore, the dimensionless threshold wavenumber $k_{c}\eta$ remains approximately unchanged from the value for isotropic turbulence when the spectrum steepens from $k^{-5/3}$ to $k^{-2}$. On the other hand, when the slope of the energy spectra is further steepened, reaching approximately that of the $k^{-3}$ power law, $k_{c}\eta$ does become smaller, as in the cases with the constant power forcing when $\Omega=1$.

To parametrise the decay rate of the synchronisation error with a physical quantity, Yoshida et al. (2005) look into the enstrophy content in the master modes. Let $H_{\omega}^{m}=\sum_{k<k_{m}}k^{2}E(k)$ be the enstrophy contained in the master modes and $H_{\omega}=\sum_{k}k^{2}E(k)$ be the enstrophy of the whole velocity field. They find that, in isotropic turbulence, the decay rate $\alpha$ is a universal function of the ratio $H_{\omega}^{m}/H_{\omega}$, and the ratio at the threshold wavenumber $k_{c}$, denoted by $H_{\omega}^{c}/H_{\omega}$, is approximately $0.35$. We plot $H_{\omega}^{c}/H_{\omega}$ as a function of $Ro_{k}$ in Fig. 11 for our simulations. For the largest $Ro_{k}$, the ratio is approximately $0.36$, which is close to the value found in Yoshida et al. (2005). The ratio consistently increases with rotation for both forcing terms, even for larger values of $Ro_{k}$ where $k_{c}\eta$ remains a constant in the flows with Kolmogorov forcing. However, the results for different cases do not collapse on a unique curve. Therefore, it seems that $H^{c}_{\omega}/H_{\omega}$ does not provide a simple way to characterise $k_{c}$ in rotating turbulence.

Another way to characterize $k_{c}$ is put forward by Leoni et al. (2020), where they observe that $k_{c}$ roughly marks the end of the inertial range. The observation is corroborated in Nikolaidis & Ioannou (2022). As the Reynolds numbers for our simulations are relatively small, this observation is not assessed here even though it is highly desirable to do so. Rather, we comment on a potential relationship between $k_{c}$ and the energy spectrum of the Lyapunov vector ${\bm{u}}^{\delta}$, which provides another perspective into the threshold wavenumbers.

Fig. 12 plots the energy spectra of ${\bm{u}}^{\delta}({\bm{x}},t)$ averaged over $t$ in the stationary stage. As the magnitude of ${\bm{u}}^{\delta}$ is irrelevant, the energy spectra have been normalised such that the total energy is unity. Also note that included in this figure are the results with $k_{m}=0$, i.e., they are the spectra of the unconditional Lyapunov vectors.

The left panel of Fig. 12 is for the flows driven by Kolmogorov forcing. First of all, the energy spectra peak at an intermediate wavenumber. That is, the perturbations with energy localised on intermediate wavenumbers are the most unstable. This observation is consistent with Ohkitani & Yamada (1989) where the Lyapunov vector for a shell model is calculated, and they find the energy spectrum of the Lyapunov vector is localised in the inertial range. Interestingly, the peaks of the spectra here are all found around $k\eta=0.2$, i.e. around the threshold wavenumber. Similar features are found in the right panel of the figure, which is the results for constant power forcing. Again in most cases the peaks are found around $k_{c}\eta$. In particular, for the two cases with $\Omega=1$, the peaks are found to shift to lower $k\eta$, consistent with Fig. 10 which shows $k_{c}\eta$ is also reduced in these two cases.

It is desirable to compare the peak wavenumbers with $k_{c}$ quantitatively. There are some challenges in extracting precise peak wavenumbers due to two factors: firstly, the spectra of ${\bm{u}}^{\delta}$ at lower wavenumbers display stronger statistical fluctuations; secondly, the gap between two data points on the spectra is $\Delta k=1$, which is fairly large and potentially introduces error into the reading of the peak wavenumbers. In order to reduce the uncertainty, we average the spectra over five realisations. We then fit a smooth curve to the spectra using cubic splines. The peak wavenumber of the fitted curve is taken to be the peak wavenumber of the spectrum.

The cubic spline fitting is conducted using scipy function UnivariateSpline, with smoothing factor $s$ chosen as $0.01\%$ of the maximum of the spectrum, which implies the 2-norm of the residue of the fitting is smaller than $s$. In other words, only a very small amount of smoothing is allowed.

The peak wavenumbers extracted in the above manner are plotted in Fig. 13 together with $k_{c}$ which has been shown in Fig. 10. The peak wavenumber obtained this way usually falls between two integer wavenumbers. These two wavenumbers are used to define the error bars in Fig. 13. The figure confirms the qualitative comments we made previously. The peak wavenumbers are slightly larger than $k_{c}$ in most cases. However they do display same trends as $k_{c}$. In particular, for flows driven by constant power forcing, the peak wavenumber clearly drops off significantly for the smallest $Ro_{k}$, despite the uncertainty in the data.

One plausible explanation of the correlation between $k_{c}$ and the peak wavenumber of the energy spectrum of ${\bm{u}}^{\delta}$ is as follows. Let the coupling wavenumber be $k_{m}$ in a synchronisation experiment. The peak wavenumber corresponds to the Fourier modes most susceptible to infinitesimal perturbations (on average). One may hypothesize that, to synchronise two flows, the perturbations to these most unstable Fourier modes should be suppressed by the coupling in the synchronisation experiments. This suggests that the coupling wavenumber $k_{m}$ should be larger than the peak wavenumber. However, even though only Fourier modes with wavenumbers up to $k_{m}$ in the two flows are coupled by design (in fact, they are exact copies of each other), the Fourier modes with wavenumbers slightly larger than $k_{m}$ are also strongly coupled, due to the fact that they are linked to the master modes through nonlinear inter-scale interactions. The coupling suppresses the growth of the synchronisation errors in these modes. Therefore, synchronisation can still be achieved even if $k_{m}$ is slightly smaller than the peak wavenumber. As a result, the threshold wavenumber $k_{c}$ could be slightly smaller than the peak wavenumber in the spectrum of ${\bm{u}}^{\delta}$.

The spectra of the Lyapunov vectors corresponding to the conditional CLEs, namely the conditional Lyapunov vectors, are given in Fig. 14 for the cases where $N=128$, and compared with the unconditional ones. For readability of the figures, only cases where $\Omega=0.1$ and $1$ with selected $k_{m}$ are included. Note that the spectra for the conditional Lyapunov vectors start from wavenumber $k_{m}+1$ as $\hat{{\bm{u}}}^{\delta}({\bm{k}},t)=0$ for $|{\bm{k}}|\leq k_{m}$. The value of $k_{m}\eta$ for each case is shown in the parentheses, which can be compared with the value of $k_{c}\eta$ to determine if synchronisation is achievable in the case. The interesting observation to note is demonstrated most clearly by both the solid and the dashed green lines with triangles in the left panel, which correspond to $k_{m}=3$ for $\Omega=0.1$ and $1$, respectively. The flows do not synchronise in these two cases while they do in all other cases depicted in the panel. The common feature of the spectra in these two cases is that the spectra peak at wavenumbers corresponding to the slaved modes. On the other hand, the spectra for the other cases all peak at $k_{m}+1$. The same trend can be observed in the right panel of the figure, and for cases with $N=192$ shown in Fig. 15. It appears that synchronisation can be achieved only when the energy spectrum of the conditional Lyapunov vector does not have a local maximum among the slave modes.

The above results for the conditional Lyapunov vectors, though are of a qualitative nature, also suggest that the threshold wavenumber $k_{c}\eta$ might be associated with the peak of the spectrum of the Lyapunov vector. Given that our simulations cover only a moderate range of Rossby numbers, with relatively low Reynolds numbers, how this observation generalises to a wider range of parameter values requires further investigation.

As previously commented, the local CLEs $\Gamma$ display significant fluctuations. In this subsection, we present statistics of $\Gamma$ for a few selected cases to highlight the qualitative trends that are shared by the other cases. The statistics in this subsection are all calculated by averaging over time as well as five independent realisations. The average is denoted by $\langle~{}\rangle$.

The variance of $\Gamma$, $Var(\Gamma)=\langle(\Gamma-\langle\Gamma\rangle)^{2}\rangle,$ is shown in Fig. 16 for the cases indicated in the figure. The variance clearly increases with rotation in all cases for a given $k_{m}$, and for a given $\Omega$, it is smaller for larger $k_{m}$. The behaviours at high rotation rates are different for the two forcing terms. The variance for constant power forcing seems to increase slower at higher rotation rates.

The PDFs of $\Gamma$, shown in Fig. 17 for the same selected cases, largely exhibit the same behaviours already shown by the mean CLEs and the variances. A common feature is that the width of the PDF increases with the rotation rate. Note that the PDFs are not normalised. The increase in the width is thus a manifestation of increased variance.

For Kolmogorov forcing (left panel), the PDFs of the unconditional $\Gamma$ (with $k_{m}=0$) do not move significantly with the rotation rate. For $k_{m}=7$, the PDF moves towards the positive values as rotation is strengthened, indicating increased mean CLE. These behaviours are consistent with Fig. 8. The behaviours of the PDFs for constant power forcing (right panel) are also consistent with Fig. 8. One notable difference with the results in the left panel is the PDFs for the unconditional Lyapunov exponents are affected more strongly by rotation in this case. For example the PDF for $\Omega=1$ is moved to the left significantly, while the same is not observed for the corresponding case in the left panel.

At higher rotation rates, the PDFs often have significant probabilities to take both positive and negative values, e.g., those for F1N192$\Omega$1K7, F1N192$\Omega$05K7 and those with constant power forcing and $k_{m}=5$. Therefore, for these cases, even if synchronisation is achieved in the long term, the synchronisation error $\Delta(t)$ may increase temporarily when the local CLE is positive. This behaviour is observed in Figs. 5 and 6 for some $k_{m}$ values near the threshold wavenumber.

Some understanding of $\Gamma$ and $\Lambda$ can be gained from Eq. (21). Our calculation shows that the contribution of the forcing term ${\mathcal{F}}$ is always negligible. In what follows we only present the results related to the production term ${\mathcal{P}}$ and the dissipation term ${\mathcal{D}}$. We use

to denote the averaged non-dimensionalised production and dissipation terms, respectively.

The values of $P$ and $D$ are shown in Figs. 18–19. For completeness, results for $N=128$ and $192$ have both been included, but we will mainly comment on those for $N=128$ as the trends are the same for $N=192$. Fig. 18 shows the results for the cases with Kolmogorov forcing. We can see that both $P$ and $D$ depend strongly on $k_{m}$, but are less sensitive to the value of $\Omega$. The production term $P$ decreases as $k_{m}$ increases, while the dissipation $D$ increases in the mean time. Thus both contribute to the decrease in $\Gamma$, hence $\Lambda$, as $k_{m}$ increases. For a given $k_{m}$, $P$ increases slightly with rotation rate $\Omega$. On the other hand, $D$ increases slightly with $\Omega$ for smaller $k_{m}$, but decreases with $\Omega$ for larger $k_{m}$. Overall, $P$ and $D$ change only slightly with $\Omega$.

For the cases with constant power forcing, Fig. 19 shows that the main impact of rotation is on the production term $P$. However, in this case $P$ decreases as rotation is increased, i.e., the trend is opposite to what is shown in Fig. 18. This trend is consistent with the previous observation that in this case synchronisation is easier when $\Omega$ is larger. Dissipation $D$ does not strongly depend on $\Omega$.

Further insights into $P$ can be obtained by looking into the alignment between ${\bm{u}}^{\delta}$ and $s_{ij}$. Let $s^{+}_{ij}\equiv\tau_{k}s_{ij}$ be the dimensionless strain rate tensor. Let ${\bm{v}}={\bm{u}}^{\delta}/\|{\bm{u}}^{\delta}\|$, and $\lambda^{s}_{\alpha}\geq\lambda^{s}_{\beta}\geq\lambda^{s}_{\gamma}$ be the eigenvalues of $s^{+}_{ij}$, with corresponding eigenvectors ${\bm{e}}_{i}$ $(i=\alpha,\beta,\gamma)$. Due to incompressibility, we have $\lambda^{s}_{\alpha}+\lambda^{s}_{\beta}+\lambda^{s}_{\gamma}=0$ with $\lambda^{s}_{\alpha}\geq 0$ and $\lambda^{s}_{\gamma}\leq 0$. In isotropic turbulence, it is well-known that $\lambda^{s}_{\beta}$ is more likely to take positive values so that the magnitude of $\lambda^{s}_{\gamma}$ tends to be the largest among the three.

Letting the angle between ${\bm{e}}_{i}$ and ${\bm{v}}$ be $\theta_{i}$, we may write

where

with $P_{\alpha}\geq 0$ and $P_{\gamma}\leq 0$. The above expressions show that $P$ is closely related to the alignment between ${\bm{u}}^{\delta}$ and the eigenvectors of $s^{+}_{ij}$, the magnitudes of the eigenvalues, and the correlations between them. As ${\bm{v}}$ is normalised, it is reasonable to expect that its magnitude is insensitive to rotation, and that rotation will mainly affect $P$ through the eigenvalues and $\cos\theta_{i}$. Since $P$ is always positive in our simulations (c.f., Figs. 18 and 19), $P_{\gamma}$ is the dominant term in $P$. As a result, we will only consider the statistics of $\cos\theta_{\gamma}$ and the eigenvalues.

The PDFs of $\cos\theta_{\gamma}$ are given in Fig. 20 for selected cases, with the left panel showing the results for Kolmogorov forcing and the right panel showing those for constant power forcing. It is evident that there is a preferable alignment between ${\bm{e}}_{\gamma}$ and ${\bm{u}}^{\delta}$ when rotation is weak, since the PDFs peak at $\cos\theta_{\gamma}=1$. Interestingly, the alignment is weaker for larger $k_{m}$, and is also weakened by rotation. These trends are consistently observed for both forcing terms. The mean values of the eigenvalues are given in Fig. 21. Here, the results for the two forcing terms display different trends. For constant power forcing, the mean eigenvalues are almost independent of the rotation rate. For Kolmogorov forcing, the magnitude of the averaged $\lambda^{s}_{\alpha}$ and $\lambda^{s}_{\gamma}$ both increase significantly with rotation.

Putting the results in Figs. 20 and 21 together, the physical picture for flows with constant power forcing appears to be simple. The production term $P$ decreases with rotation in this case, because the preferable alignment between ${\bm{e}}_{\gamma}$ and ${\bm{u}}^{\delta}$ is reduced by rotation. For the flows driven by Kolmogorov forcing, the preferable alignment is reduced by rotation, which tends to reduce $P$. However, this trend is opposed by the trend where the eigenvalues of $s^{+}_{ij}$ increase with rotation. The overall effect is that $P$ increases only slightly with rotation.

The different consequences of the two forcing terms have been made quite obvious from our analysis so far. However, the cause of the difference is not yet elucidated. The Kolmogorov forcing term is inhomogeneous whereas the constant power forcing is isotropic. However, if rotation is absent, this difference alone does not lead to significant differences in the statistics we have examined, notwithstanding the fact that the Reynolds numbers of the flows are not large. This assertion is supported by the various statistics obtained for $\Omega=0.1$, i.e., for weakest rotation. For example, Fig. 20 shows that the alignment is roughly the same for the two forces when $\Omega=0.1$. Fig. 21 shows that the mean eigenvalues are also almost the same for the two flows when $\Omega=0.1$. The same can be said about the (conditional) Lyapunov exponents as well (see, e.g., Figs. 8 and 9). Therefore, if there is no rotation, the results would be more or less independent of the forcing mechanism even if the Reynolds number is moderate, i.e. even if there is only moderate scale separation between the forced large scales and the small scales.

One may thus conclude that the drastic impacts of the forcing terms originate from the interaction between forcing and rotation. The interaction seems to alter the spectral dynamics profoundly, eluding simple phenomenological explanations (see, e.g., Dallas & Tobias (2016) and references therein). A corollary is that it is also unlikely to obtain simple mechanical explanations for the difference in the behaviours of the production term, hence those of the Lyapunov exponents, shown in previous subsections. Nevertheless, to shed some light on the physics behind the observations, we look into how different scales of the flow contribute to the production term, and how these contributions depend on the rotation rate.

To do so, we decompose the normalised strain rate tensor $s^{+}_{ij}$ into a large scale component $s^{+>}_{ij}$ and a small scale component $s^{+<}_{ij}\equiv s^{+}_{ij}-s^{+>}_{ij}$. $s^{+>}_{ij}$ is obtained by applying a low-pass filter on $s^{+}_{ij}$ (Pope, 2000). The production term calculated with $s^{+>}_{ij}$ and ${\bm{v}}$ is denoted by $P_{l}$, and that with $s^{+<}_{ij}$ and ${\bm{v}}$ is denoted by $P_{s}$, where obviously $P_{l}+P_{s}=P$. $P_{l}$ and $P_{s}$ represent the contributions from large and small scale straining, respectively, to the total production $P$.

We use the Gaussian filter (Pope, 2000) to calculate $s^{+>}_{ij}$. The filter length scale $\Delta$ is chosen to be eight times of the grid size, with the corresponding filter wavenumber $k_{\Delta}\equiv\pi/\Delta$ being $12$ for $N=192$, and $8$ for $N=128$. The values for $P$, $P_{l}$ and $P_{s}$ at different rotation rates are plotted in Fig. 22 for $k_{m}=0$. Several observations are evident. Firstly, $P_{l}$ is larger than $P_{s}$ for both forcing terms and both Reynolds numbers. That is, large scale straining makes larger contribution to the total production term, which is consistent with the fact that the spectrum of ${\bm{u}}^{\delta}$ peaks at intermediate to low wavenumbers (c.f. Fig. 12). Secondly, for constant power forcing, both $P_{l}$ and $P_{s}$ decrease as $\Omega$ increases. Both contribute roughly equally to the decrease of the total production. Thirdly, for Kolmogorov forcing, the picture again is quite different. Interestingly, $P_{l}$ barely increases (or decreases slightly) as $\Omega$ increases, whereas $P_{s}$ increases with $\Omega$ at both Reynolds numbers. Even though $P_{l}$ makes bigger contribution to $P$, the change in $P$ with $\Omega$ comes mainly from $P_{s}$.

The results in Fig. 22 are again non-trivial to interpret fully. If the impacts of forcing are confined in large scales, small scale contribution $P_{s}$ should behave in similar ways in flows driven by different forcing terms, but this is not supported by Fig. 22. If the effects of Kolmogorov forcing (coupled with rotation) on smaller scales decrease with increasing scale separation according to naive Kolmogorov phenomenology, then one should reasonably expect $P_{l}$ depends more strongly on $\Omega$ compared with $P_{s}$, which again is not supported by Fig. 22. Overall, like previous research (Dallas & Tobias, 2016), these observations suggest that large scale forcing affects the spectral dynamics of rotating turbulence in highly non-trivial ways, which is the root cause of the different synchronisability of the flows.