Non-Stationary Analysis of Plasma Turbulence-Zonal Flow Interactions Using Time-Varying Models

Abstract

Understanding the interaction between plasma turbulence and zonal flow is essential for structure formation and transport regulation in magnetically confined plasmas. Using experimental data from the PANTA linear plasma device, we reduce the high-dimensional dynamics of numerous modes into two variables: turbulence and zonal. We then identify prominent intermittent bursts localized in the zonal-flow component—a phenomenon that, to the best of our knowledge, has not been reported previously. These bursts appear as impulsive, positive-going, low-frequency excursions superimposed on rapid background fluctuations. To address this non-stationarity, we investigate two approaches: a model-driven approach and a data-driven approach. We use a generalized Predator-Prey (P-P) model identified via an adaptive extended Kalman filter to track the evolution of time-varying physical parameters, and we also apply a Time-Varying Coefficient Vector Autoregressive model (TV-VAR) to achieve high-fidelity tracking of the sharp, transient burst dynamics. A distinct reorganization of the directed coupling structure is observed following burst events, including a transition to a regime in which zonal-flow regulation of turbulence becomes transiently dominant. These results indicate that time-varying interaction measures may provide informative features for burst forecasting and offer a basis for understanding regime transitions in transient plasma dynamics.

1.  Introduction

Decomposition into spatial Fourier modes has been widely used to characterize the complex spatiotemporal dynamics of linear plasmas [1–3]. Recently, we applied multivariate time series models to the azimuthal Fourier modes of PANTA data (Plasma Assembly for Nonlinear Turbulence Analysis) and quantified inter-mode causality within a linear analysis framework [4]. Although plasmas inherently exhibit nonlinear behavior and nonlinear analysis is desirable, extending multivariate time series models to nonlinear forms often leads to the “parameter explosion”—a drastic increase in the number of parameters that significantly degrades estimation accuracy.

To overcome these difficulties, we reduce the high-dimensional dynamics of numerous modes into a system of two variables: turbulence and zonal flow (ZF). Specifically, ZF is represented by the $k$ = 0 Fourier mode energy, while turbulence is the mean energy of non-zero modes (detailed in Sec. 2). This reduction allows us to dramatically decrease the degrees of freedom while distinguishing the dynamics of each component. While bursty and intermittent phenomena in plasma turbulence have been extensively reported through theory, simulations [5–7], and experiments [8], these studies have focused on turbulence energy and transport. In the process of reducing the Fourier modes to a two-variable representation, we observe that prominent burst phenomena also occur in the zonal-flow component in this reduced representation.

To investigate the time-varying interactions and directionality between these variables, we adopt two approaches: a model-driven approach and a data-driven approach. For the model-driven approach, we focus on identifying the Predator-Prey (P-P) model [9]. While recent efforts have been made to identify P-P systems using Bayesian inference [10, 11], these studies primarily rely on numerical simulation data and focus on extracting time-invariant parameters. In contrast, characterizing the non-stationary bursts in measurements requires a framework capable of tracking parameter evolution in real time; therefore, we generalize the P-P model to incorporate time-varying parameters. To the best of our knowledge, there have been no reported cases in which time-varying parameters are estimated for this class of models. Therefore, we identify a generalized P-P model using an Adaptive Extended Kalman Filter (AEKF), which allows for time-varying parameters, such as growth rates and coupling terms, while maintaining physical interpretability.

For the data-driven approach, we utilize the Time-Varying Coefficient Vector Autoregressive (TV-VAR) model [12]. Although TV-VAR is fundamentally a linear framework, its time-varying coefficients adapt to changes in the dynamical regime at each time point, providing superior tracking and predictive performance for non-stationary data compared to conventional time-invariant coefficient models. Using these two time-varying methods, we estimate time-varying parameters in the generalized P-P model and use TV-VAR to characterize burst-time dynamics and time-dependent, directional coupling between turbulence and zonal flow in both the time and frequency domains.

2.  Data Acquisition and Preprocessing

Experiments using linearly magnetized plasma in the PANTA device enable fundamental research on plasma turbulence. The device generates plasma with a total length of 4 m and a radius of 6 cm. Depending on the neutral gas pressure ($P_{\mathrm{n}}$), various turbulence states, such as the streamer state or solitary wave state, are realized [3].

The temporal and spatial structure of turbulence can be observed with high resolution through an ultra multi-channel Langmuir probe array [13]. The ion saturation current $I_{\mathrm{is}}$ is considered a primary diagnostic for turbulent electron density (and electron temperature) fluctuations in plasma experiments. The floating potential $V_{\mathrm{f}}$ reflects fluctuations in the plasma potential and is commonly used to characterize zonal flow and large-scale electrostatic structures. In this study, the ion saturation current ($I_{\mathrm{is}}$) and floating potential ($V_{\mathrm{f}}$) data were measured under various neutral gas pressure conditions $P_{\mathrm{n}}$ from 0.1 to 0.6 Pa in steps of 0.05 Pa. These measurements were taken for 0.6 s at a sampling frequency of 10⁵ Hz using probes arranged on a ring with 32 channels.

For the analysis, we used a stable 0.4 s interval from the middle of the 0.6 s discharge. As preprocessing, we performed the following steps for each channel: (1) downsampling to 1,000 Hz with an anti-aliasing filter to focus on the low-frequency dynamics of interest, (2) linear trend removal to match the gain levels, and (3) normalization so that each signal was shifted to have a minimum of zero and scaled to unit variance. This normalization was used for the generalized P-P model estimation, where the analyzed variables were treated as non-negative energy-like quantities. For the TV-VAR estimation, we additionally applied z-score normalization so that the signals had zero mean and unit variance, in order to improve numerical stability.

To investigate the interaction between turbulence and zonal flow, we decompose the fluctuation data $\psi (x,t)$ into a zonal flow component and a turbulent component using spatial Fourier analysis. The spatial Fourier coefficient $\widehat{\psi }(k,t)$ for a wavenumber $k$ is obtained by:

  

$$\widehat{\psi }(k,t)=\sum _x\psi (x,t)e^{-ikx}.$$(1)

The zonal flow energy $U^{\mathrm{zf}}(t)$ is calculated as the power of the $k$ = 0 mode:

  

$$U^{\mathrm{zf}}(t)={\left|\widehat{\psi }(0,t)\right|}^2.$$(2)

The turbulence energy $U^{\mathrm{tb}}(t)$ is evaluated as the mean power of the non-zero wavenumber modes ($k$ ≥ 1):

  

$$U^{\mathrm{tb}}(t)={⟨{\left|\widehat{\psi }(k,t)\right|}^2⟩}_{k\ge 1}.$$(3)

We analyzed two sets of bivariate time-series data: (i) turbulence $I_{\mathrm{is}}$ vs. zonal flow $V_{\mathrm{f}}$, denoted as $I_{\mathrm{is}}^{\mathrm{tb}}$–$V_{\mathrm{f}}^{\mathrm{zf}}$, and (ii) turbulence $V_{\mathrm{f}}$ vs. zonal flow $V_{\mathrm{f}}$, denoted as $V_{\mathrm{f}}^{\mathrm{tb}}$–$V_{\mathrm{f}}^{\mathrm{zf}}$.

Figures 1(a)–(d) show the temporal evolution of the processed data at $P_{\mathrm{n}}$ = 0.3 Pa, which exhibited the most prominent burst phenomena among the measured conditions. A primary observation is the presence of intermittent low-frequency bursts. These events are characterized as impulsive, positive-going excursions superimposed on high-frequency fluctuations. Notably, we discover that these prominent bursts occur not only in $U^{\mathrm{tb}}(t)$ but also specifically within $U^{\mathrm{zf}}(t)$.

Fig. 1.  Results of the generalized Predator-Prey (P-P) model with time-varying parameters inferred using the Adaptive Extended Kalman Filter (AEKF). The left column (panels (a), (c), (e), (g), (i) and (k)) and right column (panels (b), (d), (f), (h), (j) and (l)) compare the results for the $I_{\mathrm{is}}^{\mathrm{tb}}$ – $V_{\mathrm{f}}^{\mathrm{zf}}$ and $V_{\mathrm{f}}^{\mathrm{tb}}$ – $V_{\mathrm{f}}^{\mathrm{zf}}$ datasets, respectively. (a)–(d) Comparison between the turbulence and zonal-flow signals derived from the experimental data (blue) and the AEKF-estimated states (red) based on the reduced-order physical model. The estimation reflects the model’s attempt to capture the underlying trends of the normalized turbulence energy $U^{\mathrm{tb}}$ and zonal flow energy $U^{\mathrm{zf}}$. (e)–(j) Temporal evolution of the identified physical parameters: linear growth rate ${\gamma }_{\mathrm{L}}$ (e) and (f), nonlinear coupling parameter $\alpha$ (g) and (h), and zonal flow damping rate ${\gamma }_{\mathrm{damp}}$ (i) and (j). (k) and (l) Characteristic frequency $f_c$ calculated from the identified parameters. This $f_c$ represents the intrinsic frequency inherent in the assumed physical framework, serving as a basis for discussing the discrepancy between the simplified P-P dynamics and the actual spectral features observed in the data.

The heavy-tailed probability density functions (PDFs) in Fig. 2 suggest significant deviations from Gaussian statistics, characteristic of intermittent burst phenomena. This non-stationary behavior indicates that conventional models with time-invariant parameters or coefficients are inadequate to capture the transient nature of the plasma dynamics, thereby motivating the use of the time-varying analytical frameworks described in the following sections.

Fig. 2.  Probability density functions (PDFs) of the turbulence ($I_{\mathrm{is}}$ and $V_{\mathrm{f}}$) and zonal flow ($V_{\mathrm{f}}$) energy. The prominent heavy tails observed across all distributions correspond to the intermittent, non-Gaussian low-frequency burst events discovered in the experimental data. The slight extension of the smoothed curves into the negative range reflects a smoothing effect of kernel density estimation.

3.  Model and Methodology

3.1  Identification of a generalized P-P model

3.1.1  Stochastic differential equation formulation

The P-P model was originally formulated as a deterministic system to describe the nonlinear energy exchange between turbulence and zonal flow [9]. To enable state estimation from noisy observations and real-time tracking of time-varying parameters via the Kalman filter framework, we extend the model to a stochastic differential equation (SDE) formulation.

We adopt multiplicative noise, which assumes uncertainty proportional to the state magnitude. This formulation is appropriate when the relative (rather than absolute) uncertainty is approximately constant, and it ensures that the energy-related quantities ($U^{\mathrm{tb}}$ and $U^{\mathrm{zf}}$) remain non-negative. The temporal evolution of the normalized turbulence energy $U^{\mathrm{tb}}$ and the zonal flow energy $U^{\mathrm{zf}}$ is governed by the following SDEs:

  

$$\begin{array}{cc}\hfill d{U}^{\mathrm{tb}}=& [{\gamma }_{\mathrm{L}}(t)U^{\mathrm{tb}}-\alpha (t)U^{\mathrm{tb}}{U}^{\mathrm{zf}}-{\gamma }_2(t){(U^{\mathrm{tb}})}^2]\mathit{d}t\hfill \\ & +U^{\mathrm{tb}}\sqrt{Q_{\mathrm{tb}}}\mathit{d}{W}_{\mathrm{tb}},\hfill \\ \\ \hfill d{U}^{\mathrm{zf}}=& [\alpha (t)U^{\mathrm{tb}}{U}^{\mathrm{zf}}-{\gamma }_{\mathrm{damp}}(t)U^{\mathrm{zf}}]\mathit{d}t\hfill \\ & +U^{\mathrm{zf}}\sqrt{Q_{\mathrm{zf}}}\mathit{d}{W}_{\mathrm{zf}},\hfill \end{array}$$(4)

where ${\gamma }_{\mathrm{L}}(t)$, $\alpha (t)$, ${\gamma }_2(t)$ and ${\gamma }_{\mathrm{damp}}(t)$ are time-varying parameters representing the linear growth rate, the nonlinear coupling parameter, the nonlinear saturation of turbulence, and the damping rate of the zonal flow, respectively. In this study, ${\gamma }_2$ is treated as a constant parameter based on preliminary analysis, but is retained in the general formulation of the model for notational completeness. The terms $d{W}_{\mathrm{tb}}$ and $d{W}_{\mathrm{zf}}$ denote independent Wiener processes, and $Q_{\mathrm{tb}}$ and $Q_{\mathrm{zf}}$ represent the noise intensities.

3.1.2  Numerical discretization

For the implementation of the AEKF, the deterministic part of the SDEs is discretized using the exponential Euler method. Defining the instantaneous growth rates as

  

$$\begin{array}{cc}\hfill {\mu }_{\mathrm{tb}}(t_n)=& {\gamma }_{\mathrm{L}}(t_n)-{\gamma }_2(t_n)U_n^{\mathrm{tb}}-\alpha (t_n)U_n^{\mathrm{zf}},\hfill \\ \hfill {\mu }_{\mathrm{zf}}(t_n)=& \alpha (t_n)U_n^{\mathrm{tb}}-{\gamma }_{\mathrm{damp}}(t_n),\hfill \end{array}$$(5)

the state prediction is given by

  

$$\begin{array}{cc}\hfill U_{n+1}^{\mathrm{tb}}=& U_n^{\mathrm{tb}}\exp ({\mu }_{\mathrm{tb}}(t_n)\mathrm{\Delta }t),\hfill \\ \hfill U_{n+1}^{\mathrm{zf}}=& U_n^{\mathrm{zf}}\exp ({\mu }_{\mathrm{zf}}(t_n)\mathrm{\Delta }t).\hfill \end{array}$$(6)

The stochastic contribution is incorporated through the process noise covariance matrix in the EKF covariance propagation, with

  

$$Q_n=\mathrm{diag}({(U_{n+1}^{\mathrm{tb}})}^2{Q}_{\mathrm{tb}},{(U_{n+1}^{\mathrm{zf}})}^2{Q}_{\mathrm{zf}})\mathrm{\Delta }t,$$

consistent with the multiplicative noise formulation.

3.1.3  Parameter estimation via AEKF

The system state vector is augmented to include the unknown time-varying parameter vector $\mathit{\theta }(t)={[{\gamma }_{\mathrm{L}}(t),\alpha (t),{\gamma }_2(t),{\gamma }_{\mathrm{damp}}(t)]}^{\mathsf{T}}$. In the discrete-time AEKF implementation, this parameter vector is represented as ${\mathit{\theta }}_n$, and its temporal evolution is modeled as a random walk process:

  

$${\mathit{\theta }}_{n+1}={\mathit{\theta }}_n+{\mathit{\eta }}_n,{\mathit{\eta }}_n\sim \mathcal{N}(\mathbf{0},Q_{\theta }\mathrm{\Delta }t),$$(7)

where $Q_{\theta }=\mathrm{diag}(q_{{\gamma }_{\mathrm{L}}},q_{\alpha },q_{{\gamma }_2},q_{{\gamma }_{\mathrm{damp}}})$ is the process noise intensity matrix for the time-varying parameters, with $q_{{\gamma }_{\mathrm{L}}},q_{\alpha },q_{{\gamma }_2},q_{{\gamma }_{\mathrm{damp}}}$ > 0. The AEKF recursively updates the joint state-parameter vector [14], enabling the simultaneous tracking of the latent plasma dynamics and the associated time-varying parameters from noisy experimental observations.

We employed a multi-rate AEKF, in which the model state is propagated using a finer internal time discretization than that of the measurement updates, whereas the Kalman gain is updated only every 200 prediction steps. By separating the prediction and update time scales, the filter remains numerically stable while focusing on the low-frequency dynamics described by the model.

3.1.4  Linear stability analysis and time-varying characteristic frequency

To provide a physical interpretation of the time-varying parameters inferred by the AEKF, we derive a time-dependent characteristic frequency of the coupled turbulence-zonal flow system based on the deterministic drift part of the SDEs, i.e., in the zero-noise limit. Since the model parameters, ${\gamma }_{\mathrm{L}}(t)$, $\alpha (t)$, ${\gamma }_2(t)$, and ${\gamma }_{\mathrm{damp}}(t)$, evolve in time, the local stability properties of the system are also time-dependent. We therefore adopt a frozen-time (or quasi-stationary) approximation, in which the parameters are treated as temporarily constant at each time $t$.

Let the deterministic drift vector at time $t$ be defined as

  

$$\begin{array}{cc}\hfill \mathbf{F}(U^{\mathrm{tb}},U^{\mathrm{zf}};t)=& \left(\begin{array}{c}{F}_{\mathrm{tb}}(U^{\mathrm{tb}},U^{\mathrm{zf}};t)\\ F_{\mathrm{zf}}(U^{\mathrm{tb}},U^{\mathrm{zf}};t)\end{array}\right)\hfill \\ \hfill =& \left(\begin{array}{c}\hfill {\gamma }_{\mathrm{L}}(t)U^{\mathrm{tb}}-\alpha (t)U^{\mathrm{tb}}{U}^{\mathrm{zf}}-{\gamma }_2(t){(U^{\mathrm{tb}})}^2\hfill \\ \hfill \alpha (t)U^{\mathrm{tb}}{U}^{\mathrm{zf}}-{\gamma }_{\mathrm{damp}}(t)U^{\mathrm{zf}}\hfill \end{array}\right).\hfill \end{array}$$(8)

At each time $t$, we consider the non-trivial equilibrium point $(U^{\mathrm{tb}\ast }(t),U^{\mathrm{zf}\ast }(t))$, defined by $\mathbf{F}=\mathbf{0}$ with $U^{\mathrm{zf}\ast }(t)\ne 0$, which requires $\alpha (t)\ne 0$ and $U^{\mathrm{tb}\ast }(t)\ne 0$, yielding

  

$$\begin{array}{c}{U}^{\mathrm{tb}\ast }(t)=\frac{{\gamma }_{\mathrm{damp}}(t)}{\alpha (t)},\hfill \\ U^{\mathrm{zf}\ast }(t)=\frac{{\gamma }_{\mathrm{L}}(t)-{\gamma }_2(t)U^{\mathrm{tb}\ast }(t)}{\alpha (t)}.\hfill \end{array}$$(9)

The Jacobian matrix evaluated at $(U^{\mathrm{tb}\ast }(t),U^{\mathrm{zf}\ast }(t))$ is given by

  

$$\begin{array}{cc}\hfill \mathbf{J}(t)=& \left(\begin{array}{cc}\frac{\partial F_{\mathrm{tb}}}{\partial U^{\mathrm{tb}}}& \frac{\partial F_{\mathrm{tb}}}{\partial U^{\mathrm{zf}}}\\ \frac{\partial F_{\mathrm{zf}}}{\partial U^{\mathrm{tb}}}& \frac{\partial F_{\mathrm{zf}}}{\partial U^{\mathrm{zf}}}\end{array}\right)\hfill \\ \hfill =& \left(\begin{array}{cc}-{\gamma }_2(t)U^{\mathrm{tb}\ast }(t)& -\alpha (t)U^{\mathrm{tb}\ast }(t)\\ \alpha (t)U^{\mathrm{zf}\ast }(t)& 0\end{array}\right).\hfill \end{array}$$(10)

The eigenvalues $\lambda (t)$ of $\mathbf{J}(t)$ satisfy the characteristic equation

  

$${\lambda }^2(t)-\mathrm{Tr}(\mathbf{J}(t))\lambda (t)+\det (\mathbf{J}(t))=0,$$(11)

where

  

$$\begin{array}{cc}\hfill \mathrm{Tr}(\mathbf{J}(t))=& -{\gamma }_2(t)U^{\mathrm{tb}\ast }(t),\hfill \\ \hfill \det (\mathbf{J}(t))=& {\alpha }^2(t)U^{\mathrm{tb}\ast }(t)U^{\mathrm{zf}\ast }(t).\hfill \end{array}$$(12)

Thus,

  

$$\begin{array}{cc}\hfill \lambda (t)=& \frac{-{\gamma }_2(t)U^{\mathrm{tb}\ast }(t)}{2}\hfill \\ & \pm \frac{1}{2}\sqrt{{({\gamma }_2(t)U^{\mathrm{tb}\ast }(t))}^2-4{\alpha }^2(t)U^{\mathrm{tb}\ast }(t)U^{\mathrm{zf}\ast }(t)}.\hfill \end{array}$$(13)

When the discriminant is negative, the local dynamics exhibit damped oscillatory behavior. The corresponding instantaneous angular frequency $\omega (t)$ is given by

  

$$\omega (t)=\sqrt{{\alpha }^2(t)U^{\mathrm{tb}\ast }(t)U^{\mathrm{zf}\ast }(t)-{\left(\frac{{\gamma }_2(t)U^{\mathrm{tb}\ast }(t)}{2}\right)}^2}.$$(14)

Finally, we define the time-varying characteristic frequency as

  

$$f_c(t)=\frac{\omega (t)}{2\pi }.$$(15)

3.2  Estimation of a TV-VAR model

3.2.1  TV-VAR model representation

To capture the transient and non-stationary interactions between turbulence energy and zonal flow, we employ a TV-VAR model. In this framework, the observed bivariate signal vector ${\mathbf{y}}_t={[U^{\mathrm{tb}}(t),U^{\mathrm{zf}}(t)]}^{\mathsf{T}}$ is modeled as a linear combination of its own previous states with coefficients that evolve over time:

  

$${\mathbf{y}}_t=\sum _{m=1}^p{\mathbf{A}}_{m,t}{\mathbf{y}}_{t-m}+{\mathit{ϵ}}_t,$$(16)

where $p$ is the model order, ${\mathbf{A}}_{m,t}$ are 2 × 2 time-varying coefficient matrices at lag $m$, and ${\mathit{ϵ}}_t$ is the white noise innovation vector with time-varying covariance ${\mathbf{\Sigma }}_t$. The optimal model order was selected separately for each data combination: $p$ = 3 for the $I_{\mathrm{is}}^{\mathrm{tb}}$ – $V_{\mathrm{f}}^{\mathrm{zf}}$ pair, and $p$ = 4 for the $V_{\mathrm{f}}^{\mathrm{tb}}$ – $V_{\mathrm{f}}^{\mathrm{zf}}$ pair, using the Akaike Information Criterion (AIC) [15].

3.2.2  State-space estimation via kalman filter

The time-varying coefficients ${\mathbf{A}}_{m,t}$ and covariance ${\mathbf{\Sigma }}_t$ are estimated by treating the coefficients as state variables following a random walk process. A Kalman filter is used to track dynamic changes in the coupling strength and characteristic frequencies during transient events such as bursts.

3.2.3  Time-Varying Impulse Response Function (TV-IRF)

While spectral analysis provides insight into the frequency-domain characteristics of the system, it does not directly reveal the causal and temporal propagation of perturbations between variables. To characterize these time-domain dynamics, we employ the Time-Varying Impulse Response Function (TV-IRF).

Based on the estimated TV-VAR model, the impulse response at each time step $t$ is evaluated under the assumption of local stationarity. Specifically, the time-varying autoregressive coefficients ${\{{\mathbf{A}}_{m,t}\}}_{m=1}^p$ are treated as frozen at time $t$, and the response of the linear system to an instantaneous unit impulse is computed recursively. The TV-IRF matrices ${\mathbf{\Phi }}_t(h)$ are defined as:

  

$${\mathbf{\Phi }}_t(h)=\sum _{m=1}^{\min (h,p)}{\mathbf{A}}_{m,t}{\mathbf{\Phi }}_t(h-m),$$(17)

where $h$ ≥ 1 denotes the discrete time lag (horizon), ${\mathbf{\Phi }}_t(0)=\mathbf{I}$, and $\mathbf{I}$ is the identity matrix.

Each element ${[{\mathbf{\Phi }}_t(h)]}_{ij}$ represents the response of variable $i$ at lag $h$ to a unit impulse applied to variable $j$ at time $t$. By evaluating ${\mathbf{\Phi }}_t(h)$ for all $t$, we obtain a time-resolved description of how perturbations propagate through the coupled turbulence-zonal flow system.

In contrast to conventional stationary IRF analysis, the TV–IRF explicitly captures the temporal evolution of interaction pathways and response timescales. This enables us to identify periods during which the influence of turbulence on zonal flow, or vice versa, is transiently enhanced or suppressed, providing a time-resolved measure of directional coupling in the time domain.

3.2.4  Time-varying spectral analysis (TV-PSD and TV-RPC)

The frequency-domain characteristics of the interaction between turbulence and zonal flow were examined based on the estimated TV-VAR parameters. To characterize the spectral structure and relative power contributions within the system, we employed Time-Varying Power Spectral Density (TV-PSD) and Time-Varying Relative Power Contribution (TV-RPC) analyses.

First, assuming local stationarity at each time step $t$, the time-varying transfer function matrix $\mathbf{H}(f,t)$ is defined as the inverse of the Fourier-transformed autoregressive coefficient matrix:

  

$$\mathbf{H}(f,t)={\left(\mathbf{I}-\sum _{m=1}^p{\mathbf{A}}_{m,t}{e}^{-i2\pi fm/F_s}\right)}^{\mathrm{-1}},$$(18)

where $F_s$ denotes the sampling frequency and $\mathbf{I}$ is the identity matrix.

Using this transfer function, the TV-PSD matrix $\mathbf{P}(f,t)$ is computed. A critical assumption in the standard noise contribution analysis is that the innovation processes of the variables are mutually uncorrelated. Accordingly, we approximated the noise covariance matrix ${\mathbf{\Sigma }}_t$ by its diagonal elements, i.e., ${\mathrm{\Sigma }}_{jj,t}=\mathrm{Var}[{ϵ}_{j,t}]$. Under this orthogonality assumption, the auto-PSD of the $i$-th variable, $P_{ii}(f,t)$, can be decomposed into the sum of partial power contributions from each noise source:

  

$$P_{ii}(f,t)=\sum _{j=1}^K{\left|H_{ij}(f,t)\right|}^2{\mathrm{\Sigma }}_{jj,t},$$(19)

where $H_{ij}(f,t)$ represents the element of the transfer matrix from the $j$-th input to the $i$-th output, and $K$ denotes the total number of variables.

To quantify the relative power contributions across variables in the frequency domain, we calculated the TV-RPC. We define the partial power contribution $Q_{ij}(f,t)$, which represents the portion of the power of variable $i$ induced by the intrinsic fluctuation (innovation) of variable $j$, as:

  

$$Q_{ij}(f,t)={\left|H_{ij}(f,t)\right|}^2{\mathrm{\Sigma }}_{jj,t}.$$(20)

Using this term, the TV-RPC $R_{ij}(f,t)$ is defined as the ratio of this partial contribution to the total power:

  

$$R_{ij}(f,t)=\frac{Q_{ij}(f,t)}{P_{ii}(f,t)}.$$(21)

By definition, the summation of the contributions from all variables equals unity (${\sum }_j{R}_{ij}(f,t)$ = 1). This metric provides a normalized measure of how the spectral power of each variable is distributed among the innovation sources through the time-varying linear dynamics.

4.  Results

Figure 1 presents the results of applying the generalized P-P model to the turbulence and zonal flow data using the AEKF. The left and right columns of Fig. 1 correspond to the $I_{\mathrm{is}}^{\mathrm{tb}}$ – $V_{\mathrm{f}}^{\mathrm{zf}}$ and $V_{\mathrm{f}}^{\mathrm{tb}}$ – $V_{\mathrm{f}}^{\mathrm{zf}}$ datasets, respectively. In Figs. 1(a)–(d), the turbulence and zonal flow time series (blue) are compared with the corresponding state trajectories obtained by fitting the generalized P-P model using the AEKF (red).

In both datasets, the AEKF estimates reproduce only the slowly varying components of the observed signals and fail to capture the sharp, high-amplitude burst events. As seen in Figs. 1(a)–(d), the generalized P-P model therefore tracks primarily the low-frequency envelope of the burst activity rather than its fast-time structure.

The temporal evolution of the inferred model parameters is shown in Figs. 1(e)–(j). Despite the limited ability of the generalized P-P model to reproduce individual bursts, systematic parameter variations are observed in association with the burst dynamics. Around $t$ ≃ 0.14 s, a small burst is detected in the zonal flow signal (Figs. 1(c) and (d)) without a simultaneous response in the turbulence energy (Figs. 1(a) and (b)). Following this event, the linear growth rate ${\gamma }_{\mathrm{L}}$ and the coupling parameter $\alpha$ begin to increase (Figs. 1(e)–(h)), while the zonal-flow damping rate ${\gamma }_{\mathrm{damp}}$ exhibits a decrease or a temporary suppression of its growth (Figs. 1(i) and (j)).

After the pronounced burst occurring around $t$ ≃ 0.2 s in both the turbulence and zonal flow signals (Figs. 1(a)–(d)), the model parameters exhibit distinct temporal behaviors (Figs. 1(e)–(j)). The linear growth rate ${\gamma }_{\mathrm{L}}$ and the damping rate ${\gamma }_{\mathrm{damp}}$ tend to increase following the burst, whereas the nonlinear coupling parameter $\alpha$ shows a non-monotonic evolution, with an increase prior to the burst and a step-like decrease at its onset.

The characteristic frequency $f_c$, shown in Figs. 1(k) and (l), primarily reflects variations in the nonlinear coupling parameter $\alpha$ under the present modeling assumption that ${\gamma }_2$ is fixed at 0.1 (i.e., time-invariant). This trend is consistent with enhanced turbulence–zonal flow interaction.

To further investigate the fast dynamics and causal structures that the reduced-order physical model could not fully capture, we applied the TV-VAR model to the same datasets. Fig. 3 summarizes the time series estimation and the TV-IRFs.

Fig. 3.  TV-VAR model analysis: time series estimation and impulse response functions. The left and right columns correspond to the $I_{\mathrm{is}}^{\mathrm{tb}}$ – $V_{\mathrm{f}}^{\mathrm{zf}}$ and $V_{\mathrm{f}}^{\mathrm{tb}}$ – $V_{\mathrm{f}}^{\mathrm{zf}}$ datasets, respectively. (a)–(d) Comparison between the turbulence and zonal-flow signals derived from the experimental data (blue) and the TV-VAR model estimation (red lines). (e)–(l) Temporal evolution of the impulse response functions (IRFs) derived from the model coefficients. Panels (e) and (f) show the self-response of turbulence, while (g) and (h) show the self-response of zonal flow. The cross-responses are displayed in the bottom panels: (i) and (j) represent the response of zonal flow to a turbulence impulse, and (k) and (l) represent the response of turbulence to a zonal flow impulse. The horizontal axis corresponds to the discharge time, the vertical axis represents the time lag of the response, and the color scale indicates the response intensity.

First, the time series comparison in Figs. 3(a)–(d) demonstrates the high fidelity of the TV-VAR model. In contrast to the generalized P-P model results shown in Fig. 1, the TV-VAR reconstruction (red) almost perfectly tracks the experimental data (blue), accurately reproducing not only the baseline trends but also the sharp, high-frequency burst events observed after $t$ ≃ 0.2 s. This high reproducibility indicates that the TV-VAR model captures the observed dynamics sufficiently well to support subsequent impulse-response and spectral analyses.

For the self-responses of turbulence (Figs. 3(e) and (f)), the response distribution shifts toward shorter lags after $t$ ≃ 0.2 s, while the contribution from longer lags becomes weaker. This indicates that the turbulence dynamics become more dominated by short-memory responses following the burst events.

In contrast, the self-responses of the zonal flow (Figs. 3(g) and (h)) exhibit an enhancement at relatively longer lags after $t$ ≃ 0.2 s, suggesting an increased persistence of the zonal-flow state over extended time delays.

The cross-responses show a clear asymmetry between the two directions. The response from turbulence to zonal flow (Figs. 3(i) and (j)) weakens over a broad range of lags after $t$ ≃ 0.2 s, indicating a reduction in the direct influence of turbulence fluctuations on the subsequent zonal-flow evolution. In contrast, the response from zonal flow to turbulence (Figs. 3(k) and (l)) becomes markedly enhanced after $t$ ≃ 0.2 s, with a pronounced increase in response amplitude over specific lag ranges, depending on the dataset.

Taken together, these results indicate a reorganization of the low-frequency coupling structure between turbulence and zonal flow following the burst events. After $t$ ≃ 0.2 s, the interaction becomes increasingly dominated by the influence of zonal flow on turbulence, while the opposite direction is suppressed. This behavior is consistent with a transition to a regime in which the feedback from zonal flow plays a more prominent role in regulating the turbulence dynamics.

Figures 4(e)–(h) show the TV-PSD results, where recurrent, common temporal changes are observed in both variables in association with the simultaneous bursts in turbulence and zonal flow (approximately at $t$ ≃ 0.20, 0.27, 0.34, and 0.38 s). Specifically, prior to each burst event, the low-frequency components (roughly below 50 Hz) decrease relative to the surrounding periods, and they increase again after the burst; this pattern repeats multiple times. Such a cyclic modulation of low-frequency power may indicate that the system repeatedly switches between quasi-quiescent and active phases, or that a quasi-periodic redistribution process exists in which low-frequency components are suppressed before bursts and subsequently reinforced.

Fig. 4.  Time-varying spectral analysis based on the TV-VAR model. The left and right columns correspond to the $I_{\mathrm{is}}^{\mathrm{tb}}$ – $V_{\mathrm{f}}^{\mathrm{zf}}$ and $V_{\mathrm{f}}^{\mathrm{tb}}$ – $V_{\mathrm{f}}^{\mathrm{zf}}$ datasets, respectively. (a)–(d) Comparison of time series waveforms between the experimental data (blue) and the TV-VAR model estimation (red). (e)–(h) Time-Varying Power Spectral Density (TV-PSD) of turbulence and zonal flow, where the color scale indicates the power intensity. (i)–(l) Time-Varying Relative Power Contribution (TV-RPC), representing the relative contribution of each variable to the spectral power of the other through the time-varying linear coefficients. Panels (i) and (j) show the contribution from turbulence to zonal flow, corresponding to $R_{\mathrm{zf}\leftarrow \mathrm{tb}}(f,t)$. Panels (k) and (l) show the contribution from zonal flow to turbulence, corresponding to $R_{\mathrm{tb}\leftarrow \mathrm{zf}}(f,t)$. The horizontal axis represents time, the vertical axis represents frequency, and the color scale in (i)–(l) represents the contribution ratio ranging from 0 to 1.

In terms of the TV-RPC, for the $I_{\mathrm{is}}^{\mathrm{tb}}$ – $V_{\mathrm{f}}^{\mathrm{zf}}$ dataset (Fig. 4(i)), the relative contribution from turbulence ($I_{\mathrm{is}}^{\mathrm{tb}}$ ) to zonal flow ($V_{\mathrm{f}}^{\mathrm{zf}}$), $R_{\mathrm{zf}\leftarrow \mathrm{tb}}(f,t)$, shows a gradual increase in the low-frequency range (roughly below 30 Hz) before the first burst ($t$ ≃ 0.2 s), whereas it drops sharply over a broad frequency range after $t$ ≃ 0.2 s and remains nearly absent thereafter. A similar tendency is also observed for the $V_{\mathrm{f}}^{\mathrm{tb}}$ – $V_{\mathrm{f}}^{\mathrm{zf}}$ dataset (Fig. 4(j)): before $t$ ≃ 0.2 s, non-negligible contributions appear in the bands below 10 Hz and around 50 Hz, but these contributions decrease almost instantaneously at the onset of the first burst. This behavior may indicate that, during quiescent intervals, fluctuations on the turbulence side are more readily reflected in the spectral power of the zonal-flow signal. Conversely, during burst intervals, the abrupt change in the contribution pattern implies that additional processes not fully captured by the present linear contribution decomposition may become relevant.

In contrast, the contribution from zonal flow to turbulence is broadly distributed across frequencies for the $I_{\mathrm{is}}^{\mathrm{tb}}$ – $V_{\mathrm{f}}^{\mathrm{zf}}$ dataset (Fig. 4(k)), exhibiting an instantaneous enhancement in the low-frequency range (roughly below 50 Hz) at the burst onsets after $t$ ≃ 0.2 s. Similarly, the $V_{\mathrm{f}}^{\mathrm{tb}}$ – $V_{\mathrm{f}}^{\mathrm{zf}}$ dataset (Fig. 4(l)) shows a consistent enhancement of low-frequency contributions synchronized with bursts, while no clear contributions appear in the higher frequency range (roughly above 50 Hz) regardless of the presence or absence of bursts. Overall, these time-varying reorganizations of the spectral structure and relative contributions indicate that the interaction pattern between turbulence and zonal flow changes across burst events. In particular, a regime may transiently emerge after bursts in which zonal-flow fluctuations are more strongly reflected in the turbulence dynamics, consistent with a temporarily strengthened zonal-flow regulation of turbulence.

Around $t$ ≃ 0.35 s, a back-transition is also observed in the interaction pattern, particularly in the cross-response structures in Figs. 3(i)–(l) and the relative contribution patterns in Figs. 4(i)–(l). This change is qualitatively different from the earlier burst-associated transition around $t$ ≃ 0.2 s.

5.  Discussion

The reduced P-P model, even with time-varying parameters, does not provide a high-fidelity description of the burst waveforms when fitted to the present experimental time series under the present identification setting. This limitation is intrinsic to the low-dimensional formulation: the model is designed to represent coarse-grained turbulence–zonal-flow exchange with a small number of effective coefficients rather than to reproduce fast, intermittent events. Nevertheless, the time-varying parameters do respond to burst occurrences, attempting to track the underlying changes in dynamical structure.

In particular, the nonlinear coupling parameter $\alpha$ exhibited a distinct, non-monotonic temporal evolution associated with the burst events. As shown in Figs. 1(g) and (h), $\alpha$ increases prior to the burst and undergoes a step-like decrease at its onset, indicating a transient reorganization of the interaction between turbulence and zonal flow. This non-monotonic behavior suggests that $\alpha$ reflects changes in the structure of energy transfer between turbulence and zonal flow, rather than simply a strengthening of nonlinear coupling. Accordingly, the P-P framework, while currently limited in reproducing detailed waveforms, may still offer potential for improvement through refinements such as incorporating higher-order nonlinearities or additional state variables.

In contrast, the TV-VAR framework provides an accurate empirical representation of the observed dynamics and yields time-resolved measures of interaction in the time and frequency domains. Its advantage lies in describing the data with high flexibility and in revealing when and how the linearized coupling structure reorganizes. At the same time, TV-VAR parameters are primarily descriptive, and their physical interpretation is generally non-unique, particularly in the presence of intermittent bursts and potential nonlinear effects.

The back-transition around $t$ ≃ 0.35 s suggests a further reorganization of the post-burst state rather than a simple recovery to the pre-burst condition. This may reflect burst-induced changes in the background plasma or turbulence properties. A Kelvin–Helmholtz-type instability is also possible if the zonal-flow shear is enhanced, but confirming this will require radial-structure measurements.

These distinct roles motivate future development of hybrid models that retain physical structure while achieving improved expressiveness for burst dynamics. One direction is to embed P-P-inspired physical constraints into a flexible state-space model or a time-varying linear framework, so that data-driven components account for burst-triggered deviations while the physical core preserves interpretability. Such a hybrid formulation is also a natural stepping stone toward burst prediction.

In particular, the increase in the P-P coupling parameter $\alpha$ prior to burst onset and the decrease in the low-frequency components in the TV-PSD preceding bursts suggest that these quantities may serve as candidate precursors of impending burst events. Physically constrained latent states and time-varying interaction measures could therefore be used as predictive features to estimate near-future burst probability or onset time. Evaluating predictive performance and robustness across operating conditions (including $P_{\mathrm{n}}$) is an important subject for future work.

Acknowledgements

This work was supported by JSPS KAKENHI Grant Numbers JP24K15079 and JP25K00986, and by the grant of the OML Project by the National Institutes of Natural Sciences (NINS program No. OML012511). This work was also supported in part by the Collaborative Research Program of Research Institute for Applied Mechanics, Kyushu University (2025S2-CD-3), and was performed under the auspices of the NIFS Collaboration Research program (NIFS25KISS041).

References

• [1]   T.  Kobayashi et al., Phys. Plasmas 23, 102311 (2016).
• [2]   S.  Inagaki et al., Sci. Rep. 6, 22189 (2016).
• [3]   T.  Kobayashi et al., Plasma Fusion Res. 12, 1401019 (2017).
• [4]   F.  Miwakeichi and  M.  Sasaki, Phys. Plasmas 31, 082307 (2024).
• [5]   S.  Benkadda et al., Nucl. Fusion 41, 995 (2001).
• [6]   M.A.  Malkov et al., Phys. Plasmas 8, 5073 (2001).
• [7]   E.  Kim et al., Nucl. Fusion 43, 961 (2003).
• [8]   Y.  Kosuga et al., Phys. Plasmas 29, 122301 (2022).
• [9]   P.H.  Diamond et al., Plasma Phys. Control. Fusion 47, R35 (2005).
• [10]   M.  Sasaki et al., Plasma Phys. Control. Fusion 67, 095011 (2025).
• [11]   J.C.  Huang et al., Nucl. Fusion 66, 036010 (2026).
• [12]   X.  Jiang and  G.  Kitagawa, Signal Process. 33, 315 (1993).
• [13]   T.  Yamada et al., Rev. Sci. Instrum. 78, 123501 (2007).
• [14]   R.K.  Mehra, IEEE Trans. Autom. Control 17, 693 (1972).
• [15]   H.  Akaike, IEEE Trans. Autom. Control 19, 716 (1974).