2. Kinetic diffusion equation

The starting point is the drift-kinetic equation for the gyro-averaged distribution function of particles of species j (Hinton & Hazeltine Reference Hinton and Hazeltine1976):

(2.1)\begin{equation}\partial {f_j}/\partial t + {v_\parallel }\boldsymbol{b}\boldsymbol{\cdot }\boldsymbol{\nabla }{f_j} + {Z_j}e {v_\parallel }\boldsymbol{b}\boldsymbol{\cdot }{\boldsymbol{E}_{\textrm{in}}}\partial {f_j}/\partial \varepsilon = {C_j},\end{equation}

with ${f_j} = {f_j}(r,\theta ,\phi ,\varepsilon ,\mu )$, $\varepsilon = {\textstyle{1 \over 2}}{m_j}{v^2} + {Z_j}e\varPhi$ the total energy, ${m_j}$ the mass, ${v^2} = {v_\parallel }^2 + v_ \bot ^2$, ${v_\parallel }, {v_ \bot }$ the velocities of the parallel and gyration motion respectively, ${Z_j}e$ the electric charge, $e > 0$ the elementary charge, $\varPhi$ the electrostatic potential, $\mu = {m_j}v_ \bot ^2/(2B)$ the magnetic moment, $\boldsymbol{b} = \boldsymbol{B}/B$ the magnetic unit vector, ${\boldsymbol{E}_{\textrm{in}}}$ the externally applied inductive electric field and ${C_j}$ the collision term. From now on we drop the subscript j, if not strictly necessary. As in Harvey et al. (Reference Harvey, McCoy, Hsu and Mirin1981) and Rice et al. (Reference Rice, Molvig and Helava1982), we neglect the drift velocity perpendicular to the magnetic field and we assume both the magnetic field (including its perturbation) and the electrostatic potential constant in time $(\partial B/\partial t = \partial \varPhi /\partial t = 0)$. The plasma here considered has destroyed magnetic surfaces due to the overlapping of several resonant magnetic perturbations. The total magnetic field is written as the sum of the equilibrium ${\boldsymbol{B}_0}$ and the perturbation ${\boldsymbol{B}_1}$:

(2.2a–c)\begin{equation}\boldsymbol{B} = {\boldsymbol{B}_0}(r) + {\boldsymbol{B}_1}(r,\theta ,\phi ),\quad {B_{0,r}} = 0,\quad \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{B} = \boldsymbol{\nabla }\boldsymbol{\cdot }{\boldsymbol{B}_0} = \boldsymbol{\nabla }\boldsymbol{\cdot }{\boldsymbol{B}_1} = 0.\end{equation}

We take the ratio $\delta = {B_1}/{B_0}$ to be a small parameter. It is also useful to split the magnetic field unit vector $\boldsymbol{b}$ into the equilibrium ${\boldsymbol{b}_0} = {\boldsymbol{B}_0}/{B_0}$ and perturbation components:

(2.3a,b)\begin{equation}\boldsymbol{b} = {\boldsymbol{b}_0} + \tilde{\boldsymbol{b}},\quad \tilde{\boldsymbol{b}} = {\boldsymbol{B}_1}/{B_0} - ({\boldsymbol{B}_0}\boldsymbol{\cdot }{\boldsymbol{B}_1}/B_0^2) {\boldsymbol{b}_0} + o({\delta ^2}).\end{equation}

In the presence of a stochastic magnetic field, the parallel gradient ${v_\parallel }\boldsymbol{b}\boldsymbol{\cdot }\boldsymbol{\nabla }f$ of (2.1) gives rise to a diffusive term for the distribution function, stemming from quasi-linear approximation (Harvey et al. Reference Harvey, McCoy, Hsu and Mirin1981; Rice et al. Reference Rice, Molvig and Helava1982). To prove it, let us write the distribution function as

(2.4)\begin{equation}f = \langle\, f\rangle (r) + \tilde{f}(r,\theta ,\phi ),\end{equation}

where $\langle \rangle$ denotes the average over the angles $\theta ,\phi$ and $\tilde{f}\sim o(\delta )\langle\, f\rangle$ is the perturbation component. Then, we approximate (2.1) by the following two equations:

(2.5)\begin{gather}\partial \tilde{f}/\partial t + {v_\parallel } {\boldsymbol{b}_0}\boldsymbol{\cdot }\boldsymbol{\nabla }\tilde{f} + {v_\parallel } \tilde{\boldsymbol{b}}\boldsymbol{\cdot }\boldsymbol{\nabla }\langle\, f\rangle \cong 0,\end{gather}

(2.6)\begin{gather}\partial \langle\, f\rangle /\partial t + {v_\parallel }\langle \tilde{\boldsymbol{b}}\boldsymbol{\cdot }\boldsymbol{\nabla }\tilde{f}\rangle + Ze{v_\parallel }{\boldsymbol{b}_0}\boldsymbol{\cdot }{\boldsymbol{E}_{\textrm{in}}}\partial \langle\, f\rangle /\partial \varepsilon \cong C(\langle\, f\rangle ).\end{gather}

All the terms in (2.5) are of order $({v_\parallel }/a) o(\delta )\langle\, f\rangle$. All the terms in (2.6) are assumed to be of higher order: $({v_\parallel }/a)o({\delta ^2})\langle\, f\rangle$. As far as $\partial \langle\, f\rangle /\partial t$ and $C(\langle\, f\rangle )$ are concerned, this is justified in conditions close to an equilibrium described by a Maxwellian distribution. As far as the ${\boldsymbol{E}_{\textrm{in}}}$ term is concerned, the underlying idea is that the electric field balances the diffusion ${v_\parallel }\langle \tilde{\boldsymbol{b}}\boldsymbol{\cdot }\boldsymbol{\nabla }\tilde{f}\rangle$ and the collisions. Equation (2.5) is time-integrated along the characteristics, given by the field line of the equilibrium magnetic field ${\boldsymbol{b}_0}$. Since collisions are not considered in (2.5), this case corresponds to the collision-less solution of $\tilde{f}$. The collisional correction is introduced heuristically in a second step. By inserting the $\tilde{f}$ solution into the term ${v_\parallel }\langle \tilde{\boldsymbol{b}}\boldsymbol{\cdot }\boldsymbol{\nabla }\tilde{f}\rangle$ of (2.6), one gets the following kinetic diffusion equation for $\langle\, f\rangle$:

(2.7)\begin{align}\left. \begin{array}{*{20}{c}} {\partial \langle\, f\rangle /\partial t \cong |{v_\parallel }|\dfrac{1}{r}\mathcal{L}[r {D_0}\mathcal{L} \langle\, f\rangle ] + |{v_\parallel }|L_0^{ - 1} {D_0}\mathcal{L} \langle\, f\rangle - Ze{v_\parallel }{\boldsymbol{b}_0}\boldsymbol{\cdot }{\boldsymbol{E}_{\textrm{in}}}\partial \langle\, f\rangle /\partial \varepsilon + C(\langle\, f\rangle ),}\\ \mathcal{L} = {{\left( {\dfrac{\partial }{{\partial r}}} \right)}_{\varepsilon ,\mu }},\quad {L_0} = {B_0}/(\textrm{d}{B_0}/\textrm{d}r). \end{array} \right\}\end{align}

Details of the derivation are given in Appendix A. The first two terms on the right-hand side stem from ${v_\parallel }\langle \tilde{\boldsymbol{b}}\boldsymbol{\cdot }\boldsymbol{\nabla }\tilde{f}\rangle$, and they describe radial diffusion processes due to the stochastic magnetic field. The coefficient ${D_0}$ is the magnetic diffusivity, which relates the average squared radial displacement $\langle {(\mathrm{\Delta }r)^2}\rangle$ of the field lines to the distance l covered along the lines: $\langle {(\mathrm{\Delta }r)^2}\rangle = 2{D_0}l$ (Rechester & Rosenbluth Reference Rechester and Rosenbluth1978; D'Angelo & Paccagnella Reference D'Angelo and Paccagnella1996). By defining the Fourier expansion $\tilde{\boldsymbol{b}} = \sum\nolimits_{\boldsymbol{m},\boldsymbol{n} \in \mathrm{\mathbb{Z}}} {{{\tilde{\boldsymbol{b}}}^{\boldsymbol{m},\boldsymbol{n}}}\,{\textrm{e}^{\textrm{i}(m\theta - n\phi )}}}$, ${D_0}$ is expressed in terms of the harmonics $\tilde{b}_r^{m,n}$:

(2.8a,b)\begin{equation}{D_0}(r) = \mathop \sum \limits_{\begin{array}{*{20}{c}} {m \in ,\mathrm{\mathbb{Z}}}\\ {n > 0} \end{array}} 2\mathrm{\pi }\delta (m - nq)|\tilde{b}_r^{m,n}{|^2}\left|{\frac{{R\, q}}{{{b_{0,\phi }}}}} \right|,\quad q(r) = r{B_{0,\phi }}/(R{B_{0,\theta }}).\end{equation}

Only the resonant harmonics, i.e. those satisfying the condition $q({r^{m,n}}) = m/n$ inside the plasma $({r^{m,n}} < a)$, give a contribution to ${D_0}$. Note that the sum runs over half of the harmonics: as shown in Appendix A, the delta-function includes the contribution of the complex-conjugate harmonics (i.e. the complementary half $m \in \mathrm{\mathbb{Z}},\, n < 0$). The $\delta$ function in (2.8) derives from integrating (2.5) along the unperturbed magnetic field line, and it is a typical feature of the quasi-linear theory. Such a singular function can be smoothed by exploiting nonlinear resonance broadening techniques (Ishihara, Xia & Hirose Reference Ishihara, Xia and Hirose1992), which instead consider the effect of the perturbed orbits. The development of such a method for our specific case is beyond the scope of the present paper, and it is left as future refinement of the work.

The term $|{v_\parallel }|L_0^{ - 1}{D_0}\mathcal{L}\langle\, f\rangle$ of (2.7) is not present in previous derivations (Harvey et al. Reference Harvey, McCoy, Hsu and Mirin1981; Rice et al. Reference Rice, Molvig and Helava1982), where the tokamak ordering there assumed makes it negligible: $L_0^{ - 1} = {a^{ - 1}}o({\epsilon ^2}),\epsilon = a/R$, with $\epsilon$ small. Nonetheless, this term cannot be discarded a priori, since $L_0^{ - 1} \sim {a^{ - 1}}$ for a general pinch.

Expression (2.8) corresponds to the collision-less case, since it derives from the $\tilde{f}$ obtained in this limit. Collisions hamper the free motion of the particles along the field lines, thereby reducing the stochastic transport. Following Predebon & Paccagnella (Reference Predebon and Paccagnella2007), we include their effect by an ad hoc factor in front of ${D_0}$:

(2.9a–d)\begin{equation}{D_j} = {g_j}({v_\parallel }){D_0},\quad {g_j}({v_\parallel }) = \frac{1}{{1 + {L_K}/{\lambda _j}({v_\parallel })}},\quad {\lambda _j}({v_\parallel }) = \frac{{{v_\parallel }}}{{{\nu _{0j}}}},\quad {\nu _{0j}} = \frac{{Z_j^2{e^4}\ln (\varLambda ) {n_i}}}{{4\mathrm{\pi }\varepsilon _0^2m_j^2v_\parallel ^3}},\end{equation}

where ${\lambda _j}({v_\parallel })$ is the mean free path of a test particle j in a hydrogen background Maxwellian plasma whose main ion density is ${n_i}$, computed by the fundamental collision frequency ${\nu _{0j}}$ (Trubnikov Reference Trubnikov1965) (the actual frequencies of the various relaxation processes involve extra coefficients); $\ln (\varLambda )$ is the Coulomb logarithm; and ${L_K}$ is the Kolmogorov–Lyapunov length (Rechester, Rosenbluth & White Reference Rechester, Rosenbluth and White1979; D'Angelo & Paccagnella Reference D'Angelo and Paccagnella1996), which is the characteristic length for the exponential divergence of close magnetic field lines. The collision-less limit ${L_K}/{\lambda _j} \to 0$ corresponds to (2.8). In the opposite limit ${L_K}/{\lambda _j} \to + \infty$ (as occurs for particles with ${v_\parallel } \to 0$), collisions are so frequent that particles do not experience the stochastic nature of the magnetic field: the vanishing of expression (2.9) for D in this limit means that other transport mechanisms take over the stochastic one. Actually, ${L_K}$ should be corrected by a logarithmic factor involving further characteristic lengths (Rechester & Rosenbluth Reference Rechester and Rosenbluth1978; Krommes, Oberman & Kleva Reference Krommes, Oberman and Kleva1983; Predebon & Paccagnella Reference Predebon and Paccagnella2007): these details, not crucial for our purposes, are here omitted. Hereafter, we replace ${D_0}$ with D in (2.7).

A solution of (2.7) for $\langle\, f\rangle$ is not attempted here. Instead, we take velocity moments of (2.7), and then we convert them into transport equations for macroscopic plasma quantities by approximating $\langle\, f\rangle$ with a local shifted Maxwellian.

3. Velocity moments of the kinetic diffusion equation

In order to take the velocity moments of (2.7), we perform the variable change $(\varepsilon ,\mu ) \to ({v_\parallel },{v_ \bot })$ by

(3.1a,b)\begin{equation}{v_\parallel } ={\pm} \sqrt {\frac{2}{m}(\varepsilon - \mu B - Ze\varPhi )} ,\quad {v_ \bot } = \sqrt {\frac{2}{m}\mu B} .\end{equation}

Moreover, we add a source term of particles coming from the wall with zero velocity. Equation (2.7) for $\langle\, f\rangle (r,{v_\parallel },{v_ \bot })$ becomes

(3.2)\begin{align} \partial \langle\, f\rangle /\partial t & = |{v_\parallel }|\dfrac{1}{r}\mathcal{L}[r\, D \mathcal{L} \langle\, f\rangle ] + |{v_\parallel }|L_0^{ - 1}D\mathcal{L}\langle\, f\rangle - Ze/m{\boldsymbol{b}_0}\boldsymbol{\cdot }{\boldsymbol{E}_{\textrm{in}}}\partial \langle\, f\rangle /\partial {v_\parallel } + C \nonumber\\ & \quad + \delta ({v_\parallel })\delta ({v_ \bot })/(2\mathrm{\pi }{v_ \bot })\, S(r), \end{align}

with

(3.3)\begin{equation}\mathcal{L} = {\left( {\frac{\partial }{{\partial r}}} \right)_{\varepsilon ,\mu }} = {\left( {\frac{\partial }{{\partial r}}} \right)_{{v_\parallel },{v_ \bot }}} - \frac{{{\textstyle{1 \over 2}}mv_ \bot ^2L_0^{ - 1} - Ze{E_a}}}{{m{v_\parallel }}}{\left( {\frac{\partial }{{\partial {v_\parallel }}}} \right)_{r,{v_ \bot }}} + \frac{1}{2}{v_ \bot }L_0^{ - 1}{\left( {\frac{\partial }{{\partial {v_ \bot }}}} \right)_{r,{v_\parallel }}},\end{equation}

where $S(r)$ is the radial dependence of the particle source. The form of the delta-function factors derives from the cylindrical symmetry in the velocity space introduced by the magnetic field, which gives ${\textrm{d}^3}v = 2\mathrm{\pi }{v_ \bot }\,\textrm{d}{v_ \bot }\,\textrm{d}{v_\parallel }$, with ${v_\parallel } \in [ - \infty , + \infty ]$ in the role of axial coordinate and ${v_ \bot } \in [0, + \infty ]$ in the role of radial coordinate. Since $D\sim o({\delta ^2})$, in (3.3) we approximate the radial derivative of B by means of the equilibrium component $L_0^{ - 1}$, and the radial derivative of the potential by its angular average ${E_a} ={-} \langle \partial \varPhi /\partial r\rangle$. This radial electric field ${E_a}$ is named ambipolar, because it equalizes the particle fluxes of electrons and ions, as we will see later.

Let $F({v_\parallel },{v_ \bot })$ be a generic polynomial function of the velocity components. After multiplying (3.2) by F, we integrate it over the velocity space by assuming $\langle\, f\rangle ,\mathcal{L}\langle\, f\rangle$ to vanish faster than any power of the velocity when $|{v_{{\parallel} , \bot }}|\to + \infty$. The integration over ${v_\parallel }$ can be limited to the interval $[0, + \infty ]$ by splitting any function $X({v_\parallel })$ of the parallel velocity into its even ${X_ + }({v_\parallel }) = [X({v_\parallel }) + X( - {v_\parallel })]/2$ and odd ${X_ - }({v_\parallel }) = [X({v_\parallel }) - X( - {v_\parallel })]/2$ components: $X({v_\parallel }) = {X_ + }({v_\parallel }) + {X_ - }({v_\parallel })$. Since both the operator (3.3) and the collisional factor $g({v_\parallel })$ defined in (2.9) are even in ${v_\parallel }$, and the latter also vanishes for ${v_\parallel } \to 0$, the velocity integration can be cast into the following form:

(3.4)\begin{gather} \dfrac{\partial }{{\partial t}}(\mathcal{F}_ + + \mathcal{F}_ - ) + \dfrac{\partial }{{r\partial r}}[r({\varGamma _ + } + {\varGamma _ - })] ={-} \dfrac{{Ze{E_a}}}{m}({{A}_ + } + {{A}_ - }) + \dfrac{1}{2}L_0^{ - 1}({{B}_ + } + {{B}_ - })\nonumber\\ + \dfrac{{Ze}}{m}{\boldsymbol{b}_0}\boldsymbol{\cdot }{\boldsymbol{E}_{\textrm{in}}}{\rm E} + \bar{C} + F(0,0)S(r), \end{gather}

(3.5)\begin{gather}\mathcal{F}_ \pm = 4\mathrm{\pi }\int_0^{ + \infty } {d{v_\parallel }} \int_0^{ + \infty } {\textrm{d}{v_ \bot }{v_ \bot }{F_ \pm }{{\langle\, f\rangle }_ \pm }} ,\end{gather}

(3.6)\begin{gather}{\varGamma _ \pm } ={-} 4\mathrm{\pi }\int_0^{ + \infty } {\textrm{d}{v_\parallel }{v_\parallel }} \int_0^{ + \infty } {\textrm{d}{v_ \bot }{v_ \bot }{F_ \pm }D\mathcal{L}{{\langle\, f\rangle }_ \pm }} ,\end{gather}

(3.7)\begin{gather}{{A}_ \pm } = 4\mathrm{\pi }\int_0^{ + \infty } {\textrm{d}{v_\parallel }} \int_0^{ + \infty } {\textrm{d}{v_ \bot }{v_ \bot }} \left[ {\frac{\partial }{{\partial {v_\parallel }}}{F_ \pm }} \right]D\mathcal{L}{\langle\, f\rangle _ \pm },\end{gather}

(3.8)\begin{gather}{{B}_ \pm } = 4\mathrm{\pi }\int_0^{ + \infty } {\textrm{d}{v_\parallel }} \int_0^{ + \infty } {\textrm{d}{v_ \bot }} \left[ {v_ \bot^3\frac{\partial }{{\partial {v_\parallel }}}{F_ \pm } - {v_\parallel }{v_ \bot }^2\frac{\partial }{{\partial {v_ \bot }}}{F_ \pm }} \right]D\mathcal{L}{\langle\, f\rangle _ \pm },\end{gather}

(3.9)\begin{gather}{E} = 4\mathrm{\pi }\int_0^{ + \infty } {\textrm{d}{v_\parallel }} \int_0^{ + \infty } {\textrm{d}{v_ \bot }{v_ \bot }} \left[ {{{\langle\, f\rangle }_ + }\frac{\partial }{{\partial {v_\parallel }}}{F_ - } + {{\langle\, f\rangle }_-{}}\frac{\partial }{{\partial {v_\parallel }}}{F_ + }} \right],\end{gather}

(3.10)\begin{gather}\bar{C} = 2\mathrm{\pi }\int_{ + \infty }^{ + \infty } {\textrm{d}{v_\parallel }} \int_0^{ + \infty } {\textrm{d}{v_ \bot }{v_ \bot }F\, C(\langle\, f\rangle ) } .\end{gather}

Note that the original integration interval of ${v_\parallel }$ is retained in the collision term (3.10).

4. Transport equations

As written above, we do not seek for a solution of (3.2). Instead, according to the ordering assumption previously discussed for (2.6), we approximate $\langle\, f\rangle$ by a local shifted Maxwellian:

(4.1)\begin{equation}\langle\, f\rangle \cong n{\left( {\frac{m}{{2\mathrm{\pi }T}}} \right)^{3/2}}{\textrm{e}^{ - m[{{({v_\parallel } - u)}^2} + v_ \bot ^2]/2T}},\end{equation}

with $n(r),T(r)$ the particle density and temperature (in energy units) and $u(r)$ the parallel drift velocity. All these quantities depend on the particle species. In such a way, the velocity integrals (3.5)–(3.9) can be carried out analytically, as done in previous publications (Harvey et al. Reference Harvey, McCoy, Hsu and Mirin1981; Rice et al. Reference Rice, Molvig and Helava1982; Predebon & Paccagnella Reference Predebon and Paccagnella2007), and (3.4) gives rise to transport equations for $n(r),T(r), u(r)$. Approximation (4.1) implies that collisions can establish a local thermal equilibrium. Accordingly, their effect has been included in (2.9).

Denoting ${v_{\textrm{th}}} = {(T/m)^{1/2}}$ the thermal velocity, and assuming $\zeta = {(u/{v_{\textrm{th}}})^2} = m{u^2}/T \ll 1$, we approximate the even/odd components in ${v_\parallel }$ of $\langle\, f\rangle$ with the following expansions:

(4.2)\begin{gather}{\langle\, f\rangle _ + } \cong \left( {1 - \frac{1}{2}\zeta + \frac{{m\, v_\parallel^2}}{{2T}}\zeta } \right){f_0},\quad {f_0} = n{\left( {\frac{m}{{2\pi T}}} \right)^{3/2}}{\textrm{e}^{ - m\, {v^2}/2T}},\end{gather}

(4.3)\begin{gather}{\langle\, f\rangle _ - } \approx \frac{m}{T}u\, {v_\parallel }{f_0} = {\zeta ^{1/2}}{v_\parallel }/{v_{\textrm{th}}}{f_0}.\end{gather}

Accordingly, in (3.6)–(3.8) we make use of the following relations:

(4.4)\begin{gather} \mathcal{L}{\langle\, f\rangle _ + } = \left[ \mathcal{A} + \zeta \mathcal{B} + \frac{{m{v^2}}}{{2T}}\frac{{T^{\prime}}}{T}\left( {1 - \frac{1}{2}\zeta } \right) + \zeta \left( {\frac{{m\, v_\parallel^2}}{{2T}}\mathcal{C} + \frac{{m\, v_\parallel^2}}{{2T}}\frac{{m{v^2}}}{{2T}}\frac{{T^{\prime}}}{T} - \frac{{m\, v_ \bot^2}}{{2T}}L_0^{ - 1}} \right) \right]{f_0},\end{gather}

(4.5a–c)\begin{gather}\mathcal{A} = \frac{{n^{\prime}}}{n} - \frac{3}{2}\frac{{T^{\prime}}}{T} - \frac{{Ze{E_a}}}{T},\quad \mathcal{B} ={-} \frac{{n^{\prime}}}{{2n}} + \frac{5}{4}\frac{{T^{\prime}}}{T} - \frac{{u^{\prime}}}{u} + \frac{3}{2}\frac{{Ze{E_a}}}{T},\nonumber\\ \quad \mathcal{C} = \frac{{n^{\prime}}}{n} - \frac{7}{2}\frac{{T^{\prime}}}{T} - \frac{{Ze{E_a}}}{T} + 2\frac{{u^{\prime}}}{u}, \end{gather}

(4.6)\begin{gather}\mathcal{L}{\langle\, f\rangle _ - } = \left( {\frac{{n^{\prime}}}{n} - \frac{5}{2}\frac{{T^{\prime}}}{T} + \frac{{m{v^2}}}{{2T}}\frac{{T^{\prime}}}{T} - \frac{{Ze{E_a}}}{T} + \frac{{u^{\prime}}}{u} - \frac{{{\textstyle{1 \over 2}}mv_ \bot^2L_0^{ - 1} - Ze{E_a}}}{{m{v_\parallel }^2}}} \right)\frac{m}{T}u\, {v_\parallel }{f_0}.\end{gather}

Only terms up to $o({\zeta ^{1/2}})$ have been retained in Predebon & Paccagnella (Reference Predebon and Paccagnella2007). Instead, here we push the approximation to $o(\zeta )$ to disclose velocity-related contributions in the particle and heat transport equations. Equation (3.4) is specialized to specific functions $F({v_\parallel },{v_ \bot })$ in the following subsections.

4.1. Particle transport

By taking ${F_ + } = 1, {F_ - } = 0$, we have ${{A}_ \pm } = {{B}_ \pm } = {E} = 0$ in (3.4). Since we are considering collisions that do not change the number of particles, $\bar{C} = 0$ as well. In this case (3.4) provides the following continuity equation:

(4.7)\begin{equation}\frac{{\partial n}}{{\partial t}} + \frac{\partial }{{r \partial r}}(r{\varGamma _n}) = S(r),\end{equation}

with the particle flux ${\varGamma _n}$ given by

(4.8)\begin{equation}{\varGamma _n} ={-} {D_n}\, n^{\prime}\left[ {1 + \zeta \left( {\frac{{{H_3}}}{{{H_2}}} - \frac{1}{2}} \right)} \right] + n\, {V_n},\quad ^{\prime} = \textrm{d}/\textrm{d}r,\end{equation}

where the stochastic particle diffusivity is

(4.9)\begin{equation}{D_n} = {\left( {\frac{2}{\mathrm{\pi }}} \right)^{1/2}}{v_{\textrm{th}}}{D_0}{H_2},\end{equation}

the convective velocity is

(4.10)\begin{align}{V_n} & ={-} {D_n}\left\{ \frac{{T^{\prime}}}{T}\left( {\frac{{ {H_3}}}{{{H_2}}} - \frac{1}{2}} \right) - \frac{{Ze{E_a}}}{T} + \zeta \left[ \frac{{T^{\prime}}}{T}\left( {\frac{{ {H_4}}}{{{H_2}}} - \frac{{3{H_3}}}{{{H_2}}} + \frac{3}{4}} \right) + \frac{{u^{\prime}}}{u}\left( {\frac{{2{H_3}}}{{{H_2}}} - 1} \right) \right.\right.\nonumber\\ &\quad \left.\left. +\, \frac{{Ze{E_a}}}{T}\left( {\frac{3}{2} - \frac{{ {H_3}}}{{{H_2}}}} \right) - L_0^{ - 1} \right] \right\}\end{align}

and ${H_k}$ are the following family of integral functions:

(4.11)\begin{equation}{H_k}(\alpha ) = \int_0^{ + \infty } {\textrm{d}x} \frac{{{x^k}}}{{\alpha + {x^2}}}\,{\textrm{e}^{ - x}},\end{equation}

stemming from the collision term $g({v_\parallel })$ defined in (2.9). Therefore, their argument depends on the particle species:

(4.12)\begin{equation}{\alpha _j} = \frac{1}{4}{L_K}/{\lambda _j}({v_{\textrm{th},j}}) = \frac{{Z_j^2{e^4}\ln (\varLambda ) {n_i}{L_K}}}{{4\mathrm{\pi }\varepsilon _0^2{{(2{T_j})}^2}}}.\end{equation}

We note that the combination of (4.8), (4.10) becomes (16) of Predebon & Paccagnella (Reference Predebon and Paccagnella2007) in the limit $\zeta \to 0$. The particle radial flux (4.8) is given by a diffusion term, proportional to $n^{\prime}$, plus a convective term. The integral functions (4.11) are decreasing functions of $\alpha$; hence collisions tend to decrease ${D_n}$. The collision-less case corresponds to $\alpha = 0$. Note that ${\alpha _i}\sim Z_i^2{\alpha _e}$ for the ions and ${\alpha _i}\sim {\alpha _e}$ for the common case of hydrogen isotope species. Therefore, the electron diffusivity is much larger than the ion diffusivity: ${D_{n,i}}\sim {v_{\textrm{th},i}}/{v_{\textrm{th},e}}{D_{n,e}}\sim {({m_e}/{m_i})^{1/2}}{D_{n,e}}$. Nonetheless, the particle fluxes (4.8) of ions and electrons are equalized, condition required by the plasma quasi-neutrality, by the ambipolar electric field ${E_a}$. The constraint ${\varGamma _{n,i}} = {\varGamma _{n,e}}$ implies the vanishing, apart from terms $o({({m_e}/{m_i})^{1/2}})$, of ${\varGamma _{n,e}}/{D_{n,e}}$. One gets

(4.13)\begin{align} \frac{{e{E_a}}}{{{T_e}}} & \cong{-} \left\{ \frac{{{{n^{\prime}}_e}}}{{{n_e}}} + \frac{{{{T^{\prime}}_e}}}{{{T_e}}}\left( {\frac{{{H_3}}}{{{H_2}}} - \frac{1}{2}} \right) + {\zeta_e}\left[ \frac{{{{n^{\prime}}_e}}}{{{n_e}}} + \frac{{{{T^{\prime}}_e}}}{{{T_e}}}\left( {\frac{{{H_4}}}{{{H_2}}} - \frac{{{H_3}}}{{{H_2}}} - {{\left( {\frac{{{H_3}}}{{{H_2}}}} \right)}^2}} \right)\right.\right.\nonumber\\ &\quad \left.\left. +\, \frac{{{{u^{\prime}}_e}}}{{{u_e}}}\left( {\frac{{2 {H_3}}}{{{H_2}}} - 1} \right) - L_0^{ - 1} \right] \right\}.\end{align}

In (4.13) the functions (4.11) are evaluated on the electron population: ${H_k} = {H_k}({\alpha _e})$. Note that the $L_0^{ - 1}$ term enters only as higher-order correction $(o(\zeta ))$ in (4.10) and (4.13).

4.2. Parallel momentum transport

By taking ${F_ + } = 0, {F_ - } = {m_j}{v_\parallel }$, we get from (3.4) the following equation for the parallel momentum:

(4.14)\begin{equation}\frac{\partial }{{\partial t}}(nmu) + \frac{\partial }{{r\partial r}}(r{\varGamma _u}) ={-} \frac{{Ze{E_a}}}{T}{\varTheta _u} + L_0^{ - 1}{\varLambda _u} + Zen{\boldsymbol{b}_0}\boldsymbol{\cdot }{\boldsymbol{E}_{\textrm{in}}} + {R_\parallel }.\end{equation}

The radial flux of parallel momentum ${\varGamma _u}$ on the left-hand side of (4.14) is the sum of a diffusion term and a convection term:

(4.15)\begin{gather}{\varGamma _u} = nm( - {D_u}u^{\prime} + u\, {V_u}),\end{gather}

(4.16)\begin{gather}{D_u} = 2{\left( {\frac{2}{\mathrm{\pi }}} \right)^{1/2}}{v_{\textrm{th}}}{D_0} {H_3} = 2 \frac{{{H_3}}}{{{H_2}}}{D_n},\end{gather}

(4.17)\begin{gather}{V_u} ={-} {D_u}\left[ {\frac{{n^{\prime}}}{n} + \frac{{T^{\prime}}}{T}\left( {\frac{{{H_4}}}{{{H_3}}} - \frac{3}{2}} \right) + \frac{{Ze{E_a}}}{T}\left( {\frac{{{H_2}}}{{2{H_3}}} - 1} \right) - \frac{1}{2}L_0^{ - 1}\frac{{{H_2}}}{{{H_3}}}} \right].\end{gather}

The coefficient (4.16) represents the stochastic kinematic viscosity. On the right-hand side of (4.14), the first term is a parallel force generated by the ambipolar electric field, with ${\varTheta _u}$ similar to a momentum flux:

(4.18)\begin{gather}{\varTheta _u} = n\, m\frac{{{H_2}}}{{2{H_3}}}({D_u}u^{\prime} + u\, {V_{u,\varTheta }}),\end{gather}

(4.19)\begin{gather}{V_{u,\varTheta }} = {D_u}\left[ {\frac{{n^{\prime}}}{n} + \frac{{T^{\prime}}}{T}\left( {\frac{{{H_3}}}{{{H_2}}} - \frac{3}{2}} \right) + \frac{{Ze{E_a}}}{T}\left( {\frac{{{H_1}}}{{2{H_2}}} - 1} \right) - \frac{1}{2}L_0^{ - 1}\frac{{{H_1}}}{{{H_2}}}} \right].\end{gather}

The second term in the right-hand side of (4.14) is a parallel force generated by the gradient of ${B_0}$, with ${\varLambda _u}$ similar to a momentum flux:

(4.20)\begin{gather}{\varLambda _u} = nm\frac{{{H_2}}}{{2{H_3}}}({D_u}u^{\prime} + u\, {V_{u,\varLambda }}),\end{gather}

(4.21)\begin{gather}{V_{u,\varLambda }} = {D_u}\left[ {\frac{{n^{\prime}}}{n} + \frac{{T^{\prime}}}{T}\left( {\frac{{{H_3}}}{{{H_2}}} - \frac{1}{2}} \right) + \frac{{Ze{E_a}}}{T}\left( {\frac{{{H_1}}}{{2{H_2}}} - 1} \right) - L_0^{ - 1}\frac{{{H_1}}}{{{H_2}}}} \right].\end{gather}

The function ${H_1}$, which enters (4.19), (4.21), diverges in the collision-less case $\alpha = 0$. This divergence stems from the singularity of (4.6) for ${v_\parallel } \to 0$, and it is symptomatic of the inconsistency of taking a collision-less limit in transport expressions derived by a local Maxwellian approximation of $\langle\, f\rangle$ (a distribution indeed established by collisions). Therefore, we have to keep $\alpha \ne 0$ to be consistent with our assumptions, and in this case ${H_1}$ is finite. Nonetheless, the divergence at $\alpha = 0$ is logarithmic, hence very weak: for instance, ${H_1}/{H_2}\sim 1$ for $\alpha \sim 1$, and ${H_1}/{H_2}\sim 10$ only at $\alpha \sim {10^{ - 9}}$. Moreover, the explicit computation discussed in the next paragraph estimates the two forces under discussion as almost negligible with $\alpha \sim 1$. Therefore, even the case of $\alpha \ne 0$ but very small is not problematic.

The third term on the right-hand side of (4.14) is the drive due to the applied electric field. Finally, the fourth term is the collisional transfer of parallel momentum:

(4.22)\begin{equation}{R_\parallel } = 2\mathrm{\pi }\int_{ - \infty }^{ + \infty } {\textrm{d}{v_\parallel }} \int_0^{ + \infty } {\textrm{d}{v_ \bot }{v_ \bot }m{v_\parallel }C} .\end{equation}

Following the procedure adopted in Bragiinkii (Reference Bragiinkii1965), we replace the time derivative of the density in (4.14) by the continuity equation (4.7):

(4.23)\begin{gather}nm\frac{\partial }{{\partial t}}u + m{\varGamma _n}\frac{\partial }{{\partial r}}u + \frac{\partial }{{r\partial r}}(r{\hat{\varGamma }_u}) ={-} \frac{{Ze{E_a}}}{T}{\varTheta _u} + L_0^{ - 1}{\varLambda _u} + Zen{\boldsymbol{b}_0}\boldsymbol{\cdot }{\boldsymbol{E}_{\textrm{in}}} + {R_\parallel } - muS,\end{gather}

(4.24)\begin{gather}{\hat{\varGamma }_u} = {\varGamma _u} - mu{\varGamma _n} = nm( - {D_u}u^{\prime} + u{\hat{V}_u}),\end{gather}

(4.25)\begin{gather}{\hat{V}_u} ={-} {D_u}\left[ {\frac{{n^{\prime}}}{n}\left( {1 - \frac{{{H_2}}}{{2{H_3}}}} \right) + \frac{{T^{\prime}}}{T}\left( {\frac{{{H_4}}}{{{H_3}}} - 2 + \frac{{{H_2}}}{{4{H_3}}}} \right) + \frac{{Ze{E_a}}}{T}\left( {\frac{{{H_2}}}{{{H_3}}} - 1} \right) - \frac{1}{2}L_0^{ - 1}\frac{{{H_2}}}{{{H_3}}}} \right].\end{gather}

Since $u = {\zeta ^{1/2}}{v_{\textrm{th}}}$, terms of order $\zeta$ are neglected in (4.25). Equation (4.23) is formally akin to the Braginskii motion equation ((2.2) of Bragiinkii Reference Bragiinkii1965), with an extra term related to the particle source S. In particular, two correspondences can be identified:

(4.26)\begin{equation}m{\varGamma _n}\frac{\partial }{{\partial r}}u \leftrightarrow nm\boldsymbol{V}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{V},\quad \frac{\partial }{{r\partial r}}(r{\hat{\varGamma }_u}) \leftrightarrow \boldsymbol{\nabla }\boldsymbol{\cdot }\overline{\overline{\boldsymbol{\varPi }}} ,\end{equation}

with $\boldsymbol{V}, \overline{\overline{\boldsymbol{\varPi }}}$ the Braginskii mean velocity and stress tensor, respectively.

4.3. Heat transport

Since we do not distinguish between parallel and perpendicular temperatures in the distribution function (4.1), the meaningful second-order moment is ${F_ + } = {\textstyle{1 \over 2}}m{v^2} = {\textstyle{1 \over 2}}m(v_\parallel ^2 + v_ \bot ^2),{F_ - } = 0$. Note an exact cancellation of the terms in the square bracket of expression (3.8) for ${{B}_ \pm }$: therefore, $L_0^{ - 1}$ gives no contribution in (3.4). From (3.4), one gets the following energy transport equation:

(4.27)\begin{equation}\frac{\partial }{{\partial t}}\left( {\frac{3}{2}nT + \frac{1}{2}nm{u^2}} \right) + \frac{\partial }{{r\partial r}}(r{\varGamma _E}) = Ze{E_a}{\varGamma _n} + Zenu{\boldsymbol{b}_0}\boldsymbol{\cdot }{\boldsymbol{E}_{\textrm{in}}} + {R_\parallel }u + Q,\end{equation}

where ${\varGamma _E}$ is the energy flux, not written here for the sake of brevity, and

(4.28)\begin{equation}Q = 2\mathrm{\pi }\int_{- \infty }^{ + \infty } {\textrm{d}{v_\parallel }} \int_0^{ + \infty } {\textrm{d}{v_ \bot }{v_ \bot }} \frac{1}{2}m[{({v_\parallel } - u)^2} + v_ \bot ^2]C\end{equation}

is the heat produced by collisions. Following again Bragiinkii (Reference Bragiinkii1965), we express the time derivative of ${\textstyle{1 \over 2}}nm{u^2}$ by (4.7) and (4.23) to obtain the heat transport equation:

(4.29)\begin{equation}\frac{\partial }{{\partial t}}\left( {\frac{3}{2}nT} \right) + \frac{\partial }{{r\partial r}}(r{\varGamma _Q}) ={-} {\hat{\varGamma }_u}\frac{\partial }{{\partial r}}u + Ze{E_a}{\varGamma _n} + \frac{{Ze\, {E_a}}}{T}{\varTheta _u}u - L_0^{ - 1}{\varLambda _u}u + Q + \frac{1}{2}m{u^2}S.\end{equation}

The heat flux is

(4.30)\begin{gather}{\varGamma _Q} = {\varGamma _E} - \frac{1}{2}m{u^2}{\varGamma _n} - {\hat{\varGamma }_u}u ={-} {D_T}nT^{\prime}\left[ {{H_4} + \frac{{{H_3}}}{2} + \frac{{{H_2}}}{2} + \zeta \left( {\frac{{9}}{4}{H_3} - 4{H_4} + {H_5}} \right)} \right]\nonumber\\ + nT{V_T},\end{gather}

(4.31)\begin{gather}{D_T} = {\left( {\frac{2}{\mathrm{\pi }}} \right)^{1/2}}{v_{\textrm{th}}}{D_0} = {D_n}H_2^{ - 1},\end{gather}

(4.32)\begin{gather} {V_T} ={-} {D_T}\left\{ \left( {\dfrac{{n^{\prime}}}{n} - \dfrac{{Ze{E_a}}}{T}} \right)({H_3} + {H_2}) + \zeta \left[ {\dfrac{{n^{\prime}}}{n}\left( {{H_4} - \dfrac{3}{2}{H_3}} \right) + \dfrac{{u^{\prime}}}{u}(2{H_4} - {H_2} - {H_3})} \right. \right.\nonumber\\ \left. \left. \quad + \dfrac{{Ze{E_a}}}{T}\left( {\dfrac{5}{2}{H_3} - {H_4}} \right) - L_0^{ - 1}({H_3} + {H_2}) \right] \right\},\end{gather}

with ${D_T}$ the stochastic thermal diffusivity. Since ${H_4}(\alpha ) = 2 - \alpha {H_2}(\alpha )$, we note that the flux (4.30) becomes expression (31) of Predebon & Paccagnella (Reference Predebon and Paccagnella2007) for $\zeta \to 0$.

The terms on the right-hand side of (4.29) are heating sources. In particular, by comparing (4.29) with the Braginskii heat equation ((1.23) of Bragiinkii Reference Bragiinkii1965), and taking into account the second of (4.26), we can identify the viscous heating as

(4.33)\begin{equation}{Q_{visc}} \to - {\hat{\varGamma }_u}\frac{\partial }{{\partial r}}u\textrm{.}\end{equation}

The presence of the particle source S in (4.29) is explained by the fact that zero-velocity particles from the wall are seen with energy ${\textstyle{1 \over 2}}m{u^2}$ in the plasma reference. However, this contribution is estimated as negligible (see § 5).

4.4. The Onsager reciprocal relations

The above expressions for particle, momentum and heat flux satisfy the Onsager reciprocal relations (Onsager Reference Onsager1931), a general property of the thermodynamics of the transport processes, stemming from the principle of microscopic reversibility. This result is discussed in full because it represents a verification of the model. Assuming a small deviation from thermal equilibrium due to the ‘thermodynamic forces’ ${X_i}$ (such as the gradients), the ‘thermodynamic fluxes’ ${\varGamma _i}$ are linearly related to ${X_i}$ by kinetic coefficients ${L_{ik}}$:

(4.34)\begin{equation}{\varGamma _i} = \mathop \sum \limits_k {L_{ik}}{X_k}.\end{equation}

Fluxes and forces are ‘conjugate’ if the entropy production for unit of volume can be expressed in the form

(4.35)\begin{equation}\theta = \mathop \sum \limits_i {\varGamma _i}{X_i}.\end{equation}

The reciprocal relations of the kinetic coefficients apply to conjugate fluxes and forces. In the presence of a magnetic field they are (Bragiinkii Reference Bragiinkii1965)

(4.36)\begin{equation}{L_{ik}}(\boldsymbol{B}) = {L_{ki}}( - \boldsymbol{B}),\end{equation}

when both forces ${X_i}, {X_k}$ are even or odd functions of the particle velocity,

(4.37)\begin{equation}{L_{ik}}(\boldsymbol{B}) ={-} {L_{ki}}( - \boldsymbol{B}),\end{equation}

if one force is odd and the other is even. To show that the present theory satisfies the constraints (4.36), (4.37), we first identify the pairs of conjugate thermodynamic fluxes and forces. Following Bragiinkii (Reference Bragiinkii1965), we start from the entropy per particle (j-species label omitted):

(4.38)\begin{equation}s = {\textstyle{3 \over 2}}\ln (T) - \ln (n) + \textrm{const}\textrm{.}\end{equation}

Then we compute the transport equation for the entropy per unit volume, $ns$, by exploiting the continuity equation for density (4.7), and the heat transport equation (4.29). Discarding the particle source S, one gets

(4.39)\begin{equation}\frac{\partial }{{\partial t}}(ns) + \frac{\partial }{{r\partial r}}\left[ {r\left( {s{\varGamma_n} - \frac{5}{2}{\varGamma_n} + \frac{{{\varGamma_Q}}}{T}} \right)} \right] = \theta = \mathop \sum \limits_{i = 1,4} {\varGamma _i}{X_i} + \frac{Q}{T}.\end{equation}

The entropy production $\theta$ is given by the collisional heating and by the following four pairs of conjugate fluxes and thermodynamic forces:

(4.40)\begin{equation} \left. {\begin{array}{c@{}} {{\varGamma_1} = {\varGamma_n},\quad {X_1} = \dfrac{3}{2}\dfrac{{T^{\prime}}}{T} - \dfrac{{n^{\prime}}}{n} + \dfrac{{Ze{E_a}}}{T};\quad {\varGamma_2} = {{\hat{\varGamma }}_u},\quad {X_2} ={-} \dfrac{{u^{\prime}}}{T};}\nonumber\\ {{\varGamma_3} = {\varGamma_Q},\quad {X_3} ={-} \dfrac{{T^{\prime}}}{{{T^2}}};\quad {\varGamma_4} = \dfrac{{Ze{E_a}}}{T}{\Theta_u} - L_0^{ - 1}{\Lambda_u},\quad {X_4} = \dfrac{u}{T}.} \end{array}} \right\}\end{equation}

Note that ${\varGamma _4}$ is a force density (which enters (4.23)), but we keep the term ‘thermodynamic flux’ in this context (as done in Bragiinkii (Reference Bragiinkii1965) for the collisional momentum transfer R, for instance). The forces ${X_2}, {X_3}$ are expressed by gradients of $u, T$ respectively, and they are conjugate of the parallel momentum and heat fluxes ${\varGamma _2}, {\varGamma _3}$, as expected. Instead, ${X_1}$, the conjugate of the particle flux ${\varGamma _1}$, depends on $T^{\prime}$ and ${E_a}$, in addition to $n^{\prime}$. The ambipolar electric field enters ${X_1}$ because ${\varGamma _1}Ze{E_a}$ is a heating source in (4.29). Moreover, we speculate that $T^{\prime}$ enters ${X_1}$ because particles carry energy. As written above, ${\varGamma _4}$ is actually a parallel force density, and it is conjugate of ${X_4}$ because ${\varGamma _4}{X_4}T$ is a heating source in (4.29).

By exploiting equations (4.8)–(4.10), (4.16), (4.18)–(4.21), (4.24), (4.25), (4.30)–(4.32), it is possible to establish the linear relationship (4.34) between ${\varGamma _i}$ and ${X_i}$ defined in (4.40), with the following ${L_{ik}}$. The diagonal coefficients are

(4.41)\begin{gather}\left. \begin{array}{c@{}} {{L_{11}} = n{D_n}\left[ {1 + \zeta \left( {\dfrac{{{H_3}}}{{{H_2}}} - \dfrac{1}{2}} \right)} \right];\quad {L_{22}} = 2nmT{D_n}\dfrac{{{H_3}}}{{{H_2}}};}\\ {L_{33}} = n{T^2}{D_n}\left[ {\dfrac{{{H_4}}}{{{H_2}}} + \dfrac{{2 {H_3}}}{{{H_2}}} + 2 + \zeta \left( {\dfrac{{{H_5}}}{{{H_2}}} - \dfrac{5}{2}\dfrac{{{H_4}}}{{{H_2}}}} \right)} \right];\\ \quad {L_{44}} = nmT{D_n}\dfrac{{{H_1}}}{{{H_2}}}\left[ {\dfrac{1}{2}{{\left( {\dfrac{{Ze\, {E_a}}}{T}} \right)}^2} + L_0^{ - 2} - \dfrac{{Ze\, {E_a}}}{T}L_0^{ - 1}} \right]. \end{array} \right\}\end{gather}

The off-diagonal coefficients are symmetric:

(4.42)\begin{align}\left. {\begin{array}{c@{}} {{L_{12}} = {L_{21}} = nmu{D_n}\left( {\dfrac{{2{H_3}}}{{{H_2}}} - 1} \right);\quad {L_{13}} = {L_{31}} = nT{D_n}\left[ {\dfrac{{{H_3}}}{{{H_2}}} + 1 + \zeta \left( {\dfrac{{{H_4}}}{{{H_2}}} - \dfrac{3}{2}\dfrac{{{H_3}}}{{{H_2}}}} \right)} \right];}\\ {{L_{14}} = {L_{41}} ={-} nmu{D_n}\left[ {\dfrac{{Ze{E_a}}}{T} - L_0^{ - 1}} \right];\quad {L_{23}} = {L_{32}} = nmTu{D_n}\left( {\dfrac{{2{H_4}}}{{{H_2}}} - 1 - \dfrac{{{H_3}}}{{{H_2}}}} \right);}\\ {{L_{24}} = {L_{42}} ={-} nmT{D_n}\left[ {\dfrac{{Ze{E_a}}}{T} - L_0^{ - 1}} \right];\quad {L_{34}} = {L_{43}} ={-} nmTu{D_n}\left[ {\dfrac{{Ze{E_a}}}{T}\dfrac{{{H_3}}}{{{H_2}}} - L_0^{ - 1}\left( {\dfrac{{{H_3}}}{{{H_2}}} + 1} \right)} \right].} \end{array}} \right\}\end{align}

The equalities ${L_{13}} = {L_{31}}$, ${L_{24}} = {L_{42}}$ correspond to (4.36), because both the forces involved are even or odd in the velocity. Instead, the other coefficients relate even and odd forces: the equalities ${L_{12}} = {L_{21}}$, ${L_{14}} = {L_{41}}$, ${L_{23}} = {L_{32}}$, ${L_{34}} = {L_{43}}$ indeed correspond to (4.37), because the average parallel velocity u, contained in these coefficients, is odd with $\boldsymbol{B}$. In conclusion, this theory fulfils the Onsager reciprocal relations.

The off-diagonal coefficients formalize the interdependence of the transport processes in a thermodynamic system (Onsager Reference Onsager1931). Indeed, the fluxes derived in the previous paragraph depend on all the quantities $n^{\prime}/n, T^{\prime}/T, u^{\prime}/u, Ze{E_a}/T, L_0^{ - 1}$, so that they are represented by convective velocities, in addition to diffusive components. By inspection of (4.41), (4.42), we note that the terms ${L_{12}}{X_2}, {L_{14}}{X_4}$ are $o(\zeta )$ with respect to ${L_{11}}{X_1}, {L_{13}}{X_3}$, so that the particle flux ${\varGamma _1}$ mostly depends on ${X_1}, {X_3}$. Instead, the terms ${L_{2k}}{X_k}$ share the same order, so that the parallel momentum flux ${\varGamma _2}$ depends on all the ${X_k}$ in the same way. The terms ${L_{32}}{X_2}, {L_{34}}{X_4}$ are $o(\zeta )$ with respect to ${L_{33}}{X_3}, {L_{31}}{X_1}$, so that the heat flux ${\varGamma _3}$ mostly depends on ${X_1}, {X_3}$. The terms ${L_{4k}}{X_k}$ share the same order, so that ${\varGamma _4}$ depends on all the ${X_k}$ in the same way. Of course, the actual magnitude of the various terms ${L_{ik}}{X_k}$ can be inferred only by an explicit computation, and we comment on that in the next paragraph.

As far as the collisional heating Q is concerned, we can adopt the Braginskii expressions (Bragiinkii Reference Bragiinkii1965) in the present context without violating (4.36), (4.37). Braginskii's formulas are simplified here, because: (i) the mean velocity is taken only in the direction parallel to $\boldsymbol{B}$ and (ii) due to the cylindrical geometry, the gradients parallel to $\boldsymbol{B}$ can be discarded as an effect produced by the radial magnetic field $({\nabla _\parallel }\sim {a^{ - 1}}o(\delta ))$. For electrons and main ions one gets

(4.43a,b)\begin{gather}{Q_e} = {R_{{\parallel} ,e}}({u_i} - {u_e}) - {Q_{ei}},\quad {Q_i} = {Q_{ei}},\end{gather}

(4.44a,b)\begin{gather}{R_{{\parallel} ,e}} = {(e{n_e})^2}{\eta _\parallel }({u_i} - {u_e}),\quad {Q_{ei}} \cong 1.37 \times {10^{ - 30}}\frac{{{m_e}}}{{{m_i}}}{Z_i}\ln \varLambda n_e^2({T_e} - {T_i})/T_e^{3/2},\end{gather}

with ${\eta _\parallel }$ the parallel resistivity and ${Q_{ei}}$ the heat exchanged between the two species (the ${Q_{ei}}$ expression is in SI units, except T in eV). For the electrons, (4.39), (4.43) introduce the thermodynamic force $({u_i} - {u_e})/{T_e}$ and the conjugate flux ${R_{{\parallel} ,e}}$. The same equations, summed for electrons and ions, introduce the force $({T_e} - {T_i})/({T_e}{T_i})$ and the conjugate flux ${Q_{ei}}$. These two pairs of fluxes and forces are completely decoupled by the previous ones (4.40). Therefore, (4.43), (4.44) can be included without violating the Onsager reciprocal relations.

5. Example of solution

Within the cylindrical geometry approximation, the present theory can be applied to the stochastic region of any magnetic configuration. The above equations for particle and heat transport, the latter without the terms related to the velocity u, have been already integrated into a transport model for the RFP (Predebon et al. Reference Predebon, Paccagnella, Puiatti and Valisa2009), which is currently applied to the RFX-mod device ($R = 2m$, $a = 0.459m$) (Zuin et al. Reference Zuin2017). Indeed, this theory is suitable for the RFP in the so-called multiple-helicity regime (MH), which is characterized by many $m = 1$ tearing modes with comparable amplitudes, whose overlap gives rise to a wide stochastic region. This is not true in the quasi-single-helicity regime (QSH), where the innermost resonant mode grows at the expense of the others, thus healing the magnetic chaos (Lorenzini et al Reference Lorenzini2009). The inclusion of parallel momentum transport, the novelty of the present derivation, still allows finding plausible solutions for RFX-mod in MH conditions. The unavoidable uncertainty about several input quantities, as well as the lack of the transport terms due to the drift velocity perpendicular to the magnetic field, prevent a detailed comparison of this model with the experiment. Nonetheless, the solutions are in the expected range of values. We consider a case at plasma current ${I_p}\sim 1\textrm{ MA}$ (shot 29509 as reference). We adopt SI units, but the temperature is given in eV.

First of all, at least an indicative estimate of the collision-less magnetic diffusion coefficient ${D_0}$, given by (2.8), is required. The RFX-mod sensors’ layout allows determination of the harmonics $\tilde{b}_r^{m,n}$ for the core-resonant $m = 1$ modes and for the reversal-resonant $m = 0$ modes (${r^{0,n}} \approx 0.4$ in this shot), in the interval $n = 1$–24. We focus on $m = 1$, whose amplitude is much larger than that of $m = 0$. The harmonics $\tilde{b}_r^{1,n}$ are estimated at the resonant radii ${r^{1,n}}$ by a reconstruction, from the edge measurements, of both the equilibrium (assumed force-free, ${\boldsymbol{J}_0} \propto {\boldsymbol{B}_0}$) and the perturbation radial profile (Zanca & Terranova Reference Zanca and Terranova2004). Their amplitudes are plotted in figure 1(a): even in a MH case, the innermost resonant $n = 7$ is twice as large as the following $n = 8$ (for a QSH case it would be 20 times larger). The same figure plots the Chirikov parameter $s(n,n + 1) = {\textstyle{1 \over 2}}({w^{1,n}} + {w^{1,n + 1}})/({r^{1,n + 1}} - {r^{1,n}})$ (Rechester & Rosenbluth Reference Rechester and Rosenbluth1978; D'Angelo & Paccagnella Reference D'Angelo and Paccagnella1996), with ${w^{1,n}}$ the island width (the increase of s with n is due to the resonances getting denser). The overlapping condition $s > 1$, which defines a stochastic magnetic field, is largely satisfied. As written before, we do not apply any resonance broadening method to regularize the singular functions $\delta (m - nq) = \delta (r - {r^{1,n}})/|nq^{\prime}|$. We simply average them by radial integration of (2.8) over suitable intervals centred about the resonant radii ${r^{1,n}}$. Since magnetic islands overlap $(s > 1)$, we cannot take ${w^{1,n}}$ as amplitudes of these intervals. To avoid overlapping, we adopt neighbouring intervals $\mathrm{\Delta }{r^{1,n}} = ({r^{1,n + 1}} - {r^{1,n - 1}})/2$, with extrema halfway to ${r^{1,n + 1}}$ and ${r^{1,n - 1}}$. For the innermost resonant $n = 7$ the left extremum is set at $r = 0$. The contribution of the $m = 0$ modes cannot be computed in this way, due to the impossibility of defining a radial interval to average: in fact, both resonances of $m = 1, n > 0$ and $m = 1, n < 0$ (a secondary branch resonating between the reversal and the plasma edge, not detectable by the RFX-mod sensors) become dense approaching the reversal. The average values of ${D_0}$ on these intervals give the stepwise profile plotted in figure 1(b) (black). We have verified that a more regular profile makes solving the equations easier. Therefore, we choose the profile plotted in blue $({D_0}(0) = 2.5 \times {10^{ - 4}}\ \textrm{m})$. We also keep a finite value at the edge $({D_0}(a) = 5 \times {10^{ - 5}}\ \textrm{m})$, since numerical analyses of magnetic diffusion in MH conditions give a stochastic region in the whole plasma: see figure 9(b) of Innocente et al. (Reference Innocente, Lorenzini, Terranova and Zanca2017).

Figure 1. (a) Values of $|\tilde{b}_r^{1,n}({r^{1,n}})|$ for modes $m = 1,\, n = 7$–24, averaged over a significant time interval for shot 29509 (black diamonds), and stochasticity parameter from these amplitudes (red squares). (b) Stepwise profile of ${D_0}(r)$ obtained from averaging equation (2.8) (black; the vertical segments mark the resonant radii), alongside the continuous profile taken to solve the equations (blue).

The Kolmogorov–Lyapunov length ${L_K}$, which enters ${\alpha _j}$ (equation (4.12)), is here fixed by a scaling found numerically in D'Angelo & Paccagnella (Reference D'Angelo and Paccagnella1996): ${L_K}(r) = {(0.127D_0^{0.275})^{ - 1}}$. An increasing radial profile with ${L_K}(0) \approx 80\ \textrm{m}$ is obtained. This scaling is close to the analytical result derived in Rechester et al. (Reference Rechester, Rosenbluth and White1979).

Stochastic diffusion is expected to be the dominant transport mechanism in the plasma core, but electrostatic turbulence can emerge in the edge region where ${D_0}$ becomes smaller. Indeed, the edge RFX-mod particle transport is compatible with the Bohm diffusivity ${D_{\textrm{Bohm}}} = {\textstyle{1 \over {16}}}(T/B)$ (Antoni et al. Reference Antoni1998; Spolaore et al. Reference Spolaore, Antoni, Spada, Bergsåker, Cavazzana, Drake, Martines, Regnoli, Serianni and Vianello2004); hence this additional transport mechanism should be included when solving the equations. Diffusivity ${D_{\textrm{Bohm}}}$ is typically negligible compared with the electron stochastic diffusivity, but, in the case under examination, it is of the same order of magnitude $({\sim} 10\ {\textrm{m}^2}\ {\textrm{s}^{ - 1}})$ as the ion stochastic diffusivity (see figure 2d). One option is to take a background of Bohm diffusion in the whole plasma, besides the stochastic transport, letting the equations to weigh the two different mechanisms. We now present the solution obtained by this assumption. Another option is to restrict the Bohm diffusion only to the edge region: this second case is briefly discussed later on.

Figure 2. (a) Normalized profiles obtained for density (red with diamonds), ${T_e}$ (blue with crosses) and ${T_i}$ (black). The normalized particle source is also reported (red dashed). (b) Parallel current from velocity and density solution (black) and from equilibrium reconstruction (blue with crosses). (c) Toroidal (black) and poloidal (blue with crosses) components of ${u_i}$. (d) Electron (blue with crosses) and ion (black) stochastic particle diffusivity, alongside the ion Bohm coefficient (red with diamonds). These solutions are obtained by superimposing the Bohm diffusion on the stochastic transport, in the whole plasma.

Taking stationary conditions $(\partial /\partial t = 0)$, the transport equations become second-order ordinary differential equations in the radial coordinate, for density, velocity and temperature. Regularity at $r = 0$ requires the vanishing of the on-axis radial derivative of these quantities. Therefore, one free parameter is left in each equation to make the profiles decreasing to small values towards $r = a$. The equations solved are described below.

1. (i) Equations (4.7)–(4.13) for the normalized density ${n_i}/{n_i}(0)$ of the hydrogen ion $({Z_i} = 1)$. As explained before, the diffusion term $- {D_{\textrm{Bohm},i}}{n^{\prime}_i}$ is added to the flux (4.8). Without a model for the source, we fix a plausible $S(r)/{n_i}(0)$ profile localized at the edge, but we leave a global coefficient ${S_0}$ as free parameter. We take the same normalized density profile for all the particle species, with on-axis electron density ${n_e}(0) = 2.5 \times {10^{19}}\ {\textrm{m}^{ - 3}}$.

2. (ii) Equations (4.14)–(4.22) for the parallel velocities of electrons ${u_e}$ and ions ${u_i}$, with ${\boldsymbol{E}_{\textrm{in}}} = {E_\phi }{\boldsymbol{e}_\phi },\,{E_\phi } = 25/(2\mathrm{\pi }R)\ \textrm{V}\ {\textrm{m}^{ - 1}}$. The electron collisional momentum transfer ${R_{{\parallel} ,e}}$ is expressed by (4.44) with Spitzer parallel resistivity: ${\eta _\parallel } \cong 5.2 \times {10^{ - 5}}\ln \varLambda {Z_{\textrm{eff}}}N({Z_{\textrm{eff}}})T_e^{ - 3/2}$ (Wesson Reference Wesson2004). Due to the significant dependence of ${R_{{\parallel} ,e}}$ on the plasma effective charge ${Z_{\textrm{eff}}}$, with $N({Z_{\textrm{eff}}})$ a corrective function, it is important to include impurity species, which, though present only in trace, can raise ${Z_{\textrm{eff}}}$. We consider the Carbon impurity (the material of the RFX-mod first wall), with mean charge ${Z_C} = 5$, and we take ${Z_{\textrm{eff}}} = 2.5$, a suitable value for this density level (Carraro et al. Reference Carraro, Gazza, Marrelli, Martin, Puiatti, Scarin, Valisa, Zaccaria and Zaniol2005). Equation (4.14) is considered separately for the electrons, whereas those for the main ion and the Carbon impurity are summed together, with assumption of same velocity and temperature for the two species, and exploiting the momentum conservation ${R_{{\parallel} ,i}} + {R_{{\parallel} ,C}} ={-} {R_{{\parallel} ,e}}$. The on-axis values ${u_e}(0)$, ${u_i}(0)$ are free parameters. The background Bohm diffusion terms $- {n_j}{m_j}{D_{\textrm{Bohm},j}}{u^{\prime}_j}$ are added to the momentum flux (4.15) for the different species.

3. (iii) Equations (4.29)–(4.32) for the electron ${T_e}$ and ion ${T_i}$ temperatures. The collisional heating is given by (4.43), (4.44). The background diffusion terms $- {n_j}{D_{\textrm{Bohm},j}}{T^{\prime}_j}$ are added to the heat flux (4.30). The on-axis values ${T_e}(0)$, ${T_i}(0)$ are free parameters.

All the above equations are coupled together.

To summarize, several inputs of the equations are taken from the experiments: the equilibrium magnetic field and the perturbed radial magnetic field, used to estimate ${D_0}$ and ${L_K}$, the on-axis electron density, the applied electric field and ${Z_{\textrm{eff}}}$. Instead, ${S_0}, {u_e}(0), {u_i}(0), {T_e}(0), {T_i}(0)$ are left as free parameters, tuned to have density, velocity and temperature profiles with values at $r = a$ below 15 % of the on-axis inputs. Therefore, they turn out to be an output of the solution.

The normalized solution profiles for density and temperature are shown in figure 2(a). The very steep gradient at $r \approx a$ is an artefact due to transport coefficients getting very small there (they depend on $T$). The density has the mildly hollow profile typical of RFX-mod, and the particle source $S(r)$ (red dashed line in the same figure) is very similar to the edge-localized profiles reported in figures 5 and 6 of Lorenzini et al. (Reference Lorenzini, Auriemma, Canton and Carraro2006), with a maximum about $2.4 \times {10^{23}}\ {\textrm{m}^{ - 3}}\ {\textrm{s}^{ - 1}}$. As far as the on-axis temperature is concerned, we get ${T_e}(0) \approx 335\ \textrm{eV}$, ${T_i}(0) \approx 195\ \textrm{eV}$. The ${T_e}$ value is compatible with typical measurements at ${I_p}\sim 1\textrm{ MA}$ in MH conditions (see figure 6 of Innocente et al. Reference Innocente, Alfier, Canton and Pasqualotto2009), whereas for ${T_i}$ we cannot say, lacking reliable diagnostics of this quantity in RFX-mod. The collisional heating ${Q_i} = {Q_{ei}}$ dominates the other source terms on the right-hand side of (4.29) for the ions. Hence, this model does not contain any relevant anomalous heating source for the ions. Apart from the high-density regime, ${Q_i}$ is not strong enough to equalize the ion and electron temperatures, hence our result ${T_i}(0) \approx 0.58\ {T_e}(0)$. During the magnetic reconnection events, characteristic of the RFP, it is often observed that ${T_i}> {T_e}$ (Scime et al. Reference Scime, Cekic, Den Hartog, Hokin, Holly and Watts1992; Gobbin et al. Reference Gobbin2022), supporting the existence of an important transient anomalous ion heating, outside the quasi-linear description of the present work. The identification of the process responsible for that is still a matter of investigation: see McChesney, Stern & Bellan (Reference McChesney, Stern and Bellan1987), Guo, Paccagnella & Romanelli (Reference Guo, Paccagnella and Romanelli1994) and Fiksel et al. (Reference Fiksel, Almagri, Chapman, Mirnov, Ren, Sarff and Terry2009) as examples.

In the core region the parallel velocity profiles are steeper than those of the temperature, because the ratio between the source and the stochastic diffusion coefficient is larger for the momentum than for the temperature. Moreover, the two terms with ${\boldsymbol{E}_{\textrm{in}}},{R_\parallel }$ dominate the other sources on the right-hand side of (4.14). We get ${u_e}(0) \approx 1540\ \textrm{km}\ {\textrm{s}^{ - 1}}$, ${u_i}(0) \approx 27\ \textrm{km}\ {\textrm{s}^{ - 1}}$, corresponding to Mach numbers ${\approx} 5.17,\, 0.09$ respectively (use is made of the sound speed defined in formula (6.1) below). The density and velocity solutions allow computing the parallel current density profile, ${j_\parallel } = e{n_e}({u_i} - {u_e})$, which is compared with the magnetic equilibrium reconstruction ${j_{{\parallel} ,\textrm{eq}}}$ in figure 2(b). The good deal of uncertainty about several quantities of the model (for instance ${D_0}$, ${L_k}$) would leave space for an adjustment of the same in order to get a better agreement between ${j_\parallel }$ and ${j_{{\parallel} ,\textrm{eq}}}$. Nonetheless, the result is already satisfactory. Since ${u_i} \ll {u_e}$, ${j_\parallel }$ involves mostly ${u_e}$. As far as ${u_i}$ is concerned, a comparison with the RFX-mod spectroscopic measurement of the impurity flow (Zaniol et al. Reference Zaniol, Bolzonella, Cavazzana, Innocente, Zanca and Zuin2015) is more problematic, because the present model misses the perpendicular drift velocity, particularly relevant in the edge region (Spolaore et al. Reference Spolaore, Antoni, Spada, Bergsåker, Cavazzana, Drake, Martines, Regnoli, Serianni and Vianello2004; Vianello et al. Reference Vianello, Spada, Antoni, Spolaore, Serianni, Regnoli, Cavazzana, Bergsåker and Drake2005). However, the toroidal and poloidal components of ${u_i}$ plotted in figure 2(c) are in the correct experimental range ($< 20\ \textrm{km}\ {\textrm{s}^{ - 1}}$ for impurities which emit in the region $r > a/2$). In conclusion, the solution is compatible with the experiment.

The solution so obtained weakly depends on the profile shape chosen for ${D_0}(r)$. For instance, a profile with the same on-axis and edge values as that considered so far (blue continuous line in figure 1b), but much flatter in the central region $r < a/2$, and steeper for $r > a/2$, leaves the density profile and the source $S(r)$ almost unchanged and gives about $10\,\%$ variations for the other quantities: ${T_e}(0) \approx 377\ \textrm{eV}$, ${T_i}(0) \approx 205\ \textrm{eV}$, ${u_e}(0) \approx 1750\ \textrm{km}\ {\textrm{s}^{ - 1}}$, ${u_i}(0) \approx 30\ \textrm{km}\ {\textrm{s}^{ - 1}}$.

As far as the stochastic transport coefficients are concerned, both the kinematic viscosity (4.16) and the thermal diffusivity (4.31) are proportional to the particle diffusivity (4.9). Therefore, in figure 2(d) we plot ${D_{n,e}},{D_{n,i}}$ alongside ${D_{\textrm{Bohm},i}}$ as terms of comparison: the electron stochastic diffusivity dominates the ion stochastic diffusivity, which instead is of the same order of magnitude as the Bohm diffusivity in the central region $(r < a/2)$.

The stochastic fluxes derived in § 4 are all of the form ${\varGamma _x} \propto {D_x}x^{\prime} + x{V_x}$, with ${V_x}$ mainly depending on $n^{\prime}/n$ and $T^{\prime}/T$ (taking into account expression (4.13) for the ambipolar electric field), apart from the $L_0^{ - 1}$ dependence of ${\hat{V}_u}$, and the $o(\zeta )$ corrections. Figures 3(a) and 3(b) plot the ratio between the convective and diffusive components of the particle, heat and momentum flux, for ions and electrons, respectively, in the core region (the edge region is not plotted, due the unrealistically large $n^{\prime}/n$ and $T^{\prime}/T$ obtained there). This ratio is everywhere close to unity for the electron particle flux, due to the ambipolar electric field (4.13), which makes ${\varGamma _{n,e}}/{D_{n,e}} = o({({m_e}/{m_i})^{1/2}})$. Therefore, it is not plotted. The ratio is small for the momentum, because the density and temperature profiles (which determine ${V_x}$) are rather flat in the core, while the velocity is not. Recalling the discussion below (4.42), this means that the momentum flux is mainly driven by the diagonal term ${L_{22}}{X_2}$, though all ${L_{2k}}{X_k}$ are of the same order a priori. Nonetheless, convection is not entirely negligible for the ion momentum flux. For heat flux, diffusion exceeds convection: indeed, the diagonal term ${L_{33}}{X_3}$ is larger than ${L_{31}}{X_1}$ (see figure 3c) whereas ${L_{32}}{X_2},\, {L_{34}}{X_4}$ (not plotted) are $o(\zeta )$ corrections (as stated below (4.42)). For ion particle flux, convection exceeds diffusion: indeed ${L_{13}}{X_3}$ is larger than ${L_{11}}{X_1}$ (see figure 3c), whereas ${L_{12}}{X_2},{L_{14}}{X_4}$ (not plotted) are $o(\zeta )$ corrections. This large convection produces the hollow density profile typically observed in RFX-mod. The parallel force density ${\varGamma _4}$ (not plotted) is estimated as negligible; hence the cumulative effects of the terms ${L_{4k}}{X_k}$ are very small. We cannot say whether all these properties are specific to the stochastic RFP core, or more general.

Figure 3. (a) Ratio between the convective and diffusive components of the stochastic fluxes in the core region $(r/a < 0.8)$ for the ion quantities. (b) Ratio between the convective and diffusive components of the stochastic fluxes in the core region $(r/a < 0.8)$ for the electron quantities, but for the density. (c) Some of the terms ${L_{ik}}{X_k}$ for the main ion, suitably normalized to make them comparable. All these quantities refer to the solution displayed in figure 2.

A comment is made about the terms depending on the magnetic field radial variation $L_0^{ - 1}$, which is one of the new elements of this analysis. The one entering ${\hat{V}_u}$ (equation (4.25)) is not negligible, whereas those in ${V_n}$, ${V_T}$ are higher-order $(o(\zeta ))$ corrections. As written above, the related force and heat source, on the right-hand side of (4.23), (4.29) respectively, are estimated as negligible, as those produced by ${E_a}$. Therefore, the $L_0^{ - 1}$ effect is relevant only in the convective component of the momentum flux.

Finally, we summarize what changes in the solution when restricting the Bohm diffusion to the very edge region $(r \ge 0.9a)$. As expected, ${T_e},\,{u_e}$ are almost unaffected, whereas ${T_i},\,{u_i}$ increase: ${T_i}(0) \approx 0.8 {T_e}(0)$, ${u_i}(0) \approx 54\ \textrm{km}\ {\textrm{s}^{ - 1}}$ (nonetheless, in the region $r > a/2$ we still get the correct experimental range $< 20\ \textrm{km}\ {\textrm{s}^{ - 1}}$ for the toroidal and poloidal components of ${u_i}$). The hollowness of the density profile increases (but it is still plausible for RFX-mod), due to the removal of the additional Bohm diffusion in the core. In conclusion, the ion solution is sensitive to the inclusion of anomalous transport processes different from the stochastic one.