# Nonlinear phase synchronization and the role of spacing in shell models

L. Manfredini, Ö D Gürcan

► To cite this version:

L. Manfredini, Ö D Gürcan. Nonlinear phase synchronization and the role of spacing in shell models. Physical Review E, 2026, 113 (1), pp.015101. 10.1103/2vxp-1k2t. hal-05465301

HAL Id: hal-05465301

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

Submitted on 19 Jan 2026

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
# Nonlinear phase synchronization and the role of spacing in shell models

L. Manfredini, Ö. D. Gürcan

Laboratoire de Physique des Plasmas, CNRS, Ecole Polytechnique, Sorbonne Université,
Université Paris-Saclay, Observatoire de Paris, F-91120 Palaiseau, France

A shell model can be considered as a self-similar chain of interacting triads, where each triad can be interpreted as a nonlinear oscillator that can be mapped to a spinning top. Investigating the relation between phase dynamics and intermittency in such a chain of nonlinear oscillators, it is found that synchronization is linked to increased energy transfer. In particular, our results indicate that the observed systematic increase of intermittency, as the shell spacing is decreased, is associated with strong phase alignment among consecutive triadic phases, facilitating the energy cascade. It is shown that while the overall level of synchronization can be quantified using a Kuramoto order parameter for the global phase coherence in the inertial range, a local, weighted Kuramoto parameter can be used for the detection of burst-like events propagating across shells in the inertial range. This novel analysis reveals how locally phase-locked states are associated with the passage of extreme events of energy flux. Applying this method to helical shell models (i.e. for a class of helical interactions that couple the two helicities in a non separable topology) reveals that a reduction in phase coherence correlates with suppression of intermittency. When inverse cascade scenarios are considered using two different shell models including a non local helical shell model, and a local standard shell model with a modified conservation law, it was shown that a particular phase organization is needed in order to sustain the inverse energy cascade. It was also observed that the PDFs of the triadic phases were peaked in accordance with the basic considerations of the form of the flux, which suggests that a triadic phase of $\pi/2$ and $-\pi/2$ maximizes the forward and the inverse energy cascades respectively.

## I. INTRODUCTION

One of the greatest challenges in modern physics, especially concerning statistical dynamics, is characterizing the behavior of complex systems far from thermodynamic equilibrium. These systems are not only of fundamental intellectual interest, due to the rich complexity of the dynamics that they exhibit, but also for their relevance in a wide range of engineering and environmental applications.

Many natural systems exhibit non-Gaussian fluctuations across a broad range of spatial and temporal scales, implying the presence of extreme fluctuations or rare events, whose dynamical properties and statistical signatures are essential to describe the overall behavior of the system. Examples include seismic activity, financial market crashes and extreme weather events [1–3]. An archetypical example of such complex, open system is the problem of fully developed turbulence: an out of equilibrium state of fluid or plasma flow, displaying a mixture of hierarchical self organization and chaotic behavior across a large range of scales, maintained in a statistically stationary state by continuous energy injection (e.g. at large scales) and dissipation (e.g. at small scales).

Between forcing and dissipation, turbulent systems can develop an inertial range, spanning many decades, exhibiting self-similar properties heuristically captured by the Richardson cascade picture, and the power law spectrum of Kolmogorov [4].

Since direct numerical simulations (DNS) of turbulent flows over a wide range of scales were historically too costly, various reduced models were developed in the past [5]. These models are capable of reproducing various statistical features of the full system, such as the power

law distribution of energy as a function of scale, and the scaling of higher-order structure functions. Among those, shell models [6, 7] stand out in their simplicity and in their ability to mimic certain statistical aspects of different types of turbulent systems including rotating turbulence, passive scalar advection, thermal convection etc. [8–10]. In the context of plasma physics, shell models have been used to describe magnetohydrodynamic (MHD) turbulence [11, 12], and drift waves in fusion plasmas [13, 14].

Probably the most intriguing – perhaps unexpected – feature of shell models, is their ability to produce the anomalous scaling of the structure functions similar to fully developed turbulence [15–17]. This suggests that shell models can be a useful tool to capture the statistics of extreme events [18, 19], which are manifested as burst like structures, related to finite time singularities or instantons [20, 21], playing an important role in shaping the underlying intermittency characteristics.

Shell models reduce the complexity of the nonlinear energy transfer governed by triadic interactions drastically, while preserving the key properties of the underlying system such as its scale invariance, the quadratic nature of the nonlinearity and the conservation laws - like Energy and Helicity in the case of three dimensional turbulence. This extreme reduction in the number of degrees of freedom, is a consequence of dividing the Fourier space into a set of geometrically spaced shells and retaining the nonlinear interactions that involve only consecutive shells, thus replacing the turbulent system with a one dimensional chain of connected, self-similar triads.

In order to reintroduce some spatial structure in cascade models, various more sophisticated constructions have been proposed, including nested polyhedra [22],
logarithmic discretization (LDM) [23], spiral chains [24] and logarithmic lattices [25–27]. Most of these models are self-similar by construction, and maintain a constant number of triads across different scales. However, a puzzling feature of these spatially-extended “shell” models is the apparent suppression of intermittency observed in their statistics. Similar loss of intermittency has also been observed in numerical simulations of fractally-decimated 3D turbulence [28] even keeping an overwhelming majority of the modes and dropping only a small percentage.

The effect of a fractal decimation of Fourier modes has been explored also in simpler systems, such as the Burgers equation, where a disruption of the average triadic phase coherence is correlated with a reduction of the intermittency, normally associated with shock formation [29, 30]. This suggests a picture, in which phase alignment between triads enhances energy transfer resulting in strong dissipation events [31, 32]. Extending this analysis to the Navier Stokes equations requires understanding the self-organization of triadic phases in fully developed turbulence, and the definition of a clear order parameter that identifies the interactions involved in the collective behavior responsible for intermittency, which remains a challenge. Part of the difficulty involves dealing with the huge number of degrees of freedom represented by the convolution sum in such a system. One possible way forward may be to reformulate turbulence as a network of Fourier modes [33–35], which allows us to consider reduced models, such as shell models, as reductions of a complete grid into a network of groups of nodes, in this case represented by shells. Although the individual complex phases of the dynamical variables of a shell model have no physical significance, the fact that the triadic phase coherence along the chain regulates the cascade process, suggests that triadic phase coherence among the underlying Fourier modes, lumped together to form the shells, is related to the resulting shell-to-shell energy transfer.

Thus, to clarify the relationship between intermittency and phase coherence in fully developed turbulence, we investigate this connection in shell models that preserve the same invariants as the Navier-Stokes equations and, in some cases, exhibit a level of intermittency comparable to that of 3D turbulence. Building on our earlier observation that reducing inter-shell spacing in standard shell models increases intermittency considerably [36], we demonstrate a systematic link between triadic phase synchronization (sync) and enhanced intermittency in shell models, consistently with the observations in studies of the Burger’s equation cited above.

Notice that, a shell model can be interpreted as a one dimensional chain of connected triads, each representing a “nonlinear oscillator”. Consecutive triads only differ by a scaling factor related to the inter-shell spacing. This means that by decreasing the inter-shell spacing, (while keeping the conserved quantities unchanged) both the number and the similarity of the oscillators would be

increased making it potentially easier to form a phase coherent state. To quantify global phase coherence, and its relation to different levels of anomalous scaling of the structure functions, we employ a global Kuramoto order parameter, which allows us to show that higher intermittency correlates with enhanced phase coherence.

In addition, we observe that phase coherence correlates with the existence of coherent structures in the energy cascade, appearing as localized pulses in shell models. In order to clarify the role of these localized events, we introduce a novel local order parameter –a weighted Kuramoto parameter– which enable us to observe burst-like events propagating across shells in the inertial range and analyze the alignment of triadic phases in forming an “open channel” for energy transfer during these events. We show that extreme energy bursts are consistently preceded by rapid growth in local phase coherence and that more intermittent models display stronger and longer-lasting phase-locked states, thereby facilitating larger energy transfer.

Beyond the usual linear chain, we also analyzed helical shell models, which introduce the complexity of chiral interactions, fundamental in describing transitions in turbulent cascades [37]. Helical shell models consist in two sequences of helical modes that, depending on the interaction class, may be topologically intertwined by the triads that connect them. In a particular interaction class, a reduction of the inter-shell spacing leads to a disruption of the intermittency, offering us a framework to test our diagnostic tools for quantifying phase coherence in this opposite scenario.

Finally, we consider inverse cascades in shell models, which –apart from particular examples– are usually suppressed by equipartition solutions. We consider these special cases that actually work, in order to study the role of phase coherence and to compare it with the forward cascade scenario. The results reveal that, in the forward and inverse cascade regimes, the PDFs of the triadic phases exhibit peaks at values of $\pi/2$ and $-\pi/2$ respectively, which maximize the flux term at fixed amplitude. Notably, in the case of intermittent cascades the peak at $\pi/2$ is much more marked, indicating that coherent events are clearly favored by the system; in contrast, for inverse cascades, the distribution is symmetric. This results are in parallel and complementary with recent results focusing on the relation between phase dynamics and the energy cascade direction [38].

The rest of the paper is organized as follows.

In section II we introduce the shell model framework, focusing on the dynamical equations for phases and amplitudes IIA, discussing role of phase alignments in the energy flux IIB and presenting the phase statistics for the standard forward intermittent cascade. Next in II C, we introduce the Kuramoto order parameters to quantify phase coherence in shell models, highlighting how local sync is related to extreme flux events.

Section III addresses intermittency, in particular investigating how the enhancement of the intermittency
corrections, by decreasing the shell spacing g, is associated with different statistics of the sync order parameters III A. We further characterize the dynamics of burst events - responsible for intermittency - as sync events propagating through shells in cascade process.

In Sec. IV, we extend the analysis to helical shell models, exploring how sync properties vary with shell spacing, focusing on a helical shell model exhibiting suppression of intermittency. Section V is devoted to the role of phase organization in driving inverse cascades. In particular we present an inverse energy cascade for both an elongated helical shell model and a local shell model.

Finally, Section VI summarizes conclusions, while Appendix A provides numerical details.

## II. SHELL MODEL AS A DYNAMICAL SYSTEM: EQUATIONS FOR PHASES AND AMPLITUDES

Motivated by observations of statistical isotropy and scale invariance of the turbulent cascade of fully developed turbulence without obvious sources of anisotropy, the shell model reduction consists in modeling the velocity field as a one-dimensional set of complex variables representing different scales, corresponding to a set of logarithmically spaced wave numbers $k_n = k_0 g^n$, where it is customary to take the spacing factor $g = 2$.

The evolution equation for the shell variable $u_n$, for the GOY (Gledzer-Ohkitani-Yamada) model, which is the standard example for shell models, is usually written as

$$ \frac{du_n}{dt} = i [ak_{n+1}u_{n+2}^*u_{n+1}^* + bk_n u_{n-1}^*u_{n+1}^* + ck_{n-1}u_{n-1}^*u_{n-2}^* + d_n u_n + f_n], \quad (1) $$

where the terms in the square brackets represent the reduction of the convolution sum to a local triadic interaction between three consecutive shells, $d_n = \nu k_n^{2p} + \mu k_n^{-2q}$ is the general form of (hyper/hypo)-viscosity, and $f_n$ denotes external forcing, localized at a shell $n_f$, or a small set of consecutive shells around it, corresponding to the typical energy injection scale $1/k_{n_f}$.

Here the coefficients a, b, c are chosen in such a way that they preserve typical inviscid invariants of the underlying physical system. In the case of 3D Navier-Stokes turbulence, it is the energy and the helicity that are conserved, which are defined for the shell models respectively as

$$ E = \sum_n |u_n|^2, \quad H = \sum_n (-1)^n k_n |u_n|^2. \quad (2) $$

Note the peculiar aspect of the GOY model and its variants that positive and negative helicities correspond to even and odd shells respectively.

Enforcing the conservation of these two quantities one obtains, $a = g+1$, $b = g^{-1} - g$ and $c = -1 - g$, as a function of the inter-shell spacing g. It is also common to factor out $ak_{n+1}$ in the nonlinear term of Eq. (1) so that the other two coefficients corresponding to b and c can be written as $(g^{-2} - g^{-1})$ and $-g^{-3}$ respectively. The common choice of $g=2$ would then give $-1/4$ and $-1/8$ for the second and third coefficients respectively. Another common choice in the literature is $g = \varphi$ (the golden ratio), which gains significance, from the derivation of shell models based on a decimation of Fourier space, which requires the existence of triadic interactions between the characteristic wave numbers of the shells (specifically solution to $k_{n-1} + k_n + k_{n+1} = 0$, which imposes the geometrical constrain $g \le \varphi$). While there are some examples of studies of shell models varying their nonlinear coupling coefficients [39, 40], most of these studies are carried out in such a way that varying the coefficients change the conservation law, and hence the effective dimensionality of the underlying physical model, instead of keeping the conservation laws fixed and varying the logarithmic spacing for a fixed model. Exploring different shell spacings is also important for making the connection to a perspective based on Fourier space decimation, such as log-lattices [26], LDM [23], or helical shell models with different types of elongations [41].

Note that a standard shell model with $N$ nodes has $N-2$ interacting triads. Since the equation for a single unconnected triad is equivalent to that of a spinning top, giving rise to nonlinear oscillations that can be represented in terms of Jacobi elliptic functions, each triad can be considered as a nonlinear oscillator. Furthermore since the shell model consists of self-similar triads that are scaled by the same factor, the consecutive nonlinear oscillators that are coupled are also self-similar, maintaining the same scaling factor. The energy flux through such a chain of triads then requires some sort of synchronization of the phases of those consecutive triads. Decreasing the shell spacing is thus key to understanding the synchronization behavior of shell models. Since it simultaneously allows us to increase $N$, considering the large $N$ limit, important for synchronization, and reduce the ratio of the coefficients of two consecutive triads allowing them to have closer natural nonlinear frequencies. Note that the usual GOY model with $g=2$, covers roughly seven decades in $k$-space with $N=24$, and it is not feasible to increase that number substantially without decreasing $g$, since the power law solution would start going under machine precision levels.

Thus, varying the $g$ parameter, while increasing $N$ in order to cover the same wave-number range, allows us to modify the affinity of the triadic network of the shell model to synchronization (i.e. synchronizability). Since the intermittency in shell models is a consequence of its time evolution, one would expect this to have significant impact on intermittency as well. Hence, we present below, a detailed study introducing specific methods for observing various aspects of local and global synchro-
nization, and the implications of these observations for the turbulent cascade as well as intermittency, mainly by varying the g parameter to control its sychronizability.

A. Phase Dynamics in Shell Models

In order to clarify the relevance of phase organization, we substitute $u_n = U_n e^{i\theta_n}$ in Eq. (1), and separate the system into an equation for the amplitudes

$$ \frac{dU_n}{dt} = ak_{n+1}U_{n+2}U_{n+1} \sin \varphi_{n+1} + bk_n U_{n+1}U_{n-1} \sin \varphi_n \\ + ck_{n-1}U_{n-1}U_{n-2} \sin \varphi_{n-1} \\ + d_n U_n + Re \left( f_n e^{-i\theta_n} \right) \quad (3) $$

and a corresponding equation for the phases $\theta_n$, which can be written as:

$$ \frac{d\theta_n}{dt} = ak_{n+1} \frac{U_{n+1}U_{n+2}}{U_n} \cos \varphi_{n+1} + bk_n \frac{U_{n-1}U_{n+1}}{U_n} \cos \varphi_n \\ + ck_{n-1} \frac{U_{n-1}U_{n-2}}{U_n} \cos \varphi_{n-1} \\ + Im \left( \frac{f_n}{U_n} e^{-i\theta_n} \right). \quad (4) $$

Note that the triadic phases defined as $\varphi_n = \theta_{n-1} + \theta_n + \theta_{n+1}$, govern both the amplitude and the phase evolution, as they determine in particular the signs of each nonlinear term in Eqs. (3) and (4). This means that the evolution equation are invariant under the period three transformation

$$ \left\{ \begin{array}{ll} \theta_{3n-1} & \rightarrow \theta_{3n-1} - \delta \\ \theta_{3n} & \rightarrow \theta_{3n} + 2\delta \\ \theta_{3n+1} & \rightarrow \theta_{3n+1} - \delta \end{array} \right. $$

which leaves the triadic phases invariant for any constant $\delta$. The consequence of this $U(1)$ symmetry is that single shell phases $\theta_n$ are not uniquely determined by the triadic phases, and as also numerically confirmed, they remain uniformly distributed (i.e. “random”). Thus, the appropriate quantities, whose dynamics and statistics are physically relevant are the triadic phases $\varphi_n$, rather than single shell phases $\theta_n$.

We can write the equations for the triadic phases by adding Eq. (4) for three consecutive shells, which gives:

$$ \frac{d\varphi_n}{dt} = a\Lambda_{n+1}^n \cos \varphi_{n+1} + b\Lambda_n^n \cos \varphi_n + c\Lambda_{n-1}^n \cos \varphi_{n-1} \\ + a\Lambda_n^{n-1} \cos \varphi_n + b\Lambda_{n-1}^{n-1} \cos \varphi_{n-1} + c\Lambda_{n-2}^{n-1} \cos \varphi_{n-2} \\ + a\Lambda_{n+2}^{n+1} \cos \varphi_{n+2} + b\Lambda_{n+1}^{n+1} \cos \varphi_{n+1} + c\Lambda_n^{n+1} \cos \varphi_n \quad (5) $$

where

$$ \Lambda_{\ell}^{n} \equiv k_{\ell} \frac{U_{\ell-1} U_{\ell} U_{\ell+1}}{U_{n}^{2}}. $$

Note that Eq. (5) involves the triadic phases of five consecutive triads from $n-2$ to $n+2$, which in turn include contributions from seven consecutive shells from $n+3$ to $n-3$, making this coupling longer range compared to the shell amplitude equation, which involves only the neighboring triadic phases (i.e. $\varphi_j$ with $j = n-1, n, n-1$) and contributions from shells from $n-2$ to $n+2$. Note that these couplings can also be interpreted (and generalized to more complex interactions) within the network representation of the shell models, see e.g. Ref. 34, where the evolution of a given triad can be computed as a sum over all the triads connected through shared nodes in the network topology.

B. The Role of Triadic Phase on Energy Flux

The phase dynamics of a shell model with $N$ shells can therefore be considered as being equivalent to that of a chain of $N-2$ nonlinear oscillators, each representing a triad, where the coupling terms are determined by a combination of the $N$ self-consistently evolving amplitudes and the interaction coefficients. When forcing is applied near the first few shells and dissipation occurs at high-$n$, the chain of oscillators is constrained to accommodate the energy flux from low to high $n$. In particular if the phases are initially distributed in such a way that there is no initial net flux, the amplitudes will grow until the phases will reorganize themselves in order to accommodate the spectral flux implied by forcing and dissipation, basically because of the free energy associated with such a non-uniform energy distribution. This means that the energy injection at low-$n$ drives the phases to dynamically self-organize in order to sustain the cascade mechanism, possibly reaching an eventual steady state.

In order to illustrate this steady state, let us consider the form of the flux in a shell model. The energy cascades to small scales due to nonlinear interactions, and the energy transfer in a given triad centered at shell $n$ is defined as a local spectral energy flux, expressed in terms of the three shell correlation $\chi_n = U_{n-1}U_nU_{n+1}e^{-i\varphi_n}$ as:

$$ \Pi_n^E = 2 \mathrm{Im} (ak_{n+1}\chi_{n+1} - ck_n\chi_n). \quad (6) $$

Using $a = (1+g)$ and $c = -(1+g^{-1})$ and the power law solution $|\chi_n| \propto k_n^{-1}$, we get

$$ \Pi_n^E = 2 (g^{-1} + 1) (g \sin \varphi_{n+1} + \sin \varphi_n), $$

which is maximized for $\varphi_n = \varphi_{n+1} = \pi/2$. This means that the Kolmogorov cascade with the correspondingly power law spectrum and a constant energy flux is maximized if the triads are aligned at $\varphi_n = \varphi_{n+1} = \pi/2$.
Note that since $c$ is negative, this choice maximizes the forward flux regardless of the choice of $|\chi_n|$ as well. This allows the system to deviate locally from the Kolmogorov solution. Similarly an inverse cascade is implied if $\phi_n = \phi_{n+1} = -\pi/2$ regardless of the choice of the power law.

This means that the particular arrangement of the non-linear dynamics of the consecutive triads is important in determining both the direction and magnitude of energy transfer through the “boundary” between the two triads. In the following, we try to highlight this connection, considering different types of shell models and studying the role of triadic phases, and their statistical distributions in the cascade process.

### 1. Forward cascade and intermittency

In a shell model with local interactions, that conserves the inviscid invariants typical of 3D turbulence, corresponding to Eq. (1), and with the standard shell spacing $g=2$, one gets a forward energy cascade, and an energy spectrum consistent with the Kolmogorov $-5/3$ scaling together with a level of intermittency, somewhat consistent with experiment and direct numerical simulations [15–17].

The multifractal corrections, implied by this intermittent behavior (see, e.g., Refs. 15, 16, 19, 42), violating Kolmogorov’s original self-similarity hypothesis, manifest themselves as non-gaussian and scale-dependent distributions of velocity fluctuations in the inertial range, with fatter tails at smaller scales. For shell models, this is undoubtedly a consequence of the time evolution, since these models lack any information about spatial distribution of singularities. Note that, rare events in the time evolution of shell models are linked to sudden bursts of flux that goes through the chain of local triads that constitute them. Those events require both large amplitudes and some kind of local phase synchronization so that the burst can propagate across the chain in the correct direction.

More generally, the degree of locality of the cascade processes in shell models, can be investigated through multi-shell correlation functions [43], or via the statistical properties of the Kolmogorov multipliers, defined as $w_n = \frac{U_n}{U_{n-1}}$, and the triadic phases $\phi_n$, as done in Refs. 44–46, investigating whether or not the probability distribution functions of those quantities are universal in the inertial range (i.e. independent of the shell index). It is also worth noting that the correlation of the Kolmogorov multipliers of two different shells decays faster than that of the triadic phases [45]. This is significant, because it suggests that the *longer range* interactions implied by the triadic phase equation [i.e., Eq. (5)] introduce nontrivial phase correlations that are particularly relevant for the statistics of extreme events causing intermittency.

In Figure 1, we present the results for a standard GOY model of three dimensional Navier-Stokes turbu-

<figure><img src="image_2.png" /><figcaption>Figure 1. Plot of the average energy spectrum as a function of the wave number obtained for the 3D GOY model with inter-shell spacing $g=2$, driven by large-scale forcing. The bottom-left inset shows the average spectral energy flux as a function of the shell wave number $k_n$. In the top-right inset the normalized PDFs of the triad phases $\phi_n$ are shown for different shells corresponding to the inertial range. The shell number is represented by different colors labeled by the continuous colorbar.</figcaption></figure>

lence, forced at large scales, or small $n$, which is extensively studied in the past [15, 16, 42]. The dynamics of the model exhibit a forward energy cascade with clear inertial range, with a spectral slope $|u_n|^2 \sim k_n^{-2/3}$, which is the shell model equivalent of $E(k) \propto k^{-5/3}$. In the inset plot, the normalized PDFs of the triad phases $\phi_n$ are displayed for different shells in the inertial range, showing that they collapse to a more or less scale invariant universal distribution similarly to Kolmogorov multipliers [45].

This distribution exhibits a distinct peak at the fixed point $+\pi/2$ consistent with a direct energy cascade. As we will discuss later in Sec. III, reducing the inter-shell spacing progressively disrupts scale invariance and enhances intermittency corrections.

### C. Sync and Kuramoto Order parameter

As discussed in Section II A, the cascade in a shell model happens via self-organization of the triadic phases $\phi_n$ of the consecutive triads that constitute it, in order to accommodate the flux that goes through them. For the forward cascade, this causes the triadic phases to cluster around $\pi/2$, and for the standard GOY model, this implies a scale-invariant universal distribution of their statistics in the inertial range. In the usual scenario, such statistics underlies intermittent dynamics, corresponding to rare but extreme events.
<figure><img src="image_3.png" /><figcaption>Figure 2. Normalized probability distribution functions of the Kuramoto order parameter $R$ and $\Phi$, for various shell models exhibiting different type of cascades: a classical GOY model displaying a forward energy cascade; a helical shell model with elongated triads driving an inverse cascade (see Sec. IV and V 1 for a detailed discussion and the definition of the Kuramoto parameter in the context of helical shell models); and the model investigated in Ref. 50 that drives an inverse cascade via local triads (see V 2 below).</figcaption></figure>

In order to quantify the level of coherence of the chain of triads that facilitates the energy cascade, we can define a global order parameter by considering the population of nonlinear oscillators represented by triadic phases, in order to explore the role of phase aligned triads in extreme flux events.

For a 1D system of interacting nonlinear oscillators, which are described by triadic phases, we can define the Kuramoto order parameter [47–49] as:

$$Re^{i\Phi} = \frac{1}{N_m} \sum_m e^{i\varphi_m}, \qquad (7)$$

where the sum $\sum_m$ is over all the triads in the inertial range, where a clear power low scaling is observed, independent of forcing and dissipation.

At a given time $t$, the complex phase of the Kuramoto parameter $\Phi(t) \in [0, 2\pi)$ measures the average phase of the population of triads that are considered, and the amplitude $R(t) \in [0, 1]$ quantify the phase coherence of the chain. A completely incoherent state corresponds to $R(t) \approx 0$, since all the triads are uniformly distributed on the unit circle, whereas, the case with all the triad phases being perfectly aligned corresponds to $R = 1$. Using the Kuramoto parameter for a set of triadic phases, was initially introduced in the context of the Burgers equation [30, 31, 51] in order to analyze the energy flux enhancement due to highly coherent states, related to collisions of shocks in physical space. Here, we apply it to shell models, for characterizing the dynamics of triadic phases in driving the underlying cascade and the formation of rare events that cause intermittency.

In Figure 2 we present the normalized PDF of both $R$, and $\Phi$ for the standard GOY model that drives a forward cascade, discussed above and for two other models

that are capable of generating inverse cascade (which will be discussed in more detail in Sec. V, since inverse cascade in shell models is a complicated issue). Note that the PDF of $R$ for the forward cascade peaks at around $R \approx 0.7$, it exhibits significantly higher level of coherence with respect to the two inverse cascade models, which peak at $R \approx 0.1$ and at $R \approx 0.25$ respectively. Regarding the normalized PDF of $\Phi$, a strong evidence of phase clustering is present for the forward cascade shell model showing a sharp distribution with a maxima at $\pi/2$, in contrast to the two models that result in inverse cascade displaying a broader, less pronounced distributions. Nevertheless, the two non-intermittent models display a weak clustering of the triadic phases around $-\pi/2$. However, this is so weak that it is not clearly visible for the helical inverse cascade model in the presented plot, as compared to the model of Tom and Ray [50].

## 1. Local synchronization and energy transfer

The clustering of triadic phases around $\pi/2$, and a significant value of the Kuramoto parameter, as in the forward cascade case shown in Fig. 2, provides an indicator of the global phase coherence in the entire inertial range. We can also analyze the local degree of synchronization within the triad chain and its correlation to energy flux at the corresponding shell $n$, to better understand the contribution of phase alignments in the cascade mechanism and its variations.

In order to establish the connection between coherent structures traveling across shells, extreme events, end the alignment or clustering of triadic phases through the inertial range, we propose the following local Kuramoto order parameter centered at the $n$-th shell:

$$R_n e^{i\Phi_n} = \frac{\sum_{t \in \Gamma(n)} w_t e^{i\varphi_t}}{\sum_{t \in \Gamma(n)} w_t}, \qquad (8)$$

where the summation $\sum_{t \in \Gamma(n)} w_t$ represents a shell-dependent weighted sum over a set of triads $\Gamma(n)$. A consistent choice, in order to quantify the local coherence is to restrict the summation to the set of interacting triads $\Gamma(n)$ which are connected to the central triad of the $n$-th shell. In this approach the weight of each triadic phase is equal to the number of times it appears on the right-hand side of the evolution equation for $\varphi_n$ [i.e., Eq. (5)]. While the weights of those who do not appear are set to zero. With this definition, a local coherent state ($R_n = 1$), corresponds to local phase locking, which implies a fixed point of the triadic phase evolution equation. In the case of a local shell model, this definition corresponds to a sum over a list of 5 consecutive triads with the weights given by $\{w_{n-2}, w_{n-1}, w_n, w_{n+1}, w_{n+2}\} \equiv \{1, 2, 3, 2, 1\}$.

As shown in Fig. 3, panel (a), the PDFs of the spectral energy flux $\Pi_n^E$, computed for different shells in the inertial range, display a typical form characterized by a pronounced peak, corresponding to a relatively small flux
value (with a mean value around 0.7 as can be seen in the inset of Figure 1), accompanied by heavy tails accounting for intense forward energy transfer events, which can reach values as large as 500 times the mean flux. This means that there are rare events that drive sudden and large energy transfer.

The panel (b) of the same figure, shows the PDFs of the local synchronization parameter $R_n$, computed using Eq. (8), for different shells. Notably, the distributions are scale-independent in the inertial range and imply a clear tendency of the dynamics to favor locally phase-coherent configurations over incoherent ones.

In panel (c), we consider the joint PDF of these two quantities for a given shell in the inertial range, in order to investigate if local synchronization correlates with energy transfer. The joint distribution shows that while typical flux events, close to the mean value are distributed across the whole range of $R_n$, rare but extreme flux events only occur for large values of $R_n$, where consecutive triads are partially-synchronized. In other words, while for small $R_n$ the flux events are distributed more or less symmetrically around the mean, an overwhelming majority of the large flux events, tend to have large $R_n$ values. In the same spirit, panel (d) shows the joint PDF between $\Phi_n$ and the local energy flux $\Pi_n^E$, illustrating that the extreme events are associated with a pronounced clustering of triad phases around $\pi/2$. These results show clearly that rare events go hand in hand with local synchronization of triadic phases in shell models. In the following we will study the effect of these rare events on a well-established measure of turbulence intermittency and how this behavior varies with the inter-shell spacing $g$, which changes both the synchronizability of the triadic interactions and the number of its elements.

### III. INTERMITTENCY

We have shown that extreme events in energy flux (and amplitudes $U_n$) are linked to triadic phase coherence. The statistical properties of these events determine the level of intermittency, manifesting itself in the deviation of the scaling of the structure functions, defined as $S_p(k_n) = \langle |u_n|^p \rangle$, from the scaling $S_p \propto k_n^{-p/3}$ predicted by Kolmogorov's self-similarity hypothesis. In the GOY shell model, intermittency corrections can be more accurately quantified by studying the scaling behavior of the three point correlations:

$$ \Sigma_p(k_n) = \langle |Im(\chi_n)|^{p/3} \rangle, $$

which eliminates period three oscillations of shell amplitudes [15]. The scaling exponents $\xi(p)$ of the p-th order structure function can be estimated either by performing a linear fit in a log-log plot of $S_p(k_n)$ or $\Sigma_p(k_n)$ versus $k_n$, or through a more refined method such as the Extended Self-Similarity (ESS) [52].

<figure><img src="image_4.png" /><figcaption>Figure 3. Individual and joint statistics of energy flux, and the amplitude and phase of the Kuramoto parameter. In panel (a), located at the top-left, the PDF of the local spectral energy flux $\Pi_n^E$ is shown, exhibiting the typical behavior where the distribution develops increasingly heavy tails for larger shell indices. In panel (b) at the top right, the PDF of the local Kuramoto parameter $R_n$ is presented, displaying a scale-invariant universal statistics with higher probability of coherent local $R_n$. Curve colors correspond to different shell numbers, as indicated by the continuous colorbar. In panel (c) at the bottom left, the joint PDF between $\Pi_n^E$ and $R_n$ is shown, which illustrates how extreme flux events are predominantly correlated with higher local synchronization events. In panel (d) at the bottom right, the joint PDF between $\Pi_n^E$ and $\Phi_n$ is shown, complementary to panel (c), that present the clustering of triad phases around $\pi/2$ which is associated with larger flux values.</figcaption></figure>

The ESS approach, estimates the scaling exponents $\xi_{ESS}(p)$, as the power low scaling of the structure function of order $p$, with respect to the third order one, i.e. the scaling $S_p(k_n) \sim S_3^{-\xi_{ESS}(p)} (\Sigma_p(k_n) \sim \Sigma_3^{-\xi_{ESS}(p)})$. This approach can reveal hidden scalings in turbulent structure functions (see Ref. 53 as an example in the shell model context). The improved accuracy of the ESS method has been explained as its ability to suppress sub-dominant corrections [54] and an interpretation of ESS within the multifractal formalism has been proposed in Ref. 55.

In the following analysis, we present different scaling exponents $\xi(p)$ that are extracted from numerical simulations of different shell models. In particular, for the GOY model we employ the three-point correlation $\Sigma_p(k_n)$, using the ESS procedure, restricting the fits to a range of scales corresponding to the constant flux window in the inertial range. Note that even though no hid-
<figure><img src="image_5.png" /></figure>

Figure 4. Scaling of the structure functions $S_p(k_n)$ (orders $p = 2, 3, \dots, 6, 7$) as a function of the shell wave-number $k_n$, shown for two representative values of the inter-shell spacing $g = 1.4$ and $g = 2$. The right panel ($g = 2$) presents the typical intermittency corrections, while the left panel ($g = 1.4$) presents an enhancement of intermittency. The insets display the PDFs of the triadic phases $\phi_n$ for different shells in the inertial range, with colors indicating shell index $n$ as shown in the colorbar.

The scaling exponent $\xi(p)$ is revealed by the ESS analysis (i.e. $\xi_{ESS}(p) \sim \xi(p)$); here, it helps refine our estimation of the intermittency.

## A. Shell Spacing Dependence

Although the widely adopted choice in literature, is to set the inter-shell spacing $g = 2$, having actual triadic interactions involving a Fourier grid structure with the wave-numbers given by $k_n$ requires an inter-shell ratio smaller than the golden ratio. It is therefore interesting to study how the intermittency corrections depend on the inter-shell spacing $g$. In a recent work [36], exploring a class of asymptotically self-similar shell models based on recurrent sequences -where the shell wave numbers $k_n$ are integers- we observed that the intermittency corrections increase for sequences characterized by smaller inter-shell spacings.

In order to study this phenomenon, here we provide a detailed characterization of the intermittency as a function of inter-shell spacing in regular shell models, in particular considering the role of triadic phase synchronization in the generation of rare and intense events of energy transfer that are responsible for the increased intermittency.

We first note that the intermittency corrections have a non-trivial dependence to the inter-shell spacing $g$ for the standard (GOY) shell model. Figure 5 shows the estimated scaling exponents $\xi(p,g)$ of the $p$-th order structure functions, for different values of the spacing, in the range $g = 1.1 - 2$, suggesting a progressive enhancement of intermittency, corresponding to a flattening of the scaling exponent curves as $g$ is decreased. These observations are consistent with previous studies, particularly those

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

Figure 5. The scaling exponent $\xi(p)$ of the $p$-th order structure function as a function of $p$, for different shell spacings. The dashed line represents the Kolmogorov $p/3$ scaling.

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

Figure 6. Scaling exponent $\xi_p(g)$ for structure functions of order $p$ (ranging from 4 to 10 in steps of 0.5), plotted as a function of the inter-shell spacing $g$. The marker colors correspond to different shell spacings, as in Fig. 5.

investigating the continuum limit of shell models, where $g \to 1$ leads to equations closely resembling the Burgers equation [56], known to generate shocks and an intermittent, bifractal anomalous scaling behavior [57]. Further insight is provided by Figure 6, which shows the dependence of the scaling exponents $\xi_p(g)$ as a function of the inter-shell spacing $g$. High order scaling exponents such as $p \gtrsim 8$ exhibit a sharp transition when $g$ decreases below a critical value $g_c \simeq 1.7$, suggesting a qualitative change in the scaling properties of the system. It seems that for $g \gtrsim g_c$ the scaling is roughly independent of the inter-shell spacing, while for $g \lesssim g_c$ the anomalous scaling exponents vary significantly with the inter-shell spacing, following an approximately linear scaling.

Note also how $g_c$ is just slightly above the golden ratio. While this may be a pure coincidence, there may be a link to nonlinear frequencies $\omega_n$, of the $n$-th triad being close
enough to interact more efficiently through a recurrence relation such as $\omega_n = \omega_{n-1} + \omega_{n-2}$.

## 1. Shell spacing and the level of synchronization

With the aim of highlighting the role of synchronization events in increasing intermittency corrections, we analyze the statistical properties of triadic phases across the inertial range. In Figure 4, we present the scaling of structure functions for two representative cases, $g = 2$ and $g = 1.4$, which are characterized by a standard and very high levels of intermittency, respectively. In the inset of the same figure, we present the PDFs of triadic phases for different shells within the inertial range. For the higher intermittency case ($g = 1.4$), the triadic phase distributions exhibit a pronounced scale-dependent behavior, becoming progressively narrower at smaller scales in the inertial range, a similar trend is also observed in the statistics of the Kolmogorov multipliers (not shown in the plot). This behavior is actually common for all the $g$ values below the critical value $g \lesssim g_c$ and becomes more pronounced as $g$ decreases. This suggests a tendency towards stronger phase alignment at higher wave numbers, consistent with the idea that the increase in intermittency is due to rare but large dissipation events at small scales.

In order to characterize the global level of synchronization across the inertial range, we consider the complex Kuramoto order parameter, defined in Eq. (7), and present the PDFs of the coherence parameter $R$ in Fig. 7 as a function of the inter-shell spacing $g$.

We note that, as the inter-shell spacing decreases, the PDF of the coherence parameter $R$ become increasingly skewed and a finite probability for fully synchronized events appears. These are states in which all the triads are aligned, and are associated with rare, extreme events in the forward energy flux. Remarkably, we observe that below the critical value $g \lesssim g_c$, the PDFs of $R$ are $g$-dependent, while for $g \gtrsim g_c$ they collapse to the same distribution, which is also in accordance with the behavior of intermittency corrections as a function of $g$. This suggests that the Kuramoto parameter, despite being only a measure of the phase coherence, can be used as a proxy for detecting different levels of intermittency.

The statistics of the average phase, represented by the phase $\Phi$ of the Kuramoto parameter is presented in Fig. 8. The PDFs of $\Phi$ as a function of $g$ are all centered around $\pi/2$, but with a width that decreases rapidly for $g \lesssim g_c$, indicating that for smaller $g$ the system has a stronger tendency to cluster around a common phase that maximizes the forward energy flux.

Not surprisingly, global synchronization goes hand in hand with local coherence. This can be seen clearly in Figure 9, which is equivalent to Fig. 3, but with an inter-shell spacing of $g = 1.5$, corresponding to strong intermittency. In Panel (a), the statistics of the local spectral energy flux $\Pi_n^E$ are presented, showing an increasing

<figure><img src="image_8.png" /><figcaption>Figure 7. Normalized PDFs of the Kuramoto order parameter $R$, shown for the different values of the inter-shell spacings $g$ that are considered.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 8. Normalized PDFs of the average phase $\Phi$, for different shell spacings $g$. Marker colors correspond to the different shell spacings, as in Fig. 7.</figcaption></figure>

occurrence and intensity of energy transfer along with a reduced probability of negative flux events as we move towards smaller scales. This behavior is evidenced by the progressively heavier tails of the distribution curves. The scale-dependent organization of triadic phases $\phi_n$ can be seen in the statistics of the local phase coherence parameter $R_n$, in panel (b), which show that smaller scales are associated with higher synchronization levels, and an increased probability of fully coherent states ($R_n \simeq 1$). These events are clusters of locally aligned triads with $\phi_n \simeq \pi/2$ corresponding to a fixed point of the phase dynamics.

The correlation between the energy transfer and the
<figure><img src="image_10.png" /></figure>

Figure 9. Individual and joint statistics of energy flux and the amplitude and phase of the Kuramoto parameter. In panel (a) on the top left plot, the PDF of the local spectral energy flux $\Pi_n^E$ is shown, exhibiting the typical behavior where the distributions becomes progressively more heavy-tailed for shells corresponding to smaller scales. This reflects the increasing occurrence and intensity of energy bursts as we move towards smaller scales in the inertial range. In panel (b) on the top right, the PDF of the local Kuramoto phase coherence $R_n$ is presented, where smaller scales display higher synchronization levels, with an increasing probability of local fully coherent states ($R_n \simeq 1$). Curve colors correspond to different shell numbers, as indicated by the continuous colorbar. In panel (c) on the bottom left, the joint PDF between $\Pi_n^E$ and $R_n$ is shown, which illustrates how extreme flux events are predominantly associated with high local synchronization events. In panel (d) on the bottom right, the joint PDF between $\Pi_n^E$ and $\Phi_n$ is shown, complementary to panel (c), reinforcing the observation that extreme flux events are associated with a marked clustering of triad phases around $\pi/2$. Note that this state, where the phases are locked, is highly efficient in transferring energy towards small scales.

local synchronization level along the triad chain is analyzed via the joint statistics of $\Pi_n^E$ and $R_n$, as well as $\Pi_n^E$ and $\Phi_n$, shown in panel (c) and (d) respectively. These plots demonstrate that local synchronized triad configurations - identified by $R_n \simeq 1$ and $\Phi_n \simeq \pi/2$ - are strongly correlated with extreme flux events.

## B. Burst events and temporal intermittency

So far, we have characterized the statistical properties of shell models as a function of the inter-shell spacing $g$, showing the relation between phase coherence and statistics of extreme flux events, leading to different levels of

intermittency.

In particular, we showed that the flattening of the scaling exponents $\xi(p)$ of the structure functions that is observed for decreasing values of $g$, is rather similar to what happens in Burgers equation with the formation of shocks. The equivalent process in shell models is the formation of burst-like solutions. These solutions have been addressed in the continuum limit analytically by assuming a zero background by Andersen et al. [56].

Similar burst events are also present in the fully developed turbulent state, coexisting with a background power law solution of the average spectral energy. They can be identified as coherent sequences of local maxima in the time evolution of the shell amplitudes $U_n$, propagating across shells. Such structures, sometimes referred to as *instantons*, and have been proposed as a key mechanism underlying anomalous scaling in shell models [20, 21].

A detailed characterization of the statistical properties of instantons as a function of the inter-shell spacing $g$ lies beyond the scope of the present work. This would represent an interesting future direction, particularly in the low-$g$ regime, characterized by stronger intermittency, where we observe also backward-propagation events, suggesting a scenario along the lines of a gas of colliding structures, rather than isolated, unidirectional bursts.

Nevertheless, continuing with the line of inquiry of this work, we present a phenomenological analysis of these burst like events with respect to the local synchronization properties of the system.

First we present how the dynamics of these events depend on the degree of intermittency, controlled by the inter-shell spacing $g$. In order to visualize these structures, Fig. 10 provides the temporal evolution of the compensated shell amplitudes, for different values of the inter-shell ratio $g = 2, 1.4, 1.2$. We can see rather frequent but less intense bursts for larger $g$, some of which may dissipate or annihilate within the inertial range, before reaching the smallest scales. In contrast, as $g$ is decreased, bursts become less frequent, more localized, and significantly more intense, allowing us to see their propagation across the whole inertial range. In this regime, a given shell spends most of its time in a quiescent state, interrupted by the passage of large-scale coherent bursts, which may also occasionally be accompanied by secondary events forming tree-like patterns.

In order to support our hypothesis that the dynamics of these bursts are associated with triadic phase synchronization, we also present the time evolution of the compensated amplitudes, the energy flux, and the local Kuramoto parameter $R_n$ in Figure 11, for inter-shell spacing $g = 1.5$. This demonstrates that the arrival of a burst corresponds to a rapid growth of local synchronization level, which precedes the extreme flux event. It also shows that as the flux event arrives at the smallest scales, all the shells are synchronized and they lose this synchronization slowly over time with the smallest scales becoming desynchronized first. In this sense, synchro-
<figure><img src="image_11.png" /><figcaption>Figure 10. Temporal evolution of the compensated shell amplitudes, $log(k_n^{1/3}U_n)$ along the one dimensional chain (shell index $n$ in the y-axis) shown for different inter-shell spacing $g = 2, 1.4, 1.2$. (top to bottom panels). The plots illustrate how the coherent events became more and more relevant as $g$ decreases.</figcaption></figure>

scales and propagates towards large scales.

One explanation of these observations is that the energy, being injected at large scales, tends to accumulate in low-$n$ shells, until the triadic phases can organize coherently, facilitating the release of this accumulated energy across the shells. This naturally starts at the large scales, and proceeds like a domino effect towards small scales. Therefore, as the burst propagates, the local Kuramoto parameter exhibits a sharp increase, corresponding to a local synchronization of the system. This is necessary to accommodate the passage of the large deviation in flux. After the passage of the burst, the system remains partially phase locked for some time, corresponding to a cluster of scales with $R_n \simeq 1$, maintaining an “open channel” for energy transfer (in terms of triadic phases). This synchronized configuration corresponds to a fixed point of the triadic phase equation, which provides a quasi-stationary state. It then gets desynchronized at the smallest scale first, probably because of the nonlinear time scales associated with the smallest scales are the shortest, so the decay to a desynchronized state in the absence of a synchronizing force is faster at small scales. Although we do not present the analogue of Fig. 11 for different values of the inter-shell spacing, we note that as $g$ decreases, the phase locked state appears to be more robust and stable, requiring a longer time for the system to desynchronize and return to a disordered phase configuration.

After the decay of the burst, energy injected at large scales starts to accumulate again in the low-$n$ shells, and

a new burst is eventually triggered, in a sort of cyclic dynamics of energy accumulation and release. (Interest- ing our observation echos, with the results of Ref. 58 in which avalanche-like behavior of seismic events, has been compared to energy release events and their inter-event times (waiting times). In our framework, the duration of the quiescent states — i.e., the time intervals during which energy accumulates at large scales — appears to be directly influenced by the inter-shell spacing $g$.

These observations suggest that the increased stability of phase-locked states — and consequently the longer duration of the so-called quiescent periods during which energy is accumulated in low-$n$ shells — may favor the occurrence of extreme events capable of releasing larger amounts of energy.

## IV. HELICAL SHELL MODELS AND SYNC

We have shown that the existence of phase coherent energy transport events account for the anomalous scalings of structure functions and therefore intermittency in shell models, and that the phase coherence happens through a mechanism of nonlinear synchronization. This means that the affinity of its triads to synchronization, is the key feature of a shell model that determine its intermittency. So far, we have characterized the *synchronizability* of only the simplest shell model: the classical GOY model with a one dimensional chain consisting in $N$ shells connected by $N-2$ triads.

However, the *synchronizability* of a system of triads (or, generally, a system of coupled oscillators) also depends on the way its elements are wired, which is represented by its topology. It is somewhat obvious that topological effects become relevant as we depart from a linear chain picture, as in the case of more complex model such as nested polyhedra models or log-lattices, where multiple elements interact at each scale, but even if we put aside those more complex examples, a minimal reintroduction of network topology appears in simple shell models if we try to incorporate helicity in a proper way. This happens if we consider two signs of helicity in each scale -instead of one sign in odd scales and another sign in even scales- as in helical shell models. These models are a natural generalization of the standard shell model, in which each scale is associated with two dynamical variables, $u_n^+$ and $u_n^-$, carrying positive and negative helicities, respectively. In this framework the two inviscid invariants are defined as:

$$E = \sum_n |u_n^+|^2 + |u_n^-|^2, \quad H = \sum_n k_n (|u_n^+|^2 - |u_n^-|^2). \quad (9)$$

This formulation results in two, one-dimensional sequences of shell variables, corresponding to the two helicity signs. The triadic interactions among those modes can be systematically recovered from Waleffe’s helical decomposition [59], where the nonlinear term in the Navier-
<figure><img src="image_12.png" /><figcaption>Figure 11. The left panels show colormaps of the temporal evolution of the compensated amplitudes, $log(k_n^{1/3}U_n)$, the energy flux, $log(\Pi_n^E)$, and the local Kuramoto parameter, $R_n$, across different shells (y-axis). The right panels present the temporal evolution of the same quantities during a burst event (unshaded region in the left panels), plotting distinct curves for different shells (as indicated by the colorbar) highlighting the local synchronization dynamics associated with large deviation events.</figcaption></figure>

Stokes equations is decomposed into distinct classes of helical interactions.

The general form of a helical shell model, omitting the dissipative and forcing terms, with arbitrary non-local triads connecting the three shells $n-l, n, n+m$, can be written compactly as:

$$ \frac{du_n^{s_0}}{dt} = s_0 (a_n u_{n+\ell+m}^{s_3*} u_{n+\ell}^{s_1*} + b_n u_{n-\ell}^{s_1*} u_{n+m}^{s_2*} + c_n u_{n-m}^{s_2*} u_{n-\ell-m}^{s_3*}) $$

with

$$ a_n = Q \frac{s_1}{s_0} \left( \frac{s_2}{s_0} k_{n+\ell+m} - k_{n+\ell} \right) $$

$$ b_n = Q \left( \frac{s_1}{s_0} k_{n-\ell} - \frac{s_2}{s_0} k_{n+m} \right) $$

$$ c_n = Q \frac{s_2}{s_0} \left( k_{n-m} - \frac{s_1}{s_0} k_{n-\ell-m} \right) $$

$s_3 \equiv s_0 s_1 s_2$, and where $Q = Q_{lm} (s_1/s_0, s_2/s_0)$ is the geometric factor of the triad class that is considered. The four triad classes can be defined as the 4 possible combinations of the signs of $s_1/s_0$ and $s_2/s_0$. Notice that the choice of local interactions corresponds to setting $l = 1, m = 1$, and results in $2N$ modes connected by $2N-4$ triads for each class. This represents - with the usual choice $g = 2$ - the helical shell model introduced in Ref. 60, investigated by considering the different classes of interactions separately.

However, the general formulation of the model allows us to consider arbitrarily elongated triads as investigated

by De Pietro, Biferale, and Mailybaev [61]. Moreover, one can also consider a right-hand side that constitutes a sum over different classes, each corresponding to different geometrical factors defined in terms of the shell index and the elongation of the triads [41].

Each of these choices correspond to a different network of triadic interactions, connecting two sequences of helical modes. For simplicity, we focus on models constructed from a single interaction class, and adopt the labeling SM1, SM2, etc. (see, e.g., Refs. 60, 61) following Benzi and De Pietro. There are qualitative differences among the networks represented by different classes. Note in particular that for certain models, the triadic interactions form two independent chains that do not interweave. We call those models *separable*.

This allows us to distinguish two types of topologies namely *separable* and *non separable*. The SM4 class, corresponds to $s_1/s_0 = s_2/s_0 = +$, is *separable* since it consists of homochiral triads of the type $(+, +, +)$ or $(-, -, -)$ where two helicities decouple and each helical variable evolves independently that conserve separately their corresponding positive-definite invariants, (i.e. $\sum_n k_n |u_n^+|^2$ and $\sum_n k_n |u_n^-|^2$), making it a candidate to drive an inverse cascade, which is in practice overwhelmed by equipartition. The other *separable* class is the SM1 class, corresponding to $s_1/s_0 = -, $ and $s_2/s_0 = -, $ involving heterochiral triads of the type $(-, +, -)$ or $(+, -, +)$, which naturally splits into two distinct chains of modes with alternating helicities, consisting effectively of two uncoupled GOY models.

The remaining classes SM3 and SM2, are *non separable*, since it is not possible to disentangle the triad chains
in these models into two distinct chains because of the way the interactions are coupled together. This means that for these classes the whole chain is interwoven, and the energy transfer across such a chain involves multiple interactions at each scale, which causes the synchronizability to be qualitatively different.

In this case, the right-hand side of the phase evolution equation involves contributions from up to seven different triads, in contrast to the five involved in the separable (single-chain) case, as detailed in Sec. II A. Therefore, in order to characterize the level of synchronization in a non separable chain, we define a weighted local Kuramoto parameter, constructed by summing the contributions of up to seven phases of its interacting triads.

More specifically, in a helical shell model, one can identify two distinct triads centered at the same shell, which can be written for the n-th shell as $u_{n-1}^{s_1}u_n^{s_0}u_{n+1}^{s_2} = \chi_n e^{i\varphi_n}$ and $u_{n-1}^{-s_1}u_n^{-s_0}u_{n+1}^{-s_2} = \bar{\chi}_n e^{i\bar{\varphi}_n}$. Note that the evolution equations are odd under the parity transformation, which implies that an organization of triad phase that maximizes the forward energy flux corresponds to the choice $\varphi_n = \pi/2$, $\bar{\varphi}_n = -\pi/2$. We then define the local weighted Kuramoto parameter as:

$$R_n e^{i\Phi_n} = \frac{\sum_{t \in \Gamma(n)} w_t e^{i\varphi_t} + \bar{w}_t e^{-i\bar{\varphi}_t}}{\sum_{t \in \Gamma(n)} w_t + \bar{w}_t}, \quad (10)$$

where the weights

$$\{w_{n-2}, w_{n-1}, w_n, w_{n+1}, w_{n+2}; \bar{w}_{n-2}, \bar{w}_{n-1}\bar{w}_n, \bar{w}_{n+1}, \bar{w}_{n+2}\}$$

are determined by the structure of the phase evolution equation that correspond to a set of interacting triads $\Gamma(n)$ which are connected with the central triad of the n-th shell $u_{n-1}^{s_1}u_n^{s_0}u_{n+1}^{s_2} = \chi_n e^{i\varphi_n}$. These weights, consistently with the previous definition (Eq. (8)), are equal to the multiplicity of each triad in the right-hand side of the evolution of equation for the triadic phase $\varphi_n$. In helical shell model, their values depends on the topology of the interaction class. For instance, for the triads of type SM3 the weights are $\{0, 1, 3, 1, 0; 1, 1, 0, 1, 1\}$, while for the separable class SM1, they are $\{1, 0, 3, 0, 1; 0, 2, 0, 2, 0\}$. Please note how these weights are consistent with the definition adopted for the single chain case, that we can express as $\{w_{n-2}, \bar{w}_{n-1}, w_n, \bar{w}_{n+1}, w_{n+2}\}$, since the SM1 reduces to two uncoupled GOY models.

In the following, we will characterize the level of synchronization in helical shell models with respect to their cascade properties and the presence of intermittency corrections.

### A. Shell spacing dependence

Helical shell models provide a rich framework, in which a set of interacting triads are connected in a way that mirrors the actual interactions among two sequences of

helical modes in three dimensional Navier-Stokes turbulence. Here, we consider the phase coherence of these two sequences, as a function of the inter-shell ratio $g$, in order to investigate how differences in synchronizability are related to energy transfer across scales and to compare with the single chain GOY model.

Due to their more complicated structure, helical shell model are less explored in literature. Here, we restrict our analysis to shell models constructed from a single class of interactions, as in Refs. 60, 61, focusing first on those classes that can sustain a forward energy cascade for local interactions, namely SM1, SM2 and SM3. We note that:

* The separable class SM1 behaves exactly as two uncoupled GOY models. As a result, the intermittency corrections present the same dependence with respect to the inter-shell ratio, and its synchronization properties are equivalent to those of a single chain, as characterized in detail in the first part of this work, showing a transition from the standard intermittency ($g=2$) to strong intermittency as the inter-shell spacing $g$ is decreased.

* On the other hand, SM2 is known to be non intermittent for the standard $g=2$ shell spacing, and this non-intermittent behavior is robust under variations of $g$.

* In stark contrast to all the examples that we have seen earlier, SM3 presents a level of intermittency comparable to SM1 for $g=2$ but the intermittency corrections decrease drastically in the small $g$ limit, and the model tends toward exhibiting monofractal scaling [55].

Those behaviors are presented together in Fig. 12.

The time behavior of the SM3 model for the two representative cases with inter-shell ratio $g=1.2$ and $g=2$ is presented in Fig. 13 showing the evolution of the compensated energy spectrum, that can be contrasted to Fig. 10.

In intermittent case ($g=2$), we observe bursty fluctuations around the global scaling $k_n^{-2/3}$, less defined but similar to those seen for the GOY model evolution. However no clear signature of instantons can be seen here. It is worth noting that, instantons have been observed in helical shell models when initialized from a zero-background rather than in conditions of a steady state energy cascade. In particular for SM3 class, such structures display an energy spectrum less steep than the Kolmogorov scaling, are reported in Ref. 62.

In contrast, in the non-intermittent case ($g=1.2$), while some local coherent events, that are carrying energy to smaller scales are still present, they do not seem to propagate across many shells as in the intermittent case since the compensated spectra look more noisy at smaller scales. This suggests that triadic phases lose their ability to self-organize in a way that accommodate strong flux events involving a large portion of the inertial range.
<figure><img src="image_13.png" /><figcaption>Figure 12. Plot of the Scaling exponent $\xi_p$ of the $p$-th order structure function as a function of $p$, for different Helical shell models corresponding to the classes SM1, SM2, and SM3 (indicated by different colors). The different markers styles represent the two inter-shell spacings $g = 2$ and $g = 1.1$. The scaling exponents were estimated using the ESS procedure described in Sec.III. The dashed line represents the Kolmogorov $p/3$ scaling.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 14. The plots show the decrease in phase coherence observed in the SM3 class as the inter-shell ratio is reduced, for the two representative cases $g = 1.2$ (dashed lines) and $g = 2$ (solid lines). The top panels displays local coherence showing the normalized PDFs of the global Kuramoto parameter $R$ (left), and the associated average phase $\Phi$ (right). The bottom panels display local coherence, showing the PDFs of the local Kuramoto parameter $R_n$ (left) and the associated average phase $\Phi_n$ (right) for different shell wave-numbers, as indicated by the colorbar.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 13. Temporal evolution of compensated energy, $\log(k_n^{2/3} E_n)$ across different shells (index $n$ in the $y$-axis) for the two spacing values $g = 2$ and $g = 1.2$ considered (top to bottom). The plots, corresponding to the SM3 class, illustrate the reduction of coherent structures across scales in the energy cascade as $g$ decreases.</figcaption></figure>

The SM3 class, displays a transition opposite to that of the GOY model (or SM1), showing a clear a reduction of the intermittency corrections as the inter-shell ratio $g$ is reduced. This makes SM3 suitable as a second example for investigating how the phase coherence change moving from an intermittent to a non-intermittent regime, but with an opposite dependence on $g$.

The phase coherence in the two regimes is shown in Fig. 14, where normalized PDFs of the global Kuramoto order parameter $R$, and the associated average phase $\Phi$ are presented, for the two inter-shell spacing $g = 1.2$ and $g = 2$. A clear difference in the level of synchronization is observed, as the system moves from intermittent to non-intermittent behavior, which is manifested by a marked reduction in triadic phase coherence. The monofractal case exhibits a narrow distribution at low $R$ values, in contrast to the broader distribution of $g = 2$ case, indicating that incoherent events dominate the dynamics. At the same time, the PDF of the average triadic phase $\Phi$ appears broader as the coherence is lost, even though a clustering around $\pi/2$ persist, which on average, continues to sustain a forward energy flux.

In the bottom panels, we present the PDFs of the parameters quantifying local coherence, as defined in Eq. (10), namely $R_n$ and $\Phi_n$, for shells within the inertial range. The local level of phase coherence appear to be universal, i.e. independent of the scale. Moreover, the $g = 2$ case displays a distribution of local coherence significantly different from that of the GOY model (presented in Fig. 3), also due to the differences in the definition of the local Kuramoto parameters in the two cases, even though both models show comparable levels of intermittency.

Notably, local coherence persists despite the loss of global synchronization as we can see by comparing the PDFs of $R_n$ in the two cases $g = 2$ and $g = 1.2$. When the
inter-shell spacing is decreased the position of the max-
ima shifts slightly toward lower coherence levels, and only
highly local coherent events ($R_n \simeq 0,8$) are significantly
less probable. In contrast, the statistics of both local
and global mean triadic phase $\Phi$ show that the phase dis-
tribution is more coherent for larger inter-shell spacing.
This statistical features are consistent with the observed
temporal dynamics, characterized not by structures that
travels over many shells but by events that involve neigh-
boring scales.

Regarding forward cascades, we conclude from these
results that the Kuramoto parameter, defined using the
triadic phases, provides a useful order parameter for
studying phase coherence and intermittency in a network
of interacting triads. However, it is important to stress
that it is not a unique indicator of intermittency. For in-
stance, for the same inter-shell ratio $g = 2$ the PDF of the
Kuramoto parameter differs between GOY model and the
helical model SM3, even though they present very similar
levels of intermittency. This observation suggest that the
Kuramoto parameter is also sensible to the topological
structure of the network of interacting triads. Nonethe-
less, once the topology is fixed, it proves to be a highly
sensible proxy of the modification of anomalous scaling
behavior, both when increasing and decreasing the level
of intermittency.

V. INVERSE CASCADES

Despite the success of shell models in modeling forward
cascades, attempts to extend them, in order to reproduce
the inverse energy cascade observed in 2D turbulence -
simply by changing the dimensionality of the conserved
invariants - have generally failed, mainly because of the
tendency of the inverse cascade in these models to be
overwhelmed by statistical shell-equipartition solutions
(see, e.g., Refs. 40, 63). To overcome this issue, different
generalization has been proposed, such as shell model
on hierarchical trees [64, 65], or models based on spi-
ral chains written for energy [24] or a model only for
amplitudes, or one based on a stochastic cascade mecha-
nism [66]. These results suggest that the self-consistently
evolving phase equation plays a crucial role in suppress-
ing the inverse cascade.

In this work, we will focus on the two shell model that,
to our knowledge, are able to present a genuine inverse
energy cascade: the elongated helical shell model pro-
posed by De Pietro, Biferale, and Mailybaev [61] and the
model studied by Tom and Ray [50] (see also Ref. 46
for a shell model that shows an inverse cascade driven
by resonantly interacting waves in the context of wave
turbulence). Our goal here is to analyze the statistical
properties of triadic phases and their implication in the
inverse cascade solutions produced by these models as a
complementary example to the forward cascade case.

Since inverse cascades are boringly non-intermittent
and lack coherent bursts, we do not show detailed inves-

tigations of the *synchronizability* of these models. How-
ever it is interesting to show the forms of the PDFs of
the triadic phases together with the fluxes and the re-
sulting spectra, in order to demonstrate how the triadic
phases organize statistically, to drive an inverse cascade
even when they are mostly random. We hope that, un-
derstanding how triadic phase organization works in such
reduced representations as shell models, will provide the
key to explaining how the cascade mechanism is sustained
in turbulence in general, including the inverse cascade.

1. *Inverse energy cascade in an elongated helical shell model*

A particular class of Helical interactions, represented
by triads of the type (-, +, +) or (+, -, -), i.e. the SM2,
is known to drive an inverse cascade for sufficiently elon-
gated triads. The critical elongation, at which the direc-
tion of energy transfer changes, (i.e. p/k = 0.278, assum-
ing p < k < q) was first predicted by Waleffe [59], later
confirmed with an elongated helical shell model [61], and
can apparently, also be justified using arguments based
on the conservation of pseudo-invariants [67] in a single
triad.

In this section we consider this model, representing a
unique example of an inverse cascade driven by non lo-
cal interactions as our test bed for studying the phase
coherence of elongated triads, which drive the inverse
cascade process. We recall that for non-local interac-
tions the relevant triadic phases centered at the n-th
shell can be defined as $u_{n-l}^{s_1}u_n^{s_0}u_{n+m}^{s_2} = \chi_n e^{i\varphi_n}$ and
$u_{n-l}^{-s_1}u_n^{-s_0}u_{n+m}^{-s_2} = \bar{\chi}_n e^{i\bar{\varphi}_n}$.

In Figure (15) we present the numerical results for SM2 (elongated) shell model, defined by $s_1/s_0 = -1$ and $s_2/s_0 = +1$, with an inter-shell ratio $g = \frac{1+\sqrt{5}}{2}$ and sufficiently elongated triads, specifically $l = 3$ and $m = 1$, which satisfy $g^{-l} < 0.278$. With this choice of parameters, the introduction of a small scale forcing produces an inverse cascade, as evidenced by the spectral energy flux and the presence of a clear inertial range, with an energy spectrum consistent with the Kolmogorov scaling. The probability distribution functions of non-local triadic phases, presented in the inset (top right), show a higher probability of events with triadic phases equal to $-\pi/2$, which corresponds to the choice that maximizes an inverse energy flux. The distribution appears to have a roughly sinusoidal shape and seems to be universal, in the sense that it does not seem to show any scaling with the shell number $n$.

Comparing these PDFs with the forward cascade case,
one finds a lower probability of coherent events in the in-
verse cascade. This weaker phase coherence is even more
pronounced at a global level as shown in Figure 2, where
the Kuramoto order parameter reveals a relatively inco-
herent cascade dynamics, characterized by low values of
$R$, and a clear flattening of the average phase $\Phi$ distri-
bution.
<figure><img src="image_16.png" /><figcaption>Figure 15. Plot of the energy spectrum as a function of the wave number for an elongated helical shell model (SM2), driven by small-scale forcing (forcing shells indicated in red). The inset plot on the bottom left shows the spectral flux as a function of scale, corresponding to an inverse energy flux in the inertial range. The inset plot on the top right shows the normalized PDFs of triadic phases $\varphi_n$, for the shells in the inertial range. The shell number is represented by different colors labeled by the continuous colorbar.</figcaption></figure>

## 2. Inverse energy cascades and equipartition solutions in a local shell model

For the Sabra shell model [17], a detailed phase diagram was proposed in Ref. 68, and a further investigations of the equipartition and the inverse cascade solutions have been carried out by Tom and Ray [50]. Specifically, they considered $g = 2$, by setting $a = 1$ and imposing energy conservation through the condition $c = -1-b$. This causes the second positive definite invariant (in the range $-2 < b < -1$) to have the form $H = \sum_n k_n^\alpha |u_n|^2$, with a scaling exponent $\alpha = \log_g\left(|\frac{1}{-1-b}|\right)$. Notably, the choice of enstrophy conservation, i.e. $\alpha = 2$, corresponds to $b = -5/4$ but as we already know this choice fails to drive an inverse cascade because of the issue with the equipartition in shell models. However, they showed that this shell model sustains an inverse energy cascade for values of $b$ in the range $-2 < b < b_c = -1.63$, with a spectral slope that approaches the characteristic scaling of 2D turbulence in the limit $b \to -2$. We have replicated the numerical results presented in Ref. 50, so that we can investigate the role of phase coherence in driving a local inverse cascade.

As shown in Figure 16, the model, when forced at small scales with $b = -1.8$, exhibits a clear inertial range associated with an inverse energy cascade, corresponding to a constant, on average, negative energy flux. The probability distribution of triadic phases, displayed in the top right inset, shows how the triadic phases for different shells organize around the value $-\pi/2$, which maximizes the inverse energy transfer. The distribution follows a sinusoidal shape, i.e odd with respect to $\varphi_n \to -\varphi_n$, with

<figure><img src="image_17.png" /><figcaption>Figure 16. Plot of the energy spectrum as a function of the wave number $k_n$, for the model of Tom and Ray, driven by a small-scale forcing (forcing shell indicated in red). The bottom left inset plot shows the spectral flux as a function of scale, corresponding to an inverse energy cascade. The plotted spectral slope $k_n^{-0.8}$ corresponds to the numerical estimate from previous studies cited above. The top-right inset shows the normalized PDFs of triad phases $\varphi_n$ for the shells in the inertial range. The shell number is represented by different colors labeled by the continuous colorbar.</figcaption></figure>

non universal, scale-dependent behavior, indicating that triadic phase organization is weaker, (as marked by noisier and broader curves) going to larger scales (slower time scales). In contrast (not shown), when the model exhibits energy equipartition (i.e. for $b_c < b < -1$) triadic phases appear noisy and randomly distributed.

To quantify synchronization at a global level, providing an overall view of phase coherence between forward and inverse cascades, we report in Figure 2, the global coherence order parameter $R$ and the associated average triadic phase $\Phi$. Notably, we observe that inverse cascades are generally driven by much weaker phase organization. Among the two inverse cascade cases considered, the local inverse cascade displays higher level of phase coherence compared to the one of the helical shell model driven by elongated triads. In the latter case, even though the average phase $\Phi$ appears flattened, a backward flux persists.

# VI. CONCLUSION

A detailed analysis of the dynamics of triadic phases in several shell models, interpreted as chains of interacting triads, where each triad is considered as a *nonlinear oscillator*, reveals how *synchronizability* plays a key role in the statistical properties of the energy cascade. In particular, our results show that phase-synchronized events are associated with bursts of energy transfer, which determine the intermittency of the system. Our investigation, which is in line with previous results on Burgers equation that links intermittency to phase dynamics, pro-
vides new evidence that stronger intermittency is related to enhanced triadic phase coherence, which is demonstrated first in a shell model representation of the three dimensional Navier Stokes equations.

In the classical GOY model, with fixed inter-shell spacing, $g = 2$, the system exhibits the “usual level” of intermittency. Decreasing inter-shell spacing below a critical value, while keeping the same physical quantities conserved, we systematically examine how each scaling exponent changes continuously, ultimately approaching a Burgers-like dynamical scaling in the continuum limit. This transition is accompanied by a loss of scale invariance in the statistics of triadic phases and in the Kolmogorov multipliers. Remarkably, the enhancement of intermittency correction is associated with a breaking of the statistical symmetry of both quantities: while their PDFs are scale-independent at the “usual level” of intermittency, they lose this invariance as the intermittency increases, providing an interesting challenge to hidden-symmetry theories [69, 70].

This progressive enhancement of intermittency, controlled by the inter-shell spacing, allows us to investigate the connection between synchronization and intermittency, systematically using a Kuramoto order parameter to quantify global phase coherence among the triads of the shell model. The PDFs of the order parameter $R$ reveal an increased likelihood for fully coherent events as $g$ decreases, complementary to this, the average phase $\Phi$ (i.e. the phase associated with the complex Kuramoto parameter) distribution becomes more peaked around $\pi/2$, which is incidentally the value that maximizes forward energy transfer.

In shell models, extreme energy flux events are associated with coherent pulses, propagating across the inertial range. We observe that, as $g$ decreases, those pulses becomes more localized and significantly more intense but less frequent. In order to study the correlation between the dynamics of such burst-like events and the phase coherence of the triad chain, a novel order parameter was introduced, to quantify local phase coherence. This weighted Kuramoto parameter, derived from the right-hand side of the triadic phase evolution equation, captures how a coherent pulse acts as a synchronization domino effect on consecutive triadic phases. Specifically, the passage of a coherent event is preceded by a rapid increase in the local coherence parameter that synchronizes a cluster of triads in a partially phase locked state, i.e. $R_n \simeq 1$, or $\Phi_n \simeq \pi/2$, that maximizes the energy transfer between aligned triads, suggesting that the local synchronization is accompanied by the coherent burst in energy transfer. The strong correlation between this phase locked state and large deviations in energy flux, can also be evidenced by their joint statistics.

The PDFs of the local Kuramoto parameter also show that local coherence precedes global synchronization. As intermittency increases, the coherent events generate longer-lasting phase-locked states, as demonstrated by the elevated probability of observing $R_n \simeq 1$. This suggests

that accumulation of energy at large scales due to an incoherent phase distribution at small scales, triggers the instanton, whose forced passage reorganizes the phases in such a way that the energy can continue to flow along this opened “path” for a while, which is then slowly desynchronized from small scales, and as the desynchronization propagates to large scales and becomes significant, the “path” gets closed again, which in turn causes the accumulation of energy at large scales, eventually repeating the whole cycle.

As these observations show, the local Kuramoto parameter can be a useful proxy for detecting coherent events that drive the energy cascade across range of scales. It may be possible to generalize this for detecting paths of burst events to models with spatial structure such as binary trees or more complex $k$-space structure such as nested polyhedra models, and maybe eventually to direct numerical simulations. In particular, the analysis of triadic phase dynamics in a 3D DNS could follow two different approaches: focusing on a single “path” in Fourier space, represented by a self-similar chain of connected triads in the spirit of spiral-chain models, or filtering DNS data into a set of shells and studying the triadic phase dynamics of those. However, it should be mentioned that defining a local phase coherence indicator in nonlinear systems involving many interconnected triads is a nontrivial task.

As a minimal example of the case with interconnected triads, we investigated triadic phase dynamics in helical shell models constructed from two sequences of helical variables, allowing for different classes of chiral interactions, focusing in particular on the SM3 class, which consist of two intertwined chains that are topologically non-separable. We observe that, contrary to the single chain case, here, decreasing the inter-shell spacing leads to a reduction of intermittency, with the system exhibiting self-similar monofractal scaling as it approaches the continuum limit. However, as before we find that the phase coherence correlates with intermittency, as shown by the order parameters, even though the dependence of that with $g$ is the inverse of what happens in the GOY case.

Finally, we have addressed the role of triadic phases in the inverse energy cascade, considering two specific cases, represented by a non-local helical shell model and the model investigated in Ref. 50. Our results show that inverse cascades are characterized by a PDF of triadic phases, peaked around $\phi_n = -\pi/2$, which maximizes the inverse energy transfer. However the order parameters show that the inverse cascades (at least the ones that we considered), are non-intermittent, and present significantly lower levels of phase coherence in contrast to the forward cascades. This is also in agreement with the emerging picture that phase dynamics and intermittency in cascade models are intimately connected.

In short, our findings provide new evidence that phase synchronization acts as a driver of intermittency and extreme events in shell models. We introduced novel pa-
rameters to quantify the local phase coherence in these
models, which can be extended to similar cascade mod-
els using a network perspective, and eventually to direct
numerical simulations.

ACKNOWLEDGMENT

This work has benefited from a grant managed by
the Agence Nationale de la Recherche (ANR), as part
of the program ‘Investissements d’Avenir’ under the ref-
erence (ANR-18-EURE-0014) and has been carried out
within the framework of the EUROfusion Consortium,
funded by the European Union via the Euratom Re-
search and Training Programme (Grant Agreement No
101052200 — EUROfusion) and within the framework of

the French Research Federation for Fusion Studies. The
authors would like to thank the Isaac Newton Institute
for Mathematical Sciences, Cambridge, for support and
hospitality during the programme “Anti-diffusive dynam-
ics: from sub-cellular to astrophysical scales”, where part
of the work on this paper was undertaken. This work
was supported by EPSRC grant EP/R014604/1. In par-
ticular, the authors would like to thank Dr. Santiago J.
Benavides for stimulating discussion on phase dynamics
during the program.

DATA AVAILABILITY

The data that support the findings of this article are
openly available at [71].

[1] P. Bak, K. Christensen, L. Danon, and T. Scanlon, Phys. Rev. Lett. 88, 178501 (2002).

[2] P. Gopikrishnan, V. Plerou, L. A. Nunes Amaral, M. Meyer, and H. E. Stanley, Phys. Rev. E 60, 5305 (1999).

[3] S. Lovejoy and D. Schertzer, *The weather and climate: emergent laws and multifractal cascades* (Cambridge University Press, 2013).

[4] U. Frisch, *Turbulence: The Legacy of A. N. Kolmogorov* (Cambridge University Press, Cambridge, 1995).

[5] T. Bohr, M. H. Jensen, G. Paladin, and A. Vulpiani, *Dynamical Systems Approach to Turbulence*, Cambridge Nonlinear Science Series (Cambridge University Press, 1998).

[6] L. Biferale, Ann. Rev. Fluid Mech. 35, 441 (2003).

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

[8] Y. Hattori, R. Rubinstein, and A. Ishizawa, Phys. Rev. E 70, 046311 (2004).

[9] M. H. Jensen, G. Paladin, and A. Vulpiani, Phys. Rev. A 45, 7214 (1992).

[10] A. Kumar and M. K. Verma, Phys. Rev. E 91, 043014 (2015).

[11] F. Plunian, R. Stepanov, and P. Frick, Phys. Rep. - Rev. Sect. Phys. Lett. 523, 1 (2013).

[12] A. Verdini and R. Grappin, Phys. Rev. Lett. 109, 025004 (2012).

[13] V. Berionni, P. Morel, and Ö. D. Gürcan, Physics of Plasmas 24, 122310 (2017).

[14] K. Ghantous and Ö. D. Gürcan, Phys. Rev. E 92, 033107 (2015).

[15] L. Kadanoff, D. Lohse, J. Wang, and R. Benzi, Phys. Fluids 7, 617 (1995).

[16] K. Ohkitani and M. Yamada, Progress of Theoretical Physics 81, 329 (1989).

[17] V. S. L'vov, E. Podivilov, A. Pomyalov, I. Procaccia, and D. Vandembroucq, Phys. Rev. E 58, 1811 (1998).

[18] V. S. L'vov, A. Pomyalov, and I. Procaccia, Phys. Rev. E 63, 056118 (2001).

[19] X. M. de Wit, G. Ortali, A. Corbetta, A. A. Mailybaev, L. Biferale, and F. Toschi, Phys. Rev. E 109, 055106 (2024).

[20] A. A. Mailybaev, Phys. Rev. E 86, 025301(R) (2012).

[21] A. A. Mailybaev, Phys. Rev. E 87, 053011 (2013).

[22] Ö. D. Gürcan, Phys. Rev. E 95, 063102 (2017).

[23] Ö. D. Gürcan, P. Morel, S. Kobayashi, R. Singh, S. Xu, and P. H. Diamond, Phys. Rev. E 94, 033106 (2016).

[24] Ö. D. Gürcan, S. Xu, and P. Morel, Phys. Rev. E 100, 043113 (2019).

[25] C. S. Campolina and A. A. Mailybaev, Phys. Rev. Lett. 121, 064501 (2018).

[26] C. S. Campolina and A. A. Mailybaev, Nonlinearity 34, 4684 (2021).

[27] G. Costa, A. Barral, and B. Dubrulle, Phys. Rev. E 107 (2023), 10.1103/PhysRevE.107.065106.

[28] A. S. Lanotte, R. Benzi, S. K. Malapaka, F. Toschi, and L. Biferale, Phys. Rev. Lett. 115, 264502 (2015).

[29] M. Buzzicotti, L. Biferale, U. Frisch, and S. S. Ray, Phys. Rev. E 93 (2016), 10.1103/PhysRevE.93.033109.

[30] M. Buzzicotti, B. P. Murray, L. Biferale, and M. D. Bustamante, Eur. Phys. J. E 39 (2016), 10.1140/epje/i2016-16034-5.

[31] B. P. Murray and M. D. Bustamante, Journal of Fluid Mechanics 850, 624 (2018).

[32] B. Protas, D. Kang, and M. D. Bustamante, Phys. Rev. E 109 (2024), 10.1103/PhysRevE.109.055104.

[33] Ö. D. Gürcan, Y. Li, and P. Morel, Mathematics 8, 530 (2020).

[34] Ö. D. Gürcan, Physica D: Nonlinear Phenomena 426, 132983 (2021).

[35] Ö. D. Gürcan, Reviews of Modern Plasma Physics 7, 20 (2023).

[36] L. Manfredini and O. D. Gürcan, Phys. Rev. E 111, 025103 (2025).

[37] A. Alexakis and L. Biferale, Physics Reports 767-769, 1 (2018), cascades and transitions in turbulent flows.

[38] S. J. Benavides and M. D. Bustamante, “Phase dynamics and their role determining energy flux in hydrodynamic shell models,” (2025), arXiv:2507.03397 [physics.flu-dyn].

[39] L. Biferale, A. Lambert, R. Lima, and G. Paladin, Physica D: Nonlinear Phenomena 80, 105 (1995).

[40] P. Ditlevsen and I. Mogensen, Phys. Rev. E 53, 4785 (1996).
[41] N. M. Rathmann and P. D. Ditlevsen, Phys. Rev. E **94** (2016), 10.1103/PhysRevE.94.033115.

[42] D. Pisarenko, L. Biferale, D. Courvoisier, U. Frisch, and M. Vergassola, Physics of Fluids A **5**, 2533 (1993).

[43] L. Biferale, G. Boffetta, A. Celani, and F. Toschi, Physica D: Nonlinear Phenomena **127**, 187 (1999).

[44] R. Benzi, L. Biferale, and G. Parisi, Physica D: Nonlinear Phenomena **65**, 163 (1993).

[45] L. Biferale, A. A. Mailybaev, and G. Parisi, PHYSICAL REVIEW E **95** (2017), 10.1103/PhysRevE.95.043108.

[46] N. Vladimirova, M. Shavit, and G. Falkovich, Phys. Rev. X **11**, 021063 (2021).

[47] S. H. Strogatz, Nature **410**, 268 (2001).

[48] S. H. Strogatz, R. E. Mirollo, and P. C. Matthews, Phys. Rev. Lett. **68**, 2730 (1992).

[49] Y. Kuramoto, *Chemical Oscillations, Waves, and Turbulence* (Springer–Verlag, New York, 1984).

[50] R. Tom and S. S. Ray, EPL **120** (2017), 10.1209/0295-5075/120/34002.

[51] J.-A. Arguedas-Leiva, E. Carroll, L. Biferale, M. Wilczek, and M. D. Bustamante, Phys. Rev. Res. **4**, L032035 (2022).

[52] R. Benzi, S. Ciliberto, R. Tripiccione, C. Baudet, F. Masaioli, and S. Succi, Phys. Rev. E **48**, R29 (1993).

[53] D. Banerjee, S. S. Ray, G. Sahoo, and R. Pandit, Phys. Rev. Lett. **111**, 174501 (2013).

[54] S. Chakraborty, U. Frisch, and S. S. Ray, J. Fluid Mech. **649**, 275 (2010).

[55] R. Benzi, L. Biferale, and E. Trovatore, Phys. Rev. Lett. **77**, 3114 (1996).

[56] K. Andersen, T. Bohr, M. Jensen, J. Nielsen, and P. Olesen, Physica D **138**, 44 (2000).

[57] J. Bec and K. Khanin, Physics Reports **447**, 1 (2007).

[58] R. Benzi, I. Castaldi, F. Toschi, and J. Trampert, Phil. Trans. R. Soc. A **380**, 20210074 (2022), https://royalsocietypublishing.org/doi/pdf/10.1098/rsta.2021.0074.

[59] F. Waleffe, Physics of Fluids A: Fluid Dynamics **4**, 350 (1992).

[60] R. Benzi, L. Biferale, R. M. Kerr, and E. Trovatore, Phys. Rev. E **53**, 3541 (1996).

[61] M. De Pietro, L. Biferale, and A. A. Mailybaev, Phys. Rev. E **92**, 043021 (2015).

[62] M. De Pietro, A. A. Mailybaev, and L. Biferale, Phys. Rev. Fluids **2**, 034606 (2017).

[63] E. Aurell, G. Boffetta, A. Crisanti, P. Frick, G. Paladin, and A. Vulpiani, Phys. Rev. E **50**, 4705 (1994).

[64] E. Aurell, P. Frick, and V. Shaidurov, Physica D: Nonlinear Phenomena **72**, 95 (1994).

[65] E. Aurell, E. Dormy, and P. Frick, Phys. Rev. E **56**, 1692 (1997).

[66] P. D. Ditlevsen, Physics of Fluids **24**, 105109 (2012).

[67] N. M. Rathmann and P. D. Ditlevsen, Phys. Rev. Fluids **2**, 054607 (2017).

[68] T. Gilbert, V. L'vov, A. Pomyalov, and I. Procaccia, Phys. Rev. Lett. **89** (2002), 10.1103/PhysRevLett.89.074501.

[69] A. A. Mailybaev, Phys. Rev. Fluids **7**, 034604 (2022).

[70] A. A. Mailybaev, Nonlinearity **35**, 3630 (2022).

[71] L. Manfredini and Ö. Gürcan, “Repository for nonlinear phase synchronization and the role of spacing in shell models,” (2025), zenodo dataset.

[72] C. Rackauckas, “Sciml/differentialequations.jl: v7.14.0,” https://doi.org/10.5281/zenodo.13629099 (2024).

[73] C. A. Kennedy and M. H. Carpenter, Applied Numerical Mathematics **44**, 139 (2003).

[74] J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah, SIAM Review **59**, 65 (2017), https://doi.org/10.1137/141000671.

## Appendix A: Numerical details

The dynamics of the different shell models presented have been integrated using a 4th order Implicit-Explicit (IMEX) solver [72, 73] implemented in Julia [74]. This approach allows for efficient separation of the stiff and non-stiff components in the coupled ordinary differential equations system that govern the shell model evolution.

For all the shell models considered, given the inter-shell spacing $g$, we fix a number of shells $N$ in order to cover approximately seven decades. The dissipation term is set as $d_n = \nu k_n^2$ where the viscosity $\nu$ is chosen small enough to dissipate energy only in the last few shells. For the inverse cascade case, an additional hypo-viscosity term $(\mu k_n^{-2})$ has been included to prevent energy accumulation at large scales.

The statistical analyses are performed in a steady state regime, over a time window corresponding to approximately 2,000 eddy turnover times ($\tau_E^{-1} \sim k_n u_n$) of the largest scale. The sampling time is set to $\delta t = 5 \times 10^{-5}$, ensuring a sufficiently high sampling rate to resolve the temporal dynamics within the inertial range. The initial condition typically consists of a small amplitude ($10^{-8}$) white noise. The Tables I, II and III summarize all the parameters used for the various runs that corresponding to the numerical experiments analyzed in the paper.

A special note concerns the choice of the forcing term. The common approach in the literature is to apply either a constant or random forcing. In the specific case discussed in Sec. III to investigate the role of phase self-organization, particularly in relation to phase-synchronized events approaching the continuum limit, we implemented a forcing term that does not inject any phase into the system. This is achieved by ensuring that the term $\text{Im}(\frac{f_n}{\rho_n}e^{-i\theta_n})$ on the right-hand side of Eq. (4) is identically zero, thus we fixed $f_n = A e^{i\theta_n}$, with $A = 0.5$, in all the runs corresponding to Table I. Another possible choice, $f_n = \epsilon/u_n^*$, which does not inject any phase and fix a constant energy injection $\epsilon$, as discussed in Ref. 58 in order to study the statistics of inter-event times. We note that for the result of present analysis there are no significant differences among these two choices.
GOY model Runs

<table><thead><tr><th>g</th><th>1.1</th><th>1.15</th><th>1.2</th><th>1.3</th><th>1.4</th><th>1.5</th><th>1.6</th><th>1.7</th><th>1.8</th><th>1.9</th><th>2.0</th></tr></thead><tbody><tr><th>N</th><td>164</td><td>113</td><td>88</td><td>62</td><td>50</td><td>42</td><td>38</td><td>36</td><td>33</td><td>31</td><td>29</td></tr><tr><th>ν</th><td>5 × 10<sup>-9</sup></td><td>4 × 10<sup>-9</sup></td><td>8 × 10<sup>-9</sup></td><td>14 × 10<sup>-9</sup></td><td>9 × 10<sup>-9</sup></td><td>8 × 10<sup>-9</sup></td><td>5 × 10<sup>-9</sup></td><td>9 × 10<sup>-10</sup></td><td>7 × 10<sup>-10</sup></td><td>5.5 × 10<sup>-10</sup></td><td>5.5 × 10<sup>-10</sup></td></tr></tbody></table>

Table I. Numerical parameters corresponding to the different runs of the GOY model presented is Sec. III. In all the different runs the shell wave-numbers correspond to $k_n = k_0 g^n$ with $k_0 = 1/g^3$, and the forcing act the third shell i.e. $f_3 = Ae^{i\theta_3}$ with $A = 0.5$.

Helical shell models Runs

<table><thead><tr><th>class</th><th>SM1</th><th>SM1</th><th>SM2</th><th>SM2</th><th>SM3</th><th>SM3</th><th>SM3</th></tr></thead><tbody><tr><th>g</th><td>1.1</td><td>2</td><td>1.1</td><td>2</td><td>1.1</td><td>1.2</td><td>2</td></tr><tr><th>N</th><td>164</td><td>30</td><td>140</td><td>29</td><td>140</td><td>80</td><td>30</td></tr><tr><th>ν</th><td>5 × 10<sup>-10</sup></td><td>1.5 × 10<sup>-10</sup></td><td>8 × 10<sup>-10</sup></td><td>3.5 × 10<sup>-10</sup></td><td>8 × 10<sup>-10</sup></td><td>1 × 10<sup>-9</sup></td><td>1.5 × 10<sup>-10</sup></td></tr></tbody></table>

Table II. Numerical parameters corresponding to the different runs of the local Helical shell models. $f_3 = (1 + 1j) \times 7 \times 10^{-1}$

Inverse cascade shell models Runs

<table><thead><tr><th>model</th><th>Helical-nonlocal</th><th>Sabra-local</th></tr></thead><tbody><tr><th>g</th><td>φ</td><td>2</td></tr><tr><th>N</th><td>36</td><td>28</td></tr><tr><th>ν</th><td>4 × 10<sup>-12</sup></td><td>5 × 10<sup>-12</sup></td></tr><tr><th>μ</th><td>1</td><td>1 × 10<sup>-5</sup></td></tr><tr><th>f<sub>n</sub></th><td>f<sup>+</sup><sub>32,33</sub> = 0.8, f<sup>-</sup><sub>32,33</sub> = 0.4, f<sub>23</sub> = (5 + 5j) × 10<sup>-4</sup></td><td></td></tr></tbody></table>

Table III. Numerical parameters corresponding to the different runs of the non local Helical shell model and the local sabra model discussed in Sec.V 2