Data-driven modal analysis of nonlinear quantities in turbulent plasmas using multi-field singular value decomposition

Nonlinear dynamics in the two-dimensional multi-component plasma turbulence described by the Hasegawa–Wakatani equation is investigated by using a data-driven modal analysis with the singular value decomposition (SVD). The conventional SVD is extended to ‘multi-field SVD’ which can decompose multiple turbulence fields simultaneously by a single set of orthonormal basis functions without imposing a priori scale separations. Then, in addition to the mode amplitude labeled by the singular value, the information on the phase relations in the nonlinear quantities such as a transport flux or a triad energy transfer is extracted in the mode space. Through applications to the two-dimensional plasma turbulence, it is revealed that the multi-field SVD can extract the dominant spatial structures for the turbulent transport and the nonlinear energy transfer, preserving the multi-scale nature of the original turbulent fields. It is also demonstrated that one can reduce the dimensionality or information using the multi-field SVD through comparisons with the conventional Fourier decomposition.


Introduction
In magnetically confined fusion plasmas, the drift-wave turbulence is driven by the background temperature and/or density gradient. This causes the turbulent transport of particles and energy to the outside of the confined plasmas. On the other hand, large-scale sheared flows called zonal flows are spontaneously generated through the nonlinear interactions in the drift-wave turbulence [1][2][3]. The zonal flows, which are observed not only in confined plasmas but also in neutral fluids such as the atmosphere of Earth and Jupiter [4,5], have an important role in the suppression of the turbulent transport by reducing the characteristic length of the turbulent vortices. In order to understand the physics of turbulent transport and zonal-flow generation and to achieve high-performance plasmas, many studies have been carried out both experimentally and theoretically.
In these decades, the high-resolution measurements and large-scale numerical simulations in fluids [6][7][8][9] and plasmas [10][11][12][13] including the turbulence dynamics have made remarkable progress with the development of sensing technologies and high-performance computers with improved algorithms. Nowadays, the multi-dimensional and multicomponent flow or turbulence data obtained by these measurements or simulations become so huge and complex that the analyses and understandings of their characteristics are no longer straightforward. The data-driven approach is, therefore, required to lower the dimensionality and to extract the essential information contributing to the physics from the enormous data.
To this end, the mode decomposition of flow fields is proposed and utilized [14], such as the singular value decomposition (SVD) [15][16][17], the proper orthogonal decomposition (POD) [18,19], and the dynamical mode decomposition [20,21]. These decompositions are useful to extract the predominant spatio-temporal structures of the field quantities for given data. Wavelet transform is also utilized to decompose turbulent systems, and some advantageous features are reported in comparison to the Fourier decomposition or the POD [22,23]. However, data-driven analyses such as the SVD or POD are still important in plasma systems, where global structures such as zonal flows or streamers are spontaneously generated, and it is difficult to predict what is the dominant structure a priori.
The advantageous point of the data-driven modal analysis compared to the conventional Fourier decomposition is wide applicability to inhomogeneous systems with strong gradients and/or complex geometry, as well as homogeneous fields with the periodic boundary treatment. Thanks to these characteristics, the above data-driven mode decomposition has been utilized in fluid and plasma flow dynamics, and other research areas. However, the application to strongly developed turbulent fields with higher-order correlation quantities such as the turbulent transport flux has not been established yet.
In this work, we apply the SVD to the two-dimensional plasma turbulence described by a fluid model called the Hasegawa-Wakatani equation [24,25]. As are many other fluid or turbulence systems, the Hasegawa-Wakatani equation consists of multi-component fields and nonlinear quantities such as the turbulence-driven particle flux and the triad energy transfer to the zonal flows are of quite an interest. However, the conventional SVD cannot sufficiently analyze these nonlinear quantities, because it possesses the information of only one field quantity. Then, in order to analyze the nonlinear quantities in such a multi-component turbulent system, the conventional SVD is extended to decompose multiple turbulence fields simultaneously by a single set of orthonormal basis functions without imposing a priori scale separation. This extended method is referred to as the multi-field SVD and provides us with a novel modal analysis for nonlinear quantities in turbulence. We will see that the present multi-field SVD can successfully extract both the mode amplitude and the phase relation for the multiple turbulence fields, by comparison with the conventional Fourier decomposition. The multi-field SVD enables the reduction of the dimensionality in the triad energy transfer, where the sign of the SVD-mode spectra is separated or clustered. These features are due to the inclusion of phase relation between multiple fields in the multi-filed SVD, not incorporated in the conventional SVD. Then, the spatial structures of the SVD modes to dominantly contribute to the particle flux and the triad energy transfer are clearly identified.
The remainder of the paper is organized as follows. In section 2.1, the algorithm of the conventional and the present multi-field SVD are explained. The Hasegawa-Wakatani equation and the numerical conditions in the turbulence simulation are described in section 2.2. Then, the multi-field SVD is applied to the turbulent particle flux and the triad energy transfer in the plasma turbulence in section 3, discussing some advantageous features of the multi-field SVD. Finally, the concluding remarks are given in section 4.

Multi-field SVD
The mathematical definition of SVD is briefly summarized as follows. Let F be a general real or complex M × N matrix. Then SVD is given by where the complex transpose (V T ) * is represented by V † . Here, Σ = diag(s 1 , . . . , s n ) is the M × N rectangular diagonal matrix containing the non-negative real number called 'singular value' in the diagonal components. The number of non-zero diagonal components n is given by n = min(M, N). U and V are M × M and N × N unitary matrices, respectively. When F is real, V † = V T . Let u and v be orthonormal vectors for U and V, and the SVD for F is given by In this study, SVD is applied to the third-order-tensor-like grid data of turbulent fields, where a time-dependent twodimensional field quantity f(x, y, t) is expressed as the discretized form of f(x i , y j , t k ), where 1 ⩽ i ⩽ N x , 1 ⩽ j ⩽ N y , 1 ⩽ k ⩽ N t means the indices for spatio-temporal grids. Then f(x i , y j , t k ) is rearranged to a time-dependent vector f(t k ); Subsequently, a (N x N y ) × N t matrix F is constructed as for the temporal evolution of f(x, y, t). Now we can apply SVD to equation (4), and the time-dependent two-dimensional field quantity f(x, y, t) is expressed by a series of the orthogonal functions characterized by the singular value; where the continuous limit of f( stands for the total number of the SVD modes. Each column of the matrix U corresponds to the spatial structure of each mode and we denote it as ψ i (x, y) so that ψ i (x j , y k ) = U j+Ny(k−1),i . Similarly, each column of the matrix V corresponds to the temporal evolution of each mode and we denote it as h i (t) so that h i (t j ) = V j,i . Such SVD with preconditioning for tensors is referred to as 'higher-order SVD (HOSVD)' (see e.g. reference [26]). From the orthonormality of the eigenvectors in U, the spatial structures ψ i satisfies the orthonormal condition, i.e.
In order to analyze the multiple fields correlated with each other, the above higher-order SVD is extended to decompose multiple fields simultaneously. Let f 1 (x, y, t), f 2 (x, y, t), . . . , f M (x, y, t) be M fields of interest, and they are rearranged to M time-dependent vectors f 1 (t), f 2 (t), . . . , f M (t) in the similar manner to equation (3). These vectors are aligned in the temporal direction of the matrix F on an equal footing as shown in figure 1; Then, the matrix F is decomposed in the same procedure as the 1-field case like equation (1). As a result, each field labeled by l (1 ⩽ l ⩽ M) is decomposed as where the total number of the SVD modes N is determined by N = min(M × N t , N x × N y ). As seen in this equation, the singular value s i and the spatial structure ψ i (x, y) are common to all fields, so that the correlation properties such as the phase difference are reflected in the coefficient of the temporal evolution h (l) i . Here, this extended method is referred to as 'multifield (higher-order) SVD'. The orthonormality of the spatial structures, equation (6), is preserved in the multi-field SVD since the different fields are aligned in the temporal direction. Therefore, the second-order quantity of the fields is also decomposed as It should be stressed that any scale separations are not imposed to the basis spatial structure ψ i (x, y) in the multi-field SVD, unlike the Fourier decomposition.

Hasegawa-Wakatani equation
The microscopic drift-wave turbulence and the mesoscopic zonal flows, which are spontaneously generated via the nonlinear turbulent interactions, are often modeled by the (modified) Hasegawa-Wakatani equation [24,25]: The fundamental field variables are the vorticity ζ, the electron density n, and the electrostatic potential ϕ. Here, is the Laplacian in the two-dimensional rectangular geometry, and the curly bracket denotes the Poisson bracket {A, B} = ∂ x A∂ y B − ∂ y A∂ x B.Ã = A − ⟨A⟩ is the non-zonal component of the quantity A, where ⟨A⟩ = 1/L y´d yA is the zonal average of A, and L y is the box size in the y direction. The constant parameter κ is proportional to the equilibrium density gradient and the driving source of the drift-wave instability. α is referred to as the adiabatic parameter and controls the kinetic response of electrons to the non-zonal electrostatic potential. The dissipation coefficient is denoted by D. The physical quantities in the Hasegawa-Wakatani equation are normalized as t = ω cit , x =x/ρ s , ϕ = eφ/T e , n =n/n 0 , where the overline symbol stands for the dimensional quantity. ω ci and n 0 are, respectively, the cyclotron frequency of ions and the equilibrium electron density. ρ s = T e /mω −1 ci represents the ion sound Larmor radius, where T e is the electron temperature.
The Hasegawa-Wakatani equation is numerically solved in the two-dimensional rectangular geometry, where the periodic boundary condition is applied to both x and y directions. The Fourier spectral method is used for calculating the spatial derivatives and the fourth-order Runge-Kutta method is used for time integration. The numerical and physical parameters are set as N x = N y = 100, L x = L y = 2π/0.15 ≃ 42, κ = 1, α = 1, and D = 10 −4 . Note that although a moderate spatial resolution is chosen here to suppress the matrix size for SVD, one can generally use higher spatial resolutions.
The simulation results are summarized in figure 2. The temporal evolution of the kinetic energy of the turbulent (nonzonal) component E Turb = 1/(L x L y )´dx|∇φ| 2 and the zonalflow component E ZF = 1/(L x L y )´dx(∂ x ⟨ϕ⟩) 2 are shown in figure 2(a). The zonal flow grows after the excitation of the turbulence and eventually dominates the quasi-steady state under the physical parameters chosen here. It is stressed that an important feature of the plasma turbulence governed by the Hasegawa-Wakatani equation is the existence of the turbulent transport to be expressed in detail in section 3.1. As shown in figure 2(b), the particle flux is enhanced when the turbulence develops, and then it is suppressed during the quasi-steady state dominated by the zonal flow. Figures 2(c) and (d) show the snapshot of ϕ and n, respectively. Note that the electrostatic potential ϕ exhibits a significant zonal structure, whereas the more turbulent structure appears in the electron density n. The time-averaged spectrum of the zonal-flow component of the potential, ||ϕ(k x , k y = 0)|| 2 , is shown in figure 2(e). One can see that the zonal flow is dominated by several modes around k x ∼ 0.1 and the power-law-like decay is observed in higher-k x modes. The time-averaged spectrum of the turbulent electrostatic potential, kx ||ϕ(k x , k y )|| 2 , is also shown in figure 2(f).

Multi-field SVD for the turbulent particle flux
In turbulent plasmas, the particle transport is driven by turbulent convections due to the particle drift of v = E × B/B 2 = B × ∇ϕ/B 2 . Thus, the particle flux in the x direction Γ n is expressed as where the direction of the magnetic field is assumed to be the z direction in the present situation. It should be noted that the local turbulent particle transport can generally occur both inwardly or outwardly, i.e. (∂ y ϕ)ñ < 0 or (∂ y ϕ)ñ > 0, depending on the phase relation between the two turbulent fields. First, we apply the conventional 1-field SVD to the field The dataset is obtained with the time interval of ∆t = 0.1 between t = 1500 and 2500 and thus the number of time points is N t = 10 000. Since the spatial grid numbers are N x = N y = 100, the total number of the SVD modes is N = min(N t , N x × N y ) = 10 000. The mode index is aligned in the descending order of the singular value as shown in figure 3(a). Up to the mode index of 4422, the spectrum of the singular value is continuous and indicates powerlaw-like distributions, so the extraction of a few modes to make a dominant contribution is not trivial. This clearly exemplifies the limited characteristic of the conventional SVD for the quantities in strongly developed turbulence, in contrast to the fact that global laminar flow can be represented by fewer SVD modes [15]. One also finds that the modes after 4422 can be neglected by the smallness of their singular values in comparison to those for the preceding modes. The jump of the singular value is associated with the upper bound of effective modes to represent the field quantity due to the de-aliasing in the numerical simulations. In fact, in the present Fourier spectral scheme with 100 × 100 modes, the number of the effective modes is (100 × 2/3) 2 ≃ 4444, which is almost identical to the mode index indicating the singular value jump. We also confirmed that the position of the jump scales with the effective mode number.
Due to the orthonormality of the spatial structures, equation (6), the fieldñ is also decomposed by taking the inner product with the basis function ψ i (x) obtained from the 1-field SVD of ∂ y ϕ(x, t); where ⟨A, B⟩ =´dx A B is the inner product of the fields A and B. Then, the particle flux can be decomposed as   The spectra of (a) the singular value s i and (b) the particle flux Γ n,i in the 1-field SVD for ∂yϕ.
Next, the particle flux is analyzed with the multi-field SVD. From the definition equation (13), 2-field SVD for ∂ y ϕ andñ are considered here, where the 2-field matrix is given by Then the fields are decomposed into The dataset of N t = 5000 for the time interval of ∆t = 0.1 between t = 1500 and 2000 is used. Consequently, the total number of the SVD modes N is given by N = min(2N t , N x × N y ) = 10 000, where the factor 2 of N t results from the 2-field alignment in the temporal direction in the 2-field matrix F. As shown in figure 5(a), we confirmed the power-law-like singular value spectrum, which is similar to the case of the 1-field SVD. Then the singular values for the modes after the mode index of 4422 are negligibly small. In the multi-field SVD, the obtained spatial structures ψ i (x) include the contributions from both ∂ y ϕ andñ. In order to evaluate the reproducibility of the original field, the convergence with respect to the mode index is compared to that in 1-field SVD, where the cumulative root-mean-squared error is defined by As shown in figure 5(b), the deviation in the 2-field SVD is slightly larger than that in the 1-field SVD for higher mode indices above 3000 with the value of O(10 −6 ), indicating less  impact on the reproducibility. On the other hand, the deviations in the two cases are almost identical until the mode index around 3000. The original field is, thus, reproduced with a sufficient accuracy even for applying the multiple fields in the SVD. Substituting equations (17) and (18) to equation (13), the 2-field SVD of the particle flux is given by The time-averaged particle flux of each SVD mode Γ n,i is shown in figure 6(a), where Γ n,i > 0 and Γ n,i < 0 are simultaneously displayed in red and blue symbols, respectively. In contrast to the case of the 1-field SVD, one can see that the distribution of Γ n,i is clearly separated into the positive-flux and the negative-flux regimes, i.e. Γ n,i > 0 for 1 < i < 1500 and Γ n,i < 0 for 2000 < i < 3500. The separation is attributed to that the information on the phase difference between the fields is incorporated into the multi-field SVD. The partial summations for the positive and negative fluxes are evaluated as 4.38 × 10 −2 and −1.82 × 10 −4 , respectively. The importance of this analysis is demonstrating the clear separation of the turbulent particle flux, although the negative flux is negligibly small in this situation. In fact, the coexistence of positive and negative turbulent fluxes is observed in multi-species gyrokinetic simulations [27]. The multi-field SVD is straightforwardly applicable to more complicated gyrokinetic and fluid turbulence even for non-uniform global systems in general geometries, beyond the present simplified fluid system. Then, the separation of the positive and negative turbulent fluxes by the multi-field SVD will provide deeper understanding of physical mechanisms. The temporal evolution of the coefficients h is well separated. From 2-field SVD, the physical characteristics in the particle transport can be seen from the spatial structure ψ i (x, y). As the representatives of the modes with the positive and negative fluxes, the spatial structure of mode 1, mode 5, and mode 2000 are shown in figures 7(a)-(c), respectively. In the spatial structure ψ 1 (x, y), large-scale eddies exist at −20 < x < −10 and 10 < x < 20 with small-scale eddies within them, whereas small-amplitude vertical structures are shown in −5 < x < 10. It is the advantageous characteristic of the multi-field SVD that the mode that dominates the particle flux is characterized by the dominant multi-scale spatial structure, which can not be represented by only a few monochromatic modes in the Fourier decomposition. On the other hand, the near-isotropic fine-scale eddies dominate the spatial structure ψ 2000 (x, y). This result that the particle flux is dominated by large-scale vortices is consistent with the experimental observations [28].
The results of the multi-field SVD are compared with the Fourier decomposition. The fields are decomposed into and then the particle flux is also decomposed into The time-averaged particle flux Γ n in wavenumber-space is plotted in figure 8. It is found that the Fourier modes with small wavenumbers dominantly contribute to the positive flux. However, it is difficult to find a clear separation of the positive and negative fluxes in the time-averaged twodimensional wavenumber spectrum. Indeed, the positive and the negative contribution of Γ n at each time and at each point in wavenumber-space are calculated as 4.44 × 10 −2 and −7.69 × 10 −4 , where the probability density function (PDF) is shown in figure 8(b). Note that they match the multi-field SVD case in the order estimation. From these comparisons, it is demonstrated that the multi-field SVD provides one with an useful modal analysis to separate or extract the characteristic field structures in the turbulent transport. It is also emphasized that the multi-field SVD can apply to inhomogeneous turbulence systems that can not be decomposed by the Fourier modes.
Since the multi-field SVD utilizes the time evolution data of the field, the dependence of the time points N t on the convergence of the mode spectrum should be clarified. As shown in figure 9, the negative contribution is unphysically large when N t is small. As increasing N t , the mode spectrum converges. Therefore, a sufficient number of time points are necessary in order to analyze the turbulence and transport correctly.

Multi-field SVD for the triad energy function
Next, the multi-field SVD is applied to the triad energy transfer function [29,30] S p,q k that represents the nonlinear interactions among turbulent vortices and zonal flows. The three fields, the zonal component of the potential ⟨ϕ⟩, the turbulent component of the potentialφ, and the vorticity ζ are aligned in the matrix F as Note that although the turbulent and the zonal parts can be decomposed by the SVD for the original potential field ϕ, we separately treat ⟨ϕ⟩ andφ as different fields to express the energy transfer from/to the zonal flow explicitly. Then, the multi-field SVD is applied to F and the fields are decomposed, where N t = 3400 for the time interval of ∆t = 0.1 between t = 1500 and 1840. The total number of the SVD modes is N = min(3N t , N x × N y ) = 10 000. The kinetic energy of the zonal-flow component E ZF is decomposed into the SVD mode as expressed in equation (9);  where the fields are decomposed into ⟨ϕ⟩ = k ⟨ϕ⟩ k = k s k h (ϕ) k (t)ψ k (x, y) and so on. Thus, by multiplying −1/(L x L y )⟨ϕ⟩ k to the Hasegawa-Wakatani equation (10) and integrating over the whole space, the temporal evolution of E k is expressed as where Note that the nonlinear triad energy transfer S p,q k is symmetric for superscripts p and q. From equation (25), the positive Figure 9. The spectrum of the time-averaged particle flux for Nt = 50 (red square), 100 (green circle), 500 (blue triangle), 2500 (magenta inverted triangle), and 5000 (cyan diamond). The horizontal axis is normalized to the range of [0, 1] and the vertical axis is normalized so as for the maximum value to be 1.
value of S p,q k corresponds to the net energy transfer from the turbulent SVD modes p and q to the zonal SVD mode k, i.e. the zonal flow production, and the negative S p,q k corresponds to the opposite direction, i.e. zonal-flow decay. Similar to the case of the 1-field and the 2-field SVD, we confirmed that the modes after the index of 4488 are negligible, since the singular value s k decreases by about 7 orders there.
The contour of the time-averaged triad energy transfer in the SVD mode space is shown in figure 10(a), where S p,q k is normalized by the linear driving term κΓ n . Here, the mode k is set to k = 1 which represents the zonal flow as shown in figure 12(a). As is expected by the definition, the symmetric structure in p ↔ q is observed. One finds that the SVD modes driving the positive or negative energy transfer to the zonal mode with k = 1 are simply separated in the SVD mode space.
For comparison, the time-averaged triad energy transfer S p,q k in the Fourier decomposition is plotted in figure 10(b), where the definition is given by Note that the triad condition with δ k+p+q,0 in summation operation (omitted here) is used in the second equality. It is seen that the modes with small p y indicate a positive contribution and the modes with large p y have a negative contribution to the energy transfer to the zonal mode with k = (0.6, 0). As will be shown in the following, this is consistent with the result in the multi-field SVD. It is stressed that S p,q k with the 3field SVD significantly simplifies the triad energy transfer analysis in comparison to the conventional Fourier decomposition. Indeed, as shown in figure 2(e), the zonal flow is represented by many Fourier modes including the components decaying by a power law in the Fourier decomposition. Although we selected a single wavevector k = (0.6, 0) as the representative component of zonal flows, one needs to analyze S p,q k for many zonal modes to obtain the whole structures. On the other hand, in the SVD, only a few modes around k = 1 or only one mode summing up these few modes can express the zonal flow. This provides us with a great information reduction for analyzing S p,q k , but keeping key physical pictures. In order to clarify the physical meanings of the SVD modes, we introduce quantities to characterize the fine-scale structures and their anisotropy of ψ i (x, y). For the x-direction, it is defined by For the y-direction ⟨∂ 2 y ⟩ is also defined in a similar way. Then, the anisotropy is evaluated by ⟨∂ 2 y ⟩/⟨∂ 2 x ⟩, i.e. the relatively smaller ⟨∂ 2 y ⟩/⟨∂ 2 x ⟩ means an elongated structure in the y-direction. The mode-index dependence of ⟨∂ 2 y ⟩/⟨∂ 2 x ⟩ of the present 3-field SVD is shown in figure 11. One can see the local maxima and the local minima around mode 850 and 1500.
The spatial structures of the SVD basis function ψ 850 and ψ 1506 are shown in figures 12(b) and (c), respectively, as the representative of the negative and positive contributions of S p,q k in figure 10(a). The spatial structure with S p,q k > 0, ψ 1506 , consists of vortices elongated in the y direction, whereas smallerscale vortices dominate the other case with S p,q k < 0, ψ 850 . By examining the modes around 1500, we found that, in the quasisteady state, the zonal flow with k = 1 is sustained by the balance between the net energy gain from vortices elongated in the y direction and the net energy loss to smaller-scale vortices. In the Fourier decomposition, shown in figure 10(b), the positive contribution of S p,q k exists at the low-k y region, i.e. the zonal-flow mode k = (0.6, 0) gains energy from vortices elongated in the y direction. Therefore, the qualitative physical characteristics of the zonal-flow energy transfer are consistent between the multi-field SVD and the Fourier decomposition.
Another usefulness of the present multi-field SVD is demonstrated in the context of dimensionality or information reduction. Since the positive and negative S p,q k is simply separated and aligned along the p-axis, one can integrate S p,q k by  q, i.e. q⩽p S p,q k to express the energy transfer to the mode k by the nonlinear interaction including the mode p. Figure 13 shows the one-dimensional SVD spectrum of q⩽p S p,q k , where the positive and the negative contribution are clearly separated. Such advantageous characteristics of the multi-field SVD enable us to analyze the nonlinear dynamics in turbulent systems from a different aspect than the conventional Fourier decomposition. Indeed, two or more dimensional analysis is necessary to examine the whole wavenumber-space structures of S p,q k with complicated positive-negative boundaries (shown by solid curves in figure 10(b)). It should also be emphasized that the multi-field SVD is applied to the inhomogeneous turbulence as well.

Summary
In this paper, an extended mode decomposition method, multifield SVD, to analyze the multi-component turbulence is proposed by extending the conventional SVD or higher-order SVD (HOSVD). The multi-field SVD provides us with a simultaneous decomposition of several turbulent fields with a single set of the orthonormal basis functions ψ i (x, y). Then, in addition to the mode amplitude labeled by the singular value, the phase relations in the nonlinear quantities such as a transport flux or a triad energy transfer (or more general higherorder correlations) are extracted in the mode space.
In order to analyze the turbulent particle flux Γ n and the energy transfer between the zonal flows and turbulent vortices S p,q k , the two-dimensional multi-component plasma turbulence described by the Hasegawa-Wakatani equation is examined and the multi-field SVD is applied. Since the phase relation between multiple fields is included in the multi-field SVD, the sign of the SVD-mode spectra is clearly separated for both the particle flux and the triad energy transfer. Furthermore, the spatial structures of the dominant modes are clarified, preserving the multi-scale nature of the original turbulent fields. Moreover, it is demonstrated that one can reduce the dimensionality or information using the multi-field SVD through comparisons with the conventional Fourier decomposition. This implies many capabilities of the data-driven modal analysis of various nonlinear quantities in plasma turbulence. Also, the applicability to general inhomogeneous turbulence systems with the background gradients and/or the mean flows is emphasized.
As the future aspects, the multi-field SVD can be applied to more practical issues in the fusion plasma turbulence, which is essentially a multi-component and multi-field system. In particular, a strong coupling between turbulent energy and particle transport should be clarified in multi-ion mixture plasmas [27,31,32]. There is a recent application to the experimental data analysis [33]. Also, the nonlinear stability of zonal flows and the associated energy transfer processes in fluid and kinetic turbulence [34][35][36] need to be elucidated. Moreover, the inhomogeneity of the obtained spatial structure of the SVD modes is suitable to investigate turbulence trapping due to the zonal flows [36,37]. Future works will address these issues by flexibly utilizing both the Fourier decomposition and the multi-field SVD.

Data availability statement
The data cannot be made publicly available upon publication because no suitable repository exists for hosting data in this field of study. The data that support the findings of this study are available upon reasonable request from the authors.