Gyrokinetic theory of turbulent acceleration and momentum conservation in tokamak plasmas

The intrinsic rotation can be generated by a self-organization process in turbulent atmosphere. This phenomenon is widely discovered in nature, such as the planetary atmosphere [1], solar tachocline [2], as well as in laboratory, such as the intrinsic rotation in tokamak plasmas [3]. The intrinsic rotation in tokamak plasmas has been extensively studied for decades since rotation is helpful for stabilizing magnetohydrodynamic (MHD) instabilities such as resistive wall modes [4, 5] and neoclassical tearing modes [6], and plays an important role in improving plasma confinement by suppressing micro-turbulence [7–9]. In addition, another idea to control tokamak plasma turbulence via artificial intelligence techniques with using artificial neutral networks was also proposed recently [10–12].

What the mechanisms account for the intrinsic rotation drive is still an open question. We proposed a new mechanism named turbulent acceleration for the intrinsic parallel rotation drive by electrostatic ion temperature gradient (ITG) turbulence [13]. The turbulent acceleration cannot be written as a divergence of stress term, which is different from the physics of residual stress, and so acts as a local source or sink. In other words, the turbulent acceleration is an effective volume-force, while the divergence of residual stress is a kind of surface-force. We also proposed a possible theoretical explanation for the experimental observation of electron cyclotron heating (ECH) induced decrease of co-current core toroidal rotation via the comparison between the turbulent acceleration driven by ITG turbulence and that driven by collisionless trapped electron mode (CTEM) turbulence [14]. Then, we extended our theory to electromagnetic ITG turbulence and found the importance of electromagnetic effects on intrinsic parallel rotation drive [15]. Finally, we demonstrated that the conserved quantity corresponding to axisymmetry is the total canonical momentum or total momentum carried by both particles and electromagnetic fields, and the presence of turbulent acceleration does not imply its confliction with momentum conservation [16].

The rest of this paper is organized as follows. In section 2, the new mechanism named turbulent acceleration for electrostatic ITG turbulence and the comparison of turbulent acceleration between ITG turbulence and CTEM turbulence are presented. Then, an extension to electromagnetic ITG turbulence and some discussions on relevant experiments are given in section 3. The relationship between the turbulent acceleration and the momentum conservation law is clarified in section 4. Finally, we make a summary and discuss our future work in section 5.

2.1. Electrostatic ITG turbulence

Starting from the conservative form of the nonlinear electrostatic gyrokinetic equation in the continuity form [17]

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}(F{B}_{\parallel }^{* })+{\rm{\nabla }}\cdot \left(\displaystyle \frac{{\rm{d}}{\boldsymbol{R}}}{{\rm{d}}t}F{B}_{\parallel }^{* }\right)+\displaystyle \frac{\partial }{\partial {v}_{\parallel }}\left(\displaystyle \frac{{\rm{d}}{v}_{\parallel }}{{\rm{d}}t}F{B}_{\parallel }^{* }\right)=0,\end{eqnarray} \tag{ 1 }$$

with gyrocenter equations of motion

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}{\boldsymbol{R}}}{{\rm{d}}t}={v}_{\parallel }\hat{{\boldsymbol{b}}}+\displaystyle \frac{c}{e{B}_{\parallel }^{* }}\hat{{\boldsymbol{b}}}\times (e{\rm{\nabla }}\langle \langle \delta {\phi }_{{\rm{gc}}}\rangle \rangle +\mu {\rm{\nabla }}B\\ &&\,\,+{m}_{{\rm{i}}}{v}_{\parallel }^{2}\hat{{\boldsymbol{b}}}\cdot {\rm{\nabla }}\hat{{\boldsymbol{b}}}),\end{eqnarray} \tag{ 2 }$$

and

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}{v}_{\parallel }}{{\rm{d}}t}=-\displaystyle \frac{{\boldsymbol{B}}* }{{m}_{{\rm{i}}}{B}_{\parallel }^{* }}\cdot (e{\rm{\nabla }}\langle \langle \delta {\phi }_{{\rm{gc}}}\rangle \rangle +\mu {\rm{\nabla }}B).\end{eqnarray} \tag{ 3 }$$

Here, $F=F({\boldsymbol{R}},\,\mu ,\,{v}_{\parallel },\,t)$ is the gyrocenter distribution function, $\mu$ is the gyrocenter magnetic moment and ${v}_{\parallel }={\boldsymbol{v}}\cdot \hat{{\boldsymbol{b}}}$ is the parallel velocity of ions' gyrocenter, ${B}_{\parallel }^{* }\equiv \hat{{\boldsymbol{b}}}\cdot {\boldsymbol{B}}*$ is the Jacobian of the transformation from particle phase space to the gyrocenter phase space with ${\boldsymbol{B}}* =\,{\boldsymbol{B}}+\tfrac{{m}_{{\rm{i}}}c}{e}{v}_{\parallel }{\rm{\nabla }}\times \hat{{\boldsymbol{b}}},$ and ${\boldsymbol{B}}$ is the equilibrium magnetic field with $\hat{{\boldsymbol{b}}}=\tfrac{{\boldsymbol{B}}}{B}.$ ${m}_{{\rm{i}}}$ is the ion mass, $\delta {\phi }_{{\rm{gc}}}$ is the gyrocenter electric potential fluctuation, and $\langle \langle \cdots \rangle \rangle$ denotes gyroaveraging. In the following, we denote $\langle \langle \delta {\phi }_{{\rm{gc}}}\rangle \rangle =\delta \phi$ since a long wavelength approximation ${k}_{\perp }^{2}{\rho }_{{\rm{i}}}^{2}\ll 1$ is used, where ${k}_{\perp }$ is the perpendicular wave number and ${\rho }_{{\rm{i}}}=\displaystyle \frac{{v}_{{\rm{t}}{\rm{h}}{\rm{i}}}}{{{\rm{\Omega }}}_{{\rm{i}}}}\,$ is the ion Larmor radius with ${v}_{{\rm{thi}}}=\sqrt{\tfrac{{T}_{{\rm{i}}}}{{m}_{{\rm{i}}}}}$ being the ion thermal velocity and ${{\rm{\Omega }}}_{{\rm{i}}}=\tfrac{eB}{{m}_{{\rm{i}}}}$ being the ion cyclotron frequency. In this paper, the index $\parallel$ refers to the components parallel to the equilibrium magnetic field, and the index $\perp$ refers to the components perpendicular to the equilibrium magnetic field. Note that the usual experimentally measured quantity is the ion parallel flow velocity rather than the ion parallel momentum density, so we focus on the mean parallel flow velocity equation. By taking the zeroth order and the first order moments of equation (1), decoupling the density from momentum density equation and taking flux average of parallel flow equation, we can finally obtain the evolution equation of the mean parallel flow velocity [13]

$$\begin{eqnarray}&&\displaystyle \frac{\partial \langle {U}_{\parallel }\rangle }{\partial t}+{\rm{\nabla }}\cdot {{\rm{\Pi }}}_{r,\parallel }={a}_{\parallel }.\end{eqnarray} \tag{ 4 }$$

Here, ${U}_{\parallel }$ is the ion parallel flow velocity and $\langle \cdots \rangle$ denotes flux average. ${{\rm{\Pi }}}_{r,\parallel }=\langle \delta {{\bf{v}}}_{E\times B,r}\delta {U}_{\parallel }\rangle$ is the usual parallel Reynolds stress with $\delta {{\bf{v}}}_{E\times B}=\tfrac{c\hat{{\boldsymbol{b}}}\times {\rm{\nabla }}\delta \phi }{B}$ being the fluctuating ${\boldsymbol{E}}\times {\boldsymbol{B}}$ drift velocity and $\delta {U}_{\parallel }$ being the parallel flow fluctuation, which consists of three terms including diffusion, convection, and residual stress ${{\rm{\Pi }}}_{r,\parallel }^{{\rm{res}}}.$ The parallel Reynolds stress has been intensively investigated [18, 19], ${a}_{\parallel }=\tfrac{1}{{m}_{{\rm{i}}}{n}_{0}}\langle \delta \hat{n}\hat{{\boldsymbol{b}}}\cdot {\rm{\nabla }}\delta {P}_{\parallel {\rm{i}}}\rangle$ is the turbulent acceleration with $\delta \hat{n}=\tfrac{\delta n}{{n}_{0}}$ being the normalized density fluctuation and $\delta {P}_{\parallel {\rm{i}}}$ being the ion pressure fluctuation with ${P}_{\parallel {\rm{i}}}=2\pi \displaystyle \int {\rm{d}}{v}_{\parallel }{\rm{d}}\mu {B}_{\parallel }^{* }F{v}_{\parallel }^{2}={n}_{{\rm{i}}}{T}_{{\rm{i}}}.$ The turbulent acceleration cannot be written as a divergence of the parallel Reynolds stress, and plays a role of local source/sink of parallel rotation. Its existence does not depend on the magnetic geometry and it is independent of the mean parallel velocity or the gradient of mean parallel velocity. Therefore, it is a new candidate mechanism for the origin of intrinsic rotation.

For electrostatic ITG turbulence, an adiabatic electron response is assumed, i.e., $\delta \hat{n}=\delta {\hat{\phi }}_{k}$ with $\delta {\hat{\phi }}_{k}=\tfrac{e\delta {\phi }_{k}}{{T}_{{\rm{e}}}}.$ By linearizing the ion pressure equation we can obtain the ion pressure fluctuation

$$\begin{eqnarray}&&\delta {P}_{\parallel {\rm{i}}}={P}_{\parallel {\rm{i}}}\displaystyle \frac{{\omega }_{* {\rm{e}}}}{{\omega }_{{\rm{r}}}}\left(1-{\rm{i}}\displaystyle \frac{| {\gamma }_{k}| }{{\omega }_{{\rm{r}}}}\right)\delta {\hat{\phi }}_{k},\end{eqnarray} \tag{ 5 }$$

where ${\omega }_{{\rm{r}}}$ and ${\gamma }_{k}$ are real frequency and linear growth rate, respectively. Then, a quasilinear expression for the turbulent acceleration can be written as [13]

$$\begin{eqnarray}&&{a}_{\parallel }=\displaystyle \frac{1}{{m}_{{\rm{i}}}{n}_{0}}\langle \delta \hat{n}\hat{{\boldsymbol{b}}}\cdot {\rm{\nabla }}\delta {P}_{\parallel {\rm{i}}}\rangle \approx \tau (1+{\eta }_{{\rm{i}}})\displaystyle \frac{1}{{L}_{n}}\displaystyle \sum _{k}{{\rm{\Pi }}}_{{\rm{E}}{\rm{S}},k},\end{eqnarray} \tag{ 6 }$$

where $\tau =\tfrac{{T}_{{\rm{i}}}}{{T}_{{\rm{e}}}},$ ${\eta }_{{\rm{i}}}=\tfrac{{L}_{n}}{{L}_{T{\rm{i}}}}$ with ${L}_{n}$ being the density gradient scale length and ${L}_{T{\rm{i}}}$ being the ITG scale length. ${{\rm{\Pi }}}_{r,\parallel }^{{\rm{res}}}\,=\displaystyle \sum _{k}{{\rm{\Pi }}}_{{\rm{E}}{\rm{S}},k}=\displaystyle \sum _{k}{\rho }_{{\rm{s}}}{c}_{{\rm{s}}}^{3}\tfrac{| {\gamma }_{k}| }{{\omega }_{{\rm{r}}}^{2}}{k}_{\theta }{k}_{\parallel }{| \delta {\hat{\phi }}_{k}| }^{2}$ is electrostatic residual stress which is consistent with previous result [19], ${c}_{{\rm{s}}}=\sqrt{\tfrac{{T}_{{\rm{e}}}}{{m}_{{\rm{i}}}}}$ is the ion acoustic velocity and ${\rho }_{{\rm{s}}}=\tfrac{{c}_{{\rm{s}}}}{{{\rm{\Omega }}}_{{\rm{i}}}}.$ By taking the parallel symmetry breaking induced by intensity gradient [19], i.e., ${I}_{k}(x)={| \delta {\hat{\phi }}_{k}| }^{2}(x)={I}_{k}(0)+x\left(\tfrac{\partial {I}_{k}}{\partial x}\right)$ and ${k}_{\parallel }=\tfrac{{k}_{\theta }x\hat{{\rm{s}}}}{q{R}_{0}},$ where $\hat{{\rm{s}}}$ is the magnetic shear, $q$ is the safety factor, ${R}_{0}$ is the major radius, $x={r}_{m,n}-r,$ with ${r}_{m,n}$ being the radial location of the resonant surface, we can obtain the turbulent acceleration

$$\begin{eqnarray}&&{a}_{\parallel }\approx \tau (1+{\eta }_{{\rm{i}}})\displaystyle \frac{1}{{L}_{n}}\displaystyle \frac{\hat{{\rm{s}}}}{q{R}_{0}}\displaystyle \sum _{k}{\rho }_{{\rm{s}}}{c}_{{\rm{s}}}^{3}\displaystyle \frac{| {\gamma }_{k}| }{{\omega }_{{\rm{r}}}^{2}}{k}_{\theta }^{2}{x}^{2}\displaystyle \frac{\partial {I}_{k}}{\partial x}.\end{eqnarray} \tag{ 7 }$$

It is obvious that the turbulent acceleration can provide a co-current intrinsic rotation drive for positive magnetic shear $(\hat{s}\gt 0)$ and positive intensity gradient $\left(\tfrac{\partial {I}_{k}}{\partial x}\gt 0\right).$ For convenience of comparison, the divergence of the residual stress obtained in [19] is written as

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot {{\rm{\Pi }}}_{r,\parallel }^{{\rm{res}}}\approx \pm \displaystyle \frac{1}{L}\displaystyle \frac{\hat{{\rm{s}}}}{q{R}_{0}}\displaystyle \sum _{k}{\rho }_{{\rm{s}}}{c}_{{\rm{s}}}^{3}\displaystyle \frac{| {\gamma }_{k}| }{{\omega }_{{\rm{r}}}^{2}}{k}_{\theta }^{2}{x}^{2}\displaystyle \frac{\partial {I}_{k}}{\partial x}.\end{eqnarray} \tag{ 8 }$$

Here, $L$ is the scale length of the radial variation of the residual stress. Thus, the ratio of the turbulent acceleration to the divergence of the residual stress is $\tau (1+{\eta }_{{\rm{i}}})\tfrac{L}{{L}_{n}},$ which means the turbulent acceleration can be comparable to the divergence of the residual stress, depending upon $\tau (1+{\eta }_{{\rm{i}}})$ and $\tfrac{L}{{L}_{n}}.$ Hence, the new mechanism of turbulent acceleration is qualitatively different from the usual Reynolds stress, but it is quantitatively comparable to the divergence of residual stress.

2.2. Electrostatic CTEM turbulence

For electromagnetic turbulence, the shear component of magnetic perturbation, i.e., $\delta {A}_{\parallel },$ is considered [15]. Starting from the nonlinear electromagnetic gyrokinetic equation with gyrocenter equations of motion in the symplectic formulation (i.e., ${v}_{\parallel }$ representation)

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}{\boldsymbol{R}}}{{\rm{d}}t}={v}_{\parallel }\hat{{\boldsymbol{b}}}* +\displaystyle \frac{c}{eB}\hat{{\boldsymbol{b}}}\times e{\rm{\nabla }}\langle \langle \delta \phi \rangle \rangle \end{eqnarray} \tag{ 9 }$$

and

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}{{\boldsymbol{v}}}_{\parallel }}{{\rm{d}}t}=-\displaystyle \frac{{\boldsymbol{e}}}{{m}_{{\rm{i}}}}\left(\hat{{\boldsymbol{b}}}* \cdot {\rm{\nabla }}\langle \langle \delta \phi \rangle \rangle +\displaystyle \frac{1}{c}\displaystyle \frac{\langle \langle \partial \delta {A}_{\parallel }\rangle \rangle }{\partial t}\right).\end{eqnarray} \tag{ 10 }$$

Here, uniform equilibrium magnetic field ${\boldsymbol{B}}$ is assumed. ${\boldsymbol{B}}* =\,{\boldsymbol{B}}+\delta {{\boldsymbol{B}}}_{\perp },$ $\delta {{\boldsymbol{B}}}_{\perp }=-\hat{{\boldsymbol{b}}}\times {\rm{\nabla }}\delta {A}_{\parallel }.$ $\hat{{\boldsymbol{b}}}* =\,\hat{{\boldsymbol{b}}}+\delta {\hat{{\boldsymbol{B}}}}_{\perp }$ with $\delta {\hat{{\boldsymbol{B}}}}_{\perp }=\tfrac{\delta {{\boldsymbol{B}}}_{\perp }}{B}.$ By taking the process similar to electrostatic case, we can obtain the mean field equation of the parallel flow velocity

$$\begin{eqnarray}&&\displaystyle \frac{\partial {U}_{\parallel }}{\partial t}+{\rm{\nabla }}\cdot \langle \delta {{\bf{v}}}_{E\times B,r}\delta {U}_{\parallel }\rangle +\displaystyle \frac{1}{{m}_{{\rm{i}}}{n}_{0}}{\rm{\nabla }}\cdot \langle \delta {\hat{{\boldsymbol{B}}}}_{r}\delta {P}_{\parallel {\rm{i}}}\rangle \\ &&\,+\displaystyle \frac{e}{{m}_{{\rm{i}}}}{\rm{\nabla }}\cdot \langle \delta {\hat{{\boldsymbol{B}}}}_{r}\delta \phi \rangle =\displaystyle \frac{1}{{m}_{{\rm{i}}}{n}_{0}}[\langle \delta \hat{n}\hat{{\boldsymbol{b}}}\cdot {\rm{\nabla }}\delta {P}_{\parallel {\rm{i}}}\rangle \\ &&\,+\langle \delta \hat{n}\delta {\hat{{\boldsymbol{B}}}}_{r}\cdot {\rm{\nabla }}{P}_{\parallel {\rm{i}}}\rangle ].\end{eqnarray} \tag{ 11 }$$

Here, $\langle \delta {\hat{{\boldsymbol{B}}}}_{r}\delta {P}_{\parallel {\rm{i}}}\rangle$ is known as the kinetic stress [26, 27] and it is somewhat analogous to the kinetic dynamo physics which is a transport process of electron parallel momentum (current) [27, 28]. The kinetic stress was experimentally identified to drive the intrinsic parallel plasma flow in MST [26]. $\langle \delta {\hat{{\boldsymbol{B}}}}_{r}\delta \phi \rangle$ is the cross Maxwell stress, which was also presented in [29]. Both kinetic stress and cross Maxwell stress are surface force since they enter the parallel flow equation via their divergence. This is similar to the usual Reynolds stress $\langle \delta {{\bf{v}}}_{E\times B,r}\delta {U}_{\parallel }\rangle .$ While, the turbulent acceleration on the right hand side of equation (11) is an effective volume-force, because it cannot be recast into the divergence of stress term. $\tfrac{1}{{m}_{{\rm{i}}}{n}_{0}}\langle \delta \hat{n}\hat{{\boldsymbol{b}}}\cdot {\rm{\nabla }}\delta {P}_{\parallel {\rm{i}}}\rangle$ is the turbulent acceleration driven by parallel gradient of ion pressure fluctuation, which has been investigated for electrostatic turbulence as mentioned above. $\tfrac{1}{{m}_{{\rm{i}}}{n}_{0}}\langle \delta \hat{n}\delta {\hat{{\boldsymbol{B}}}}_{r}\rangle \cdot {\rm{\nabla }}{P}_{\parallel {\rm{i}}}$ is a new turbulent acceleration related to the correlation between density and magnetic fluctuations and the equilibrium ion pressure gradient as well. Combining these two terms, the total electromagnetic turbulent acceleration can be rewritten as $\tfrac{1}{{m}_{{\rm{i}}}{n}_{0}}\langle \delta \hat{n}\,\delta (\hat{{\boldsymbol{b}}}* \cdot {\rm{\nabla }}{P}_{\parallel {\rm{i}}})\rangle ,$ which represents the correlation between density fluctuation and perturbed pressure gradient along the total magnetic field.

Quasilinear estimates of all the intrinsic rotation drive can be calculated for the typical parameters of the pedestal top region on DIII-D where the dominant instability is electromagnetic ITG turbulence [30]. The total stress force can be written as

$$\begin{eqnarray}&&{{\rm{\Lambda }}}_{{\rm{tot}}}=\mp \displaystyle \frac{1}{L^{\prime} }\displaystyle \sum _{k}{{\rm{\Pi }}}_{{\rm{E}}{\rm{S}},k}(\mathrm{1-23}{\hat{\beta }}_{e}),\end{eqnarray} \tag{ 12 }$$

and the total turbulent acceleration is

$$\begin{eqnarray}&&{a}_{{\rm{tot}}}=\displaystyle \frac{1}{{L}_{n}}\displaystyle \sum _{k}{{\rm{\Pi }}}_{{\rm{E}}{\rm{S}},k}16\left(1+\displaystyle \frac{17}{4}{\hat{\beta }}_{e}\right).\end{eqnarray} \tag{ 13 }$$

Here, $L^{\prime}$ is the scale length of the radial variation of the stresses and $\displaystyle \sum _{k}{{\rm{\Pi }}}_{{\rm{E}}{\rm{S}},k}$ is the electrostatic residual stress mentioned in section 2.1. ${\hat{\beta }}_{{\rm{e}}}={\beta }_{{\rm{e}}}\tfrac{{q}^{2}{R}_{0}^{2}}{{L}_{n}^{2}}\approx 0.144$ for the typical parameters of the pedestal top region on DIII-D with ${\beta }_{{\rm{e}}}=\tfrac{8\pi {n}_{0}{T}_{{\rm{e}}}}{{B}^{2}}\approx 0.4 \% ,$ and the terms related to ${\hat{\beta }}_{{\rm{e}}}$ come from electromagnetic effects. From equations (12) and (13), we can see that the electromagnetic effects can reduce the electrostatic stress force and even reverse it, but enhance the turbulent acceleration, and the total turbulent acceleration is comparable to the total stress force. Therefore, taking into account the electromagnetic effects and the turbulent acceleration is important for more accurate understanding and prediction of intrinsic rotation in pedestal region.

Now, we discuss some possible relevant experimental observations. In the MST, the kinetic stress related to density fluctuations $\langle \delta \hat{n}\delta {\hat{{\boldsymbol{B}}}}_{r}\rangle$ was directly measured [26], so the turbulent acceleration related to the kinetic stress (i.e., $\tfrac{1}{{m}_{{\rm{i}}}{n}_{0}}\langle \delta \hat{n}\delta {\hat{{\boldsymbol{B}}}}_{r}\rangle \cdot {\rm{\nabla }}{P}_{\parallel {\rm{i}}}$) can be also easily obtained there. This turbulent acceleration $\tfrac{1}{{m}_{{\rm{i}}}{n}_{0}}\langle \delta \hat{n}\delta {\hat{{\boldsymbol{B}}}}_{r}\rangle \cdot {\rm{\nabla }}{P}_{\parallel {\rm{i}}}$ can provide counter-current rotation drive for normal ion pressure profile $({\rm{\nabla }}{P}_{\parallel {\rm{i}}}\lt 0)$ and positive kinetic stress. Taking this turbulent acceleration into account in the intrinsic flow drive of MST will lead to the reversal point of flow drive closer to that of the plasma parallel flow [26], which may result in better agreement between the intrinsic flow drive and the observed parallel flow profile.

As mentioned above, turbulent acceleration is obtained by decoupling density from the gyrocenter parallel momentum density equation and then taking flux average of flow velocity, i.e.,

$$\begin{eqnarray}&&\displaystyle \frac{\partial \langle {U}_{\parallel }\rangle }{\partial t}=\left\langle \displaystyle \frac{1}{n}\left[\displaystyle \frac{\partial (n{U}_{\parallel })}{\partial t}-{U}_{\parallel }\displaystyle \frac{\partial n}{\partial t}\right]\right\rangle .\end{eqnarray} \tag{ 14 }$$

Besides, there is an alternative method to derive the mean parallel flow velocity equation by taking flux average of the gyrocenter parallel momentum density first, and then decoupling flow velocity and density as follows

$$\begin{eqnarray}&&\displaystyle \frac{\partial \langle {U}_{\parallel }\rangle }{\partial t}=\displaystyle \frac{1}{{n}_{0}}\left[\displaystyle \frac{\langle \partial n{U}_{\parallel }\rangle }{\partial t}-\langle {U}_{\parallel }\rangle \displaystyle \frac{\partial {n}_{0}}{\partial t}-\left\langle \delta {U}_{\parallel }\displaystyle \frac{\partial }{\partial t}\delta n\right\rangle \right.\\ &&\,\,\,\,\,-\left.\left\langle \delta n\displaystyle \frac{\partial }{\partial t}\delta {U}_{\parallel }\right\rangle \right].\end{eqnarray} \tag{ 15 }$$

Through equations (14) and (15), we can obtain the same mean field equation of the parallel flow velocity (i.e., equation (11)) regardless of the sequence of taking flux average and decoupling the parallel flow from parallel momentum density [16]. Therefore, the presence of the turbulent acceleration does not result from regrouping terms as one may wonder.

By taking summation of gyrocenter momentum density equation over both species and using quasi-neutrality condition for electron density ${n}_{{\rm{e}}}$ and ion gyrocenter density ${\bar{n}}_{{\rm{i}}},$ i.e., ${n}_{{\rm{e}}}-{\bar{n}}_{{\rm{i}}}={n}_{{\rm{pol}}}$ with ${n}_{{\rm{pol}}}={{\rm{\nabla }}}_{\perp }\cdot \left(\tfrac{{c}^{2}{m}_{{\rm{i}}}{n}_{0}}{e{B}^{2}}{{\rm{\nabla }}}_{\perp }\delta \phi \right)$ being ion polarization density, we can obtain the gyrocenter momentum conservation equation after some algebra [16]

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}\displaystyle \sum _{s}\langle {\bar{p}}_{\parallel s}\rangle +{{\rm{\nabla }}}_{r}\cdot \langle \delta {{\rm{v}}}_{E\times B,r}\delta ({m}_{{\rm{i}}}{\bar{n}}_{{\rm{i}}}{\bar{U}}_{\parallel {\rm{i}}})+\delta {\hat{{\boldsymbol{B}}}}_{r}\delta \bar{P}\rangle \\ &&\,-{{\rm{\nabla }}}_{r}\cdot \left\langle \displaystyle \frac{{c}^{2}{m}_{{\rm{i}}}{n}_{0}}{{B}^{2}}{{\rm{\nabla }}}_{r}\delta \phi \hat{{\boldsymbol{b}}}\cdot {\rm{\nabla }}\delta \phi \right\rangle =0,\end{eqnarray} \tag{ 16 }$$

where s includes ions and electrons, $\displaystyle {\sum }_{s}\langle {\bar{p}}_{\parallel s}\rangle$ = $\left\langle {m}_{{\rm{i}}}{\bar{n}}_{{\rm{i}}}{\bar{U}}_{\parallel {\rm{i}}}-{n}_{{\rm{pol}}}\displaystyle \frac{e}{c}\delta {A}_{\parallel }\right\rangle$ with ${\bar{p}}_{\parallel s}={m}_{s}{\bar{n}}_{s}{\bar{U}}_{\parallel s}+{\bar{n}}_{s}\displaystyle \frac{{e}_{s}}{c}\delta {A}_{\parallel }$ being the gyrocenter canonical parallel momentum of one species including the ordinary gyrocenter kinematic momentum and the gyrocenter magnetic vector momentum. The ordinary kinematic momentum carried by electrons was neglected because of small electron mass as compared to ion mass $({m}_{{\rm{e}}}\ll {m}_{{\rm{i}}}).$ $-{n}_{{\rm{pol}}}\tfrac{e}{c}\delta {A}_{\parallel }$ is the total gyrocenter magnetic vector momentum resulting from the ion polarization density ${n}_{{\rm{pol}}}.$ Then, we obtain the conservation equation of the total gyrocenter canonical parallel momentum carried by ions and electrons. The terms under the divergence are radial flux of the ion gyrocenter parallel kinematic momentum, total kinetic stress and electric Maxwell stress, respectively. By using the relationship between total gyrocenter magnetic vector momentum and the electromagnetic fields momentum

$$\begin{eqnarray}\displaystyle \frac{\partial \left\langle -{n}_{{\rm{pol}}}\displaystyle \frac{e}{c}\delta {A}_{\parallel }\right\rangle }{\partial t} & = & \displaystyle \frac{\partial \left\langle {m}_{{\rm{i}}}{n}_{0}\displaystyle \frac{{c}^{2}}{{B}^{2}}{(\delta {\boldsymbol{E}}\times \delta {\boldsymbol{B}})}_{\parallel }\right\rangle }{\partial t}\\ & & -{\rm{\nabla }}\cdot \left\langle \displaystyle \frac{{c}^{2}{m}_{{\rm{i}}}{n}_{0}}{{B}^{2}}{{\rm{\nabla }}}_{r}\delta \phi \displaystyle \frac{1}{c}\displaystyle \frac{\partial \delta {A}_{\parallel }}{\partial t}\right\rangle ,\end{eqnarray} \tag{ 17 }$$

equation (16) can be rewritten as

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}\left\langle {m}_{{\rm{i}}}{\bar{n}}_{{\rm{i}}}{\bar{U}}_{\parallel {\rm{i}}}+{m}_{{\rm{i}}}{n}_{0}\displaystyle \frac{{c}^{2}}{{B}^{2}}{(\delta {\boldsymbol{E}}\times \delta {\boldsymbol{B}})}_{\parallel }\right\rangle \\ &&\,+{{\rm{\nabla }}}_{r}\cdot \langle \delta {{\rm{v}}}_{E\times B,r}\delta ({m}_{{\rm{i}}}{\bar{n}}_{{\rm{i}}}{\bar{U}}_{\parallel {\rm{i}}})+\delta {\hat{{\boldsymbol{B}}}}_{r}\delta \bar{P}\rangle \\ &&\,-{{\rm{\nabla }}}_{r}\cdot \left\langle \displaystyle \frac{{c}^{2}{m}_{{\rm{i}}}{n}_{0}}{{B}^{2}}{{\rm{\nabla }}}_{r}\delta \phi \left(\hat{{\boldsymbol{b}}}\cdot {\rm{\nabla }}\delta \phi +\displaystyle \frac{1}{c}\displaystyle \frac{\partial \delta {A}_{\parallel }}{\partial t}\right)\right\rangle =0.\end{eqnarray} \tag{ 18 }$$

This is another conservative format of momentum which is the total gyrocenter parallel momentum density including the ion gyrocenter kinematic momentum and electromagnetic fields momentum. The pieces under the divergence are radial flux of the ion gyrocenter kinematic momentum, total kinetic stress, and electric Maxwell stress including both parallel electrostatic field and inductive electric field.

These two gyrocenter momentum conservation equations mentioned above (i.e., equations (16) and (18)) are equivalent. We can see that the conserved quantity can be either the total gyrocenter parallel canonical momentum carried by both species or the total gyrocenter parallel momentum including the ion gyrocenter kinematic momentum and the electromagnetic fields momentum, but not the ion parallel kinematic momentum or the ion parallel flow velocity. Hence, the turbulent acceleration-the local source or sink of the mean ion parallel flow does not imply that the momentum conservations is broken [16].

A new mechanism for the origin of intrinsic rotation named turbulent acceleration has been discovered. It cannot be written as a divergence of stress term and acts as the local source or sink of parallel rotation. The order of magnitude of turbulent acceleration is comparable to that of the divergence of residual stress which is also thought to be the origin of intrinsic rotation in previous works. Therefore, the turbulent acceleration is significant for intrinsic rotation. The co-current acceleration driven by electrostatic ITG turbulence vanishes after mode transition to CTEM turbulence, which qualitatively agrees with experimental observations. An extension to electromagnetic ITG turbulence has also been investigated, showing that the electromagnetic effects are important for both stress force and turbulent acceleration. We have derived two equivalent equations of momentum conservation and discussed the conserved quantity corresponding to the axisymmetry of tokamak plasmas. We clarified that the turbulent acceleration does not contradict with momentum conservation law.

Our ongoing work is extending to toroidal geometry, such as considering kinetic ballooning mode, which is the dominant electromagnetic mode in the peaked gradient region of DIII-D pedestal plasmas [30]. Moreover, inspired by the investigation of intrinsic rotation (which is mainly carried by ions) driven by turbulence, we are also working on the intrinsic non-inductive current (which is mainly carried by electrons) driven by turbulence. This could be also important for tokamak performance, since the modification of current profile may affect the confinement and a variety of MHD instabilities.

These works are supported by National Natural Science Foundation of China (NSFC) under Contract Nos. 11675059 and 11305071, and the Ministry of Science and Technology of China under Contract No. 2013GB112002.