Derivation

Brief outline of the derivation

Braginskii equations

The derivation starts from Braginskii equations in the limit $\Omega_e\tau_e\gg 1$, where we only keep parallel friction, $$\begin{align} \frac{\partial}{\partial t}n&+\nabla\cdot(n\mathbf{v}_e)=0, \label{braginskii_electroncontinuity} \\ % \frac{\partial}{\partial t}n&+\nabla\cdot(n\mathbf{v}_i)=0, \label{braginskii_ioncontinuity} \\ % m_e n\left[\frac{\partial}{\partial t}+\mathbf{v}_e\cdot\nabla\right]\mathbf{v}_e&=-\nabla p_e-en\left(\mathbf{E}+\frac{1}{c}\mathbf{v}_e\times\mathbf{B}\right)+R_\parallel\mathbf{b}, \label{braginskii_electronmomentum} \\ % M_i n\left[\frac{\partial}{\partial t}+\mathbf{v}_i\cdot\nabla\right]\mathbf{v}_i&=-\nabla p_i -\nabla\cdot \mathbf{P}_i+en\left(\mathbf{E}+\frac{1}{c}\mathbf{v}_i\times\mathbf{B}\right)-R_\parallel\mathbf{b}, \label{braginskii_ionmomentum} \\ % \left[\frac{\partial}{\partial t}+\mathbf{v}_e\cdot\nabla\right]T_e&+\frac{2}{3}T_e\nabla\cdot\mathbf{v}_e=-\frac{2}{3n}\nabla\cdot\mathbf{q}_e+\frac{2}{3n}(Q_e-Q_\Delta) \label{braginskii_electrontemperature} \\ % \left[\frac{\partial}{\partial t}+\mathbf{v}_i\cdot\nabla\right]T_i&+\frac{2}{3}T_i\nabla\cdot\mathbf{v}_i=-\frac{2}{3n}\nabla\cdot\mathbf{q}_i - \frac{2}{3n} \mathbf{P}_i : \nabla\mathbf{v}_i + \frac{2}{3n}Q_\Delta, \label{braginskii_iontemperature} \end{align}$$

with, $$\begin{align} R_\parallel = & en\eta_{\parallel}j_\parallel-0.71n\nabla_\parallel T_e, \\ \mathbf{q}_e = & -0.71\frac{T_e}{e}j_\parallel\mathbf{b} -\frac{5}{2}\frac{cT_en}{eB^2}\mathbf{B}\times\nabla T_e -\chi_\parallel^{e}\nabla_\parallel T_e, \\ \mathbf{q}_i = & \frac{5}{2}\frac{cnT_i}{eB^2}\mathbf{B}\times\nabla T_i-\chi_\parallel^{i}\nabla_\parallel T_i, \\ Q_e =&\eta_\parallel j_\parallel^2-0.71\frac{j_\parallel}{e}\nabla_\parallel T_e, \\ Q_\Delta = &\frac{3m_e}{M_i}\frac{n}{\tau_e}(T_e-T_i) \end{align}$$

A brief discussion of the ion stress tensor $\mathbf{P}_i$ is given in \ref{sec:neocl-visc-corrections}. A more comprehensive discussion on the ion stress tensor $\mathbf{P}_i$ can be found e.g. in the original Braginskii article of 1965 or in Zeiler's habilitation, which is for us also the main reference for the drift-reduction of it $\mathbf{P}_i$. The electron stress tensor is neglected as it will drop out later anyway based on the assumption $\omega\ll\Omega_e$. We note that we assumed quasi-neutrality, i.e. $n:=n_e=Zn_i$, and one of the continuity equations can be dropped in favour of the quasi-neutrality condition $\nabla\cdot\mathbf{j}=\nabla\cdot\left(en\mathbf{v}_i-en\mathbf{v}_e\right)=0$. For charge state $Z>1$, the thermal force coefficient 0.71, resistivity $\eta_\parallel$, electron heat conductivity $\chi_\parallel^e$ and collision time $\tau_e$ have to be adapted.

Drift approximation

Under the drift approximation, the electric field is given as $\mathbf{E}=-\nabla\phi-\frac{1}{c}\frac{\partial}{\partial t}A_\parallel\mathbf{b}$. Multiplying the momentum equations (\ref{braginskii_electronmomentum}) and (\ref{braginskii_ionmomentum}) with $\frac{c}{B^2en}\mathbf{B}\times$ yields for the perpendicular velocities $(\mathbf{v}_{\perp}:=-\mathbf{b}\times\mathbf{b}\times\mathbf{v})$

$$\begin{align} \mathbf{v}_{e\perp}=&\frac{c}{B^2}\mathbf{B}\times\nabla\phi-\frac{c}{enB^2}\mathbf{B}\times\nabla p_e-\frac{1}{\Omega_e}\mathbf{b}\times\frac{d_e}{dt}\mathbf{v}_e, \\ \mathbf{v}_{i\perp}= &\frac{c}{B^2}\mathbf{B}\times\nabla\phi+\frac{c}{enB^2}\mathbf{B}\times\nabla p_i+\frac{1}{\Omega_i}\mathbf{b}\times\frac{d_e}{dt}\mathbf{v}_i+\frac{c}{enB^2}\mathbf{B}\times\nabla\cdot \mathbf{P}_i, \end{align}$$

where we define here the $E\times B$ drift $\mathbf{v}_E:=c/B^2\mathbf{B}\times\nabla\phi$ and the diamagnetic drift $\mathbf{v}_*^{e,i}:=\mp c/(enB^2)\mathbf{B}\times\nabla p_{e,i}$ for the electrons $(e,-)$ and ions $(i,+)$, respectively. We take the limit of low frequencies $d/dt\ll\Omega_i\ll\Omega_e$, where we develop the electron velocity to lowest order and the ion velocity to first order,

$$\begin{align} \mathbf{v}_{e\perp}=\mathbf{v}_E+\mathbf{v}_*^e \\ \mathbf{v}_{i\perp}=\mathbf{v}_E+\mathbf{v}_*^i+\mathbf{u}_{pol}, \end{align}$$ The reason for the latter is that the divergence of the lowest order drifts is comparable to the divergence of the first order drifts.

We obtain the polarisation velocity from inserting the lowest order drift back in as $$\begin{align} \mathbf{u}_{pol}=\frac{1}{\Omega_i}\mathbf{b}\times\frac{d_i}{dt}\left(\mathbf{v}_E+\mathbf{v}_*^i\right)+\frac{1}{\Omega_i}\mathbf{b}\times \left(\mathbf{v}_*^i\cdot\nabla\right)\left(\mathbf{v}_E+\mathbf{v}_*^i\right) + \frac{c}{enB^2}\mathbf{B}\times\left(\nabla\cdot \mathbf{P}_i\right), \label{eq:v_pol} \end{align}$$ where the second (diamagnetic advection) term cancels with a part of the third term (FLR-gyroviscosity), such that the diamagnetic velocity does not appear in the advective derivative, defined as: $$\begin{align} \frac{d_e}{dt}=&\frac{\partial}{\partial t}+\left(\mathbf{v}_E+v_\parallel\mathbf{b}\right)\cdot\nabla, \\ \frac{d_i}{dt}=&\frac{\partial}{\partial t}+\left(\mathbf{v}_E+u_\parallel\mathbf{b}+\mathbf{u}_{pol}\right)\cdot\nabla, \end{align}$$ where in the code the polarisation velocity is neglected in the advective derivative due to practical reasons.

The polarisation velocity will (only) be of importance if it appears under a divergence, $$\begin{align} \nabla\cdot\left(n\mathbf{u}_{pol}\right) = \nabla\cdot\left[-\frac{nM_ic^2}{eB^2}\frac{d_i}{dt}\left(\nabla_\perp\phi+\frac{\nabla_\perp p_i}{en}\right)\right]-\frac{1}{6e}\mathcal{C}(G), \end{align}$$ where we commuted the magnetic field with the advective derivative based on the assumption $L_B^{-1}\gg k_\perp$. We also introduced the ion stress function $$\begin{align} G=-\eta_i\left[\frac{2}{B^{3/2}}\nabla\cdot\left( u_\parallel B^{3/2} \mathbf{b} \right) - \frac{\mathcal{C}(\phi)}{2}-\frac{\mathcal{C}(p_i)}{2en}\right], \label{eq:ion_stress_function} \end{align}$$ with the coefficient $\eta_i=0.96nT_i\tau_i\propto T_i^{5/2}$ and $\tau_i$ the ion collision time. The ion stress tensor and related derivations will be discussed in more detail below.

The curvature operator was introduced as $$\begin{align} \mathcal{C}(f):=-\nabla\cdot\left(\frac{c}{B^2}\mathbf{B}\times\nabla f\right) = -\left(\nabla\times\frac{c\mathbf{B}}{B^2}\right)\cdot\nabla f. \end{align}$$ This simplifies the semantics of the following expressions: $$\begin{align} \nabla\cdot\mathbf{v}_E=&-\mathcal{C}(\phi), \\ \nabla\cdot\left(n\mathbf{v}_*^e\right)=&\frac{1}{e}\mathcal{C}(p_e), \\ \nabla\cdot\left(n\mathbf{v}_*^i\right)=&-\frac{1}{e}\mathcal{C}(p_i). \end{align}$$

Magnetic fluctuations

Additionally to the above drift approximation, we take the small beta limit $$\begin{align} \beta_e = \frac{4\pi p_e}{B^2} = \frac{c_s^2}{v_A^2} \ll 1, \end{align}$$ where $v_A = B / \sqrt{4\pi n_i M_i}$ is the Alfven speed. Note that $\beta_e$ is the dynamic electron beta, and not the full plasma beta. This limit allows to neglect magnetic field disturbances $\mathbf{\tilde{B}}_\perp$ in the drift operator $(c/B^2)\mathbf{B}\times$, since $|\mathbf{\tilde{B}}_\perp|\ll|\mathbf{B}|$. This is fine for the perpendicular motion, but the parallel motion is fast, and therefore it can contribute to perpendicular transport through $\mathbf{\tilde{B}}_\perp$ nonetheless.

To take this into account, the perturbed magentic field is computed as $\mathbf{\tilde{B}}_\perp=\nabla\times A_1\mathbf{b}_\mathrm{tor}$ and the unit vector as $\mathbf{\tilde{b}}_\perp=\mathbf{\tilde{B}}_\perp / B_\mathrm{tor}$, in which $A_1$ is defined as the fluctuating part of $A_\parallel$ (see magnetic shift). We treat $\mathbf{\tilde{B}}_\perp$ carefully to ensure the exact divergence-free $\nabla \cdot \mathbf{\tilde{B}}_\perp = 0$. $$\begin{align} &\mathbf{\tilde{B}}_\perp=\nabla\times A_1\mathbf{b}_\mathrm{tor}=-\mathbf{B}_\mathrm{tor} \times \nabla \frac{A_1}{B_\mathrm{tor}} + \frac{A_1}{B_\mathrm{tor}} \nabla \times \mathbf{B}_\mathrm{tor}, \\ &=-\mathbf{B}_\mathrm{tor} \times \nabla \frac{A_1}{B_\mathrm{tor}}. \label{magnetic_fluctations} \end{align}$$ The reduction will be exact as long as the magnetic equilibrium has $\nabla \times \mathbf{B}_\mathrm{tor}=0$, namely the vacuum field assumption.

The parallel gradient is then defined as $$\begin{align} &\nabla_\parallel f = \nabla_\parallel^{(0)} f + \nabla_\parallel^{(1)} f, \\ &\nabla_\parallel^{(1)} f = - \mathbf{b}_\mathrm{tor} \cdot \left( \nabla \frac{A_1}{B_\mathrm{tor}} \times \nabla f \right), \\ \label{eq:gradpar2} \end{align}$$ The zero order term $\nabla_\parallel^{(0)} f$ is along the equilibrium field. $\nabla_\parallel^{(1)} f$ can be assumed to be the leading order correction due to the magnetic perturbation (Zeiler's habil. eq. (2.78), Scott's habil. eq. (14.7)).

For the divergence, the result is very similar because the perturbed magnetic field is also divergence free, $$\begin{align} \nabla \cdot \mathbf{b}_\perp f = \nabla \cdot \mathbf{\tilde{B}}_\perp \frac{f}{B_\mathrm{tor}} = \mathbf{\tilde{B}}_\perp \cdot \nabla \frac{f}{B_\mathrm{tor}} = B_\mathrm{tor} \nabla_\parallel^{(1)} \frac{f}{B_\mathrm{tor}} \end{align}$$

In the following, we will write $\mathbf{b}$ instead of $\mathbf{b}_\mathrm{tot}$ for the total magnetic field unity vector, and for the equilibrium part $\mathbf{b}_0$ instead of $\mathbf{b}$.

Gyroviscous stress tensor

We follow here the discussion of Zeiler's habilitation. The ion viscous stress tensor consists of a collisionless FLR part: $$\begin{align*} \nabla\cdot\mathbf{P}_i^{FLR} = -M_in\mathbf{v}_*^i\cdot\nabla\mathbf{v}_i + p_i\nabla\times\frac{\mathbf{b}}{\Omega_i}\cdot\nabla\mathbf{v}_i + \nabla_\perp\left(\frac{p_i}{2\Omega_i}\nabla\cdot\mathbf{b}\times\mathbf{v}_i\right) + \mathbf{b}\times\nabla\left(\frac{p_i}{2\Omega_i}\nabla_\perp\cdot\mathbf{v}_i\right), \end{align*}$$ where the first term cancels the diamagnetic advection at the computation of the polarisation drift in eq. (\ref{eq:v_pol}) and lateron also in the parallel momentum equation. On the other terms there are partly some subtle approximations applied, for which we refer to the discussion in Zeilers habilitation.

The viscous part of the stress tensor is: $$\begin{align*} \mathbf{P}_i^{vis} = \left(\mathbf{b}\mathbf{b}-\mathbf{I}/3\right)G, \end{align*}$$ where for the ion stress function Zeiler concludes with the following approximate expression $$\begin{align*} G = -\eta_i\left(2\nabla_\parallel u_\parallel -\mathbf{\kappa}\cdot\mathbf{v}_{\perp,i0}\right), \end{align*}$$ where here and in all other expressions for the gyroviscous stress tensor the zeroth order velocities, i.e. $\mathbf{v}_{\perp,i0} = \mathbf{v}_E+\mathbf{v}_*^i + u_\parallel\mathbf{b}$ is used, which ensures $P_i:\nabla\mathbf{v}=0$ for whichever vector $\mathbf{u}$ is used to construct $\mathbf{P}_i$ (see Scott's habilitation). In GRILLIX the slightly modified expression (\ref{eq:ion_stress_function}) is used (possibly with (see transcollisional modifications). Again the divergence of the viscous ion stress tensor is of relevance: $$\begin{align*} \nabla\cdot\mathbf{P}_i^{vis} = G\mathbf{\kappa} - \nabla G/3 +\mathbf{B}\nabla_\parallel\frac{G}{B}. \end{align*}$$

Finally, the following expression is of importance for computations related with the energy theorem: $$\begin{align} \nabla\cdot\left(fn\mathbf{u}_{pol}^{vis}\right) := \nabla\cdot\left(\frac{cf}{eB^2}\mathbf{B}\times\left(\nabla\cdot\mathbf{P}_i^{vis}\right)_\perp\right) = f\nabla\cdot\left(\frac{c}{eB^2}\mathbf{B}\times\left(\nabla\cdot\mathbf{P}_i^{vis}\right)_\perp\right) + \nabla f \cdot\left(\frac{c}{eB^2}\mathbf{B}\times\left(\nabla\cdot\mathbf{P}_i^{vis}\right)_\perp\right)., \end{align}$$ where $f$ is some arbitrary field. Then the first term can be can be evaluated using the result for the polarisation current derived e.g. in Zeiler's habilitation and the second term can be reformulated into a drift-advection like term and a curvature term acting on $G$: $$\begin{align} \nabla\cdot\left(fn\mathbf{u}_{pol}^{vis}\right) = -\frac{f}{6e}\mathcal{C}(G) + \frac{1}{3e}\left(\frac{c\mathbf{B}}{B^2}\times\nabla f\right)\cdot\nabla G - \frac{G}{2e}\mathcal{C}(f) \end{align}$$