Distinguishing Chaos from Periodicity in Fusion Plasma Fluctuations: A Volterra-Kernel Approach with BIC Model Selection

Distinguishing Chaos from Periodicity in Fusion Plasma Fluctuations: A Volterra-Kernel Approach with BIC Model Selection Christopher K. Merrill With computational collaboration by Claude (Anthropic) and ChatGPT (OpenAI) December 24, 2025 Abstract Distinguishing chaotic from periodic dynamics in plasma fluctuation measurements is essential for selecting appropriate modeling and control strategies. We present a method based on Bayesian Information Criterion (BIC) model selection between linear and quadratic polynomial autoregressive models. The discriminant statistic ∆BIC quantifies whether nonlinear terms improve prediction sufficiently to justify their complexity cost. Crucially, this criterion detects chaos— dynamics requiring nonlinear modeling—rather than nonlinearity in the governing equations: limit cycle oscillations from nonlinear predator-prey systems yield ∆BIC <0 (linear sufficient), while deterministic chaos yields ∆BIC ≫0 (nonlinear required). We validate on canonical systems including the logistic map (∆BIC = +87,548), H´enon map (+130,589), and predator-prey oscillators (−6,867), demonstrating correct classification and robustness to measurement noise up to 30%. The method detects dynamical transitions, showing elevated ∆BIC near Hopf bifurcations relevant to confinement transitions. Application to Hasegawa-Wakatani turbulence simulations reveals regime-dependent behavior: nonlinear dynamics in the hydrodynamic limit, linear in the adiabatic limit. Complete algorithmic specification enables independent implementation. The method provides fusion plasma diagnostics with a simple, interpretable criterion answering: “Does this signal require nonlinear modeling for prediction?” Keywords: plasma turbulence, chaos detection, zonal flows, model selection, nonlinear dynamics, Bayesian Information Criterion 1 Introduction Turbulent transport remains a central challenge in magnetic confinement fusion, determining the energy and particle losses that limit plasma performance [Horton, 1999, Diamond et al., 2005]. The interaction between drift-wave turbulence and self-generated zonal flows plays a 1 crucial regulatory role: zonal flows, driven by Reynolds stress from the turbulence itself, act to shear apart turbulent eddies and suppress transport [Diamond et al., 2005, Fujisawa, 2009]. This feedback system has been described using predator-prey models, where turbulence (prey) drives zonal flow growth (predator), which in turn suppresses the turbulence [Kim and Diamond, 2003, Malkov and Diamond, 2009]. Experimental observations in tokamaks have revealed limit cycle oscillations (LCOs) consistent with predator-prey dynamics, particularly near the low-to-high confinement (L-H) transition [Schmitz et al., 2012, Conway et al., 2011, Estrada et al., 2011]. These oscillations, occurring at frequencies of 1–10 kHz, show characteristic phase relationships where zonal flow shear lags turbulence intensity by approximately 90 ° . Understanding and predicting these dynamics is essential for controlling access to high-confinement regimes in future devices including ITER. A fundamental question arises when analyzing fluctuation time series from plasma diagnostics: are the observed dynamics chaotic or periodic? This distinction has practical consequences. Chaotic dynamics, characterized by sensitive dependence on initial conditions, are fundamentally unpredictable beyond short horizons and require nonlinear modeling approaches. Periodic or quasi-periodic dynamics, even when arising from nonlinear mechanisms, are predictable using simpler linear methods. Misidentifying the nature of the dynamics leads to either unnecessary model complexity or inadequate predictive capability. Existing methods for characterizing nonlinear dynamics include Lyapunov exponent estimation [Wolf et al., 1985], surrogate data hypothesis testing [Theiler et al., 1992], recurrence quantification analysis [Marwan et al., 2007], and correlation dimension calculations [Grassberger and Procaccia, 1983]. While powerful, these techniques often require long, clean time series and careful parameter tuning. They also typically answer the question “is this system nonlinear?” rather than the more operationally relevant question “does this system require nonlinear modeling for prediction?” We present a complementary approach based on model selection between linear and quadratic polynomial autoregressive models. Using the Bayesian Information Criterion (BIC) to penalize model complexity, we identify when quadratic (nonlinear) terms provide sufficient predictive improvement to justify their inclusion. The key insight is that this criterion detects chaos—dynamics where nonlinear terms are essential for prediction—rather than nonlinearity in the governing equations per se. This distinction is crucial for fusion plasma applications. Predator-prey limit cycles, despite arising from nonlinear ODEs with terms like NU and N2, produce quasi-periodic time series that linear autoregressive models predict well. The BIC correctly identifies these as “linear” because adding quadratic terms does not improve prediction enough to offset the complexity penalty. Conversely, turbulent cascade dynamics with broadband fluctuations and intermittent bursts show strong preference for quadratic models. The contributions of this paper are: 1. A simple, interpretable criterion (∆BIC) for distinguishing chaotic from periodic dynamics in scalar time series 2. Complete algorithmic specification enabling independent implementation 2 3. Validation on canonical nonlinear systems including maps, flows, and predator-prey oscillators 4. Demonstration of bifurcation detection near dynamical transitions 5. Application to Hasegawa-Wakatani plasma turbulence simulations 6. Explicit characterization of limitations and failure modes The remainder of this paper is organized as follows. Section 2 presents the mathematical framework, BIC model selection criterion, noise titration protocol, and complete algorithm pseudocode. Section 3 validates the method on canonical systems and applies it to plasma turbulence. Section 4 interprets the results, discusses limitations, and suggests applications. Section 5 summarizes our findings. A reference implementation exists and is under active development; we welcome methodological correspondence. 2 Methods We present a model selection approach to distinguish chaotic dynamics from periodic or stochastic behavior in scalar time series. The method compares linear and quadratic autoregressive models using the Bayesian Information Criterion (BIC), where preference for quadratic terms indicates deterministic nonlinear structure characteristic of chaos. 2.1 Polynomial Autoregressive Framework Consider a scalar time series {yt}N t=1 sampled at uniform intervals ∆t. We model the dynamics using delay-coordinate embedding with maximum lag L, constructing feature vectors: xt= (yt−1, yt−2, . . . , yt−L)⊤∈RL(1) The linear autoregressive model of order L(AR(L)) predicts: ˆy(1) t=β0+ L X j=1 βjyt−j(2) where superscript (1) denotes first-order (linear) terms. This model has k1=L+1 parameters including the intercept. The quadratic autoregressive model augments Eq. (2) with second-order polynomial terms: ˆy(2) t=β0+ L X j=1 βjyt−j+ L X j=1 L X k=j γjkyt−jyt−k(3) The quadratic model has k2= 1 + L+L(L+ 1)/2 parameters. For L= 20, this yields k1= 21 and k2= 231 parameters. 3 Both models are fit using ridge regression to handle potential multicollinearity: ˆ θ= arg min θ(N X t=L+1 (yt−ˆyt)2+λ∥θ∥2 2)(4) where λ > 0 is the regularization parameter. We use λ= 10−3throughout, which provides numerical stability without significantly biasing coefficient estimates for well-conditioned problems. 2.2 BIC Model Selection The Bayesian Information Criterion balances goodness-of-fit against model complexity: BIC = nln RSS n+kln(n) (5) where n=N−Lis the number of predictions, RSS = PN t=L+1(yt−ˆyt)2is the residual sum of squares, and kis the number of model parameters. Lower BIC indicates a preferred model. We define the discriminant statistic: ∆BIC = BIClinear −BICquadratic (6) The interpretation is:  ∆BIC >0: Quadratic model preferred despite complexity penalty →nonlinear/chaotic  ∆BIC <0: Linear model preferred →periodic or stochastic  |∆BIC|<2: Inconclusive (models approximately equivalent) This criterion is conservative: the quadratic model must overcome a substantial complexity penalty (additional ∼L2/2 parameters) to be selected. In practice, chaotic systems yield ∆BIC ∼104–105, while periodic systems yield ∆BIC ∼ −102–−103. Critical insight: A positive ∆BIC indicates that quadratic terms provide predictive value beyond what linear terms capture. For deterministic systems, this occurs when the dynamics exhibit sensitive dependence on initial conditions (chaos), not merely when the governing equations are nonlinear. Limit cycles from nonlinear ODEs produce quasi-periodic time series that linear AR models predict well. 2.3 Noise Titration Protocol To assess robustness and quantify the “degree of nonlinearity,” we employ noise titration [Poon and Barahona, 2001]. The protocol adds calibrated Gaussian noise to the time series and measures how ∆BIC degrades: 1. Compute baseline ∆BIC0on original series {yt} 2. For noise levels σadd ∈ {0.05,0.10,0.15,...,0.50} × std(y): 4 (a) Generate corrupted series: ˜yt=yt+ϵt, where ϵt∼ N(0, σ2 add) (b) Compute ∆BIC(σadd) (c) Repeat M= 20 times with different noise realizations (d) Record mean and standard deviation 3. Identify critical noise level σ∗where ∆BIC crosses zero The critical noise level σ∗quantifies predictability: systems with larger σ∗have stronger deterministic structure that persists despite measurement noise. For chaotic systems, we typically observe σ∗>0.3×std(y); for periodic systems, ∆BIC remains negative for all noise levels. 2.4 Algorithm Summary Algorithm 1 provides complete pseudocode for the chaos detection procedure. 5 Algorithm 1 Volterra-Kernel Chaos Detection Require: Time series y= (y1, . . . , yN), max lag L, regularization λ Ensure: Classification ∈ {chaotic,periodic/stochastic}, ∆BIC 1: Preprocessing: 2: Remove mean: yt←yt−¯y 3: (Optional) Standardize: yt←yt/std(y) 4: Construct design matrices: 5: n←N−L ▷ Number of samples 6: for t=L+ 1 to Ndo 7: x(1) t←(1, yt−1, yt−2, . . . , yt−L)▷Linear features 8: x(2) t←(x(1) t, y2 t−1, yt−1yt−2, . . . , y2 t−L)▷Add quadratic 9: end for 10: X(1) ∈Rn×(L+1),X(2) ∈Rn×(1+L+L(L+1)/2) 11: ytarget ←(yL+1, . . . , yN)⊤ 12: Fit linear model (ridge regression): 13: ˆ β←(X(1)⊤X(1) +λI)−1X(1)⊤ytarget 14: RSS1← ∥ytarget −X(1) ˆ β∥2 2 15: BIC1←nln(RSS1/n)+(L+ 1) ln(n) 16: Fit quadratic model (ridge regression): 17: ˆ θ←(X(2)⊤X(2) +λI)−1X(2)⊤ytarget 18: RSS2← ∥ytarget −X(2)ˆ θ∥2 2 19: k2←1 + L+L(L+ 1)/2 20: BIC2←nln(RSS2/n) + k2ln(n) 21: Compute discriminant: 22: ∆BIC ←BIC1−BIC2 23: Classify: 24: if ∆BIC >2then 25: return chaotic, ∆BIC 26: else if ∆BIC <−2then 27: return periodic/stochastic, ∆BIC 28: else 29: return inconclusive, ∆BIC 30: end if 2.5 Parameter Selection Guidelines Maximum lag L: Should capture the characteristic correlation time of the dynamics. For oscillatory systems with period T, use L≳T/∆t. Typical values: L= 10–30. Excessively large Lincreases computational cost (O(L2) features) and may cause overfitting. Regularization λ: Values 10−4–10−2provide numerical stability without significant bias. Increase λfor noisy data or when n/k2<10. Minimum sample size: We recommend n > 10 ×k2for reliable BIC estimates. For L= 20, this requires N > 2500 samples. Stationarity: The method assumes stationary or quasi-stationary dynamics. For non6 stationary data, apply windowed analysis with segments of length ∼5–10 ×L. 2.6 Computational Complexity Feature construction requires O(nL2) operations. Ridge regression via the normal equations costs O(k3 2) = O(L6) for matrix inversion, though iterative methods reduce this for large L. In practice, analysis of N= 5000 samples with L= 20 completes in <1 second on standard hardware. 3 Validation Results We validate the method on canonical dynamical systems with known properties, then apply it to plasma turbulence simulations. All analyses use L= 15–20 lags and λ= 10−3unless otherwise noted. 3.1 Canonical Chaotic Systems Table 1 summarizes results for well-characterized chaotic maps and flows. Table 1: Chaos detection in canonical nonlinear systems. Positive ∆BIC correctly identifies chaotic dynamics; negative values indicate periodic or quasi-periodic behavior predictable by linear models. System Parameters N∆BIC Classification Chaotic systems Logistic map r= 3.9 4800 +87,548 Nonlinear Logistic map r= 3.7 4800 +52,341 Nonlinear H´enon map (x)a= 1.4, b = 0.3 4800 +130,589 Nonlinear H´enon map (y)a= 1.4, b = 0.3 4800 +98,205 Nonlinear Periodic/quasi-periodic systems Sine wave ω= 1 4000 −995 Linear Predator-prey LCO See text 4750 −6,867 Linear Predator-prey LCO + 15% noise 4750 −256 Linear Stochastic systems White noise N(0,1) 5000 −342 Linear AR(1) process ϕ= 0.8 5000 −287 Linear The logistic map xn+1 =rxn(1−xn) at r= 3.9 exhibits fully developed chaos with positive Lyapunov exponent λLyap ≈0.5. The method yields ∆BIC = +87,548, unambiguously indicating that quadratic terms dramatically improve prediction. At r= 3.7 (weaker chaos), ∆BIC remains strongly positive but reduced. The H´enon map (xn+1, yn+1) = (1 −ax2 n+yn, bxn) produces even larger ∆BIC values, reflecting the map’s explicit quadratic structure. Both components are correctly classified as chaotic. 7 For the predator-prey limit cycle oscillator (Section 3.2), ∆BIC is strongly negative despite the underlying nonlinear ODEs. This confirms the critical insight: the method detects chaos (unpredictability requiring nonlinear modeling), not merely nonlinearity in the governing equations. Quasi-periodic oscillations, regardless of their mechanistic origin, are well-predicted by linear autoregressive models. Pure stochastic processes (white noise, AR(1)) yield negative ∆BIC, correctly indicating no deterministic nonlinear structure. 104103102 0 102103104105 BIC (signed log scale) Logistic (r=3.9) Hénon map Predatorprey Sine wave White noise Chaotic (nonlinear required) Periodic (linear sufficient) (a) Classification by system type Figure 1: Classification by ∆BIC: chaotic systems (logistic, H´enon) yield strongly positive values indicating nonlinear terms are required for prediction, while periodic systems (predator-prey, sine) and stochastic processes yield negative values indicating linear models suffice. 3.2 Predator-Prey Dynamics To model turbulence-zonal flow interactions observed in fusion plasmas [Diamond et al., 2005, Kim and Diamond, 2003], we examine a predator-prey system: dN dt =γNN−αNU −σNN2(7) dU dt =−γUU+βNU −σUU2(8) where Nrepresents turbulence intensity (prey) and Urepresents zonal flow shear (predator). Parameters are chosen to produce stable limit cycle oscillations: γN= 5000, α= 8000, γU= 2000, β= 6000, σN= 1000, σU= 500 (units: s−1). Despite the nonlinear coupling terms NU,N2, and U2, the resulting time series are quasi-periodic with well-defined frequency f≈500–1000 Hz (Figure 2). Linear AR models capture this periodicity efficiently, yielding ∆BIC <0 for both N(t) and U(t) components. 3.3 Bifurcation Detection We test the method’s sensitivity to dynamical transitions using the Malkov-Diamond model [Malkov and Diamond, 2009], which exhibits a Hopf bifurcation as the heating power P 8 0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0 Time (ms) 0.2 0.4 0.6 0.8 Energy (a.u.) (a) Time series Turbulence N Flow U 0.34 0.36 0.38 0.40 0.42 Turbulence energy N 0.54 0.56 0.58 0.60 Flow energy U (b) Phase space (limit cycle) Logistic (chaotic) Hénon (chaotic) Pred-prey (periodic) 4 2 0 2 4 BIC (signed log) +87,548 +130,589 -6,867 (c) Model selection 2.0 1.5 1.0 0.5 0.0 0.5 1.0 1.5 2.0 Lag (ms) 1 0 1 2 Cross-correlation (d) N - U correlation Figure 2: Predator-prey dynamics analysis. (a) Time series of turbulence intensity Nand flow shear Ushowing quasi-periodic oscillations. (b) Phase space trajectory forming a closed limit cycle. (c) ∆BIC comparison showing predator-prey classified as “linear” despite nonlinear governing equations. (d) Cross-correlation function confirming phase relationship between Nand U. 9