# Minimal model of quasi-cyclic behaviour in turbulence driven by Taylor-Green forcing

Ryo Araki, Wouter Bos, Susumu Goto

► To cite this version:

Ryo Araki, Wouter Bos, Susumu Goto. Minimal model of quasi-cyclic behaviour in turbulence driven by Taylor-Green forcing. Fluid Dynamics Research, 2023, 55 (3), pp.035507. 10.1088/1873-7005/acdff7. hal-04271757

HAL Id: hal-04271757

https://hal.science/hal-04271757v1

Submitted on 6 Nov 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Minimal model of quasi-cyclic behaviour in turbulence driven by Taylor–Green forcing

Ryo Araki

Univ Lyon, École Centrale de Lyon, CNRS, Univ Claude Bernard Lyon 1, INSA
Lyon, LMFA, UMR5509, 69130, Écully, France

Graduate School of Engineering Science, Osaka University, 1-3 Machikaneyama,
Toyonaka, Osaka 560-8531, Japan

E-mail: araki.ryo@ec-lyon.fr

Wouter J. T. Bos

Univ Lyon, École Centrale de Lyon, CNRS, Univ Claude Bernard Lyon 1, INSA
Lyon, LMFA, UMR5509, 69130, Écully, France

Susumu Goto

Graduate School of Engineering Science, Osaka University, 1-3 Machikaneyama,
Toyonaka, Osaka 560-8531, Japan

March 2023

**Abstract.** We attempt to formulate the simplest possible model mimicking turbulent dynamics, such as quasi-cyclic behaviour (QCB), using only three variables. To this end, we first conduct direct numerical simulations of three-dimensional flow driven by the steady Taylor–Green forcing to find a similarity between a stable periodic orbit (SPO) at a small Reynolds number (Re) and turbulent QCB at higher Re. A close examination of the SPO allows the heuristic formulation of a three-equation model, representing the evolution of Fourier modes in three distinct scales. The model reproduces the continuous bifurcation from SPO to turbulence with QCB when Re is varied. We also demonstrate that, by changing model parameters, the proposed model exhibits a discontinuous transition from steady to chaotic solutions without going through an SPO.

## 1. Introduction

The dynamics of turbulent flows is determined by the collective behaviour of a large number of interacting modes. The very large number of triadic interactions between these modes, even in moderately turbulent flows, prevent us from understanding the global flow features through the direct analysis of the interactions [1, 2, 3]. To gain such understanding, the complexity of the description needs to be drastically reduced.
Systematic approaches for that purpose have been applied to turbulence research, such as proper orthogonal decomposition (POD) [4, 5, 6] or Galerkin truncation [7]. For a list of the key publications in projection-based reduced-order modelling, see a recent review article [8]. In these approaches, the complex dynamics are dissected by projecting the Navier-Stokes equations on a low-dimensional basis of eigenfunctions.

A more heuristic approach to reduce complexity is the direct modelling of the dynamics by a small number of variables retaining a number of constraints (such as energy or helicity conservation). A well-known example of such an approach is the development of shell models [9, 10], as first proposed by Obukhov [11]. This approach bypasses the definition of the basis function by directly modelling the dynamics of an ensemble of modes.

In the current study, we combine observations of direct numerical simulations (DNS) and heuristic modelling. We assess the detailed dynamics of a numerical simulation and investigate the interactions between small groups of Fourier modes. Inspired by the form of the Navier-Stokes equations, we represent the full dynamics by an ODE system of three interacting variables, yielding a sort of shell model with both dyadic and triadic interactions between the groups of modes.

The turbulent flow we characterise is incompressible Navier-Stokes turbulence, driven by a large-scale steady forcing in a spatially periodic domain. At low Reynolds numbers, the considered flow becomes temporally periodic. We will show that, even for this specific periodic flow, retaining the modes governing both energy and enstrophy in the flow considered in the present investigation leads to a subset of several dozens of complex-valued Fourier modes. This results in a system of an important number of coupled ODE, which will not allow analytical treatment. Therefore, using a more heuristic approach, we analyse the periodic flow, identify the key interactions between scales, and formulate the simplest model which retains these interactions and the forcing and dissipation mechanisms. This approach allows us to formulate a model containing only three degrees of freedom, reproducing certain characteristics of the investigated fluid flow. In particular, one feature we want to reproduce with our model is quasi-cyclic behaviour (QCB).

A number of laminar and turbulent flows display QCB. An illustrative example is vortex shedding behind an obstacle. For low Reynolds number (Re), the so-called von Kármán vortex street behind a cylinder is perfectly periodic, which corresponds to a stable periodic solution (SPO) in phase space. Even when the flow becomes fully turbulent at higher Re, this periodicity is still present, though the stochastic nature of turbulence motion prevents the system from being perfectly periodic. This close-to-periodic motion, embedded in turbulent fluctuations, is what we will call QCB.

Another important example of QCB is the temporal behaviour of turbulent channel flow, where a self-sustaining process governs the dynamics [12, 13, 14]. In particular, in small channel flow domains (the so-called minimal flow unit), close to periodic behaviour is observed (Fig. 6 of Ref. [15]). The simplified descriptions of this phenomenon are specific to channel flow or the simplified case of Waleffe flow [16, 17, 18, 19, 20]. Non-
### Minimal model of quasi-cyclic behaviour in turbulence driven by Taylor-Green forcing

trivial QCB was also observed [21] in a confined cylindrical flow between two counter-rotating disks (the so-called von Kármán flow). Furthermore, QCB is observed in periodic box flow with steady forcing [22, 23, 24]. These observations suggest that such dynamics might be more general than wall-bounded flow or flow behind obstacles. Moreover, a vast amount of recent research is dedicated to identifying unstable periodic orbits (UPO) embedded in turbulent flows [25, 26, 27, 28, 29, 30, 31]. Recently, the periodic orbit-aided reduced-order model was discussed [32].

In the present study, we construct a three-mode model which exhibits turbulent QCB. Our strategy is as follows. First, we conduct DNS of turbulence in a periodic cube to find an SPO at low Re, which resembles the turbulent QCB observed at higher Re (see § 2). Secondly, in § 3, we construct a minimal model by carefully examining nonlinear interactions in the SPO. Then, in § 4.1, we demonstrate that an SPO of the model bifurcates to a chaotic solution that indeed shows QCB. In § 4.2, we also demonstrate that another route, via subcritical transition, can be possible for another set of model parameters to show the potential ability of the proposed model to be extended to other types of turbulence.

## 2. Observation of Quasi-cyclic behaviour

To illustrate the features we want to reproduce and guide the formulation of a minimal model reproducing these features, we conduct numerical simulations of both turbulent and temporally periodic flows with the same type of forcing. More precisely, we conduct DNS of three-dimensional incompressible flow governed by the Navier-Stokes equations,

$$ \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla)\mathbf{u} = -\nabla p + \nu \nabla^2 \mathbf{u} + \mathbf{f}, \quad (1) $$

with a steady forcing of the two-dimensional Taylor-Green type [22, 23, 24],

$$ \mathbf{f} = (-f_0 \sin x \cos y, f_0 \cos x \sin y, 0), \quad (2) $$

and the continuity equation, $\nabla \cdot \mathbf{u} = 0$. Here, $\mathbf{u}$, $p$, and $\mathbf{f}$ are the velocity, pressure, and forcing fields, respectively. The forcing amplitude $f_0$ is set to unity. The only control parameter is the kinematic viscosity $\nu$. We employ a pseudo-spectral method in a $(2\pi)^3$ periodic box. See Appendix A for details of the DNS. We define the Reynolds number and the characteristic timescale of large-scale flow as

$$ \mathrm{Re} \equiv \frac{\sqrt{f_0}}{|\mathbf{k}_f|^{3/2}\nu} \quad \text{and} \quad T \equiv \frac{1}{\sqrt{|\mathbf{k}_f|f_0}} = 0.840, \quad (3) $$

respectively. Here, $\mathbf{k}_f = (\pm 1, \pm 1, 0)$ is the wavevector of the forcing (2).

Figure 1 shows the temporal evolution of the energy input rate $P(t) \equiv \langle \mathbf{f} \cdot \mathbf{u} \rangle$ against the energy dissipation rate $\epsilon(t)$ given by $\nu \langle |\omega|^2 \rangle$ for various Reynolds numbers. Here, $\langle \cdot \rangle$ denotes the spatial average and $\omega \equiv \nabla \times \mathbf{u}$.
<figure><img src="image_2.png" /><figcaption>**Figure 1.** (a) Parametric plots of the instantaneous values of the energy dissipation rate $\epsilon(t)$ and the energy input rate $P(t)$ for the turbulent flow at Re = 29.7 for 50T [See Fig. B1 (b) in Appendix B]. (b) Phase-averaged values $\langle P \rangle_{\text{phase}}$ and $\langle \epsilon \rangle_{\text{phase}}$ at Re = 29.7 [same as in (a)] for 20T. (c) Parametric plots of the phase averaged values $\langle P \rangle_{\text{phase}}$ and $\langle \epsilon \rangle_{\text{phase}}$ at four different values of Re (29.7, 11.9, 8.49, and 6.61). Note that the orbit in panel (b) at Re = 29.7 is re-plotted in panel (c). (d) The SPO at Re = 5.83 is shown with a coloured line. Note that this orbit is also shown in panel (c). The gap between two consecutive dots corresponds to 5T for all panels. In panels (a), (b), and (d), the time evolves from dark to light colours. Four black cross symbols in panel (d) denote the instances shown in Fig. 2.</figcaption></figure>

Figure 1 (a) shows the turbulent time series, where the time-averaged Taylor scale-based Reynolds number $\langle \text{Re}_\lambda \rangle_t$ is about 90. Note that $\langle \cdot \rangle_t$ denotes the time average. We also show the snapshot of this flow in Fig. B1 (a) in Appendix B. The time series exhibits QCB in a counter-clockwise direction behind the chaotic fluctuations. This time delay between the large- and small-scale representatives (i.e. $P(t)$ and $\epsilon(t)$) reflects the causal nature of the energy cascade.

We apply a phase average to the complex time series of $P(t)$ and $\epsilon(t)$ conditioned on the local maxima of $P(t)$ in order to extract smooth, time-delayed oscillations shown in Fig. 1 (b). See Appendix B for the detailed procedure. We denote the phase-averaged quantities by $\langle \cdot \rangle_{\text{phase}}$. These results suggest that the QCB of turbulent flow driven by
Minimal model of quasi-cyclic behaviour in turbulence driven by Taylor-Green forcing⁵

the steady body force (2) is robust. Such QCB is also shown in Fig. 12 of Ref. [23] for two different forcing types at even higher Re. The physical origin of QCB is rooted in the energy cascading process from larger to smaller scales. Since the coherent structures at these scales are composed of a large number of Fourier modes, to describe the QCB in terms of Fourier modes, we need to understand the underlying nonlinear interactions among them. However, identifying the direct cause of QCB from tens of thousands of excited Fourier modes seems illusory. Thus, we decrease Re to reduce the complexity of the flow.

In Fig. 1 (c), we show the phase-averaged plots of the parametric time series of $P(t)$ and $\epsilon(t)$ for four different values of Re. The change in the shape of the parametric plots is gradual, suggesting that the quasi-cyclic orbit in the turbulent flow is continuously connected to an SPO at $Re \approx 5.83$, which is also shown in Fig. 1 (d) for comparison. As will be shown in Fig. 2 below, the SPO at $Re = 5.83$ is not the laminar solution of the system which corresponds to a purely two-dimensional structure resulting from a balance between viscous stress and the forcing (2). We emphasise that this SPO plays a key role in constructing our model.

We find that the amplitude and the period of the periodic and quasi-cyclic flows monotonically increase when we decrease Re from 29.7 to 5.83. This does not prove that the dynamics are identical, but the turbulent QCB and periodic flow seem to share the same driving mechanism. Note that in high-Reynolds-number turbulence beyond $Re = 30$, the amplitude and period seem to saturate to values of the same order as in Fig. 1 (b) (See Fig. 12 of Ref. [23]).

In Fig. 2, we visualise the periodic flow (i.e. SPO) discussed in Fig. 1 (d), similar to the three-dimensional periodic solution reported in [24, Fig. 5]. We distinguish four large-scale columnar vortices associated with the Taylor-Green force (2) and counter-rotating pairs of smaller vortices perpendicular to them. Note that we do not perform low-pass filtering [See Fig. B1 (a) of Appendix B] since there is no significant scale separation in the periodic flow. Nevertheless, we can observe a one-step energy cascading process from the four large-scale columnar vortices to smaller-scale lateral vortices. More concretely, we observe only large-scale vortices at $t = T_0$ [Fig. 2 (a)], then the energy cascade starts to create smaller-scale vortices [Fig. 2 (b), $t = T_0+9.42T$], while the large-scale vortices get weaker [Fig. 2 (c), $t = T_0+17.0T$]. Afterwards, the energy dissipation dominates to weaken smaller-scale vortices, and then the entire system becomes calm [Fig. 2 (d), $t = T_0+19.3T$]. When small-scale vortices disappear, energy input by the external force exceeds dissipation to reestablish the large-scale vortices, and the system returns to the initial state [Fig. 2 (a)]. We emphasise that this periodic behaviour is similar to turbulent QPB observed at higher Re (Fig. B1 of Appendix B and Figs. 12-17 of Ref. [23]). This similarity manifests itself in the continuous change between the SPO and turbulence seen in Fig. 1.

In the next section, we analyse the SPO to unveil the essential physics behind QCB. Even though we have not rigorously shown the connection between the SPO and turbulence, we hope to obtain new insights into QCB in Navier-Stokes flow by dissecting
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Visualisation of vortical structures of the SPO (at Re = 5.83) at four instances with isosurface of $|\omega| = 5$. See Fig. 1 (d) for the corresponding instances.</figcaption></figure>

the SPO.

### 3. Three-equation model

#### 3.1. Construction of the model

Our objective is to construct the simplest possible model capable of reproducing QCB while retaining a close connection with the structure of the Navier–Stokes equations (1). For this purpose, we recall that in a Fourier representation of (1), the individual modes $q_i$ for the $i$th wavevector $\mathbf{k}_i$ are governed by [33, 34],

$$ \left( \frac{\partial}{\partial t} + \nu |\mathbf{k}_i|^2 \right) q_i = \sum_{j,m} A_{ijm} q_j q_m + f_i, \quad (4) $$

where $f_i$ is the forcing applied to the $i$th mode, and $A_{ijm}$ are the coupling constants resulting from the advection and pressure terms of (1). The nonlinear term associated with triad interactions rapidly yields an overwhelming complexity when the number of retained modes increases. Even in our SPO, a large number of modes are dynamically active. In order to develop an analytically tractable model, we use a coarse-graining
approach where we group subsets of Fourier modes and represent each group by a single variable, leading to a sort of shell-model [11, 10].

The shells or groups used in our model are not regrouping modes as a function of scale using a rigorous criterion but as a function of the type of nonlinear interactions and energetic content. Indeed, we investigate the Fourier decomposition of the SPO to find that only Fourier modes with wavevectors ($k_x, k_y, k_z$) of

$$ (\pm1, \pm1, 0), (\pm1, 0, 0), (0, \pm1, 0), (0, 0, \pm2), (\pm1, \pm1, \pm2), (\pm2, 0, \pm2), \text{ and } (0, \pm2, \pm2) \tag{5} $$

are responsible for 98% of its energy. See Appendix C for details of these energetic modes. Figure 3 (a) illustrates that the time evolution of the kinetic energy is closely reproduced, retaining only these modes.

A close inspection of the seven modes shows that all the nonlinear interactions involve the forced mode and two of the six other modes (See Fig. C1 in Appendix C). In the following, $X \in \mathbb{R}$ denotes the characteristic velocity of the forced modes $k = (\pm1, \pm1, 0)$ and $Y \in \mathbb{R}$ corresponds to that of the remaining modes in (5). At this point, we suppose that there are only these two classes of modes and that we represent each class by a single, real variable. Furthermore, we assume (4) to govern the interaction of these two variables, $X$ and $Y$, yielding,

$$ \begin{aligned} \frac{dX}{dt} &= -AY^2 - \nu K_X^2 X + F, \\ \frac{dY}{dt} &= +AXY - \nu K_Y^2 Y, \end{aligned} \tag{6} $$

with a coefficient $A > 0$, typical wavenumbers $K_\alpha > 0$ with $\alpha \in \{X, Y\}$, and a steady force $F > 0$. The first term on the RHS of each equation represents the nonlinear coupling between $X$ and $Y$. This interaction conserves the global energy, $(X^2 + Y^2)/2$. Note that since we model the triadic nonlinear term of (4) by regrouping the modes into two families (See Fig. C1 in Appendix C), the resulting interactions which appear in the model (6) are dyadic. For notation, we employ both $X$ and $Y$ as the principal variables of our model and as subscripts to denote quantities associated with these variables.

An extensive parameter scan of the two-equation model shows that the model always converges to a steady solution, and we do not observe an SPO or QCB. In fact, linear stability analysis of the fixed points of (6) shows that there are only stable steady solutions (See Appendix D). Thus, retaining only this simple interaction between the forced and most energetic modes seems insufficient to reproduce QCB via supercritical bifurcations. Results of the parameter scan further suggest that the subcritical route to QCB is not present either.

The additional ingredient for QCB turns out to be a small-scale representative and its associated triad interaction terms. Figure 3 (a) shows the time series of energy $E(t)$ and energy dissipation rate $\epsilon(t)$ in the SPO along with partial energy $E_{X+Y} \equiv \langle |\mathbf{u}_X|^2 \rangle / 2 + \langle |\mathbf{u}_Y|^2 \rangle / 2$ and partial energy dissipation rate $\epsilon_{X+Y} \equiv \nu (\langle |\omega_X|^2 \rangle + \langle |\omega_Y|^2 \rangle)$ contained by the forced and primary modes. While the energy is almost entirely
<figure><img src="image_4.png" /><figcaption>**Figure 3.** (a) Time series of energy $E(t)$ and energy dissipation rate $\epsilon(t)$ computed from all modes (solid lines) and those of the forced plus the primary energetic modes, denoted by $\langle \cdot \rangle_{X+Y}$ (dashed lines). (b) Schematic of three different scales: “forced”, “primary”, and “secondary”. We visualise $|\omega|$ distributions of typical Fourier modes in each scale. The forced scale corresponds to $k_f = (\pm 1, \pm 1, 0)$. In the primary scale, we visualise $k = (0, 0, \pm 2)$ and $(0, \pm 2, \pm 2)$ modes. For the secondary scale, we visualise $k = (\pm 3, \pm 1, 0)$ and $(\pm 2, \pm 2, \pm 2)$ modes for example. Note that the contributions of modes with all the possible sign combinations $(\pm k_x, \pm k_y, \pm k_z)$ are gathered in the visualisations. Triangles denote triad interactions between different scales. For details of interactions between the forced and primary scales, see Fig. C1 in Appendix C.</figcaption></figure>

contained in $E_{X+Y}$, there is a visible difference between the full and partial energy dissipation rates. This reveals that the rest of the Fourier modes contribute significantly to the dynamics of the energy dissipation, representing the small scales. We denote the ensemble of these residual modes by $Z$. The essential nonlinear interactions of $Z$ form triads with one mode of the $Y$-ensemble and another mode from either the $Z$-ensemble or the forced mode $X$. These observations lead to a refined three-equation model,

$$ 
\begin{aligned}
\frac{dX}{dt} &= -A_1 Y^2 + A_3 Y Z - \nu K_X^2 X + F, \\
\frac{dY}{dt} &= +A_1 XY - A_2 Z^2 + A_4 X Z - \nu K_Y^2 Y, \\
\frac{dZ}{dt} &= + A_2 Y Z - (A_3 + A_4)XY - \nu K_Z^2 Z,
\end{aligned}
\qquad (7) $$

which is represented by a schematic in Fig. 3 (b). Here, $A_1, A_2 > 0$ and $A_3, A_4 \in \mathbb{R}$ are triad coefficients which retain the discrete Navier–Stokes structure (4). We choose the
### Minimal model of quasi-cyclic behaviour in turbulence driven by Taylor-Green forcing⁹

signs and the values of the triad coefficients such that the detailed balance holds in the energy transfer between the three scales. The signs of $A_1$ and $A_2$ are defined so that energy cascades towards small scales: from X to Y and Y to Z. This two-step energy cascade (for $A_3 = A_4 = 0$) is similar to the Obukhov two-stage cascade model [11]. The triads with coefficients $A_3$ and $A_4$ represent the “non-local” interactions involving all three scales. Note that this system is different from, but is of the same level of complexity, as the well-known Lorenz [35] or Rössler models [36, 37]. An important difference is that each variable denotes a Fourier mode in the Lorenz model, while in our model, it represents a group of modes.

## 3.2. Determination of the parameters

The model (7) is a simplified representation of the SPO, where all Fourier modes are sorted into three scales; the forced mode X, the energetic modes Y directly draining energy from X through the $A_1$ interaction, and the small scale modes Z which couple through the local direct cascade interaction $A_2$ with Y. There are also scale non-local interactions represented by $A_3$ and $A_4$. Even though such a representation of the flow discards details of the actual flow obtained by the DNS, we will fit the model parameters to the DNS data to assess how the model can reproduce actual flow properties.

We can fit six out of eight model constants in (7) by comparing them to the DNS of the periodic flow: $A_i$ with $i = 1, 2, 3, 4$, $K_\alpha^2$, where $\alpha \in \{X, Y, Z\}$, and $F$. To do so, we use the energy equations associated with (7),

$$ 
\begin{aligned}
\frac{dE_X}{dt} &= T_X - \epsilon_X + P, \\
\frac{dE_Y}{dt} &= T_Y - \epsilon_Y, \\
\frac{dE_Z}{dt} &= T_Z - \epsilon_Z.
\end{aligned}
\qquad (8) $$

Here, $E_\alpha \equiv \alpha^2/2$ is the energy,

$$ 
\begin{aligned}
T_X &\equiv -A_1XY^2 + A_3XYZ, \\
T_Y &\equiv +A_1XY^2 - A_2YZ^2 + A_4XYZ, \\
T_Z &\equiv +A_2YZ^2 - (A_3 + A_4)XYZ.
\end{aligned}
\qquad (9) $$

are the energy transfer terms, $\epsilon_\alpha \equiv 2\nu K_\alpha^2 E_\alpha$ is the energy dissipation rate, and $P \equiv FX$ is the energy input rate. The model parameters are determined by their corresponding quantities of the SPO obtained by DNS. The resulting values are

$$ A_1 = 0.4, A_2 = 4, F = 0.7, K_X^2 = 2, K_Y^2 = 5, \text{ and } K_Z^2 = 15. \qquad (10) $$

We describe the parameter determination procedure in Appendix E. The energy flux coefficients $A_1$ and $A_2$ are determined by the energy transfer terms $T_\alpha$ in (9) while ignoring the nonlocal coefficients $A_3$ and $A_4$ (E.4-E.5). The forcing coefficient $F$
### Minimal model of quasi-cyclic behaviour in turbulence driven by Taylor-Green forcing¹⁰

is evaluated by $P$ and the $X$-scale energy $E_X$ (E.1). The squared characteristic wavenumber $K_\alpha^2$ is set by $E_\alpha$ and $\epsilon_\alpha$ in each scale (E.2). Note that $K_X = \sqrt{2}$ of the model parameter can be related to $|\mathbf{k}_f| = \sqrt{2}$ of the forcing (2) of the DNS. We remark here that our parameter choice (10) supports the energy cascade picture with $T_X(t) < 0$: the forced scale $X$ transfers its energy to smaller scales ($Y, Z$) on average. And $T_Y(t), T_Z(t) > 0$ means that the smaller scales receive energy from the larger scales. The undetermined parameters of the model are the scale non-local interaction coefficients $A_3$ and $A_4$, which can be freely chosen. The only control parameter is $\mathrm{Re} \equiv 1/\nu$. We numerically integrate the model with a fourth-order Runge-Kutta scheme and $\Delta t = 0.01$ starting from random initial conditions. See Ref. [38] for the solver information. Our numerical simulations seem to indicate that no periodic solutions exist without the complete non-local interactions: $A_3 = 0, A_4 = 0$, or $A_3 + A_4 = 0$. Conversely, periodic behaviour is observed for a wide range of values when $A_3 \neq 0, A_4 \neq 0$, and $A_3 + A_4 \neq 0$. This observation emphasises the importance of non-local triad interactions for periodic behaviour.

### 3.3. Comparison between the model and the DNS result

Figure 4 compares the SPO obtained by the model and the DNS. Figure 4 (a) shows the time series of the model with the parameters (10) and $(A_3, A_4) = (0.5, -0.95)$. Since the definitions of $\mathrm{Re}$ are different in the model and DNS, we have chosen a Reynolds number in the model, which allows qualitatively reproducing the DNS results. We compute two quantities. One is $E_{X-X_0} \equiv (X - X_0)^2/2$, which is the fluctuating energy of the forced mode around the laminar base flow $X_0 \equiv F \mathrm{Re}/K_X^2$. The other quantity $E_{Y+Z} \equiv Y^2/2 + Z^2/2$ is the energy of the rest of the modes. We compare them to the corresponding quantities in the DNS of the SPO [Fig. 4 (b)], where the base flow is $\mathbf{u}_0 \equiv \mathbf{f}/2\nu |\mathbf{k}_f|^2$, the forced-mode fluctuating energy is $E_{X-X_0} \equiv \langle |\mathbf{u}_X - \mathbf{u}_0|^2 \rangle /2$, and $E_{Y+Z}$ is defined by the energy possessed by the non-forced modes. We can observe similar periodic behaviour of $E_{X-X_0}$ and $E_{Y+Z}$ in the model (7) and in the SPO driven by the steady forcing (2). In particular, there are predator-prey-like exponential growth and decay in both systems. Although fast oscillations are observed in the model but not in the DNS, a close analysis (See Appendix C) of the DNS of the SPO reveals the presence of rapid oscillations in specific Fourier modes. These oscillations are compensated by modes that display the same energy oscillations with an opposite phase and do not appear in Fig. 4 (b). We stress that this SPO is independent of the exact amplitude of the initial conditions because the present model is a dissipative system.

## 4. Dynamics of the model

### 4.1. Bifurcation from SPO to chaos with QCB

We observe a chaotic state of the model by varying $\mathrm{Re}$ from 14.05 to 14.1 while keeping the model parameters as in Fig. 4 (a). Figure 5 (a) shows the orbits in phase space
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Time series of fluctuating energy $E_{X-X_0}(t)$ of the forced scale and residual energy $E_{Y+Z}(t)$ of periodic solutions of (a) model (7) at Re = 14.05 and (b) the Navier-Stokes equations (1) at Re = 5.83. Parameters of the model are (10) and $(A_3, A_4) = (0.5, -0.95)$. Note that time in panel (b) is normalised by T.</figcaption></figure>

for both the periodic (at Re = 14.05) and chaotic (at Re = 14.1) cases. The chaotic solution remains close to the SPO as it shows chaotic QCB and is permanent as in the turbulence investigated in § 2. Thus, the same model reproduces SPO and chaotic QCB. Incidentally, the SPO resembles a Shilnikov homoclinic orbit [39]. We also plot a simpler periodic orbit at Re = 12.5 in this figure. It is almost two-dimensional as opposed to the complex three-dimensional periodic and chaotic orbits, suggesting a possible connection with a two-dimensional periodic orbit in the same forcing configuration [24, Fig. 4]. However, we do not focus on this orbit as it is not directly connected to a chaotic one.

To further assess the behaviour of the system, we draw the bifurcation diagram in Fig. 5 (b) with the same parameter set as in Fig. 4 (a) and Fig. 5 (a). We observe a supercritical transition from periodic to chaotic solutions at a critical Reynolds number $Re_{cr} \in [14.060, 14.061]$ [inset of Fig. 5 (b)], and, as observed in Fig. 5 (a), the chaotic orbit remains close to the SPO. We note that the solution becomes periodic again when we further increase Re beyond the range of Fig. 5 (b), probably because the model contains only a small number of degrees of freedom. The inset of Fig. 5 (b) shows that
<figure><img src="image_6.png" /><figcaption>**Figure 5.** (a) Simple periodic (Re = 12.50), complex periodic (Re = 14.05) [Fig. 4 (a)], and chaotic (Re = 14.1) orbits of the model (7). The parameters are the same as in Fig. 4 (a). The chaotic orbit is tracked over 100 periods. The arrow indicates the direction of the orbit. (b) Bifurcation diagram of the model with changing Re for the same parameters as in Fig. 4 (a). We plot the local extrema of Y. We have determined the periodicity by Poincaré analysis. The black vertical dotted line corresponds to Re = 12.50 used in panel (a). Inset: close-up in the range shown by the red rectangle in the main plot. The red vertical dashed line corresponds to Re = 14.05 used for panel (a) and Fig. 4 (a). For both panels, black, purple, orange, and yellow data denote steady, simple periodic, complex periodic, and chaotic solutions, respectively.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 6.** (a) Time series of (X, Y, Z) of the model (7) with parameters (10) and (A₃, A₄) = (0.4, -0.5) at Re = 32. A random initial condition is used. (b) The bifurcation diagram for the same parameter set. The red vertical dashed line denotes Re = 32, which is used for Fig. 6 (a). Inset: close-up of the diagram in the range shown by the red rectangle in the main plot.</figcaption></figure>

there is a hysteresis in the range Re ∈ [13.82, 14.03], below Recr, which corresponds to a subcritical bifurcation from a periodic solution to another periodic solution shown in Fig. 4 (a). The appearance of the multiplicity of local extrema corresponds to the spiralling behaviour of the orbit in phase space. Thus, although the bifurcations from the trivial steady solution to the SPO are rather complicated, that from the SPO to chaos with QCB is simple. Although there is no clear scenario for the route to turbulence with QCB, the present model results may give us a hint to describe the route in real turbulence.

## 4.2. Subcritical bifurcation to chaos

Since it is well known that, in some cases, turbulence appears via a subcritical transition, here we demonstrate that our model also expresses such a route to chaos. We stress that we cannot use the strategy above to determine the model parameters since there
<figure><img src="image_8.png" /><figcaption>**Figure 7.** (a) Survival probability $P_{\text{Re}}(t)$ of the transient chaos of the model (7) evaluated from 10,000 samples for each Re. The parameter set is the same as in Fig. 6. Dashed line denotes exponential fitting by (11) using $0.01 \le P_{\text{Re}}(t) \le 0.9$. (b) The escape rate $1/\tau$ as a function of Re. The dashed line denotes the exponential fitting by (12).</figcaption></figure>

is no SPO in such a system. Instead, by varying the undetermined parameters of the model, we observe transient chaos at $(A_3, A_4) = (0.4, -0.5)$ as shown in Fig. 6 (a). The corresponding bifurcation diagram in Fig. 6 (b) shows a subcritical bifurcation between steady and chaotic solutions around $\text{Re} \approx 32.3$. There are bi-stable states for $33 \le \text{Re}(\le 35)$ of steady and chaotic solutions. The inset of Fig. 6 (b) shows that there are multiple windows of periodic solutions in the chaotic regime, probably due to the limited number of degrees of freedom of the model (7).

The transient behaviour in Fig. 6 (a) reminds us of the sudden relaminarisation observed in a linearly forced turbulence [40], turbulent Kolmogorov flow [41], pipe flow [42], and even in the Lorenz system [43, 44]. We evaluate the survival probability $P_{\text{Re}}(t)$, representing how likely the solution remains in a chaotic regime at a given time $t$, to investigate this phenomenon. To evaluate $P_{\text{Re}}(t)$, we identify the relaminarisation time $t_r$ by the first time when the local maxima of oscillating energy $E_{Y-Y_0} \equiv$
Minimal model of quasi-cyclic behaviour in turbulence driven by Taylor-Green forcing15

$(Y - Y_0)^2/2$ becomes smaller than a threshold $\delta = 1 \times 10^{-3}$. Here, $Y_0$ is the stable and steady solution. Then, the probability $P_{\text{Re}}(t)$ for given $t$ can be evaluated by the ratio of a number of samples with $t_r < t$ against the number of the whole sample. We plot $P_{\text{Re}}(t)$ in Fig. 7 (a) to find that an exponential scaling,

$$P_{\text{Re}}(t) \propto \exp\left[-\frac{t}{\tau(\text{Re})}\right], \qquad (11)$$

fits the data. The characteristic time scale $\tau$ in Fig. 7 (b) also displays an exponential scaling,

$$\tau(\text{Re}) \propto \exp[a \text{Re}], \qquad (12)$$

against Re. Although the scaling (11) of $P_{\text{Re}}(t)$ is consistent with the observations in the previous studies [40], the exponential scaling (12) of $\tau(\text{Re})$ differs from a super-exponential behaviour observed in Ref. [40]. This qualitative difference may also be caused by the minimal number of degrees of freedom in the model.

Note that the Taylor-Green forcing (2) in the DNS does not permit such a transition since the laminar base flow $u_0 \equiv f/2\nu|k_f|^2$ is linearly unstable. However, the steady Kolmogorov forcing with a linearly stable laminar base flow exhibits sudden relaminarisations [41]. Thus, we can speculate that the model can reflect different forcing set-ups applied to the Navier-Stokes equations by varying the parameters ($A_3, A_4$).

## 5. Conclusion

The present investigation attempts to construct a minimal model of turbulence with quasi-cyclic behaviour (QCB) in a steady-force driven flow while keeping the structure of the Navier-Stokes equations. First, through the DNS of Navier-Stokes turbulence, we show that QCB in high-Re turbulence is continuously connected to an SPO at small Re by extracting the intrinsic periodicity of QCB via a phase averaging technique (§ 2). Next, we conduct a mode-by-mode analysis of the SPO to identify the flow's forced, primary energetic, and secondary scales. We propose the three-equation model (7) describing the evolution of such three distinct scales (§ 3.1). By adjusting the model parameters, we observe that the model reproduces an SPO similar to that of the DNS (§ 3.2). We emphasise that scale non-local nonlinear interactions (interactions involving three separate scales) are mandatory for reproducing these dynamics. Then, we conduct a bifurcation analysis to show that the model also exhibits chaotic QCB via a supercritical bifurcation, which is continuously connected to the SPO (§ 4.1). Thus, we conclude that the proposed model reproduces turbulent QCB and its relation to an SPO using a minimum number of degrees of freedom.

Further analysis of the model by varying the undetermined parameters yields transient chaos with sudden relaminarisation, which is also observed in turbulent flow with different forcing set-ups (§ 4.2). Thus, we speculate that the present model can be a minimal model for certain features of turbulence.
*Minimal model of quasi-cyclic behaviour in turbulence driven by Taylor-Green forcing* 16

An outstanding open question is how QCB survives in spatially extended flows. How will the global dynamics change when the forcing is applied to scales smaller than the domain size? In other words, how will the modes larger than the forced scale alter QCB turbulence, and how can we model it? Investigating the relation between space and scale locality and temporal dynamics of turbulence is left for further research.

## Acknowledgments

All DNS were carried out using the facilities of the PMCS2I (École Centrale de Lyon). R.A. is supported by the Takenaka Scholarship Foundation, and S.G. is funded by the JSPS Grant-in-Aid for Scientific Research 20H02068. The authors thank Dr. Y. Duguet for his valuable comments and suggestions on the preprint. The authors also thank the anonymous referees for many valuable suggestions to improve the manuscript. R.A. deeply thanks Dr. T. Tanogami for his insightful comments and discussions and Dr. T. Yasuda for his helpful comments on an early draft. For the purpose of Open Access, a CC-BY public copyright licence has been applied by the authors to the present document and will be applied to all subsequent versions up to the Author Accepted Manuscript arising from this submission.

## References

[1] Kraichnan R H 1971 *Journal of Fluid Mechanics* **47** 525-535

[2] Domaradzki J A and Rogallo R S 1990 *Physics of Fluids A: Fluid Dynamics* **2** 413-426

[3] Waleffe F 1992 *Physics of Fluids A: Fluid Dynamics* **4** 350-363

[4] Bakewell Jr H P and Lumley J L 1967 *The Physics of Fluids* **10** 1880-1889

[5] Berkooz G, Holmes P and Lumley J L 1993 *Annual Review of Fluid Mechanics* **25** 539-575

[6] Holmes P, Lumley J L and Berkooz G 1996 *Turbulence, Coherent Structures, Dynamical Systems and Symmetry* (Cambridge University Press)

[7] Rempfer D 2000 *Theoretical and Computational Fluid Dynamics* **14** 75-88

[8] Ahmed S E, Pawar S, San O, Rasheed A, Iliescu T and Noack B R 2021 *Physics of Fluids* **33** 091301

[9] Biferale L 2003 *Annual Review of Fluid Mechanics* **35** 441-468

[10] Ditlevsen P D 2010 *Turbulence and Shell Models* (Cambridge University Press)

[11] Obukhov A M 1971 *Academy of Sciences, USSR, Izvestiya, Atmospheric and Oceanic Physics* **7** 471-475

[12] Waleffe F 1995 *Studies in Applied Mathematics* **95** 319-343

[13] Hamilton J M, Kim J and Waleffe F 1995 *Journal of Fluid Mechanics* **287** 317-348

[14] Panton R L 2001 *Progress in Aerospace Sciences* **37** 341-383

[15] Jiménez J and Moin P 1991 *Journal of Fluid Mechanics* **225** 213-240

[16] Waleffe F 1997 *Physics of Fluids* **9** 883-900

[17] Thomas V L, Lieu B K, Jovanović M R, Farrell B F, Ioannou P J and Gayme D F 2014 *Physics of Fluids* **26** 105112

[18] Thomas V L, Farrell B F, Ioannou P J and Gayme D F 2015 *Physics of Fluids* **27** 105104

[19] Alizard F and Biau D 2019 *Journal of Fluid Mechanics* **864** 221-243

[20] Cavalieri A V G 2021 *Physical Review Fluids* **6** 034610

[21] Araki R and Goto S 2021 *Physical Review Fluids* **6** 084603

[22] Yasuda T, Goto S and Kawahara G 2014 *Fluid Dynamics Research* **46** 061413
[23] Goto S, Saito Y and Kawahara G 2017 *Physical Review Fluids* **2** 064603

[24] van Veen L, Kawahara G and Yasuda T 2018 *The European Physical Journal Special Topics* **227** 463–480

[25] Auerbach D, Cvitanović P, Eckmann J P, Gunaratne G and Procaccia I 1987 *Physical Review Letters* **58** 2387

[26] Kawahara G and Kida S 2001 *Journal of Fluid Mechanics* **449** 291–300

[27] Toh S and Itano T 2003 *Journal of Fluid Mechanics* **481** 67–76

[28] Kawasaki M and Sasa S 2005 *Physical Review E* **72** 037202

[29] Kawahara G, Uhlmann M and van Veen L 2012 *Annual Review of Fluid Mechanics* **44** 203–225

[30] Lucas D and Kerswell R 2017 *Journal of Fluid Mechanics* **817**

[31] van Veen L, Vela-Martín A and Kawahara G 2019 *Physical Review Letters* **123** 134502

[32] Yalniz G, Hof B and Budanur N B 2021 *Physical Review Letters* **126** 244502

[33] Kraichnan R H 1958 *Physical Review* **109** 1407

[34] Kraichnan R H 1988 *Journal of Statistical Physics* **51** 949–963

[35] Lorenz E N 1963 *Journal of Atmospheric Sciences* **20** 130–141

[36] Rössler O E 1976 *Physics Letters A* **57** 397–398

[37] Rössler O E 1979 *Physics Letters A* **71** 155–157

[38] Rackauckas C and Nie Q 2017 *Journal of Open Research Software* **5** 15

[39] Shilnikov L P 1965 A case of the existence of a denumerable set of periodic motions *Doklady Akademii Nauk* **160** 558–561

[40] Linkmann M F and Morozov A 2015 *Physical Review Letters* **115** 134502

[41] van Veen L and Goto S 2016 *Fluid Dynamics Research* **48** 061425

[42] Hof B, Westerweel J, Schneider T M and Eckhardt B 2006 *Nature* **443** 59–62

[43] Yorke J A and Yorke E D 1979 *Journal of Statistical Physics* **21** 263–277

[44] Maslennikov O V and Nekorkin V I 2013 *Chaos* **23** 023129

[45] Delache A, Cambon C and Godeferd F 2014 *Physics of Fluids* **26** 025104

[46] Rogallo R S 1981 *Numerical Experiments in Homogeneous Turbulence* (National Aeronautics and Space Administration)
# Appendix

## Appendix A. Direct Numerical Simulations

This appendix describes the detail of the DNS condition. We use an in-house parallelised code [45] to conduct DNS. It employs a pseudo-spectral method with the 2/3 dealiasing rule for spatial discretisation and the Adams-Bashforth scheme in the time domain. The initial condition is generated in Fourier space by Rogallo's method [46].

We perform DNS in a $(2\pi)^3$ triply periodic box. We focus on two distinct flows: three-dimensional periodic and turbulent. The SPO is obtained by the DNS with $64^3$ Fourier modes by adjusting the viscosity to $\nu = 0.102$. This corresponds to the value of the Reynolds number (3) of $Re = 5.83$. We use $128^3$ Fourier modes to simulate turbulent flow at $\nu = 0.02$ (Re = 29.7). For the phase averaging procedure (See Appendix B), we use the time series in the interval $1.28 \times 10^3 \le t \le 1.04 \times 10^4$ to guarantee statistical convergence. Note that we discard the transient part from the analysis. This interval is approximately $1.08 \times 10^4 T$ with $T$ defined in (3).

## Appendix B. Detailed procedure of the phase averaging

This appendix explains the detailed procedure of the phase averaging shown in Fig. 1 (b). We first describe the flow observed at $Re = 29.7$ ($Re_\lambda$ is about 90) in Fig. B1. Figure B1 (a) shows isosurfaces of $|\omega|$ capturing small-scale structures, whereas the forcing-induced columnar vortices emerge by visualising the isosurfaces of $|\omega^{<}|$. Here, $\omega^{<} \equiv \nabla \times u^{<}$, which is obtained by applying a low-pass filter to the velocity field, defined as $u^{<}(x) \equiv \int dr G(r/r_0)u(x+r)$ with $G$ being the Gaussian function. We set $r_0 = 2$. Figure B1 (b) shows the temporal evolution of the energy input rate $P(t)$ and the energy dissipation rate $\epsilon(t)$. Both time signals exhibit significant fluctuations with a clear time-delayed correlation. See Fig. 1 (a) for the 2D projection of the same time series.

To conduct the phase-averaging, first, we pick up the local maxima of $P(t)$ [Fig. B2 (a)] with the following two criteria: (i) It must be larger than $\langle P \rangle_t + \sigma(P)$ where $\langle \cdot \rangle_t$ and $\sigma(\cdot)$ denote the time average and the standard deviation, respectively. The horizontal pink line indicates this value in Fig. B2 (a). (ii) The temporal gap between two consecutive local maxima must be larger than $\tau_{max}/2$ where $\tau_{max}$ is the time for the second peak of the autocorrelation function of $P$ (note that the first peak is at $\tau = 0$). We denote the identified local maximum of $P(t)$ and the corresponding time by $P_0$ and $t_0$, respectively. Second, the segments of the time series of $P(t)$ around the local maximum $P_0$ are overlapped, as shown in Fig. B2 (b). We normalise the segments by $P_0$ to avoid overestimation due to huge intermittent peaks. Third, we compute the average over the overlapped and normalised time series to obtain the phase averaged time series $\langle P \rangle_{\text{phase}}$ shown in Fig. B2 (c).

We apply a similar procedure to $\epsilon(t)$. However, the time is shifted for $t_0$, and $\epsilon(t)$
<figure><img src="image_9.png" /><figcaption>**Figure B1.** (a) A snapshot of the turbulent flow at $Re = 29.7$ ($Re_λ ≈ 90$). Isosurfaces of $|ω| = 20$ (blue) and low-pass filtered $|ω| < = 4$ (red) are visualised. (b) Time series of energy input rate $P(t)$ and energy dissipation rate $ε(t)$ of the same flow. The red rectangle denotes the time interval examined in Fig. 1 (a).</figcaption></figure>

is normalised by $P_0$ so that we can evaluate the time delay and the relative amplitude difference between the two quantities. The pink vertical dashed line shows the time delay in Fig. B2 (c), which is 2.80T. Figure 1 (b) is a parametric plot of Fig. B2 (c).

## Appendix C. Primary energetic modes of the periodic flow

This appendix shows the detailed behaviour of the seven most energetic modes. Figure C1 shows $|ω|$ distributions of these seven modes. Triangles indicate combinations of different modes where energy transfer via triad interactions is possible. Figure C1 also compares isosurfaces of $|ω|$ of the sum of these seven primary energetic modes and that of all modes. We find similar principal structures: the large columnar vortices and the small counter-rotating pairs of vortices. Here, we denote the velocity field consisting of the forced mode by $u_X$ and the other six primary energetic modes by $u_Y$. The corresponding vorticity fields are denoted by $ω_X$ and $ω_Y$, respectively.
<figure><img src="image_10.png" /><figcaption>**Figure B2.** (a) Time series of $P(t)$ with its local maxima $P_0$ denoted by dots. The pink horizontal line corresponds to the threshold of the magnitude $\langle P \rangle_t + \sigma(P)$. (b) Overlapped segments of the time series of $P(t)$ around $t_0$ and normalised by $P_0$. (c) Phase averaged time series of $P(t)$ and $\epsilon(t)$. Shaded region represents $\langle f \rangle_t \pm \sigma(f)$ where $f$ is $P(t)$ or $\epsilon(t)$. The pink vertical dashed line indicates the maximum of $\langle \epsilon \rangle_{\text{phase}}$, denoting the average time delay between $P(t)$ and $\epsilon(t)$.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure C1.** Schematic of forced (centre) plus six primary energetic (surrounding) Fourier modes in the SPO at $Re = 5.83$. Visualisations show distributions of $|\omega|$ at the same instance. The three-digit numbers on the visualisations indicate three components of wavevector $k_x k_y k_z$. Note that all the possible sign combinations $(\pm k_x, \pm k_y, \pm k_z)$ are gathered. Triangles denote the possible triad interactions. On top, we compare the isosurfaces of $|\omega_X + \omega_Y|$ with $|\omega|$ of the full flow. Here, $\omega_X$ and $\omega_Y$ denote the vorticity of the forced and the primary energetic modes, respectively.</figcaption></figure>

We plot the time series of the energy of the forced and six primary modes of the SPO at $Re = 5.83$ in Fig. C2. Although the energy $E_{110}(t)$ of the forced mode dominates, which is approximately equal to the total energy $E(t)$, we observe a distinctive difference between $E(t)$ and $E_{110}(t)$ when the primary scale energies are excited. By summing up the contributions of these seven modes, we obtain $E_{X+Y}$ shown in Fig. 3 (a). We also note that there are fast oscillations in $E_{100}(t)$ and $E_{010}(t)$. However, these two modes are compensated with each other, and such rapid dynamics are not visible in $E_{X+Y}(t)$ [Fig. 3 (a)]. This observation explains why there are no fast oscillations in the time series shown in Fig. 4 (b).
<figure><img src="image_12.png" /><figcaption>**Figure C2.** Time series of energy of the forced and six primary energetic modes in the SPO. $E_{k_x k_y k_z}$ denotes the energy summed-up for the wavevectors ($\pm k_x, \pm k_y, \pm k_z$). Total energy $E(t)$ is also shown for reference.</figcaption></figure>

## Appendix D. Linear stability analysis of the two-equation model

In this appendix, we show the results of the linear stability analysis of the fixed points of (6). There are two kinds of fixed points: namely,

$$ \overline{\mathbf{X}}_1 = \left( \frac{F}{\nu K_X^2}, 0 \right) \quad \text{and} \quad \overline{\mathbf{X}}_2 = \left( \frac{\nu K_Y^2}{A}, \pm \frac{1}{A} \sqrt{AF - \nu^2 K_X^2 K_Y^2} \right), \quad (D.1) $$

where $\overline{\mathbf{X}} = (\overline{X}, \overline{Y})$. Note that the fixed points $\overline{\mathbf{X}}_2$ exist only for $\nu < \sqrt{AF}/K_X K_Y$. The perturbation $(x, y)$ in the vicinity of the fixed points $(\overline{X}, \overline{Y})$ obeys

$$ \frac{dx}{dt} = -A(\overline{Y}^2 + 2\overline{Y}y) - \nu K_X^2 (\overline{X} + x) + F, \quad (D.2) $$

$$ \frac{dx}{dt} = +A(\overline{X}\overline{Y} + \overline{X}y + x\overline{Y}) - \nu K_Y^2 (\overline{Y} + y), \quad (D.3) $$

where we have neglected second-order terms $x^2, y^2$, and $xy$. The Jacobian matrix is then expressed as

$$ \mathbf{J} = \begin{pmatrix} -\nu K_X^2 & -2A\overline{Y} \\ A\overline{Y} & A\overline{X} - \nu K_Y^2 \end{pmatrix}, \quad (D.4) $$

whose eigenvalues are

$$ \lambda = -\frac{1}{2} [-A\overline{X} + \nu(K_X^2 + K_Y^2)] \pm \frac{1}{2} \sqrt{[A\overline{X} + \nu(K_X^2 - K_Y^2)]^2 - 8A^2\overline{Y}^2}. \quad (D.5) $$

The eigenvalues for $\mathbf{X}_1$ are

$$ \lambda_1^{(\mathbf{X}_1)} = -\nu K_X^2, \quad \lambda_2^{(\mathbf{X}_1)} = \frac{AF}{\nu K_X^2} - \nu K_Y^2, \quad (D.6) $$
<figure><img src="image_13.png" /><figcaption>**Figure E1.** Time series of (a) forcing coefficient $F(t)$ and (b) scale coefficients $K_\alpha^2(t)$ in the DNS of the SPO.</figcaption></figure>

which are both negative for $\nu > \sqrt{AF}/K_X K_Y$. Therefore, $\bar{X}_1$ is stable for $\nu > \sqrt{AF}/K_X K_Y$, and a pitchfork bifurcation takes place at $\nu = \sqrt{AF}/K_X K_Y$. Then, for $\nu < \sqrt{AF}/K_X K_Y$, $\bar{X}_2$ exists, which is stable irrespective of $\nu$ because the eigenvalues are

$$ \lambda_{1,2}^{(\mathbf{X}_2)} = -\frac{\nu K_X^2}{2} \pm \frac{\sqrt{-8AF + \nu^2 K_X^2 (1 + 8K_Y^2)}}{2}. \qquad (D.7) $$

## Appendix E. Detailed procedure of the parameter fitting

This appendix discusses the detailed procedure of the parameter fitting (10) of the three-equation model (7). Figure E1 (a) shows the time evolution of the forcing coefficient defined by

$$ F(t) \equiv \frac{P}{\sqrt{2E_X}}. \qquad (E.1) $$

We estimate the model parameter $F = 0.7$, since the time average $\langle F(t) \rangle_t = 0.696$. The periodic drops of $F(t)$ are associated with a phase-desynchronisation between the forcing and the forcing-induced velocity field, $u_X$.
<figure><img src="image_14.png" /><figcaption>**Figure E2.** Time series of (a) energy transfer terms $T_\alpha(t)$ and (b) transfer coefficients $A_i(t)$ in the DNS of the SPO.</figcaption></figure>

We also compute the scale factors

$$K_{\alpha}^{2}(t) \equiv \frac{\epsilon_{\alpha}}{2\nu E_{\alpha}} \quad (\alpha \in \{X, Y, Z\}). \qquad (E.2)$$

Figure E1 (b) shows their temporal evolutions. The forced scale factor $K_X^2(t) = 2$ is constant, since it corresponds to $k_f = (\pm 1, \pm 1, 0)$ mode. On the other hand, $K_Y^2(t)$ and $K_Z^2(t)$ fluctuate, reflecting the competition of different Fourier modes in these scales. We estimate the model parameters by $K_Y^2 = 5$ and $K_Z^2 = 15$, since $\langle K_Y^2(t) \rangle_t = 4.97$ and $\langle K_Z^2(t) \rangle_t = 15.4$, respectively.

To obtain rough estimates of the scale local coefficients $A_1$ and $A_2$, we compute the average energy transfer rate from X to Y and Y to Z while ignoring the scale non-local interactions by setting $A_3 = A_4 = 0$. In this way, the energy transfer terms (9) of the energy equation (8) of the model are approximated by

$$ \begin{aligned} T_X(t) &\approx -A_1XY^2, \\ T_Y(t) &\approx +A_1XY^2 - A_2YZ^2, \\ T_Z(t) &\approx + A_2YZ^2. \end{aligned} \qquad (E.3) $$
Minimal model of quasi-cyclic behaviour in turbulence driven by Taylor-Green forcing25

Figure E2 (a) shows their time series by the DNS of the SPO. $T_X(t) < 0$ supports the energy cascade picture; the forced scale $X$ is transferring energy to smaller scales ($Y, Z$) on average. Similarly, $T_Y(t), T_Z(t) > 0$ means that these smaller scales receive energy from the larger scales. We then evaluate the time-dependent coefficients,

$$A_1(t) \approx -\frac{T_X}{XY^2} = -\frac{1}{2\sqrt{2}} \frac{T_X}{\sqrt{E_X}E_Y}, \qquad (E.4)$$

$$A_2(t) \approx \frac{T_Z}{YZ^2} = \frac{1}{2\sqrt{2}} \frac{T_Z}{\sqrt{E_Y}E_Z}. \qquad (E.5)$$

Again, we neglect the scale non-local interactions ($A_3 = A_4 = 0$) in these expressions. The result is shown in Fig. E2 (b), and we estimate $A_1 = 0.4$ and $A_2 = 4$ as the model parameters from the time-averaged values $\langle A_1(t) \rangle_t = 0.440$ and $\langle A_1(t) \rangle_t = 4.04$, respectively.

The above argument allows us to determine the model parameters in (10). The non-local interaction coefficients $A_3$ and $A_4$ are left to be determined. In § 3.2, we vary these two parameters to investigate the model properties.