This paper explores the fundamental characteristics of electron-temperature-gradient (ETG)-driven turbulence in the tokamak pedestal. The extreme gradients in the pedestal produce linear instabilities and nonlinear turbulence that are distinct from the corresponding ETG phenomenology in the core plasma. The linear system exhibits multiple (greater than ten) unstable eigenmodes at each perpendicular wave vector, representing different toroidal and slab branches of the ETG instability. Proper orthogonal decomposition of the nonlinear fluctuations reveals no clear one-to-one correspondence between the linear and nonlinear modes for most wave vectors. Moreover, nonlinear frequencies deviate strongly from those of the linear instabilities, with spectra peaking at positive frequencies, which is opposite the sign of the ETG instability. The picture that emerges is one in which the linear properties are preserved only in a narrow range of k-space. Outside this range, nonlinear processes produce strong deviations from both the linear frequencies and eigenmode structures. This is interpreted in the context of critical balance, which enforces alignment between the parallel scales and fluctuation frequencies. We also investigate the nonlinear saturation processes. We observe a direct energy cascade from the injection scale to smaller scales in both perpendicular directions. However, in the bi-normal direction, there is also nonlocal inverse energy transfer to larger scales. Neither streamers nor zonal flows dominate the saturation.

## I. INTRODUCTION

The electron-temperature-gradient instability is a microinstability present in fusion plasmas at electron (gyroradius) scales. It is responsible for driving microturbulence that leads to transport at small scales. One manifestation of this instability is of a class of “bad curvature” instabilities—akin to the Rayleigh–Taylor instability. This curvature-driven ETG mode has been the most widely studied because it produces transport in the core of tokamak plasmas. Here, we study turbulence driven by the slab version of the ETG instability, which is manifest in edge transport barriers where the temperature gradients can become extremely steep.

Historically, ETG turbulence was thought to be negligible due to the low transport levels predicted by simple mixing length estimates: $ \chi e = \rho e 2 v e , t h / L T e$, which is lower than the ion-scale counterpart by a factor of the square root of the mass ratio (*χ _{e}* is the electron thermal diffusivity,

*ρ*is the electron gyroradius, $ v e , t h$ is the electron thermal velocity, and

_{e}*L*is the electron temperature gradient scale length). This notion was dispelled by nonlinear gyrokinetic simulations,

_{Te}^{1–3}demonstrating that the naive mixing length estimate was violated by radially extended structures, called streamers, that form in the nonlinear turbulent state, thus establishing the potential for ETG as a major transport mechanism. Subsequent work observed saturation via inverse energy transfer and demonstrated that transport levels are independent from streamer scale lengths, calling into question the mixing-length conceptualization of ETG turbulence.

^{4}

Recent research has discovered another scenario (not relying on streamers) in which ETG can contribute to tokamak transport. In the tokamak H-mode,^{5} a transport-barrier forms at the plasma boundary, resulting in a region with steep temperature and density gradients and a corresponding dramatic boost to plasma confinement. This region of steep gradients, located at the edge of the plasma (radially outward from the core), is known as the pedestal. ETG turbulence in the pedestal is largely slab-like, which typically produces fluctuations that are much more isotropic in comparison with core, streamer-dominated, turbulence.^{6,7} In this system, the transport is indeed quite low, scaling roughly as the gyroBohm prediction.^{8} In the pedestal, however, this does not relegate ETG to irrelevance; low transport is, after all, the hallmark of a transport barrier. In fact, pedestal ETG turbulence has been demonstrated in multiple studies to produce experimentally relevant transport levels.^{6,9–15} This occurs due to the extreme pedestal gradients (larger than core values by an order of magnitude) compensating for the low diffusivity: $ Q = n \u2207 T \chi $. Transport can be increased further due to “stiffness”: the propensity for transport to increase rapidly as the driving gradient surpasses its threshold. Consequently, ETG can be a major transport mechanism in the pedestal.

Pedestal ETG turbulence exhibits several unique characteristics. The magnetic drift resonances that are typical of the core are weakened by the discrepancy between scale lengths of magnetic field in comparison with those of the density and temperature. This high-gradient scenario favors^{16,17} slab resonances,^{18} which result in high parallel wavenumbers (far surpassing the inverse connection length $ 1 / q R$) and demand extreme resolution in the parallel direction for numerical studies. The slab-like nature of the turbulence has further implications: instabilities with little or no sensitivity to ballooning angle, producing much more isotropic fluctuations and many unstable eigenmodes (representing a spectrum of $ k \u2225$) at each wavenumber. The toroidal drift resonances that do persist are also unique, destabilized at the high radial wavenumbers necessary to compensate for the relatively small magnetic scale lengths.^{19,20}

Although several studies have established experimental relevance, less has been done to study the fundamental nature of turbulence in this new regime. Questions like the following remain unanswered: (1) What structures persist in the nonlinear turbulent state? (2) What kind of connections exist between the linear eigenmodes and these nonlinear structures? (3) In what manner do the fluctuations self-organize to achieve a saturated state—i.e., how is energy injected, redistributed among scales, and ultimately dissipated? Such investigations are the major goal of this work.

We find that at the most unstable wavenumbers, the primary unstable modes in the linear system are highly ballooned (i.e., localized at the outboard midplane) slab modes. However, in the nonlinear system, there are significant contributions to the heat flux from multiple nonlinear modes, and proper orthogonal decomposition (POD) reveals no clear correspondence between the linear and nonlinear modes for most wave vectors. The nonlinear structures that develop in the turbulent system have very fine parallel structure and high frequencies well beyond the linear frequencies. This is true throughout k-space—even in the peak drive range of the system. Notably, the frequency spectrum often peaks at positive frequencies for many wavenumbers even though the ETG frequency is negative. We find that the linear properties persist only in a narrow range of k-space and even then, only for a single POD mode among many modes with significant amplitude. We interpret these observations in the context of critical balance theory, which enforces alignment between parallel scales and fluctuation frequencies.

Investigation into the nonlinear saturation process shows that energy is injected at low radial wavenumbers and high bi-normal wavenumbers. From this injection scale, the energy is directly transferred to small scales in the radial direction. In the bi-normal direction, we identify both direct (to small scales) and nonlocal (inversely to large scales^{4,21}) transfer.

Overall, this study provides insights into the unique characteristics of ETG turbulence in the tokamak pedestal and its nonlinear behavior. In addition to the fundamental scientific interest, this study also informs the construction of reduced models, which are being formulated for this important transport mechanism.^{8,12,14}

The paper is organized as follows: Section II describes the equations solved and the numerical methods used. Section III presents results from linear initial value simulations and eigenvalue calculations. Section IV presents results from nonlinear initial value simulations and compares them to those from Sec. III. Section V concludes with a discussion of the results. Finally, two Appendixes contain definitions used in the paper ( Appendix A) and a brief explanation of the proper orthogonal decomposition used in this paper ( Appendix B).

## II. METHOD

In the strong-guide-field plasma of a tokamak, the gyrokinetic model^{22} is appropriate^{23} and is employed here. The gyrokinetic code, gene,^{2} is used to numerically solve the gyrokinetic equation. The small scales characteristic of ETG turbulence are additionally amenable to the local flux-tube approximation with adiabatic ions^{24} and kinetic electrons, which allows the radial and bi-normal directions of the domain to be decomposed into a Fourier basis and periodic boundary conditions to be applied.^{25}

*j*may be written as

*g*, which is sometimes called the modified gyrocenter distribution function, is related to the total first-order fluctuating distribution function,

*f*, as $ g j ( k x , k y , z , \mu , v \u2225 ) = f j ( k x , k y , z , \mu , v \u2225 ) + 2 q j m j v T j v \u2225 A \xaf \u2225 j F 0 j$, where

*k*is the radial wavenumber,

_{x}*k*is the bi-normal wavenumber,

_{y}*z*is the parallel (to the magnetic field) coordinate (parameterized by the poloidal angle),

*μ*is the magnetic moment, and $ v \u2225$ is the parallel velocity. The linear terms are

^{26}Curvature terms in Eq. (2) are $ K x = \u2212 a B 0 ( \u2202 y B 0 + \gamma 2 \gamma 1 \u2202 z B 0 ) , \u2009 K y = a B 0 ( \u2202 x B 0 \u2212 \gamma 3 \gamma 1 \u2202 z B 0 ) , \u2009 \Gamma x , y = \u2202 x , y g + \u2212 e T 0 e F 0 \u2202 x , y \chi $, and $ \Gamma z = \u2202 z g + \u2212 e T 0 e F 0 \u2202 z \chi + \u2212 e v T e T 0 e v \u2225 \mu F 0 A \u2225 \u2202 z B 0$, where $ \gamma 1 = g 11 g 22 \u2212 g 21 g 12 , \u2009 \gamma 2 = g 11 g 23 \u2212 g 21 g 13 , \u2009 \gamma 3 = g 12 g 23 \u2212 g 22 g 13$, and $ g \alpha \beta $ are metric coefficients. The distribution function and fields are normalized as $ g \u2192 \rho s a n 0 e 2 c s 3 v T e 3 g , \u2009 F 0 \u2192 n 0 e 2 c s 3 v T e 3 F 0 , \u2009 \Phi \u2192 \rho s a T 0 e e \Phi $, and $ A \u2225 \u2192 \rho s a B 0 \rho s A \u2225$.

The equilibrium parameters for the magnetic field *B*_{0}, temperature $ T e 0$, and density $ n e 0$, are summarized in Table I, along with the scale length *a* and reference (ion) mass *m _{i}*. Other physical parameters are summarized in Table II, and they are the radial location

*x*

_{0}, safety factor

*q*

_{0}, magnetic shear $ s \u0302$, density gradient $ \omega n \u2261 a / L n$, electron temperature gradient $ \omega T \u2261 a / L T$, ratio of gradients $ \eta = \omega T / \omega n$, and plasma $ \beta = 2 \mu 0 n e 0 T e 0 / B 0 2$. The reference velocity is the sound speed $ c s = T e 0 / m i$, so that time is in units of $ a / c s$. Note that values reported in this article, where not explicitly stated, will be non-dimensionalized with respect to these reference values. For example, frequencies will be in units of $ c s / a$, and the heat flux will be given in gyroBohm units $ Q GB = c s \rho s 2 n e 0 T e 0 / a 2$. With $ \beta > 0$, electromagnetic perturbations are included in the model. For additional details on the numerical implementation, see, e.g., Told

^{7}and Merz.

^{27}

B_{0} | $ 2.03 \u2009 T$ |

$ T e 0$ | $ 3.67 \xd7 10 \u2212 1 \u2009 \u2009 keV$ |

$ n e 0$ | $ 4.14 \xd7 10 \u2212 19 \u2009 m \u2212 3$ |

a | $ 7.69 \xd7 10 \u2212 1 \u2009 \u2009 m$ |

m _{i} | $ 2.00 \u2009 m p$ |

B_{0} | $ 2.03 \u2009 T$ |

$ T e 0$ | $ 3.67 \xd7 10 \u2212 1 \u2009 \u2009 keV$ |

$ n e 0$ | $ 4.14 \xd7 10 \u2212 19 \u2009 m \u2212 3$ |

a | $ 7.69 \xd7 10 \u2212 1 \u2009 \u2009 m$ |

m _{i} | $ 2.00 \u2009 m p$ |

x_{0}
. | q_{0}
. | $ s \u0302$ . | $ a / L n$ . | $ a / L T$ . | η
. | β
. |
---|---|---|---|---|---|---|

0.973 | 4.57 | 3.22 | 22.72 | 47.12 | 2.07 | 0.14% |

x_{0}
. | q_{0}
. | $ s \u0302$ . | $ a / L n$ . | $ a / L T$ . | η
. | β
. |
---|---|---|---|---|---|---|

0.973 | 4.57 | 3.22 | 22.72 | 47.12 | 2.07 | 0.14% |

Initial value simulations (linear and nonlinear) and eigenvalue simulations are studied in this work.

Parameters were taken from DIIID discharge 162 940. This pedestal was also studied by Guttenfelder *et al.*^{12} and Hassan *et al.*^{28,29} As in these references, we adapt the background gradients to enhance ETG transport by increasing (decreasing) the electron temperature (density) gradient by 20%, which is within experimental uncertainties.

## III. LINEAR GROWTH RATES AND EIGENMODES

We begin this study with an investigation of the linear eigenmodes in the ETG range of the spectrum $ 1 < k y \rho s < 200$. The growth rates and frequencies from linear initial value (LIV) calculations are plotted in Fig. 1 along with the growth rates found from the gene eigenvalue (EV) solver.

We note the deep spectrum of unstable modes at each wavenumber. This is in stark contrast with the core plasma, where typically a small number of modes (one to three) is unstable at each wavenumber and contributes to the transport.^{1,30,31} This raises the question of how many modes are active in the nonlinear system, which will be addressed later in Sec. IV D.

The parallel mode structures of the electrostatic potential $\Phi $ of a sample of these most unstable eigenmodes is plotted in Fig. 2. There is a small peak of the ETG growth rate near $ k y \rho s \u2248 5$. These modes at small *k _{y}* ( $ k y \rho s \u2272 10$) appear to be toroidal ETG modes, described in detail in Ref. 19, which are characterized by relatively long parallel wavelengths, low

*k*wavenumbers, and large radial

_{y}*k*wavenumbers. Modes in a mid range of $ 10 \u2272 k y \rho s \u2272 120$ have finer

_{x}*z*structure characteristic of slab ETG modes (note the narrow range plotted in

*z*, where

*z*is the poloidal angle, which parameterizes the parallel direction). Modes with high

*k*( $ 120 \u2272 k y \rho s$) appear to also be slab in nature, with very small parallel scales adapted to minimize finite-Larmor-radius (FLR) effects, which rise precipitously along the field line at these wavenumbers: $ k x = z s \u0302 k y$.

_{y}Near the peak of the growth rate spectrum ( $ k y \rho s \u2248 128$), the eigenvalue solver reveals two mode branches with similar growth rates. See Fig. 3, which shows the first few eigenmodes for four *k _{y}* near the peak. The singly peaked mode that is subdominant at lower

*k*is ultimately the dominant branch for all $ k y \rho s \u2265 128$. While stable modes were not found with the eigenvalue solver at these wavenumbers, it is possible that the other branches stabilize at high

_{y}*k*.

_{y}These eigenmodes will be compared with the nonlinear structures extracted from the turbulence in Sec. IV.

## IV. NONLINEAR SIMULATIONS

We now turn to an investigation of the nonlinear state. This section begins with a description of system setup and convergence checks. Then an overview of the system in the nonlinear state is given by examining time traces and spectra. Subsequently, the production, transfer, and dissipation of energy is traced through the system followed by an examination of the nonlinear mode structures. Finally, in Sec. IV E, timescales and parallel scale lengths are discussed in the context of critical balance.

### A. Setup and convergence

For this study, the numerical parameters for the nominal setup were chosen as follows. The nominal box size was $ l x \rho s \u2212 1 = 0.776 \u2009 657$ and $ k y , min \rho s = 8$ in physical space^{32} and $ v \u2225 \u2208 [ \u2212 3 v T e , 3 v T e ]$ and $ \mu \u2208 [ 0 , 9 T e 0 / B 0 ]$ in velocity space. The number of grid points in each direction was *n _{x}* = 256, $ n k y = 64$,

*n*= 2048, $ n v \u2225 = 40$, and $ n \mu = 8$.

_{z}Energy is injected into the system by the electron temperature gradient term as defined in Eq. (8). The local flux tube formalism ensures that the driving gradients are fixed throughout the simulation time.

Convergence checks of the nonlinear setup were carried out in parallel resolution *n _{z}*, and box size

*l*and $ k y , min$. Additionally, velocity-space convergence was checked by increasing the resolution for the fastest-growing mode in the system. Shown in Fig. 4 are two measures of convergence with respect to parallel resolution: heat flux and fractional collisional dissipation.

_{x}^{33}is given by

*n*are the fluctuating moments. The hat denotes normalization such that the heat flux is in gyroBohm units. The fractional collisional dissipation is defined as $ D coll / ( D coll + D hypz )$, where $ D coll$ and $ D hypz$ are the dissipation due to collisions and parallel hyperviscosity, respectively.

The heat flux requires at least 1024 parallel grid points to resolve, as seen by the squares in Fig. 4. A parallel resolution of $ n z > 512$ is required for physical dissipation to dominate over a hyperviscosity (circles in the figure). While these extreme resolutions are not necessary for rough estimates of overall heat flux, they are important for accurately tracking the energetics.

Convergence with respect to the box size was checked by doubling the nominal box size of $ l x \rho s \u2212 1 = 0.776 , \u2009 k y , min \rho s = 8$ (case 1), resulting in case 6. The heat flux and fractional collisional dissipation for both of these setups are plotted in Fig. 4. The values are the same (overlapping in the figure), indicating that total heat flux is converged in box size. Chen *et al.*^{4} notes the requirement that ETG transport be converged in radial box size in order to justify the flux-tube approximation, a condition which is satisfied here. For these same two cases, the heat flux spectra are the same (not shown)—except for the addition of values at lower *k* because of the larger box—indicating convergence by this metric as well. A contour plot of $\Phi $ (Fig. 5) shows the larger of these two cases, with structures that are very fine scale in *y* and extended in *x*. Even though all other setups in this work are half as large in the *x* and *y* dimensions, the correlation lengths are smaller than the box dimensions in all cases. Specifically, they are $ ( C x , C y ) \rho s \u2212 1 = ( 0.12 , 0.013 )$. Notably, this turbulence is more anisotropic than observed in other ETG pedestal turbulence,^{6,7,14} with structures that are elongated in the radial direction ( $ C x / C y \u2248 9$).

Simulations were performed in both double and single precision, with no significant differences observed between the two using the measures described above (heat flux spectra, etc.) Therefore, the single-precision setups will be used for analysis because they have longer integration times for better steady-state statistics.

All setups analyzed in this work are summarized in Table III.

Case . | n
. _{x} | $ n k y$ . | n
. _{z} | $ l x \rho s \u2212 1$ . | $ k y , min \rho s$ . | Prec. . | $ t max c s a$ . | $ Q es / Q GB$ . |
---|---|---|---|---|---|---|---|---|

1 | 128 | 64 | 256 | 0.776 | 8 | S | 1.00 | 25.2 ± 0.1 |

2 | 128 | 64 | 512 | 0.776 | 8 | S | 1.25 | 28.1 ± 0.1 |

3 | 128 | 64 | 1024 | 0.776 | 8 | D | 0.28 | 31.4 ± 0.2 |

4 | 128 | 64 | 2048 | 0.776 | 8 | D | 0.53 | 31.9 ± 0.3 |

5 | 128 | 64 | 2048 | 0.776 | 8 | S | 0.68 | 32.3 ± 0.1 |

6 | 256 | 128 | 256 | 1.553 | 4 | S | 1.14 | 25.3 ± 0.1 |

Case . | n
. _{x} | $ n k y$ . | n
. _{z} | $ l x \rho s \u2212 1$ . | $ k y , min \rho s$ . | Prec. . | $ t max c s a$ . | $ Q es / Q GB$ . |
---|---|---|---|---|---|---|---|---|

1 | 128 | 64 | 256 | 0.776 | 8 | S | 1.00 | 25.2 ± 0.1 |

2 | 128 | 64 | 512 | 0.776 | 8 | S | 1.25 | 28.1 ± 0.1 |

3 | 128 | 64 | 1024 | 0.776 | 8 | D | 0.28 | 31.4 ± 0.2 |

4 | 128 | 64 | 2048 | 0.776 | 8 | D | 0.53 | 31.9 ± 0.3 |

5 | 128 | 64 | 2048 | 0.776 | 8 | S | 0.68 | 32.3 ± 0.1 |

6 | 256 | 128 | 256 | 1.553 | 4 | S | 1.14 | 25.3 ± 0.1 |

### B. Time traces and spectra

Time traces for several quantities are shown in Fig. 6, exhibiting an initial transient linear phase followed by a quasi-steady turbulent state. The heat flux saturates earlier ( $ t \u2248 0.1 a / c s$) than the moments, which saturate at $ \u2248 0.4 a / c s$. With $ t max \u2248 0.7 a / c s$ for the highest resolution run, it might appear that the averaging time $ \Delta t \u2248 0.3 a / c s$ is not sufficient. However, typical nonlinear correlation times $ \tau \u2248 10 \u2212 3 a / c s$, so that averages are taken over 30–50 nonlinear times.

*k*and $ k y \rho s \u2248 100$ (Fig. 7), very near the wavenumber at which the linear growth rate peaks (Fig. 1). The energy spectra, defined as

_{x}*k*spectrum, this is due to a greater number of energy-containing $ k x \u2260 0$ modes at

_{y}*k*close to the flux peak, as can be seen in Fig. 9. However, that figure also shows that the modes with the most energy are at lower

_{y}*k*than where the energy is injected.

_{y}The discrepancy between the peak heat flux wavenumbers and those of the fluctuation energy is related to the different saturation time scales observed in the time traces; the heat flux saturates early at small scales, and an inverse transfer drives the growth of fluctuation energy (and fluctuating moments) at larger scales on the slower timescale. The heat flux is driven directly by the linear instabilities. The low-*k _{y}* fluctuation energy is driven by an inverse transfer, which has no special inclination to enforcing a flux-producing phase relation. We present a more detailed discussion of the nonlinear saturation mechanisms below in the context of nonlinear transfer functions in Sec. IV C.

Considering the frequencies in the nonlinear system, shown in Fig. 8, it is notable that there are both positive and negative frequency branches. In that figure, the frequency spectrum of $\Phi $ is plotted for a range of *k _{y}*. The black line denotes the frequencies of the most unstable linear modes. The most unstable modes found both from LIV and EV simulations all have negative frequencies, which is expected from ETG modes. Nearly all of the subdominant eigenmodes also have negative frequency. The exceptions are modes with very low frequencies ( $ \omega a / c s < 10$) found at $ 80 \u2264 k y \rho s \u2264 104$. However, in Fig. 8, frequencies of the positive band are in the range of $ 100 \u2272 \omega a / c s \u2272 500$, clearly not these subdominant eigenmodes. Apparently, they do not correspond to any single eigenmode or set of linear eigenmodes. It will be seen in Sec. IV D that the part of the frequency spectrum at large, positive values of

*ω*corresponds to particular nonlinear modes.

### C. Saturation

**k**from interaction with all other modes. It is given by

^{34}

**k**).

Figure 9 displays the three dominant terms in Eq. (7), as well as *E* itself.^{36} Since the energy injection is proportional to the heat flux, the dominant flux scales approximately match the regions of high growth rates, as discussed above in Sec. IV. Nonlinear energy transfer is the main mechanism balancing the energy drive at these **k**. While Fig. 9 does not show the exact path of this energy transfer, general trends can be inferred. The net nonlinear transfer is from modes near $ k \rho s = ( 0 , 96 )$ to other $ k \u2032$. A particular net receiver of energy from the transfer are modes with the smallest, nonzero *k _{x}* and small

*k*[red region in Fig. 9(d)]. An apparent mismatch between the total positive and negative contributions in $ d E d t | NL$ is due to the broad transfer of energy to many high $ | k |$ modes (the range of

_{y}**k**is higher than shown in Fig. 9). One can infer from Fig. 9 that there is an inverse transfer

^{37}of energy in

*k*to wavenumbers smaller than the driving range. This is in addition to an apparent direct transfer in

_{y}*k*to wavenumbers higher than the driving range and a direct transfer to higher

_{y}*k*. Subsequent detailed analysis of the nonlinear transfer for individual wavenumbers will confirm this to be the case.

_{x}In contrast to the energy drive, the collisional dissipation is distributed across a broad range of scales in both *k _{x}* and

*k*and is mostly proportional to the energy, see Fig. 9. However, the rate of energy dissipation ( $ 1 E d E d t | coll$) at wavenumbers with the most energy is actually less than high

_{y}*k*,

_{x}*k*wavenumbers (not shown).

_{y}The substantial dissipation in the drive range is in stark contrast with the standard picture of fluid turbulence (e.g., high Reynolds number Navier–Stokes), which exhibits distinct scale separation between injection and dissipation. The lack of scale separation between drive and dissipation has been investigated in several studies. Hatch *et al.*^{38,39} note that, since the collision operator is primarily enhanced by small scales in velocity space rather than real space, dissipation can occur at all spatial scales, with a propensity toward smaller scales only as collisionality is decreased and/or nonlinear energy transfer is enhanced. In Refs. 31 and 40, this phenomenon is explained in terms of damped eigenmodes that are excited in the same scale range as the instability. See, e.g., Refs. 41 and 42 for this phenomenon in a fluid context.

In addition to the overall nonlinear transfer at some scale *k*, one can look for the individual wavenumbers involved in the nonlinear interaction. In the cases of other instabilities, particular modes have had important roles in saturation and steady state. For example, zonal modes^{43} (*k _{y}* = 0) have been found to play a dominant role in the saturation phase of the instability for core ITG (both toroidal and slab), but they have a weaker role in the steady-state energy transfer.

^{34,35}During saturation, they are net energy sinks relative to the non-zonal modes. While in steady state, the zonal modes mediate the transfer of energy from non-zonal modes at some

*k*to non-zonal modes at some higher $ k x \u2032$. On the other hand, in core toroidal ETG saturation, radial streamers (

_{x}*k*= 0) play the dominant role.

_{x}^{3,35}Others have seen the presence of zonal flows in the long-time evolution of ETG-driven systems,

^{44}but that is not observed here for a case that was run until $ t c s / a \u2248 6.2$ ( $ t v T e / a \u2248 514$). In the current work, we show that neither zonal flows nor streamers dominate the energy transfer. The energy transfer is rather complex, exhibiting a direct forward cascade in

*k*at large wavenumbers and a nonlocal inverse cascade in

_{y}*k*in the drive range. In

_{y}*k*, the energy transfer is to higher

_{x}*k*.

_{x}The energy balance at each wavenumber **k** is given by the sum of injection, dissipation, and nonlinear terms. Nonlinear transfer functions, defined in Eq. (11), are used to single out particular wavenumbers in the nonlinear redistribution process. Considering the nonlinearity in the perpendicular direction, the nonlinear transfer function $ T k , k \u2032$ gives the transfer of energy to modes with **k** from the interaction of two modes with $ k \u2032$ and $ k \u2033$ satisfying the three-wave interaction condition $ k \u2212 k \u2032 \u2212 k \u2033 = 0$.

Figure 10 shows the time-averaged nonlinear transfer functions $ \u27e8 T k , k \u2032 \u27e9$ plotted for select wavenumbers (*k _{x}*,

*k*). In that plot, the green dot marks the wavenumber

_{y}**k**labeled in the figure, and the color bar indicates whether

**k**is receiving energy from (positive values) or giving energy to (negative values) the mode $ k \u2032$.

^{45}It is worth noting that while overall the nonlinear transfer terms must sum to zero—they cannot create or destroy energy, only redistribute it—that is, an individual wavenumber may have a net positive or negative transfer of energy. i.e., $ d E k d t | NL \u2260 0$, but $ \u2211 k d E k d t | NL , = 0$.

Beginning at the high *k _{y}* = 128 range shown in the top row of Fig. 10, we see that the wavenumbers

*k*,

_{x}*k*(green dots) receive energy from a band of $ k x \u2032 , k y \u2032$ (the red region) just below

_{y}*k*,

_{x}*k*(the green dots). This is indicative of a relatively local direct cascade at high

_{y}*k*.

_{y}Moving now to the driving range of the instability, that is, **k** with $ | k x | \rho s \u2264 8.09$ and $ 72 \u2264 | k y | \rho s \u2264 96$, the energy transfer is more complex but can still be traced through the system. Examining Fig. 10, a picture forms of energy transported from the driving range to higher and lower $ k y \u2032$, both locally and non-locally, and to higher $ k x \u2032$ locally. For example, from $ ( k x , k y ) \rho s = ( 8.09 , 96 )$, the energy is predominantly transferred to a range of modes at a single $ k x \u2032$ ( $ = 16.2 \rho s \u2212 1$) but different $ k y \u2032$ ( $ 32 < k y \u2032 \rho s < 120$). In fact, one of the strongest transfers is non-locally in *k _{y}* to $ k \u2032 \rho s = ( 16.2 , 40 )$.

Through this non-local transfer, energy in the driving range can move to low *k _{y}* through only a few interactions. Once at low (

*k*,

_{x}*k*), the power budget is balanced between dissipation and the nonlinear transfer, which are comparable in amplitude. From here, the transfer is from modes with similarly low

_{y}*k*to those with higher

_{y}*k*. Primarily, at low

_{x}*k*, energy goes to $ k x \u2032 = k x + k x , min$. It is both local and direct (i.e., to $ k x \u2032 > k x$). It is only at

_{x}*k*ρ

_{x}_{s}≥ 48.5 and

*k*= 0 that energy transfer to $ k x \u2032$ significantly greater than

_{y}*k*is seen. This is seen in Fig. 11. As

_{x}*k*increases, there is a greater range of modes with $ k x \u2032$ that receive a significant portion of the energy from the mode $ ( k x , 0 )$.

_{x}At higher *k _{x}*, it is also notable that the amplitudes of the nonlinear transfer functions (Figs. 10 and 11) are reduced compared to the driving range. This is explained by the observation described above, where the dissipation is comparable to the nonlinear transfer at many of the wavenumbers through which the energy travels to arrive at high

*k*, low

_{x}*k*. In other words, there does not exist an “inertial range” in the classic sense of hydrodynamic turbulence; dissipation is present at significant amplitudes at all scales, sapping energy as it moves through the system.

_{y}Compared to the core ITG and ETG, e.g., Ref. 35, this pedestal system does not exhibit zonal modes or streamers playing a dominant role in steady-state energy transfer. What would be expected if radial streamers, for example, were important in energy transfer? Either they would receive a large share of energy compared to other non-streamer wavenumbers, or they (*k _{x}* = 0 modes) would represent the wavenumbers $ k \u2033$ that predominantly interact with

**k**to transfer energy to/from $ k \u2032$.

Consider the first scenario. In Fig. 10, it is observed that modes with $ k x \rho s = 0$ and $ k x \rho s = 8.09$ (at the same *k _{y}*) have similar magnitudes of energy transfer. This is also seen in the bottom-right subfigure of Fig. 9. Figure 10 shows that modes $ k = ( 0 , k y )$ in the driving range are responsible for moving energy to other higher- $ k x \u2032$ wavenumbers. As one moves to wavenumbers where the drive is less strong, then the

*k*= 0 modes gain as much energy as they lose, consistent with the fact that the drive and dissipation terms are nearly equal at these wavenumbers. Importantly, this behavior of $ k x = 0$ wavenumbers is not exclusive and is observed in wavenumbers with $ k x \u2260 0$ as well.

_{x}In the second scenario, one would expect to see a vertical banded structure in Fig. 10 at wavenumbers where $ k x \u2032 = k x$, indicating that $ k x \u2032 \u2032 = 0$ is the dominant wavenumber completing the triplet. This is clearly not the case, since the regions of high energy transfer (bright red or blue) are never vertically aligned with the green dot (i.e., the $ k x \u2032$ where the transfer peaks is never equal to *k _{x}*).

Similar arguments can be made about zonal flows by simply exchanging *k _{x}* and

*k*.

_{y}A direct cascade in *k _{x}* along with an inverse cascade in

*k*for ETG was also observed in Refs. 4 and 21, based on a combined gyrokinetic and analytic analysis. Similarly to that result, we do observe long wavelengths in the spectrum of $ \Phi ( k y )$ generated earlier than shorter wavelengths. Also, nonlinear transfer functions show interactions between driving modes with

_{y}*k*interacting with modes with $ k \u2032 y \u2248 k y$. However, some differences between the works warrant noting. First, in addition to the inverse energy transfer, we also see direct transfer for very large wavenumbers. Also, the previous work was in the context of streamers. Although the fluctuations are not perfectly isotropic, we do not observe streamers to dominate the saturation as discussed above. In the toroidal mode coupling work, three-wave resonant interactions are negligible. Here, three-wave coupling is an intrinsic part of the nonlinear saturation.

_{y}In summary, it is observed that for $ k y \rho s \u2273 100$, the energy is predominantly transferred to modes $ k \u2032$ with $ k \u2032 y > k y$, indicating a direct transfer of energy to smaller scales. For modes with $ k y \rho s \u2272 80$, the energy is predominantly transferred to modes with $ k \u2032 y < k y$, indicating the nonlinear transfer of energy to large scales. In *k _{x}*, a direct transfer to larger $ k \u2032 x$ is observed at all scales. While

*k*= 0 modes and

_{x}*k*= 0 modes are present in the steady state and important to the energy transfer, they do not dominate over other $ k x \u2260 0$ and $ k y \u2260 0$ modes, indicating that neither streamers nor zonal flows are dominant in the saturation.

_{y}### D. Nonlinear modes

^{46–49}The POD decomposes the nonlinear data into two orthonormal bases: one in time, another in space, ordered by their associated amplitude. The bases are optimal in the sense that they reproduce more of the original data at a given truncation level than any other possible basis of this form. By performing a POD in the time and

*z*domains, the data can be written as

*σ*are the singular values (POD amplitudes). Here,

_{i}*ξ*is a vector that contains data from multiple fields, $ \xi = [ \Phi \u0303 , T \u0303 \u22a5 , T \u0303 \u2225 , n \u0303 ] T$ such that $ \xi i ( z ) = [ \Phi \u0303 i , T \u0303 \u22a5 , i , T \u0303 \u2225 , i , n \u0303 i ] T ( z )$. This decomposition is performed at a particular

*k*and

_{y}*k*, with all connected

_{x}*k*,

_{x}^{50}extending the

*z*domain. See Appendix B for a detailed explanation.

*k*to account for the heat flux. Figure 12 shows the singular values

_{y}*σ*and electrostatic heat fluxes

_{i}*Q*for each POD mode

_{i}*i*. Also plotted is the cumulative fraction of each quantity for a given number of POD modes, defined as

*h*related to POD mode

_{i}*i*, where

*N*is the total number of POD modes. From the figures, it is clear that $ CF [ Q ] ( n ) > CF [ \sigma ] ( n )$, i.e., a greater amount of the total heat flux is contributed by a given number of modes than amount of the amplitude of the data. For example, 91% of the heat flux is contributed by the first ten POD modes, while they contribute only 34% of the total “energy” of the data. This is somewhat to be expected from the quadratic form of the heat flux, $ Q i \u221d \sigma i 2$. Even so, it is clear that several modes contribute to the heat flux at a given wavenumber; apparently, no single POD mode establishes dominance in the nonlinear state. This is in contrast with the typical core scenario. There, a small number of modes is unstable, and a small number of structures contribute to the heat flux.

^{31,40,51}

One might expect the high number of active modes to complicate quasilinear approaches to reduced modeling for this scenario. However, Hatch *et al.*^{8} demonstrate that an accurate reduced model can in fact be constructed from the most unstable eigenmode alone (albeit with an additional parameter dependence on $ \eta = L n / L T$). Apparently, this serves as a reasonable proxy for the contribution of a large number of modes. The question of the quasilinear modeling also motivates a comparison of the nonlinear modes with their linear counterparts.

It would be useful to establish what signatures, if any, of the linear eigenmodes appear in the nonlinear data. This is somewhat limited because linear eigenmodes are not easily distinguished in the nonlinear data. Non-orthogonality of eigenmodes means that there is not a one-to-one correspondence between any single eigenmode and a corresponding POD mode. Nevertheless, some similar mode structures do appear. Consider the POD modes for $ k y \rho s = 96$ shown in Fig. 13. Compared with eigenmodes seen in Fig. 3 for $ k y \rho s = 100$,^{52} the first two POD modes are similar to the second and first eigenmodes, respectively. However, these POD modes do not have the same frequencies as the eigenmodes. In fact, they have large positive frequencies, which can be seen in Fig. 14. In Sec. IV E, we show that the bulk nonlinear fluctuations have parallel scale lengths much smaller than the linear eigenmodes in most regions of k-space. It is only at the peak of the spectrum ( $ k y \rho s \u223c 96$) where the parallel scale lengths and frequencies of the POD modes are similar to the linear eigenmodes.

As noted above, the frequency spectrum exhibits a lobe of positive frequencies, while the linear instability has negative frequencies. Using the POD analysis, certain frequencies can be isolated to particular POD modes. Figure 14 shows the frequency spectrum $ | \sigma i u \u0302 i ( \omega ) |$, for selected **k**, plotted as a function of POD mode number *i* and *ω*, where $ u \u0302 i ( \omega )$ is the Fourier transform of the time POD mode $ u i ( t )$ defined above.

Near the peak of the heat flux spectrum at $ k y \rho s = 96$, the dominant frequencies alternate between positive and negative. The first POD mode has positive frequencies, while the second through fourth have mixed or negative frequencies. The positive frequency modes have no corresponding eigenmode that has been found. They appear to be nonlinear creations, excited through triplet interactions. To see how that might occur, it possible to link individual POD modes at different *k _{y}* through triplet interactions, with the conditions $ k \u2212 k \u2032 \u2212 k \u2033 = 0$ and $ \omega \u2212 \omega \u2032 \u2212 \omega \u2033 = 0$. As seen in Fig. 14, many POD modes have frequency spectra with well-defined peaks that can be identified as the dominant frequency of the mode. For example, consider the frequency of the first POD mode of $ ( k x , k y ) \rho s = ( 0 , 96 )$: $ \omega k , 1 \u2248 325$ (Fig. 14, center). This matches with $ ( k x \u2032 , k y \u2032 ) \rho s = ( 8.09 , 64 )$: $ \omega k \u2032 , 1 \u2248 350$ (Fig. 14, right) and $ ( k x \u2033 , k y \u2033 ) \rho s = ( \u2212 8.09 , 32 )$: $ \omega k \u2033 , 1 \u2248 \u2212 25$ (Fig. 14, left). Therefore, these modes all meet the condition to interact nonlinearly. This demonstrates how the large, positive frequency modes might receive energy despite the eigenvalue solver finding no corresponding unstable modes.

Another possibility is that these modes are in the pseudospectrum of the linear operator. That is, they may be eigenmodes of a perturbed linear operator, and the nonlinearity provides just such a perturbation.^{53}

It is worth discussing the POD modes described here in the context of the model of Chen, Zonca, and Lin.^{4} The POD modes encompass the entire spectrum of fluctuations, including the nonlinear structures that correspond to the linear eigenmodes. These POD modes (those corresponding to the instabilities) clearly are not connected to the quasimodes. However, our sense is that some subset of the POD modes is indeed similar to the quasimodes, i.e., those that are nonlinearly excited. One possible distinction is that many of the POD modes—even those dissimilar to the linear modes—are observed to contribute substantially to the heat flux. On the other hand, our expectation is that the quasimodes probably do not contribute to the heat flux.^{21} In short, we think it likely that some subset of the POD modes—particularly those in the low-*k _{y}* region of high energy and low heat-flux—are likely very similar to the quasimodes.

### E. Linear and nonlinear timescales

While saturation of the nonlinear state is the consequence of the balance of the linear injection with the nonlinear transfer at some scale, there is another question of how the rest of the scales relate to one another. Prior gyrokinetic studies have found evidence in similar systems for the critical-balance conjecture that relates the characteristic timescales over a range of spatial scales.^{54} Here, we provide empirical numerical evidence that this system satisfies the condition of comparable parallel streaming times and nonlinear perpendicular decorrelation times over a wide range of perpendicular spatial scales.

To compare the parallel streaming time and perpendicular decorrelation time, we compute the characteristic parallel wavenumber $ k \u2225 , rms$ and nonlinear frequency, respectively. The streaming time is then defined as $ \tau \u2225 \u2212 1 \u2261 k \u2225 , rms a v T e / c s$. These are computed by weighting with $\Phi $, as described in Appendix A, and the results are shown in Fig. 15. In this figure, values of $ \tau \u2225 \u2212 1$ and $ \omega rms$ computed at *k _{x}* = 0 for each

*k*are labeled “NL.” Also plotted are values for the first two POD modes (separately labeled) and the linear frequencies of the unstable eigenmodes.

_{y}It is observed that while both $ \tau \u2225 \u2212 1$ and $ \omega rms$ grow with *k _{y}*, their ratio (bottom of Fig. 15) remains fairly constant. In fact, it is very near unity, providing evidence for critically balanced fluctuations. This happens despite a system where the drive is dominant over nonlinear transfer for $ k y \rho s \u2272 200$, and where an inertial range—where both drive and dissipation terms are subdominant to nonlinear transfer—does not exist.

This observation of critically balanced fluctuations enables a broad interpretation of many of the observations of this study. It is interesting to observe (see Fig. 15) that the only wavenumber at which the time scales of the dominant nonlinear fluctuations match those of the linear eigenmodes is near $ k y \rho s \u2248 100$. Notably, this is also the wavenumber range at which the nonlinear spectrum peaks (recall Fig. 7). Based on these observations, we posit the following picture. The linear eigenmodes are predominant only in one narrow wavenumber range: $ k y \rho s \u2248 100$. In this wavenumber range, nonlinear interaction stimulates additional fluctuations at *k _{z}* larger than that of the linear eigenmodes (evidenced by the large values of $ \tau \u2225 \u2212 1$ in the top plot of Fig. 15). Critical balance then compels the nonlinear time scales to increase to match the higher

*k*fluctuations resulting in the superlinear time scales observed in the nonlinear frequency spectrum (evidenced by the large values of $ \omega rms$ in the middle plot of Fig. 15). At scales $ k y \rho s < 100$, the fluctuations are nonlinearly driven by interaction with the dominant linear wavenumber ( $ k y \rho s \u2248 100$). This is confirmed by the nonlinear transfer functions shown in Fig. 10, where there is a red lobe at the dominant wavenumber ( $ k \u2032 y \rho s \u2248 100$) for all of the smaller wavenumbers $ k y \rho s < 100$ shown in the figure. The fast (superlinear) time scales and large

_{z}*k*persist at these scales because they are predominantly nonlinearly driven and adapt in combination with one another to satisfy critical balance.

_{z}## V. CONCLUSION

This paper reports on an investigation into the fundamental nature of ETG turbulence in the pedestal as modeled with gyrokinetic simulations. We employ frequency spectra, free energy balance, nonlinear transfer functions, and POD to probe the system. The main results can be summarized as follows. Linear physics dominates only in a narrow range of k-space, $ k y \rho s \u2248 100$, close to the peak in the linear spectrum. In this region, the *k _{z}* and

*ω*of the linear eigenmodes closely match those of an important nonlinear mode extracted via POD. In all other regions of k-space, this does not pertain: the linear frequencies and parallel wavenumbers are far smaller than those of the nonlinear fluctuations. Even in the range where the linear physics persists, additional, nonlinear high-

*k*, high-

_{z}*ω*(superlinear) fluctuations make substantial contributions to the free energy and heat flux. We posit that a nonlinear cascade drives these high-

*k*fluctuations in this peak region of k-space. Subsequently, these fluctuations couple with all other regions of k-space to drive intrinsically nonlinear fluctuations that dominate over the linear instabilities. This process is mediated by critical balance, which is manifest across the entire domain despite enormous variation in linear drive, parallel fluctuation scales, and nonlinear time scales. This picture is supported by (1) frequency spectra, which exhibit frequencies far surpassing the linear frequencies at all wavenumbers, even exhibiting high amplitude at positive frequencies (the opposite sign of the ETG instability); (2) nonlinear transfer functions, which demonstrate an inverse cascade mediated by scales in the driving region $ k y \rho s \u2248 100$; and (3) the ratio $ k z v T e / \omega rms$ calculated from the nonlinear fluctuations, which is bounded between 1.5 and 2 in all regions of k-space (critical balance predicts order unity). Although this work substantially advances our understanding of this turbulent system, some questions remain, e.g., what precise nonlinear mechanism is responsible for exciting the small fluctuation scales in the drive range.

_{z}This system exhibits distinct characteristics in comparison with core turbulence phenomenology, which generally exhibits a single dominant linear eigenmode spanning the entire drive range (with some subdominant stable mode activity). In the conventional core scenario, this unstable eigenmode dominates the turbulent transport and clear signatures of its frequency and parallel length scales are observed in the fluctuations. Moreover, in this pedestal system, neither streamers nor zonal flows play a dominant role in the saturation process.

This work has implications for formulating reduced models of pedestal transport. Despite the fact that intrinsically nonlinear structures dominate the fluctuations at most scales, the persistence of linear physics in a single narrow region of k-space is consistent with the efficacy of a generalized quasilinear model, which uses the linear eigenmode at a single wavenumber to predict nonlinear heat flux.^{8}

## ACKNOWLEDGMENTS

The authors would like to thank Ben Chapman-Oplopoiou, Gabriele Merlo, Ian Abel, Adrian Fraser, and M.J. Pueschel for useful discussions.

This work was supported by the U.S. DOE Office of Fusion Energy Sciences Scientific Discovery through Advanced Computing (SciDAC) program under Award No. DE-SC0018429 and by the U.S. DOE under Contract No. DEFG02–04ER54742. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Justin Walker:** Conceptualization (equal); Data curation (lead); Investigation (lead); Methodology (equal); Software (equal); Supervision (equal); Visualization (lead); Writing – original draft (lead); Writing – review & editing (equal). **David R. Hatch:** Conceptualization (equal); Funding acquisition (lead); Methodology (equal); Project administration (lead); Resources (lead); Supervision (equal); Writing – original draft (supporting); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### APPENDIX A: DEFINITIONS

*J*(

*z*) is the Jacobian, and the geometry factors are related by $ B 2 = C 2 [ g x x g y y \u2212 ( g x y ) 2 ]$ and $ J ( z ) = det ( g i j ( z ) )$, where

*g*is the metric. Note that we have dropped the explicit

*k*,

_{x}*k*dependence for readability.

_{y}### APPENDIX B: PROPER ORTHOGONAL DECOMPOSITION

In this section, a short description of proper orthogonal decomposition (POD) and the details of our particular process are presented.

##### 1. Introduction

POD allows one to construct orthonormal bases for the data in time and space. The basis vectors are ordered such that the first *r* vectors are the optimal projection of the data onto a rank–*r* subspace. With a POD, one can check the contributions that each mode makes to various quantities, such as the heat flux *Q*. The POD is related to the singular value decomposition (SVD) and can be computed that way.

*M*×

*N*. However, this paper will use the continuous variables

*t*and

*z*to label equations. By performing an SVD on the function, one obtains

*U*,

*V*are unitary matrices and $ ( \xb7 ) \u2020$ denotes a conjugate transpose. The columns of $ U M \xd7 M$ are the left singular vectors of

*f*, the columns of $ V N \xd7 N$ are the right singular vectors of

*f*, and $ \Sigma M \xd7 N$ is a diagonal matrix with singular values $ { \sigma i | i \u2208 { 1 , \u2026 , min ( M , N )}}$ along the diagonal. Thus, one can think of the columns $ v i ( z )$ of

*V*as the “spatial modes” of

*f*and the columns $ u i ( t )$ of

*U*as the “temporal modes” of

*f*, when weighted by the singular value

*σ*. One can rewrite Eq. (B1) as

_{i}*f*as a sum of the outer product of the “temporal” and “spatial” basis vectors.

##### 2. Application to gyrokinetic system

*M*×

*N*dimensional matrix,

*ζ*would then be an $ M \xd7 4 N$ dimensional matrix. It must be noted that the

*z*used here includes all

*k*connected through the parallel boundary condition. See, e.g., Told.

_{x}^{7}

*z*with the modes from a POD of

*ζ*, the Jacobian prevents the modes from being orthogonal. Therefore, the modes must be re-scaled to be orthogonal with respect to the Jacobian norm. By choosing to re-scale the fields by $ J ( z )$, a new vector can be defined,

##### 3. $ \omega rms$ and $ k \u2225 , rms$ of POD modes

## REFERENCES

^{56}explicit tests demonstrate that they have a negligible impact in our simulations.

*j*will be dropped in the remainder of this article, since electrons are the only kinetic species.

*k*range of the toroidal ETG modes, as studied in Ref. 19, which is a multiscale problem beyond the scope of the work.

_{y}*k*,

_{x}*k*).

_{y}*k*, the same eigenmode branches are seen for $\n\nk\ny\n\n\rho \ns\n=\n80$, and can be inferred to exist at $\n\nk\ny\n\n\rho \ns\n=\n96$.

_{y}