Percolation transition in the kinematics of nonlinear resonance broadening in Charney-Hasegawa-Mima model of Rossby wave turbulence

We study the kinematics of nonlinear resonance broadening of interacting Rossby waves as modelled by the Charney-Hasegawa-Mima equation on a biperiodic domain. We focus on the set of wave modes which can interact quasi-resonantly at a particular level of resonance broadening and aim to characterise how the structure of this set changes as the level of resonance broadening is varied. The commonly held view that resonance broadening can be thought of as a thickening of the resonant manifold is misleading. We show that in fact the set of modes corresponding to a single quasi-resonant triad has a nontrivial structure and that its area in fact diverges for a finite degree of broadening. We also study the connectivity of the network of modes which is generated when quasi-resonant triads share common modes. This network has been argued to form the backbone for energy transfer in Rossby wave turbulence. We show that this network undergoes a percolation transition when the level of resonance broadening exceeds a critical value. Below this critical value, the largest connected component of the quasi-resonant network contains a negligible fraction of the total number of modes in the system whereas above this critical value a finite fraction of the total number of modes in the system are contained in the largest connected component. We argue that this percolation transition should correspond to the transition to turbulence in the system.


I. INTRODUCTION
The phenomenon of dispersive wave propagation is fundamental to our understanding of a wide variety of spatially extended physical systems. In such systems, the frequency, ω k , of each wave mode is a nonlinear function of its wave-vector, k [1]. Examples include gravitycapillary waves on fluid interfaces [2], flexural waves in thin elastic plates [3], drift waves in strongly magnetized plasmas [4] and Rossby waves in planetary oceans and atmospheres [5]. In this article, we will be interested in Rossby waves as modeled by the Charney-Hasegawa-Mima (CHM) equation [4] on the β-plane: This is the simplest two-dimensional model of the large scale dynamics of a shallow layer of fluid on the surface of a strongly rotating spehere. The surface of the sphere is approximated locally by a plane, x ∈ R 2 , with x varying in the longitudinal (meridional) direction and the y varying in the latitudinal (zonal) direction. The field ψ(x, t) is the geopotential height, β is the Coriolis parameter measuring the variation of the Coriolis force with latitude, F is the inverse of the square of the deformation radius and J [f, g] denotes the Jacobian of two functions, f and g. This equation admits harmonic solutions, ψ(x, t) = Re [A k exp(ik · x − iω k t)] with k ∈ R 2 . These solutions, known as Rossby waves, have the anisotropic dispersion relation Since Eq.(1) is nonlinear, modes with different wavevectors couple together and exchange energy. If the non-linearity is weak, one finds that this energy exchange is generally quite slow and occurs most efficiently between groups of modes which are in resonance. For the CHM equation, such resonances involve three modes since the nonlinearity is quadratic. Four modes would be involved in the case of systems with cubic nonlinearity. Three wave-vectors (k 1 , k 2 , k 3 ) satisfying the resonance conditions, are referred to as a resonant triad. If one projects the spectral representation of the wave equation onto a resonant triad, one obtains a set of ordinary differential equations for the coupled time evolution of the amplitudes of the constituent modes. Such systems of equations appeared as basic models of nonlinear mode coupling in a variety of physical systems including plasma physics [6], nonlinear optics [7] and oceanic internal waves [8]. An advantage of such models is that the equations of motion for a resonant triad are simple enough that explicit formulae can be obtained for both amplitudes and phases of the resonant modes [9][10][11][12][13]. A disadvantage, however, is that such triads are generally not closed. Even if energy is initially mostly restricted to a single triad, other resonant modes can be generated which are not in the original triad. This process can repeat in a cascade-like fashion and result in the excitation of a large number of modes. If a large number of degrees of freedom are excited, a statistical description of energy transfer between modes is preferable. Such a description is provided by the theory of wave turbulence [2,14]. This theory provides a kinetic description of energy transfer in ensembles of weakly interacting disper-sive waves in which conserved quantities are redistributed along the resonant manifolds. See [15,16] for a review.
For an infinite system, in which wave modes are indexed by a continuous wave-vector, the theory of weakly nonlinear wave turbulence becomes asymptotically exact in the weakly nonlinear limit. For finite sized systems, in which the wave modes are indexed by a discrete wavevector, some subtleties arise. The simplest case, which is particularly relevant to numerical studies of wave turbulence, is a bi-periodic box. In this case, k, is restricted to a periodic lattice with a minimum spacing, ∆k, between modes. Modulo this spacing, the components of k must be integer valued. This is an issue because if the components of k are integers, then the resonance conditions, Eq.(3), become a problem of Diophantine analysis. Such problems typically have far fewer solutions than their real-valued counterparts and it is generally quite difficult to find them. A complete enumeration of all solutions for the case of Eq.(2) with F = 0 was recently provided in [17]. This sparseness of solutions means that, in discrete systems, resonant triads can exist in isolation or in finite groups of triads known as resonant clusters. Two triads belong to the same cluster if they share at least one mode. The dynamics of small clusters consisting of two triads has been studied in considerable detail in [18]. Small clusters have attracted some interest in the context of atmospheric dynamics as a possible explanation of the unusual periods of certain observed atmospheric oscillations [19]. Depending on the dispersion relation, there may or may not exist large clusters capable of distributing energy over a large range of scales in a discrete system. In the case of Rossby waves on a sphere with infinite deformation radius, such a large cluster does exist [20]. On the other hand, for capillary waves there are no exactly resonant triads at all [21]. For the dispersion relation Eq.(2), numerical explorations indicate that for general values of F large exactly resonant clusters are rare. Thus, in discrete systems, it is often necessary to rely on approximate resonances to account for energy transfer.
Approximate resonance is possible due to the phenomenon known as nonlinear resonance broadening. This is an effect whereby the frequency of a wave acquires a correction to its linear value which depends on the amplitude (see Chaps. 14 & 15 of [1]). Triads which are not exactly in resonance can then interact at finite amplitude if the frequency mismatch is less than this correction. Such triads are known as quasi-resonant triads and satisfy the broadened resonance conditions where δ is a characteristic value for the resonance broadening taken to be positive. Although Eq.(4) provides only a kinematic description of resonance broadening, the analogous dynamical effect can be very strikingly visualised in linear stability analyses of weakly nonlinear waves [22,23] where it is found that the set of unstable perturbations lie in a neighbourhood around the set of exactly resonant perturbations. This set of quasi-resonant modes is often pictured as a "thickened" or broadened version of the exactly resonant manifold. For weakly nonlinear systems, this broadening is expected to be small since amplitudes are small. It may nevertheless be large enough to overcome frequency mismatches which arise when wave-vectors are restricted to a discrete grid and prevent discreteness from impeding the cascade of energy. A striking example of this effect is observed for capillary wave turbulence in a biperiodic box. For this system, since there are no exact resonances, as the nonlinearity is decreased the resonance broadening eventually becomes smaller than the frequency mismatches due to the grid. Direct numerical simulations illustrate that the cascade of energy to small scales stops entirely when the level of nonlinearity gets sufficiently small leading to the phenomenon of "frozen turbulence" [24][25][26].
The interplay between exactly resonant and quasiresonant clusters means that wave turbulence in discrete systems is nowadays believed to exhibit several regimes. If the typical resonance broadening, δ, is small enough that effectively only exactly resonant clusters can interact, the dynamics are referred to as discrete wave turbulence [20,27,28]. If δ is larger than the typical spacing between modes then effectively all triads can interact at least quasi-resonantly and the classical statistical theory is expected to be valid. In between is a regime consisting of a mixture of exactly resonant and quasi-resonant clusters which has been termed mesoscopic wave turbulence [28,29]. In this intermediate regime, it has been suggested [30] that forced systems could exhibit some aspects of self-organised criticality. This suggestion is motivated by the idea that the forcing will cause the characteristic value of δ to increase until it is large enough for a large quasi-resonant cluster to form which will then facilitate an "avalanche" of energy transfer to the dissipation scale thereby reducing wave amplitudes and the corresponding value of δ.
In this paper we develop a kinematic concept of criticality in quasi-resonant interactions in the Eq.(1). Specifically, we address the question of how a large quasiresonant cluster emerges in the CHM equation as δ is increased. Inspired by the theory of percolation on random networks [31], we take "large cluster" to mean a cluster that consists of a finite fraction of all modes in the system. We begin by analytically characterising the shape of the quasi-resonant set defined by Eq.(4) as a function of δ for a single triad in Sec. II. By expressing the boundary of the quasi-resonant set in terms of the intersection of a pair of quadratic forms, we find some surprises. In particular, we find that the area of the set diverges at a finite value of δ illustrating that the common perception of the quasi-resonant set as a "thickened" version of the exact resonant manifold is potentially quite misleading. In Sec. III we numerically construct the set of quasi-resonant clusters as a function of δ for various system sizes. We show that a percolation transition occurs at a critical value, δ * , of the resonance broadening as δ is increased. At this critical value, the size of the largest cluster rapidly goes from containing a negligible fraction of the modes in the system to containing a finite fraction of them. The value of δ * decreases as the inverse cube of the system size, a fact which we trace to quasi-resonant interactions between small scale meridional modes and large scale zonal modes. We finish with a short summary and discussion about what conclusions can be drawn about Rossby wave turbulence from our results.

II. CHARACTERISATION OF THE QUASI-RESONANT SET FOR A SINGLE TRIAD
In what follows, we shall take k 3 to be fixed with k 2 = k 3 − k 1 . The δ-detuned quasi-resonant set of k 3 is the set of modes, k 1 , which satisfy the inequality This section is devoted to determining the structure of this set as a function of the detuning, δ. The boundaries of this set are given by the pair of curves We begin by finding these curves. Clearly it suffices to solve Eq.(6) since the second boundary can be obtained from this by setting δ → −δ. To fix notation, let us write For the CHM dispersion relation, Eq.(2), the boundary of the quasi-resonant set, Eq.(6), then corresponds to the curve in the (x, y) plane implicitly defined by Subsequent formulae will be more compact if we shift the origin to the centre of symmetry of the curve, by defining x = r − p/2, y = s − q/2. Also, we rescale δ by setting β = 1 from here on. The curve we wish to study is therefore Gathering these terms together with a common denominator we obtain the quartic curve c(x, y) = 0, where c(x, y) = a 1 + a 2 (x 2 + y 2 ) 2 + a 3 x 2 + a 4 y 2 − a 5 xy, (9)  11) and (13) respectively. and the coefficients are given by If we introduce new variables u = x 2 , v = y 2 and w = xy, then it becomes clear that the boundary curve, Eq. (8), corresponds to the intersection of 2 quadratic surfaces Notice that the surface w 2 = uv is a cone, which is singular at the origin x = y = 0. We have performed a full analysis of the properties of the intersection curves for general values of the parameters p, q, F and δ. This is a technical exercise which is not very illuminating. In the interests of clarity, we will restrict ourselves here to illustrating the essential qualitative features of these curves accompanied by some explicit examples.
A. The curve typically exists only for a finite range of δ. Examining Eq. (10), it is clear that since u and v are both positive, it is generally not possible to balance the quadratic terms if u 2 + v 2 → ∞. Hence the curve, if it exists is bounded. The single exception to this occurs when the coefficient of the quadratic terms vanishes. The curve may therefore diverge at the single critical value of detuning given by: Note that this corresponds to δ = ω(k 3 ). At δ = δ 1 , some algebra shows that the curve is given by the hyperbola It follows that the area of the quasi-resonant set diverges as δ tends to δ 1 from below. This result illustrates that the common conception that the quasi-resonant set looks like a "thickened" version of the exact resonant manifold is a misconception. The special case p → 0 corresponding to the case of k 3 becoming zonal is discussed separately below. C. The curve may self-intersect only at a single critical value δ.
From Eq.(10), self-intersection is only possible if the curve passes through (0, 0) in the (u, v) plane. This can only happen if the coefficient a 1 = 0. Thus we identify a second critical value of δ where self-intersection may occur: Note that a 1 = 0 does not necessarily mean the curve selfintersects. It is possible that at δ = δ 2 , the curve reduces to a single point. To probe what happens at δ = δ 2 , we consider the surface z = c(x, y) defined by Eq.(9) when δ = δ 2 in the neighbourhood of the origin. Calculation of the partial derivatives indicate that this surface has a critical point at the origin. After some tedious algebra, we find that the determinant of the matrix of second derivatives is If ∆ > 0 then the critical point at (0, 0) is a maximum or a minimum. The curve c(x, y) = 0 is then an isolated point. On the other hand, if ∆ < 0 the critical point is a saddle and the curve c(x, y) = 0 has a self-intersection at (0, 0). The condition for self-intersection is therefore We note that for F = 0 we always have a self-intersection. These qualitative features of the boundary of the quasiresonant set are illustrated graphically in Figs. 1 and 2 which show the shape of the curve for different values of δ. Generic parameters, specified in the captions, have been chosen with no particular symmetries. Hence the curves shown in these figures are representative of what is found from the complete analysis of Eq.(8). Fig. 1 shows an example in which no self-intersection occurs at δ = δ 2 , while Fig. 2 shows an example in which a selfintersection occurs. The hyperbolic curves identified at δ = δ 1 are clearly visible in both cases.
A couple of special cases are worth noting: (i) p = 1, q = 0 and F = 0 This corresponds to a meridional mode. In this case, a 5 = 0 so the intersection of quadratic forms given by Eq.(10) lies entirely in the w = 0 plane and reduces to  We can immediately identify the critical points. The curve self-intersects at δ = 3 and diverges at δ = −1. The point δ = 1 is also noteworthy, being the complementary boundary to the divergent case. For this value of δ the curve is a perfect circle.
(ii) p = 0, q = 1 and F = 0 This corresponds to a zonal mode. In this case the curve simplifies to The only special value of δ for this case is δ = 0. We then recover the exact resonant manifold of a zonal mode, x y = 0. This consists of the two coordinate axes. It is now less surprising that the boundary of the quasi-resonant set can diverge for finite δ once we appreciate that the exact resonant manifold of a zonal mode is unbounded. The divergence of the boundary of the quasi-resonant set of the non-zonal modes in some sense reflects the presence of this structure in the dispersion relation.
Remark. Notice that the exactly resonant manifold in case (ii) consists of non-generic triads, i.e., triads where one or more interaction coefficients are identically zero. A point on the x-axis corresponds to a so-called catalytic interaction (|k 1 | = |k 2 |), for which the zonal mode k 3 does not change its energy but influences the energy exchange between the other two modes in the triad. A point on the y-axis corresponds to a 'spurious' triad, formed by purely zonal modes, which do not interact at all: the interaction coefficients are all identically zero because the three modes are collinear. In our computations of the network of quasi-resonant modes (Section III), non-generic triads are discarded from the start. Having given a fairly complete qualitative description of the behaviour of the curves which define the boundaries of the resonant sets as the detuning is varied, it now remains to assemble these boundary curves for postive and negative values of δ to determine the interior and exterior of the quasi-resonant set. This is again best accomplished by illustration. Figs. 3 and 4 illustrate the quasi-resonant set for the two special cases discussed above for a range of increasing values of δ. The shaded areas in these figures correspond to the quasi-resonant sets. The exact resonant curve is also shown for reference. Fig.3 illustrates one of the key points of this article: the common picture of the quasi-resonant set as a thickened version of the exact resonant manifold is appropriate only for small values of the broadening (eg Fig.3(b)). One might counter this with the observation that it is only in the weakly nonlinear regime that it makes sense to be discussing quasi-resonant interactions in the first place and in this regime, the broadening is necessarily small. Eq.(11) tells us, however, that no matter how small the broadening, there are always modes close to the zonal axis, whose quasi-resonant set diverges.

III. STRUCTURE OF THE NETWORK OF QUASI-RESONANT MODES
In a turbulent system many modes are excited. The fact that a mode can be a member of more than one quasi-resonant set leads to overlap between sets. This allows modes to join together to form a network of quasiresonant clusters analogous to the exactly resonant clusters which have been extensively studied in the literature to date. In this section we study the structure of this quasi-resonant network as the typical amount of broadening, δ, in the system is varied. We consider a finite system of size 2π × 2π and truncate the wavenumber space such that there are M modes in the k x and k y directions, M 2 modes in all. The spacing between modes is 2π/M . The quasi-resonant clusters were obtained by an exhaustive numerical search speeded up by incorporating symmetries of the triads. As explained in the remark below equation (16), we have eliminated from our list of triads the so-called non-generic triads, i.e., the triads for which one or more interaction coefficients are zero. This includes triads formed out of collinear modes and triads where any two wave-vectors have the same modulus.
We first measure the total number of clusters normalised by the maximum number of possible isolated clusters which is approximately M 4 /3. This is shown as a function of δ for several different system sizes, M , in Fig. 5. For each system size, as δ increases, the total number of clusters first increases, then reaches a maximum at a particular value of δ, which we shall denote by δ * , and then decreases. The value, δ * , at which the maximum occurs decreases as the system size, M , is increased. In order to understand these results one must realise that for the CHM dispersion relation, the number of exactly resonant clusters which live on the grid of discrete wavevectors is rather small with most discrete triads necessarily exhibiting a small amount of detuning enforced by the geometry. The initial growth in the number of clusters is explained by the addition of isolated triads to the list of quasi-resonant clusters as the broadening grows large enough to incorporate the geometrical detuning generated by the grid. While one would expect this rate of increase to slow down as clusters start joining to form larger clusters, the attainment of a maximum and subsequent decrease in the total number of clusters becomes much easier to understand when we measure the size of largest cluster, ρ, which is simply the number of modes in the largest cluster normalised by the system size, M 2 . This is shown in Fig. 6. We see that as δ approaches δ * the largest cluster quickly makes a transition from including a negligible fraction of the total number of modes in the system to including almost all of them. Thus the network of quasi-resonant modes undergoes a percolation transition at δ = δ * . While the results shown in Figs. 5 and 6 were obtained for β = F = 1, we checked that the results are not very sensitive to this choice. Although this is an entirely kinematic study, we expect that the emergence of a percolating quasi-resonant cluster is the explanation for the dynamical transition from discrete wave turbulence to Kolmogorov turbulence in the CHM equation as the nonlinearity of the system is increased.
The decrease of the percolation threshold, δ * , as the system size, M , increases is clearly due to the fact that as M increases, the spacing between discrete modes decreases. The amount of broadening required to overcome the geometrical detuning is therefore smaller and hence clusters can form more easily. We would like to quantify this decrease. We have already learned that the area of the quasiresonant set of any particular mode, k, diverges at a finite value of δ = ω(k) as indicated in Eq. (11). It seems plausible that as soon as the broadening is sufficiently large to allow any mode in the system to have a divergent quasi-resonant set, then the largest cluster must be of the size of the system. Thus one might esti- By this argument, Eq.(2) tells us that for large systems we should see the scaling, δ * ∼ M −2 . A numerical measurement of how δ * scales with system size for a range of values of F , is plotted in figure 7. We see that the value of F becomes redundant for large M as already mentioned. Each of the curves converge to the straight line indicated by the scaling with exponent σ 1 = 3.00 with a 95% confidence intervals if (2.97, 3.03), which was obtained by bias-corrected bootstrapping. The connected component therefore forms much more easily than the naive argument above would suggest. The origin of this M −3 scaling can be traced to the fact that zonal modes require very little detuning in order to interact with high k almost-meridional modes. Take F = 0 for simplicity, although the following argument works for any F > 0. Let us consider the largestscale zonal mode in the system, k = (0, 1). We recall that the exactly resonant manifold xy = 0 of this zonal mode, figure 4(a), gives rise to non-generic triads which do not interact efficiently and therefore must be discarded. So we ask what is the minimum value of the detuning so that the quasiresonant set of the mode k = (0, 1) contains some new modes. From Fig. 4 obtain For large M , this has solution δ * ∼ M −3 , in agreement with Fig. 7 and the inset of Fig. 6. This suggests that the percolation transition is driven by interactions between large scale zonal modes and small scale meridional modes. This is consistent with the known scale-nonlocality of wave turbulence in the CHM equation (see [32] and the references therein) and provides further evidence, albeit circumstantial, that the percolation transition is associated with the onset of turbulence in the CHM model. In order to further illustrate this point, we can ask which modes have the greatest or least tendency to become part of a quasi-resonant cluster. This is done in figure (8), where we have coloured each mode according to the minimal amount of resonance broadening required for this mode to join a cluster of any size. Blue modes become active very easily whereas red modes are resilient to becoming part of any cluster. Rather than appearing homogeneous and random, we see the appearance of definite structure. We see a circular region containing modes with a strong propensity to join clusters including a narrow large-scale region of zonal modes accumulating near the k x = 0 axis with very low interaction thresholds as expected from the discussion above. We also remark upon the group of large-scale meridional scales reluctant to form any quasi-resonant connections. The relative reluctance of large-scale meridional modes to exchange energy has already been remarked up in the literature and suggested as an explanation of the inherent anisotropy of Rossby waves turbulence. For example in [33], a wave-turbulence boundary was computed by comparing the CHM dispersion relation to the inverse of the eddy-turnover time. It is then argued that inside this region Rossby waves dominate, but with a frequency incommensurate with that of the surrounding turbulence so that energy cannot penetrate into this region. The boundary thus obtained seems similar in position and shape to the dark region in Fig. 8. This anisotropic energy distribution at large scales has been documented in numerical simulations of Rossby wave turbulence (see [34] and the references therein). It is somewhat surprising to see it emerging again here from purely kinematic considerations. Finally we might ask whether the density of the giant cluster illustrated in Fig. 6 shows the generic profile for a second order phase transition, and, if so, what is the exponent z. Notice that the system size M has been absorbed after appropriate rescalings of δ and δ * by the factor M 3 . Eq.(19) contains 3 adjustable parameters, δ * , c and z, which makes it difficult to unambiguously determine the exponent z. As pointed out in [35], if Eq. (19) holds, then for δ > δ * . Plotting this quantity against δ produces an easier fitting problem because the amplitude c, cancels out, the fitting becomes linear and the value of δ * corresponds to the point at which the fitted line crosses zero. The values of dρ dδ would ideally be obtained independently from the values for ρ. In our case, this was not possible so we obtained them by locally interpolating the measured values of ρ using M athematica and differentiating the result. The outcome of this analysis is shown in Fig. 9. We can see that a case can be made for a linear fit in the neighbourhood of δ * . Taking δ * = 0.815 (the point at which the straight line fit obviously starts to fail), and fitting the data over the range [δ * : 2] gives the fit shown in the figure. The slope gives a value of z ≈ 1 2 . This is standard Landau value for the mean field theory of a second order phase transition of a scalar field. These results, while suggestive, are far from definitive. A more detailed numerical study in the vicinity of δ * will be required before we can start putting estimates of uncertainty on these values. For the purposes of comparision, Eq. (19) with the best fit values of the parameters, is plotted with the original data in the inset of Fig.6 (solid line).

IV. CONCLUSIONS AND OUTLOOK
To conclude, we have presented a kinematic analysis of the properties of quasi-resonant triads in the CHM equation. We described the analytic form of the quasiresonant set defined by the quasi-resonance conditions, Eq. (4) as a function of the resonance broadening, δ. We found that they have non-trivial geometric shape and are not well described as simply thickened versions of the exact resonant manifold as is commonly assumed. In particular, we found that the quasi-resonant set becomes unbounded for above a critical value of δ and that this can occur for arbitrarily small values of δ as we consider modes approaching the zonal axis. We then conducted an in-depth numerical study of the structure of quasiresonant clusters as a function of δ and identified a percolation transition as δ is increased. At the transition, a large cluster is formed which contains a finite fraction of all the modes in the system. For a system containing M 2 modes, the value of the percolation threshold decreases as M −3 . This scaling results from the ease with which large scale zonal modes interact with small scale meridional modes, a reflection of the nonlocality of Rossby wave turbulence.
We speculate that the percolation transition corresponds dynamically to the transition from mesoscopic to classical wave turbulence. The fact that a percolation transition exists is consistent with earlier work on capillary waves [26]. In fact, we believe that this transition is not a special feature of the CHM dispersion relation and is generic [36]. Furthermore our results provide circumstantial support for the sandpile picture of mesoscopic wave turbulence suggested in [30] since a small change in resonance broadening in the vicinity of δ * can trigger a transition from a state which cannot support an energy cascade to one which can. In order to better understand these issues, we believe that it is important to move beyond the kinematic picture of resonance broadening and attempt to devise methods of studying these effects dynamically.