Turbulent transport and heating of trace heavy ions in hot, magnetized plasmas

Heavy ions are present in hot, magnetized plasmas both in laboratory experiments and in nature. These heavy ions are often trace, i.e., their densities are small enough that they have only a small direct effect on the bulk plasma dynamics. Nonetheless, trace heavy ions are important in numerous contexts: main ion properties are often inferred from heavy ion measurements because heavy ions radiate more readily islerPPCF94 ; accumulation of heavy ions leads to dilution and increased radiative energy losses in magnetic confinement fusion meserveyNF76 ; tokarNF97 ; and temperature measurements of minority ions in space and astrophysical plasmas indicate the existence of a novel heating mechanism schmidtGRL80 ; collierGRL96 ; kohlSP97 .

Considerable effort has gone into understanding the particle transport of trace heavy ions, or impurities, in the context of magnetized toroidal plasmas for fusion. In particular, the scaling with charge number $Z$ and mass number $A$ of the impurity particle flux were predicted with a quasilinear fluid model and found to be in relatively good agreement with numerical and experimental results angioniPRL06 ; angioniPPCF09 . However, little to no work has been done on impurity momentum and energy fluxes or for turbulent heating of impurities. The latter may play a role not only in fusion plasmas, but also in the context of astrophysical plasmas, where the temperature of minority ions has been observed to increase with increasing ion mass schmidtGRL80 ; collierGRL96 ; kohlSP97 . Cyclotron heating cranmerApJ99 and stochastic heating via large-amplitude fluctuations chandranApJ10a have been proposed as possible explanations for this mass dependence. The turbulent heating mechanism described here provides an alternative explanation for the mass dependence of the minority ion temperature that is present even for low frequency, low amplitude fluctuations.

In this Letter we use local, nonlinear, $\delta f$-gyrokinetic theory cattoPP78 ; friemanPoF82 ; howesApJ06 to provide scaling predictions for trace heavy ion particle, momentum, and energy fluxes, as well as turbulent heating in hot, magnetized plasmas. This approach has already proven successful in determining scalings of temperature-gradient driven turbulence in tokamaks barnesPRL11b . We consider an inhomogeneous, axisymmetric plasma rotating toroidally at angular frequency $\omega_{\phi}$, immersed in a curved, inhomogeneous magnetic field. To simplify our analysis, we restrict our attention to a region of plasma with rotation speed well below the ion sound speed but with a strong rotation gradient. We also consider only moderate values of $\beta=8\pi p/B^{2}\lesssim 1$, where $p$ is the mean plasma pressure and $B$ is the mean magnetic field magnitude. This is directly applicable to toroidal confinement experiments in magnetic confinement fusion, but the scaling laws we obtain are general: they do not change for a stationary, homogeneous plasma slab and therefore also pertain to various space and astrophysical plasmas.

The $\delta f$-gyrokinetic theory is obtained by performing an asymptotic expansion in the small ratio of the Larmor radius, $\rho$, to system size, $L$, and averaging over the fast Larmor motion of particles. It is valid for low-amplitude turbulence with time scales long compared to the Larmor frequency, $\Omega$, and spatial scales comparable to $\rho$ and $L$ in the directions across and along the mean magnetic field, respectively. While initially developed for magnetic confinement fusion plasmas, $\delta f$-gyrokinetics can also be applied to small-scale turbulence in the solar wind, solar corona, accretion disks, and galaxy clusters  howesJGR08 ; schekApJ09 .

We use ($\mathbf{R},\mu,\varepsilon$) as our coordinate system, where $\mathbf{R}$ is the position of the center of a particle’s Larmor orbit, $\varepsilon=mv^{2}/2$ its kinetic energy, and $\mu=mv_{\perp}^{2}/2B$ its magnetic moment, with $m$ its mass and $v$ its speed. The subscripts $\perp$ and $\parallel$ are used to denote the components perpendicular and parallel to the mean magnetic field, respectively, with the magnetic field magnitude given by $B$. With this choice of coordinates, the electromagnetic gyrokinetic equation governing the evolution of the fluctuating piece of the distribution function, $\delta f_{s}$, is

where $g_{s}=\delta f_{s}+Z_{s}eF_{M,s}(\Phi-{\left<\chi\right>}_{s})/T_{s}$, ${\left<.\right>}_{s}$ denotes an average over Larmor angle at fixed $\mathbf{R}_{s}$, ${\left<\chi\right>}_{s}={\left<\Phi-v_{\parallel}\delta A_{\parallel}/c+\int_{0}^{\mu_{s}}d\mu_{s}^{\prime}\delta B_{\parallel}/Z_{s}e\right>}_{s}$ ¹¹1The fields $\Phi$, $\delta A_{\parallel}$, and $\delta B_{\parallel}$ are independent of Larmor angle at fixed particle position, $\mathbf{r}$, but not at fixed $\mathbf{R}=\mathbf{r}+\mathbf{v}_{\perp}\times\mathbf{\hat{b}}/\Omega$. Thus care must be taken to specify which spatial coordinate is held fixed for velocity integration. The $\mu$-integral contained in ${\left<\chi\right>}$ is performed at fixed $\mathbf{R}$, but all other velocity integrals in this Letter are performed at fixed $\mathbf{r}$., $\Phi$ is the fluctuating electrostatic potential, $\delta A_{\parallel}$ and $\delta B_{\parallel}$ are the parallel components of the fluctuating magnetic vector potential and magnetic field, respectively, $Z_{s}$ is the charge number, $e$ the proton charge, $c$ the speed of light, $T_{s}$ the mean temperature, $F_{M,s}$ is a stationary Maxwellian distribution of velocities in the frame rotating with velocity $\mathbf{u}=R^{2}\omega_{\phi}\nabla\phi$, $\phi$ is the toroidal angle, $R$ the plasma major radius, $D/Dt=\partial/\partial t+\mathbf{u}\cdot\nabla$, $\dot{\mathbf{R}}_{s}=\mathbf{v}_{\parallel}+\mathbf{v}_{M,s}+{\left<\mathbf{v}_{\chi}\right>}_{s}$, with $\mathbf{v}_{M,s}=\mathbf{\hat{b}}/\Omega_{s}\times(v_{\parallel}^{2}\mathbf{\hat{b}}\cdot\nabla\mathbf{\hat{b}}+v_{\perp}^{2}\nabla B/2B)$ the drift velocity due to a curved, inhomogeneous mean magnetic field and ${\left<\mathbf{v}_{\chi}\right>}_{s}=c\mathbf{\hat{b}}\times\nabla{\left<\chi\right>}_{s}/B$ the drift due to the fluctuating electromagnetic fields, $\mathbf{\hat{b}}$ the unit vector along the mean magnetic field, $\Omega_{s}=Z_{s}eB/m_{s}c$ the Larmor frequency, and $C$ describes two-particle Coulomb interactions. Plasma species is indicated by the subscript $s$, which we henceforth drop unless it is needed to avoid ambiguity.

By definition, the trace ions considered here do not contribute to the fields. They are instead determined solely by the electron and main ion dynamics through the low-frequency Maxwell’s equations, supplemented by the quasineutrality constraint:

where $\Phi$ enters Eqs. (2-4) through the definition for $\delta f$ given below Eq. (1).

With $g$ and $\{\Phi,\delta A_{\parallel},\delta B_{\parallel}\}$ specified by Eqs. (1-4), one can evaluate the turbulent heating,

and the turbulent fluxes,

where $r$ labels surfaces of constant mean pressure, $\left<a\right>_{\Lambda}=\int d^{3}\mathbf{r}\int d^{3}\mathbf{v}\ a/\int d^{3}\mathbf{r}$ is an integral over all velocity space and over a volume of width $w$ ($\rho\ll w\ll L$) encompassing the mean magnetic field line of interest, and $\Gamma$, $\Pi$, and $Q$ are the particle, toroidal angular momentum, and energy fluxes, respectively. Note that the momentum flux defined in Eq. (8) does not include each species’ contribution to the Maxwell stress.

To obtain scaling laws for the turbulent fluxes and heating of trace heavy ions, we take $Z\sim A\gg 1$, $d\omega_{\phi}/dr\sim v_{ti}/R^{2}$, and expand $g=g_{0}+g_{1}+...$ in powers of $A^{1/2}$. Here $v_{ti}$ is the main ion thermal speed. We restrict our attention to $\beta=8\pi p/B^{2}\lesssim 1$, and assume the collisional mean free path is sufficiently long that collisions may be neglected in our analysis. In what follows, we keep $Z$ and $A$ dependences separate so that we can consider the subsidiary expansion $A^{1/2}\ll Z\ll A$.

Because the heavy ions are trace, their space and time scales are those of the bulk plasma turbulence. Thus, $Z$ and $A$ only enter Eq. (1) through explicit factors of $m$, $v\sim v_{t}$, and $Z$, as well as through $g$ itself. In what follows, we assume the ratio of the heavy ion to proton temperature is much smaller than $A$, giving $v_{t}\sim A^{-1/2}$. The two lowest order equations in our expansion are thus

where $\mathbf{v}_{E}=c\mathbf{\hat{b}}\times\nabla\Phi/B$ and $\mathbf{v}_{A}=v_{\parallel}\mathbf{\hat{b}}\times\nabla\delta A_{\parallel}/B$.

There are two possible scalings for both $g_{0}$ and $g_{1}$ due to a competition between terms with different $A$ and $Z$ dependences in Eqs. (9) and (10). In particular, $g_{0}\propto A^{1/2},\ Z/A^{1/2}$, and $g_{1}\propto 1,\ Z/A$. By considering the limits $|d\omega_{\phi}/dr|\ll(Z/A)v_{ti}/R^{2}$ and $|d\omega_{\phi}/dr|\gg(Z/A)v_{ti}/R^{2}$, the number of such scalings is reduced. The $A$- and $Z$-scalings for $g_{0}$ and $g_{1}$ in these limits, as well as for the general case, are summarized in Table 1.

If $g_{0}(v_{\parallel})$ is a solution to Eq. (9), then $-g_{0}(-v_{\parallel})$ is also a solution. Thus, $\overline{\int_{-\infty}^{\infty}dv_{\parallel}g_{0}\{\Phi,\delta A_{\parallel},\delta B_{\parallel}\}}=0$, where the overline denotes a statistical average. As a result, $g_{0}$ does not contribute to the lowest order (i.e., electrostatic) heating or particle and heat fluxes, Eqs. (5)-(7), whose integrands are otherwise even functions of $v_{\parallel}$. Conversely, the lowest order momentum flux integrand has a component proportional to $mv_{\parallel}$, so $\Pi\sim mv_{t}g_{0}\propto A^{1/2}g_{0}$. Using our scalings for $g_{0}$, we see that $\Pi$ has competing terms scaling as $Z$ and $A$, respectively.

Note that Eq. (10) has a $v_{\parallel}$ symmetry that is opposite that of Eq. (9): if $g_{1}(v_{\parallel})$ is a solution, then $g_{1}(-v_{\parallel})$ is also a solution. Furthermore, for all higher order equations, one can show that the symmetry in $v_{\parallel}$ alternates between that of Eqs. (9) and (10). As a result, the only components of $g$ that contribute to the particle and heat fluxes and heating are $g_{1}$, $g_{3}$, etc. Using Eqs. (6) and (7), we have $\{\Gamma,Q\}\sim g_{1}$, which in the general case has competing terms scaling as $Z/A$ and $1$ (no $Z$ or $A$ dependence), respectively.

The first term in the heating expression, (5), is the Joule heating and is scaled up by an explicit factor of $Z$ (arising from the current), while the second term is viscous heating. At lowest order, the Joule heating term gives $H\propto(Zv_{\parallel}g_{0},\ Zg_{1})$, giving $H\propto(Z,\ Z^{2}/A)$. The viscous heating is proportional to $\Pi\propto(Z,\ A)$. $H$ thus has competing terms scaling as $Z^{2}/A$, $Z$, and $A$, respectively. The scalings of the various fluxes and heating are summarized in Table 1.

Integrating Eq. (5) by parts in time and using Eq. (1), the heating can be expressed as krommesPoP94 ; howesApJ06 ; abelRPP12

Our scalings indicate that $H_{s}$ increases in magnitude with $A$ or $Z$, but $\Gamma_{s}$ and $Q_{s}$ do not. The first term in Eq. (11) must thus dominate for $A$ or $Z$ large. This term, which we identify as the collisional entropy generation, is positive definite when summed over species. We argue that it is also positive species by species for the low collisionalities considered here.

The collision operator, $C$, consists of a test-particle piece, which is a diffusion operator in velocity space, and a field-particle piece, which is an integral operator helander . Both contributions are inversely proportional to the collisional mean free path and thus small, except at small scales in the velocity space where large derivatives in the test-particle operator compensate abelPoP08 ; barnesPoP09 ; schekApJ09 . The test-particle operator should thus dominate in weakly collisional plasmas, and its diffusive nature ensures that its contribution to entropy generation is positive-definite.

Consequently, trace heavy ions must be heated by turbulence instead of cooled. For this heating process to subside, the trace ion temperature must become large enough to interfere with our large $A$ expansion. This happens when the heavy ion temperature exceeds the main ion temperature by a factor of $A\sim Z$. In this limit, the turbulent heating $H$ becomes comparable to the heat flux $Q$ so that $H$ is no longer required to be positive definite. Our theory thus predicts that heavy ions will be hotter than light ions by a factor of $A\sim Z$ – but only if turbulent heating is larger than collisional temperature equilibration.

The collisional temperature equilibration of the main ions, $i$, and a trace heavy ion species, $s$, is $\mathcal{E}_{s}\equiv(8/3\sqrt{\pi})(Z_{s}^{2}/A_{s})n_{s}\Delta T_{s}v_{ti}/\lambda_{\textnormal{mfp}}$, where $\Delta T_{s}=T_{s}-T_{i}$, and $\lambda_{\textnormal{mfp}}$ is the mean free path for collisions between the main ions. From Eq. (5), we estimate $H_{s}\sim Sn_{s}T_{i}(\delta n_{i}/n_{i})^{2}v_{ti}/L$, where $S$ is the scaling of $H$ with $A$ and $Z$ given in Table 1, and we have assumed $e\Phi/T_{i}\sim\delta n_{i}/n_{i}\sim\delta n_{s}/n_{s}$. The ratio of turbulent heating to collisional temperature equilibration is thus $H/\mathcal{E}\sim S(A/Z^{2})(\lambda_{\textnormal{mfp}}/L)(\delta n_{i}/n_{i})^{2}\sim S(A/Z^{2})\lambda_{\textnormal{mfp}}/v_{ti}\tau_{E}$, with $\tau_{E}$ the characteristic time scale over which the equilibrium density and temperature vary.

To test our predictions for the scalings of turbulent transport and heating, we employ the local, $\delta f$-gyrokinetic code GS2 dorlandPRL00 . We consider an axisymmetric system with sheared magnetic field lines mapping out nested toroidal surfaces with circular cross sections (known as the Cyclone Base Case dimitsPoP00 and parametrized using the Miller local equilibrium model millerPoP98 ). Each simulation is electrostatic and includes kinetic electrons, as well as kinetic main and trace heavy ions with a wide range of $Z$ and $A$ values. The turbulence is driven by gradients in the mean ion and electron densities and temperatures, with $R_{0}(d\ln n/dr)=2.2$ for the electrons and main ions, and $R_{0}(d\ln T/dr)=6.9$ for all species, with $R_{0}$ the major radius at the center of the constant pressure surface. The collision frequency is chosen small, $R_{0}/\lambda_{\textnormal{mfp}}=0.003$, so that heavy ion collisions do not affect our scalings.

Two sets of simulations were carried out: one with a stationary plasma ($d\omega_{\phi}/dr=0$) and one with a differentially rotating plasma ($d\omega_{\phi}/dr=4.67v_{ti}/R^{2}$). The simulation results are shown in Figs. (1)-(4). Data points for fluxes and heating at various $Z$ and $A$ values are plotted as solid circles and fit using a least-squares analysis with the predicted lowest order $Z$ and $A$ dependences, as well as the first order correction. In each case, the predicted scalings fit the data well. It should be noted that the momentum flux for $d\omega_{\phi}/dr=0$ is zero for all species due to a fundamental symmetry of the $\delta f$-gyrokinetic equation parraPoP11 .

We now discuss the implications of the trace heavy ion scalings derived in this Letter. First, the preferential heating of heavy ions should lead to large temperature disparities between different ion species in nearly collisionless plasmas. Many space and astrophysical plasmas are weakly collisional enough (i.e., $(\delta n/n)^{2}>L/\lambda_{\textnormal{mfp}}$) that turbulent heating should dominate over collisional equilibration, and preferential heating of heavy ions is indeed observed schmidtGRL80 ; collierGRL96 . However, for such low collisionalities the equilibrium can deviate strongly from the isotropic Maxwellian assumed in our analysis, which cannot consequently address the large $T_{\perp}/T_{\parallel}$ values observed in coronal holes and the fast solar wind kohlSP97 .

Magnetic confinement fusion plasmas typically have $(\delta n/n)^{2}<L/\lambda_{\textnormal{mfp}}$ so that collisional temperature equilibration dominates over turbulent heating and all ions have the same temperature. However, for rotating plasmas our results indicate that the turbulent heating is enhanced by an additional factor of $(A/Z)^{2}$ relative to the equilibration. It may therefore be possible for heavy, partially ionized impurities to be heated by turbulence to temperatures significantly larger than the main ions.

Because the momentum transport of heavy ions is enhanced by $A$, it can generate flows of order the ion thermal speed for densities as small as $n_{i}/A$. In this limit, heavy ions could thus significantly alter bulk plasma momentum transport.

We thank S. C. Cowley, E. Quataert, and A. A. Schekochihin for useful discussions. M.B. was supported by a US DoE FES Postdoctoral Fellowship, F.I.P. was supported by US DoE Grant No DE-FG02-91ER-54109, and computing time was provided by HPC-FF (Jülich).