Ultra Long Turbulent Eddies, Magnetic Topology, and the Triggering of Internal Transport Barriers in Tokamaks

Understanding and reducing the turbulence that transports energy out of the plasma is crucial to designing an economically competitive magnetic confinement fusion power plant [1]. In this Letter, we will present high-fidelity local (flux tube) gyrokinetic simulations of tokamak plasmas demonstrating that nonlinear turbulent eddies become much longer than previously appreciated at very low magnetic shear $\hat{s}\approx 0$ (i.e. weak variation of the field line pitch between magnetic surfaces). While eddies are typically assumed to have lengths comparable to the connection length [2, 3] (i.e. the distance along the magnetic field between the outboard and inboard midplanes of the torus), when $\hat{s}\approx 0$ we show that they can be hundreds of times longer. This can be seen as a nonlinear analogue of “giant electron tails” — extended parallel structures observed in linear modes that result from kinetic electron physics [4, 5].

When eddies become hundreds of poloidal turns long, they can wrap around the device many times and “bite their own tail.” As we will see, the resulting parallel self-interaction can dramatically alter the transport properties of the turbulence. Hence, to accurately simulate real experiments with low $\hat{s}$, it becomes necessary to faithfully model the correct field line topology. This has rarely been studied. Even the original paper formulating the flux tube [6] recommended always increasing the length of the simulation domain in order to entirely eliminate self-interaction. More recent works have investigated the impact of field line topology, but for standard values of shear $\hat{s}\sim 1$ [7, 8, 9, 10, 11, 12, 13, 14]. While some global gyrokinetic simulations can correctly model the topology, due to their high computational cost, there have been few studies that included both $\hat{s}\approx 0$ and kinetic electrons [15]. An important exception is [16], which investigated Internal Transport Barriers (ITBs). ITBs are an improved confinement regime characterized by a localized steepening of the pressure profile in the core [17, 18, 19]. Importantly, experimental observations have revealed that ITBs are enabled by low $\hat{s}$ and preferentially form near surfaces with rational values of the safety factor $q$ (i.e. the number of toroidal turns field lines make around the torus per poloidal turn) [20, 21]. Accordingly, [16] found large corrugations in temperature, density, and flow profiles associated with $\hat{s}\approx 0$ rational surfaces, but with limited insight into the physics at play.

In this Letter, we employ gyrokinetics, a high fidelity kinetic model of magnetized plasma turbulence [22, 23, 24, 25], and solve the coupled system of gyrokinetic Vlasov and Maxwell’s equations numerically using the spectral code GENE [26]. The results that follow were obtained with collisionless, electrostatic simulations and require a full kinetic treatment of both ions and electrons. We specify the equilibrium with Miller geometry [27] and use a local domain [6]. A local treatment takes advantage of the anisotropic nature of the fluctuations and constructs a long, slender domain along magnetic field lines where all equilibrium quantities are Taylor expanded around a magnetic surface of interest and only vary along the magnetic field.

We will study standard Ion Temperature Gradient (ITG)-driven turbulence at zero magnetic shear $\hat{s}=0$ using two cases: the standard Cyclone Base Case (CBC) [28] and a “pure ITG” (pITG) case. The latter is simply the former without the density and electron temperature gradients, so that the ion temperature gradient is the only source of free energy driving turbulence. We choose these two cases as we will find them to be dominated by different ITG instability branches. We will employ both linear and nonlinear local gyrokinetic simulations and vary two parameters, $N_{pol}$ and $\Delta y$, to understand the impact of topology. The parameter $N_{pol}$ is the length of the simulation domain in number of poloidal turns, which controls the distance that an eddy must stretch before it can return to its starting parallel location. The parameter $\Delta y$ is the shift in the binormal direction that all field lines experience as they pass through the parallel boundary of the domain (i.e. in the direction perpendicular to the field line, but within the magnetic surface). Mathematically, this shift enters the parallel boundary condition as [11]

where $\delta A$ represents any turbulent quantity, $L_{y}$ is the binormal domain width, $x$ and $y$ are the minor radial and the binormal coordinates respectively, and $z$ is the poloidal angle (i.e. the parallel coordinate). Statistically motivated periodic boundary conditions are used in the x and y directions, so applying binormal periodicity removes $N_{*}L_{y}$ (given that $N_{*}\in\mathbb{N}$). Thus, $\Delta y=0$ creates field lines that all exactly close on themselves after one time through the domain and, for example, setting $\Delta y=L_{y}/2$ will let field lines pass through the parallel boundary twice before closing. At finite shear $\Delta y$ corresponds to a simple radial shift of the lowest order rational surfaces, so without loss of generality $\Delta y$ can be set to zero, but this is no longer true when $\hat{s}=0$ [6]. Accordingly, we have implemented the shifted parallel boundary condition in GENE and benchmarked the changes against analytical ITG and parallel velocity gradient (PVG) instabilities in the cold ion limit [29, 30]. Including $\Delta y$ is essential to properly model physically possible magnetic topologies.

While appropriately choosing $N_{pol}$ and $\Delta y$ given an arbitrary physical equilibrium is quite complicated [30], many cases of practical importance are straightforward. For example, a field line on a $q=3$ surface exactly closes on itself after one poloidal turn, so one should take $N_{pol}=1$ and $\Delta y=0$. Similarly, a $q=5/2$ field line closes after two turns, meaning one can take $N_{pol}=2$ and $\Delta y=0$. Finite $\Delta y$ is needed for surfaces that are very close to, but not exactly rational. For example, since eddies are $\sim 10\rho_{i}$ wide in the binormal direction, if a field line comes within four ion gyroradii $\rho_{i}$ of closing on itself after one poloidal turn, we should set $N_{pol}=1$ and $\Delta y=4\rho_{i}$. This ensures that we capture the complexity of an eddy partially (but not exactly) biting its own tail. Given the scale separation of gyrokinetics, a field line that comes within a few ion gyroradii of being integer would correspond to a physical value of $q$ that is asymptotically close to, but not actually a rational, e.g. $q=3+O(\rho_{i}/a)$ where $a$ is the device minor radius.

In linear simulations with $N_{pol}=1$ and $\Delta y=0$ we observe unique behavior as the magnetic shear drops below $|\hat{s}|\lesssim 0.1$. The ballooning structure of eigenmodes become broader and toroidal ITG (tITG) weakens relative to slab ITG (sITG), causing a transition in the dominant linear mode in both the CBC and pITG cases. We identify the tITG $\rightarrow$ sITG mode transition based on qualitative agreement with previous analytical results [31, 32], a large discontinuous increase in the real frequency, and a sharp broadening of the eigenmode in ballooning space. This broadening at low magnetic shear suggests self-interaction is becoming increasingly important.

Accordingly, Fig. 1 shows a scan in $N_{pol}$ at $\Delta y=0$ and $\hat{s}=0$. Physically, this corresponds to studying how linear stability is affected by changing the order of the rational surfaces (e.g. integer, half-integer, third-integer). We see that the sITG mode is independent of the order of the rational surface, whereas the tITG mode is stabilized by low order rational surfaces. In a long domain, sITG modes are dominant for pITG gradients and tITG modes for CBC gradients. However, as the domain is shortened the tITG instability is stabilized and transitions to sITG at $N_{pol}=1$ for CBC as well.

Figure 2 shows how $\Delta y$ changes the linear drive when $\hat{s}=0$ and $N_{pol}=1$, indicating the impact of magnetic field lines not exactly closing on themselves. We see that small changes of order $\Delta y\approx\rho_{i}$ can cause the dominant instability to transition between sITG and tITG and change the growth rate by a factor of 2 (finite $\Delta y$ increases the growth rate for CBC and decreases it for pITG). Though not shown here, in some cases changing $\Delta y$ can even control if the mode is linearly unstable or not. We also see that the growth rate for sITG peaks around $\Delta y=0$, i.e. when magnetic field lines exactly close, whereas the tITG growth rate is maximum at $\Delta y\approx\rho_{i}$. Thus, we can expect that sITG will tend to be strongest exactly on rational surfaces and tITG away from them.

Finally, using Floquet-Bloch theory [33, 34], the linear growth rates from a $\Delta y$ scan with $N_{pol}=1$ (e.g. Fig. 2) can be used to analytically calculate growth rates for any value of $N_{pol}>1$ and $\Delta y$ [30]. This is because each poloidal turn is identical due to axisymmetry. For instance, a simulation with $N_{pol}=4$ and a given $\Delta y|_{N_{pol}=4}$ is equivalent to 4 identical linked $N_{pol}=1$ simulations, each with a binormal shift $\Delta y|_{N_{pol}=1}$ that, when repeated 4 times, corresponds to $\Delta y|_{N_{pol}=4}$. Thus, the allowed $\Delta y|_{N_{pol}=1}$ values are $\Delta y|_{N_{pol}=1}=(\Delta y|_{N_{pol}=4}+pL_{y})/4$, where $p\in\mathbb{N}$ is a free parameter. Each of these allowed $\Delta y|_{N_{pol}=1}$ has a different growth rate, the fastest growing of which will dominate and set the growth rate of the $N_{pol}=4$ simulation. Generalizing this to domains longer than $N_{pol}=4$, we see that as $N_{pol}\rightarrow\infty$ the binormal shift of the domain ceases to matter and the free parameter $p$ always enables the simulation to find the maximum growth rate from the $\Delta y$ scan with $N_{pol}=1$.

Analogously to the linear study of Fig. 1, we carried out a nonlinear scan in $N_{pol}$ with $\hat{s}=0$. Figure 3 shows how the total electrostatic heat flux $Q_{es}$ and $C_{||}(z_{1},z_{2})$ vary with the parallel domain length, where $C_{||}(z_{1},z_{2})$ is the two-point parallel correlation function [11] of the non-zonal electrostatic potential averaged over time, $x$, and $y$. We see that the heat flux requires $N_{pol}>100$ for convergence and parallel locations more than 50 poloidal turns apart remain highly correlated — direct evidence of ultra long turbulent eddies. In fact, we believe their length is only limited by critical balance — the distance that electrons can travel along the field line within the lifetime of a turbulent eddy [2]. Accordingly, simulations confirm that if the electron thermal velocity is artificially reduced (by increasing their mass), the parallel correlation length shrinks proportionally. In contrast, simulations with an adiabatic approximation for electrons exhibit parallel correlation spanning no more than a couple poloidal turns, reflecting the distance an ion can travel in an eddy’s lifetime. Thus, fully resolving kinetic electron effects (in particular with the physical $m_{e}/m_{i}$ mass ratio) is crucial to properly model low magnetic shear surfaces.

At small $N_{pol}$ we expect the heat flux in Fig. 3 to follow the linear trends of Fig. 1. Indeed, the heat flux has little variation for the pITG case as expected. However, for CBC, while the small jump in the heat flux at $N_{pol}=3$ could be consistent with the increase in the linear growth rate, a larger absolute change was expected. Hence, another $N_{pol}$ scan was performed with $\Delta y\approx 0.6\rho_{i}$ for the pITG case, which shows a sharp increase with $N_{pol}$ as is expected from linear results. At large $N_{pol}>40$, both parameter sets exhibit a similar decrease in the heat flux and parallel correlation, suggesting a common mechanism. This appears to be a result of the parallel domain being long enough to accommodate multiple eddies, instead of just one. When there is just a single correlated eddy, it can organize internally to maximize transport. However, when there is more than one eddy along the parallel direction they can organize within themselves, but there will always be transition regions where the different eddies meet and the radial transport is reduced.

Additionally, in the CBC simulations we observed long parallel waves that appeared when the domain exceeded $N_{pol}\approx 20$. As with the eddy length, their wavelength is proportional to the thermal electron velocity. These waves cause the oscillations in the parallel correlation profile shown in Fig. 3 (bottom). Their appearance is associated with the $\sim 30\%$ increase in the particle and heat transport as shown in Fig. 3 (top). Thus, they seem important, but a detailed study of them is left for future work.

Lastly, we scanned the binormal shift $\Delta y$ for domains with $N_{pol}=1$ and $\hat{s}=0$, which is equivalent to considering near integer $q$. We observe dramatic, fine-scale variation in the heat flux with $\Delta y$ as seen in Fig. 4. While this is surprising, independent evidence of this can be found in [35]. We see that low but finite $\Delta y\leq\rho_{i}$ (circle markers) can double the heat flux or completely stabilize turbulence, which can be explained by the variation of the linear drive (thin red line). At somewhat larger values of $\Delta y$ (triangle markers), the heat flux time traces become bursty, alternating between two distinct values of heat flux. During the peak of the bursts for pITG case (top panel), turbulence shows strong parallel correlation across the domain $C_{||}(z_{1}=0,z_{2}=\pi)\approx 0.8$, indicating that the turbulence has exactly bitten its tail. During the troughs of the bursts, the parallel correlation $C_{||}(z_{1}=0,z_{2})\approx 0.4$ is much lower, suggesting that the eddies don’t exactly bite their own tails. This trend flips for CBC simulations (as expected from Fig. 2)— turbulence shows weak parallel correlation during the peak of the burst and high correlation in the troughs. Thus, at these values of $\Delta y$ (triangle markers) the eddies stochastically connect/disconnect from themselves and alternate between the two states due to chaotic nature of the system. At larger values of $\Delta y$ corresponding to high-order rational surfaces (square markers), the eddies pass through the domain several times and then exactly bite their own tails. Since ultra long eddies can span hundreds of poloidal turns, on some high-order rational surfaces a single eddy can cover the full flux surface, even in a real device. To study this, we can set $L_{y}$ to correspond to the full flux surface of a given device. For example, on the eighth-order rational surface shown in Fig. 4 (bottom), we see that a single eddy tiles the entire surface. More important, this eight-fold periodicity does not give the eddy enough space to be its preferred size to maximise transport. Instead the lack of space in the poloidal direction “squeezes” the eddy, restricting its size through perpendicular self-interaction and reduces its ability to transport energy. Such poloidal eddy squeezing presents a novel strategy for improving confinement. This is accomplished by tailoring the safety factor profile to have particular high-order rational values where $\hat{s}=0$. It is attractive because, unlike integer surfaces, high-order rational surfaces do not enable magnetohydrodynamic (MHD) instabilities [36]. Lastly, at particular high-order rational values of $\Delta y$ (diamond markers) it is also possible to flip this effect and somewhat “stretch” eddies, thereby slightly increasing turbulent transport.

We have observed turbulent eddies that span hundreds of poloidal turns in gyrokinetic simulations of tokamaks with low magnetic shear $\hat{s}\approx 0$. Such ultra long eddies mean that the topology of the magnetic field becomes very important. A small binormal shift of $\Delta y\approx\rho_{i}$ at the parallel boundary (corresponding to a small change in the safety factor away from being rational) was found to either completely stabilise turbulence or increase the heat flux by a factor of two. This demonstrates the importance of faithfully constructing the topology of the device being modeled. Moreover, it may explain the experimental observation that ITBs are easier to trigger in the vicinity of rational surfaces. Additionally, we have seen that an individual ultra long eddy may have the ability to entirely cover a flux surface, even in large devices. This enables a novel strategy to improve confinement — by carefully selecting a high-order rational value of the safety factor, the eddy can box itself in and squeeze itself to a smaller perpendicular size. Several important considerations like plasma shaping (which can cause strong local magnetic shear), collisions (which may shorten the parallel correlation length), and electromagnetic effects (which allow the turbulence itself to modify the magnetic field geometry) are left for future work. Additionally, stellarators (which commonly have very weak global magnetic shear) have exhibited improved confinement around rational magnetic surfaces [37], but are also left for future study.