1. Introduction
Turbulence generally evolves through distinct stages, beginning with the growth of small-amplitude fluctuations driven by linear instabilities and progressing toward a statistically stationary, broadband state. In the early, or developing, phase, energy is injected into a limited set of unstable modes, producing coherent spectral peaks and a time-dependent redistribution of power as nonlinear interactions gradually become important. As the fluctuation amplitudes increase, nonlinear mode coupling transfers energy away from these unstable components, broadening the spectrum and forming a cascade toward damped scales. Eventually, a balance is established between linear growth, nonlinear transfer and dissipation, leading to a fully developed turbulent regime in which the statistical properties of the fluctuations become approximately stationary (Ritz, Powers & Bengtson Reference Ritz, Powers and Bengtson1989; Horton Reference Horton1999).
In fluid and plasma turbulence, the time evolution of fluctuating physical quantities, such as plasma density and electrostatic potential, are often approximated as a superposition of statistically independent wave components. Under this assumption, the fluctuations can be partially characterised through their power spectra (Smith et al. Reference Smith, Powers and Caldwell1974a ). During the 1970s, linear spectral analysis techniques based on second-order statistics (Smith & Powers Reference Smith and Powers1973; Smith et al. Reference Smith, Powers and Caldwell1974c ), including auto-power and cross-power spectra, were introduced as diagnostic tools for plasma waves and instabilities. These methods have since been widely and successfully used to investigate wave dynamics and spectral properties in laboratory and space plasmas. However, linear spectral methods inherently assume statistical independence between modes and, therefore, cannot capture the phase-coherent interactions generated by nonlinear dynamics. As a result, they provide limited insight into how energy is exchanged among spectral components. To overcome this limitation, higher-order spectral techniques have been developed to reveal nonlinear coupling processes. In particular, bispectral methods, based on third-order statistics, allow the identification of phase-correlated wave triplets and the assessment of nonlinear wave–wave interactions that remain hidden in conventional power spectra. Using the third-order cumulant, it has been demonstrated in both fluids and plasmas that three-wave coupling processes play an essential role in turbulent dynamics (Kim & Powers Reference Kim and Powers1979; Miksad, Jones & Powers Reference Miksad, Jones and Powers1983; Elgar & Guza Reference Elgar and Guza1985). Similar three-wave coupling interactions between coherent modes and broadband turbulence have also been reported in magnetised plasmas, such as the coupling between geodesic acoustic modes and broadband turbulence in the T-10 tokamak (Melnikov et al. Reference Melnikov, Eliseev, Lysenko, Ufimtsev and Zenin2017).
While the bispectrum is effective in identifying phase-coherent triadic interactions among waves, they primarily quantify the presence of nonlinear coupling and do not, by themselves, specify the direction of energy flow. Subsequent studies have shown that if there are only a few coherent modes, then direction of energy transfer can be computed as discussed by Kim et al. (Reference Kim, Beall, Powers and Miksad1980). However, this approach is not readily applicable to more complex situations such as developing or fully developed plasma turbulence, where a large number of interacting modes are present. To address this limitation, Ritz & Powers (Reference Ritz and Powers1986) and Ritz et al. (Reference Ritz, Powers, Miksad and Solis1988, Reference Ritz, Powers and Bengtson1989) introduced the concept of a power transfer function followed by an energy transfer function, which quantifies the net nonlinear energy exchange among spectral modes. This formulation provides both qualitative and quantitative insight into how fluctuation energy is redistributed across the spectrum. However, this method does not provide very accurate results for general fluctuation data because it approximates the fourth-order moments using the square of second-order moments (Ritz et al. Reference Ritz, Powers and Bengtson1989). Moreover, it produces large, non-physical damping coefficients when applied to fully developed turbulent systems, both in simulations (Shen et al. Reference Shen, Li, Li and Li2020) and in measured fluctuation data (Ritz & Powers Reference Ritz and Powers1986). Consequently, Ritz’s method performs reasonably well to estimate quadratic nonlinearity only for developing turbulence models, whereas it fails for fully developed turbulent systems (Shen et al. Reference Shen, Li, Li and Li2020). To address limitations of the Ritz method, Kim et al. (Reference Kim, Durst, Fonck, Fernandez, Ware and Terry1996) introduced a modified formulation that retains fourth-order cumulants and separates the data into ideal and non-ideal components. In this framework, only the ideal part is assumed to participate in linear growth (or damping) and three-wave coupling processes. Further it is observed that, in models of fully developed turbulence, unphysical behaviour in the estimated damping coefficients can arise if non-ideal contributions – such as higher-order nonlinear couplings, systematic errors or measurement noise – are neglected (Ritz et al. Reference Ritz, Powers and Bengtson1989; Shen et al. Reference Shen, Li, Li and Li2020). To mitigate these effects, they imposed a stationarity condition by enforcing equality between the input and output ideal powers, thereby reducing the influence of non-ideal components on the inferred energy transfer (Kim et al. Reference Kim, Durst, Fonck, Fernandez, Ware and Terry1996). Ritz’s method is generally not well suited for fully developed turbulent systems, whereas Kim’s method can be applied to both developing and fully developed turbulence, as discussed in Shen et al. (Reference Shen, Li, Li and Li2020). There are, however, situations in which neither Ritz’s nor Kim’s method performs satisfactorily. For example, in modelling nonlinear drift oscillations of moored vessels under non-Gaussian random sea-wave excitation, Ritz’s method becomes ineffective due to system and measurement errors, while Kim’s method fails in this case due to a strong power imbalance between the input and output spectra, which exhibit different spectral trends, leading to the loss of important spectral information and inaccurate transfer estimates (Kim et al. Reference Kim, Powers, Ritz, Miksad and Fischer1987). In such cases, bispectral analysis can still be applied effectively by employing the complete iterative method, as discussed by Shen et al. (Reference Shen, Li, Li and Li2020).
The physical origin of the three-wave coupling processes mainly comes from the quadratic
$\boldsymbol{E}\times \boldsymbol{B}$
convective nonlinearity inherent in drift-wave (DW) dynamics. Drift waves arise due to the background density gradient
$\boldsymbol{\nabla }n_0$
, where electrostatic potential fluctuations generate an
$\boldsymbol{E}\times \boldsymbol{B}$
velocity that advects the equilibrium density through the term
$\boldsymbol{v}_E {\cdot }\boldsymbol{\nabla }n_0$
(Hasegawa & Wakatani Reference Hasegawa and Wakatani1983; Anderson et al. Reference Anderson, Kim, Hnat and Rafiq2020). In contrast, the Rayleigh–Taylor/interchange mode develops when an effective gravitational (fluids) or magnetic curvature (toroidal plasmas) or centrifugal (linear magnetised plasma) force acts on the same unfavourable density gradient, producing flute-like perturbations that extract energy from the equilibrium (Krall & Trivelpiece Reference Krall and Trivelpiece1973). In both cases, the nonlinear interaction appears through the quadratic term
$\boldsymbol{v}_E {\cdot }\boldsymbol{\nabla }F$
, where,
$F$
represents advected plasma fields such as density or vorticity, which reduces in two-dimensional geometry to the Poisson bracket
$\{\phi , F\}$
and is explicitly quadratic in fluctuation amplitudes (Terry & Diamond Reference Terry and Diamond1985). In Fourier space, this nonlinear term generates triadic interactions satisfying the resonance conditions
$\boldsymbol{k} = \boldsymbol{k}_1 + \boldsymbol{k}_2$
and
$f = f_1 + f_2$
. Finite collisionality further enhances this process by breaking adiabatic electron shielding, increasing the density-potential phase shift, and broadening the unstable spectral bandwidth; this leads to larger fluctuation amplitudes and strengthens the nonlinear
$\boldsymbol{E}\times \boldsymbol{B}$
coupling (Hasegawa & Wakatani Reference Hasegawa and Wakatani1983).
In this work we present the nonlinear three-wave interaction between DW and Rayleigh–Taylor (RT) modes, identified using the digital bispectrum, which directly captures phase coupling arising from the quadratic Poisson-bracket nonlinearity. The efficiency of this nonlinear coupling is enhanced when the DW and RT eigenfunctions spatially coincide, allowing stronger nonlinear coupling and sustained phase locking, as evidenced by the finite bispectrum (Kim & Powers Reference Kim and Powers1979). The existing computational methods for a single-field turbulence model such as the Ritz method and Kim method have been explored to understand the turbulence dynamics. We observe that the Ritz method performs well for developing turbulent systems, its accuracy depends strongly on the statistical properties of the data. In particular, it is reliable for fluctuations that are close to a Gaussian nature, characterised by low skewness and kurtosis. When the kurtosis is large, the method deteriorates because Ritz’s formulation approximates fourth-order moments using second-order statistics, whereas kurtosis directly reflects fourth-order contributions. Consequently, strong departures from Gaussianity introduce errors in the estimated transfer functions. This issue has been systematically investigated in this work. We examine the applicability of the Ritz method to both simulation data and experimental fluctuation measurements, with emphasis on the role of statistical properties in developing plasma turbulence. We show that Ritz’s method provides reliable transfer functions only when the data are close to Gaussian, and that beyond a threshold level of kurtosis it fails to recover physically meaningful transfer functions, the result verified using both simulated and measured fluctuations from the inverse mirror plasma experimental device (IMPED). In contrast, Kim’s method performs robustly regardless of the statistical character of the data. Even for signals with large kurtosis, Kim’s formulation yields consistent and physically correct transfer functions, as demonstrated in both simulations and experiments. Although Kim’s method was originally developed for fully developed, stationary turbulence where input and output powers balance, we apply it here to developing turbulent regimes in which the spectra are non-stationary and not fully coincident but showing similar trends. Despite this, Kim’s method continues to provide accurate transfer functions for both simulated and experimental fluctuation data from IMPED.
The flow of the paper is organised as follows. In § 2 we introduce the spectral analysis tools required to compute linear, quadratic and nonlinear energy transfer functions. Section 3 describes both Ritz’s and Kim’s method, including simulation tests used to validate Ritz’s approach and to examine its applicability with respect to the statistical properties of the data. This section also presents Kim’s method, along with simulation studies demonstrating its performance and its ability to address the limitations of Ritz’s formulation. In § 4 we apply both methods to experimental fluctuation data, including nonlinear transfer function and growth rate analyses at two radial locations, thereby assessing the validity of each approach in experimental data. Finally, §§ 5 and 6 present the discussion and conclusions, respectively, while § 7 presents the outlook and future directions of the present work.
2. Spectral analysis framework for nonlinear energy transfer function estimation
This section establishes the spectral analysis framework required for estimating nonlinear energy transfer functions. Linear interactions are characterised using second-order spectral quantities such as the power spectral density, which describe the distribution of energy across frequencies. Nonlinear interactions, arising from quadratic phase coupling and triadic interactions, are captured using higher-order spectral measures. We consider time-series signals divided into
$M$
ensembles, each containing
$N$
data points, and estimate Fourier-based spectral quantities through ensemble averaging. In addition to auto-spectral measures, cross-spectral quantities are introduced to characterise inter-signal coupling and associated energy transfer mechanisms. The definitions of these spectral quantities are provided below.
2.1. Auto-power spectral density
It describes how the signal power is distributed across frequencies. It is defined as (Smith et al. Reference Smith, Powers and Caldwell1974a )
\begin{align} P_{a}(f) = \frac {1}{M} \sum _{l=1}^M X_{\kern-2pt f}^{*l} X_{\kern-2pt f}^l , \end{align}
where
$X_f^l$
denotes the Fourier transform of the
$l$
th ensemble and
$*$
represents complex conjugation.
2.2. Cross-power spectral density
For two time-series signals
$x(t)$
and
$y(t)$
with Fourier transforms
$X_f$
and
$Y_f$
, respectively, the cross-power spectral density is defined as (Smith et al. Reference Smith, Powers and Caldwell1974a
)
\begin{align} P_{c}(f) = \frac {1}{M} \sum _{l=1}^M X_{\kern-2pt f}^{*l} Y_{\kern-2pt f}^l . \end{align}
The cross-power spectral density quantifies the correlation and phase relationship between the two signals at each frequency and is generally complex valued.
2.3. Auto-bispectrum
The auto-bispectrum is a third-order spectral quantity that captures nonlinear interactions and phase coupling between different frequency components. It is defined as (Kim & Powers Reference Kim and Powers1979)
\begin{align} B_a(f_1,f_2) = \frac {1}{M}\sum _{l=1}^M X_{\kern-2pt f_1}^l \, X_{\kern-2pt f_2}^l \, X_{\kern-2pt f_1+f_2}^{*l}. \end{align}
The auto-bispectrum quantifies phase-coupled interactions among frequency triplets
$(f_1, f_2, f_1+f_2)$
and is sensitive to quadratic nonlinearities in the signal.
2.4. Cross-bispectrum
For two signals
$x(t)$
and
$y(t)$
, the cross-bispectrum is defined as (Kim & Powers Reference Kim and Powers1979)
\begin{align} B_{c}(f_1,f_2) = \frac {1}{M}\sum _{l=1}^M X_{\kern-2pt f_1}^l \, X_{\kern-2pt f_2}^l \, Y_{f_1+f_2}^{*l}. \end{align}
The cross-bispectrum characterises quadratic phase coupling between frequency components of one signal and their nonlinear interaction with another signal.
2.5. Auto-bicoherence
To obtain a normalised measure of phase coupling that is independent of signal amplitude, the bicoherence is defined as (Kim & Powers Reference Kim and Powers1979)
\begin{align} b_a(f_1,f_2) = \frac {\big| ({1}/{M})\sum _{l=1}^M X_{\kern-2pt f_1}^l X_{\kern-2pt f_2}^l X_{\kern-2pt f_1+f_2}^{*l} \big|} {\big(({1}/{M})\sum _{l=1}^M |X_{\kern-2pt f_1}^l X_{\kern-2pt f_2}^l|^2\big)^{1/2} \big(({1}/{M})\sum _{l=1}^M |X_{f_1+f_2}^l|^2\big)^{1/2}} . \end{align}
The auto-bicoherence measures the normalised strength of phase coupling between
$f_1$
,
$f_2$
and
$f_1+f_2$
, with values ranging from
$0$
(no coupling) to
$1$
(perfect phase coupling).
2.6. Cross-bicoherence
To obtain a normalised measure of quadratic phase coupling between two different signals that is independent of signal amplitude, the cross-bicoherence is defined as (Kim & Powers Reference Kim and Powers1979)
\begin{align} b_{c}(f_1,f_2) = \frac {\big| ({1}/{M})\sum _{l=1}^M X_{\kern-2pt f_1}^l X_{\kern-2pt f_2}^l Y_{f_1+f_2}^{*l} \big|} {\big(({1}/{M})\sum _{l=1}^M |X_{\kern-2pt f_1}^l X_{\kern-2pt f_2}^l|^2\big)^{1/2} \big(({1}/{M})\sum _{l=1}^M |Y_{f_1+f_2}^l|^2\big)^{1/2}} . \end{align}
The cross-bicoherence measures the normalised strength of phase coupling between frequency components
$f_1$
and
$f_2$
of one signal and their nonlinear interaction with the component
$f_1+f_2$
of another signal. Its values range from
$0$
(no phase coupling) to
$1$
(perfect quadratic phase coupling).
For random (uncorrelated) signals, the expected bicoherence is not exactly zero due to finite ensemble averaging. Therefore, the statistical significance level for the bicoherence is given by (Kim & Powers Reference Kim and Powers1979; Nikias & Petropulu Reference Nikias and Petropulu1993)
\begin{align} b_{\mathrm{sig}} \approx \frac {1}{\sqrt {M}}. \end{align}
The nonlinear coupling is considered statistically significant when
$b(f_{1},f_{2}) \gt {1}/{\sqrt {M}}$
. Also, for squared bicoherence, the statistical significance level is written as:
$b_{\mathrm{sig}}^{2} \approx {1}/{{M}}$
(Nikias & Petropulu Reference Nikias and Petropulu1993).
The spectral and higher-order spectral quantities introduced above provide a framework for characterising linear and nonlinear interactions among frequency components. Power and cross-power spectra describe energy distribution and linear coupling, while bispectral measures capture triadic phase-coupled interactions responsible for nonlinear energy exchange. These quantities form the basis for constructing the nonlinear energy transfer functions developed in the following section.
3. Methods for estimating nonlinear energy transfer function
Using the spectral framework developed in the previous section, we describe two approaches for estimating nonlinear energy transfer between spectral components: (i) the Ritz method, based on the work of Ritz & Powers (Reference Ritz and Powers1986), Ritz et al. (Reference Ritz, Powers, Miksad and Solis1988), Ritz et al. (Reference Ritz, Powers and Bengtson1989); and (ii) the Kim method, based on the work of Kim et al. (Reference Kim, Durst, Fonck, Fernandez, Ware and Terry1996). The formulation and implementation of these methods are presented in the following subsections.
3.1. Ritz method
Based on the nonlinear energy transfer formalism introduced by Ritz et al. (Reference Ritz, Powers and Bengtson1989) and subsequent studies (Ritz & Powers Reference Ritz and Powers1986; Ritz et al. Reference Ritz, Powers, Miksad and Solis1988), the evolution of turbulent spectral components can be described using a nonlinear wave coupling equation. This formulation separates the dynamics into linear and quadratic interaction terms, providing a natural framework to analyse the direction of nonlinear energy transfer among interacting modes. Although Ritz’s original formulation was developed in wavenumber space, we apply it in the frequency domain using Taylor’s frozen-flow hypothesis (Taylor Reference Taylor1938; Howes, Klein & TenBarge Reference Howes, Klein and TenBarge2014), which enables a direct transformation between wavenumber and frequency (
$k \rightarrow f$
,
$f \rightarrow k$
). This hypothesis assumes that the fluctuations are effectively ‘frozen’ and are advected past the probe by the background flow without undergoing significant evolution during the measurement interval. This hypothesis is expected to be valid when the advection velocity exceeds the intrinsic phase velocity of the fluctuations. Physically, this implies that the temporal variations measured by a stationary probe arise predominantly from the convection of spatial structures rather than from their intrinsic temporal evolution. The applicability of the hypothesis may be reduced when the intrinsic phase velocity becomes comparable to the advection velocity or when the fluctuations evolve significantly during the measurement interval (Howes et al. Reference Howes, Klein and TenBarge2014; Klein, Howes & TenBarge Reference Klein, Howes and TenBarge2014). In the present experiment, the measured
$E\times B$
flow velocity exceeds the phase velocity of the dominant modes at the investigated radial locations, thereby supporting the applicability of Taylor’s hypothesis for the frequency-to-wavenumber conversion employed in the energy transfer analysis. This approach simplifies the experimental procedure: wavenumber-space analysis typically requires multi-probe measurements, whereas frequency-space analysis can be performed with a two-probe set-up. The governing equation in the frequency domain can then be expressed as (Ritz et al. Reference Ritz, Powers, Miksad and Solis1988)
\begin{align} \frac {\partial \phi ({f},x)}{\partial x} = \varLambda _{\kern-1pt f}^{L} \phi ({f},x) + \frac {1}{2} \sum _{\substack {{f}_1,\, {f}_2 \\ {f}_1 + {f}_2 = {f}}} \varLambda _{\kern-1pt f}^{Q}({f}_1, {f}_2)\, \phi ({f}_1,x)\, \phi ({f}_2,x) ,\end{align}
where
$ \phi ({f},x)$
is the spatial Fourier spectrum of the fluctuating field
$ \varphi (t,x)$
defined by
\begin{align} \phi ({f},x) = \sum _{x} \varphi (t,x) {\rm e}^{-i\kern-1pt {f}\kern-1pt t} .\end{align}
Here, we assume that
$\phi (t,x)$
is temporally stationary and has a zero mean value. Temporal stationarity implies that the statistical properties of the signal are invariant in time; in particular, the power associated with each frequency remains unchanged across different time segments. Equation (3.1) represents the nonlinear wave coupling equation governing the evolution of the Fourier spectrum through linear and quadratic nonlinear interactions. These linear and nonlinear processes are characterised by the linear transfer coefficient (
$\varLambda _f^{L}$
) and the quadratic nonlinear coupling coefficient (
$\varLambda _{{f}}^{Q}$
), respectively. The linear coefficient accounts for physical effects such as linear growth (or damping) and dispersion, whereas the quadratic nonlinear coupling coefficient characterises the strength of the three-wave interaction. It governs the transfer of energy among coupled modes through three-wave processes, whereby a parent wave (
$f$
) can decay into two daughter waves (
$f_1$
) and (
$f_2$
), or conversely, two waves can merge to form a single mode. The quadratic coupling coefficient arises from convective nonlinearities, including the nonlinear polarisation drift and the
$\boldsymbol{E}\times \boldsymbol{B}$
advection (Hasegawa & Mima Reference Hasegawa and Mima1978). These nonlinear interactions facilitate the transfer of energy between different modes through wave merging and decay processes, thereby playing a significant role in the spectral evolution of the plasma fluctuations. The linear transfer coefficient is defined as
$ \varLambda _f^{L} = \gamma _f + ik_f$
, with
$ \gamma _f$
and
$ k_f$
being the linear growth rate and the dispersion relation, respectively. Here
$ \varLambda _{{f}}^{Q}$
is the nonlinear coupling coefficient, giving the strength of mode coupling between
$ {f}_1$
and
$ {f}_2$
.
The nonlinear wave coupling equation arises from the retention of convective nonlinearities in the governing dynamical equations. For neutral fluids, this coupling originates from the nonlinear advective term
$(\boldsymbol{v}{\cdot }\boldsymbol{\nabla })\boldsymbol{v}$
in the Navier–Stokes equations, whose Fourier representation naturally leads to triadic mode interactions. In plasma systems, an analogous coupling appears in DW turbulence through the nonlinear polarisation drift and
$\boldsymbol{E}\times \boldsymbol{B}$
convection, and is formally described by the Hasegawa–Mima equation (Hasegawa & Mima Reference Hasegawa and Mima1978). Similar three-wave coupling dynamics have been invoked to model turbulence in diverse physical systems, including shallow-water surface waves (Freilich & Guza Reference Freilich and Guza1984) and atmospheric Rossby waves (Hasegawa et al. Reference Hasegawa, Maclennan and Kodama1979), providing a common framework for describing nonlinear energy transfer across scales.
The method itself consists in computing the fast Fourier transforms,
$ \phi ({f},x)$
and
$ \phi ({f}, x + \Delta x)$
, from the data of two time series measured at a fixed spatial position
$ x$
and at another different position
$ x + \Delta x$
(i.e.
$ \varphi (t,x)$
and
$ \varphi (t,x + \Delta x)$
). We can consider the first signal
$ \phi ({f},x)$
as the input and the space separated signal
$ \phi ({f},x + \Delta x)$
as the output, such that the output signal can be expressed in terms of the input signal via the linear and quadratic transfer functions, in which the linear and nonlinear coupling coefficients are involved.
We start from (3.1) to establish the relation between the output and input signals for a change of the spectrum along
$ x$
. Considering that the Fourier spectrum of
$ \phi (f,x)$
is a complex quantity
$ \phi (f,x) = |\phi (f,x)|\mathrm{e}^{i\varTheta (f,x)}$
, where
$ |\phi (f,x)|$
is the amplitude and
$ \varTheta (f,x)$
is the phase, we can estimate how it evolves along
$ x$
by a difference approach:
\begin{align} \frac {\partial \phi ({f},x)}{\partial x} &= \frac {\partial \,|\phi ({f},x)|}{\partial x}\, \mathrm{e}^{i\varTheta ({f},x)} + |\phi ({f},x)|\,\frac {\partial \mathrm{e}^{i\varTheta ({f},x)}}{\partial x} \notag \\[7pt] &= \frac {\partial \,|\phi ({f},x)|}{\partial x}\, \mathrm{e}^{i\varTheta ({f},x)} + \frac {\partial \,\varTheta ({f},x)}{\partial x}\, i\,\mathrm{e}^{i\varTheta ({f},x)}|\phi ({f},x)|. \end{align}
Using the forward finite difference formulas, i.e.
\begin{align} \begin{aligned} \frac {\partial \,|\phi ({f},x)|}{\partial x} &\approx \frac {|\phi ({f},x+\Delta x)| - |\phi ({f},x)|}{\Delta x}, \\[7pt] \frac {\partial \,\varTheta ({f},x)}{\partial x} &\approx \frac {\varTheta ({f},x+\Delta x) - \varTheta ({f},x)}{\Delta x}, \end{aligned} \end{align}
respectively, in (3.1), and rearranging the terms, we obtain
\begin{align} \phi ({f}, x + \Delta x) &= \frac {\varLambda _{{f}}^{L} \Delta x + 1 - i \left [ \varTheta ({f}, x + \Delta x) - \varTheta ({f}, x) \right ]} {\mathrm{e}^{-i\left [ \varTheta ({f}, x + \Delta x) - \varTheta ({f}, x) \right ]}} \, \phi ({f}, x) \notag \\[7pt] &\quad + \frac {1}{2} \sum _{\substack {{f}_1,\, {f}_2 \\ {f}_1 + {f}_2 = {f}}} \frac {\varLambda _{{f}}^{Q}({f}_1, {f}_2) \, \Delta x} {\mathrm{e}^{-i\left [ \varTheta ({f}, x + \Delta x) - \varTheta ({f}, x) \right ]}} \, \phi ({f}_1, x)\, \phi ({f}_2, x) .\end{align}
If we define the input and output Fourier spectra as
\begin{align} \begin{aligned} X &= \phi ({f}, x) \quad\text{(input)},\\[4pt] Y &= \phi ({f}, x + \Delta x)\quad\text{(output)}, \end{aligned} \end{align}
and the linear and quadratic transfer functions as
\begin{align} L_{\kern-1pt f} &= \frac {\varLambda _{{f}}^{L} \Delta x + 1 - i \left [ \varTheta ({f}, x + \Delta x) - \varTheta ({f}, x) \right ]} {\mathrm{e}^{-i\left [ \varTheta ({f}, x + \Delta x) - \varTheta ({f}, x) \right ]}} ,\\[-10pt]\nonumber\end{align}
\begin{align} Q_{{f}}^{{f}_1 {f}_2} &= \frac {\varLambda _{{f}}^{Q}({f}_1, {f}_2)\, \Delta x} {\mathrm{e}^{-i\left [ \varTheta ({f}, x + \Delta x) - \varTheta ({f}, x) \right ]}} ,\end{align}
then, (3.5) can be written concisely as
\begin{align} Y_{\kern-2pt f} = L_{\kern-1pt f} X_{\kern-1.5pt f} + \sum _{\substack {{f}_1 \geq {f}_2 \\ {f} = {f}_1 + {f}_2}} Q_{{f}}^{{f}_1 {f}_2} X_{{\kern-2pt f}_{\kern-1.5pt 1}}\, X_{{\kern-2pt f}_{\kern-1.5pt 2}} .\end{align}
3.1.1. Estimation of linear and quadratic transfer functions
We now proceed to determine
$L_{\kern-1pt f}$
and
$Q_f^{f_1 f_2}$
, and apply them within the framework of the expressions derived above. To obtain the expressions of
$ L_{\kern-1pt f}$
and
$ Q_{{f}}^{{f}_1 {f}_2}$
, one can multiply (3.9) with the complex conjugate of
$X_{\kern-1.5pt f}$
, and
$ X_{{\kern-2pt f}_{\kern-1.5pt 1}} X_{{\kern-2pt f}_{\kern-1.5pt 2}}$
, respectively, and average over many statistically similar realisations. The results are
\begin{align} \langle Y_{\kern-2pt f} X_{\kern-1.5pt f}^* \rangle = L_{\kern-1pt f} \langle X_{\kern-1.5pt f} X_{\kern-1.5pt f}^* \rangle + \sum _{\substack {{f}_1 \geq {f}_2 \\ {f} = {f}_1 + {f}_2}} Q_{{f}}^{{f}_1, {f}_2} \big\langle X_{{\kern-2pt f}_{\kern-1.5pt 1}} X_{{\kern-2pt f}_{\kern-1.5pt 2}} X_{\kern-1.5pt f}^* \big\rangle,\end{align}
\begin{align} \big\langle Y_{\kern-2pt f} X_{{f}_1^{\prime }}^{*} X_{{f}_2^{\prime }}^{*} \big\rangle = L_{\kern-1pt f} \big\langle X_{\kern-1.5pt f} X_{{f}_1^{\prime }}^{*} X_{{f}_2^{\prime }}^{*} \big\rangle + \sum _{\substack {{f}_1 \geq {f}_2 \\ {f} = {f}_1 + {f}_2}} Q_{{f}}^{{f}_1, {f}_2} \big\langle X_{{\kern-2pt f}_{\kern-1.5pt 1}} X_{{\kern-2pt f}_{\kern-1.5pt 2}} X_{{f}_1^{\prime }}^{*} X_{{f}_2^{\prime }}^{*} \big\rangle ,\end{align}
where
$ {f} = {f}_1 + {f}_2 = {f}_1^{\prime } + {f}_2^{\prime }$
. The notation
$\langle {\cdot }\rangle$
denotes the ensemble average over realisations.
The evaluation of the fourth-order moment
$\langle X_{f_1} X_{f_2} X_{f_1'}^{*} X_{f_2'}^{*} \rangle$
is computationally demanding. To reduce the complexity, the Ritz method employs Millionshchikov’s approximation (Millionshchikov Reference Millionshchikov1941; Ritz & Powers Reference Ritz and Powers1986), whereby contributions from terms satisfying
$(f_1, f_2) \neq (f_1', f_2')$
are neglected. A similar approach has been employed in various theoretical models of strong turbulence (Terry & Diamond Reference Terry and Diamond1985; Terry et al. Reference Terry, Diamond, Shaing, Garcia and Carreras1986). Consequently, (3.11) reduces to
\begin{align} \big\langle Y_{\kern-2pt f} X_{{\kern-2pt f}_{\kern-1.5pt 1}}^{*} X_{{\kern-2pt f}_{\kern-1.5pt 2}}^{*} \big\rangle = L_{\kern-1pt f} \big\langle X_{\kern-1.5pt f} X_{{\kern-2pt f}_{\kern-1.5pt 1}}^{*} X_{{\kern-2pt f}_{\kern-1.5pt 2}}^{*} \big\rangle + Q_{{f}}^{{f}_1 {f}_2} \big\langle |X_{{\kern-2pt f}_{\kern-1.5pt 1}} X_{{\kern-2pt f}_{\kern-1.5pt 2}}|^2 \big\rangle. \end{align}
Operating on this equation, we obtain
\begin{align} Q_{{f}}^{{f}_1 {f}_2} = \frac {\big\langle Y_{\kern-2pt f} X_{{\kern-2pt f}_{\kern-1.5pt 1}}^{*} X_{{\kern-2pt f}_{\kern-1.5pt 2}}^{*} \big\rangle - L_{\kern-1pt f} \big\langle X_{\kern-1.5pt f} X_{{\kern-2pt f}_{\kern-1.5pt 1}}^{*} X_{{\kern-2pt f}_{\kern-1.5pt 2}}^{*} \big\rangle } {\big\langle \big|X_{{\kern-2pt f}_{\kern-1.5pt 1}} X_{{\kern-2pt f}_{\kern-1.5pt 2}}\big|^2 \big\rangle }. \end{align}
Substituting this expression of
$ Q_{{f}}^{{f}_1 {f}_2}$
into (3.12), we get
\begin{align} L_{\kern-1pt f} = \frac { \langle Y_{\kern-2pt f} X_{\kern-1.5pt f}^{*}\rangle - \sum _{\substack {{f}_1 \geq {f}_2 \\ {f} = {f}_1 + {f}_2}} \big({ \big\langle X_{\kern-1.5pt f}^{*} X_{{\kern-2pt f}_{\kern-1.5pt 1}} X_{{\kern-2pt f}_{\kern-1.5pt 2}} \big\rangle \big\langle Y_{\kern-2pt f} X_{{\kern-2pt f}_{\kern-1.5pt 1}}^{*} X_{{\kern-2pt f}_{\kern-1.5pt 2}}^{*} \big\rangle }\big)/\left({\left\langle |X_{{\kern-2pt f}_{\kern-1.5pt 1}} X_{{\kern-2pt f}_{\kern-1.5pt 2}}|^2 \right\rangle }\right) } { \langle X_{\kern-1.5pt f} X_{\kern-1.5pt f}^{*} \rangle - \sum _{\substack {{f}_1 \geq {f}_2 \\ {f} = {f}_1 + {f}_2}} \big({ \big|\big\langle X_{\kern-1.5pt f}^{*} X_{{\kern-2pt f}_{\kern-1.5pt 1}} X_{{\kern-2pt f}_{\kern-1.5pt 2}} \big\rangle \big|^2 }\big)/\big({\big\langle \big|X_{{\kern-2pt f}_{\kern-1.5pt 1}} X_{{\kern-2pt f}_{\kern-1.5pt 2}}\big|^2 \big\rangle }\big) }, \end{align}
which is needed to calculate the linear growth rate (
$ \gamma _{{f}}$
) and dispersion relation
$(k_f)$
. The value of the quadratic transfer function
$ Q_{{f}}^{{f}_1 {f}_2}$
can be obtained by substituting
$ L_{\kern-1pt f}$
into (3.13) together with the ‘averaged products’ of Fourier transforms of the measured signals. In our notation of spectral quantities described in § 2, we can rewrite (3.13) and (3.14) as
\begin{align} L_{\kern-1pt f} = \frac { P_c(f) - \sum _{\substack {{f}_1 \geq {f}_2 \\ {f} = {f}_1 + {f}_2}} ({ B_a(f_1, f_2) B_c^*(f_1, f_2) })/({E_{12}}) } { P_a(f) - \sum _{\substack {{f}_1 \geq {f}_2 \\ {f} = {f}_1 + {f}_2}} ({ |B_a(f_1, f_2)|^2 })/({E_{12}}) } ,\end{align}
\begin{align} Q_{{f}}^{{f}_1 {f}_2} = \frac {B_c^*(f_1, f_2) - L_{\kern-1pt f} B_a^*(f_1, f_2)} {E_{12}} .\end{align}
Here, we set
$E_{12} = \langle |X_{{\kern-2pt f}_{\kern-1.5pt 1}} X_{{\kern-2pt f}_{\kern-1.5pt 2}}|^2 \rangle$
. These transfer function definitions are used to estimate the corresponding quantities in the following sections.
3.1.2. Estimation of power and energy transfer functions
One can write a wave kinetic equation for the power spectrum
$ P_{{f}} = \langle X_{\kern-1.5pt f} X_{\kern-1.5pt f}^* \rangle$
in terms of the coupling coefficients
$ \varLambda _{{f}}^L$
and
$ \varLambda _{{f}}^Q$
, by multiplying (3.1) with
$ \phi ^*({f}, x)$
(Ritz et al. Reference Ritz, Powers and Bengtson1989):
\begin{align} \frac {\partial \left [ \phi ({f}, x)\, \phi ^*({f}, x) \right ]}{\partial x} = \frac {\partial \phi ({f}, x)}{\partial x} \, \phi ^*({f}, x) + \frac {\partial \phi ^*({f}, x)}{\partial x} \, \phi ({f}, x) .\end{align}
Using (3.1) in (3.17) and ensemble averaging leads to the following wave kinetic equation (Kim et al. Reference Kim, Durst, Fonck, Fernandez, Ware and Terry1996), i.e.
\begin{align} \frac {\partial P_{{f}}}{\partial x} = 2\gamma _{{f}} P_{{f}} + T_f, \quad \text{where}\quad T_f = \sum _{\substack {{f}_1 , {f}_2 \\ {f} = {f}_1 + {f}_2}} T_{{f}}({f}_1, {f}_2) .\end{align}
Here the power transfer function is written as
\begin{align} T_{{f}}({f}_1, {f}_2) = \mathrm{Re} \left [ \varLambda _{{f}}^{Q}({f}_1, {f}_2)\, \langle X_{\kern-1.5pt f}^{*} X_{{\kern-2pt f}_{\kern-1.5pt 1}} X_{{\kern-2pt f}_{\kern-1.5pt 2}} \rangle \right ] = \mathrm{Re} \left [ \varLambda _{{f}}^{Q}({f}_1, {f}_2)\, B_a(f_1, f_2) \right ], \end{align}
where
$\varLambda _{{f}}^{Q}({f}_1, {f}_2)$
can be obtained from (3.8).
Once the power transfer function is calculated then the energy transfer function can be easily calculated as
\begin{align} W^{f}_{\mathrm{NL}} = T_f = \sum _{\substack {f_1, f_2 \\ f = f_1 + f_2}} T_{f}(f_1, f_2) .\end{align}
To check the reliability of the estimated nonlinear energy transfer, we use spectral energy conservation as a physical consistency condition, i.e.
$\sum _f W_{\mathrm{NL}}^{f} = 0$
. Any deviation from this relation indicates an imbalance in the evaluated transfer. To quantify this imbalance, we define the energy mismatch parameter
$ W_{\mathrm{mis}}$
(Xia & Shats Reference Xia and Shats2004) as
\begin{align} W_{\mathrm{mis}} = \frac {\sum _f W_{\mathrm{NL}}^{f}}{\sum _f \lvert W_{\mathrm{NL}}^{f} \rvert }. \end{align}
This parameter measures how far the computed nonlinear energy transfer departs from exact conservation. For a perfectly balanced transfer,
$ W_{\mathrm{mis}} = 0$
, while non-zero values indicate residual imbalance in the spectral energy budget.
3.1.3. Computational algorithm for Ritz method
In this subsection we outline the computational procedure employed to estimate the linear and quadratic transfer functions using the Ritz method. The algorithm is implemented in the frequency domain and is based on ensemble-averaged spectral quantities obtained from multiple realisations of the input–output data. The steps of the algorithm are as follows.
-
(i) Step 1 (original time-series data): the input and output signals are given as time-domain sequences,
$ X(t)$
and
$ Y(t)$
, respectively. -
(ii) Step 2 (pre-processing).
-
(a) Step 2.1 (ensemble and data points): define the number of ensemble realisations (M) and the number of data points per realisation (N).
-
(b) Step 2.2 (mean removal and detrending): apply mean removal or detrending to the time-series data as appropriate. Mean removal is performed when the data exhibit a DC offset, whereas detrending is applied in the presence of low-frequency trends.
-
(c) Step 2.3 (overlapping and window function): No overlapping or window functions are applied to the time-series data. This is because highly correlated samples do not significantly improve statistical convergence when sufficient independent realisations are available. Furthermore, windowing can distort phase relationships and bias higher-order statistics.
-
(d) Step 2.4 (fast Fourier transform): apply the fast Fourier transform to the input and output time-series data to obtain their frequency-domain representations. The resulting spectral components of the input and output signals are denoted by
$ X_i(f)$
and
$ Y_i(f)$
, respectively, for
$ i = 1, 2, \ldots , M$
.
-
-
(iii) Step 3 (computation of valid frequency pairs and second- and third-order cumulants): for each frequency index
$ f = - {N}/{2}, \ldots , {N}/{2}-1$
, all frequency pairs
$ (f_1, f_2)$
satisfying the resonance condition
$ f_1 + f_2 = f$
are calculated. Using these valid triads, the relevant statistical quantities are evaluated, including the power spectral densities
$ P_a(f)$
and
$ P_c(f)$
, the bispectra
$ B_a(f_1, f_2)$
and
$ B_c(f_1, f_2)$
, and the corresponding bicoherences
$ b_a(f_1, f_2)$
and
$ b_c(f_1, f_2)$
. -
(iv) Step 4 (computation of transfer functions): calculate the linear and quadratic transfer functions, i.e.
$ L_{\kern-1pt f}$
and
$ Q_f^{f_1, f_2}$
. -
(v) Step 5 (computation of coupling coefficients): calculate the linear and quadratic coupling coefficients, i.e.
$ \varLambda _f^L$
and
$ \varLambda _f^Q(f_1, f_2)$
, respectively, by multiplying the phase difference between input and output data to the transfer functions. -
(vi) Step 6 (computation of power and energy transfer functions): calculate the power transfer function
$ (T_f(f_1, f_2))$
and energy transfer function (
$ W_{NL}^f$
) by summing
$ T_f(f_1, f_2)$
along the
$ f = f_1 + f_2$
axis. -
(vii) Step 7 (outputting results): output the following quantities:
$L_{\kern-1pt f}, Q_f^{f_1, f_2},\varLambda _f^L,$
$\varLambda _f^Q(f_1, f_2), T_f(f_1, f_2), W_{NL}^f$
.
We use the above algorithm for the Ritz method in case of both simulated data and experimental data. However, applying the Ritz method to experimental data requires a data selection process that is discussed in the following sections.
3.1.4. Simulation test and analysis
To verify the validity of the method described in the previous section, we performed a computer simulation. First, we define the analytical form of linear and quadratic transfer functions, then we generate an output signal from a given Gaussian input. Next, we estimate the linear and quadratic transfer functions using the method described in § 3.1.1, and then we compare the results of the estimated transfer function with the analytical ones.
3.1.4.1. Simulation data generation
To carry out the simulation test, we defined a ‘nonlinear black-box model’ in which we can generate a non-Gaussian output from a Gaussian input signal following (3.9). This model consists of a linear transfer function (
$L_{\kern-1pt f}$
) that makes the input and output spectra similar in shape as expected from a stationary state and a quadratic transfer function (
$Q_{f}^{f_1,f_2}$
), which incorporates the nonlinear coupling among the spectral components of the input signal. The schematic diagram of this nonlinear black-box model is shown in figure 1.
Nonlinear black-box model.
The analytical linear and quadratic transfer functions are defined as (Ritz & Powers Reference Ritz and Powers1986)
\begin{align} L_{\kern-1pt f} = 1.0 - 0.4\frac {f^2}{f_{\text{Nyq}}^2} + i \, 0.8 \frac {f}{f_{\text{Nyq}}}, \end{align}
\begin{align} Q_{f}^{f_1,f_2} = \frac {i}{5 f_{\text{Nyq}}^4} \frac { f_1 f_2 \left ( f_2^2 - f_1^2 \right )}{1 + \dfrac {f^2}{f_{\text{Nyq}}^2}} , \end{align}
where
$f= f_1 + f_2$
and
$f_{\text{Nyq}}$
is the Nyquist frequency.
The choice of these functions is arbitrary, but their values are of the same order of magnitude as those predicted by the Hasegawa–Mima equation (Hasegawa & Mima Reference Hasegawa and Mima1978). These analytical transfer functions closely resemble those reported by Shen et al. (Reference Shen, Li, Li and Li2020) for a nonlinear developing turbulence model, in which the input and output power spectra differ, a feature also observed in shallow-water ocean waves. Now, to carry out a simulation, we need to generate nonlinearity in our input and output signal as we are looking for the nonlinear energy transfer function. A white Gaussian noise as an input signal shows zero bispectrum indicating there is no nonlinearity in the signal. The nonlinearity can be introduced to the signal by passing it through a ‘nonlinear black box’ of the type shown in figure 1. Here we have taken zero mean white Gaussian noise as an input signal and passed it through a series of five identical ‘nonlinear black boxes’ that are connected in series. So, in this case a Gaussian white noise signal is fed into the first black box. Due to the nonlinear behaviour of this system, the output is no longer Gaussian and is passed as the input to the next black box. This process is repeated sequentially and, for the simulation, we consider the input–output pair of the fifth black box as our final input (
$X_f$
) and output (
$Y_f$
) in frequency domain. We can always do an inverse Fourier transform to get the input and output time-series signals.
As we have time-series signals, we can now study the power spectrum, bispectrum and bicoherence, which are already defined in § 2. The auto-power spectra
$(P_a(f))$
of the input and output signals after the fifth iteration (i.e. after the fifth nonlinear black box) are shown in figure 2.
Power spectrum of input and output signals after the fifth iteration.
From the power spectra shown in figure 2, it is evident that the input and output spectra do not completely coincide. This mismatch is more pronounced at higher frequencies, a feature commonly observed in shallow-water surface waves in the ocean (Freilich & Guza Reference Freilich and Guza1984), and it is indicative of a developing turbulence regime (Shen et al. Reference Shen, Li, Li and Li2020).
The auto-bispectrum
$(B_a(f_1, f_2))$
and cross-bispectrum
$(B_c(f_1, f_2))$
for the input signal after the fifth iteration is shown in figure 3.
(a) Auto-bispectrum, and (b) cross-bispectrum after the fifth iteration.
3.1.4.2. Convergence test
In the Ritz method, obtaining reliable estimates of the linear and quadratic transfer functions requires a sufficiently large number of ensembles. These quantities are constructed from auto- and cross-power spectra as well as higher-order cumulants, namely the auto- and cross-bispectra, which are known to converge slowly with respect to the number of realisations. When only a small number of ensembles is used, the estimates are dominated by statistical fluctuations, resulting in large variance and poor accuracy. Increasing the ensemble size reduces these fluctuations through averaging, allowing the transfer functions to converge toward their true statistical values. To quantify convergence, we examine the evolution of the auto- and cross-power spectra and the auto- and cross-bispectra as functions of the number of ensembles. The results are considered to be convergent when the relative change of these quantities between successive ensemble sizes falls below 5 %. This criterion ensures that both second- and third-order statistics are sufficiently stable for subsequent analysis of the linear and quadratic transfer functions. The corresponding convergence behaviour is illustrated in figure 4.
Convergence plot of (a) auto- and cross-power spectra, (b) auto- and cross-bispectra.
From the above plots, we observe that the amplitudes of the auto- and cross-power spectra as well as the auto- and cross-bispectra exhibit a pronounced oscillatory behaviour for small ensemble sizes
$M$
, indicating strong statistical fluctuations when only a limited number of realisations are available. As
$M$
increases, these fluctuations are progressively suppressed and the estimates approach stable values. Beyond
$M \approx 2000$
, the relative change in these quantities reduces to about
$0.5\,\%$
, which is already well below the adopted convergence threshold of
$5\,\%$
. Furthermore, increasing the ensemble size improves the statistical accuracy even more. As seen in figure 4, the relative change falls below
$0.1\,\%$
once
$M \gtrsim 3500$
, indicating near-complete convergence of both second- and third-order statistics. For the present simulations, we therefore use
$M = 2500$
realisations of a Gaussian signal, each consisting of
$N = 64$
points, which is sufficient to obtain reliable estimates of the linear and quadratic transfer functions. However, for experimental data, where additional noise and non-idealities are present, we employ ensemble sizes larger than
$3500$
to ensure higher statistical robustness and a more reliable transfer function estimation.
3.1.4.3. Simulation validation
To validate our formulation, we benchmark the estimated linear and quadratic transfer functions against their corresponding analytical expressions. Figures 5 and 6 shows the comparison for
$L_{\kern-1pt f}$
and
$Q^{f_1,f_2}_f$
with their respective analytical functions.
Benchmarking of
$L_{\kern-1pt f}$
with analytical function.
Benchmarking of
$Q_f^{f_1,f_2}$
with analytical function.
As discussed earlier, we use the input and output data at iteration
$5$
, for which the estimated transfer functions show good agreement with the analytical results. A small deviation is observed in the linear transfer function near zero frequency, and the contour plot of the estimated quadratic transfer function appears slightly noisier than its analytical counterpart. This behaviour primarily arises from the use of the Millionshchikov’s approximation in the present method, which affects the estimation of higher-order statistics. Such limitations do not arise in Kim’s method, where this approximation is not used, as discussed in the subsequent sections.
3.1.5. Applicability of Ritz method
In the previous section we discussed that if we feed a white Gaussian noise at the input of the nonlinear black box then it will gain some nonlinearity at the output due to the nonlinear nature of the black box. This signature of the nonlinearity can be identified by the non-zero bispectrum of the signal. On the other hand, we can identify the nonlinearity by looking at the distribution function of the signal or, more specifically, the nonlinearity can be identified by skewness and kurtosis of the signal’s distribution function. The skewness of a signal
$\{x_i\}_{i=1}^N$
is defined as (Bendat & Piersol Reference Bendat and Piersol2010; Spiegel & Stephens Reference Spiegel and Stephens2018)
\begin{align} \text{Skewness} = \frac {1}{N}\sum _{i=1}^N \left ( \frac {x_i - \mu }{\sigma } \right )^3, \end{align}
where
$\mu = ({1}/{N})\sum _{i=1}^N x_i$
is the mean and
$\sigma$
is the standard deviation of the signal. Skewness corresponds to the third-order central moment and it measures the asymmetry of the probability distribution. A positive skewness indicates a distribution with a longer tail on the right, while a negative skewness indicates a longer tail on the left. For a Gaussian distribution, skewness equals 0 (Bendat & Piersol Reference Bendat and Piersol2010).
The kurtosis of the signal is defined as (Bendat & Piersol Reference Bendat and Piersol2010; Spiegel & Stephens Reference Spiegel and Stephens2018)
\begin{align} \text{Kurtosis} = \frac {1}{N}\sum _{i=1}^N \left ( \frac {x_i - \mu }{\sigma } \right )^4, \end{align}
which is the fourth-order central moment, measuring the peakedness and tail heaviness of the distribution. A higher kurtosis indicates heavier tails and sharper peak compared with a Gaussian distribution, while lower kurtosis indicates lighter tails and a flatter peak. For a Gaussian distribution, kurtosis equals 3. The excess kurtosis is defined as the kurtosis measured relative to its Gaussian benchmark value of 3 (Bendat & Piersol Reference Bendat and Piersol2010).
Skewness, being a third-order statistic, is directly related to the bispectrum (Ritz & Powers Reference Ritz and Powers1986), which is the Fourier transform of the third-order cumulant. On the other hand, kurtosis, as a fourth-order statistic, is connected to fourth-order cumulants, thereby providing insight into higher-order correlations. For purely Gaussian signals, all higher-order cumulants beyond the second order vanish; hence, any non-zero skewness or excess kurtosis serves as a clear indicator of non-Gaussianity and nonlinearity in the signal. Thus, skewness and kurtosis serve as compact statistical measures of nonlinearity and deviations from Gaussian behaviour in signal data, complementing the more detailed spectral tools such as bispectrum and the fourth-order cumulants.
In the previous section we discussed the nonlinear black-box model in which a Gaussian signal was passed through a number of black boxes connected in series. As an effect of this, we obtained some amount of nonlinearity in the output signal at every iteration. Here, we can observe the signature of nonlinearity by calculating skewness and kurtosis of the output signal at every iteration. The distribution functions before the first iteration (i.e. first black box) and after the fifth iteration (i.e. fifth black box) along with the skewness and kurtosis of the distribution are shown in figure 7.
(a) Probability distribution of input time series having
$\mu = 0.001$
,
$\sigma = 1.00$
, skewness
$=0.001$
, excess kurtosis
$=0.01$
. (b) Probability distribution of output time series having
$\mu = 0.002$
,
$\sigma = 0.52$
, skewness
$= 0.15$
, excess kurtosis
$=0.40$
at the fifth iteration.
From the figure 7 we observe that before the first iteration (i.e. first black box) the distribution of the input time series (Gaussian white noise) gives the expected mean and standard deviation, but we find small non-zero skewness and kurtosis values are still there that are nothing but the consequence of not taking a long signal length. On the other hand, the distribution at the output of the fifth iteration (i.e. the fifth black box) exhibits a noticeable increase in skewness and excess kurtosis, indicating that substantial nonlinearity has been introduced into the signal after passing through the fifth nonlinear black box (i.e. iteration
$5$
). By observing the output across successive iterations, we find that both skewness and excess kurtosis increase with the number of iterations. This trend is illustrated in figure 8(a), where results up to 10 iterations are shown.
Evolution of statistical parameters with nonlinear iterations to examine the validity of Millionshchikov’s approximation.
From figure 8(a) it is clear that kurtosis increases more rapidly than skewness with iteration. So, we see that with the increase of nonlinearity, skewness and kurtosis of the distribution increases – here comes the main issue of applicability of the Ritz method. In the Ritz method a key assumption is that the fourth-order statistical moments of the system can be approximated by the square of the second-order moments. This assumption is known as Millionshchikov’s approximation (Millionshchikov Reference Millionshchikov1941; Ritz & Powers Reference Ritz and Powers1986). However, this approximation has a limitation. The quantity called kurtosis is a direct statistical measure of the fourth-order central moment of a distribution. When kurtosis remains close to its Gaussian value, the approximation works reasonably well, but as kurtosis increases significantly, the higher-order statistics deviate strongly from Gaussian behaviour. In such cases, the fourth-order moment can no longer be represented simply as the square of the second-order moment. Furthermore, the inapplicability of Millionshchikov’s approximation is also supported by figure 8(b), where the square of the standard deviation (
$\sigma ^2$
) and the excess kurtosis are plotted as a function of iteration number. Since the standard deviation (
$\sigma$
) is related to the second-order central moment, the fourth-order moments can be approximated by the square of the second-order moments only when
$\sigma ^2$
remains comparable to the excess kurtosis, i.e. when Millionshchikov’s approximation is valid. From figure 8(b) it can be seen that up to iteration
$5$
the excess kurtosis and
$\sigma ^2$
remain relatively close. However, beyond this point, the excess kurtosis increases to significantly higher values, resulting in a substantial deviation between the two quantities and indicating the breakdown of Millionshchikov’s approximation. Therefore, when the kurtosis becomes large, the fundamental assumption behind the Ritz method breaks down, and the method becomes invalid. This marks the key limitation in the applicability of the Ritz method. Hence, the Ritz method cannot be applied to highly nonlinear signals.
In this study, using artificial data, we observed that the Ritz method fails to correctly estimate the transfer functions once the excess kurtosis exceeds approximately 0.40, which occurs after five iterations. As shown earlier, both the linear and quadratic transfer functions exhibit good agreement with the analytical functions at iteration
$5$
in figures 5 and 6. However, at iteration 7, where the kurtosis reaches approximately 0.55, the transfer function plots exhibit clear deviations from the analytical functions, as shown in figure 9 and figure 10. Figures 9 and 10 show a significant deviation of the estimated transfer functions from the analytical ones. These results are also consistent with figure 8(b), which shows that Millionshchikov’s approximation breaks down beyond iteration
$5$
. Consequently, the assumptions underlying the Ritz method are no longer satisfied, indicating that the Ritz method is not applicable for iterations greater than 5. In the next section we discuss the Kim method and validate it using simulation data, as done for the Ritz method, and also examine its performance at higher iterations.
Plot of
$L_{\kern-1pt f}$
with analytical function at iteration 7.
Plot of
$Q_f^{f_1,f_2}$
with analytical function at iteration 7.
3.2. Kim method
As discussed earlier, the Ritz method suffers from several inherent limitations. One major source of inaccuracy arises from the use of the Millionshchikov approximation (Millionshchikov Reference Millionshchikov1941), in which fourth-order moments are approximated by the square of second-order moments. As a consequence, the Ritz method loses accuracy for fluctuation data exhibiting significant non-Gaussianity or high kurtosis, as demonstrated in the previous section. In addition, the application of the Ritz method to experimentally measured fluctuations often results in unrealistically large growth or damping coefficients across the entire spectrum. According to Kim et al. (Reference Kim, Durst, Fonck, Fernandez, Ware and Terry1996), these unphysical results originate from the fact that the Ritz method implicitly assumes that all measured fluctuations satisfy the nonlinear wave coupling equation (3.1). In practice, experimental signals contain non-ideal fluctuations that do not obey this equation, such as contributions from systematic errors, higher-order nonlinear interactions beyond quadratic coupling and instrumental or environmental noise. To overcome these limitations, Kim et al. (Reference Kim, Durst, Fonck, Fernandez, Ware and Terry1996) proposed a modified formalism in which the measured spectra
$(X_f, Y_f)$
are decomposed into an ideal component
$(\beta _f, \alpha _f)$
that participates in the linear and quadratic interaction processes described by (3.1), and a non-ideal component
$(X_f^{ni}, Y_f^{ni})$
that does not contribute to these processes. That is,
\begin{align} X_f = \beta _f + X_f^{\mathrm{ni}}, \qquad Y_f = \alpha _f + Y_f^{\mathrm{ni}}. \end{align}
Furthermore, unlike the Ritz method, Kim et al. (Reference Kim, Durst, Fonck, Fernandez, Ware and Terry1996), avoided the Millionshchikov’s approximation and explicitly retained the fourth-order moments in the calculation of the nonlinear transfer functions, leading to a more accurate and physically consistent description of nonlinear energy transfer in experimental turbulence data. Using (3.26) in (3.9), and expanding the wave coupling equation into a hierarchy of moment equations by multiplying by
$X_f^{}$
,
$Y_f^{}$
and
$X_{f_1'}^{} X_{f_2'}^{}$
, we obtain
\begin{align} \left . \begin{aligned} \langle Y_f X_f^{*} \rangle &= L_{\kern-1pt f} \langle \beta _f \beta _f^{*} \rangle + \sum _{\substack{f_1 \ge f_2 \\ f = f_1 + f_2}} Q_{f}^{f_1,f_2}\, \big\langle X_{f_1} X_{f_2} X_f^{*} \big\rangle , \\[4pt] \langle \alpha _f \alpha _f^{*} \rangle &= L_{\kern-1pt f} \langle X_f Y_f^{*} \rangle + \sum _{\substack{f_1 \ge f_2 \\ f = f_1 + f_2}} Q_{f}^{f_1,f_2}\, \big\langle X_{f_1} X_{f_2} Y_f^{*} \big\rangle , \\[4pt] \big\langle Y_f X_{f_1}^{*} X_{f_2}^{*} \big\rangle &= L_{\kern-1pt f} \big\langle X_f X_{f_1}^{*} X_{f_2}^{*} \big\rangle + \sum _{\substack{f_1' \ge f_2' \\ f = f_1' + f_2' = f_1 + f_2}} Q_{f}^{f_1',f_2'}\, \big\langle X_{f_1'} X_{f_2'} X_{f_1}^{*} X_{f_2}^{*} \big\rangle .\end{aligned} \right \} \end{align}
Here
$\langle {\cdot }\rangle$
denotes ensemble averaging. In contrast to the Millionshchikov’s approximation, the fourth-order cumulant
$\langle X_{f_1'} X_{f_2'} X_{f_1}^{} X_{f_2}^{} \rangle$
is explicitly retained in the moment equation (3.27). This results in a more general nonlinear energy transfer formalism. Moreover, retaining the fourth-order term improves the accuracy of the energy transfer estimates, albeit at the cost of increased computational effort. Inspection of the resulting moment equations shows that there are four unknown quantities, namely,
$L_{\kern-1pt f}$
,
$Q_{f}^{f_1,f_2}$
,
$\langle \beta _f \beta _f^{} \rangle$
and
$\langle \alpha _f \alpha _f^{} \rangle$
, while only three independent moment equations are available.
To close the system, Kim et al. (Reference Kim, Durst, Fonck, Fernandez, Ware and Terry1996) imposed the condition that the power spectrum is fully saturated by linear and three-wave coupling processes. This assumption corresponds to spatial stationarity of the turbulence and leads to the constraint
\begin{align} \frac {\partial P_f}{\partial x} \approx \frac {\langle \alpha _f \alpha _f^{*} \rangle - \langle \beta _f \beta _f^{*} \rangle }{\Delta x} = 2\gamma _f P_f + T_f = 0 , \end{align}
where
$P_f$
is the ideal power spectrum,
$\gamma _f$
is the linear growth (or damping) rate and
$T_f$
denotes the nonlinear energy transfer term. From (3.28), the fourth closure relation immediately follows as
\begin{align} \langle \alpha _f \alpha _f^{*} \rangle = \langle \beta _f \beta _f^{*} \rangle . \end{align}
Consequently, the set of moment equations (3.27) can be written in matrix form, and from that we can easily calculate the linear and quadratic transfer functions as
\begin{align} L_{\kern-1pt f} = \frac { \langle Y_f X_f^{*} \rangle - (\boldsymbol{B}^{*})^{T} \cdot \boldsymbol{F}^{-1} \cdot \boldsymbol{A} }{ \langle \beta _f \beta _f^{*} \rangle - (\boldsymbol{A}^{*})^{T} \cdot \boldsymbol{F}^{-1} \cdot \boldsymbol{A} }, \end{align}
\begin{align} \boldsymbol{Q}_f = \bigl [ (\boldsymbol{B}^{*})^{T} - L_{\kern-1pt f} (\boldsymbol{A}^{*})^{T} \bigr ] {\cdot }\boldsymbol{F}^{-1}, \end{align}
where
$\boldsymbol{A}$
and
$\boldsymbol{B}$
denote vectors composed of third-order moments, corresponding to the auto- and cross-bispectra, respectively, while
$\boldsymbol{F}$
represents the matrix of fourth-order moments (Kim et al. Reference Kim, Durst, Fonck, Fernandez, Ware and Terry1996). Here the expression for the power transfer function
$T_f(f_1,f_2)$
is identical to that used in the Ritz method, since both approaches are based on the same wave kinetic equation (3.18). Consequently, the power transfer function in Kim’s method is given by
\begin{align} T_{{f}}({f}_1, {f}_2) = \mathrm{Re} \left [ \varLambda _{{f}}^{Q}({f}_1, {f}_2)\, \langle X_{\kern-1.5pt f}^{*} X_{{\kern-2pt f}_{\kern-1.5pt 1}} X_{{\kern-2pt f}_{\kern-1.5pt 2}} \rangle \right ] = \mathrm{Re} \left [ \varLambda _{{f}}^{Q}({f}_1, {f}_2)\,{\cdot }A \right ].\end{align}
Here, the expression for
$\varLambda _f^{Q}(f_1,f_2)$
follows from (3.8), provided that (3.31) is used for the quadratic transfer function estimation. Also, from the above definition of the power transfer function, we can easily estimate the nonlinear energy transfer function (
$W_{\text{NL}}^f$
) from (3.20).
Furthermore, the Kim method allows for a decomposition of the measured fluctuations into ideal and non-ideal components. The non-ideal component accounts for contributions arising from experimental noise, measurement uncertainties and interactions with physical processes that lie outside the applicability of the underlying model. In the present work, this decomposition is not performed and all observed fluctuations are treated as ideal. With this prescription, we use
$ X_f$
and
$Y_f$
instead of
$\beta$
and
$\alpha$
, respectively, in (3.27) and (3.29) and subsequently obtain the linear transfer function as
\begin{align} L_{\kern-1pt f} = \frac { \langle Y_f X_f^{*} \rangle - (\boldsymbol{B}^{*})^{T} \cdot \boldsymbol{F}^{-1} \cdot \boldsymbol{A} }{ \langle X_f X_f^{*} \rangle - (\boldsymbol{A}^{*})^{T} \cdot \boldsymbol{F}^{-1} \cdot \boldsymbol{A} } \end{align}
and the quadratic transfer function follows accordingly from (3.31). In this treatment, we can write (3.28) for the spatial stationarity
$(\text{i.e.}\,({\partial P_f}/{\partial x})=0)$
as
\begin{align} \frac {\langle Y_f Y_f^{*} \rangle - \langle X_f X_f^{*} \rangle }{\Delta x} = 2\gamma _f P_f + T_f = 0 .\end{align}
Here,
$P_f$
represents the total power associated with the frequency
$f$
. Furthermore, the closure relation reduces to
\begin{align} \langle X_f X_f^{*} \rangle = \langle Y_f Y_f^{*} \rangle .\end{align}
In this formulation, the linear growth rate is expressed as
$\gamma _f = \mathrm{Re} (\varLambda _f^L)$
, with the corresponding linear transfer coefficient
$(\varLambda _f^L)$
extracted from the linear transfer function
$(L_{\kern-1pt f})$
defined in (3.7). However, the power transfer function and the energy transfer function follow the same form as in (3.32) and (3.20).The formulation developed above is applicable to experimental fluctuation data; however, its reliability must first be established. Accordingly, we validate the method using numerical simulations with analytically prescribed transfer functions, following an approach similar to that employed for the Ritz method. The details of this validation procedure and the corresponding results are presented in the subsequent section.
3.2.1. Computational algorithm for Kim method
For the Kim method, the initial computational process is the same as that of the Ritz method, i.e. we follow the same procedures from step 1 to step 2. From step 3 onward, the procedure is modified and presented as follows.
-
(i) Step 3: for
$f = - {N}/{2}$
to
$ {N}/{2}-1$
, execute the following steps.-
(a) Step 3.1: identify all frequency pairs
$(f_1, f_2)$
that satisfy
$f_1 + f_2 = f$
, and store the number of such pairs in Npairs. -
(b) Step 3.2: compute the arrays
$A_f$
and
$B_f$
, and construct the matrix
$F_f$
. -
(c) Step 3.3: evaluate the linear and quadratic transfer functions,
$L_{\kern-1pt f}$
and
$Q_f$
. -
(d) Step 3.4: calculate the linear transfer coefficient and store it in the array
$\varLambda _f^L(f)$
. Compute the quadratic coupling coefficients for all pairs
$(1 \ldots$
Npairs) and assign them to the array
$\varLambda _f^Q$
.
-
Subsequently, to compute the power transfer function
$T_f(f_1, f_2)$
and the nonlinear energy transfer function
$W_{\text{NL}}^f$
, steps 6 and 7 from the Ritz method are used. We use the above algorithm for the Kim method in case of both simulated data and experimental data, which is discussed in the subsequent sections.
3.2.2. Simulation test and analysis
To validate the method, we employ the same analytical linear and quadratic transfer functions introduced in § 3.1.4 ((3.22) and (3.23)). The input and output data are generated using the same nonlinear black-box model adopted in the Ritz-method validation. For the present simulations, the input–output pair from the fifth black box (i.e. iteration
$5$
) is taken as the final input
$X_f$
and output
$Y_f$
in the frequency domain. The simulations consist of
$M = 2500$
realisations, each containing
$N = 64$
data points.
To benchmark the implementation, the estimated linear coefficient
$L_{\kern-1pt f}$
and nonlinear coupling coefficient
$Q_f^{f_1,f_2}$
are compared with their corresponding analytical expressions, and the results are presented in the following figures. From figure 11 it is evident that the estimated transfer functions closely match their corresponding analytical expressions, thereby confirming the correct implementation of the Kim method. Furthermore, the results obtained using the Kim method are noticeably smoother and more accurate than those obtained from the Ritz method, as demonstrated by the benchmarking plots. This improvement arises from the explicit retention of fourth-order moments in the Kim formulation. In the following section we therefore employ the Kim method to address the limitations encountered in the Ritz approach.
Simulation validation of the Kim method. Panels (a,b) show the real and imaginary parts of the linear transfer function
$L_{\kern-1pt f}$
. Panels (c,d) compare the analytical and estimated magnitudes of the quadratic coupling coefficient
$\lvert Q_f(f_1,f_2) \rvert$
.
Figure 11. Long description
Panel A: A line graph shows the real part of the linear transfer function with frequency (Hz) on the x-axis and the real part on the y-axis. The true and estimated values are plotted in blue and red, respectively. Panel B: A line graph shows the imaginary part of the linear transfer function with frequency (Hz) on the x-axis and the imaginary part on the y-axis. The true and estimated values are plotted in blue and red, respectively. Panel C: A heatmap shows the analytical magnitude of the quadratic coupling coefficient with frequency f1 (Hz) on the x-axis and frequency f2 (Hz) on the y-axis. The color scale indicates the magnitude values. Panel D: A heatmap shows the estimated magnitude of the quadratic coupling coefficient with frequency f1 (Hz) on the x-axis and frequency f2 (Hz) on the y-axis. The color scale indicates the magnitude values.
Panels (a,b) show the real and imaginary parts of the linear transfer function
$L_{\kern-1pt f}$
. Panels (c,d) compare the analytical and estimated magnitudes of the quadratic coupling coefficient
$\lvert Q_f(f_1,f_2) \rvert$
at iteration 10.
3.2.3. Addressing the limitations of the Ritz method through Kim method
The primary limitation of the Ritz method is its inability to adequately handle data exhibiting high kurtosis. As discussed in § 3.1.5, kurtosis is a direct statistical measure of the fourth-order central moment of a distribution. However, the Ritz method relies on Millionshchikov’s approximation (Millionshchikov Reference Millionshchikov1941), in which the fourth-order cumulant is replaced by the square of second-order cumulants. This approximation becomes invalid for strongly non-Gaussian or highly nonlinear data.
In contrast, the Kim method does not invoke this approximation and explicitly retains the fourth-order moments, enabling it to more accurately capture nonlinear interactions and to remain robust for data with high kurtosis. This can be easily justified from the benchmarking plots of
$L_{\kern-1pt f}$
and
$Q_f^{f_1,f_2}$
at higher iterations. The plots for 10 iterations is shown in figure 12.
From figure 12 it is evident that the Kim method continues to yield results in good agreement with the analytical solutions even at iteration 10, corresponding to a kurtosis value of 1.04. In contrast, the Ritz method shows noticeable deviations beginning at iteration 7, where the kurtosis reaches approximately 0.55. These observations indicate that the Kim method is more robust to increasing non-Gaussianity and is therefore more suitable than the Ritz method for the analysis of general nonlinear data. In the analysis of experimental fluctuation data, we primarily employ the Kim method due to its robustness in handling strong nonlinearities. However, in the present work, both the Kim and Ritz methods are applied to the same experimental dataset in order to quantitatively assess the deviation between the two approaches when applied to real measurements. This comparison and its implications are discussed in detail in the next section.
3.2.4. Discussion on growth rate and energy transfer function
As discussed earlier, we can estimate the linear growth rate
$\gamma _f$
and the nonlinear energy transfer function
$T_f$
by using the Kim method. These two quantities describe different but complementary aspects of fluctuation dynamics in frequency space. The growth rate
$\gamma _f$
indicates whether the fluctuation at frequency
$f$
grows or decays in time. A positive value
$\gamma _f\gt 0$
means that the mode at frequency
$f$
is unstable and its amplitude increases with time. In this case, the mode extracts energy from the system and acts as a source of energy. A negative value
$\gamma _f\lt 0$
means that the mode is stable and damped, so its amplitude decreases and the mode loses energy unless it is supported by nonlinear interactions with other modes. However, the energy transfer function
$T_f$
describes how energy is redistributed among different frequencies through nonlinear three-wave (triad) interactions satisfying
$f=f_1+f_2$
. With the present sign convention,
$T_f\gt 0$
represents a net nonlinear transfer of energy into the mode at frequency
$f$
from other frequencies, while
$T_f\lt 0$
represents a net transfer of energy out of the mode at
$f$
toward other frequencies.
The relation between these two quantities is given by (3.34), which shows that growth and transfer play opposite roles. If
$\gamma _f\gt 0$
then
$T_f\lt 0$
, meaning that an unstable mode generates energy and redistributes it to other frequencies through nonlinear coupling. If
$\gamma _f\lt 0$
then
$T_f\gt 0$
, meaning that a damped mode can persist only by receiving energy from other, typically dominant, modes. Therefore, the growth rate identifies where energy is created or removed in the spectrum, while the energy transfer function determines how that energy is transported across frequencies. Their combined behaviour controls the spectral evolution and saturation of the fluctuations.
Growth rate
$\gamma _f$
and nonlinear energy transfer
$T_f/2P_f$
as a function of frequency.
The growth rate and energy transfer function are plotted in figure 13. It can be observed that
$\gamma _f$
and
$T_f/2P_f$
do not sum to zero, indicating that the spatial stationarity condition given in (3.34) is not satisfied within this method. From figure 2 it is evident that the power spectra of the input and output data do not perfectly coincide, which leads to a mismatch in the complementary behaviour of
$\gamma _f$
and
$T_f/2P_f$
. This deviation is not unphysical; rather, it arises from the choice of the analytical linear and quadratic transfer functions adopted from Ritz & Powers (Reference Ritz and Powers1986). As discussed in § 3.1.4, these transfer functions correspond to a developing turbulence regime, in which the growth rate and nonlinear transfer function are not balanced and the input and output power spectra are unequal. A different choice of transfer functions, corresponding to fully developed turbulence (Ritz et al. Reference Ritz, Powers and Bengtson1989) and more consistent with the simulated dynamics, may improve the spatial stationarity of the data. This would better satisfy (3.34) and yield a closer complementary behaviour between the nonlinear transfer function and the linear growth rate.
Here we introduce a parameter called the energy imbalance parameter
$(\varDelta _{\mathrm{im}})$
, defined as the frequency-averaged spectral energy balance between growth rate and the nonlinear energy transfer function:
\begin{align} \varDelta _{\mathrm{im}} = \left \langle 2\gamma _f P_f + T_f \right \rangle _f . \end{align}
For an ideally spatially stationary system, the local balance is satisfied at each frequency, yielding
$\varDelta _{\mathrm{im}} = 0$
. A finite value of
$\varDelta _{\mathrm{im}}$
therefore measures the net excess or deficit of energy when linear growth and nonlinear redistribution do not fully compensate, indicating that the spectrum is not perfectly stationary. Equivalently,
$\varDelta _{\mathrm{im}} = 0$
corresponds to a fully developed turbulent regime, whereas deviations from zero indicate a developing turbulent state.
In our case, for the simulated data analysed using the Kim method, we obtain
$\varDelta _{\mathrm{im}} = 0.007$
, which is very small but indicates a developing turbulent state. This allows (3.34) to be written approximately as
$2\gamma _f P_f + T_f \approx 0$
, or
$\gamma _f \approx -T_f/(2P_f)$
, consistent with the behaviour observed in figure 13. Consequently, (3.34) implies that
$\langle Y_f Y_f^{*} \rangle \approx \langle X_f X_f^{*} \rangle$
, which supports the closure equation, i.e. (3.35). Therefore, although the turbulence is still developing, the spatial stationarity condition is approximately satisfied and the Kim method can be applied with reasonable accuracy.
Typical plasma and operational conditions of the IMPED device.
Table 1. Long description
A table with two columns and sixteen rows. The first column is labeled Parameter and the second column is labeled Range. The table lists various parameters and their corresponding ranges. Row 1: Magnetic field, 100-1200. Row 2: Operating pressure (Argon), 2 times 10 to the power of negative 5-1 times 10 to the power of negative 3. Row 3: Peak electron density, 8 times 10 to the power of 11. Row 4: Electron temperature, 2-5. Row 5: Ion temperature, 0.025. Row 6: Density scale length, 2-5. Row 7: Electron plasma frequency, 0.3-9. Row 8: Electron gyro frequency, 0.3-3.4. Row 9: Ion gyro frequency, 3.8-45.6. Row 10: Electron gyro radius, 0.04-0.3. Row 11: Ion gyro radius, 1-13. Row 12: Electron-neutral collision frequency, 6.6-213. Row 13: Ion-neutral collision frequency, 1.18-23.7. Row 14: Electron-ion collision frequency, 0.6-800. Row 15: Ion-neutral mean free path, 2.9-58. Row 16: Plasma radius, 8.
4. Application on experimental data and analysis
In the present study, experiments are carried out in a cylindrical plasma device known as the IMPED (Bose et al. Reference Bose, Chattopadhyay, Ghosh, Sengupta, Saxena and Pal2015a
,
Reference Bose, Kaur, Chattopadhyay, Ghosh and Saxenab
). The plasma is axially uniform over a length of approximately 2 m in the main chamber. The radius of the plasma column in the IMPED device is
$8\,\mathrm{cm}$
, and the corresponding operating conditions and plasma parameters are summarised in table 1. The IMPED has an additional control parameter,
$R_m$
, defined as the ratio of the magnetic field in the main chamber to that in the source chamber. Previously, the IMPED has been operated in different experimental regimes to investigate nonlinear interactions among primary instabilities such as DW, RT and Kelvin–Helmholtz (KH) modes. At a magnetic field of 550 G and a neutral pressure of
$1 \times 10^{-4}\,\text{mbar}$
, it has been demonstrated that variations in the density scale length and velocity shear scale length can lead to a transition from DW instability to KH instability, along with their nonlinear interactions (Roy et al. Reference Roy, Karmakar, Lachhvani, Chattopadhyay, Sen and Bose2025a
). With a further reduction in the magnetic field, reduced magnetic tension facilitates the onset of RT instability along with different KH instability modes (Roy et al. Reference Roy, Karmakar, Lachhvani, Khodiyar, Chattopadhyay and Sen2025b
). An increase in ion-neutral collisions introduces nonlinearities in both density and velocity fields, rendering the system more susceptible to nonlinear mode coupling. Recently, a study conducted by our group at 350 G and a pressure of
$4 \times 10^{-4}\,\text{mbar}$
provided a detailed identification of the primary modes, namely DW instability and RT instability, along with their nonlinear interactions and energy transfer mechanisms, which is submitted for publication. In the present work, for the same parameter, we perform a comprehensive analysis of nonlinear energy transfer between DW and RT modes using both the Ritz and Kim methods at two radial locations exhibiting distinct statistical properties. These locations are selected to assess the validity and applicability of each method with respect to the underlying statistical character and the spatial stationarity of the measured fluctuations, as discussed in the following sections. Long time-series measurements of plasma density fluctuations (
$\tilde {n}$
) and floating potential fluctuations (
$\tilde {\phi }_f$
), each of 10 s duration, are obtained using a three-tip Langmuir probe to have a large number of statistical realisations such that the error in the higher-order moments can be minimised as discussed in § 3.1.4. The fluctuations
$\tilde {n}$
and
$\tilde {\phi }_f$
are used as proxies for ion saturation current fluctuations (
$\tilde {I}_{sat}$
) and plasma potential fluctuations (
$\tilde {\phi }_p$
), respectively, assuming that electron temperature fluctuations (
$\tilde {T}_e$
) are negligible (Burin et al. Reference Burin, Tynan, Antar, Crocker and Holland2005; Roy et al. Reference Roy, Karmakar, Lachhvani, Chattopadhyay, Sen and Bose2025a
,
Reference Roy, Karmakar, Lachhvani, Khodiyar, Chattopadhyay and Senb
). Each probe tip (
$C_1$
,
$C_2$
and
$C_3$
) has a length of 4 mm and a diameter of 0.5 mm. The
$C_1$
and
$C_2$
probe tips are used to measure density fluctuations (
$\tilde {n}$
), whereas the
$C_3$
tip measures floating potential fluctuations (
$\tilde {\phi }_f$
). A detailed description of the probe configuration can be found in figure 2 of Roy et al. (Reference Roy, Karmakar, Lachhvani, Khodiyar, Chattopadhyay and Sen2025b
). The data are acquired at a sampling frequency of 1 MHz, corresponding to a sampling rate of 1 MS/s. For a data acquisition duration of 10 s, the total number of acquired samples is therefore
$10^7$
. The cross-phase between density and potential fluctuations measured using probes
$C_1$
and
$C_3$
, along with the normalised fluctuation level of the corresponding modes, is used to identify the underlying instability mechanisms. For the estimation of the nonlinear energy transfer function, the ‘two-probe method’ is employed, in which two Langmuir probes are spatially separated by a distance
$\Delta x$
to measure fluctuations at
$x$
,
$\phi (x,t)$
, and at
$x+\Delta x$
,
$\phi (x+\Delta x,t)$
. In the present experiment, density fluctuation data (
$\tilde {n}$
) are obtained using two poloidally separated Langmuir probes with a fixed separation of
$\Delta x = 6.4$
mm. In this framework, the fluctuation signal measured at position
$x$
,
$\phi (x,t)$
, is treated as the input signal (
$X$
), while the signal measured at
$x+\Delta x$
,
$\phi (x+\Delta x,t)$
, is considered as the output signal (
$Y$
), as discussed earlier. In the experiment, only two spatially separated signals are available. To identify which signal should be treated as the input and which as the output, a cross-correlation function (CCF) analysis (Bendat & Piersol Reference Bendat and Piersol2010) is performed. The CCF provides information on the direction of flow or wave propagation, thereby allowing the identification of propagation from
$X$
(input signal) to
$Y$
(output signal) or
$Y$
to
$X$
. The details of this analysis are discussed in the next section.
The CCF at a radial position of 5.76 cm showing a positive time delay between the input (
$X$
) and output (
$Y$
) signals.
4.1. Estimation of direction of flow propagation
To determine the direction of propagation, the CCF (Bendat & Piersol Reference Bendat and Piersol2010) is evaluated and is defined as
\begin{align} C_{XY}(\tau ) = \sum _{\tau =-\infty }^{\infty } X(t)\,Y(t+\tau ) {\rm d}\tau, \end{align}
where
$X(t)$
and
$Y(t)$
represent the time series of the two measured potential signals and
$\tau$
denotes the time lag. If a finite time delay exists between the two signals, the CCF peak is shifted away from zero lag. A positive lag indicates propagation from
$X$
to
$Y$
, validating the initial assumption, whereas a shift toward negative lag implies that the propagation direction is reversed. Figure 14 shows the CCF at a radial position of 5.76 cm. The cross-correlation between
$X$
and
$Y$
exhibits a distinct peak at a positive time lag, indicating a finite propagation delay. This confirms that the wave propagates from
$X$
to
$Y$
, justifying the identification of
$X$
as the input signal and
$Y$
as the output signal. The same analysis is performed in other radial locations before undergoing the nonlinear energy transfer analysis.
4.2. Experimental data selection
To ensure the validity and applicability of the Ritz and Kim methods on experimental data, density fluctuations (
$\tilde {n}$
) are selected at two distinct radial locations based on their statistical parameters (mainly skewness and kurtosis) and spatial stationarity. These criteria and their implications are discussed in detail in the following paragraphs.
4.2.1. Statistical parameters
From the discussion in the preceding sections, it is evident that the Ritz method cannot be applied to arbitrary time-series data. The applicability of the Ritz method is restricted to signals that are approximately Gaussian in nature, which implies that the data should exhibit low skewness and kurtosis. In contrast, the Kim method does not impose such constraints on the underlying data distribution. The statistical properties of the fluctuations, including skewness and kurtosis, are evaluated over a range of radial positions from 0 to 7.04 cm to guide the appropriate data selection for the application of the Ritz and Kim methods, as shown in figure 15(a).
Statistical properties of the fluctuation signal at different radial positions. Panel (a) shows radial variation of statistical parameters (mean, standard deviation, skewness and excess kurtosis) of the measured density fluctuations (
$\tilde {n}$
), while panel (b) compares the variance (
$\sigma ^2$
) and excess kurtosis to assess the validity of Millionshchikov’s approximation.
From figure 15(a) it can be observed that the data for the Ritz method can be selected below a radial position of 4 cm to maintain the excess kurtosis within
$0.40$
, consistent with the simulation results. Furthermore, figure 15(b) shows that, for radial positions up to approximately
$r = 4$
cm, the values of
$\sigma ^2$
and excess kurtosis remain relatively close. This suggests that Millionshchikov’s approximation may be applicable in this region, and consequently, the Ritz method is expected to provide reliable estimates. Moreover, at certain radial locations within
$r = 4$
cm, the agreement between
$\sigma ^2$
and excess kurtosis is particularly good, indicating conditions under which the Ritz method is most likely to be valid. In the present analysis, the absolute values of the statistical parameters are used, as their signs do not influence the interpretation of the results. Based on the above analysis, the energy transfer function analysis is performed at two distinct radial locations: one closer to the plasma centre and the other farther away. In this study, these locations correspond to
$r = 2.24$
cm and
$r = 5.76$
cm, respectively. The probability distribution functions of the input data at these radial locations are shown in figure 16. The smooth curves represent the Gaussian kernel density estimation, which provides a continuous estimate of the probability density function (PDF) by smoothing the discrete data using Gaussian kernels. From figure 16(b) it is observed that at
$r = 2.24$
cm the distribution is nearly Gaussian, with a skewness of
$0.042$
and an excess kurtosis of
$0.093$
. In contrast, at
$r = 5.76$
cm the distribution deviates significantly from Gaussian behaviour, exhibiting a skewness of
$0.205$
and an excess kurtosis of
$0.961$
, as shown in figure 16(a). Accordingly, Ritz’s method is expected to yield reliable results at
$r = 2.24$
cm but fail at
$r = 5.76$
cm. This is consistent with figure 15(b), where
$\sigma ^2$
remains close to the excess kurtosis at
$r = 2.24$
cm, but a significant deviation between the two quantities is observed at
$r = 5.76$
cm. Because Kim’s method does not depend explicitly on the skewness or kurtosis of the data, it is expected to provide consistent results at both radial locations. A detailed comparison of both approaches is presented in the following sections.
Probability distribution functions of the measured density fluctuations (
$\tilde {n}$
) at two radial locations. Gaussian kernel density estimation (KDE) of the PDF of the measured fluctuations.
The PDF at
$r=5.76$
cm exhibits a multimodal structure with more than one prominent peak, indicating a significant departure from a single Gaussian distribution. This behaviour is consistent with the enhanced nonlinear activity at this radial location, as evidenced by the higher-order statistical moments and the comparatively higher bicoherence values presented in the following sections. The multimodal nature of the PDF may arise from nonlinear intermittency and active energy transfer between coupled modes, leading to multiple preferred fluctuation-amplitude states.
4.2.2. Spatial stationarity
Spatial stationarity is defined as the condition in which the power associated with a given frequency component of a signal remains invariant in space. This implies that the power spectrum preserves both its shape and magnitude across different spatial locations. In this study, the fluctuation signals are considered spatially stationary when the power spectra of the input (
$X$
) and output (
$Y$
), separated by
$\Delta x$
, are equal and exhibit the same spectral shape as indicated by (3.34).
Previous studies (Ritz et al. Reference Ritz, Powers and Bengtson1989; Shen et al. Reference Shen, Li, Li and Li2020) have shown that Ritz’s method yields accurate results when applied to developing plasma turbulence, but fails for fully developed turbulence. Although Kim’s method was originally formulated for fully developed turbulence, it has also been shown to provide reliable results for developing turbulence models (Shen et al. Reference Shen, Li, Li and Li2020). In fully developed turbulence the input and output power spectra typically coincide, indicating good spatial stationarity and broadband characteristics, whereas in developing turbulence the spectra do not fully coincide, suggesting poorer spatial stationarity and that the turbulent cascade has not yet reached a fully developed state. However, Kim’s method becomes inapplicable when the input and output spectra do not follow similar trends (Kim et al. Reference Kim, Powers, Ritz, Miksad and Fischer1987). Consequently, Kim’s method remains suitable for developing turbulence cases in which the input and output spectra exhibit similar trends, even if they are not completely coinciding.
In the experimental data, coherent structures appear at low frequencies superimposed on a broadband turbulent background, indicating a developing turbulence regime. Although the input and output spectra do not fully coincide, they exhibit similar spectral trends, under which Kim’s method performs well. The Ritz method is also applicable under these developing turbulence conditions. As described earlier, the data are analysed at two radial locations,
$r = 2.24$
cm and
$r = 5.76$
cm. The power spectra of the input (
$X$
) and output (
$Y$
) signals at these two positions are shown in figure 17. Here
$P_X$
and
$P_Y$
denote the power spectra of the density fluctuations corresponding to the input signal
$X(t)$
and output signal
$Y(t)$
, respectively. The input and output signals are obtained from the
$C_2$
and
$C_1$
probe tips, which measure density fluctuations at different poloidal positions. At
$r = 5.76$
cm, the power spectra of
$X$
and
$Y$
nearly coincide, indicating good spatial stationarity. In contrast, at
$r = 2.24$
cm the spectra do not coincide exactly, although in both cases the input and output spectra exhibit similar trends. Since the power spectra of the input and output signals at both radial locations exhibit similar spectral trends, we argue that the Kim method is applicable at both positions. These arguments are verified through nonlinear energy transfer function analysis, which is discussed in the following sections.
Comparison of the input power spectrum (
$P_X$
) and output power spectrum (
$P_Y$
) of density fluctuations corresponding to the input (
$X(t)$
) and output (
$Y(t)$
) signals, measured using the
$C_2$
and
$C_1$
probe tips, respectively, at
$r = 5.76$
cm and
$r = 2.24$
cm. The probe tips are located at different poloidal positions.
4.3. Analysis and results
We now present the nonlinear transfer function analysis of the experimental data using both the Kim and Ritz methods. As discussed earlier, the selected radial locations in this study correspond to
$r = 2.24\,\mathrm{cm}$
and
$r = 5.76\,\mathrm{cm}$
, respectively. As discussed in the previous section (see figure 15), the data at
$r = 2.24$
cm exhibit low kurtosis, making them suitable for application of Ritz’s method. In contrast, Kim’s method does not depend explicitly on the skewness or kurtosis of the data and is therefore expected to provide consistent results at both radial locations. Moreover, from the stationarity viewpoint, the input and output spectra at both positions follow similar trends, supporting the applicability of Kim’s method at each location. To assess the validity of Ritz’s method, its results are systematically compared with those obtained using Kim’s method. Agreement between the two approaches indicates the applicability of Ritz’s method, whereas discrepancies signify its breakdown at a given radial position. Detailed analyses using both methods are presented in the following sections.
(a) Auto-bicoherence analysis. (b) Density fluctuation spectra for
$\tilde {n}$
of DW and RT modes at
$r = 5.76$
cm.
In the present analysis, the frequency range is restricted to
$40\,\mathrm{kHz}$
because the dominant spectral modes lie within this band. Although modes inside this range can exchange energy with modes outside the measured spectrum, such interactions cannot be directly resolved. As a result, energy transferred from the measured frequencies to unobserved higher-frequency modes may appear as linear damping, even though it actually originates from nonlinear spectral redistribution. Accordingly, the computed
$W_{\mathrm{mis}}$
does not represent the energy mismatch over the entire spectrum, but only the mismatch within the limited
$40\,\mathrm{kHz}$
band considered in the analysis.
For the estimation of the power spectra, bicoherence and energy transfer functions, the number of data points per segment is chosen as
$N = 2048$
and the number of ensemble realisations as
$M = 4882$
, corresponding to a total data length of
$10\,\mathrm{s}$
sampled at
$1\,\mathrm{MHz}$
. This choice yields a frequency resolution of approximately
$244\,\mathrm{Hz}$
. A large ensemble size is employed to ensure accurate estimation of higher-order moments, at the expense of spectral resolution. However, the resulting resolution of
$244\,\mathrm{Hz}$
is adequate for the present analysis, since the dominant modes are well resolved within this bandwidth.
4.3.1. Nonlinear energy transfer analysis at radial position
$r = 5.76$
cm
In this section the nonlinear energy transfer analysis at the radial location
$r = 5.76$
cm is presented using both the Kim and Ritz methods, as described in the following paragraphs. Before discussing the energy transfer results, we first present the power spectrum and auto- bicoherence of the density fluctuations (
$\tilde {n}$
), shown in figure 18.
The auto-power spectrum and auto–bicoherence are computed using (2.1) and (2.5), respectively, as described in § 2. The auto–bicoherence is evaluated over both positive and negative frequency regions by exploiting the symmetry relations discussed by Kim & Powers (Reference Kim and Powers1979). From figure 18(b), three dominant frequency modes are observed, which are listed in table 2, The instability types and poloidal wavenumbers corresponding to these modes are also listed in table 2. From figure 18(a) the frequency components at
$ f_1 = 8.3\,\mathrm{kHz}$
,
$ f_2 = 2.5\,\mathrm{kHz}$
and
$ f = 10.8\,\mathrm{kHz}$
are observed to be nonlinearly coupled, as indicated in the positive-bicoherence region. This interaction satisfies the frequency resonance condition
$ f_1 + f_2 = f$
as well as the corresponding wavenumber matching condition
$k_{\theta ,1} + k_{\theta ,2} = k_{\theta }$
, (see table 2). Since the wave propagation is predominantly in the poloidal direction, the resonance condition is evaluated using only the poloidal wavenumber component. However, a small deviation from the exact wavenumber matching condition can arise due to the finite spectral width of turbulent modes, plasma inhomogeneity, experimental uncertainty in wavenumber estimation and resonance broadening associated with collisional and nonlinear effects. Therefore, nonlinear coupling is considered satisfied within the experimental uncertainty range. In the negative region, additional interactions consistent with the relations
$ f - f_1 = f_2$
and
$ f - f_2 = f_1$
are also observed. The observed finite bicoherence confirms coherent phase coupling among these modes, with all identified interactions showing bicoherence values above the statistical significance level (2.7),
$b_{\mathrm{sig}}^{2} = {1}/{M} = {1}/{4882} \approx 2 \times 10^{-4}$
. Furthermore, the frequency mode at
$21\,\mathrm{kHz}$
(where the 21 kHz mode is identified as the second harmonic of the 10.8 kHz mode, within the spectral resolution of 244 Hz) is observed to interact with a broad band extending from the low-frequency region up to
$21\,\mathrm{kHz}$
in the negative-bicoherence domain. However, the bicoherence values indicate that the nonlinear interaction is strongest for the coupling between the 21 kHz mode and the 8–12.5 kHz frequency band in the negative-bicoherence region. In the positive-bicoherence region, frequency modes spanning from the low-frequency range up to
$10.8\,\mathrm{kHz}$
interact with the band between
$10.8\,\mathrm{kHz}$
and
$21\,\mathrm{kHz}$
, as illustrated in figure 18(a).
Measured mode frequencies with corresponding poloidal (
$k_{\theta }$
) and radial (
$k_r$
) wavenumbers (absolute values), and identified instability type.
Table 2. Long description
A table with three rows and three columns. The columns are labeled Frequency (kHz), kθ (cm⁻¹), and Mode type. The rows are as follows: Row 1: Frequency (kHz), 2.5; kθ (cm⁻¹), 0.36; Mode type, Drift-wave. Row 2: Frequency (kHz), 8.3; kθ (cm⁻¹), 1.12; Mode type, Mix mode of DW and RT. Row 3: Frequency (kHz), 10.8; kθ (cm⁻¹), 1.49; Mode type, Rayleigh-Taylor.
While these results clearly indicate the presence of three-wave coupling among the modes, bicoherence alone does not provide information about the direction or magnitude of the associated energy exchange. To address this limitation, we employ nonlinear energy transfer analysis to quantify the direction and rate of energy flow among the interacting modes, as discussed in the following paragraphs.
4.3.1.1. Analysis using Kim method
Figure 19(a,b) shows the power transfer function
$(T_f(f_1,f_2))$
and the corresponding nonlinear energy transfer function
$(W^{f}_{\mathrm{NL}})$
computed from two density fluctuation signals separated by 6.4 mm in the poloidal direction at the radial position
$r=5.76$
cm. In the two-dimensional contour map of
$ T_f(f_1,f_2)$
shown in figure 19(a), positive (red) contour regions correspond to energy gain at the coupled frequency
$ f = f_1 + f_2$
due to contributions from
$ f_1$
and
$ f_2$
, whereas negative (blue) contour regions indicate energy transfer from
$ f$
back to its coupling frequencies. Figure 19(b) presents the nonlinear energy transfer function
$ W^{f}_{\mathrm{NL}}$
, obtained by integrating the power transfer function along the resonance line
$ f = f_1 + f_2$
, where positive peaks indicate the energy gain and negative peaks indicate energy loss from the individual modes.
Results obtained using the Kim method at
$r = 5.76\,{\rm cm}$
. (a) Power transfer function
$T_f(f_1,f_2)$
. (b) Energy transfer function
$W^{f}_{\mathrm{NL}}(f)$
for
$\tilde {n}$
fluctuations associated with DW and RT modes. The coordinates
$(f_1,f_2)$
(in kHz) corresponding to the dominant nonlinear interaction regions are marked in panel (a), while the corresponding frequencies
$f=f_1+f_2$
are indicated in panel (b). Positive (red) regions in
$T_f(f_1,f_2)$
denote triads for which energy is transferred from the interacting modes
$(f_1,f_2)$
to the mode at
$f=f_1+f_2$
, producing positive peaks in
$W^{f}_{\mathrm{NL}}(f)$
. Negative (blue) regions denote triads for which energy is transferred from the mode at
$f=f_1+f_2$
to the interacting modes
$(f_1,f_2)$
, resulting in negative peaks (dips) in
$W^{f}_{\mathrm{NL}}(f)$
. For better visualisation, the colour scale in panel (a) is limited to the
$99.9$
th percentile of
$T_f(f_1,f_2)$
; values exceeding this threshold are saturated in the colour map. The contribution at
$f = 31.8$
kHz is relatively weak but remains above the noise floor; therefore, a zoomed inset is included in panel (b).
Figure 19. Long description
Panel A: A heatmap graph shows the power transfer function. The x-axis is labeled f1 (kHz) and the y-axis is labeled f2 (kHz). The color scale ranges from blue to red, indicating different values of the power transfer function. Specific coordinates corresponding to dominant nonlinear interaction regions are marked with arrows and labels. Panel B: A line graph displays the energy transfer function for fluctuations. The x-axis is labeled f = f1 + f2 (kHz) and the y-axis is labeled W_NL(f). The graph includes several peaks and dips, with annotations indicating specific frequencies. A zoomed inset is included to better visualize a weak contribution at a specific frequency.
From the two-dimensional contour plot of
$ T_f(f_1,f_2)$
as shown in figure 19(a), it is observed that the DW mode at 2.5 kHz is gaining energy from the mix mode and RT mode at 8.3 kHz and 10.8 kHz, respectively, following the interaction
$(f_1 = 10.8\,\text{kHz}) + (f_2 = -8.3\,\text{kHz})\rightarrow (f = 2.5\,\text{kHz})$
. The one-dimensional spectrum
$ W_{\mathrm{NL}}^{f}$
in figure 19(b) further supports this observation showing a positive peak near 2.5 kHz and negative peaks near 8.3 kHz and 10.8 kHz, indicating a net transfer of energy to the 2.5 kHz (DW) mode and from the 8.3 kHz (mix mode) and 10.8 kHz (RT) modes through nonlinear coupling. Furthermore, figure 19(a) shows that the RT mode at
$10.8\,\mathrm{kHz}$
loses energy by transferring it to modes at
$21\,\mathrm{kHz}$
,
$10.2\,\mathrm{kHz}$
and
$31.8\,\mathrm{kHz}$
. These transfers occur through the triad interactions
$(f = 10.8\,\mathrm{kHz}) \rightarrow (f_1 = 21\,\mathrm{kHz}) + (f_2 = -10.2\,\mathrm{kHz})$
and
$(f = 10.8\,\mathrm{kHz}) \rightarrow (f_1 = 31.8\,\mathrm{kHz}) + (f_2 = -21\,\mathrm{kHz})$
, which are identified by the negative contour (blue) regions in the two-dimensional power transfer function. Furthermore, energy is transferred from the 10.8 kHz and 10.2 kHz modes to the 21 kHz mode via the triad interaction
$ (f_1 = 10.8\,\mathrm{kHz}) + (f_2 = 10.2\,\mathrm{kHz}) \rightarrow (f = 21\,\mathrm{kHz}),$
as indicated by the positive (red) contour region in figure 19(a). The
$W_{\mathrm{NL}}^{f}$
spectrum in figure 19(b) further supports these observations, exhibiting a negative peak at
$10.8\,\mathrm{kHz}$
and positive peaks at
$21\,\mathrm{kHz}$
. The positive peaks at
$31.8\,\mathrm{kHz}$
and
$10.2\,\mathrm{kHz}$
are comparatively weak, indicating that a net transfer of energy occurs away from the
$10.8\,\mathrm{kHz}$
mode toward the
$21\,\mathrm{kHz}$
mode. In addition, the
$21\,\mathrm{kHz}$
mode interacts with a band of frequencies in the range
$10.4$
–
$12.5\,\mathrm{kHz}$
and transfers energy to the
$8.5$
–
$10.6\,\mathrm{kHz}$
band, as indicated by the positive contour regions (red and yellow) in figure 19(a). This observation is consistent with figure 19(b), which exhibits a broad positive peak spanning approximately
$8.5$
–
$10.6\,\mathrm{kHz}$
. The remaining peaks in the
$W_{\mathrm{NL}}^{f}$
spectrum correspond to relatively weak nonlinear interactions.
Results obtained using the Ritz method at
$r=5.76\,\mathrm{cm}$
. (a) Power transfer function
$T_f(f_1,f_2)$
. (b) Energy transfer function
$W^{f}_{\mathrm{NL}}(f)$
for
$\tilde {n}$
fluctuations associated with DW and RT modes. The coordinates
$(f_1,f_2)$
(in kHz) corresponding to the dominant nonlinear interaction regions are marked in panel (a), while the corresponding frequencies
$f=f_1+f_2$
are indicated in panel (b). Positive (red) regions in
$T_f(f_1,f_2)$
denote triads for which energy is transferred from the interacting modes
$(f_1,f_2)$
to the mode at
$f=f_1+f_2$
, producing positive peaks in
$W^{f}_{\mathrm{NL}}(f)$
. Negative (blue) regions denote triads for which energy is transferred from the mode at
$f=f_1+f_2$
to the interacting modes
$(f_1,f_2)$
, resulting in negative peaks (dips) in
$W^{f}_{\mathrm{NL}}(f)$
. For better visualisation, the colour scale in panel (a) is limited to the
$99.9$
th percentile of
$T_f(f_1,f_2)$
; values exceeding this threshold are saturated in the colour map.
4.3.1.2. Analysis using Ritz method
Here the nonlinear energy transfer analysis of the same data at the radial location
$ r = 5.76\,\mathrm{cm}$
is performed using the Ritz method. Figure 20(a,b) presents the power transfer function
$ T_f(f_1,f_2)$
and the corresponding nonlinear energy transfer rate
$ W^{f}_{\mathrm{NL}}$
, computed from the same density fluctuation data used for the Kim method analysis in the previous section. From the two-dimensional contour plot of
$T_f(f_1,f_2)$
shown in figure 20(a), the RT mode at
$10.8\,\mathrm{kHz}$
is observed to gain energy from the DW mode at
$2.5\,\mathrm{kHz}$
and the mix mode at
$8.3\,\mathrm{kHz}$
through the triad interaction
$(f_1 = 2.5\,\mathrm{kHz}) + (f_2 = 8.3\,\mathrm{kHz}) \rightarrow (f = 10.8\,\mathrm{kHz})$
indicated by the positive contour (red) region. In addition, the mode at
$21\,\mathrm{kHz}$
interacts with a frequency band from
$9.5\,\mathrm{kHz}$
to
$13.5\,\mathrm{kHz}$
, producing positive contour (red) regions in the power transfer function. The mode at
$12.7\,\mathrm{kHz}$
transfers energy to
$21\,\mathrm{kHz}$
and
$8.3\,\mathrm{kHz}$
via the triad interaction
$(f = 12.7\,\mathrm{kHz}) \rightarrow (f_1 = 21\,\mathrm{kHz}) + (f_2 = -8.3\,\mathrm{kHz})$
, which appears as negative contours (blue) in the transfer map. Furthermore, the RT mode at
$10.8\,\mathrm{kHz}$
loses energy to it’s harmonics at
$21\,\mathrm{kHz}$
and
$31.8\,\mathrm{kHz}$
modes through the triad interaction
$(f = 10.8\,\mathrm{kHz}) \rightarrow (f_1 = 31.8\,\mathrm{kHz}) + (f_2 = -21\,\mathrm{kHz})$
, also indicated by negative contour (blue) regions. The one-dimensional spectrum
$W_{\mathrm{NL}}^{f}$
in figure 20(b) supports these trends, showing negative peaks near
$2.5\,\mathrm{kHz}$
and
$12.7\,\mathrm{kHz}$
, and positive peaks in the
$7.5$
–
$11.5\,\mathrm{kHz}$
band and near
$31.8\,\mathrm{kHz}$
. The remaining peaks correspond to relatively weak nonlinear interactions.
The accuracy of the nonlinear energy transfer estimation is assessed using the energy mismatch parameter
$ W_{\mathrm{mis}}$
(shown in the legends of figures 19
b and 20
b). For the density fluctuation
$ \tilde {n}$
at the radial position
$ r = 5.76\,\mathrm{cm}$
, the mismatch is
$ W_{\mathrm{mis}} = 0.456 \approx 0.46$
for the Kim method, whereas for the Ritz method it is
$ W_{\mathrm{mis}} = -0.199 \approx -0.20$
. This implies that about 46 % of the net transfer is unbalanced for the Kim method and nearly 20 % for the Ritz method. In the Kim method,
$W_{\mathrm{mis}}$
is positive, indicating that the net nonlinear transfer into the resolved modes slightly exceeds the transfer out of them. This imbalance likely arises from additional energy input from the background plasma or from couplings with other fluctuations (e.g. potential fluctuations) that are not included in the present analysis. In contrast, for the Ritz method,
$W_{\mathrm{mis}}$
is negative, implying that the estimated loss marginally exceeds the corresponding gain. This suggests that part of the transferred energy is exchanged with the background plasma or with unaccounted modes outside the scope of the present analysis.
From the nonlinear energy transfer analysis using both methods at the radial location 5.76 cm, it is evident that the results differ significantly and, in some cases, even show opposite trends. From figure 20(a,b) it is observed that the contour plots obtained using the Ritz method are comparatively noisier than those obtained using the Kim method. Consistent with this observation, the corresponding one-dimensional profile of
$W_{NL}$
exhibits broader peaks around the 2.5 kHz and 21 kHz modes, along with several spurious oscillations that are not expected from the physical analysis. This discrepancy can be understood from the data selection discussion in § 4.2, where the kurtosis is found to be
$0.96$
(figure 15), indicating a strong deviation from Gaussian statistics. From the simulation study, it was also observed that the Ritz method becomes unreliable for kurtosis values beyond
$\sim 0.40$
. Therefore, at the radial location
$ r = 5.76\,\mathrm{cm}$
, the Ritz method is not valid and its results are not physically reliable. In contrast, the Kim method does not rely on such statistical approximations and remains applicable for non-Gaussian data. Moreover, at this radial location we observe good spatial stationarity (figure 17
a), and the input and output power spectra exhibit similar spectral trends; despite the developing turbulence regime indicated by the spectra, the Kim method provides a more reliable description of nonlinear energy transfer at this radial position. In summary, at this radial location we find that the RT mode at
$10.8\,\mathrm{kHz}$
transfers energy to the DW mode at
$2.5\,\mathrm{kHz}$
and mix mode at
$8.3\,\mathrm{kHz}$
and also to other modes at
$10.2\,\mathrm{kHz}$
and
$21\,\mathrm{kHz}$
. In the next section we present the nonlinear energy transfer analysis using both the Ritz and Kim methods at a radial location close to the plasma centre
$(r = 2.24\,\text{cm})$
, where the data exhibit low kurtosis, in order to examine whether the results obtained from the Ritz method agree with those from the Kim method.
Measured mode frequencies with corresponding poloidal (
$k_{\theta }$
) and radial (
$k_r$
) wavenumbers (absolute values), and identified instability type.
Table 3. Long description
The table presents data on measured mode frequencies with corresponding poloidal and radial wavenumbers and identified instability types. It has three rows and three columns. The columns are labeled Frequency (kHz), kθ (cm⁻¹), and Mode type. Row 1: Frequency (kHz), 2.5; kθ (cm⁻¹), 0.95; Mode type, Drift-wave. Row 2: Frequency (kHz), 8.8; kθ (cm⁻¹), 1.94; Mode type, Mix mode of DW and RT. Row 3: Frequency (kHz), 11.3; kθ (cm⁻¹), 3.15; Mode type, Rayleigh-Taylor.
(a) Auto-bicoherence analysis. (b) Density fluctuation spectra for
$\tilde {n}$
of DW and RT modes at
$r = 2.24$
cm.
4.3.2. Nonlinear energy transfer analysis at radial position
$r = 2.24$
cm
We analyse the nonlinear energy transfer at the radial position
$ r = 2.24\,\mathrm{cm}$
using both the Kim and Ritz approaches. As a first step, the spectral properties of the density fluctuations
$ \tilde {n}$
are examined through the power spectrum and auto-bicoherence, presented in figure 21, before discussing the corresponding energy transfer results.
The spectral resolution and data length are kept the same as in the case of the radial position 5.76 cm. From the auto-power spectrum plot in figure 21(b), three dominant frequency modes are observed (see table 3). The instability types, radial and poloidal wavenumbers corresponding to these modes are also listed in table 3. From figure 21(a) it is evident that the frequency components at
$ f_1 = 8.8\,\text{kHz}$
,
$ f_2 = 2.5\,\text{kHz}$
and
$ f = 11.3\,\text{kHz}$
nonlinearly interact as indicated in the positive-bicoherence region. This interaction satisfies the frequency resonance conditions
$ f_1 + f_2 = f$
. In the negative region, additional interactions consistent with the relations
$ f - f_1 = f_2$
and
$ f - f_2 = f_1$
are also observed. Similarly, we observe another relatively weak interaction:
$ f_1 = 6.3\,\text{kHz}$
,
$ f_2 = 2.5\,\text{kHz}$
and
$ f = 8.8\,\text{kHz}$
nonlinearly interact, satisfying the resonance condition
$ f_1 + f_2 = f$
, as evident in the positive-bicoherence region. Similarly, additional interactions are observed in the negative region as discussed earlier. The observed finite bicoherence confirms coherent phase coupling among these modes, with all identified interactions showing bicoherence values above the statistical significance level,
$b_{\mathrm{sig}}^{2} = {1}/{4882} \approx 2 \times 10^{-4}$
. However, the wavenumber resonance condition between the poloidal wavenumbers
$(k_{\theta })$
(see table 3) is not exactly satisfied, which may be attributed to relatively weak nonlinear three-wave coupling (indicated by low bicoherence magnitude) at this radial location and uncertainties in wavenumber estimation. Since the wave propagation is predominantly poloidal, the resonance condition is evaluated using only the poloidal wavenumber component. Small deviations from exact matching can arise due to the finite spectral width of turbulent fluctuations, plasma inhomogeneity and nonlinear or collisional resonance broadening. Therefore, the nonlinear coupling condition is considered satisfied within experimental uncertainty.
We now determine the direction of energy transfer among these frequency components using nonlinear energy transfer analysis based on both the Kim and Ritz methods, as discussed below.
4.3.2.1. Analysis using Kim method
Figure 22(a,b) presents the power transfer function
$ T_f(f_1,f_2)$
and the associated nonlinear energy transfer function
$ W^{f}_{\mathrm{NL}}$
, evaluated from two density fluctuation signals separated by 6.4 mm in the poloidal direction at the radial location
$ r = 2.24$
cm. The computation uses the same frequency resolution as adopted in the bicoherence analysis.
Results obtained using the Kim method at
$r = 2.24\,\text{cm}$
. (a) Power transfer function
$T_f(f_1,f_2)$
. (b) Energy transfer function
$W^{f}_{\mathrm{NL}}(f)$
for
$\tilde {n}$
fluctuations associated with DW and RT modes. The coordinates
$(f_1,f_2)$
(in kHz) corresponding to the dominant nonlinear interaction regions are marked in panel (a), while the corresponding frequencies
$f=f_1+f_2$
are indicated in panel (b). Positive (red) regions in
$T_f(f_1,f_2)$
denote triads for which energy is transferred from the interacting modes
$(f_1,f_2)$
to the mode at
$f=f_1+f_2$
, producing positive peaks in
$W^{f}_{\mathrm{NL}}(f)$
. Negative (blue) regions denote triads for which energy is transferred from the mode at
$f=f_1+f_2$
to the interacting modes
$(f_1,f_2)$
, resulting in negative peaks (dips) in
$W^{f}_{\mathrm{NL}}(f)$
. For better visualisation, the colour scale in panel (a) is limited to the
$99.9$
th percentile of
$T_f(f_1,f_2)$
; values exceeding this threshold are saturated in the colour map.
Figure 22. Long description
Panel A: A heatmap graph shows the power transfer function. The x-axis is labeled f1 (kHz) and the y-axis is labeled f2 (kHz). The color scale ranges from blue to red, indicating different values of the power transfer function. Specific coordinates corresponding to dominant nonlinear interaction regions are marked. Panel B: A line graph displays the energy transfer function. The x-axis is labeled f = f1 + f2 (kHz) and the y-axis represents the energy transfer function. Peaks and dips in the graph indicate positive and negative energy transfers, respectively. Specific frequencies corresponding to dominant nonlinear interactions are marked.
From the two-dimensional contour plot of
$ T_f(f_1,f_2)$
shown in figure 22(a), it is observed that the RT mode at 11.3 kHz loses energy through nonlinear coupling to the mix mode at 8.8 kHz and the DW mode at 2.5 kHz, satisfying the interaction
$(f = 11.3\,\text{kHz}) \rightarrow (f_1 = 8.8\,\text{kHz}) + (f_2 = 2.5\,\text{kHz})$
. Another relatively weak interaction is also evident,
$(f_1 = 6.3\,\text{kHz}) + (f_2 = 2.5\,\text{kHz}) \rightarrow (f = 8.8\,\text{kHz})$
, indicating that the mix mode at 8.8 kHz gains energy through nonlinear coupling with the 6.3 kHz mode and DW mode at 2.5 kHz. Another triad interaction
$(f_1 = 8.8\,\text{kHz}) + (f_2 = -6.3\,\text{kHz}) \rightarrow (f = 2.5\,\text{kHz})$
shows that the DW mode at 2.5 kHz gains energy from the modes at 6.3 kHz and 8.8 kHz. The one-dimensional spectrum
$ W_{\mathrm{NL}}^{f}$
in figure 22(b) supports these observations, showing a negative peak near 11.3 kHz and prominent positive peaks near 8.8 kHz and 2.5 kHz. The weak positive peak at 2.5 kHz arises because this mode gains energy in one interaction and loses energy in another, leading to a small net gain. In the case of the 8.8 kHz mode, energy is gained through two nonlinear interactions, while energy is lost through only a single interaction, as discussed earlier. Consequently, this mode exhibits a prominent net energy gain, which appears as a pronounced positive peak in the
$W_{\text{NL}}^f$
spectrum. The results indicate a net nonlinear energy transfer from the 11.3 kHz RT mode toward the 8.8 kHz mix mode and the 2.5 kHz DW mode, while the faint contribution near 6.3 kHz reflects comparatively weak nonlinear coupling, as evidenced by the low bicoherence associated with this mode. Other peaks in
$ W_{\mathrm{NL}}^{f}$
correspond to minor interactions. Now we analyse the same data with the Ritz method, which is presented in the next section.
4.3.2.2. Analysis using Ritz method
Here, the nonlinear energy transfer analysis of the same density fluctuation data at the radial location
$ r = 2.24\,\mathrm{cm}$
is performed using the Ritz method. The two-dimensional contour plot of
$ T_f(f_1,f_2)$
shown in figure 23(a) is qualitatively identical to that obtained using the Kim method in figure 22(a), indicating consistent nonlinear interactions. The dominant coupling processes are the same as those discussed for the Kim method. The one-dimensional spectrum
$ W_{\mathrm{NL}}^{f}$
in figure 23(b) further supports these observations and shows the same trend as in the Kim method, leading to the same conclusion: a net nonlinear energy transfer from the 11.3 kHz RT mode to the 8.8 kHz mix mode and the 2.5 kHz DW mode. The faint contribution near 6.3 kHz indicates relatively small coupling as evidenced by the low bicoherence associated with this mode, while the remaining peaks in
$ W_{\mathrm{NL}}^{f}$
correspond to minor interactions.
Results obtained using the Ritz method at
$r = 2.24\,\text{cm}$
. (a) Power transfer function
$T_f(f_1,f_2)$
. (b) Energy transfer function
$W^{f}_{\mathrm{NL}}(f)$
for
$\tilde {n}$
fluctuations associated with DW and RT modes. The coordinates
$(f_1,f_2)$
(in kHz) corresponding to the dominant nonlinear interaction regions are marked in panel (a), while the corresponding frequencies
$f=f_1+f_2$
are indicated in panel (b). Positive (red) regions in
$T_f(f_1,f_2)$
denote triads for which energy is transferred from the interacting modes
$(f_1,f_2)$
to the mode at
$f=f_1+f_2$
, producing positive peaks in
$W^{f}_{\mathrm{NL}}(f)$
. Negative (blue) regions denote triads for which energy is transferred from the mode at
$f=f_1+f_2$
to the interacting modes
$(f_1,f_2)$
, resulting in negative peaks (dips) in
$W^{f}_{\mathrm{NL}}(f)$
. For better visualisation, the colour scale in panel (a) is limited to the
$99.9$
th percentile of
$T_f(f_1,f_2)$
; values exceeding this threshold are saturated in the colour map.
The accuracy of the nonlinear energy transfer estimation is evaluated using the energy mismatch parameter
$ W_{\mathrm{mis}}$
, shown in the legends of figures 22(b) and 23(b). For the density fluctuation
$ \tilde {n}$
at
$ r = 2.24\,\mathrm{cm}$
, the mismatch is
$ W_{\mathrm{mis}} = -0.1837 \approx -0.18$
for the Kim method and
$ W_{\mathrm{mis}} = -0.056 \approx -0.06$
for the Ritz method. This indicates that about 18 % and 6 % of the net transfer remain unbalanced for the Kim and Ritz methods, respectively. In both cases, the negative sign of
$ W_{\mathrm{mis}}$
implies that the estimated loss slightly exceeds the corresponding gain, suggesting that part of the transferred energy is exchanged with background plasma or other unaccounted modes not included in the present analysis.
From the nonlinear energy transfer analysis at the radial location
$ r = 2.24\,\mathrm{cm}$
using both the Kim and Ritz methods, it is evident that the results obtained from the two approaches agree closely. In this region, the Ritz method performs reliably because the kurtosis of the data is low (
$ \approx 0.09$
), which is well below the critical limit (
$\approx 0.40$
) identified from the simulation study for the validity of the Ritz method. This good agreement with the Kim method indicates that the Ritz method can be safely applied in regions characterised by low-kurtosis statistics. However, from figure 17(b) we find that the spatial stationarity at this location is poor, although the input and output power spectra exhibit similar spectral shapes, which supports the validity of the Kim method at this location. We also find that at this location, the Ritz method yields a smaller value of the mismatch parameter
$W_{\text{mis}}$
than the Kim method, indicating a smaller residual energy imbalance in the estimated nonlinear energy budget. It should be noted that
$W_{\text{mis}}$
is not used as a direct criterion for method applicability, since its value can be influenced by factors such as finite frequency-band truncation, transfer to unmeasured fields, transfer to mean profiles and statistical uncertainties. In summary, at this radial location the dominant interaction shows that the RT mode at
$11.3\,\mathrm{kHz}$
transfers energy to the DW mode at
$2.5\,\mathrm{kHz}$
as well as to the mix mode at
$8.8\,\mathrm{kHz}$
. In the next section we present the growth rate analysis at the radial location
$r = 2.24\,\mathrm{cm}$
, where both the Kim and Ritz methods perform reliably.
4.3.3. Discussion on growth rate
In the previous section we identified which modes transfer energy to other modes. In this section we present a growth rate analysis that allows us to examine the growth or damping of individual modes and how these behaviours are connected to the energy transfer processes. Specifically, we analyse the growth rates at the radial position
$ r = 2.24\,\mathrm{cm}$
(since at this location both methods are applicable), as shown in figure 24.
Growth rate spectrum
$\gamma _f$
as a function of frequency. The dashed line indicates
$\gamma _f = 0$
.
Here we provide a qualitative discussion of the growth rate plot rather than a quantitative analysis. This approach allows us to focus on the overall trends, such as which modes are growing or damping and how their behaviour correlates with the energy transfer between modes, without delving into exact numerical values or rates. From figure 24, the mode near
$ f = 10.8\,\mathrm{kHz}$
has a small positive growth rate and is therefore unstable, whereas the modes near
$ f = 2.5\,\mathrm{kHz}$
and
$ f = 8.3\,\mathrm{kHz}$
are strongly damped. This implies that the
$10.8\,\mathrm{kHz}$
mode grows linearly and the energy gained by this mode is redistributed to other modes or to the background equilibrium due to nonlinear interaction, while the lower-frequency modes are damped and receive energy through nonlinear interactions following the spatial stationarity condition (3.34). The energy transfer analysis reveals the nonlinear redistribution of energy, demonstrating that the
$10.8\,\mathrm{kHz}$
mode acts as a donor, transferring energy to the
$8.3\,\mathrm{kHz}$
and
$2.5\,\mathrm{kHz}$
modes, which is clearly reflected in the energy transfer function shown in figure 22. Thus, although the growth rate spectrum identifies which modes are linearly driven or damped, the direction of spectral energy flow is governed by nonlinear interactions captured by the energy transfer analysis. Furthermore, the positive growth rates in the higher-frequency band indicate that these modes are unstable, and the energy contained in these frequencies is transported to lower-frequency bands or to the background equilibrium, consistent with the damping observed in the low-frequency modes. In addition, we evaluate the energy imbalance parameter
$\varDelta _{\mathrm{im}}$
from (3.36) at this location and obtain a value of
$0.06$
. Although small, this finite value indicates a developing turbulent state. Nevertheless, the small magnitude of
$\varDelta _{\mathrm{im}}$
suggests that the Kim method remains valid and applicable at this position.
5. Discussion
The behaviour of nonlinear interactions among various instability modes, namely RT and drift modes, observed in the measured plasma density fluctuations in the IMPED, are investigated in detail. Prior to applying the methods to experimental data, the computational implementations of both the Ritz and Kim methods are benchmarked against analytical transfer functions obtained from a developing plasma turbulence model. The existing computational methods are further explored with suitable improvisations to overcome their inherent limitations and constraints. It is observed that the Ritz method performs well for fluctuations that are close to Gaussian, with low skewness and kurtosis. However, for cases with large kurtosis, the method deteriorates because its formulation approximates fourth-order moments, thereby introducing errors in the computation of the energy transfer functions. In contrast, Kim’s method performs much better for fluctuation data with larger kurtosis and yields consistent estimates of the energy transfer functions, because the method does not rely on the Millionshchikov approximation and instead retains the fourth-order moments.
Both Ritz and Kim methods are tested to compute the energy transfer among the modes observed in the experimental fluctuation data of the IMPED. The applicability of these methods with respect to spatial stationarity is also examined. We find that Kim’s method yields reliable estimates of the energy transfer function even under weak spatial stationarity in developing turbulence, provided the spectral shapes of the power spectra are identical. These trends are clearly observed in the experimental data at the two radial locations,
$r = 5.76\,\mathrm{cm}$
and
$r = 2.24\,\mathrm{cm}$
. In this work, we observe that the Kim method performs well at both radial locations, whereas the Ritz method yields reliable results only at
$r = 2.24\,\mathrm{cm}$
. This behaviour can be explained by the statistical nature of the data and the degree of spatial stationarity. In particular, the mismatch parameter
$W_{\text{mis}}$
for the Ritz method is relatively small at
$r = 2.24\,\mathrm{cm}$
(about
$0.06$
), indicating a relatively small residual energy mismatch at this location. The applicability of the Ritz method is assessed primarily from the statistical properties of the fluctuations, particularly the validity of the Millionshchikov approximation and the low excess kurtosis observed at this radial position.
6. Conclusion
The nonlinear energy transfer among DW and RT modes in the IMPED plasma has been investigated using the Ritz and Kim methods. Benchmarking against a developing plasma turbulence model showed that the Ritz method provides reliable estimates for near-Gaussian fluctuations, whereas its performance deteriorates with increasing kurtosis. In contrast, the Kim method remains robust for both Gaussian and non-Gaussian fluctuations. Consistent trends are observed in the experimental data. The Kim method yields reliable estimates at both
$r = 5.76\,\mathrm{cm}$
and
$r = 2.24\,\mathrm{cm}$
, whereas the Ritz method performs reliably only at
$r = 2.24\,\mathrm{cm}$
. These results highlight the importance of fluctuation statistics, particularly excess kurtosis, in assessing the suitability of energy transfer estimation methods. The energy transfer analysis further reveals that, at
$r = 5.76\,\mathrm{cm}$
, the RT mode at
$10.8\,\mathrm{kHz}$
and the mixed mode at
$8.3\,\mathrm{kHz}$
primarily transfer energy to the DW mode at
$2.5\,\mathrm{kHz}$
. In contrast, at
$r = 2.24\,\mathrm{cm}$
, the RT mode at
$11.3\,\mathrm{kHz}$
transfers energy to both the DW mode at
$2.5\,\mathrm{kHz}$
and the mixed mode at
$8.8\,\mathrm{kHz}$
. Thus, the analysis indicates a net transfer of energy from RT-dominated fluctuations toward a comparatively lower-frequency DW mode at both radial locations. However, since the nonlinear interaction is significantly stronger at
$r = 5.76\,\mathrm{cm}$
, the overall direction of energy transfer is inferred primarily from the results at this radial location. Therefore, the dominant nonlinear energy transfer process can be interpreted as the transfer of energy from an RT mode and a mixed DW–RT mode to a comparatively low-frequency DW mode.
7. Outlook
In the present study, the energy transfer analysis is restricted to a single-field framework, where the transfer functions are estimated using the Kim and Ritz methods. Although single-field analysis provides useful insight into nonlinear energy transfer within an individual fluctuating quantity, it cannot capture the full complexity of plasma turbulence, in which multiple coupled fields such as density, potential, temperature and velocity interact nonlinearly (Xu et al. Reference Xu, Tynan, Holland, Yan, Muller and Yu2009, Reference Xu, Tynan, Holland, Yan, Muller and Yu2010). In DW turbulence experiments, where the density-potential cross-correlation is well below unity, key processes such as DW zonal-flow coupling and the exchange between internal and kinetic energy channels cannot be fully described within a single-field formulation (Xu et al. Reference Xu, Tynan, Holland, Yan, Muller and Yu2009). A multi-field turbulence approach is therefore required to quantify nonlinear coupling and energy transfer among different fluctuating fields, as demonstrated in experimental studies and theoretical treatments of plasma turbulence (Manz et al. Reference Manz, Ramisch, Stroth, Naulin and Scott2008; Xu et al. Reference Xu, Tynan, Holland, Yan, Muller and Yu2009; Itoh et al. Reference Itoh, Itoh, Nagashima and Kosuga2017). Our future work will involve extending the current methodology to a multi-field framework, building upon the validated single-field results achieved here through the Kim and Ritz methods.
Acknowledgements
The authors acknowledge Bhoomi Khodiyar, Jignesh Patel, Kalpesh Doshi, Kirti Mahajan, Manisha Bhandankar, Praveenlal, Minsha Shah and Vijay Patel for their experimental support.
Editor Cary Forest thanks the referees for their advice in evaluating this article.
Declaration of interests
The authors report no conflict of interest.
Lf
Qff1,f2
μ=0.001
σ=1.00
=0.001
=0.01
μ=0.002
σ=0.52
=0.15
=0.40
Lf
Qff1,f2
Lf
|Qf(f1,f2)|
Lf
|Qf(f1,f2)|
γf
Tf/2Pf
X
Y
n~
σ2
n~
PX
PY
X(t)
Y(t)
C2
C1
r=5.76
r=2.24
n~
r=5.76
kθ
kr
r=5.76cm
Tf(f1,f2)
WNLf(f)
n~
(f1,f2)
f=f1+f2
Tf(f1,f2)
(f1,f2)
f=f1+f2
WNLf(f)
f=f1+f2
(f1,f2)
WNLf(f)
99.9
Tf(f1,f2)
f=31.8
r=5.76cm
Tf(f1,f2)
WNLf(f)
n~
(f1,f2)
f=f1+f2
Tf(f1,f2)
(f1,f2)
f=f1+f2
WNLf(f)
f=f1+f2
(f1,f2)
WNLf(f)
99.9
Tf(f1,f2)
kθ
kr
n~
r=2.24
r=2.24cm
Tf(f1,f2)
WNLf(f)
n~
(f1,f2)
f=f1+f2
Tf(f1,f2)
(f1,f2)
f=f1+f2
WNLf(f)
f=f1+f2
(f1,f2)
WNLf(f)
99.9
Tf(f1,f2)
r=2.24cm
Tf(f1,f2)
WNLf(f)
n~
(f1,f2)
f=f1+f2
Tf(f1,f2)
(f1,f2)
f=f1+f2
WNLf(f)
f=f1+f2
(f1,f2)
WNLf(f)
99.9
Tf(f1,f2)
γf
γf=0
Open access
