Key impact of phase dynamics and diamagnetic drive on Reynolds stress in magnetic fusion plasmas

The expression of the poloidal Reynolds stress emerges from the charge conservation equation, which governs the dynamics of the generalized vorticity. Its complete derivation and expression in the gyrokinetic framework, including all radial currents responsible for its time evolution, can be found in [21, 22], and more recently in [23]. In the present work, we only focus on the component governed by the advection by the E × B drift—then putting all other contributions into the right-hand side (RHS) term. Then, the dimensionless gyrokinetic equation reads as follows:

$$\begin{align} \partial_t \bar{f} + \{J[\phi],\, \bar{f}\,\,\} = \mathrm{RHS} \end{align} \tag{ 1 }$$

with the Poisson brackets defined by ⁴ $\{g,\, f \,\,\} = \frac{1}{r} \left[ \partial_r(g \partial_\theta f\,\,) - \partial_\theta(g \partial_r f\,\,) \right]$. Here, $\bar{f}$ is the distribution function of the gyro-centers, φ is the electric potential and J[φ] stands for the gyro-average operator applied to φ (strictly speaking, the gyro-average operator is $\langle J[\phi] \rangle_v$, with $\langle X \rangle_v = \int \textrm{d}^3v X = \int_0^{+\infty}2\pi\textrm{d}\mu \int_{-\infty}^{+\infty} \textrm{d} v_\parallel X$ the integral over the velocity space). For simplicity, the inhomogeneity of the magnetic field is not retained in the first stage, so that grad-B effects leading to turbulence equipartition are not accounted for. B is then simply replaced by B₀, which is absorbed in normalizations. The latter are the Larmor radius $\rho_0 = m_pv_{T0}/eB_0$ for length scales, the cyclotron period $\omega_{c0}^{-1} = m_p/eB_0$ for time scale and thermal velocity $v_{T0} = (T_0/m_p)^{1/2}$ for velocities, where m_p is the proton mass, e is the elementary charge and T₀ is a reference temperature. The electric potential is normalized by T₀/e. In the long-wavelength limit, J can be approximated by:

$$\begin{align} J \approx 1 + \frac{\mu}{2}\nabla_\perp^2 \end{align} \tag{ 2 }$$

where µ is the adiabatic invariant, namely the magnetic moment and $\nabla_\perp^2 = \frac{1}{r}\partial_r(r\partial_r \;) + \frac{1}{r^2}\partial_\theta^2$. Let us apply the gyro-average operator J to equation (1) and integrate it over the velocity space $\langle J[$(1)$] \rangle_v$. At this stage, one can use the quasi-neutrality constraint which can be formulated as follows in the Boussinesq approximation:

$$\begin{align} n_0 \Omega = \langle {J[\bar{\,f\,}]} \rangle_v \approx \langle \bar{\,f\,} \rangle_v + \frac{1}{2} \nabla_\perp^2 p_\perp \end{align} \tag{ 3 }$$

where $\Omega = \phi - \langle \phi \rangle - \nabla_\perp^2 \phi$ is the potential vorticity, with $\langle \phi \rangle \approx \int \!\! \int \phi\, \textrm{d}\theta \textrm{d}\varphi/4\pi^2$ the flux surface averaged potential. The first term on the RHS is the gyro-center density, and the second one accounts for finite Larmor radius effects in the long-wavelength limit, with $p_\perp = \langle \mu\bar{f}\,\, \rangle_v$ the transverse pressure. One then obtains up to higher-order terms in $(k_\perp\rho_i)^2$:

$$\begin{align} \partial_t \Omega & - \left\{\phi, \nabla_\perp^2 \phi + \frac{1}{2} \nabla_\perp^2 p_\perp \right\} - \frac{1}{2} \left\{p_\perp, \nabla_\perp^2 \phi \right\}\nonumber\\ & - \frac{1}{2} \nabla_\perp^2 \left\{p_\perp, \phi \right\} = \langle J[\mathrm{RHS}] \rangle_v. \end{align} \tag{ 4 }$$

Noticing that $\nabla_\perp^2 \{\phi, p_\perp \} = \{\nabla_\perp^2 \phi, p_\perp \} + \{\phi, \nabla_\perp^2 p_\perp \} + 2 \{\nabla_{\perp,i} \phi, \nabla_{\perp,i} p_\perp \}$, with $\nabla_{\perp,i}$ the derivative along the ith Cartesian coordinate in the transverse plane and taking the sum of doubled indexes, equation (4) can be recast as follows:

$$\begin{align} \partial_t \Omega - \left\{\phi + p_\perp, \nabla_\perp^2 \phi \right\} - \{\nabla_{\perp,i} p_\perp, \nabla_{\perp,i} \phi \} = \langle J[\mathrm{RHS}] \rangle_v. \end{align} \tag{ 5 }$$

Equivalently, equation (5) also admits a conservative form:

$$\begin{align} \partial_t \Omega - \nabla_{\perp,i} \left\{\phi + p_\perp, \nabla_{\perp,i} \phi \right\} = \langle {J[\mathrm{RHS}]} \rangle_v. \end{align} \tag{ 6 }$$

Here, the Poisson bracket corresponds to the flux of vorticity. Further taking the flux surface average and integrating over the radial direction provides the evolution equation of the E × B poloidal flow $\langle u_{E\theta} \rangle = \partial_r \langle \phi \rangle$:

$$\begin{align} \frac{\partial \langle u_{E\theta} \rangle}{\partial t} + \frac{1}{r^2} \frac{\partial}{\partial r} \left[ r^2 (\pi+\pi^\ast) \right] = \mathrm{RHS} \end{align} \tag{ 7 }$$

where the flux-surface averaged Reynolds stress (hereafter simply called Reynolds stress in short) and its diamagnetic component are given by:

$$\begin{align} \pi & = \langle \tilde{u}_{Er} \, \tilde{u}_{E\theta} \rangle \end{align} \tag{ 8 }$$

$$\begin{align} \pi^\ast & = \langle \tilde{u}_r^\ast \, \tilde{u}_{E\theta} \rangle. \end{align} \tag{ 9 }$$

The tilde refers to non-axisymmetric components, while $\tilde{u}_{Er} = -\frac{1}{r} \partial_\theta \tilde{\phi}$ and $\tilde{u}_r^\ast = -\frac{1}{r} \partial_\theta \tilde{p}_\perp$ denote the radial components of the electric and diamagnetic drifts, respectively. As expected, the flux of vorticity is related to the Reynolds stress, which is a restatement of Taylor's theorem [24]. Also, it appears that finite Larmor radius effects, carried by the J operator, lead to an additional contribution to the Reynolds stress, namely $\pi^\ast$, as already acknowledged in early publications [15, 16]. Appendix A shows the consistency between this calculation within the gyrokinetic framework and the fluid moment approach in the adiabatic limit. In particular, the diamagnetic origin of $\pi^\ast$ is assessed.

Hereafter, π and $\pi^\ast$ are named Reynolds stress and diamagnetic Reynolds stress, respectively.

The structure and dynamics of the total Reynolds stress are explored in a flux-driven gyrokinetic simulation of the GYSELA code [18]. Ion temperature gradient (ITG)-driven turbulence is modeled with an adiabatic electron response.

The simulation is similar to the one already reported in [25] (it is actually a re-run so as to get three-dimensional (3D) data which are required in the present analysis). Its radial domain covers $0\leqslant r/a\leqslant 1.3$, with a toroidal limiter treated as an immersed boundary, where relaxation toward a low-temperature Maxwellian is imposed by means of an additional Krook operator [26]. In the simulated scrape-off layer $r/a\geqslant1$, the electron density relative fluctuation is clamped to the floating potential φ = 3T_e /e. The $\rho^\ast = \rho_i/a$ parameter is equal to $\rho^\ast \approx 1/151$ and $\nu^\ast \approx 7.5\times10^{-2}$ at r/a = 0.5.

We have also performed the same analysis on another well-resolved simulation for which 3D output data were available. Its radial domain is $0\leqslant r/a\leqslant 1$, the distribution function driven toward a Maxwellian at the outer edge, characterized by e-folding density and temperature profiles which mimic those expected in the scrape-off layer of tokamak plasmas. Staircases are observed in the fully developed turbulent regime. In this case $\rho^\ast \approx 1/244$ at r/a = 0.5 where the collisionality is $\nu^\ast \approx 0.23$. Results regarding the present issue are quite similar.

Figure 1 portrays the 2-dimensional (2D) behavior, in space and time, of the flux surface averaged E × B poloidal flow $\langle u_{E\theta} \rangle$ and of the total Reynolds stress $\pi+\pi^\ast$. Two main observations can be made. First, the dynamics essentially features avalanche-like events, which draw an array made of diagonal structures. These avalanches propagate both inwards and outwards, roughly at the same speed. Second, there appears to be some correlation between the two images, attesting to the role of the total Reynolds stress in driving zonal flows. Deriving a vorticity equation which is tractable numerically is difficult in the gyrokinetic framework (see, for instance, the expressions proposed in [21–23]). Hereafter, we focus on the total Reynolds stress, bearing in mind its likely dominant contribution to the generation of dynamical corrugated features of the zonal radial electric field.

Zoom In Zoom Out Reset image size

The respective contributions of the electric and diamagnetic flows to the Reynolds stress π and $\pi^\ast,$ respectively are displayed in figure 2. Two main observations can be made. First, π and $\pi^\ast$ are approximately in phase, with minima and maxima occurring at roughly the same time and radial locations. Second, π* is about twice the value of π. This is shown in figure 3, where each point corresponds to the values of π and $\pi^\ast$ at any given radial position $0.3\leqslant r/a \leqslant0.9$ (the radial domain excludes the central region where a non-vanishing heat source is applied and the outer buffer region) and time $335\,700\leqslant \omega_c t \leqslant 387\,900$ (which corresponds to $347.88 \lesssim t\, v_{Ti}/R_0 \lesssim 402.02$, so well after the overshoot, into the nonlinear saturation phase). The important conclusion is that FLR corrections to the Reynolds stress tensor cannot be ignored, and may even prove dominant in magnitude.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

At this point, it is worth noticing that this result neither implies nor means that $(k_\perp\rho_i)^2$ terms are of order unity; they are actually not, as attested to by the strong cut-off of the turbulence spectra above $k_\perp\rho_i\sim1$ (see e.g. [27]). The phase shift plays a critical role there. Accounting for higher-order terms in $(k_\perp\rho_i)^{2m}$ ($m\geqslant2$) would involve cross-correlations between fluctuations of the electric drift velocity and of high-order derivatives of high-order fluid moments.

The fact that $\pi^\ast$ is found to be larger than π is in marked difference with previous findings. To our knowledge, the only other observation in a gyrokinetic simulation was done by Dimits et al [17]. They report a negligible role of the diamagnetic Reynolds stress in an ITG turbulence simulation for the cyclone base case. Apart from the possible issue with the definition of $\pi^\ast$, there are several differences between this simulation and the ones reported in our paper, which may result in these different findings. First of all, the simulation by Dimits is gradient-driven, so that the equilibrium pressure profile is kept constant over time, conversely to our flux-driven simulations. It remains unknown how this constraint may impact the fluctuation dynamics and magnitude of the transverse pressure fluctuations and their gradient, which governs the diamagnetic Reynolds stress $\pi^\ast$. Second, cyclone base case gradient-driven simulations are close to the so-called Dimits upshift [28], i.e. close to the nonlinear threshold for turbulent transport. As reported e.g. in [29], they tend to be characterized by a secular growth of zonal flows, rendering early conclusions difficult to interpret. In this respect, and finally, Dimits' result is obtained in the early phase of the simulation: the time window which we have considered for the analysis starts at a time $t = 33.57\times10^4\,\omega_c^{-1} \approx 2.10^3 \, L_T/v_{Ti}$ roughly equal to twice the final time of the simulation considered by Dimits. In the well-saturated turbulent regime that we consider, Reynolds stress temporal fluctuations are much larger than their almost vanishing mean, conversely to what happens just after the turbulence overshoot in figure 2(a) of [17]. In this same figure, the final state actually looks closer to our finding, especially in the sense that both Reynolds stress components tend to oscillate around a much smaller mean than their standard deviation, and that these standard deviations seem to reach comparable magnitudes, possibly in phase (although this latter point is more difficult to claim from the time traces). Incidentally, the analysis of the second simulation where we have 3D data from the linear phase reveals that the diamagnetic Reynolds stress is always twice as large as the Reynolds stress, even in the early turbulence regime from the overshoot onwards.

The approximate proportionality $\pi^\ast \approx 2\pi$ also indicates that the transverse pressure is not simply advected by the $\mathbf{E} \times \mathbf{B}$ flow. Indeed, assuming $\partial_t p_\perp + \mathbf{u}_E\cdot\boldsymbol{\nabla} p_\perp = 0$, the linear analysis would lead to the following relation in the 3D Fourier space: $\hat p_{\perp,mn\omega} = -(\omega^\ast_p/\omega)\hat \phi_{mn\omega}$, with $\omega^\ast_p = k_\theta \partial_r \langle p_\perp \rangle$ the diamagnetic frequency. Here, m and n are the poloidal and toroidal wave numbers, $k_\theta = m/r$ is the poloidal wave vector and $\langle p_\perp \rangle$ is the time and flux surface average of $p_\perp$. In this case, $\pi^\ast = - \sum_{\omega,k} k_rk_\theta (\omega^\ast_p/\omega) |\hat\phi_{mn\omega}|^2$ is expected to be out of phase with respect to $\pi = \sum_{\omega,k} k_rk_\theta |\hat\phi_{mn\omega}|^2$. This result could have important implications when interpreting experimental data (see e.g. [30]).

In order to look for the origin of the close to zero-phase shift between π and $\pi^\ast$, the properties of the Fourier components of pressure and electric potential are investigated; see figure 4. Only resonant modes are selected, i.e., such that n + m/q(r)≈ 0, since they are the main contribution to the Reynolds stress. The top panel shows that the amplitudes of the modes peak approximately at the diamagnetic frequency $\omega^\ast_p$ (dot–dash black line). The bottom graph is a cut at m =−45. One recovers the ratio of about two between $|\hat p_{\perp mn\omega}|$ and $|\hat\phi_{mn\omega}|$. The relationship between $\hat p_{\perp mn\omega}$ and $\hat \phi_{mn\omega}$ is investigated on a systematic basis in figure 5 which displays all data points of $\Re(\hat p_{\perp mn\omega})$ as a function of $-(\omega^\ast_p/\omega)\Re(\hat\phi_{mn\omega})$, $\Re$ denoting the real part. These are 3D arrays, functions of ω, m and r. The range of m values has been restricted to $-200\leqslant m\leqslant0$. The cloud of points does not exhibit any clear structure. Consequently, there is no particular trend indicating some form of correlation. If anything, the sign between both quantities rather appears to be negative. In agreement with the result that π and $\pi^\ast$ are in phase, one can notice that $\Re(\hat p_{\perp mn\omega})$ and $\Re(\hat\phi_{mn\omega})$ have the same sign at large amplitude. Since the points are not aligned on the diagonal (black dashed line), this means that the fluctuations of $p_\perp$ are not simply advected by the $\mathbf{E}\times \mathbf{B}$ flow.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

From the theoretical point of view, several mechanisms can explain the departure from the passive scalar hypothesis. Within the linear framework and from the fluid perspective, it was already noticed [15] in equation (27) that compressibility could break the simple relationship $\hat p_{\perp,mn\omega} = -(\omega^\ast_p/\omega)\hat \phi_{mn\omega}$. Also, considering the drift-kinetic (i.e., neglecting finite Larmor radius corrections) linear response, one obtains for the Fourier modes $\mathbf{k}, \omega$ ($\mathbf{k}$ standing for m, n and k_r , the radial wave vector):

$$\begin{align} \hat p_{\perp \mathbf{k}, \omega} = - \left\langle \frac{\omega^\ast - \omega_d - k_\parallel v_\parallel }{\omega - \omega_d - k_\parallel v_\parallel} \; \mu B\; F_M \right\rangle_v \hat \phi_{\mathbf{k},\omega} \end{align} \tag{ 10 }$$

with $\omega_d \approx \mathbf{k}.\mathbf{v}_d$ the frequency-operator associated to the vertical drift $\mathbf{v}_d \approx (mv_\parallel^2+\mu B) (\mathbf{b}\times \boldsymbol{\nabla} B)/eB^2$ and F_M the equilibrium Maxwellian distribution function. The brackets $\langle ... \rangle_v$ denote velocity-space integration. It readily appears that the transverse pressure and electric potential deviate from the above-mentioned simple relationship due to both parallel dynamics $k_\parallel v_\parallel$ and compressibility ω_d . The simulation results reported here show that the usual assumption that their effect is negligible does not hold.

Looking for possible leading terms in the above linear relationship, and/or whether it remains valid in the turbulent regime is certainly an important issue worth addressing in future works.

More insight can be gained into the structure of the Reynolds stress equations (8) and (9) by looking at their expression in the Fourier space. The electric potential φ and transverse pressure $p_\perp$ are decomposed in Fourier modes along the periodic directions, the poloidal θ and toroidal ζ angles (no Fourier transform in time is performed there):

$$\begin{align*} \phi & = \sum_{m,n} \hat\phi_{mn} \textrm{e}^{i(m\theta+n\zeta)} = \sum_{m,n} |\hat\phi_{mn}| \textrm{e}^{i\varphi^\phi_{mn}} \; \textrm{e}^{i(m\theta+n\zeta)}\\ p_\perp & = \sum_{mn} \hat p_{\perp,mn} \textrm{e}^{i(m\theta+n\zeta)} = \sum_{m,n} |\hat p_{\perp,mn}| \textrm{e}^{i\varphi^p_{mn}} \; \textrm{e}^{i(m\theta+n\zeta)} \end{align*}$$

with $\varphi^\phi_{mn}(r,t)$ and $\varphi^p_{mn}(r,t)$ the phase of the modes. Then, π and $\pi^\ast$ can be recast as follows:

$$\begin{align} \pi & = -\frac{2}{r} \sum_n\sum_{m\gt0} m\; \Im( \hat\phi_{mn}^* \partial_r \hat\phi_{mn} ) \nonumber\\ & = -\frac{2}{r} \sum_n\sum_{m\gt0} m |\hat\phi_{mn}|^2\; \partial_r \varphi^\phi_{mn} \end{align} \tag{ 11 }$$

and

$$\begin{align} \pi^\ast& = -\frac{2}{r} \sum_n\sum_{m\gt0} m\; \Im( \hat p_{\perp,mn}^* \partial_r \hat\phi_{mn} ) \nonumber \\ & = -\frac{2}{r} \sum_n\sum_{m\gt0} m |\hat\phi_{mn}| |\hat p_{\perp,mn}|\\ &\quad \times \left[ \cos\Delta\varphi_{mn} \partial_r \varphi^\phi_{mn} + \sin\Delta\varphi_{mn} \partial_r (\ln |\hat\phi_{mn}|) \right] \nonumber \end{align} \tag{ 12 }$$

with $\Delta\varphi_{mn} = \varphi^\phi_{mn} - \varphi^p_{mn}$ the phase shift between the electric potential and transverse pressure fluctuations. $\Im$ denotes the imaginary part, and 'g*' is the complex conjugate of g. Expression (11) highlights that both phase gradient and mode amplitude are key to having a non-vanishing Reynolds stress π, as already emphasized in early works [14]. Conversely, $\pi^\ast$ can be generated by the gradient of either the phase or the amplitude depending on the phase shift Δϕ_mn between the two fields, which then governs the respective weight of each contribution.

To highlight the role of the phase in the dynamics of the Reynolds stress, we first consider here the reduced nonlinear model proposed in [31] and recalled in appendix B. It derives from the interchange instability, further considering the flute modes $k_\parallel = 0$ and a single poloidal wave vector denoted k_y . It extends the previous 1D model that features avalanche-like transport events [32] by keeping track of the phase of the density (or pressure since temperature is assumed constant) and electric potential fluctuations. Note that the model is derived in the limit of cold ions, so that the diamagnetic Reynolds stress is not accounted for. On the basis of the results discussed in the previous section, the magnitude of the total Reynolds stress $\pi+\pi^\ast$—and possibly of the zonal flows—is expected to be enhanced when relaxing this assumption T_i ≈ 0. For the time being, given the considered limit, the poloidal momentum equation reduces to $\partial_t V_\mathrm{eq} = - \partial_x \pi$, with $V_\mathrm{eq} = \langle u_{E\theta} \rangle$ the poloidal equilibrium flow and x standing for the radial coordinate. The Reynolds stress π takes the simple following form:

$$\begin{align} \pi = - \Im( 2k_y\, \phi_k^* \, \partial_x\phi_k ) \end{align} \tag{ 13 }$$

with φ_k given equation (B.4). As expected from Taylor's identity, the Reynolds force is equal to the vorticity flux: $\partial_x \pi = - \Im(2k_y\ \phi_k^*\, \partial_x^2\phi_k)$. An alternative form of π can be obtained by introducing the magnitude and phase of the Fourier mode of the electric potential fluctuations: $\phi_k(x,t) = |\phi_k(x,t)|\, \exp[i\varphi_\phi(x,t)]$, with $\varphi_\phi \in \mathbb{R}$ the phase of the φ_k field. The Reynolds stress then reads

$$\begin{align} \pi = -2k_y\ |\phi_k|^2\, \partial_x \varphi_\phi. \end{align} \tag{ 14 }$$

The time evolution of the complex fluctuating fields n_k (or p_k since temperature is assumed constant) and $\phi_k \in \mathbb{C}$ can be re-expressed in terms of phase and amplitude dynamics as detailed in appendix B.2.

Here, we compare two simulations which only differ by the magnitude S₀ of the source term. The values of the various parameters of equation (B.6) have been chosen to be close to the ones used in [31]: N_x  = 256, L_x  = 1, Δt = 2.10⁻³, k_y  = 2π, g = 1 and $D_0 = D_1 = \nu_0 = \nu_1 = 6.5\times10^{-2}$. The source term has a Gaussian shape:

$$\begin{align} S = S_0\exp\left\{- \frac{(x-0.1)^2}{0.01}\right\} \end{align} \tag{ 15 }$$

with S₀∈{20, 30}. The boundary conditions are Dirichlet at x = {0, L} for all fields, except for $n_\mathrm{eq}$ at x = 0 where we use Neumann: $\partial_x n_\mathrm{eq}(x = 0,t) = 0$.

The time dynamics between the two cases S₀ = 20 and S₀ = 30 are different, as evidenced on figure 6. For the small source magnitude, fluctuations first grow exponentially in time during the linear phase, and then saturate at a steady-state value. The system reaches a low confinement regime with large amplitude fluctuations and no equilibrium flow. Conversely, for S₀ = 30, the system bifurcates toward an improved confinement regime at t ≈ 75. It appears that this transition is due to the generation of an equilibrium flow $V_\mathrm{eq}$.

Zoom In Zoom Out Reset image size

Most interestingly, it turns out that the build-up of $V_\mathrm{eq}$ results from an instability of the gradient of the phase of the electric potential $\partial_x \varphi_\phi$. Indeed, as shown on figure 7, the phase gradient starts growing exponentially in time when the magnitude of the fluctuations reaches the first saturated regime. In fact, this governs the exponential growth of the Reynolds stress π up to a second saturation regime, at lower amplitude $|\phi_k|$. Conversely, $\partial_x \varphi_\phi$ does not develop any instability and remains vanishing for S₀ = 20, resulting in a vanishing Reynolds stress.

Zoom In Zoom Out Reset image size

Finding the expression of the growth rate of the phase instability is, however, not obvious. The expression of the time dynamics of the phase gradient $\partial_x \varphi_\phi$ is cumbersome, involving the equilibrium flow $V_\mathrm{eq}$ and the cross-phase between the various fluctuating fields. Yet, one can gain some insight about the physics at work by first noticing that the gradient of the vorticity phase $\partial_x \varphi_w$ evolves according to $k_y\partial_x V_\mathrm{eq}$ among other terms, equation (B.18): $\partial_t (\partial_x \varphi_w) = -k_y\partial_x V_\mathrm{eq} + \ldots$. This relation is reminiscent of the main expected effect of the velocity shear on turbulence eddies, namely to increase their radial wave vector k_x : $\textrm{d} k_x/\textrm{d} t = -k_y\ \partial_x V_\mathrm{eq}$, with $k_x \approx \partial_x \varphi_\phi$ [4]. During the exponential growth, it appears that all phase gradients exhibit the same growth rate, so that one can also infer $\partial_t (\partial_x \varphi_\phi) = -k_y\partial_x V_\mathrm{eq} + \ldots$. Then, accounting for the momentum balance equation with the Reynolds stress, a given equation (14) leads to:

$$\begin{align} \partial_t^2 (\partial_x \varphi_\phi) \approx -2k_y^2\partial_x^2 (|\phi_k|^2)\, \partial_x \varphi_\phi. \end{align} \tag{ 16 }$$

Dirichlet boundary conditions impose $|\phi_k|^2$ to be concave on average. In fact, with the chosen parameters, $|\phi_k|$ is broad scale and exhibits a single arch of sinusoid, so that $\partial_x^2 (|\phi_k|^2)$ is negative everywhere. The positive expression $-2k_y^2\partial_x^2 (|\phi_k|^2)$ then provides a possible explanation for the observed growth rate, which would scale like $\gamma \approx \sqrt{2}|k_y| \ [\partial_x^2 (|\phi_k|^2)]^{1/2} \sim \sqrt{2}\ |k_yk_x \ \phi_k|$. In this framework, the phase instability results from the self-reinforced shearing of turbulent eddies ($\partial_t (\partial_x \varphi_\phi)$ term) by the Reynolds stress. In the course of the above derivation of the growth rate, we have assumed that the curvature of the Reynolds stress is the dominant term, i.e. $\partial_x^2 (|\phi_k|^2 \partial_x \varphi_\phi ) \approx \partial_x^2 (|\phi_k|^2)\, \partial_x \varphi_\phi$. Although the detailed analysis of the simulation results reveals that this approximation is marginally valid, it turns out that $\partial_x^2 (|\phi_k|^2 \partial_x \varphi_\phi )$ and $\partial_x^2 (|\phi_k|^2)\, \partial_x \varphi_\phi$ exhibit a similar shape, hence providing confidence in the qualitative expression of γ.

This instability is also known as the tilting instability, which transforms convection into sheared flow in 2D ideal fluids [19]. It is responsible for the generation of zonal flows by large amplitude drift (or Rossby) waves governed by the Charney–Hasegawa–Mima equation. In this case, it was shown to exist above a finite amplitude threshold of the electric potential [33]. Consistently, the different dynamical regimes characterizing the 1D model are in agreement with the existence of a threshold, with the instability only developing above a critical magnitude $S_{0,\mathrm{crit}}$ of the driving source term S₀, in between $20\lt S_{0,\mathrm{crit}} \lt30$ (cf figure 7). More recently, the possible role of the tilting instability in the bifurcation toward the H-mode in tokamak plasmas has been addressed [34].

In the next section, the evidence of such a phase instability is reported in the saturated regime of flux-driven ITG turbulence.

In order to identify the respective role of phase gradient and mode amplitude in the time dynamics of the Reynolds stress π, equation (11) is recast as follows:

$$\begin{align} \pi(r,t) & = \mathcal{I}(r,t) \; \varphi^{^{\prime}}(r,t) \; \Delta(r,t) \end{align} \tag{ 17 }$$

$$\begin{align} \textrm{with}\;\; \mathcal{I}(r,t)& = \sum_n\sum_{m\gt0} \frac{2m}{r} |\hat\phi_{mn}|^2 \nonumber \\ \varphi^{^{\prime}}(r,t) & = - \sum_n\sum_{m\gt0} W_{mn}\, \partial_r \varphi^\phi_{mn} \end{align} \tag{ 18 }$$

where $\mathcal{I}$ captures the turbulence intensity. $W_{mn}(r) = \langle|\pi_{mn}|\rangle_t/ \langle\pi\rangle_t$—with π_mn defined by the implicit expression $\pi(r,t) = \sum_n\sum_{m\gt0}\pi_{mn}(r,t)$—is the normalized time-averaged radial profile of the (m, n) spectrum of the Reynolds stress π. This weight aims at discarding contributions $\partial_r \varphi^\phi_{mn}$ to the phase gradient $\varphi^{^{\prime}}(r,t)$ of (m, n) modes that are far from the resonance condition $k_\parallel = 0$, i.e. n + m/q ≈ 0. The weight W_mn (r) at r ≈ 180 is displayed on figure 8. The width of W_mn partly accounts for the fact that GYSELA uses the geometrical angle θ as poloidal coordinate, so that the resonant condition slightly departs from $n+m/q = 0$. Finally, Δ(r, t) accounts for the mismatch between π and the product $\mathcal{I} \varphi^{^{\prime}}$. Note that $|\Delta|$ can be larger than unity due to the use of the weight function in the definition of $\varphi^{^{\prime}}$.

Zoom In Zoom Out Reset image size

The question is whether π is primarily governed by the turbulence intensity $\mathcal{I}$ or by the phase gradient $\varphi^{^{\prime}}$. To this end, we perform the cross-correlation in time between the Reynolds stress π and $\varphi^{^{\prime}}$ on the one hand, and π and $\mathcal{I}$ on the other hand. Both are plotted at each radial location in figure 9. These are computed during the same time window of saturated turbulence reported in the previous section 2. It appears that $\varphi^{^{\prime}}$ exhibits a large correlation with π, up to 0.77 at vanishing time lag. Conversely, at most of the radial positions, the cross-correlation between π and $\mathcal{I}$ remains small or even negative ⁵ . An example of the time dynamics of these three—rescaled—quantities is plotted on figure 10, at r ≈ 180 (r/a ≈ 0.72). We obtain confirmation that both the increase and decay of π are well correlated with those of the phase gradient $\varphi^{^{\prime}}$. Even sometimes, e.g., at t ∼ 35 × 10⁴ and t ∼ 36.7 × 10⁴, π is growing while the amplitude contribution is decreasing. In this case, π is only driven by the growth of the phase gradient.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The critical role of the gradient of the phase in the triggering of a finite poloidal flow has already been acknowledged in the literature [35]. There, from a heuristic model, the authors emphasize that a finite phase curvature alone—i.e. even if turbulence is homogeneous (which translates into $\partial_r |\hat\phi_{mn}|^2 = 0$ in our notations)—can generate a finite Reynolds force. Such an effect was later reported in experimental measurements in a linear plasma device [20]. The properties of the Reynolds stress in the gyrokinetic simulations that we report here extend these predictions and observations in the temporal space, in the sense that the time evolution of π can be entirely governed by that of $\varphi^{^{\prime}}$, even if the turbulence intensity $\mathcal{I}$ remains constant in time.

This result has important implications regarding experimental measurements. Indeed, one of its consequences is that the time dynamics of the Reynolds stress cannot be inferred from that of the amplitude of density—or electric potential—fluctuations only. Indeed, our analysis has shown that π and $\mathcal{I}$ are weakly correlated in general, and can event exhibit anti-correlation. In particular, even though the magnitude of the electric potential fluctuations would be roughly constant, this would not necessarily imply that the magnitude of the Reynolds stress—hence possibly of the zonal flows—remains constant either. Conversely, both can actually exhibit opposite trends, as sometimes observed in figure 10. Our results provide the first evidence from gyrokinetic simulations that the time variation of the Reynolds stress is mainly governed by the phase gradient dynamics, rather than by the magnitude of turbulent fluctuations.

We start from the following minimal system for electrostatic interchange turbulence, reduced to two dimensions by considering flute modes $k_\parallel = 0$. It involves density n and electric potential φ, coupled together via the continuity and charge conservation equations:

$$\begin{align} \partial_t n + [\phi,n] & = S_n + \boldsymbol{\nabla}_\perp.(D \boldsymbol{\nabla}_\perp n) \end{align} \tag{ B.1 }$$

$$\begin{align} \partial_t w + [\phi,w] + g\; \partial_y \log n & = \boldsymbol{\nabla}_\perp.(\nu \boldsymbol{\nabla}_\perp w) \end{align} \tag{ B.2 }$$

$$\begin{align} w & = \nabla_\perp^2\phi. \end{align} \tag{ B.3 }$$

Distances are normalized to the ion thermal Larmor radius $\rho_c = (m_iT)^{1/2}/eB$ ((x, y)→(x, y)/ρ_c ), and time to the ion cyclotron frequency $\omega_c = eB/m_i$ (t → ω_c t), with e the elementary charge, T and B being the constant electron temperature (cold ions are considered) and magnetic field. Finally, g ∼ ρ_c /R accounts for the average curvature of the magnetic field lines [37]. Each field is split into equilibrium and fluctuating parts. A single poloidal (y direction) mode k is then retained for the fluctuating fields:

$$\begin{align} n(x,y,t) & = n_{eq}(x,t) + n_k(x,t)\; \textrm{e}^{iky} + n_k^*(x,t)\; \textrm{e}^{-iky} \nonumber\\ \phi(x,y,t) & = \phi_{eq}(x,t) + \phi_k(x,t)\; \textrm{e}^{iky} + \phi_k^*(x,t)\; \textrm{e}^{-iky}\nonumber\\ w(x,y,t) & = w_{eq}(x,t) + w_k(x,t)\; \textrm{e}^{iky} + w_k^*(x,t)\; \textrm{e}^{-iky}. \end{align} \tag{ B.4 }$$

Equilibrium (y-averaged) quantities $n_\mathrm{eq}$, $\phi_\mathrm{eq}$ and $w_\mathrm{eq}$ are real ($\in \mathbb{R}$), while n_k , φ_k and w_k are complex fields ($\in \mathbb{C}$). This approach is in marked contrast with the derivation of [32], where the phase shift between density and potential fluctuations was assumed to be clamped to π/2, hence maximizing turbulent transport. Introducing the equilibrium velocity:

$$\begin{align} V_\mathrm{eq} = \partial_x \phi_\mathrm{eq} \end{align} \tag{ B.5 }$$

one finally obtains the following 4-field system of equations:

$$\begin{align} \partial_t n_\mathrm{eq} & = 2k\; \partial_x \left[\Im( n_k\phi_k^*)\right] + \partial_x(D_0 \partial_x n_\mathrm{eq}) + S \nonumber\\ \partial_t V_\mathrm{eq} & = 2k\; \partial_x \left[\Im( \phi_k^* \, \partial_x\phi_k )\right] + \partial_x(\nu_0 \partial_x V_\mathrm{eq}) \nonumber\\ \partial_t n_k & = i k \left(\phi_k\, \partial_xn_\mathrm{eq} - V_\mathrm{eq}\,n_k \right) - D_1k^2n_k + \partial_x(D_1 \partial_x n_k)\\ \partial_t w_k & = i k \left(\phi_k\, \partial_x^2V_\mathrm{eq} - V_\mathrm{eq}\, w_k - g\frac{n_k}{n_\mathrm{eq}} \right) - \nu_1k^2w_k\nonumber\\ & \quad + \partial_x(\nu_1 \partial_x w_k)\nonumber\\ w_k & = \left(\partial_x^2 - k^2\right) \phi_k. \nonumber \end{align} \tag{ B.6 }$$

The turbulent particle flux and the Reynolds stress then read:

$$\begin{align} \Gamma & = -2k\ \Im( n_k\phi_k^*) \end{align} \tag{ B.7 }$$

$$\begin{align} \pi & = -2k\ \Im( \phi_k^* \, \partial_x\phi_k ). \end{align} \tag{ B.8 }$$

Alternative forms of the Reynolds force may prove more appropriate for numerical schemes: $\partial_x \pi = -2k\ \Im(\phi^*_k w_k) = -2k\ \Im(\phi_k^*\, \partial_x^2\phi_k)$.

The system (B.6) was initially derived in [31]. Notice that the system proposed in [38] differs from this one by neglecting the radial curvature of the potential with respect to its poloidal one, i.e., it assumes $|k^2 \phi_k| \gg |\partial_x^2 \phi_k|$. In this case, vorticity and potential are in phase. Actually, this latter system proves to be numerically unstable, leading to the transient excitation of fluctuations at spatial scales down to the grid size.

An alternative form of the previous system can be derived by introducing the magnitude and phase of the retained Fourier component k of the fluctuations.

$$\begin{align*} \hat{f}(x,t) = A_f(x,t)\, \ e^{i\varphi_f(x,t)} \end{align*}$$

with $A_f \equiv |\,\,\hat f\,| \geqslant 0$ and $\varphi_f \in \mathbb{R}$. One further introduces the following phase shifts (for any f, g ∈ n, φ, w):

$$\begin{align} \Delta\varphi_{fg} \equiv \varphi_f - \varphi_g. \end{align} \tag{ B.9 }$$

System (B.6) can then be recast in terms of amplitudes and phases. The equations for the mean quantities read:

$$\begin{align} \partial_t n_\mathrm{eq} & = 2 k\ \partial_x \left(A_nA_\phi\, \sin \Delta\varphi_{n\phi}\right) + \partial_x(D_0\, \partial_x n_\mathrm{eq}) + S_n \end{align} \tag{ B.10 }$$

$$\begin{align} \partial_t V_\mathrm{eq} & = 2 k\ \partial_x \left(A_\phi^2\, \partial_x \varphi_\phi \right) + \partial_x(\nu_0\, \partial_x V_\mathrm{eq}) -\mu\left( V_\mathrm{eq} - V_0 \right). \end{align} \tag{ B.11 }$$

Notice that the nonlinear particle flux and the Reynolds stress take the following simple forms:

$$\begin{align} \Gamma & = -2 k\ A_nA_\phi\, \sin \Delta\varphi_{n\phi} \end{align} \tag{ B.12 }$$

$$\begin{align} \pi & = -2k\ A_\phi^2\, \partial_x \varphi_\phi. \end{align} \tag{ B.13 }$$

The time evolution of the fluctuations is given by:

$$\begin{align} \partial_t A_n & = k A_\phi\, \partial_x n_\mathrm{eq}\, \sin \Delta\varphi_{n\phi} - (D_1k^2+\alpha_n A_n^2)\, A_n \nonumber \\[4pt] & \quad + \partial_x( D_1 \partial_x A_n) - D_1A_n (\partial_x \varphi_n)^2 \end{align} \tag{ B.14 }$$

$$\begin{align} \partial_t A_w & = k A_\phi \partial_x^2 V_\mathrm{eq}\, \sin \Delta\varphi_{w\phi} + gk\frac{A_n}{n_\mathrm{eq}}\sin \Delta\varphi_{nw}\nonumber \\[5pt] & \quad - (\nu_1k^2+\alpha_w A_w^2)\, A_w + \partial_x(\nu_1 \partial_x A_w)-\nu_1A_w(\partial_x \varphi_w)^2 \end{align} \tag{ B.15 }$$

with the following relationship between A_w and $A_\phi$:

$$\begin{align} A_w = \left\{\partial_x^2A_\phi - \left[k^2 + (\partial_x\varphi_\phi)^2 \right] A_\phi \right\} \ [\cos \Delta\varphi_{w\phi} + \tan \Delta\varphi_{w\phi}]. \end{align} \tag{ B.16 }$$

The phases are governed by the following equations:

$$\begin{align} \partial_t \varphi_n & = k\frac{A_\phi}{A_n}\partial_x n_\mathrm{eq} \cos \Delta\varphi_{n\phi} -kV_\mathrm{eq} \nonumber\\[3pt] & \quad + \partial_x(D_1 \partial_x \varphi_n) + 2D_1\partial_x (\ln A_n) \partial_x\varphi_n \end{align} \tag{ B.17 }$$

$$\begin{align} \partial_t \varphi_w & = k\frac{A_\phi}{A_w}\partial_x^2 V_\mathrm{eq} \cos \Delta\varphi_{w\phi}-gk\frac{A_n}{n_\mathrm{eq}A_w} \cos \Delta\varphi_{nw}-kV_\mathrm{eq} \nonumber\\[4pt] &\quad + \partial_x(\nu_1 \partial_x \varphi_w) + 2\nu_1\partial_x (\ln A_w) \partial_x\varphi_w \end{align} \tag{ B.18 }$$

with

$$\begin{align} \tan \Delta\varphi_{w\phi} = \frac{2\partial_xA_\phi\partial_x\varphi_\phi + A_\phi\partial_x^2 \varphi_\phi}{ \partial_x^2A_\phi - \left[k^2 + (\partial_x\varphi_\phi)^2 \right] A_\phi} \end{align} \tag{ B.19 }$$