From Locking to Collective Chaos: A Reproducible Review of Kaneko's Circle Maps, Torus Instabilities, Coupled Map Lattices and Globally Coupled Maps

2.1 Circle map and Arnold tongues #

The circle map is defined by

Figure 2 shows the full tongue structure computed on a $2000\times 1000$ grid in the $(\mathrm{\Omega },K)$ plane. Each tongue emanates from a rational $\mathrm{\Omega }=p/q$ at $K=0$ and widens as $K$ increases, with the width proportional to $K^q$ for small $K$ (Arnol’d 1961). The dominant tongues $0/1$, $1/3$, $1/2$, $2/3$, and $1/1$ are visible; between them, the Farey-mediant tongues $1/4$, $2/5$, $3/7,\dots$ nest self-similarly. Below the critical line $K_c=1/(2\pi )\approx 0.1592$, the map is a diffeomorphism of the circle and tongues do not overlap: between any two tongues lies a gap of quasi-periodic (irrational $\rho$) motion. Above $K_c$, tongues overlap, the map ceases to be invertible, and chaotic dynamics are commonly observed in the overlap regions.

sec02_circle_map/arnold_tonguescodeinteractcopy linkenlarge

How this figure was computed

# Reduced preview (grid 10 x 5). Full figure: dynachaos run sec02_circle_map
import numpy as np

from dynachaos.maps.circle_map import rotation_number


def is_locked(r: float, max_denom: int = 5, tol: float = 1e-2) -> bool:
    return any(abs(r - p / q) < tol for q in range(1, max_denom + 1) for p in range(q + 1))


Omega_values = np.linspace(0.0, 1.0, 10)
K_values = np.linspace(0.0, 0.3, 5)

for K in K_values:
    row = [rotation_number(K, D=Omega, n_transient=200, n_iter=1000) for Omega in Omega_values]
    locked = sum(is_locked(r) for r in row)
    print(f"K={K:.2f}  locked points (q<=5): {locked}/{len(Omega_values)}")

copy

Artifact: arnold_tongues.npz (Omega, K, rho)

2.2 Devil’s staircase and period-adding #

Kaneko (1982) studied the circle map (1) at fixed rotation frequency $\mathrm{\Omega }=0.25$ and swept the nonlinearity $K$. The resulting plot (Figure 3): within the pre-chaotic regime ($K\lesssim 0.182$), $\rho (K)$ is a nonincreasing function with broad rational plateaus (the frequency lockings) separated by quasi-periodic intervals. Although the standard “complete devil’s staircase” is more often plotted as $\rho (\mathrm{\Omega })$ at fixed critical $K$, Kaneko’s $\mathrm{\Omega }=0.25$ scan isolates the same Farey-locking mechanism along a physically important one-parameter path.

Between any two rational plateaus $p/q$ and $p^{\prime }/q^{\prime }$ lies another plateau at the Farey mediant $(p+p^{\prime })/(q+q^{\prime })$. This nesting repeats at every scale, which makes the staircase a fractal and is why it is called a devil’s staircase. At the noninvertibility threshold $K_c=1/(2\pi )\approx 0.1592$, mode-locked plateaus cover full Lebesgue measure; the complement is a Cantor set of measure zero (Jensen et al. 1983, 1984). Below $K_c$, the quasi-periodic windows have positive measure; above it, they shrink to zero. Kaneko (1984e) analysed the supercritical behaviour ($K>K_c$) of disordered orbits in detail. The universality of this transition was analysed independently by Feigenbaum et al. (1982) and by Rand et al. (Rand et al. 1982) and Ostlund et al. (Ostlund et al. 1983).

For the specific scan at $\mathrm{\Omega }=0.25$, the dominant accumulation is the $1/5$ locking. Starting from $\rho =1/4$ at $K=0$ ($n=1$), the staircase passes through the period-adding sequence $2/9$, $3/14$, $4/19$, $\dots$, with general term $n/(5n-1)$ (Kaneko 1982, 1983b), before the $1/5$ plateau opens fully. The sequence accumulates at $K_{\infty }\approx 0.1567168$, determined here by bisection (dynachaos run sec02_circle_map) on the lift locking condition ${\theta }_{n+5}-{\theta }_n=1$ after ${10}^5$ transient iterates as the parameter where the boundary of the $1/5$ locking tongue reaches $\mathrm{\Omega }=0.25$. Between successive members of this sequence lie narrower Farey-descendant windows, so the local structure still reflects the standard devil’s-staircase hierarchy.

2.3 Diagnostic spotlight: rotation number and Lyapunov exponent #

Figure 3 shows the staircase computed with $2\times {10}^5$ parameter values (compared to Kaneko’s $\sim {10}^2$), together with the Lyapunov exponent $\lambda$ in a separate diagnostic panel. The rotation-number plot isolates the period-adding sequence approaching the $1/5$ locking, while the Lyapunov panel separates the invertible quasi-periodic regime ($\lambda \approx 0$), the stable locked windows ($\lambda <0$), and the chaotic regime that appears after the $1/5$ orbit loses stability.

The Lyapunov exponent is the natural diagnostic for the circle map: it directly distinguishes locked motion ($\lambda <0$) from quasi-periodic motion ($\lambda =0$) and chaos ($\lambda >0$), providing a one-number summary of the dynamics at each parameter value.

sec02_circle_map/devils_staircasecodeinteractcopy linkenlarge

How this figure was computed

# Reduced preview (20 values). Full figure: dynachaos run sec02_circle_map
import numpy as np

from dynachaos.maps.circle_map import bisect_k_inf, rotation_number

Omega = 0.25
K_values = np.linspace(0.0, 1.0, 20)
rho = [rotation_number(K, D=Omega, n_transient=500, n_iter=2000) for K in K_values]

for K, r in zip(K_values[::4], rho[::4], strict=False):
    print(f"K={K:.2f}  rho={r:.4f}")

# K_infty, the onset of the 1/5 lock, by bisection on the lift (the figure marks it; the
# full run uses 50 steps and a 100,000-step transient).
k_inf, low, high = bisect_k_inf(n_steps=12, n_transient=20_000)
print(f"K_inf = {k_inf:.6f} (bracket width {high - low:.1e})")

copy

Artifact: devils_staircase.npz (A, rho, lam)

The self-similar locking hierarchy becomes clearer when the staircase is resolved near the onset of the $1/5$ plateau. Figure 4 zooms into the $K\in [0.125,0.158]$ region containing the dominant period-adding sequence $2/9\to 3/14\to 4/19\to \cdots \to 1/5$, computed with ${10}^5$ parameter values. At this scale the successive rational plateaus converging on $1/5$ are cleanly separated, and narrower subwindows appear between them according to the Farey-tree construction. The point of the zoom is therefore not a generic “between $1/3$ and $1/2$” staircase, but the specific accumulation sequence Kaneko analysed in the $\mathrm{\Omega }=0.25$ problem (Kaneko 1982, 1983b).

sec02_circle_map/staircase_zoomcodeinteractcopy linkenlarge

How this figure was computed

# Reduced preview (K in [0.125, 0.158], 10 values). Full figure: dynachaos run sec02_circle_map
# Below chaos onset (K ~ 0.182) lambda <= 0 here: mode-locked, regular motion.
import numpy as np

from dynachaos.maps.circle_map import lyapunov_exponent, rotation_number

Omega = 0.25
K_values = np.linspace(0.125, 0.158, 10)
for K in K_values:
    rho = rotation_number(K, D=Omega, n_transient=500, n_iter=2000)
    lam = lyapunov_exponent(K, D=Omega, n_transient=500, n_iter=2000)
    print(f"K={K:.3f}  rho={rho:.4f}  lambda={lam:+.4f}")

copy

Artifact: staircase_zoom.npz (A, rho)

3.1 The delayed logistic map (canonical model) #

Kaneko (1984c) studied torus instabilities in their purest form using the two-dimensional delayed logistic map:

3.2 Torus oscillation #

As $D$ increases beyond $D_c$, the torus initially grows smoothly. But at a critical value of $D$, the torus begins to oscillate: the closed curve develops wiggles along its length. Kaneko (1984c) traced this to the damped oscillation of the unstable manifold of a periodic saddle near a frequency locking, and read it as an instability along the amplitude direction of the torus, distinct from both frequency locking and torus doubling.

Kaneko describes the progression as $D$ increases (his Fig. 1 (Kaneko 1984c)): oscillation of the torus, then frequency lockings, then chaos emerging via a locking, then a widening belt-like chaotic attractor, until the unstable fixed point becomes a snap-back repeller. Figure 5 shows attractors along this sequence at $\alpha =0.3$, and Figure 7 adds two measurements: panel (a) shows the Lyapunov spectrum (${\lambda }_1\approx 0$ in the torus regime and ${\lambda }_1>0$ in the chaotic regime); panel (b) shows the power spectrum of $x_n$ (dynachaos run sec05_oscillation) at $D=1.70$ and $D=1.86$. At both values every peak is a harmonic $k{f}_0$ (mod 1, folded into $[0,0.5]$) of the single rotation frequency $f_0\approx 0.225$: an invariant circle of a two-dimensional map carries one frequency. From $D=1.70$ to $D=1.86$ the higher harmonics grow by about one decade, which records the wrinkling of the circle. No independent modulation frequency appears, so in this map the oscillation of the torus is seen in its shape and in the growth of the harmonics, not as sidebands at a second frequency.

sec05_oscillation/attractorscodeinteractenlarge

How this figure was computed

# Reduced preview (3 D values). Full figure: dynachaos run sec05_oscillation
# Attractor portraits of the delayed logistic map at alpha=0.3.
from dynachaos.maps.delayed_logistic import compute_attractor

alpha = 0.3
for D, label in ((1.55, "torus"), (1.94, "chaos"), (2.16, "developed chaos")):
    traj = compute_attractor(alpha, D, n_transient=3000, n_plot=2000)
    print(f"D={D:.2f} ({label})  y range=[{traj[:, 1].min():.3f}, {traj[:, 1].max():.3f}]")

copy

Artifact: attractors.npz (D_values, A)

The transition from locking to chaos in panels (f)–(h) of Figure 5 compresses a rich micro-structure into three frames. Figure 6 zooms into $D\in [1.86,1.95]$ with eight closely-spaced $D$ values, each plotted with ${10}^5$ iterates. The sequence makes clear that the transition is not a single bifurcation but a compressed interval containing near-locking, relocking, and only then the loss of smooth invariant curves. Beyond $D\approx 1.93$, the attractor develops fractal filaments before settling into a fully chaotic cloud at the far end of the window.

sec05_oscillation/locking_sequencecodeenlarge

How this figure was computed

# Reduced preview (3 D values). Full figure: dynachaos run sec05_oscillation
# Dynamics across the locking-to-chaos transition, D in [1.86, 1.95] at alpha=0.3;
# lambda_1 turns positive as the frequency-locked torus gives way to chaos.
import numpy as np

from dynachaos.diagnostics.lyapunov import lyapunov_spectrum
from dynachaos.maps.delayed_logistic import delayed_logistic, delayed_logistic_jac

alpha = 0.3
for D in (1.860, 1.930, 1.950):
    fp = (np.sqrt(1.0 + 4.0 * D) - 1.0) / (2.0 * D)
    x0 = np.array([fp + 0.01, fp - 0.01])

    def f(state, D=D):
        return delayed_logistic(state, alpha, D)

    def jac(state, D=D):
        return delayed_logistic_jac(state, alpha, D)

    spectrum = lyapunov_spectrum(f, jac, x0, n_iter=4000, n_transient=2000)
    print(f"D={D:.3f}  lambda_1={spectrum[0]:+.4f}")

copy

Artifact: locking_sequence.npz (D_values, A)

sec05_oscillation/lyapunov_vs_Dcodeinteractcopy linkenlarge

How this figure was computed

# Reduced preview (3 D values, 2 power spectra). Full figure: dynachaos run sec05_oscillation
# Panel (a): Lyapunov spectrum of the delayed logistic map vs D at alpha=0.3.
# Panel (b): Power spectrum S(f) harmonics/sidebands (f0 approx 0.225) at D=1.70 and D=1.86.
import numpy as np

from dynachaos.diagnostics.lyapunov import lyapunov_spectrum
from dynachaos.maps.delayed_logistic import compute_psd, delayed_logistic, delayed_logistic_jac

alpha = 0.3
for D in (1.6, 1.95, 2.2):
    fp = (np.sqrt(1.0 + 4.0 * D) - 1.0) / (2.0 * D)
    x0 = np.array([fp + 0.01, fp - 0.01])

    def f(state, D=D):
        return delayed_logistic(state, alpha, D)

    def jac(state, D=D):
        return delayed_logistic_jac(state, alpha, D)

    spectrum = lyapunov_spectrum(f, jac, x0, n_iter=4000, n_transient=2000)
    print(f"D={D:.2f}  lambda_1={spectrum[0]:+.4f}  lambda_2={spectrum[1]:+.4f}")

# Power spectrum peak / harmonic check at D = 1.70 and D = 1.86 (rotation f0 approx 0.225)
for d_val in (1.70, 1.86):
    freqs, psd = compute_psd(alpha, d_val, n_transient=2000, n_record=2048)
    print(f"D={d_val:.2f}  dominant harmonic={freqs[np.argmax(psd)]:.3f}")

copy

Artifact: lyapunov_vs_D.npz (D, spectra, spectra_err, A, freqs_170, psd_170, freqs_186, psd_186)

3.3 Torus doubling (Maps I and IV) #

Period doubling, the cascade $1\to 2\to 4\to 8\to \cdots$ whose universality Feigenbaum established (Feigenbaum 1978), is the best-known route to chaos in one-dimensional maps. Arnéodo, Coullet, and Arnéodo et al. (1983) discovered that an invariant torus can undergo an analogous cascade of period doublings in continuous-time flows. Kaneko (1983a) then demonstrated that discrete maps can exhibit the same torus-doubling cascade. Just as a fixed point can bifurcate into a period-2 orbit (two alternating points), a smooth invariant circle (${\mathbb{𝕋}}^1$ in the map; “torus” in Kaneko’s flow language) can bifurcate into a doubled torus, on which the trajectory alternates between two closed curves on successive iterates. A further doubling produces a $4\times$-torus, then an $8\times$-torus, and so on until the tori become so convoluted that they merge into a chaotic attractor.

Kaneko studied four maps exhibiting torus doubling (Kaneko 1983a). Maps (II) and (III) are intermediate cases that interpolate between the three-dimensional Map (I) and the four-dimensional Map (IV); we focus on the two extremes. Map (I) is a three-dimensional delayed logistic:

The parameter $D$ controls how strongly the delayed variable drives the system. As $D$ increases at fixed $\alpha$, the torus undergoes successive doublings: $1\times \to 2\times \to 4\times \to \cdots$, accumulating at a critical value $D_{\infty }$ beyond which the motion is chaotic.

Kaneko found that the doubling cascade obeys a scaling relation analogous to Feigenbaum’s:

sec04_doubling/map_I_attractorscodeinteractcopy linkenlarge

How this figure was computed

# Reduced preview (3 D values, 2000 points). Full figure: dynachaos run sec04_doubling
# Torus doubling in Map (I), Eq. (8), at alpha = 0.4: 1x torus, 2x torus, and chaos (Dc = 2.163).
# The static figure uses 10^5 iterates after a 20,000-step transient.
import numpy as np

from dynachaos.maps.torus_doubling import iterate_map, map_I

alpha = 0.4
x0 = np.array([0.5, 0.5, 0.5])
for D, label in ((2.11, "torus"), (2.16, "2x torus"), (2.19, "chaos")):
    traj = iterate_map(map_I, x0, alpha, D, n_transient=3000, n_plot=2000)
    print(f"D={D:.2f} ({label})  X range=[{traj[:, 0].min():.3f}, {traj[:, 0].max():.3f}]")

copy

Artifact: map_I_attractors.npz (D_values, A, scan_D, scan_lyap)

sec04_doubling/map_IV_lyapunovcodeinteractcopy linkenlarge

How this figure was computed

# Reduced preview (3 D values). Full figure: dynachaos run sec04_doubling
# Lyapunov spectrum of Map (IV), the 4D delayed logistic map, vs D (alpha=0.3).
import numpy as np

from dynachaos.diagnostics.lyapunov import lyapunov_spectrum
from dynachaos.maps.torus_doubling import map_IV, map_IV_jac

alpha = 0.3
x0 = np.array([0.5, 0.45, 0.52, 0.48])
for D in (1.50, 1.5206, 1.5212):

    def f(state, D=D):
        return map_IV(state, alpha, D)

    def jac(state, D=D):
        return map_IV_jac(state, alpha, D)

    spectrum = lyapunov_spectrum(f, jac, x0, n_iter=3000, n_transient=2000)
    print(f"D={D}  lambda_1={spectrum[0]:+.4f}  lambda_2={spectrum[1]:+.4f}")

copy

Artifact: map_IV_lyapunov.npz (D, spectra)

3.4 Fractalization of torus #

At the boundary between quasi-periodic motion and chaos, Kaneko (1984b) found no single breaking event but a continuous process: the torus fractalizes.

On a smooth torus, the attractor is a one-dimensional curve (a closed orbit in the Poincaré section). Its fractal dimension is exactly 1. As the control parameter approaches the chaos threshold, the curve develops wrinkles, fine-scale oscillations that fold it at ever smaller length scales. The fractal dimension increases continuously from 1 toward $\sim 1.3$–$1.5$ at the onset of chaos.

At this stage, the attractor is a wrinkled invariant set whose maximal Lyapunov exponent has not yet turned positive. Classic strange nonchaotic attractors (SNAs) arise in quasiperiodically forced systems (Grebogi et al. 1984; Kuznetsov et al. 1995); here, because the delayed logistic map is an autonomous two-dimensional system and whether ${\lambda }_1\le 0$ holds strictly on the wrinkled set is not established here, we avoid designating it an SNA and treat it simply as a fractal invariant set.

sec07_fractalization/fractal_attractorscodeenlarge

How this figure was computed

# Reduced preview (3 D values, 2000 points). Full figure: dynachaos run sec07_fractalization
# Attractor portraits of the delayed logistic map (alpha=0.3) showing progressive fractalization.
from dynachaos.maps.fractalization import iterate

alpha = 0.3
for D, label in ((1.75, "smooth torus"), (1.92, "fractal torus"), (1.945, "chaos")):
    traj = iterate(alpha, D, n_transient=2000, n_record=2000)
    x_lo, x_hi = traj[:, 0].min(), traj[:, 0].max()
    print(f"D={D} ({label})  points={traj.shape[0]}  x range=[{x_lo:.3f}, {x_hi:.3f}]")

copy

Artifact: fractal_attractors.npz (D_values, A)

3.5 Diagnostic spotlight: correlation dimension and SALI #

We quantify the fractalization using the Grassberger–Procaccia correlation dimension $D_2$ (Grassberger and Procaccia 1983). The correlation integral

Figure 11 plots $D_2$ versus $D$. Across the main fractalization window $D\in [1.92,1.95]$, the mean correlation dimension is $D_2=1.24\pm 0.07$ (mean $\pm$ standard deviation across parameter values in this window, computed by dynachaos run sec07_fractalization). Individual fits that yield $D_2<1$ are identified as failed fits arising from finite-sample and scaling-window sensitivity and are excluded from the window average: 3 of the 101 fits fail ($D=1.700$, $1.913$, $1.940$, with $D_2=0.99$, $0.98$, $0.84$), one of them inside the window, so the window mean uses 9 of 10 values. The quoted uncertainty is the scatter over the window, not a fit error; the individual regression errors are $\approx 0.01$. Below the window ($D<1.91$) the curve sits at $D_2\approx 1.03$ (a smooth closed curve) apart from one outlier ($D_2=1.20$ at $D=1.799$), and above it rises to $D_2\approx 1.3$ near $D=2$.

The SALI (Smaller Alignment Index) (Skokos 2001) provides an alignment diagnostic based on two tangent vectors under the linearised dynamics: $\mathrm{\text{SALI}}=\min (\parallel {\hat{v}}_1+{\hat{v}}_2\parallel ,\parallel {\hat{v}}_1-{\hat{v}}_2\parallel )$. The generalised index ${\mathrm{\text{GALI}}}_k$ (Skokos et al. 2007) is the volume of the parallelepiped spanned by $k$ unit deviation vectors; ${\mathrm{\text{GALI}}}_2$ and SALI are functions of each other. In a dissipative map, phase-space contraction pulls deviation vectors onto the neutral and expanding directions of the tangent map, even for regular orbits (Moges et al. 2025). What follows depends on how many such directions exist. Write ${\lambda }_1\ge {\lambda }_2\ge \cdots$ for the Lyapunov spectrum. For random initial vectors, $\ln {\mathrm{\text{GALI}}}_k$ falls at the rate ${\sum }_{i=2}^k({\lambda }_i-{\lambda }_1)$ per iterate, and SALI falls at the rate of ${\mathrm{\text{GALI}}}_2$, ${\lambda }_2-{\lambda }_1$, whenever that is negative. We test this on the four-dimensional coupled delayed logistic map of Section 5.1.

On a 2-torus of the map, ${\lambda }_1={\lambda }_2=0$. The contracting directions (${\lambda }_3$, ${\lambda }_4<0$) die out, and the two deviation vectors end in the two-dimensional tangent plane of the torus. There they need not align, so SALI and ${\mathrm{\text{GALI}}}_2$ stay bounded, while ${\mathrm{\text{GALI}}}_3$ and ${\mathrm{\text{GALI}}}_4$ decay at the rates ${\lambda }_3$ and ${\lambda }_3+{\lambda }_4$. On partial locking (${\lambda }_1=0>{\lambda }_2$) and on chaos (${\lambda }_1>0$) the first two vectors align, and SALI decays at the rate ${\lambda }_2-{\lambda }_1$.

Figure 12 tests these predictions over ${10}^6$ iterates on two orbits with two zero exponents, one partially locked orbit, and one chaotic orbit. On the two tori, SALI and ${\mathrm{\text{GALI}}}_2$ stay between about ${10}^{-3}$ and $1$. They fluctuate over about three decades and return to values near $1$; they do not converge to a constant. The fitted slope of $\ln {\mathrm{\text{GALI}}}_2$ over $n={10}^5$ to ${10}^6$ is $+5\times {10}^{-8}$ and $-1.8\times {10}^{-8}$ per iterate (seed 0). On the same orbits ${\mathrm{\text{GALI}}}_3$ decays at $-0.2161$ and $-0.2146$ per iterate (predicted $-0.2161\pm 1.2\times {10}^{-4}$ and $-0.2146\pm 1.1\times {10}^{-4}$), and ${\mathrm{\text{GALI}}}_4$ at $-0.6252$ and $-0.6150$ (predicted $-0.6253\pm 1.7\times {10}^{-4}$ and $-0.6151\pm 1.4\times {10}^{-4}$). On the partially locked orbit ${\mathrm{\text{GALI}}}_2$ decays at $-1.8\times {10}^{-3}$ per iterate (predicted $-1.7745\times {10}^{-3}\pm 6.7\times {10}^{-6}$), and on the chaotic orbit at $-0.0370$ (predicted $-0.0368\pm 2.5\times {10}^{-4}$). SALI and ${\mathrm{\text{GALI}}}_2$ decay on both of these orbits, so SALI decay alone does not separate partial locking from chaos; the decay rate does ($-1.8\times {10}^{-3}$ against $-0.037$ per iterate). Over ${10}^6$ iterates the run cannot exclude a decay of SALI and ${\mathrm{\text{GALI}}}_2$ on the tori slower than about ${10}^{-5}$ per iterate; the rates predicted from the spectrum, $-9\times {10}^{-6}\pm {10}^{-5}$ and $-9.6\times {10}^{-6}\pm 1.2\times {10}^{-5}$, are consistent with zero. The result is for one map ($\alpha =0.4$, $\epsilon =5\times {10}^{-3}$), two tori, one locked orbit, and one chaotic orbit.

sec07_fractalization/correlation_dimensioncodeinteractcopy linkenlarge

How this figure was computed

# Reduced preview (3 D values, 3000 points). Full figure: dynachaos run sec07_fractalization
# Correlation dimension D_2(D) of the delayed logistic map (alpha = 0.3).
# Full figure uses 101 values; window [1.92, 1.95] gives D_2 = 1.24 +/- 0.07 (excluding D_2 < 1).
from dynachaos.diagnostics.correlation import correlation_dimension
from dynachaos.maps.fractalization import iterate

alpha = 0.3
for D in (1.75, 1.90, 1.945):
    traj = iterate(alpha, D, n_transient=2000, n_record=3000)
    res = correlation_dimension(traj, n_r=30, max_pairs=200_000)
    D2 = res[0]
    print(f"D={D}  D2={D2:.3f}")

copy

Artifact: correlation_dimension.npz (D, D2, D2_err, A, valid_mask, failed_mask, mean_window_D2, std_window_D2, window_min, window_max)

sec11_diagnostics/alignment_three_toruscodeenlarge

How this figure was computed

# Reduced preview (2 orbits, 20 000 iterates). Full figure: dynachaos run sec11_diagnostics
# GALI_k decays at the rate sum_{i=2..k} (lambda_i - lambda_1) set by the Lyapunov spectrum.
import numpy as np

from dynachaos.diagnostics.alignment_three_torus import predicted_rates
from dynachaos.diagnostics.lyapunov import lyapunov_spectrum_blocks
from dynachaos.diagnostics.sali_gali import alignment_indices
from dynachaos.maps.coupled_delayed import coupled_delayed, coupled_delayed_jac

alpha, eps = 0.4, 5e-3
x0 = np.array([0.5, 0.5, 0.3, 0.3])
n_iter, n_transient = 20_000, 5000

for DB, label in ((2.37, "partial locking"), (2.55, "chaos")):
    DA = DB + 0.1

    def f(s, DA=DA, DB=DB):
        return coupled_delayed(s, alpha, DA, DB, eps)

    def jac(s, DA=DA, DB=DB):
        return coupled_delayed_jac(s, alpha, DA, DB, eps)

    spectrum, se = lyapunov_spectrum_blocks(f, jac, x0, n_iter, n_transient)
    gali = alignment_indices(f, jac, x0, k=4, n_iter=n_iter, n_transient=n_transient)
    rate, _ = predicted_rates(spectrum, se)
    fit = (gali[-1, 1:] - gali[n_iter // 2, 1:]) * np.log(10.0) / (n_iter - n_iter // 2)
    print(f"DB={DB} ({label})  lambda={np.round(spectrum, 4)}")
    print(f"  GALI_2..4 rate  fitted={np.round(fit, 4)}  predicted={np.round(rate, 4)}")

copy

Artifact: alignment_three_torus.npz (DB_values, regimes, n_samples, spectrum, spectrum_se, log10_gali, log10_sali, decade_min, decade_max, rate_fit, rate_fit_early, rate_pred, rate_pred_se)

4.1 The coupled logistic map (canonical model) #

To study two-frequency quasi-periodicity (an invariant circle ${\mathbb{𝕋}}^1$ in this 2D return map), Kaneko introduced the coupled logistic map (Kaneko 1983c):

4.2 Phase diagram and symmetry-broken attractors #

Kaneko showed that as $a$ increases at fixed $\epsilon$, the attractor can pass through a sequence of symmetry-preserving and symmetry-broken regimes:

1. Symmetric regular dynamics: fixed or periodic orbits with $x\approx y$.

2. Locked windows: regular periodic windows embedded in the quasiperiodic route.

3. Symmetry breaking: the $x\leftrightarrow y$ symmetry breaks, creating two coexisting attractors related by exchange.

4. Chaotic bands: irregular sets with positive finite-time expansion rate.

Figure 13 therefore serves as a locator map: panel (a) uses the symmetry order parameter $m=⟨|x-y|⟩$, while panel (b) uses a finite-time estimate of the largest Lyapunov exponent ${\lambda }_1$. The two fields separate regions where the dynamics remains close to the invariant diagonal from regions where symmetry is broken and the motion becomes irregular. The dashed line at $\epsilon =0.1$ marks the slice used in Figure 14. That attractor gallery follows a representative branch through the route $2T\to 4T\to 8T\to 8C\to 4C$; for the first two panels we deliberately start off the diagonal to expose the coexisting symmetry-broken tori rather than the synchronized cycles.

sec03_transition/phase_diagramcodeinteractcopy linkenlarge

How this figure was computed

# Reduced preview (grid 4 x 2). Full figure: dynachaos run sec03_transition
# Symmetry breaking <|x-y|> for the coupled logistic map over an (a, eps) grid.
import numpy as np

from dynachaos.maps.coupled_logistic import coupled_logistic

for eps in (0.05, 0.1):
    for a in (1.0, 1.2, 1.35, 1.5):
        state = np.array([0.1, 0.2])
        for _ in range(5000):
            state = coupled_logistic(state, a, eps)
        asym = 0.0
        for _ in range(2000):
            state = coupled_logistic(state, a, eps)
            asym += abs(state[0] - state[1])
        asym /= 2000
        print(f"eps={eps:.2f}  a={a:.2f}  <|x-y|>={asym:.4f}")

copy

Artifact: phase_diagram.npz (A, D, asym, lyap, schema_version)

sec03_transition/attractorscodeinteractcopy linkenlarge

How this figure was computed

# Reduced preview (3 states, 2000 points). Full figure: dynachaos run sec03_transition
# Representative attractor portraits at eps=0.1 along the broken-symmetry route (2T -> 4T -> 4C).
# Panels (a)-(b) use initial condition (0.1, 0.6); panels (c)-(f) use (0.1, 0.2).
import numpy as np

from dynachaos.maps._iter import trajectory_after_transient
from dynachaos.maps.coupled_logistic import coupled_logistic

eps = 0.1
cases = (
    (1.10, "2T", np.array([0.1, 0.6])),
    (1.25, "4T", np.array([0.1, 0.6])),
    (1.373, "4C", np.array([0.1, 0.2])),
)

for a, label, init_state in cases:
    traj = trajectory_after_transient(
        init_state,
        lambda state, a=a: coupled_logistic(state, a, eps),
        n_transient=5000,
        n_record=2000,
    )
    print(f"a={a:.4f} ({label})  x range=[{traj[:, 0].min():.3f}, {traj[:, 0].max():.3f}]")

copy

Artifact: attractors.npz (A_values, labels, initial_states, x_limits, y_limits, D, schema_version)

4.3 Self-similar basin structure #

Takesue and Kaneko (1984) showed that the basins of attraction for the two symmetry-broken attractors display an alternating stripe structure: the boundary between basins is a fractal. Figure 15 reproduces this structure at $a=1.35344$ on an $800\times 800$ grid. At this parameter the map has two coexisting asymmetric 4-cycles in Kaneko’s reduced notation (period-32 cycles for the original map, where the reduced notation counts the four macroscopic clusters visited across the rapid underlying 8-cycle dynamics) together with the invariant diagonal $x_0=y_0$, which flows to the synchronized cycle. The two zoom panels show stripe-like interleaving of the two basins near the diagonal at two successive magnifications. Classification evaluates the minimum Euclidean distance to the two asymmetric reference orbits over a full 32-step period rather than relying on the single-step sign of $x-y$; points on the invariant diagonal remain exactly equidistant and stay grey (dynachaos run sec03_transition).

sec03_transition/basinscodeenlarge

How this figure was computed

# Reduced preview (9 initial states). Full figure: dynachaos run sec03_transition
# Basins of attraction for two mirror-image asymmetric period-32 orbits at a = 1.35344, eps = 0.1.
# Orbits are classified by distance to the two reference cycles; the diagonal x = y is symmetric.
import numpy as np

from dynachaos.maps.coupled_logistic import _find_reference_orbit, coupled_logistic

a, eps = 1.35344, 0.1
ref_A = _find_reference_orbit(a, eps, 0.1, 0.6, n_transient=50_000, period=32)
ref_B = ref_A[:, ::-1]

for x0 in np.linspace(-0.05, 0.05, 9):
    state = np.array([x0, -x0])
    for _ in range(5000):
        state = coupled_logistic(state, a, eps)
    orbit = [state]
    for _ in range(31):
        state = coupled_logistic(state, a, eps)
        orbit.append(state)
    orbit = np.array(orbit)
    dA = min(np.min(np.sum((orbit - r) ** 2, axis=1)) for r in ref_A)
    dB = min(np.min(np.sum((orbit - r) ** 2, axis=1)) for r in ref_B)
    if np.isclose(dA, dB, atol=1e-12):
        basin = "diagonal"
    else:
        basin = "A" if dA < dB else "B"
    print(f"x0={x0:+.3f}  basin={basin}  dist_A={dA:.2e}  dist_B={dB:.2e}")

copy

Artifact: basins.npz (x, y, basin, A, D, x_zoom1, y_zoom1, basin_zoom1, rect_zoom1, x_zoom2, y_zoom2, basin_zoom2, rect_zoom2)

4.4 Diagnostic spotlight: the 0-1 test for chaos #

The 0-1 test (Gottwald and Melbourne 2004, 2009) provides a binary answer from a scalar time series alone: the test statistic $K_{01}$ approaches 0 for regular dynamics and 1 for chaotic dynamics. The method projects the time series onto a 2D extension $(p_n,q_n)$ via

Figure 16 shows $K_{01}$ versus $a$ for the logistic map $f(x)=1-a{x}^2$. The periodic windows (where $K_{01}\approx 0$) are sharply delineated from the chaotic bands (where $K_{01}\approx 1$), agreeing with the sign of the Lyapunov exponent $\lambda (a)$ but without requiring any knowledge of the map’s equations or derivatives. Frequencies $c$ are sampled randomly from $(\pi /5,4\pi /5)$ with $N_c=50$ to avoid false chaos from resonances with low-period orbits (such as $c=2\pi /3$ resonant with period-3 cycles), taking the median correlation value as $K_{01}$.

sec11_diagnostics/test01_sweepcodeinteractcopy linkenlarge

How this figure was computed

# Reduced preview (4 a values, 5000 points, n_c=50). Full figure: dynachaos run sec11_diagnostics
# 0-1 test statistic K_01 across the logistic map's route to chaos.
import numpy as np

from dynachaos.diagnostics.zero_one_test import zero_one_statistic
from dynachaos.maps.primitives import logistic

for a in (1.2, 1.6, 1.9, 2.0):
    x = 0.1
    for _ in range(2000):
        x = logistic(x, a)
    series = np.empty(5000)
    for i in range(5000):
        x = logistic(x, a)
        series[i] = x
    K = zero_one_statistic(series, n_c=50)
    print(f"a={a:.2f}  K_01={K:.3f}")

copy

Artifact: test01_sweep.npz (a, K)

5.1 The coupled delayed logistic map #

A three-torus ${\mathbb{𝕋}}^3$ is a quasi-periodic attractor with three independent, incommensurate frequencies ${\omega }_1,{\omega }_2,{\omega }_3$. It is the dynamical object that appears in Landau’s route to turbulence (Landau 1944) after three Hopf bifurcations. Ruelle and Takens (1971) argued that ${\mathbb{𝕋}}^3$ should be structurally unstable: generic perturbations destroy it and replace it with a strange attractor. Newhouse et al. (1978) proved this for certain perturbations. Kaneko’s question was which objects replace ${\mathbb{𝕋}}^3$.

To study ${\mathbb{𝕋}}^3$, Kaneko constructed a four-dimensional map (Kaneko 1984a):

5.2 Lyapunov counting for $k$-tori in maps versus flows #

The four-dimensional Lyapunov spectrum $({\lambda }_1\ge {\lambda }_2\ge {\lambda }_3\ge {\lambda }_4)$ distinguishes the dynamical regimes:

• ${\lambda }_1={\lambda }_2=0$, ${\lambda }_3,{\lambda }_4<0$: three-frequency quasi-periodicity (a ${\mathbb{𝕋}}^2$ invariant set in the map, corresponding to three-torus motion in the underlying flow).

• One zero exponent: a single resonant frequency locking condition reduces the neutral directions by one; the remaining two-frequency motion lives on an invariant circle ${\mathbb{𝕋}}^1$ in the map.

• All exponents negative: fully locked periodic orbit.

• ${\lambda }_1>0$: chaos (the torus has collapsed).

Figure 17 shows the Lyapunov spectrum versus $D_B$ for three values of the coupling $\epsilon$. The corresponding $(x_n,z_n)$ projections (Figure 18) provide a geometric gallery at $\epsilon =5\times {10}^{-3}$: a near-three-frequency curve, resonance webs and locked circles, weak chaos, a narrow periodic window, and developed chaos.

5.3 Double devil’s staircase #

The double devil’s staircase is Kaneko’s name for the locking structure of the partially locked torus: when the three-frequency motion partially locks, the rotation number ${\rho }_{\theta }$ takes values of the form $(p+sC)/r$ (where $p,s,r\in \mathbb{\mathbb{Z}}$), with $C=(\sqrt{5}-1)/2$ the irrational modulation frequency. The staircase is “double” because its plateaus are built from two incommensurate frequencies ($1$ and $C$), in contrast to the single-generator structure of the ordinary devil’s staircase (Kaneko 1984a). Our sweep computes the rotation numbers directly; it agrees with the two-generator picture but does not decompose every plateau into its rational parts.

5.4 Modulated circle map #

To demonstrate the double staircase cleanly, Kaneko introduced the modulated circle map (Kaneko 1984a):

sec06_three_torus/lyapunov_vs_DBcodeinteractcopy linkenlarge

How this figure was computed

# Reduced preview (3 D_B values). Full figure: dynachaos run sec06_three_torus
# Lyapunov spectrum of the 4D coupled delayed logistic map vs D_B (alpha=0.4, eps=5e-3).
import numpy as np

from dynachaos.diagnostics.lyapunov import lyapunov_spectrum
from dynachaos.maps.coupled_delayed import coupled_delayed, coupled_delayed_jac

alpha, eps = 0.4, 5e-3
x0 = np.array([0.5, 0.5, 0.3, 0.3])
for DB in (2.37, 2.478, 2.55):
    DA = DB + 0.1

    def f(state, DA=DA, DB=DB):
        return coupled_delayed(state, alpha, DA, DB, eps)

    def jac(state, DA=DA, DB=DB):
        return coupled_delayed_jac(state, alpha, DA, DB, eps)

    spectrum = lyapunov_spectrum(f, jac, x0, n_iter=2000, n_transient=1000)
    print(f"D_B={DB:.3f}  lambda_1={spectrum[0]:+.4f}  lambda_2={spectrum[1]:+.4f}")

copy

Artifact: lyapunov_vs_DB.npz (DB, eps_values, eps_0.001_spectra, eps_0.005_spectra, eps_0.01_spectra)

sec06_three_torus/xz_projectionscodeinteractcopy linkenlarge

How this figure was computed

# Reduced preview (2 D_B values, 1000 points). Full figure: dynachaos run sec06_three_torus
# (x_n, z_n) projections of the 4D coupled delayed logistic map (alpha=0.4, eps=5e-3).
import numpy as np

from dynachaos.maps._iter import trajectory_after_transient
from dynachaos.maps.coupled_delayed import coupled_delayed

alpha, eps = 0.4, 5e-3
for DB, label in ((2.37, "near-T^3"), (2.55, "developed chaos")):
    DA = DB + 0.1
    traj = trajectory_after_transient(
        np.array([0.5, 0.5, 0.3, 0.3]),
        lambda state, DA=DA, DB=DB: coupled_delayed(state, alpha, DA, DB, eps),
        n_transient=2000,
        n_record=1000,
        project_fn=lambda state: state[[0, 2]],
    )
    print(f"D_B={DB} ({label})  x range=[{traj[:, 0].min():.3f}, {traj[:, 0].max():.3f}]")

copy

Artifact: xz_projections.npz (DB_values, labels, render_modes, schema_version)

sec06_three_torus/double_staircasecodeinteractcopy linkenlarge

How this figure was computed

# Reduced preview (6 points). Full figure: dynachaos run sec06_three_torus
# Double devil's staircase: rotation number rho_theta(D) of the modulated circle map.
import numpy as np

from dynachaos.maps.modulated_circle import C_GOLDEN, rotation_numbers

A, eps = 0.10, 0.05
for D in np.linspace(0.0, 1.0, 6):
    rho_theta, rho_phi = rotation_numbers(A, C_GOLDEN, D, eps, n_transient=500, n_iter=2000)
    print(f"D={D:.2f}  rho_theta={rho_theta:.4f}  rho_phi={rho_phi:.4f}")

copy

Artifact: double_staircase.npz (D, rho_theta, rho_phi, A, C, eps)

The partial locking structure becomes clearer when zoomed into narrower windows of $D$. Figure 20 shows two such windows. Panel (a) resolves a narrow plateau near ${\rho }_{\theta }\approx 1/4$; panel (b) shows a broader plateau near ${\rho }_{\theta }\approx C$. Together they illustrate that the staircase contains both ordinary rational lockings and plateaus tied to the irrational modulation frequency.

sec06_three_torus/double_staircase_zoomcodeenlarge

How this figure was computed

# Reduced preview (8 values in 1/4 plateau). Full figure: dynachaos run sec06_three_torus
# Zoom into the 1/4 mode-locked plateau of the double devil's staircase.
import numpy as np

from dynachaos.maps.modulated_circle import C_GOLDEN, rotation_numbers

A, eps = 0.10, 0.05
for D in np.linspace(0.2631, 0.2654, 8):
    rho_theta, _ = rotation_numbers(A, C_GOLDEN, D, eps, n_transient=2000, n_iter=8000)
    print(f"D={D:.4f}  rho_theta={rho_theta:.5f}")

copy

5.5 Diagnostic spotlight: permutation entropy and complexity-entropy plane #

Permutation entropy (Bandt and Pompe 2002) quantifies the complexity of a time series by examining the relative ordering of $d$ consecutive values. Each window of $d$ values defines an ordinal pattern (permutation of $\{0,1,\dots ,d-1\}$). The normalised Shannon entropy of the pattern distribution gives

For an embedding dimension $d=5$, a periodic orbit of period $p$ supports at most $p$ distinct ordinal permutations among the $5!=120$ possible patterns. When each pattern appears with equal frequency $1/p$, the normalised permutation entropy takes the exact value $H_{\mathrm{\text{PE}}}=\ln (p)/\ln (5!)$, producing discrete plateaus at $H_{\mathrm{\text{PE}}}\approx 0.145$ for period 2 ($p=2$), $\approx 0.290$ for period 4 ($p=4$), and $\approx 0.434$ for period 8 ($p=8$). These plateaus explain the discrete steps visible in Figure 21, rather than an immediate drop to zero. More irregular dynamics populates many patterns, pushing $H_{\mathrm{\text{PE}}}$ upward. Combined with the Jensen–Shannon statistical complexity $C_{\mathrm{\text{JS}}}$ (Rosso et al. 2007), the pair $(H_{\mathrm{\text{PE}}},C_{\mathrm{\text{JS}}})$ locates the dynamics in a complexity-entropy plane (Figure 22). Noise occupies the high-$H$, low-$C$ corner; chaos occupies an intermediate region; and periodic motion clusters near the lower boundary.

As a compact benchmark, we evaluate these diagnostics on scalar observables from the logistic and delayed logistic maps (Figure 21). This keeps the interpretation simple while showing how entropy-based diagnostics track the shift from regular to more irregular dynamics without requiring Jacobian estimates or explicit phase-space reconstructions.

sec11_diagnostics/permutation_entropycodeinteractcopy linkenlarge

How this figure was computed

# Reduced preview (4 a values). Full figure: dynachaos run sec11_diagnostics
# Permutation entropy H_PE(a) for the logistic map (d=5).
import numpy as np

from dynachaos.diagnostics.permutation import permutation_entropy
from dynachaos.maps.primitives import logistic

for a in (1.2, 1.6, 1.9, 2.0):
    x = 0.1
    for _ in range(2000):
        x = logistic(x, a)
    series = np.empty(2000)
    for i in range(2000):
        x = logistic(x, a)
        series[i] = x
    H = permutation_entropy(series, d=5)
    print(f"a={a:.2f}  H_PE={H:.3f}")

# Full figure: dynachaos run sec11_diagnostics

copy

Artifact: permutation_entropy.npz (a, H_logistic, D, H_delayed)

sec11_diagnostics/complexity_entropy_planecodeinteractcopy linkenlarge

How this figure was computed

# Reduced preview (4 a values). Full figure: dynachaos run sec11_diagnostics
# Complexity-entropy plane (H, C) locating logistic-map regimes (d=5).
import numpy as np

from dynachaos.diagnostics.permutation import complexity_entropy
from dynachaos.maps.primitives import logistic

for a in (1.2, 1.6, 1.9, 2.0):
    x = 0.1
    for _ in range(2000):
        x = logistic(x, a)
    series = np.empty(2000)
    for i in range(2000):
        x = logistic(x, a)
        series[i] = x
    H, C = complexity_entropy(series, d=5)
    print(f"a={a:.2f}  H={H:.3f}  C={C:.3f}")

# Full figure: dynachaos run sec11_diagnostics

copy

Artifact: complexity_entropy_plane.npz (a, H_logistic, C_logistic, D, H_delayed, C_delayed)

6.1 CML definition and three models #

A coupled map lattice (CML) is a discrete-time, discrete-space dynamical system (Kaneko 1985):

Kaneko studied three CML models with different local dynamics:

1. A piecewise map, $f(x)=x+x^2+a$ for $x<c$ and $f(x)=-3(x-c)+1+a$ for $x\ge c$, with $c=(\sqrt{5}-1)/2$, $a=-0.01$ and $g=f$.

2. Coupled circle maps with $f(\theta )=\theta +K\sin (2\pi \theta )+\mathrm{\Omega }$ and $g(\theta )=\sin (2\pi \theta )$.

3. Coupled logistic maps with $f(x)=1-a{x}^2$ and $g=f$.

6.2 Spatiotemporal intermittency and directed percolation #

All three models exhibit spatiotemporal intermittency (STI): the spacetime diagram (Figure 23) shows alternating laminar (ordered) and burst (disordered) regions. As the coupling $\epsilon$ is tuned, the fraction of bursts changes continuously, analogous to a phase transition.

Kaneko reduced the laminar/burst dynamics to a binary (cellular automaton) description (Kaneko 1985) and characterised the transition through burst-propagation speed $v_B\propto (\epsilon -{\epsilon }_c{)}^{\nu }$ with $\nu \approx 0.6$–$0.9$, and the growth of positive Lyapunov exponents. Pomeau (1986) suggested a connection to the universality class of directed percolation (DP), the standard absorbing-state transition of non-equilibrium statistical mechanics (Hinrichsen 2000). Chaté and Chaté and Manneville (1988a, 1988b) investigated this connection and found DP-like scaling features in deterministic spatiotemporal intermittency, while noting that the critical exponents may depend on whether the laminar state is strictly absorbing (Chaté and Manneville 1989b, 1989a); Kaneko himself noted that “the universality class on critical exponents seems to be different” (Kaneko 1989b).

To complement the qualitative spacetime view with quantitative measures, we extract an STI “spine” from the Kaneko Model-A CML (dynachaos run sec12_intermittency; Figure 24): the turbulent fraction ${\varphi }_{\mathrm{\text{turb}}}$ as a function of the coupling $\epsilon$, together with the laminar-domain size distribution, where an active site is defined by the local burst departure threshold $|x_n(i+1)-x_n(i)|>0.05$. The turbulent fraction is zero up to $\epsilon =0.050$, rises between $\epsilon =0.050$ and $0.065$, and reaches a burst-dominated plateau ${\varphi }_{\mathrm{\text{turb}}}\approx 0.88$ for $\epsilon \ge 0.07$. At $\epsilon =0.050$ the whole lattice is laminar, so there is no laminar-length distribution. At $\epsilon =0.060$ the distribution is broad (lengths up to 50 sites); among a power law, an exponential and a truncated power law, the power law has the lowest BIC ($\alpha \approx 1.16$), but the Kolmogorov–Smirnov distances, evaluated at each integer length, reject all three: 0.074 for the power law and the truncated power law and 0.14 for the exponential, against a 5 % critical value of $1.36/\sqrt{n}\approx 0.010$ for $n\approx 1.7\times {10}^4$ lengths (the critical value assumes independent lengths). At $\epsilon =0.080$ laminar runs are at most six sites long and the exponential fits (rate $\approx 2.0$, distance 0.002 against a critical value of 0.017). The data therefore locate the onset but do not determine the form of the laminar-length distribution there. Testing whether the system falls in the directed percolation (DP) class requires finite-size scaling at the critical coupling ${\epsilon }_c$, which is not performed here. Critical exponents of directed percolation have been measured in coupled circle-map lattices (Janaki et al. 2003; Jabeen and Gupte 2005, 2006) and in a ferrofluid experiment (Rupp et al. 2003). What remains open is a rigorous mathematical proof and whether the Model-A lattice studied here belongs to the DP universality class.

sec08_sti/spacetime_diagramscodeinteractenlarge

How this figure was computed

# Reduced preview (N = 60, 50 steps). Full figure: dynachaos run sec08_sti
# Spacetime diagrams of a coupled map lattice (Kaneko 1985 Model C, logistic local map).
from dynachaos.cml.spatiotemporal import model_C_f, simulate_cml

N = 60
for eps in (0.16, 0.20, 0.30):
    spacetime = simulate_cml(model_C_f, model_C_f, eps, N=N, n_transient=200, n_record=50)
    print(f"eps={eps:.2f}  final-row var={spacetime[-1].var():.4f}")

# Full figure: dynachaos run sec08_sti

copy

Artifact: spacetime_diagrams.npz (A_eps_0.06, B_eps_0.02, C_eps_0.16)

sec12_intermittency/sti_spinecodeenlarge

How this figure was computed

# Reduced preview (N = 256, 100 steps). Full figure: dynachaos run sec12_intermittency
# Quantitative STI spine for Kaneko Model-A CML: turbulent fraction vs coupling eps, the critical
# coupling eps_c, and maximum-likelihood fits of the laminar-length tail (power law, exponential,
# truncated power law).
# Full figure: N = 512, phi_turb at 19 couplings in [0.03, 0.12], fits at eps_c, 0.060 and 0.080.
import numpy as np

from dynachaos.cml.spatiotemporal import model_A_f, simulate_cml
from dynachaos.cml.sti_spine_figure import (
    _burst_mask,
    _periodic_run_lengths,
    fit_tail_models,
    locate_critical_coupling,
)

N = 256
sweep_eps = np.linspace(0.03, 0.12, 10)
fractions = []
for idx, eps in enumerate(sweep_eps):
    rng = np.random.default_rng(42 + idx)
    spacetime = simulate_cml(
        model_A_f,
        model_A_f,
        float(eps),
        N=N,
        n_transient=500,
        n_record=100,
        x0=rng.uniform(0, 1, N),
    )
    fractions.append(_burst_mask(spacetime).mean())
    print(f"eps={eps:.3f}  turbulent_fraction={fractions[-1]:.3f}")
print(f"eps_c = {locate_critical_coupling(sweep_eps, fractions):.3f}")

lengths = _periodic_run_lengths(~_burst_mask(spacetime))
fit = fit_tail_models(lengths, x_min=2)
print(f"laminar tail at eps=0.12: n={fit['n_tail']}  preferred={fit['preferred']}")

copy

Artifact: sti_spine.npz (schema_version, source_file, seed, model_a_parameter, display_eps, sweep_eps, turbulent_fraction, eps_c, fit_eps_values, spacetime, turbulent_mask, laminar_values_0, laminar_probs_0, laminar_values_1, laminar_probs_1, laminar_values_2, laminar_probs_2, fit_x_min, fit_n_tail, fit_alpha_pl, fit_ll_pl, fit_bic_pl, fit_ks_pl, fit_lam_exp, fit_ll_exp, fit_bic_exp, fit_ks_exp, fit_alpha_tr, fit_lam_tr, fit_ll_tr, fit_bic_tr, fit_ks_tr, fit_preferred)

6.3 Pattern dynamics (five phases) #

Kaneko (1989b) systematically classified the dynamics of the CML (13) with the logistic map $f(x)=1-a{x}^2$ and $g=f$. The canonical taxonomy in the $(a,\epsilon )$ plane contains five phases:

1. Frozen random ($a\lesssim 1.5$, moderate $\epsilon$): the lattice settles into a static, spatially disordered pattern.

2. Pattern selection ($a\approx 1.6$–$1.7$): the lattice spontaneously selects a spatially periodic pattern, breaking translational symmetry.

3. Defect turbulence ($a\approx 1.7$–$1.8$): the selected pattern is disrupted by defects (domain walls) that nucleate, propagate, and annihilate chaotically.

4. Pattern competition intermittency ($a\approx 1.8$–$1.9$): large coherent domains compete for territory, punctuated by intermittent bursts of turbulence.

5. Fully developed turbulence ($a\gtrsim 1.9$): all spatial coherence is lost.

Figure 25 shows a spatial activity map computed on a $200\times 160$ grid (versus Kaneko’s $\sim 50\times 50$): the mean difference between neighbouring sites, $⟨|x_i-x_{i-1}|⟩$, averaged over the lattice and over time. It separates spatially smooth from spatially rough states, so the figure should be read as a locator map consistent with the five-phase taxonomy, not as an automatic phase classifier. Figure 26 then shows space-amplitude plots at $\epsilon =0.1$ for five values of $a$ marked on that map.

sec09_pattern/phase_diagramcodeinteractcopy linkenlarge

How this figure was computed

# Reduced preview (3 (a, eps) states, N = 100, 500-step transient, 200 samples).
# Full figure: dynachaos run sec09_pattern (200 x 160 grid, 2000 samples after 5000 transient).
# Spatial activity <|x_i - x_{i-1}|> of the logistic coupled map lattice
# (Kaneko 1989) separates spatially smooth from spatially rough states.
import numpy as np

from dynachaos.cml.primitives import cml_step_logistic as cml_step

N = 100
rng = np.random.default_rng(42)
for a, eps in ((1.50, 0.10), (1.72, 0.10), (1.90, 0.10)):
    x = rng.uniform(-1, 1, N)
    for _ in range(500):
        x = cml_step(x, a, eps)
    activity = 0.0
    for _ in range(200):
        x = cml_step(x, a, eps)
        activity += np.mean(np.abs(x - np.roll(x, 1)))
    print(f"a={a:.2f} eps={eps:.2f}  spatial activity={activity / 200:.4f}")

copy

Artifact: phase_diagram.npz (a, eps, lam, spatial_activity)

sec09_pattern/space_amplitudecodeinteractcopy linkenlarge

How this figure was computed

# Reduced preview (2 a states, N = 100). Full figure: dynachaos run sec09_pattern
# Space-amplitude snapshots x(i) vs site i of a logistic coupled map lattice (N = 100).
import numpy as np

from dynachaos.cml.primitives import cml_step_logistic as cml_step

N = 100
rng = np.random.default_rng(42)
for a, eps, label in ((1.50, 0.10, "frozen random"), (1.90, 0.10, "fully developed turbulence")):
    x = rng.uniform(-1, 1, N)
    for _ in range(500):
        x = cml_step(x, a, eps)
    snapshot = x.copy()
    print(f"a={a:.2f} ({label})  amplitude mean={snapshot.mean():.4f} std={snapshot.std():.4f}")

copy

Artifact: space_amplitude.npz (params, schema_version)

6.4 Co-moving Lyapunov exponent and information flow #

Kaneko (1986b) introduced the co-moving Lyapunov exponent $\lambda (v)$ to quantify how perturbations propagate through a CML. Consider a localised perturbation $\delta x_n(i)$ evolving under the linearised CML dynamics:

When $\lambda (v)>0$, a perturbation travelling at velocity $v$ grows exponentially; the propagation speed $v_p$ is the velocity at which $\lambda (v)$ crosses zero. Kaneko (1992c) showed that $v_p$ is finite for spatiotemporally chaotic CMLs, so disturbances spread at a well-defined speed, a discrete-map precursor of the convective/absolute instability distinction in spatiotemporal chaos.

Figure 27 shows $\lambda (v)$ for three CML regimes: pattern selection ($a=1.70$), pattern competition intermittency ($a=1.85$), and fully developed turbulence ($a=1.95$). The zero crossings give propagation speeds $v_p\approx 0.44$–$0.45$ for $a=1.85$ and $v_p\approx 0.47$–$0.48$ for $a=1.95$; the range is the change when the averaging time is doubled from $T=1000$ to $2000$. For $a=1.70$ the curve is positive near $v=0$ ($\lambda (0)\approx 0.07$) and crosses zero between $v\approx 0.15$ and $0.20$, but the crossing moves with $T$ and is not converged.

sec08_sti/comoving_lyapunovcodeinteractcopy linkenlarge

How this figure was computed

# Reduced preview (windowed estimator, N = 400, T = 150, 4 velocities).
# Full figure: dynachaos run sec08_sti
# Windowed co-moving Lyapunov exponent lambda(v) for a logistic coupled map lattice (eps = 0.3):
# growth of the perturbation norm over a window of 2 w + 1 sites that moves with velocity v.
# Static figure: N = 3000, T = 1000, w = 15, 51 velocities on [0, 1]; zero crossings give v_p.
import numpy as np

from dynachaos.diagnostics.comoving_lyapunov import (
    comoving_lyapunov_windowed_logistic,
    comoving_zero_crossings,
)

v_values = np.linspace(0.0, 1.0, 4)
for a in (1.70, 1.85, 1.95):
    lam_v = comoving_lyapunov_windowed_logistic(
        a=a, eps=0.3, N=400, T=150, v_values=v_values, n_transient=500, w=5, i0=50
    )
    crossings = comoving_zero_crossings(v_values, lam_v)
    print(f"a={a:.2f}  lambda(v)={np.round(lam_v, 3)}  zero crossings={np.round(crossings, 2)}")

copy

Artifact: comoving_lyapunov.npz (v_values, a_values, eps, N, T, w, i0, lambda_a1.70, lambda_a1.85, lambda_a1.95, zc_a1.70, zc_a1.85, zc_a1.95)

6.5 Thermodynamics of spatiotemporal chaos #

For a CML with $N$ sites, the full Lyapunov spectrum $\{{\lambda }_i{\}}_{i=1}^N$ grows with system size. Kaneko (1989d) introduced the concept of intensive thermodynamic quantities for spatiotemporal chaos: the KS entropy density $h_{\mathrm{\text{KS}}}=H_{\mathrm{\text{KS}}}/N$, the Lyapunov dimension density $d_{\mathrm{\text{KY}}}/N$, and the density of positive Lyapunov exponents. When the spatial correlation length $\xi$ is much smaller than the system size $N$, these densities converge to well-defined intensive values, and the CML behaves like a thermodynamic system with local degrees of freedom.

This convergence breaks down near phase transitions (where $\xi \to \infty$) and in the pattern selection phase (where long-range spatial order imposes correlations across the lattice). Figure 28 combines the spatial correlation function $C(r)$ with the Kolmogorov–Sinai entropy density $h_{\mathrm{\text{KS}}}(L)/L$ (dynachaos run sec08_sti). Rather than using an intensive ratio like ${\lambda }_{\max }/L$ (which decreases as $1/L$ by construction), extensivity is tested by computing the sum of all positive Lyapunov exponents of the $L$-site subsystem Jacobian divided by $L$. For $a=1.85$ and $1.95$, $h_{\mathrm{\text{KS}}}(L)/L$ stays within about 5 % of a constant over $L=10$–$100$ ($\approx 0.076$ and $\approx 0.095$), consistent with extensive chaos. For $a=1.50$ it is also nearly constant ($\approx 0.03$); for $a=1.70$ it falls from $0.007$ to $0.002$ over the same range, so extensivity is not seen there.

sec08_sti/correlation_decaycodeenlarge

How this figure was computed

# Reduced preview (N = 60, 4 a values). Full figure: dynachaos run sec08_sti
# Spatial autocorrelation C(r)/C(0) of a logistic coupled map lattice (eps = 0.3).
# 4 a values (1.50, 1.70, 1.85, 1.95); static figure adds panel (b), the KS entropy density
# h_KS(L)/L, the sum of the positive Lyapunov exponents of an L-site subsystem divided by L.
import numpy as np

from dynachaos.cml.correlation_figure import compute_ks_entropy_density
from dynachaos.cml.primitives import cml_step_logistic as cml_step

N = 60
eps = 0.3
rng = np.random.default_rng(42)
for a in (1.50, 1.70, 1.85, 1.95):
    x = rng.uniform(-0.5, 0.5, N)
    for _ in range(500):
        x = cml_step(x, a, eps)
    snapshots = np.empty((200, N))
    for s in range(200):
        x = cml_step(x, a, eps)
        snapshots[s] = x
    fluct = snapshots - snapshots.mean()
    fft_f = np.fft.rfft(fluct, axis=1)
    power = (fft_f * np.conj(fft_f)).real
    full_corr = np.fft.irfft(power, n=N, axis=1).mean(axis=0) / N
    c_norm = full_corr / full_corr[0]
    print(f"a={a:.2f}  C(r=0..3)/C(0)={np.round(c_norm[:4], 3)}")

# h_KS(L)/L on a small lattice (N = 40, L = 10 and 20, 400 steps): the chaotic states carry
# entropy per site, the pattern-selection state (a = 1.70) very little.
a_values, L_values, h_density = compute_ks_entropy_density(
    N=40, a_values=(1.70, 1.95), L_values=(10, 20), n_transient=300, n_iter=400
)
for a, row in zip(a_values, h_density, strict=True):
    print(f"a={a:.2f}  h_KS(L)/L for L={L_values.tolist()}: {np.round(row, 4)}")

# Full figure: dynachaos run sec08_sti

copy

Artifact: correlation_decay.npz (a_corr, r_vals, all_corr, xi_values, a_lyap, L_vals, h_density)

6.6 Diagnostic spotlight: RQA and spacetime recurrence #

Recurrence plots (Eckmann et al. 1987; Marwan et al. 2007) visualise the times at which a trajectory returns to a previously visited neighbourhood. From the binary recurrence matrix $R_{ij}=\mathrm{\Theta }(\delta -\parallel {\mathit{𝒙}}_i-{\mathit{𝒙}}_j\parallel )$ (where $\delta$ is the recurrence threshold), one extracts quantitative measures:

• RR (recurrence rate): fraction of recurrent points.

• DET (determinism): fraction in diagonal lines ($\ge l_{\mathrm{\text{min}}}$).

• LAM (laminarity): fraction in vertical lines ($\ge v_{\mathrm{\text{min}}}$).

• ENTR: Shannon entropy of the diagonal line length distribution.

RQA tracks recurrent states and temporal structures: the laminarity measure LAM tracks the fraction of near-recurrent states, while ENTR captures the structural complexity of diagonal recurrence lines. We evaluate these measures on the single-site delayed logistic map benchmark (Section 3.1) rather than across the full CML lattice, avoiding embedding-dimension ambiguities. Figure 29 shows these four measures versus $D$. The transition from torus to chaos manifests as a jump in ENTR and a restructuring of DET and LAM.

sec11_diagnostics/rqa_measurescodeinteractcopy linkenlarge

How this figure was computed

# Reduced preview (3 D values). Full figure: dynachaos run sec11_diagnostics
# RQA measures (RR, DET, LAM) along the delayed logistic map (alpha=0.3).
import numpy as np

from dynachaos.diagnostics.recurrence import recurrence_matrix, rqa
from dynachaos.maps.primitives import delayed_logistic

alpha = 0.3
for D in (1.6, 1.9, 2.1):
    state = np.array([0.51, 0.49])
    for _ in range(2000):
        state = delayed_logistic(state, alpha, D)
    traj = np.empty((500, 2))
    for i in range(500):
        state = delayed_logistic(state, alpha, D)
        traj[i] = state
    R, _ = recurrence_matrix(traj, percentile=5)
    stats = rqa(R, l_min=2, v_min=2)
    print(f"D={D:.2f}  RR={stats['RR']:.3f}  DET={stats['DET']:.3f}  LAM={stats['LAM']:.3f}")

# Full figure: dynachaos run sec11_diagnostics

copy

Artifact: rqa_measures.npz (D, RR, DET, LAM, ENTR)

6.7 Diagnostic spotlight: the temporal-intermittency taxonomy #

The laminar/burst dichotomy that organises spatiotemporal intermittency has a classical low-dimensional counterpart. Pomeau and Manneville (1980) classified the intermittent routes to chaos by the way a periodic orbit loses stability: Type I through a tangent (saddle–node) bifurcation, Type II through a subcritical Hopf bifurcation, and Type III through a subharmonic (inverse period-doubling) bifurcation. A distinct mechanism, on–off intermittency, arises when an invariant subspace loses transverse stability at a blowout bifurcation (Platt et al. 1993; Heagy et al. 1994). We apply laminar-phase detection and heavy-tail estimators (maximum-likelihood exponent with bootstrap confidence interval and a goodness-of-fit $p$-value) to canonical realisations of each route, as a check that the diagnostics behave correctly where the underlying scaling laws are known.

Type I (Figure 30) is realised by the logistic map just below its period-3 window, where three tangent channels open between the third iterate and the diagonal. The laminar-length distribution follows a power law ($\alpha \approx -1.5$, goodness-of-fit $p=0.45$) and the mean laminar length obeys the normal-form scaling $⟨\ell ⟩\propto {\epsilon }^{-1/2}$; the Lorenz system at $\rho =166.2$ provides a continuous-time exemplar of the same near-diagonal reinjection channel.

Type III (Figure 31) is illustrated through the flip channel of a second-iterate return map, with the reinjection probability density $M(x)$ characterised in the formulation of Elaskar et al. (2011). A power-law fit to the laminar distribution gives goodness-of-fit $p=0.03$, which rejects a pure power law at the 5% significance level.

Type II (Figure 32) is shown only as a normal-form demonstration: the spiral growth and reinjection of the subcritical-Hopf normal form yield a laminar-length distribution with an exponential (not power-law) envelope, and a power-law fit is rejected (goodness-of-fit $p=0$). No clean physical Type II example is included (the p–n diode and the forced jet are the usual citations), so no claim is made beyond the normal form.

On–off intermittency (Figure 33) uses a skew logistic system tuned near its blowout onset, where the transverse Lyapunov exponent ${\lambda }_{⟂}$ crosses zero. The off-time (laminar) distribution has a fitted exponent $\alpha =-1.66$ (95% CI $[-2.62,-1.40]$), with goodness-of-fit $p=0.02$, which also rejects a pure power law at the 5% level. The mean off-time follows the expected scaling $⟨\tau ⟩\propto |{\lambda }_{⟂}{|}^{-1}$. The burst-amplitude distribution is not a power law (goodness-of-fit $p=0$); we report the rejection and fit nothing else. These temporal routes give the spatiotemporal case a vocabulary, not a decomposition. Each laminar interval in the coupled-map-lattice intermittency of Section 6.2 is, locally in time, a site lingering near a marginally stable state before a burst, and the classical types list the ways such local stability can be lost. No single type accounts for the full lattice dynamics, which couple these local events across space.

sec12_intermittency/type_i_intermittencycodeenlarge

How this figure was computed

# Type-I intermittency: tangent bifurcation and laminar-phase statistics.
# Full figure: dynachaos run sec12_intermittency
from dynachaos.diagnostics.intermittency import detect_laminar_phases
from dynachaos.maps.intermittency import LOGISTIC_TYPE_I_ONSET, logistic_type_i_oracle

series = logistic_type_i_oracle(20_000, x0=0.2, r=LOGISTIC_TYPE_I_ONSET - 1e-4)
_, lengths = detect_laminar_phases(series, method="period", period=3, percentile=70.0)
print(
    f"r_c={LOGISTIC_TYPE_I_ONSET:.4f}  n_laminar_phases={len(lengths)}  "
    f"mean_length={lengths.mean():.2f}"
)

copy

Artifact: type_i_intermittency.npz (schema_version, source_file, seed, logistic_mechanism_r, logistic_tail_r, logistic_laminar_lengths, type_i_tail_alpha, type_i_vuong_z, normal_form_beta, normal_form_eps, normal_form_mean_lengths, logistic_f3_return_points, logistic_f3_channel_slope, lorenz_return_points, lorenz_channel_slope, type_i_tail_alpha_ci, type_i_tail_gof_p)

sec12_intermittency/type_iii_intermittencycodeenlarge

How this figure was computed

# Type-III intermittency (flip bifurcation): escape-episode statistics.
# Full figure: dynachaos run sec12_intermittency
import numpy as np

from dynachaos.maps.intermittency import pm_type_iii_oracle

eps, a, escape_threshold = 2e-3, 1.0, 0.35
for x0 in (2e-6, 2e-4, 2e-3):
    orbit = pm_type_iii_oracle(8000, x0=x0, eps=eps, a=a)
    escaped = np.abs(orbit) >= escape_threshold
    if escaped.any():
        status = f"{int(np.argmax(escaped))} steps before escape"
    else:
        status = f">={len(orbit)} (no escape)"
    print(f"x0={x0:.0e}  laminar_length={status}")

copy

Artifact: type_iii_intermittency.npz (schema_version, source_file, seed, eps, a, escape_threshold, return_grid, f2_return_points, f2_linear_slope, f2_cubic_coefficient, reinjection_points, series, laminar_mask, laminar_lengths, laminar_tail_alpha, laminar_tail_alpha_ci, laminar_tail_gof_p, rpd_thresholds, rpd_conditional_means, rpd_slope, rpd_intercept, rpd_alpha, rpd_rvalue)

sec12_intermittency/type_ii_intermittencycodeenlarge

How this figure was computed

# Type-II intermittency (subcritical Hopf): normal-form escape-time computation.
# Full figure: dynachaos run sec12_intermittency
import warnings

import numpy as np

from dynachaos.maps.intermittency import pm_type_ii_oracle

eps, a, theta, escape_threshold = 2e-3, 1.0, 0.17, 0.35
with warnings.catch_warnings():
    warnings.filterwarnings("ignore", category=RuntimeWarning)
    for r0 in (1e-3, 1e-2, 2e-2):
        orbit = pm_type_ii_oracle(2500, x0=r0, y0=0.0, eps=eps, a=a, theta=theta)
        radius = np.linalg.norm(orbit, axis=1)
        escaped = radius >= escape_threshold
        if escaped.any():
            status = f"{int(np.argmax(escaped))} steps before escape"
        else:
            status = f">={len(orbit)} (no escape)"
        print(f"r0={r0:.0e}  laminar_length={status}")

copy

Artifact: type_ii_intermittency.npz (schema_version, source_file, seed, eps, a, theta, escape_threshold, spiral_orbit, spiral_radius, spiral_escape_index, reinjection_radii, laminar_lengths, laminar_histogram_edges, laminar_histogram_density, laminar_tail_alpha, laminar_tail_alpha_ci, laminar_tail_gof_p, exponential_rate, exponential_intercept, exponential_rvalue)

sec12_intermittency/on_off_intermittencycodeenlarge

How this figure was computed

# On-off intermittency across blowout onset (lambda_perp crossing zero).
# Full figure: dynachaos run sec12_intermittency
import numpy as np

from dynachaos.maps.intermittency import on_off_skew_logistic_oracle


def mean_run_length(mask):
    lengths, run = [], 0
    for flag in mask:
        if flag:
            run += 1
        elif run > 0:
            lengths.append(run)
            run = 0
    if run > 0:
        lengths.append(run)
    return np.mean(lengths) if lengths else 0.0


for eps in (0.48, 0.50, 0.51):
    skew = on_off_skew_logistic_oracle(20_000, x0=0.217, y0=1e-2, eps=eps)
    y = np.abs(skew[1000:, 1])
    off_mask = y <= np.percentile(y, 90)
    lambda_perp = np.log(2.0 * eps)
    print(
        f"eps={eps:.3f}  lambda_perp={lambda_perp:+.4f}  "
        f"mean_off_run_length={mean_run_length(off_mask):.1f}"
    )

copy

Artifact: on_off_intermittency.npz (schema_version, source_file, seed, benchmark_eps, benchmark_lambda_perp, benchmark_series, benchmark_laminar_mask, benchmark_laminar_lengths, benchmark_burst_lengths, benchmark_burst_amplitudes, benchmark_threshold_percentile, off_time_alpha, off_time_alpha_ci, off_time_gof_p, burst_amplitude_alpha, burst_amplitude_alpha_ci, burst_amplitude_gof_p, scaling_eps_values, lambda_abs_values, mean_off_lengths, mean_off_beta, skew_driver_series, skew_transverse_series)

7.1 Violation of the law of large numbers #

Kaneko’s claim is that collective chaos violates the law of large numbers but not the central limit theorem: $P(h)$ approaches Gaussian form, yet its variance stays macroscopic (Kaneko 1990b).

For $N$ independent, identically distributed random variables, the variance of the mean decreases as $1/N$. Kaneko showed that for the GCM at certain parameter values, the mean-square deviation (MSD)

7.2 Collective chaos and cluster states #

The nontrivial mean-field dynamics indicates that the $N$ oscillators, although not synchronised, maintain correlations that resist averaging. Kaneko showed that the system self-organises into cluster states: groups of oscillators that share the same trajectory, with the number and sizes of clusters depending on $a$ and $\epsilon$. These clusters are not imposed by the coupling topology (which is all-to-all and thus featureless) but emerge from the dynamics, a form of dynamical symmetry breaking.

sec10_gcm/gcm_msdcodeenlarge

How this figure was computed

# Reduced preview (N in 50..5000, a = 1.99, single seed). Full figure: dynachaos run sec10_gcm
# Mean-square deviation (variance) of the mean field h_n of a globally coupled logistic map vs N.
# Full figure averages 8 seeds across 5 a values with saturation fit MSD(N) = A/N + B (B > 0).
from dynachaos.cml.globally_coupled import compute_msd_sweep

result = compute_msd_sweep(
    a_values=(1.99,),
    N_values=(50, 500, 5000),
    n_seeds=1,
    n_transient=300,
    n_sample=400,
    output_path=None,
)
for N, msd in zip(result["N_values"], result["msd_mean"][0], strict=True):
    print(f"N={N:5d}  Var(h)={msd:.3e}")

copy

Artifact: gcm_msd.npz (a_values, N_values, eps, n_seeds, msd_all, msd_mean, msd_std, msd_stderr, fits_A, fits_B, fits_B_stderr, fits_r2)

sec10_gcm/gcm_distributioncodeenlarge

How this figure was computed

# Reduced preview (N in (50, 5000), a=1.99). Full figure: dynachaos run sec10_gcm
# Distribution P(h) of the mean field h_n of a globally coupled logistic map.
import numpy as np

from dynachaos.cml.primitives import gcm_step
from dynachaos.maps.primitives import logistic

a = 1.99
eps = 0.1
rng = np.random.default_rng(42)
for N in (50, 5000):
    x = rng.uniform(-1, 1, N)
    for _ in range(300):
        x = gcm_step(x, a, eps)
    h_series = np.empty(400)
    for t in range(400):
        x = gcm_step(x, a, eps)
        h_series[t] = np.mean(logistic(x, a))
    print(f"N={N:5d}  mean(h)={h_series.mean():.4f}  std(h)={h_series.std():.4f}")

# Full figure: dynachaos run sec10_gcm

copy

Artifact: gcm_results.npz (N_values, a, eps)

7.3 Clustered attractors and switching #

In the GCM, a cluster state is a partition of the $N$ sites into $k$ groups with nearly identical instantaneous values at a given iterate. The partition is characterised by the cluster multiplicities $(N_1,N_2,\dots ,N_k)$ with ${\sum }_l{N}_l=N$. Kaneko (1990a) showed that at moderate coupling ($\epsilon \approx 0.1$), the GCM spontaneously selects a small number of clusters ($k=2$–$5$), and that the partition coding provides a discrete “label” for each macroscopic state.

In two-cluster attractors, coded by the populations $N_{+}$ and $N_{-}$ of the two clusters, Kaneko found posi-nega switching (Kaneko 1989a): an input on a single site moves it from one cluster to the other, so the system switches among the coded attractors. When $N_{+}$ is pushed beyond a threshold $N_{\mathrm{t}\mathrm{h}\mathrm{r}}$, the system passes through a chaotic transient and every $+$ site changes to $-$ and vice versa. Figure 36 visualises the cluster structure as a spacetime diagram after sorting the site values at each time step, so that coherent clusters appear as broad plateaus whose widths encode their populations.

sec10_gcm/gcm_clusterscodeenlarge

How this figure was computed

# Cluster count of a globally coupled logistic map (a=1.55, N=100) at eps=0.1 and eps=0.3.
# Full figure: dynachaos run sec10_gcm
from dynachaos.cml.gcm_clusters import compute_clusters

for eps in (0.1, 0.3):
    result = compute_clusters(
        a=1.55,
        eps=eps,
        n_sites=100,
        n_transient=2000,
        n_record=10,
        output_path=None,
    )
    n_clusters = len(set(result["cluster_labels"][-1]))
    print(f"eps={eps:.2f}  clusters at final step={n_clusters}")

copy

Artifact: gcm_clusters.npz (cluster_labels, x_record, a, eps, N)

7.4 Collective degrees of freedom #

Shibata and Kaneko (1998) introduced a collective Lyapunov exponent ${\lambda }_C$ for the mean field $h_n=N^{-1}{\sum }_{j}f(x_n(j))$ of a GCM whose elements have a distributed nonlinearity. A perturbation of $h_n$ observed at scales above a microscopic scale ${\mathrm{\Delta }}_m$ grows at the rate ${\lambda }_C$, which is much smaller than the largest microscopic exponent ${\lambda }_m$; when ${\mathrm{\Delta }}_m$ decreases as $N^{-1/2}$, the mean field has its own low-dimensional chaotic dynamics in the limit $N\to \infty$.

The partially ordered (PO) phase, between full synchronisation and full turbulence, is characterised by clusters that form, merge, and split in a marginally stable regime. Kaneko (1994a) showed that in this regime, information cascades through the cluster hierarchy with marginal stability: the system sits at the edge of a bifurcation, generating long-range temporal correlations. Takeuchi, Ginelli, and Chaté (Takeuchi et al. 2009) later identified collective Lyapunov modes that capture macroscopic fluctuations in large coupled chaotic systems.

Figure 37 shows the largest Lyapunov exponent ${\lambda }_{\mathrm{\text{max}}}$ of the GCM (16) with identical elements at $\epsilon =0.1$, computed from the tangent map (dynachaos run sec10_gcm). It becomes positive, by more than two standard errors over eight seeds, from $a\approx 1.64$ and reaches $0.45$ at $a=2$, with no dependence on $N$ between 500 and 2000. This is the exponent of the microscopic chaos. An infinitesimal perturbation of the mean field grows at the same rate, so ${\lambda }_{\mathrm{\text{max}}}>0$ alone says nothing about collective chaos (Shibata and Kaneko 1998).

Shibata and Kaneko therefore keep the perturbation finite and measure the finite-size Lyapunov exponent (Aurell et al. 1996) of the mean field, ${\lambda }_{{\delta }_0}(\mathrm{\Delta })=\ln (\mathrm{\Delta }/{\delta }_0)/⟨\tau (\mathrm{\Delta })⟩$, where $\tau (\mathrm{\Delta })$ is the time for $|h_n^{\prime }-h_n|$ to grow from ${\delta }_0$ to $\mathrm{\Delta }$. The normalised time $t(\mathrm{\Delta })=⟨\tau (\mathrm{\Delta })⟩+\ln {\delta }_0/{\lambda }_m$ grows as $\ln \mathrm{\Delta }/{\lambda }_m$ below ${\mathrm{\Delta }}_m$ and as $\ln (\mathrm{\Delta }/{\mathrm{\Delta }}_m)/{\lambda }_C+\ln {\mathrm{\Delta }}_m/{\lambda }_m$ above it (their Eq. 6). Figure 38 reproduces their chaotic case, $a(i)=1.9+0.05(2i-N)/(2N)$, $\epsilon =0.098$ (dynachaos run sec10_gcm). The perturbed copy starts at $x+{10}^{-3}{\delta }_0\sigma$ with $\sigma$ uniform in $[-1,1]$, and $\tau$ is counted from the first step at which $|h_n^{\prime }-h_n|\ge {\delta }_0$, so that the perturbation has aligned with the growing direction. Below $\mathrm{\Delta }\approx {10}^{-3}$ the exponent equals ${\lambda }_m=0.410$ (0.41 in (Shibata and Kaneko 1998)). Fits of Eq. 6 for $\mathrm{\Delta }\le 0.05$ give ${\lambda }_C=0.029$, $0.021$, $0.024$ and $0.022$ for $N={10}^4$, ${10}^5$, ${10}^6$ and ${10}^7$ (fit errors at most $0.0013$), against $0.02$ in (Shibata and Kaneko 1998), and ${\mathrm{\Delta }}_m$ decreases as $N^{-0.46\pm 0.03}$, with ${\mathrm{\Delta }}_m=7.7\times {10}^{-4}$ at $N={10}^6$ ($8.6\times {10}^{-4}$ in (Shibata and Kaneko 1998)). The values depend on the upper end of the fit range: from $\mathrm{\Delta }=0.01$ to $0.1$, ${\lambda }_C$ at $N\ge {10}^6$ lies between 0.020 and 0.026, and the exponent of ${\mathrm{\Delta }}_m(N)$ between $-0.33$ and $-0.48$. For identical elements ($a=1.6962$, $\epsilon =0.008$) the curves $t(\mathrm{\Delta })$ for $N={10}^5$ and ${10}^6$ do not separate, as Shibata and Kaneko found: there the scale of the mean-field motion does not separate from the microscopic one.

sec10_gcm/gcm_lyapunov_maxcodeenlarge

How this figure was computed

# Reduced preview (N = 200, 3 a values, 1 seed). Full figure: dynachaos run sec10_gcm
# Largest Lyapunov exponent lambda_max of the identical-element GCM (eps = 0.1) by the tangent map:
# v -> (1 - eps) f'(x) v + eps <f'(x) v>, renormalised every step. This is the microscopic
# (site) exponent, not the collective exponent of the mean field.
# Full figure: 8 seeds for N = 500, 1000, 2000 over a in [1.40, 2.00].
from dynachaos.cml.gcm_lyapunov_max import lyapunov_max

for a in (1.50, 1.80, 1.95):
    lam = lyapunov_max(a, 0.1, 200, n_transient=500, n_measure=2000, seed=42)
    print(f"a={a:.2f}  lambda_max={lam:+.4f}")

copy

Artifact: gcm_lyapunov_max.npz (a_values, N_values, eps, n_seeds, n_transient, n_measure, seed_base, runs_500, runs_1000, runs_2000, mean_500, mean_1000, mean_2000, se_500, se_1000, se_2000)

sec10_gcm/gcm_collective_fslecodeenlarge

How this figure was computed

# Reduced preview (N = 1000, 3 perturbations). Full figure: dynachaos run sec10_gcm
# Finite-size Lyapunov exponent (FSLE) of the mean field of the heterogeneous GCM
# a(i) = 1.9 + 0.05 (2i - N) / (2N), eps = 0.098 (Shibata & Kaneko 1998, Figs. 3-5):
# lambda(Delta) = log(Delta / delta0) / <tau(Delta)>, tau = first time |h' - h| reaches Delta.
# Full figure: N = 1e4, 1e5, 1e6, 100 perturbations; Eq. (6) fits give lambda_C, Delta_m.
# Larger sizes: python -m dynachaos.cml.gcm_collective_fsle --add 10000000
import numpy as np

from dynachaos.cml.gcm_collective_fsle import (
    CHAOS,
    DELTAS,
    a_profile,
    first_passage,
    lambda_m,
    step,
)

N = 1000
delta0 = 1e-7
a = a_profile(CHAOS["a0"], CHAOS["delta_a"], N)
rng = np.random.default_rng(N)
x = rng.uniform(-1.0, 1.0, N)
for _ in range(2000):
    x, _ = step(x, a, CHAOS["eps"])

taus = []
for _ in range(3):
    tau, x = first_passage(x, a, CHAOS["eps"], delta0, rng, 3000)
    taus.append(tau)
mean_tau = np.mean(taus, axis=0)

lam_m = lambda_m(a, CHAOS["eps"], x, n_transient=500, n_measure=2000)
print(f"lambda_m = {lam_m:.3f} (tangent map)")
for D in (1e-5, 1e-4, 1e-3, 1e-2):
    k = int(np.argmin(np.abs(np.log(DELTAS / D))))
    if mean_tau[k] <= 3000:
        print(
            f"Delta={DELTAS[k]:.1e}  lambda(Delta)={np.log(DELTAS[k] / delta0) / mean_tau[k]:.3f}"
        )

copy

Artifact: gcm_collective_fsle.npz (deltas, Ns, n_max, n_samples, chaos_a0, chaos_delta_a, chaos_eps, control_a0, control_delta_a, control_eps, chaos_tau_N10000, chaos_lam_m_N10000, chaos_tau_N100000, chaos_lam_m_N100000, chaos_tau_N1000000, chaos_lam_m_N1000000, chaos_tau_N1000000_d1e-06, chaos_tau_N1000000_d1e-05, control_tau_N100000, control_lam_m_N100000, control_tau_N1000000, control_lam_m_N1000000)

8.1 Mechanistic interpretation #

The mechanism-based ordering also explains why the same diagnostic cannot carry the whole paper. The appropriate observable is fixed by the invariant object being destabilised. In the phase-locking sections, that object is an invariant circle with a well-defined winding rate, so the rotation number is the natural observable. For torus oscillation and doubling, the change is not only in frequency content but in the stability of transverse directions, which makes the Lyapunov spectrum decisive. The same logic fixes the remaining pairings: basin and order-parameter views for symmetry breaking, where exchange-related attractors coexist; spatial diagnostics for CMLs, where perturbations propagate through laminar and burst regions; and finite-size and mean-field statistics for GCMs, where persistent correlations among nominally identical elements survive averaging.

Higher resolution makes these distinctions visible. Narrow windows resolve locking transitions organised by rational approximants, long trajectories reduce misclassification of weak chaos and nearly neutral tori, and finite-size sweeps test whether collective chaos persists as mean-field averaging becomes more effective. A diagnostic is decisive only when the degree of freedom that goes unstable is one it can see.

8.2 Toric chaos #

Yamagishi and Kaneko (2020) showed that chaos can exist on a high-dimensional torus and not only through its collapse. In systems of globally coupled phase oscillators, the attractor can be a high-dimensional torus ${\mathbb{𝕋}}^k$ (with $k\gg 1$) on which the dynamics is chaotic, which they call “toric chaos.” This adds a fate to Kaneko’s list for ${\mathbb{𝕋}}^3$ and beyond: a higher-dimensional torus can carry chaos without collapsing.

The coupled delayed logistic map (Section 5.1) supports three-torus motion that collapses to chaos as $D_B$ increases at fixed coupling (Figure 17). At the transition, two Lyapunov exponents are near zero while one becomes positive. This is the classical torus collapse, not toric chaos in the sense of Yamagishi and Kaneko (2020), which requires the torus to persist as a manifold even while the on-torus dynamics is chaotic.

8.3 Biological relevance #

Kaneko extended the dynamical clustering framework to biological networks. Kaneko (1994b) showed that the partition structure of GCM cluster states provides a natural encoding of cell-type diversity. In coupled maps with growth and death dynamics, clusters of oscillators that share a common time trace can then be reinterpreted as cell types, with cluster formation playing the role of cell differentiation (Kaneko 1997a). Ecological applications of coupled logistic maps appeared early (Satoh and Aihara 1990; Satoh 1991). In a recent preprint, Shimada and Kaneko (2026) continue this direction by connecting high-dimensional chaos and intermittent switching to coexistence in cross-feeding microbial-community models. In these models the macroscopic state is carried by partitions that form dynamically, not by fixed species labels or imposed network modules, which is why the mean-field and clustering mechanisms are still in use.

8.4 Diagnostic synergies #

The four modern diagnostics distributed across the preceding sections provide complementary views of Kaneko’s transitions:

• The 0-1 test (Section 4) is the most practical for binary classification (chaos vs. regular) from scalar data because it does not require equations or tangent maps.

• SALI (Section 3) requires the linearised flow; in the dissipative maps examined here, phase-space contraction aligns deviation vectors, so SALI stays bounded on a torus and decays on partial locking and on chaos, and its decay alone does not separate those two (cf. Moges et al. (Moges et al. 2025)); the decay rate ${\lambda }_2-{\lambda }_1$ does.

• Permutation entropy (Section 5) works from the ordinal patterns of a scalar observable alone, tracking regular windows and chaotic bands without a Jacobian.

• RQA (Section 6) responds to laminar residence times and repeated local patterns; we test it here on the single-site torus-to-chaos transition.

These diagnostics are not interchangeable. Each responds to a different way in which order is lost, and the limit of SALI on the dissipative maps above is one reason a single diagnostic cannot carry the paper.

8.5 Computational limitations: then and now #

Kaneko’s original parameter scans typically used $\sim 100$ points with $\sim 50,000$ iterations per point. Our scans use from 500 (Figure 17) to 200 000 (Figure 3) parameter values (Table B1). This increase matters because several transitions are organised by narrow windows, thin strands, or fractal boundaries; at high resolution, the tongue boundaries, period-adding plateaus, and basin stripes can be read directly off the plot.

8.6 What persisted after the 1990s #

Three later lines of work extended the framework beyond spatiotemporal chaos. First, the collective chaos discovered in the GCM led to a systematic study of Milnor attractor networks (Kaneko 1997b, 1998): in systems with many coexisting attractors, noise-induced transitions are governed by the Milnor attractor structure and not by basin volume. Second, chaotic itinerancy, a slow, directed wandering among partially stable attractor ruins proposed by Tsuda (1991) and developed with Kaneko (1992b, 1995; Kaneko and Tsuda 2001), became an organising idea for high-dimensional chaos. Third, dynamic clustering was reinterpreted as a model of cell differentiation: coupled maps with growth and death dynamics produce spontaneous diversification of cell types (Kaneko 1997a). Separately, renormalisation-group analysis clarified the scaling structure of torus doubling and strange nonchaotic attractors (Kuznetsov et al. 1995; Kuznetsov et al. 1998; Nishikawa and Kaneko 1996).