TY - JOUR
AU - Norman, C, A
AB - ABSTRACT We reconsider the theory of particle transport in a scattering medium, allowing for relativistic flow velocities. The theory uses a mixed set of variables, with position and time measured in the laboratory frame, and particle energy and momentum measured in the fluid rest frame, the reference frame where scattering is assumed to be elastic. We give a new derivation for the fictitious force terms in the equation of motion that are present if the fluid rest frame is not an inertial frame. By using a 3 + 1 notation we discuss the physical interpretation of the different terms in the fictitious force. It is shown that different approaches to the problem of particle propagation in a magnetized medium due to Skilling and Kulsrud are largely equivalent. We extend known results for non-relativistic flows to include the effects of cross-field diffusion for cosmic rays in a magnetized plasma. We also carefully consider the correct form of the diffusion approximation for scattering, and show that the resulting equations can be cast in conservation form. diffusion, MHD, relativistic processes, scattering, cosmic rays 1 INTRODUCTION The transport theory of charged particles or photons in a scattering and moving medium has many applications in astrophysics. In the case of energetic charged particles (cosmic rays) the scattering is due to random magnetic fields (magnetohydrodynamic, MHD waves), as was realized a long time ago by Parker (1965), Dolginov & Toptygin (1967), and Gleeson & Axford (1967). A very systematic and general derivation of the cosmic ray transport equation was given by Skilling (1975), which included the effect of the wave motion with respect to a non-relativistic flow and did allow for relativistic particles with particle momentum p ≫ mc. More or less in parallel Kulsrud (1983) formulated a theory in the limit of the vanishing gyro-radius for the kinetic description of charged particles in a moving magnetized plasma, assuming non-relativistic particles and slow motion (bulk velocity much less than the speed of light). The application of these two approaches (i.e. Skilling’s formulation versus Kulsrud’s approach) was mostly limited to the astrophysics and the plasma physics community, respectively, even though they are (in the appropriate limit) essentially equivalent in physical content, as we will show. A generalization of Skilling’s approach to a relativistic flow was first formulated by Webb (1985). However, as was the case with the calculations of both Skilling (1975) and Kulsrud (1983), the final equations of Webb (1985) are obtained by taking an average over the gyro-phase, the phase of the rapid circular motion of a charged particle in the plane perpendicular to the magnetic field. This is allowed if the gyro-frequency Ω =qB/γmc of a charge q of mass m and Lorentz factor γ in a magnetic field B is very large so that the gyration period is usually much smaller than any of the hydrodynamical time-scales. This averaging procedure is mathematically convenient and allows for a relatively simple final form of the transport equation. But as already noted by Skilling (1975), the gyro-average eliminates physical processes such as diffusion across the magnetic field and viscous interactions in a shear flow which both depend on anisotropies in the gyro-phase. Diffusion along the magnetic field is unaffected by the gyro-averaging and is obtained in a simple fashion. The importance of viscous effects of cosmic rays that are coupled to a flow by frequent scattering was already noted by Jokipii & Williams (1992) and Williams & Jokipii (1993). For photons, radiation transport in a moving medium has been considered by many authors. The classical results and principles of radiation hydrodynamics can be found in the books of Pomraning (1973), Mihalas & Weibel Mihalas (1984), and Castor (2004). The viscous effects of a radiation fluid were considered by Weinberg (1971) and more recently by Coughlin & Begelman (2014). In this paper we reconsider the general theory with an emphasis on charged particles in a magnetized medium. We extend existing calculations to include anisotropies in the gyro-phase, a necessary step for the calculation of the full diffusion tensor and of viscous effects. We will derive the fundamental equations and show the relation to known results for non-relativistic flows. We also cast the equations in a more intuitive form, employing a 3 + 1 split for space–time. This approach allows an easy comparison with known results for non-relativistic flows. We use the results to derive a transport equation for the particles, showing the effects of diffusion and drift with respect to the bulk fluid. We also consider (and prove) the general covariance of the final equations. The treatment of viscous effects and a detailed discussion of the form of the stress tensor for charged particles in a magnetized medium are deferred to the companion paper (Paper 2). Radiation transport is treated in Paper 3. 2 BASIC EQUATIONS: USE OF MIXED VARIABLES In what follows we assume a flow with four-velocity |$U^{\mu } = (\Gamma\,, {\boldsymbol{U}}) = (\sqrt{1 + U^2}\,, {\boldsymbol{U}})$| in the laboratory frame (LF). This flow scatters particles (or photons). Scattering is assumed to be elastic in the fluid rest frame (FRF). We employ units with c = 1 throughout the paper, but reinstate c when it clarifies results. The assumption of elastic scattering in the FRF is a simplification in the case of scattering of cosmic rays by MHD waves: It neglects the phase velocity of the waves with respect to the fluid. As shown by Skilling (1975), this becomes important when the intensity of the waves responsible for the scattering (in particular Alfvén waves) is not equal in the two possible propagation directions: along and against the large-scale magnetic field. Since scattering takes the simplest form in the fluid rest frame, it is mathematically convenient to use a set of mixed variables, where position and time are measured in the LF but particle momentum and energy are measured in the FRF. In what follows we will employ the following notation: The LF momentum four-vector is denoted using a Greek index notation and the four-momentum in the FRF is denoted using a hatted Roman index notation (e.g. |$\hat{a}, \, \hat{b}, \, \ldots$|⁠): \begin{eqnarray*} p^{\mu } = ({\cal E}, \, {\boldsymbol{p}}), p^{\hat{a}} = (\bar{{\cal E}}, \, \bar{{\boldsymbol{p}}}) . \end{eqnarray*} (1) The particle energy |$\bar{{\cal E}}$| and the momentum three-vector |$\bar{{\boldsymbol{p}}}$| in the FRF are barred in order to distinguish them from the corresponding quantities in the LF. The same convention is used when dealing with other non-covariant quantities such as three-vectors. For most of the paper we employ a 3 + 1 split of the equations. This makes a direct comparison with previous results more easy. We will show however that the final results can always be expressed covariantly in terms of scalars, four-vectors, and four-tensors. That means that, even though the calculations are done neglecting gravity, they are straightforwardly generalized to the general-relativistic case. As long as we neglect the effects of gravity (i.e. in special relativity), pμ and |$p^{\hat{a}}$| are connected by Lorentz transformations. The case with gravity is considered briefly in Appendix A. For technical reasons (again, see Appendix A) we employ the covariant momentum four-vectors |$p_{\mu } = ({\cal E}\ ,, \, - {\boldsymbol{p}})$| and |$p_{\hat{a}} = (\bar{{\cal E}}, - \bar{{\boldsymbol{p}}})$|⁠, defined formally through |$p_{\mu } = \eta_{\mu \nu } \, p^{\nu }$|⁠, |$p_{\hat{a}} = \eta_{\hat{a}\hat{b}} \, p^{\hat{b}}$| with |$\eta_{\mu \nu } = \eta_{\hat{a}\hat{b}} = {\rm diag}(1\,, -1\,, -1\,, -1)$| the Minkowski tensor, the metric tensor of flat space–time. Using the Einstein summation convention for double indices we can write: \begin{eqnarray*} p_{\mu } = {\rm {E}}_{\mu }^{\, \hat{a}} \, p_{\hat{a}}\ ,, \quad p_{\hat{a}} = {\rm {E}}_{\, \hat{a}}^{\mu } \, p_{\mu } . \end{eqnarray*} (2) Here |$\rm {E}_{\mu }^{\, \hat{a}}$| and |$\rm {E}_{\, \hat{a}}^{\mu }$| are Lorentz transformations (see for instance Landau & Lifshitz 1975, ch. 1) between the two frames. Represented as 4 × 4 matrices \begin{eqnarray*} \rm {E}_{\mu }^{\, \hat{a}} = \left(\begin{array}{@{}c@{\quad}c@{\quad}c@{}}\Gamma & & - {\boldsymbol{U}} \\ & & \\ - {\boldsymbol{U}} & & {\displaystyle {\mathbf{I}} + \frac{{\boldsymbol{U}} \,{\boldsymbol \otimes}\,{\boldsymbol{U}}}{\Gamma + 1}} \\ \end{array} \right)\,, \quad \rm {E}_{\, \hat{a}}^{\mu } = \left(\begin{array}{@{}c@{\quad}c@{\quad}c@{}}\Gamma & & {\boldsymbol{U}} \\ & & \\ {\boldsymbol{U}} & & {\displaystyle {\mathbf{I}} + \frac{{\boldsymbol{U}} \,{\boldsymbol \otimes}\,{\boldsymbol{U}}}{\Gamma + 1}} \\ \end{array} \right). \end{eqnarray*} (3) They are each other's inverse, or alternatively |${\rm E}^{\mu }_{\hat{a}}({\boldsymbol{U}}) = {\rm E}_{\mu }^{\hat{a}}(- {\boldsymbol{U}})$|⁠. The notation |${\boldsymbol{U}} \,{\boldsymbol \otimes}\,{\boldsymbol{U}}$| is the direct product, with components |$({\boldsymbol{U}} \,{\boldsymbol \otimes}\,{\boldsymbol{U}})_{ij} = U^{i} U^{j}$|⁠. We use unhatted roman upper indices |$i\,, j\,, \ldots$|⁠, with |$(i\,, j\,, \ldots) = 1\,, 2\,, 3$|⁠, when considering the spatial components of four-vectors in the LF. The following relations hold: \begin{eqnarray*} {\rm E}_{\mu }^{\hat{a}} \, {\rm E}^{\mu }_{\hat{b}} = \eta^{\hat{a}}_{\hat{b}} = \delta^{\hat{a}}_{\hat{b}}, \quad {\rm E}_{\mu }^{\hat{a}} \, {\rm E}^{\nu }_{\hat{a}} = \eta^{\nu }_{\mu } = \delta^{\nu }_{\mu }, \quad \eta_{\mu \nu } \, {\rm E}^{\mu }_{\hat{a}}{\rm E}^{\nu }_{\hat{b}} = \eta_{\hat{a}\hat{b}}, \quad \eta_{\hat{a}\hat{b}} \, {\rm E}_{\mu }^{\hat{a}} {\rm E}_{\nu }^{\hat{b}} = \eta_{\mu \nu } . \end{eqnarray*} (4) The last two relations guarantee the invariance of the scalar product of two four-vectors, in particular |$p_{\mu } p^{\mu } = {\cal E}^2 - p^2 = p_{\hat{a}} p^{\hat{a}} = \bar{{\cal E}}^2 -\bar{p}^2 = m^2$|⁠. 3 THE FICTITIOUS FORCE Generally, the FRF is not an inertial frame. The LF on the other hand is assumed to be an inertial frame. The use of the FRF four-momentum |$p_{\hat{a}}$| implies the appearance of, in classical terms, a fictitious force term in the equation of motion. In relativity the equivalence principle implies that this can also be thought of as space–time curvature in the accelerated frame. We will return to this point in Section 7. 3.1 Special-relativistic calculation In special relativity the form of the fictitious force can be derived by using the generalization of Newton’s law of inertia for a free particle. In the absence of gravity or other forces the LF four-momentum is conserved. This means we can write \begin{eqnarray*} \frac{{\rm d} p_{\mu }}{{\rm d} \lambda } = \frac{{\rm d}}{{\rm d} \lambda } \left({\rm E}_{\mu }^{\hat{a}} \, p_{\hat{a}} \right) = 0 . \end{eqnarray*} (5) Here λ is the affine parameter along the orbit of a particle. This yields \begin{eqnarray*} {\rm E}_{\mu }^{\hat{a}} \, \frac{{\rm d} p_{\hat{a}}}{{\rm d} \lambda } + \frac{{\rm d} {\rm E}_{\mu }^{\hat{a}}}{{\rm d} \lambda } \, p_{\hat{a}} = 0 . \end{eqnarray*} (6) Using the first relation in equation (4) one finds, after multiplication by |${\rm E}^{\mu }_{\hat{c}}$| and re-ordering: \begin{eqnarray*} {\rm E}^{\mu }_{\hat{c}} {\rm E}_{\mu }^{\hat{a}} \, \frac{{\rm d} p_{\hat{a}}}{{\rm d} \lambda } = \frac{{\rm d} p_{\hat{c}}}{{\rm d} \lambda } = - {\rm E}^{\mu }_{\hat{c}} \, \frac{{\rm d} {\rm E}_{\mu }^{\hat{a}}}{{\rm d} \lambda } \, p_{\hat{a}} . \end{eqnarray*} (7) As d/dλ is the proper time derivative along the particle’s orbit it equals \begin{eqnarray*} \frac{{\rm d}}{{\rm d} \lambda } \equiv \frac{p^{\nu }}{m} \, \frac{\partial }{\partial x^{\nu }} = {\rm E}^{\nu }_{\, \hat{b}} \, \frac{p^{\hat{b}}}{m} \, \frac{\partial }{\partial x^{\nu }} . \end{eqnarray*} (8) One finds: \begin{eqnarray*} \frac{{\rm d} p_{\hat{c}}}{{\rm d} \lambda } = - {\rm E}^{\mu }_{\hat{c}} \, {\rm E}^{\nu }_{\hat{b}} \, \frac{\partial {\rm E}_{\mu }^{\hat{a}}}{\partial x^{\nu }} \, \frac{p_{\hat{a}} p^{\hat{b}}}{m} \equiv K_{\hat{c}} . \end{eqnarray*} (9) This gives the fictitious force on a free particle, now expressed entirely in terms of FRF momentum components. Additional force terms can now be simply added. 3.2 General-relativistic calculation If one allows for space–time curvature with a metric gμν, equation (5) is replaced by (see for instance Carroll 2004, ch.3 and defining uα = pα/m and |${\rm d}/{\rm d} \lambda = u^{\mu } \, (\partial /\partial x^{\mu })$|⁠) \begin{eqnarray*} \frac{{\rm d} p_{\mu }}{{\rm d} \lambda } = \Gamma^{\alpha }_{\mu \beta } \, u_{\alpha } \, p^{\beta } = - \frac{1}{2m} \, \left(\frac{\partial g^{\alpha \beta }}{\partial x^{\mu }} \right) \, p_{\alpha } p_{\beta } . \end{eqnarray*} (10) In this case the FRF coordinates can be put in quasi-Minkowskian form (so that locally|$g_{\hat{a}\hat{b}} = \eta_{\hat{a}\hat{b}}$|⁠) by using the tetrad formalism (see Appendix A). The transformations |${{\rm E}}^{\hat{a}}_{\, \mu }$| and the inverse |${{\rm E}}_{\mu }^{\, \hat{a}}$| now satisfy relations analogous to those in equation (4), where one replaces ημν (ημν) by gμν (gμν). In particular one has |$\eta_{\hat{a}\hat{b}} \, {\rm E}^{\hat{a}}_{\, \mu } \, {\rm E}^{\hat{b}}_{\, \nu } = g_{\mu \nu }$| and |$g_{\mu \nu } \, {\rm E}^{\mu }_{\, \hat{a}} \, {\rm E}^{\nu }_{\, \hat{b}} = \eta_{\hat{a}\hat{b}}$|⁠. Then expressing all LF momentum variables in equation (10) in terms of FRF momentum variables one arrives at \begin{eqnarray*} \frac{{\rm d} p_{\hat{c}}}{{\rm d} \lambda } = {\rm E}^{\mu }_{\, \hat{c}} \, \left(\frac{\partial {\rm E}_{\nu }^{\, \hat{a}}}{\partial x^{\mu }} - \frac{\partial {\rm E}_{\mu }^{\, \hat{a}} \,}{\partial x^{\nu }} \right) \, \frac{p_{\hat{a}} \, p^{\hat{b}}}{m} . \end{eqnarray*} (11) A formal proof of this result is given in Appendix A, using a variational principle. The extra term in equation (11) with respect to special-relativistic result (9) (i.e. the first term inside the bracket) vanishes in flat space–time. In this paper we will use another route to obtain results that can be applied in general relativity, namely by proving the covariance of the final equations so that this transition is one from ordinary derivatives to covariant derivatives. 3.3 Fictitious force calculation for a non-inertial fluid rest frame In flat space–time all variations of |${\rm E}^{\mu }_{\hat{a}}$| and its inverse stem from changes in the velocity |${\boldsymbol{V}}$| of the fluid in the laboratory frame. This means that we can write the derivatives that appear in equation (9) in terms of the derivatives with respect to the spatial components Ui = ΓVi (with |$i = 1\ ,, 2\ ,, 3$|⁠) of the four-velocity of the FRF in the LF: \begin{eqnarray*} \frac{{\rm d} {\rm E}_{\mu }^{\hat{a}}}{{\rm d} \lambda } = \frac{\partial {\rm E}_{\mu }^{\hat{a}}}{\partial U^{k}} \, \frac{{\rm d} U^{k}}{{\rm d} \lambda }, \quad \frac{\partial {\rm E}_{\mu }^{\hat{a}}}{\partial x^{\nu }} = \frac{\partial {\rm E}_{\mu }^{\hat{a}}}{\partial U^{k}} \, \frac{\partial U^{k}}{\partial x^{\nu }} \end{eqnarray*} (12) The following relations hold: \begin{eqnarray*} \frac{\partial {\rm E}^{\hat{0}}_{0}}{\partial U^{k}} = \frac{\partial \Gamma }{\partial U^{k}} = \frac{U^{k}}{\Gamma } = V^{k}, \quad \frac{\partial {\rm E}^{\hat{i}}_{0}}{\partial U^{k}} = - \frac{\partial U^{i}}{\partial U^{k}} = - \delta_{ik}, \quad \frac{\partial {\rm E}^{\hat{j}}_{i}}{\partial U^{k}} = \frac{1}{\Gamma + 1} \, \left(\delta_{ik} U^{j} + \delta_{jk} U^{i} - \frac{U^{i} U^{j} U^{k}}{\Gamma (\Gamma + 1)} \right) . \end{eqnarray*} (13) We are somewhat sloppy with the index notation, replacing |$(\hat{i}\,, \hat{j})$| by |$(i\,, j)$| when no confusion can arise. All unhatted roman indices |$i\,, j\,, \ldots$| run from 1 to 3. 3.3.1 Energy change The change of FRF energy due to the fictitious force follows from |$p_{0} = \bar{{\cal E}}$| and equation (9) with |$\hat{c} = \hat{0}$|⁠. Using equations (12) and (13), together with |${\boldsymbol{U}} = \Gamma \, {\boldsymbol{V}}$| and |$\Gamma = \sqrt{1 + |{\boldsymbol{U}}|^2}$|⁠, one finds after a significant amount of algebra (see Appendix B) \begin{eqnarray*} \frac{{\rm d} \bar{{\cal E}}}{{\rm d} \lambda } = - (\Gamma + 1) \, \left[ \bar{{\boldsymbol{p}}}{\boldsymbol{\cdot}}\frac{{\rm d}}{{\rm d} \lambda } \left(\frac{{\boldsymbol{U}}}{\Gamma + 1} \right) \right]. \end{eqnarray*} (14) 3.3.2 Momentum change In a similar fashion, the momentum change due to the fictitious force is found to be (again see Appendix B for details) \begin{eqnarray*} \frac{{\rm d} \bar{{\boldsymbol{p}}}}{{\rm d} \lambda } = - (\Gamma + 1) \, \bar{{\cal E}}\, \frac{{\rm d}}{{\rm d} \lambda } \left(\frac{{\boldsymbol{U}}}{\Gamma + 1} \right) + \bar{{\boldsymbol{p}}}{\boldsymbol{\times}}{\boldsymbol \omega }, \end{eqnarray*} (15) where \begin{eqnarray*} {\boldsymbol \omega } \equiv \frac{{\boldsymbol{U}}}{\Gamma +1} {\boldsymbol{\times}}\frac{{\rm d} {\boldsymbol{U}}}{{\rm d} \lambda } . \end{eqnarray*} (16) The second term on the right-hand side of equation (15) is reminiscent of the Coriolis term in a rotating reference frame in classical mechanics. Indeed it stems from the change in the direction of |${\boldsymbol{U}}$|⁠. In addition it is a purely relativistic effect. Writing |${\boldsymbol{U}} = U \, {\boldsymbol{\hat{\ell}}}$|⁠, with |${\boldsymbol{\hat{\ell}}}$| a unit vector along |${\boldsymbol{U}}$| and |$U = | {\boldsymbol{U}} |$|⁠, one immediately finds using U2 = Γ2 − 1 \begin{eqnarray*} {\boldsymbol \omega } = \frac{U^2}{\Gamma + 1} \left({\boldsymbol{\hat{\ell}}}{\boldsymbol{\times}}\frac{{\rm d} {\boldsymbol{\hat{\ell}}}}{{\rm d} \lambda } \right) = (\Gamma - 1) \, \left({\boldsymbol{\hat{\ell}}}{\boldsymbol{\times}}\frac{{\rm d} {\boldsymbol{\hat{\ell}}}}{{\rm d} \lambda } \right). \end{eqnarray*} (17) For non-relativistic flows with U ≪ 1 we have Γ − 1 ≃ V2/2c2 ≪ 1, where we reinstated c for clarity, and the |$\bar{{\boldsymbol{p}}}{\boldsymbol{\times}}{\boldsymbol \omega }$| term can be neglected. It is easily checked that results of equations (14) and (15) maintain the Lorentz invariance of |$\bar{{\cal E}}^2 - \bar{p}^2$| since \begin{eqnarray*} \bar{{\cal E}}\, \frac{{\rm d} \bar{{\cal E}}}{{\rm d} \lambda } - \bar{{\boldsymbol{p}}}{\boldsymbol{\cdot}}\frac{{\rm d} \bar{{\boldsymbol{p}}}}{{\rm d} \lambda } = 0 . \end{eqnarray*} (18) 3.3.3 The origin of the |$\bar{{\boldsymbol{p}}}{\boldsymbol{\times}}{\boldsymbol \omega }$| term The appearance of this term in the fictitious force is related to the fact that two Lorentz boosts L1 and L2, where the boost directions are not equal, do not commute, so L1L2 ≠ L2L1. The difference L1L2 − L2L1 corresponds to a spatial rotation. The same effect is, for instance, responsible for the Thomas precession and the similar appearance of a cross-product in the theory of Fermi–Walker transport of a tetrad of quasi-Minkowskian unit vectors carried by an accelerated observer (e.g. Misner, Thorne, & Wheeler 1973, ch. 6) when the acceleration is not uniform. A simple physical argument for the appearance and form of this term goes as follows: Consider a free particle with constant LF energy |${\cal E}$| and momentum |${\boldsymbol{p}}$|⁠. The FRF momentum can be expressed using the Lorentz transformation from the LF to the FRF: \begin{eqnarray*} \bar{{\boldsymbol{p}}}= {\boldsymbol{p}} - U \, {\cal E}\, {\boldsymbol{\hat{\ell}}}+ \left(\Gamma - 1 \right) \, ({\boldsymbol{p}} {\boldsymbol{\cdot}}{\boldsymbol{\hat{\ell}}}) \, {\boldsymbol{\hat{\ell}}}. \end{eqnarray*} (19) Now consider an infinitesimal change in the direction of |${\boldsymbol{U}} \equiv U \, {\boldsymbol{\hat{\ell}}}$| with |$U \equiv | {\boldsymbol{U}} |$|⁠: \begin{eqnarray*} {\boldsymbol{\hat{\ell}}}\, \Rightarrow {\boldsymbol{\hat{\ell}}}+ \delta {\boldsymbol{\hat{\ell}}}, \, {\boldsymbol{\hat{\ell}}}{\boldsymbol{\cdot}}\delta {\boldsymbol{\hat{\ell}}}= 0 \,, \end{eqnarray*} (20) so that the vector |${\boldsymbol{\hat{\ell}}}+ \delta {\boldsymbol{\hat{\ell}}}$| remains a unit vector to first order. The induced linear change in |$\bar{{\boldsymbol{p}}}$| with U, |${\cal E}$|⁠, and |${\boldsymbol{p}}$| constant is \begin{eqnarray*} \delta \bar{{\boldsymbol{p}}}= -U \, {\cal E}\, \delta {\boldsymbol{\hat{\ell}}}+ \left(\Gamma - 1 \right) \, \left[ \, ({\boldsymbol{p}} {\boldsymbol{\cdot}}\delta {\boldsymbol{\hat{\ell}}}) \, {\boldsymbol{\hat{\ell}}}+ ({\boldsymbol{p}} {\boldsymbol{\cdot}}{\boldsymbol{\hat{\ell}}}) \, \delta {\boldsymbol{\hat{\ell}}}\, \right]. \end{eqnarray*} (21) Next we write |${\cal E}$| and |${\boldsymbol{p}}$| in terms of the (unperturbed) |$\bar{{\cal E}}$| and |$\bar{{\boldsymbol{p}}}$| using the Lorentz transformation back to the LF: \begin{eqnarray*} {\cal E}= \Gamma \, \bar{{\cal E}}+ U \, (\bar{{\boldsymbol{p}}}{\boldsymbol{\cdot}}{\boldsymbol{\hat{\ell}}}), \quad {\boldsymbol{p}} = \bar{{\boldsymbol{p}}}+ U \, \bar{{\cal E}}\, {\boldsymbol{\hat{\ell}}}+ \left(\Gamma - 1 \right) \, (\bar{{\boldsymbol{p}}}{\boldsymbol{\cdot}}{\boldsymbol{\hat{\ell}}}) \, {\boldsymbol{\hat{\ell}}}. \end{eqnarray*} (22) Substituting this into relation (21) yields, again to first order in |$\delta {\boldsymbol{\hat{\ell}}}$|⁠: \begin{eqnarray*} \delta \bar{{\boldsymbol{p}}}& = & - U \, \bar{{\cal E}}\, \delta {\boldsymbol{\hat{\ell}}}+ \left(\Gamma - 1 \right) \, \left[ \, (\bar{{\boldsymbol{p}}}{\boldsymbol{\cdot}}\delta {\boldsymbol{\hat{\ell}}}) \, {\boldsymbol{\hat{\ell}}}- (\bar{{\boldsymbol{p}}}{\boldsymbol{\cdot}}{\boldsymbol{\hat{\ell}}}) \, \delta {\boldsymbol{\hat{\ell}}}\, \right] = - U \, \bar{{\cal E}}\, \delta {\boldsymbol{\hat{\ell}}}+ \left(\Gamma - 1 \right) \, \left[ \, \bar{{\boldsymbol{p}}}{\boldsymbol{\times}}\left({\boldsymbol{\hat{\ell}}}{\boldsymbol{\times}}\delta {\boldsymbol{\hat{\ell}}}\, \right) \, \right]. \end{eqnarray*} (23) The |$(\Gamma - 1) \, [ \bar{{\boldsymbol{p}}}{\boldsymbol{\times}}\left({\boldsymbol{\hat{\ell}}}{\boldsymbol{\times}}\delta {\boldsymbol{\hat{\ell}}}\, \right)]$| term is the source of the cross-product term involving |${\boldsymbol \omega }$| in equation (15). 4 VLASOV EQUATION Consider charged particles with rest mass m and charge q. The Vlasov equation for such particles in the laboratory frame reads (Klimontovich 1960; Hakim 2011; see also Debbasch & van Leeuwen 2009a, b, who discuss the general-relativistic case): \begin{eqnarray*} \frac{p^\mu }{m} \, \frac{\partial {\cal F}}{\partial x^{\mu }} + \frac{q}{m} {\rm \mathit{ F}}_{\mu }^{\, \nu } p_{\nu } \, \frac{\partial {\cal F}}{\partial p_{\mu }} = {\rm St}({\cal F}) . \end{eqnarray*} (24) Here |${\cal F}(x\,, p)$| is the distribution function in terms of the eight phase-space variables |$(x^{\mu }\,, p_{\mu })$|⁠, defined in such a way that (see Klimontovich 1960) |${\rm d} N^{\mu } = {\cal F}(x\,, p) \, (p^{\mu }/m) {\rm d}^4 x \, {\rm d}^4 p$| is the number flux four-vector of particles residing in the eight-dimensional phase-space element. It is related to the more familiar distribution |$f({\boldsymbol x }\,, t\,, {\boldsymbol{p}})$| in seven-dimensional phase space through \begin{eqnarray*} {\cal F}(x^{\mu }, \quad p_{\mu }) = 2 m \Theta (p_{0}) \, f({\boldsymbol x }, t, {\boldsymbol{p}}) \, \delta (p^2 - m^2) \,, \end{eqnarray*} (25) with |$\Theta (p_{0}) = \Theta ({\cal E})$| the Heaviside step function and |$p^2 = p^{\mu } p_{\mu } = {\cal E}^2 - | {\boldsymbol{p}} |^2$|⁠. The delta function keeps particles on mass shell while the step function Θ(p0) ensures positive energy. The tensor |${\rm \mathit{ F}}_{\mu }^{\, \nu }$| is the electromagnetic Faraday tensor, see for instance Jackson (1998), and |$(q/m) \, {\rm \mathit{ F}}_{\mu }^{\, \nu } \, p_{\nu } = \gamma ({\boldsymbol E } {\boldsymbol{\cdot}}{\boldsymbol{V}}\,, - [{\boldsymbol E } + {\boldsymbol{V}} {\boldsymbol{\times}}{\boldsymbol B }])$| is the covariant Lorentz four-force with |$({\boldsymbol E }\,, {\boldsymbol B })$| the electromagnetic field in the laboratory frame and |${\boldsymbol{V}} = {\boldsymbol{p}}/\gamma m$| the particle velocity. |${\rm St}({\cal F})$| is a scattering term that will be specified below. If we change momentum variables from LF four-momentum pμ to FRF four-momentum |$p_{\hat{a}}$|⁠, we have to use \begin{eqnarray*} \frac{\partial }{\partial x^{\mu }} \Longrightarrow \frac{\partial }{\partial x^{\mu }} + \frac{\partial p_{\hat{a}}}{\partial x^{\mu }} \, \frac{\partial }{\partial p_{\hat{a}}} = \frac{\partial }{\partial x^{\mu }} + \frac{\partial {\rm E}^{\nu }_{\hat{a}}}{\partial x^{\mu }}\, p_{\nu } \, \frac{\partial }{\partial p_{\hat{a}}}, \end{eqnarray*} (26) where ∂/∂xμ is now with |$(\bar{{\cal E}}, \bar{{\boldsymbol{p}}})$| fixed. Also \begin{eqnarray*} \frac{\partial }{\partial p_{\mu }} = \frac{\partial p_{\hat{a}}}{\partial p_{\mu }} \, \frac{\partial }{\partial p_{\hat{a}}} = {\rm E}^{\mu }_{\hat{a}} \, \frac{\partial }{\partial p_{\hat{a}}} . \end{eqnarray*} (27) The Jacobian of the momentum transformation is |${\rm det}({\rm E}_{\hat{a}}^{\mu }) = 1$|⁠, so |${\cal F}$| is invariant under this particular change of momentum variables. Equation (24) then becomes \begin{eqnarray*} {\rm E}^{\mu }_{\hat{a}} \, \left(\frac{p^{\hat{a}}}{m} \right) \, \frac{\partial {\cal F}}{\partial x^{\mu }} + \left(\frac{p^\mu }{m} \right) \, \frac{\partial {\rm E}^{\nu }_{\hat{a}}}{\partial x^{\mu }}\, {\rm E}_{\nu }^{\hat{b}} \, p_{\hat{b}} \, \frac{\partial {\cal F}}{\partial p_{\hat{a}}} + \frac{q}{m} {\rm E}_{\nu }^{\hat{b}} {\rm E}^{\mu }_{\hat{a}} \,{\rm \mathit{ F}}_{\mu }^{\, \nu } \, p_{\hat{b}} \, \frac{\partial {\cal F}}{\partial p_{\hat{a}}} = {\rm St}({\cal F}) . \end{eqnarray*} (28) Here we have used transformation rules (equation 2) repeatedly. The Faraday tensor in the FRF is |${\rm \mathit{ F}}^{\hat{b}}_{\hat{a}} = {\rm E}_{\nu }^{\hat{b}}{\rm E}^{\mu }_{\hat{a}} \,{\rm \mathit{ F}}_{\mu }^{\, \nu }$|⁠. The first relation in equation (4) allows us to write the factor in front of |$\partial {\cal F}/\partial p_{\hat{a}}$| in the second term of equation (28) as: \begin{eqnarray*} \left(\frac{p^\mu }{m} \right) \, \frac{\partial {\rm E}^{\nu }_{\hat{a}}}{\partial x^{\mu }} \, {\rm E}_{\nu }^{\hat{b}} \, p_{\hat{b}} = - \left(\frac{p^\mu }{m} \right) \, {\rm E}^{\nu }_{\hat{a}} \, \frac{\partial {\rm E}_{\nu }^{\hat{b}}}{\partial x^{\mu }} \, p_{\hat{b}} = - {\rm E}^{\mu }_{\hat{c}} \, {\rm E}^{\nu }_{\hat{a}} \, \frac{\partial {\rm E}_{\nu }^{\hat{b}}}{\partial x^{\mu }} \, \frac{p^{\hat{c}} p_{\hat{b}}}{m} . \end{eqnarray*} (29) A comparison with definition (9) shows immediately that this is the fictitious force |$K_{\hat{a}}$| (apart from a trivial relabelling of indices) as is expected on physical grounds. This gives the final form of the Vlasov equation, now in mixed LF/FRF variables: \begin{eqnarray*} {\rm E}^{\mu }_{\hat{a}} \, \left(\frac{p^{\hat{a}}}{m} \right) \, \frac{\partial {\cal F}}{\partial x^{\mu }} + \left(K_{\hat{a}} + \frac{q}{m} {\rm \mathit{ F}}^{\hat{b}}_{\hat{a}} \, p_{\hat{b}} \right) \, \frac{\partial {\cal F}}{\partial p_{\hat{a}}} = {\rm St}({\cal F}) . \end{eqnarray*} (30) Here |${\cal F}$| should now be interpreted as |${\cal F} \equiv {\cal F}(x^{\mu }\,, p_{\hat{a}}) = 2 m \Theta (\bar{{\cal E}}) \, f({\boldsymbol x }\,, t\,, \bar{{\boldsymbol{p}}}) \, \delta (\bar{{\cal E}}^2 - |\bar{{\boldsymbol{p}}}|^2 - m^2)$|⁠. Using the properties of the δ-function and performing an integration over FRF energy |$\bar{{\cal E}}= p_{\hat{0}}$| of the form |$\int {\rm d} \bar{{\cal E}}\, \bar{{\cal E}}\, \ldots$| of the whole equation one can express equation (30) in terms of |$f({\boldsymbol x }\,, t\,, \bar{{\boldsymbol{p}}})$| as \begin{eqnarray*} {\rm E}^{\mu }_{\hat{a}} \, \frac{p^{\hat{a}}}{m} \, \frac{\partial f}{\partial x^{\mu }} + \left\lbrace \bar{\gamma } q \left({\boldsymbol {\bar{E}} } + \bar{{\boldsymbol{V}}}{\boldsymbol{\times}}{{\boldsymbol {\bar{B} }}} \right) - (\Gamma + 1) \, \bar{{\cal E}}\, \frac{{\rm d}}{{\rm d} \lambda } \left(\frac{{\boldsymbol{U}}}{\Gamma + 1} \right) + \bar{{\boldsymbol{p}}}{\boldsymbol{\times}}{\boldsymbol \omega } \right\rbrace {\boldsymbol{\cdot}}\frac{\partial f}{\partial \bar{{\boldsymbol{p}}}} = {\rm St}(f) . \end{eqnarray*} (31) Here St(f) is the collision term acting on the ordinary Vlasov distribution |$f({\boldsymbol x }\,, t\,, \bar{{\boldsymbol{p}}})$|⁠. For pure scattering that only changes the direction of particle momenta we have we have |${\rm St}({\cal F}) = 2 m \Theta (p_{0}) \, \delta (p^2 - m^2) \, {\rm St}(f)$| (cf. equation 25). Also |$\bar{\gamma } = \bar{{\cal E}}/m$| is the Lorentz factor of a particle in the FRF with velocity |$\bar{{\boldsymbol{V}}}= \bar{{\boldsymbol{p}}}/\bar{\gamma } m$|⁠, and where |$\bar{{\cal E}}$| must be interpreted as |$\bar{{\cal E}}\equiv \sqrt{\bar{p}^2 + m^2}$|⁠. The three-vectors |$({\boldsymbol {\bar{E}} }\ ,, {\boldsymbol {\bar{B}} })$| are the electromagnetic fields in the FRF. Definition (8) for d/dλ together with the definition of |${\rm E}^{\mu }_{\hat{a}}$| allows us to write \begin{eqnarray*} \frac{{\rm d}}{{\rm d} \lambda } = {\rm E}^{\mu }_{\hat{a}} \, \frac{p^{\hat{a}}}{m} \, \frac{\partial }{\partial x^{\mu }} = \bar{\gamma } \, \left(\frac{{\rm d}}{{\rm d} \tau } + \bar{{\boldsymbol{V}}}{\boldsymbol{\cdot}}{\boldsymbol {\tilde{{\nabla}}}}\right) = \bar{\gamma } \, \left\lbrace \Gamma \, \left[ \, 1 + (\bar{{\boldsymbol{V}}}{\boldsymbol{\cdot}}{\boldsymbol{V}}) \, \right]\, \frac{\partial }{\partial t} + (\bar{{\boldsymbol{V}}}+ {\boldsymbol{U}}) {\boldsymbol{\cdot}}{\nabla}+ \frac{(\bar{{\boldsymbol{V}}}{\boldsymbol{\cdot}}{\boldsymbol{U}})}{\Gamma +1} \, ({\boldsymbol{U}} {\boldsymbol{\cdot}}{\nabla}) \right\rbrace . \end{eqnarray*} (32) Here we use |${\boldsymbol{U}} = \Gamma {\boldsymbol{V}}$| and |$\bar{{\boldsymbol{V}}}\equiv \bar{{\boldsymbol{p}}}/\bar{\gamma }m$| and introduce the notation \begin{eqnarray*} \frac{{\rm d}}{{\rm d} \tau } = \frac{\partial }{\partial x^{\hat{0}}} = U {\boldsymbol{\cdot}}\nabla = \Gamma \, \frac{\partial }{\partial t} + ({\boldsymbol{U}} {\boldsymbol{\cdot}}{\nabla}) \,, \end{eqnarray*} (33) for the comoving proper time derivative associated with the bulk plasma, and \begin{eqnarray*} \tilde{{\nabla}} = \frac{\partial }{\partial {\boldsymbol {\bar{x}} }} = {\nabla}+ {\boldsymbol{U}} \, \left(\frac{\partial }{\partial t} + \frac{({\boldsymbol{U}} {\boldsymbol{\cdot}}{\nabla})}{\Gamma +1} \right) \end{eqnarray*} (34) for the FRF gradient operator expressed in LF time and spatial derivatives. 4.1 Collision operator To describe the effect of ‘collisions’ (for instance scattering by magnetic irregularities) we choose the Anderson–Witting form for the relaxation operator (Anderson & Witting 1974) with |$p {\boldsymbol{\cdot}}U \equiv p_{\mu } U^{\mu } = \bar{{\cal E}}= \bar{\gamma }m$|⁠: \begin{eqnarray*} {\rm St}(f) \equiv \bar{\gamma } \, \overline{{\rm St}}(f) = - \frac{p {\boldsymbol{\cdot}}U}{m} \, \frac{f - f_{0}}{\bar{\tau }_{\rm r}} = - \, \frac{\bar{\gamma } \left(f - f_{0}\right)}{\bar{\tau }_{\rm r}} . \end{eqnarray*} (35) This operator is the proper relativistic version of the well-known Krook collision operator (see Bhatnagar, Gross, & Krook 1954). Here f0 is an equilibrium distribution function (usually chosen to be isotropic in the FRF) and |$\bar{\tau }_{\rm r}$| the relaxation time in the FRF. For our present purposes this is the most convenient choice, but other choices (such as a Fokker–Planck type diffusion operator in momentum that describes scattering) are possible. In Section 11 below we briefly discuss the effect of using an alternative choice for the relaxation operator: diffusive (Fokker–Planck type) scattering. If we use the Anderson–Witting relaxation operator the Vlasov equation becomes \begin{eqnarray*} \Gamma \, \left[ \, 1 + (\bar{{\boldsymbol{V}}}{\boldsymbol{\cdot}}{\boldsymbol{V}}) \, \right]\, \frac{\partial f}{\partial t} + (\bar{{\boldsymbol{V}}}+ {\boldsymbol{U}}) {\boldsymbol{\cdot}}{\nabla}f + \frac{(\bar{{\boldsymbol{V}}}{\boldsymbol{\cdot}}{\boldsymbol{U}})}{\Gamma +1} \, ({\boldsymbol{U}} {\boldsymbol{\cdot}}{\nabla}) f + \left\lbrace {\boldsymbol {\bar{K}} }_{\rm f} + \bar{{\boldsymbol{p}}}{\boldsymbol{\times}}{\boldsymbol {\bar{\omega}}} + q \left({\boldsymbol {\bar{E}} } + \bar{{\boldsymbol{V}}}{\boldsymbol{\times}}{\boldsymbol {\bar{B}} } \right) \right\rbrace {\boldsymbol{\cdot}}\frac{\partial f}{\partial \bar{{\boldsymbol{p}}}} = - \frac{f - f_{0}}{\bar{\tau }_{\rm r}} . \end{eqnarray*} (36) Here we have cancelled a common factor |$\bar{\gamma }$| and defined the fictitious three-force as |$({\rm d} \bar{{\boldsymbol{p}}}/{\rm d}t)_{\rm fict} \equiv {\boldsymbol {\bar{K}} }_{\rm f} + \bar{{\boldsymbol{p}}}{\boldsymbol{\times}}{\boldsymbol {\bar{\omega}}}$|⁠, where \begin{eqnarray*} {\boldsymbol {\bar{K}} }_{\rm f} \equiv \frac{1}{\bar{\gamma }} \, \frac{{\rm d} \bar{{\boldsymbol{p}}}}{{\rm d} \lambda } = - (\Gamma + 1) \, \frac{{\rm d}}{{\rm d} \lambda } \left(\frac{m \, {\boldsymbol{U}}}{\Gamma + 1} \right) = - \bar{\gamma } m \, (\Gamma + 1) \, \left\lbrace \Gamma \, \left[ \, 1 + (\bar{{\boldsymbol{V}}}{\boldsymbol{\cdot}}{\boldsymbol{V}}) \, \right]\, \frac{\partial }{\partial t} + (\bar{{\boldsymbol{V}}}+ {\boldsymbol{U}}) {\boldsymbol{\cdot}}{\nabla}+ \frac{(\bar{{\boldsymbol{V}}}{\boldsymbol{\cdot}}{\boldsymbol{U}})}{\Gamma +1} \, ({\boldsymbol{U}} {\boldsymbol{\cdot}}{\nabla}) \right\rbrace \left(\frac{{\boldsymbol{U}}}{\Gamma + 1} \right) \end{eqnarray*} (37) and \begin{eqnarray*} {\boldsymbol {\bar{\omega}}} \equiv \frac{{\boldsymbol \omega }}{\bar{\gamma }} = \frac{\Gamma - 1}{\bar{\gamma }} \, \left({\boldsymbol{\hat{\ell}}}{\boldsymbol{\times}}\frac{{\rm d} {\boldsymbol{\hat{\ell}}}}{{\rm d} \lambda } \right) = (\Gamma - 1) \, {\boldsymbol{\hat{\ell}}}{\boldsymbol{\times}}\left\lbrace \Gamma \, \left[ \, 1 + (\bar{{\boldsymbol{V}}}{\boldsymbol{\cdot}}{\boldsymbol{V}}) \, \right]\, \frac{\partial {\boldsymbol{\hat{\ell}}}}{\partial t} + (\bar{{\boldsymbol{V}}}+ {\boldsymbol{U}}) {\boldsymbol{\cdot}}{\nabla}{\boldsymbol{\hat{\ell}}}+ \frac{(\bar{{\boldsymbol{V}}}{\boldsymbol{\cdot}}{\boldsymbol{U}})}{\Gamma +1} \, ({\boldsymbol{U}} {\boldsymbol{\cdot}}{\nabla}) {\boldsymbol{\hat{\ell}}}\right\rbrace . \end{eqnarray*} (38) Note that our definition for |${\boldsymbol {\bar{K}} }_{\rm f}$| agrees with the usual definition of the Minkowski force in special relativity (see for instance Goldstein, Poole, & Safko 2002). Defining |${\boldsymbol {\bar{K}} } = {\boldsymbol {\bar{K}} }_{\rm f} + \bar{{\boldsymbol{p}}}{\boldsymbol{\times}}{\boldsymbol {\bar{\omega}}}$|⁠, the change of the FRF four-momentum |$p^{\hat{a}}=(\bar{\cal {E}} \,,\quad \bar{{\boldsymbol{p}}})$| can be written in standard form: \begin{eqnarray*} K^{\hat{a}} \equiv \frac{{\rm d} p^{\hat{a}}}{{\rm d} \lambda } = \bar{\gamma } \, \left(\bar{{\boldsymbol{V}}}{\boldsymbol{\cdot}}{\boldsymbol {\bar{K}} }\,, \quad {\boldsymbol {\bar{K}} } \right)= \bar{\gamma } \left(\bar{{\boldsymbol{V}}}{\boldsymbol{\cdot}}{\boldsymbol {\bar{K}} }_{\rm f}\,, \quad {\boldsymbol {\bar{K}} }_{\rm f}+ \bar{{\boldsymbol{p}}}{\boldsymbol{\times}}{\boldsymbol {\bar{\omega}}} \right) . \end{eqnarray*} (39) Here we used |$\bar{{\boldsymbol{V}}}{\boldsymbol{\cdot}}(\bar{{\boldsymbol{p}}}{\boldsymbol{\times}}{\boldsymbol {\bar{\omega}}})=0$|⁠. With definitions (33) and (34), Vlasov equation (36) and relations (37) and (38) can be written more compactly: \begin{eqnarray*} \frac{{\rm d}f}{{\rm d} \tau } + (\bar{{\boldsymbol{V}}}{\boldsymbol{\cdot}}{\boldsymbol {\tilde{{\nabla}}}})f + \left\lbrace {\boldsymbol {\bar{K}} }_{\rm f} + \bar{{\boldsymbol{p}}}{\boldsymbol{\times}}{\boldsymbol {\bar{\omega}}} + q \left({\boldsymbol {\bar{E}} } + \bar{{\boldsymbol{V}}}{\boldsymbol{\times}}{\boldsymbol {\bar{B}} } \right) \right\rbrace {\boldsymbol{\cdot}}\frac{\partial f}{\partial \bar{{\boldsymbol{p}}}} = - \frac{f - f_{0}}{\bar{\tau }_{\rm r}} \,, \end{eqnarray*} (40) with \begin{eqnarray*} {\boldsymbol {\bar{K}} }_{\rm f} = - \bar{\gamma } m \, (\Gamma + 1) \, \left(\frac{{\rm d}}{{\rm d} \tau } + \bar{{\boldsymbol{V}}}{\boldsymbol{\cdot}}{\boldsymbol {\tilde{{\nabla}}}}\right) \left(\frac{{\boldsymbol{U}}}{\Gamma + 1} \right), \quad {\boldsymbol {\bar{\omega}}} = (\Gamma - 1) \, {\boldsymbol{\hat{\ell}}}{\boldsymbol{\times}}\left(\frac{{\rm d} {\boldsymbol{\hat{\ell}}}}{{\rm d} \tau } + (\bar{{\boldsymbol{V}}}{\boldsymbol{\cdot}}{\boldsymbol {\tilde{{\nabla}}}}) {\boldsymbol{\hat{\ell}}}\right) . \end{eqnarray*} (41) Version (40) of the Vlasov equation is the equation that we will use in our further analysis, as it allows direct comparison with known (non-relativistic) results. 4.2 Liouville’s theorem By introducing the seven-dimensional state vector in phase-space |$\tilde{{\boldsymbol Z }}$| as |$\tilde{{\boldsymbol Z }} \equiv \left(x^{\mu }\,, - p_{\hat{i}} \right) = \left(t\,, {\boldsymbol x }, \bar{{\boldsymbol{p}}}\right)$|⁠, equation (40) can be written as \begin{eqnarray*} \dot{\tilde{{\boldsymbol Z }}} {\boldsymbol{\cdot}}\frac{\partial f(\tilde{{\boldsymbol Z }})}{\partial \tilde{{\boldsymbol Z }}} = \overline{{\rm St}}(f) = - \left(\frac{f - f_{0}}{\bar{\tau }_{\rm r}} \right) . \end{eqnarray*} (42) Here we define \begin{eqnarray*} \dot{\tilde{{\boldsymbol Z }}} = \frac{1}{\bar{\gamma }} \, \frac{{\rm d} \tilde{{\boldsymbol Z }}}{{\rm d} \lambda } = \left(\dot{t}\,, \quad \dot{{\boldsymbol x }}\,, \quad \dot{\bar{{\boldsymbol{p}}}} \right) \end{eqnarray*} (43) with \begin{eqnarray*} \dot{t} & = & \frac{1}{\bar{\gamma }} \, \frac{{\rm d}t}{{\rm d} \lambda } = \Gamma \, \left[ \, 1 + (\bar{{\boldsymbol{V}}}{\boldsymbol{\cdot}}{\boldsymbol{V}}) \, \right] = \frac{\gamma }{\bar{\gamma }} \,, \nonumber \\ \dot{{\boldsymbol x }} & = & \frac{1}{\bar{\gamma }} \, \frac{{\rm d} {\boldsymbol x }}{{\rm d} \lambda } = \bar{{\boldsymbol{V}}}+ {\boldsymbol{U}} + \frac{(\bar{{\boldsymbol{V}}}{\boldsymbol{\cdot}}{\boldsymbol{U}})}{\Gamma +1} \, {\boldsymbol{U}} = \frac{{\boldsymbol{p}}}{\bar{\gamma }m} \,, \nonumber \\ \dot{\bar{{\boldsymbol{p}}}} & = & \frac{1}{\bar{\gamma }} \, \frac{{\rm d} \bar{{\boldsymbol{p}}}}{{\rm d} \lambda }= {\boldsymbol {\bar{K}} }_{\rm f} + \bar{{\boldsymbol{p}}}{\boldsymbol{\times}}{\boldsymbol {\bar{\omega}}}+ q \left({\boldsymbol {\bar{E}} } + \bar{{\boldsymbol{V}}}{\boldsymbol{\times}}{\boldsymbol {\bar{B}} } \right) \equiv {\boldsymbol {\bar{K}} } + q \left({\boldsymbol {\bar{E}} } + \bar{{\boldsymbol{V}}}{\boldsymbol{\times}}{\boldsymbol {\bar{B}} } \right) . \end{eqnarray*} (44) In order to ensure particle conservation after integration over the invariant phase-space element |${\rm d}^3 \bar{{\boldsymbol{p}}}/\bar{{\cal E}}= {\rm d}^3 \bar{{\boldsymbol{p}}}/\bar{\gamma }m$| one should be able to write equation (42) in conservation form: \begin{eqnarray*} \frac{\partial }{\partial \tilde{{\boldsymbol Z }}} {\boldsymbol{\cdot}}\left(\dot{\tilde{{\boldsymbol Z }}} \, f(\tilde{{\boldsymbol Z }}) \right) = - \left(\frac{f - f_{0}}{\bar{\tau }_{\rm r}} \right) . \end{eqnarray*} (45) This implies that the flow in phase space is incompressible (Liouville’s theorem): \begin{eqnarray*} \frac{\partial }{\partial \tilde{{\boldsymbol Z }}} {\boldsymbol{\cdot}}\dot{\tilde{{\boldsymbol Z }}} \equiv {\rm div}_{Z} \dot{\tilde{{\boldsymbol Z }}} = 0 . \end{eqnarray*} (46) For the original (covariant) equation (24) proving the equivalence of the two formulations is simple. The equations of motion in phase space can be derived from a Hamiltonian |${\cal H}(x\,, p)$| and the Vlasov equation has a symplectic structure, meaning that the transition from the non-conservative form to the conservation form generates a term that vanishes trivially: \begin{eqnarray*} - \left(\frac{\partial^2 {\cal H}}{\partial x^{\mu } \, \partial p_{\mu }} - \frac{\partial^2 {\cal H}}{\partial p_{\mu } \, \partial x^{\mu }} \right) \, {\cal F}(x\,, p) = 0 \end{eqnarray*} (47) When one uses mixed phase-space coordinates |$\tilde{{\boldsymbol Z }}= \left(t\,, {\boldsymbol x }\,, \bar{{\boldsymbol{p}}}\right)$|⁠, the proof is more difficult. It is therefore shown in detail. 4.2.1 Liouville’s theorem for mixed coordinates The space–time part of |${\rm div}_{Z} \dot{\tilde{{\boldsymbol Z }}}$| is \begin{eqnarray*} \left({\rm div}_{Z} \dot{\tilde{{\boldsymbol Z }}} \right)_{\rm st} = \frac{\partial \dot{x}^{\mu }}{\partial x^{\mu }} = \frac{\partial \Gamma }{\partial t} + {\nabla}{\boldsymbol{\cdot}}{\boldsymbol{U}} + \left(\frac{\partial }{\partial t} + \frac{{\boldsymbol{U}} {\boldsymbol{\cdot}}{\nabla}}{\Gamma +1} \right)\, \left(\bar{{\boldsymbol{V}}}{\boldsymbol{\cdot}}{\boldsymbol{U}} \right) + \left(\bar{{\boldsymbol{V}}}{\boldsymbol{\cdot}}{\boldsymbol{U}} \right) \, \left[ \, {\nabla}{\boldsymbol{\cdot}}\left(\frac{{\boldsymbol{U}}}{\Gamma +1} \right)\, \right]. \end{eqnarray*} (48) Using the definition for |${\boldsymbol {\tilde{{\nabla}}}}$|⁠, the relation \begin{eqnarray*} \frac{\partial \Gamma }{\partial t} + {\nabla}{\boldsymbol{\cdot}}{\boldsymbol{U}} = \left(\Gamma + 1 \right) {\boldsymbol {\tilde{{\nabla}}}}{\boldsymbol{\cdot}}\left(\frac{{\boldsymbol{U}}}{\Gamma +1} \right) \end{eqnarray*} (49) and the relation valid for an infinitesimal change |$\delta {\boldsymbol{U}}$| \begin{eqnarray*} {\boldsymbol{U}} {\boldsymbol{\cdot}}\delta \left(\frac{{\boldsymbol{U}}}{\Gamma +1} \right) = \frac{\delta \Gamma }{\Gamma +1} \ ,, \end{eqnarray*} (50) one finds that one can write: \begin{eqnarray*} \left({\rm div}_{Z} \dot{\tilde{{\boldsymbol Z }}} \right)_{\rm st} = \left(1 + \frac{(\bar{{\boldsymbol{V}}}{\boldsymbol{\cdot}}{\boldsymbol{U}})}{\Gamma +1} \right) \, \left\lbrace \left(\Gamma + 1 \right) {\boldsymbol {\tilde{{\nabla}}}}{\boldsymbol{\cdot}}\left(\frac{{\boldsymbol{U}}}{\Gamma +1} \right)\right\rbrace + (\Gamma + 1) \, \left(\frac{\partial }{\partial t} + \frac{{\boldsymbol{U}} {\boldsymbol{\cdot}}{\nabla}}{\Gamma +1} \right)\left(\frac{(\bar{{\boldsymbol{V}}}{\boldsymbol{\cdot}}{\boldsymbol{U}})}{\Gamma +1} \right). \end{eqnarray*} (51) The momentum part of |${\rm div}_{Z} \dot{\tilde{{\boldsymbol Z }}}$| is (in an obvious notation) \begin{eqnarray*} \left({\rm div}_{Z} \dot{\tilde{{\boldsymbol Z }}} \right)_{\rm mom} = \frac{\partial \dot{p}_{\hat{i}}}{\partial p_{\hat{i}}} = \frac{\partial }{\partial \bar{{\boldsymbol{p}}}} {\boldsymbol{\cdot}}\left({\boldsymbol {\bar{K}} }_{\rm f} + \bar{{\boldsymbol{p}}}{\boldsymbol{\times}}{\boldsymbol {\bar{\omega}}} \right) . \end{eqnarray*} (52) The term in |$\dot{\bar{{\boldsymbol{p}}}}$| involving the Lorentz force does not contribute: It’s momentum divergence vanishes trivially. We write |${\boldsymbol {\bar{\omega}}}$| as \begin{eqnarray*} {\boldsymbol {\bar{\omega}}} = {\boldsymbol{U}} {\boldsymbol{\times}}\left(\frac{{\rm d}}{{\rm d} \tau } + \bar{{\boldsymbol{V}}}{\boldsymbol{\cdot}}{\boldsymbol {\tilde{{\nabla}}}}\right) \left(\frac{{\boldsymbol{U}}}{\Gamma +1} \right)= \left(\frac{{\boldsymbol{U}}}{\bar{\gamma }} \right) {\boldsymbol{\times}}\frac{{\rm d}}{{\rm d} \lambda } \left(\frac{{\boldsymbol{U}}}{\Gamma +1} \right) \end{eqnarray*} (53) and expand the double cross-product: \begin{eqnarray*} \bar{{\boldsymbol{p}}}{\boldsymbol{\times}}{\boldsymbol {\bar{\omega}}} = m \, \left[ \bar{{\boldsymbol{V}}}{\boldsymbol{\cdot}}\left(\frac{{\rm d}}{{\rm d} \lambda } \left(\frac{{\boldsymbol{U}}}{\Gamma +1} \right)\right) \, \right] {\boldsymbol{U}} - m \, \left(\bar{{\boldsymbol{V}}}{\boldsymbol{\cdot}}{\boldsymbol{U}} \right) \, \frac{{\rm d}}{{\rm d} \lambda } \left(\frac{{\boldsymbol{U}}}{\Gamma +1} \right). \end{eqnarray*} (54) One then finds \begin{eqnarray*} \left({\rm div}_{Z} \dot{\tilde{{\boldsymbol Z }}} \right)_{\rm mom} & = & - \left(1 + \frac{(\bar{{\boldsymbol{V}}}{\boldsymbol{\cdot}}{\boldsymbol{U}})}{\Gamma +1} \right) \, \left\lbrace \left(\Gamma + 1 \right) {\boldsymbol {\tilde{{\nabla}}}}{\boldsymbol{\cdot}}\left(\frac{{\boldsymbol{U}}}{\Gamma +1} \right)\right\rbrace + \bar{{\boldsymbol{V}}}\, {\boldsymbol{U}} {{\boldsymbol \, : \,}}{\boldsymbol {\tilde{{\nabla}}}}\left(\frac{{\boldsymbol{U}}}{\Gamma +1} \right)- (\Gamma + 1) \, \bar{{\boldsymbol{V}}}{\boldsymbol{\cdot}}\frac{{\rm d}}{{\rm d} \tau } \left(\frac{{\boldsymbol{U}}}{\Gamma +1} \right). \end{eqnarray*} (55) The first term in this expression cancels the corresponding term in the space–time part |$\left({\rm div}_{Z} \dot{\tilde{{\boldsymbol Z }}} \right)_{\rm st}$|⁠. For the two remaining terms we can proceed as follows. Use the definitions |${\rm d}/{\rm d}\tau = \Gamma \, \partial /\partial t + {\boldsymbol{U}} {\boldsymbol{\cdot}}{\nabla}$|⁠, |${\rm d}/{\rm d} \lambda = \bar{\gamma } \, {\rm d}/{\rm d} \tau + \bar{\gamma } \, \bar{{\boldsymbol{V}}}{\boldsymbol{\cdot}}{\boldsymbol {\tilde{{\nabla}}}}$| together with the fact that all derivatives are with |$\bar{{\boldsymbol{p}}}$| (and therefore |$\bar{{\boldsymbol{V}}}$|⁠) kept constant. This allows us to write \begin{eqnarray*} \bar{{\boldsymbol{V}}}\, {\boldsymbol{U}} {{\boldsymbol \, : \,}}{\boldsymbol {\tilde{{\nabla}}}}\left(\frac{{\boldsymbol{U}}}{\Gamma +1} \right)= {\boldsymbol{U}} {\boldsymbol{\cdot}}{\boldsymbol {\tilde{{\nabla}}}}(\bar{{\boldsymbol{V}}}{\boldsymbol{\cdot}}{\boldsymbol{U}})\,, \quad \bar{{\boldsymbol{V}}}{\boldsymbol{\cdot}}\frac{{\rm d}}{{\rm d} \tau } \left(\frac{{\boldsymbol{U}}}{\Gamma +1} \right)= \frac{{\rm d}}{{\rm d} \tau } \, \left(\frac{(\bar{{\boldsymbol{V}}}{\boldsymbol{\cdot}}{\boldsymbol{U}})}{\Gamma +1} \right) . \end{eqnarray*} (56) The definition of |${\boldsymbol {\tilde{{\nabla}}}}$| allows one to combine the second and third term in expression (55) and manipulate them into \begin{eqnarray*} \bar{{\boldsymbol{V}}}\, {\boldsymbol{U}} {{\boldsymbol \, : \,}}{\boldsymbol {\tilde{{\nabla}}}}\left(\frac{{\boldsymbol{U}}}{\Gamma +1} \right)- (\Gamma + 1) \, \bar{{\boldsymbol{V}}}{\boldsymbol{\cdot}}\frac{{\rm d}}{{\rm d} \tau } \left(\frac{{\boldsymbol{U}}}{\Gamma +1} \right)= - (\Gamma + 1) \, \left(\frac{\partial }{\partial t} + \frac{{\boldsymbol{U}} {\boldsymbol{\cdot}}{\nabla}}{\Gamma +1} \right)\left(\frac{(\bar{{\boldsymbol{V}}}{\boldsymbol{\cdot}}{\boldsymbol{U}})}{\Gamma +1} \right) . \end{eqnarray*} (57) Substituting this into equation (55) quickly establishes that \begin{eqnarray*} \left({\rm div}_{Z} \dot{\tilde{{\boldsymbol Z }}} \right)_{\rm mom} = - \left({\rm div}_{Z} \dot{\tilde{{\boldsymbol Z }}} \right)_{\rm st} \, \Longleftrightarrow \, {\rm div}_{Z} \dot{\tilde{{\boldsymbol Z }}} = 0 \end{eqnarray*} (58) and Liouville’s theorem holds for the Vlasov equation in mixed variables. The fact that the Vlasov equation in mixed variables can be written in conservation form will prove useful later, in particular in Paper 2. 5 GYRO-AVERAGED EQUATION In order to compare our results with the results obtained for a non-relativistic bulk flow by Skilling (1975; see also Jones 1990) and Kulsrud (1983) we first consider the gyro-averaged Vlasov equation. The basic idea is that the magnetic term in the Lorentz force (due to a large-scale regular magnetic field) is dominant. If the field equals |${\boldsymbol {\bar{B}} }_{0}$| in the FRF, the corresponding term in the Lorentz force is, defining the gyro-frequency |$\bar{\Omega } = q \bar{B}_{0}/\bar{\gamma } m$|⁠, including sign of charge, and the gyro-angle |$\bar{\phi }$| , \begin{eqnarray*} \left(\bar{{\boldsymbol{V}}}{\boldsymbol{\times}}{\boldsymbol {\bar{B}} }_{0} \right) {\boldsymbol{\cdot}}\frac{\partial f}{\partial \bar{{\boldsymbol{p}}}} = - \bar{\Omega } \, \frac{\partial f}{\partial \bar{\phi }} . \end{eqnarray*} (59) If the gyro-period |$P_{\rm g} = 2 {\pi} /|\bar{\Omega }|$| is much shorter than all other relevant time-scales and the gyro-radius |$r_{\rm g} \sim \bar{v}/|\bar{\Omega }|$| is much smaller than all relevant length-scales, one may assume that the distribution is gyrotropic with |$\partial \langle f\rangle /\partial \bar{\phi }$| vanishingly small to lowest order in |$1/\bar{\Omega }$|⁠. One can then average the Vlasov equation over gyro-angle. The averaging is most easily performed by erecting a Cartesian coordinate system |$(x\,, y\,, z)$| with unit base vectors |${{\boldsymbol {\hat{x}}}}$|⁠, |${{\boldsymbol {\hat{y}}}}$|⁠, and |${{\boldsymbol {\hat{z}}}}$| in the FRF, choosing |${\boldsymbol {\bar{B}} }_{0}$| locally along that z-axis and writing the FRF three-momentum as \begin{eqnarray*} \bar{{\boldsymbol{p}}}= \bar{p}_{\perp }\, \left(\cos \bar{\phi } \, {{\boldsymbol {\hat{x}}}}+ \sin \bar{\phi } \, {{\boldsymbol {\hat{y}}}}\right) + \bar{p}_{\parallel }\, {{\boldsymbol {\hat{z}}}}. \end{eqnarray*} (60) The momentum pitch angle |$\bar{\theta }$| is defined in the usual manner by |$\bar{{\boldsymbol{p}}}{\boldsymbol{\cdot}}\mbox{$\boldsymbol{\hat{b}}$}= \bar{p}_{\parallel }= \bar{p}\, \bar{\mu }$|⁠, with |$\bar{\mu }= \cos \bar{\theta }$|⁠, so that |$\bar{p}_{\perp }= \bar{p}\, \sin {\bar{\theta }} =\bar{p}\, \sqrt{1 - \bar{\mu }^2}$|⁠. The gyro-averaged distribution is \begin{eqnarray*} \bar{f}({\boldsymbol x }\,, t\,, \bar{p}, \bar{\mu }) \equiv \frac{1}{2 {\pi} } \int_{0}^{2 {\pi} } {\rm d} \bar{\phi } \, f({\boldsymbol x }\,, t\,, \bar{{\boldsymbol{p}}}) . \end{eqnarray*} (61) The variable change, from |$\bar{{\boldsymbol{p}}}$| to |$\bar{p}$|⁠, |$\bar{\mu }$|⁠, and |$\bar{\phi }$|⁠, implies that one has to take account of the variation of |$\mbox{$\boldsymbol{\hat{b}}$}= {\boldsymbol {\bar{B}} }_{0}/\bar{B}_{0}$|⁠, the unit vector along the FRF magnetic field. A change in |$\mbox{$\boldsymbol{\hat{b}}$}$| affects the definition of the pitch angle cosine |$\bar{\mu }$|⁠. In particular: \begin{eqnarray*} \left. \frac{\partial }{\partial t} \right)_{{\bar{{\boldsymbol{p}}}}} \, \Longrightarrow \, \left. \frac{\partial }{\partial t} \right)_{\bar{p},\mu } + \left(\frac{\bar{{\boldsymbol{p}}}}{\bar{p}} {\boldsymbol{\cdot}}\frac{\partial \mbox{$\boldsymbol{\hat{b}}$}}{\partial t} \right) \, \frac{\partial }{\partial \bar{\mu }}, \left. {\nabla}\right)_{{\bar{{\boldsymbol{p}}}}} \, \Longrightarrow \, \left. {\nabla}\right)_{\bar{p},\mu } + \frac{ (\bar{{\boldsymbol{p}}}{\boldsymbol{\cdot}}{\nabla}) \mbox{$\boldsymbol{\hat{b}}$}}{\bar{p}} \, \frac{\partial }{\partial \bar{\mu }} . \end{eqnarray*} (62) The gyro-averaged Vlasov equation for |$\bar{f}({\boldsymbol x }\,, t\,, \bar{p}, \bar{\mu })$| reads: \begin{eqnarray*} \frac{{\rm d} \bar{f}}{{\rm d} \tau } + \bar{v}\bar{\mu } (\mbox{$\boldsymbol{\hat{b}}$}{\boldsymbol{\cdot}}{\boldsymbol {\tilde{{\nabla}}}}) \bar{f} + \dot{\bar{\mu }} \, \frac{\partial \bar{f}}{\partial \bar{\mu }} + \dot{\bar{p}} \, \frac{\partial \bar{f}}{\partial \bar{p}} = - \frac{\bar{f} - f_{0}}{\bar{\tau }_{\rm r}} . \end{eqnarray*} (63) Here \begin{eqnarray*} \dot{\bar{\mu }} & \equiv & \frac{1 - \bar{\mu }^2}{2} \, \left\lbrace \bar{v}({\boldsymbol {\tilde{{\nabla}}}}{\boldsymbol{\cdot}}\mbox{$\boldsymbol{\hat{b}}$}) - \frac{2(\Gamma + 1)}{\bar{v}} \, \left[ \, \mbox{$\boldsymbol{\hat{b}}$}{\boldsymbol{\cdot}}\frac{{\rm d}}{{\rm d}\tau } \left(\frac{{\boldsymbol{U}}}{\Gamma +1} \right)\, \right] -(\Gamma +1) \, \bar{\mu } \, \left[ \, 3 \mbox{$\boldsymbol{\hat{b}}$}\, \mbox{$\boldsymbol{\hat{b}}$}{{\boldsymbol \, : \,}}{\boldsymbol {\tilde{{\nabla}}}}\left(\frac{{\boldsymbol{U}}}{\Gamma +1} \right) - {\boldsymbol {\tilde{{\nabla}}}}{\boldsymbol{\cdot}}\left(\frac{{\boldsymbol{U}}}{\Gamma +1} \right) \, \right] \right. \nonumber \\ && \left. +\, (\Gamma - 1) \, \bar{v}\, \left[ \, (\mbox{$\boldsymbol{\hat{b}}$}{\boldsymbol{\cdot}}{\boldsymbol{\hat{\ell}}})({\boldsymbol {\tilde{{\nabla}}}}{\boldsymbol{\cdot}}{\boldsymbol{\hat{\ell}}}) - \mbox{$\boldsymbol{\hat{b}}$}\, {\boldsymbol{\hat{\ell}}}{{\boldsymbol \, : \,}}{\boldsymbol {\tilde{{\nabla}}}}{\boldsymbol{\hat{\ell}}}\, \right] \right\rbrace \end{eqnarray*} (64) and \begin{eqnarray*} \dot{\bar{p}} \equiv (\Gamma +1) \, \left\lbrace - \frac{\bar{\mu }}{\bar{v}} \, \left[ \, \mbox{$\boldsymbol{\hat{b}}$}{\boldsymbol{\cdot}}\frac{{\rm d}}{{\rm d} \tau } \left(\frac{{\boldsymbol{U}}}{\Gamma +1} \right) \, \right] + \frac{1 - 3 \bar{\mu }^2}{2} \, \left[ \, \mbox{$\boldsymbol{\hat{b}}$}\, \mbox{$\boldsymbol{\hat{b}}$}{{\boldsymbol \, : \,}}{\boldsymbol {\tilde{{\nabla}}}}\left(\frac{{\boldsymbol{U}}}{\Gamma +1} \right) - \frac{1}{3}\, {\boldsymbol {\tilde{{\nabla}}}}{\boldsymbol{\cdot}}\left(\frac{{\boldsymbol{U}}}{\Gamma +1} \right) \, \right] - \frac{1}{3}\, \left[ \, {\boldsymbol {\tilde{{\nabla}}}}{\boldsymbol{\cdot}}\left(\frac{{\boldsymbol{U}}}{\Gamma +1} \right) \, \right] \right\rbrace \, \bar{p}. \end{eqnarray*} (65) We employ the notation |$\mbox{$\boldsymbol{\hat{b}}$}\, \mbox{$\boldsymbol{\hat{b}}$}{{\boldsymbol \, : \,}}{\boldsymbol {\tilde{{\nabla}}}}{\boldsymbol A } = \hat{b}^{i} \hat{b}^{j} \tilde{\nabla }_{j} A^{i}$| for a double contraction between three-vectors and the gradient operator. We assumed for simplicity that there are no electromagnetic fields save for the large-scale magnetic field |${\boldsymbol {\bar{B}} }_{0}$| in the FRF. This corresponds to the limit of magnetohydrodynamics (see below). This result generalizes the result of Skilling (1975) to a relativistic flow. His result can be obtained by taking the limit |$|{\boldsymbol{U}}| \simeq |{\boldsymbol{V}}| \ll 1$|⁠, Γ ≃ 1 and neglecting all terms of order |$|{\boldsymbol{V}}|^2$| and higher. This implies \begin{eqnarray*} {\boldsymbol {\tilde{{\nabla}}}}\simeq {\nabla}+ {\boldsymbol{V}} \, \frac{\partial }{\partial t}\,, \quad \frac{{\rm d}}{{\rm d} \tau } \simeq \frac{\partial }{\partial t} + {\boldsymbol{V}} {\boldsymbol{\cdot}}{\nabla}\equiv \frac{{\rm d}}{{\rm d}t} . \end{eqnarray*} (66) Using these approximations in the expressions (64) and (65) and substituting the results into (63) leads to \begin{eqnarray*} & & \left\lbrace 1 + \bar{v}\bar{\mu } \, \left(\mbox{$\boldsymbol{\hat{b}}$}{\boldsymbol{\cdot}}{\boldsymbol{V}} \right) \right\rbrace \, \frac{\partial \bar{f}}{\partial t} + \left({\boldsymbol{V}} + \bar{v}\bar{\mu } \, \mbox{$\boldsymbol{\hat{b}}$}\right) {\boldsymbol{\cdot}}{\nabla}\bar{f} \nonumber \\ & & + \left\lbrace \bar{v}\, \left({\nabla}{\boldsymbol{\cdot}}\mbox{$\boldsymbol{\hat{b}}$}+ {\boldsymbol{V}} {\boldsymbol{\cdot}}\frac{\partial \mbox{$\boldsymbol{\hat{b}}$}}{\partial t} \right) + \bar{\mu } \, ({\nabla}{\boldsymbol{\cdot}}{\boldsymbol{V}}) + \frac{1 - 3 \bar{\mu }^2}{2} \, \left(\mbox{$\boldsymbol{\hat{b}}$}\mbox{$\boldsymbol{\hat{b}}$}{{\boldsymbol \, : \,}}{\nabla}{\boldsymbol{V}} \right) - \frac{2}{\bar{v}} \, \left(\mbox{$\boldsymbol{\hat{b}}$}{\boldsymbol{\cdot}}\frac{{\rm d} {\boldsymbol{V}}}{{\rm d}t } \right) \right\rbrace \frac{1 - \bar{\mu }^2}{2} \,\frac{\partial \bar{f}}{\partial \bar{\mu }} \nonumber \\ & & + \left\lbrace \frac{1 - 3 \bar{\mu }^2}{2} \, \left(\mbox{$\boldsymbol{\hat{b}}$}\mbox{$\boldsymbol{\hat{b}}$}{{\boldsymbol \, : \,}}{\nabla}{\boldsymbol{V}} \right) - \frac{1 - \bar{\mu }^2}{2} \, \left({\nabla}{\boldsymbol{\cdot}}{\boldsymbol{V}} \right) - \frac{\bar{\mu }}{\bar{v}} \, \left(\mbox{$\boldsymbol{\hat{b}}$}{\boldsymbol{\cdot}}\frac{{\rm d} {\boldsymbol{V}}}{{\rm d}t } \right) \right\rbrace \, \bar{p}\frac{\partial \bar{f}}{\partial \bar{p}} = - \frac{\bar{f} - f_{0}}{\bar{\tau }_{\rm r}} . \end{eqnarray*} (67) This agrees exactly with Skilling’s result (equation 4 in Skilling 1975), except for a different choice for the scattering operator on the right-hand side of this equation. In the term in front of |$\bar{p}\, (\partial \bar{f}/\partial \bar{p})$| Skilling makes the additional approximation |${\rm d} {\boldsymbol{V}}/{\rm d}t \simeq \partial {\boldsymbol{V}}/\partial t$|⁠. While this is formally correct since the convective |$({\boldsymbol{V}} {\boldsymbol{\cdot}}{\nabla}) {\boldsymbol{V}}$| term in |${\rm d} {\boldsymbol{V}}/{\rm d}t \equiv \partial {\boldsymbol{V}}/\partial t + ({\boldsymbol{V}} {\boldsymbol{\cdot}}{\nabla}) {\boldsymbol{V}}$| is of order |$|{\boldsymbol{V}}|^2$|⁠, both terms in |${\rm d} {\boldsymbol{V}}/{\rm d}t$| are usually of similar magnitude in realistic (magneto)hydrodynamic flows. Therefore, in practical applications both should be retained (or discarded). 6 RELATIVISTIC DRIFT KINETIC EQUATION Kulsrud (1983) employs |$\bar{p}_{\perp }$| and |$\bar{p}_{\parallel }$| as fundamental momentum variables. In addition he allows for a weak quasi-static electric field |$\bar{E}_{\parallel }$| along the magnetic field in the FRF, corresponding to a small deviation from the MHD approximation for the bulk flow. The gyro-averaged equation for |$\bar{f}(\bar{p}_{\perp }, \bar{p}_{\perp }\,, {\boldsymbol x }\,, t)$| in these variables formally reads (with |$\bar{v}_{\parallel } = \bar{v}\bar{\mu }$|⁠) \begin{eqnarray*} \frac{{\rm d} \bar{f}}{{\rm d} \tau } + \bar{v}_{\parallel } \, (\mbox{$\boldsymbol{\hat{b}}$}{\boldsymbol{\cdot}}{\boldsymbol {\tilde{{\nabla}}}}) \bar{f} + \dot{\bar{p}}_{\perp } \, \frac{\partial \bar{f}}{\partial \bar{p}_{\perp }} + \left(\dot{\bar{p}}_{\parallel } + q \bar{E}_{\parallel } \right) \, \frac{\partial \bar{f}}{\partial \bar{p}_{\parallel }} = - \frac{\bar{f} - f_{0}}{\bar{\tau }_{\rm r}} . \end{eqnarray*} (68) We have separated out the Coulomb force |$q \bar{E}_{\parallel }$| in the parallel force term for clarity. The quantities |$\dot{\bar{p}}_{\perp }$| and |$\dot{\bar{p}}_{\parallel }$| follow from equations (64) and (65) as \begin{eqnarray*} \dot{\bar{p}}_{\perp } & = & \sqrt{1 - \bar{\mu }^2}\, \dot{\bar{p}} - \frac{\bar{p}\bar{\mu } \, \dot{\bar{\mu }}}{\sqrt{1 - \bar{\mu }^2}} = (\Gamma + 1) \, \frac{\bar{p}_{\perp }}{2} \, \left\lbrace \mbox{$\boldsymbol{\hat{b}}$}\, \mbox{$\boldsymbol{\hat{b}}$}{{\boldsymbol \, : \,}}{\boldsymbol {\tilde{{\nabla}}}}\left(\frac{{\boldsymbol{U}}}{\Gamma +1} \right) - {\boldsymbol {\tilde{{\nabla}}}}{\boldsymbol{\cdot}}\left(\frac{{\boldsymbol{U}}}{\Gamma +1} \right) \right\rbrace \nonumber \\ && -\, \frac{\bar{p}_{\perp } \, \bar{v}_{\parallel }}{2} \, \left\lbrace \left({\boldsymbol {\tilde{{\nabla}}}}{\boldsymbol{\cdot}}\mbox{$\boldsymbol{\hat{b}}$}\right) + (\Gamma - 1) \, \left[ \, (\mbox{$\boldsymbol{\hat{b}}$}{\boldsymbol{\cdot}}{\boldsymbol{\hat{\ell}}})({\boldsymbol {\tilde{{\nabla}}}}{\boldsymbol{\cdot}}{\boldsymbol{\hat{\ell}}}) - \mbox{$\boldsymbol{\hat{b}}$}\, {\boldsymbol{\hat{\ell}}}{{\boldsymbol \, : \,}}{\boldsymbol {\tilde{{\nabla}}}}{\boldsymbol{\hat{\ell}}}\, \right] \right\rbrace . \end{eqnarray*} (69) and \begin{eqnarray*} \dot{\bar{p}}_{\parallel } & \equiv & \bar{\mu } \, \dot{\bar{p}} + \bar{p}\, \dot{\bar{\mu }} = - (\Gamma +1) \, \bar{\gamma } m \, \left\lbrace \, \mbox{$\boldsymbol{\hat{b}}$}{\boldsymbol{\cdot}}\left[ \, \frac{{\rm d}}{{\rm d} \tau } \left(\frac{{\boldsymbol{U}}}{\Gamma +1} \right) + \bar{v}_{\parallel } \, (\mbox{$\boldsymbol{\hat{b}}$}{\boldsymbol{\cdot}}{\boldsymbol {\tilde{{\nabla}}}}) \left(\frac{{\boldsymbol{U}}}{\Gamma +1} \right) \, \right] \right\rbrace \nonumber \\ &&+\, \frac{\bar{p}_{\perp } \, \bar{v}_{\perp }}{2} \, \left\lbrace \left({\boldsymbol {\tilde{{\nabla}}}}{\boldsymbol{\cdot}}\mbox{$\boldsymbol{\hat{b}}$}\right) + (\Gamma - 1) \, \left[ \, (\mbox{$\boldsymbol{\hat{b}}$}{\boldsymbol{\cdot}}{\boldsymbol{\hat{\ell}}})({\boldsymbol {\tilde{{\nabla}}}}{\boldsymbol{\cdot}}{\boldsymbol{\hat{\ell}}}) - \mbox{$\boldsymbol{\hat{b}}$}\, {\boldsymbol{\hat{\ell}}}{{\boldsymbol \, : \,}}{\boldsymbol {\tilde{{\nabla}}}}{\boldsymbol{\hat{\ell}}}\, \right] \right\rbrace . \end{eqnarray*} (70) These results can be simplified if we assume that – to lowest order – relativistic MHD applies when describing the large-scale magnetic field. In that case one can derive an equation for the evolution of |${{\rm B}}= |\bar{{\boldsymbol B }}_{0}|$|⁠, the strength of the magnetic field in the FRF. This is in fact a Lorentz-invariant quantity that follows from |${\rm F}^{\mu \nu } \, {\rm F}_{\mu \nu } = 2(| {\boldsymbol E } |^2 - | {\boldsymbol B } |^2) \simeq - 2 {{\rm B}}^2$| as long as |$|\bar{E}_{\parallel }| \ll {{\rm B}}$| as is assumed here. As shown in Appendix C we have \begin{eqnarray*} \frac{1}{{{\rm B}}} \, \frac{{\rm d} {{\rm B}}}{{\rm d} \tau } = (\Gamma +1) \, \left\lbrace \mbox{$\boldsymbol{\hat{b}}$}\, \mbox{$\boldsymbol{\hat{b}}$}{{\boldsymbol \, : \,}}{\boldsymbol {\tilde{{\nabla}}}}\left(\frac{{\boldsymbol{U}} }{\Gamma +1} \right) - {\boldsymbol {\tilde{{\nabla}}}}{\boldsymbol{\cdot}}\left(\frac{{\boldsymbol{U}} }{\Gamma +1} \right) \right\rbrace \,, \end{eqnarray*} (71) and \begin{eqnarray*} {\boldsymbol {\tilde{{\nabla}}}}{\boldsymbol{\cdot}}\mbox{$\boldsymbol{\hat{b}}$}= - \frac{(\mbox{$\boldsymbol{\hat{b}}$}{\boldsymbol{\cdot}}{\boldsymbol {\tilde{{\nabla}}}}) {{\rm B}}}{{{\rm B}}} + (\Gamma + 1) \left\lbrace \mbox{$\boldsymbol{\hat{b}}$}{\boldsymbol{\cdot}}\frac{{\rm d}}{{\rm d}\tau } \left(\frac{{\boldsymbol{U}}}{\Gamma +1} \right) \right\rbrace . \end{eqnarray*} (72) The second term in this last expression is a relativistic correction that occurs when the FRF is not an inertial frame. This allows us to write \begin{eqnarray*} \dot{\bar{p}}_{\perp } = \frac{\bar{p}_{\perp }}{2 {{\rm B}}} \, \left(\frac{{\rm d} {{\rm B}}}{{\rm d} \tau } + \bar{v}_{\parallel } \, (\mbox{$\boldsymbol{\hat{b}}$}{\boldsymbol{\cdot}}{\boldsymbol {\tilde{{\nabla}}}}) {{\rm B}}\right) - \nu_{\rm p} \, \frac{\bar{p}_{\perp } \, \bar{v}_{\parallel }}{2} \,, \end{eqnarray*} (73) and \begin{eqnarray*} \dot{\bar{p}}_{\parallel } = - \frac{\bar{p}_{\perp } \, \bar{v}_{\perp }}{2 {{\rm B}}} \, (\mbox{$\boldsymbol{\hat{b}}$}{\boldsymbol{\cdot}}{\boldsymbol {\tilde{{\nabla}}}}) {{\rm B}}- (\Gamma +1) \, \bar{\gamma } m \, \mbox{$\boldsymbol{\hat{b}}$}{\boldsymbol{\cdot}}\left\lbrace \frac{{\rm d}}{{\rm d} \tau } \left(\frac{{\boldsymbol{U}}}{\Gamma +1} \right) + \bar{v}_{\parallel } \, (\mbox{$\boldsymbol{\hat{b}}$}{\boldsymbol{\cdot}}{\boldsymbol {\tilde{{\nabla}}}}) \left(\frac{{\boldsymbol{U}}}{\Gamma +1} \right) \right\rbrace + \nu_{\rm p} \, \frac{\bar{p}_{\perp } \, \bar{v}_{\perp }}{2} . \end{eqnarray*} (74) Here νp is a frequency that is essentially a relativistic correction, given by \begin{eqnarray*} \nu_{\rm p} \equiv (\Gamma + 1) \left[ \, \mbox{$\boldsymbol{\hat{b}}$}{\boldsymbol{\cdot}}\frac{{\rm d}}{{\rm d} \tau } \left(\frac{{\boldsymbol{U}}}{\Gamma +1} \right) \right] + (\Gamma - 1) \, \left[ \, (\mbox{$\boldsymbol{\hat{b}}$}{\boldsymbol{\cdot}}{\boldsymbol{\hat{\ell}}})({\boldsymbol {\tilde{{\nabla}}}}{\boldsymbol{\cdot}}{\boldsymbol{\hat{\ell}}}) - \mbox{$\boldsymbol{\hat{b}}$}\, {\boldsymbol{\hat{\ell}}}{{\boldsymbol \, : \,}}{\boldsymbol {\tilde{{\nabla}}}}{\boldsymbol{\hat{\ell}}}\, \right]. \end{eqnarray*} (75) 6.1 Non-relativistic limit and conservation of magnetic moment If we now take the limit of a non-relativistic flow and assume non-relativistic particles so that |$\bar{v}\ll 1$|⁠, we can make the replacements (equation 66) and in addition neglect the terms involving νp. In that case equation (68) takes the form derived by Kulsrud (1983), in his equation (51), after substituting the non-relativistic equivalent of (72), \begin{eqnarray*} {\nabla}{\boldsymbol{\cdot}}\mbox{$\boldsymbol{\hat{b}}$}= - \frac{(\mbox{$\boldsymbol{\hat{b}}$}{\boldsymbol{\cdot}}{\nabla}) {{\rm B}}}{{{\rm B}}} \,, \end{eqnarray*} (76) and adding the scattering term. One finds \begin{eqnarray*} \frac{{\rm D} \bar{f}}{{\rm D}t} + \left(\frac{\bar{p}_{\perp }}{2 {{\rm B}}} \, \frac{{\rm D} {{\rm B}}}{{\rm D} t} \right) \, \frac{\partial \bar{f}}{\partial \bar{p}_{\perp }} + \left(q \, \bar{E}_{\parallel } - \frac{\bar{p}_{\perp } \, \bar{v}_{\perp }}{2 {{\rm B}}} \, (\mbox{$\boldsymbol{\hat{b}}$}{\boldsymbol{\cdot}}{\nabla}) {{\rm B}}- m \, \mbox{$\boldsymbol{\hat{b}}$}{\boldsymbol{\cdot}}\frac{{\rm D} {\boldsymbol{V}}}{{\rm D} t} \right) \, \frac{\partial \bar{f}}{\partial \bar{p}_{\parallel }} = - \frac{\bar{f} - f_{0}}{\bar{\tau }_{\rm r}} . \end{eqnarray*} (77) Here the operator \begin{eqnarray*} \frac{{\rm D}}{{\rm D}t} \equiv \frac{\partial }{\partial t} + \left({\boldsymbol{V}} + \bar{v}_{\parallel } \, \mbox{$\boldsymbol{\hat{b}}$}\right) {\boldsymbol{\cdot}}{\nabla} \end{eqnarray*} (78) is the comoving time derivative along the orbit of the guiding centre of a particle, which includes the motion along a field line. The fully relativistic version can be cast in an analogous form by defining the relativistic generalization of (78): \begin{eqnarray*} \frac{{\rm D}}{{\rm D}\tau } \equiv \frac{{\rm d}}{{\rm d} \tau } + \bar{v}_{\parallel } \, (\mbox{$\boldsymbol{\hat{b}}$}{\boldsymbol{\cdot}}{\boldsymbol {\tilde{{\nabla}}}}) . \end{eqnarray*} (79) Equation (68) then becomes \begin{eqnarray*} \frac{{\rm D}\bar{f}}{{\rm D} \tau } + \left\lbrace \frac{\bar{p}_{\perp }}{2 {{\rm B}}} \, \frac{{\rm D} {{\rm B}}}{{\rm D} \tau } - \nu_{\rm p} \, \frac{\bar{p}_{\perp } \, \bar{v}_{\parallel }}{2} \right\rbrace \, \frac{\partial \bar{f}}{\partial \bar{p}_{\perp }} + \left\lbrace - \frac{\bar{p}_{\perp } \, \bar{v}_{\perp }}{2 {{\rm B}}} \, (\mbox{$\boldsymbol{\hat{b}}$}{\boldsymbol{\cdot}}{\boldsymbol {\tilde{{\nabla}}}}) {{\rm B}}- (\Gamma +1) \, \bar{\gamma } m \, \mbox{$\boldsymbol{\hat{b}}$}{\boldsymbol{\cdot}}\frac{{\rm D}}{{\rm D} \tau } \left(\frac{{\boldsymbol{U}}}{\Gamma +1} \right) + \nu_{\rm p} \, \frac{\bar{p}_{\perp } \, \bar{v}_{\perp }}{2} + q \bar{E}_{\parallel } \right\rbrace \, \frac{\partial \bar{f}}{\partial \bar{p}_{\parallel }} = - \frac{\bar{f} - f_{0}}{\bar{\tau }_{\rm r}} . \nonumber\\ \end{eqnarray*} (80) In the special case that the FRF is an inertial frame so that νp = 0 and |${\rm D}({\boldsymbol{U}}/\Gamma + 1)/{\rm D} \tau = 0$| and in the absence of scattering equation (80) simplifies to \begin{eqnarray*} \frac{{\rm D}\bar{f}}{{\rm D} \tau } + \left(\frac{\bar{p}_{\perp }}{2 {{\rm B}}} \, \frac{{\rm D} {{\rm B}}}{{\rm D} \tau } \right) \, \frac{\partial \bar{f}}{\partial \bar{p}_{\perp }} + \left(q \bar{E}_{\parallel } - \frac{\bar{p}_{\perp } \, \bar{v}_{\perp }}{2 {{\rm B}}} \, (\mbox{$\boldsymbol{\hat{b}}$}{\boldsymbol{\cdot}}{\boldsymbol {\tilde{{\nabla}}}}) {{\rm B}}\right) \, \frac{\partial \bar{f}}{\partial \bar{p}_{\parallel }} = 0 . \end{eqnarray*} (81) This is equivalent with \begin{eqnarray*} \frac{{\rm D} {\bar{p}}_{\perp }}{{\rm D} \tau } = \frac{\bar{p}_{\perp }}{2 {{\rm B}}} \, \frac{{\rm D} {{\rm B}}}{{\rm D} \tau } \, \Longleftrightarrow \, \frac{\bar{p}_{\perp }^2}{2 {{\rm B}}} = \mbox{constant}, \quad \frac{{\rm D} {\bar{p}}_{\parallel }}{{\rm D} \tau } = q \bar{E}_{\parallel } - \frac{\bar{p}_{\perp } \, \bar{v}_{\perp }}{2 {{\rm B}}} \, (\mbox{$\boldsymbol{\hat{b}}$}{\boldsymbol{\cdot}}{\boldsymbol {\tilde{{\nabla}}}}) {{\rm B}}. \end{eqnarray*} (82) The first relation respectively corresponds with the adiabatic invariance of |$\bar{p}_{\perp }^2/2 {{\rm B}}$| if gyro-motion is sufficiently rapid as assumed here. The second term on the right-hand side of the second relation is the magnetic mirror force that occurs under the same circumstances if the strength of the large-scale field varies (see for instance Melrose 1989, ch. 12). 7 LINK WITH WEBB’s APPROACH The above derivation of the transport equation (40) is more intuitive than that of Webb (1985). In addition there are a number of notational differences with Webb’s paper that are designed to make a direct comparison with the previous results of Skilling (1975) and Kulsrud (1983) more simple. Other than that our results are equivalent with Webb’s, as we show briefly here. In both cases the calculation is based on the evaluation of the fictitious force (9). That expression (for flat space–time) can be written as \begin{eqnarray*} \frac{{\rm d} p_{\hat{c}}}{{\rm d} \lambda } = \omega^{\hat{a}}_{\hat{b}\hat{c}}\, \frac{p_{\hat{a}} p^{\hat{b}}}{m}, \omega^{\hat{a}}_{\hat{b}\hat{c}}\equiv - {\rm E}^{\mu }_{\hat{c}} \, {\rm E}^{\nu }_{\hat{b}} \, \frac{\partial {\rm E}_{\mu }^{\hat{a}}}{\partial x^{\nu }} \,, \end{eqnarray*} (83) which is equivalent to equation (6.4) in Webb (1985), who uses the notation |$\Gamma^{\hat{a}}_{\hat{c}\hat{b}}$| for |$\omega^{\hat{a}}_{\hat{b}\hat{c}}$|⁠, essentially taking the approach mentioned in Section 3 that treats the effects of acceleration as space–time curvature in the accelerated frame. The |$\omega^{\hat{a}}_{\hat{b}\hat{c}}$| indeed play the role of the connection coefficients, as follows rigorously if one derives these results using the so-called tetrad formalism (see Appendix A). Using |$p_{\hat{i}} = - p^{\hat{i}}$| for |$\hat{i} = 1\ ,, 2\ ,, 3$| with the |$p^{\hat{i}}$| the components of |$\bar{{\boldsymbol{p}}}$|⁠, a straightforward calculation using equation (32) or equations (14) and (15) yields \begin{eqnarray*} \omega^{\hat{0}}_{\hat{0}\hat{0}} & = & 0 \,, \nonumber \\ \omega^{\hat{i}}_{\hat{0}\hat{0}} = \omega^{\hat{0}}_{\hat{0}\hat{i}} & = & \left(\Gamma + 1 \right) \, \frac{{\rm d}}{{\rm d} \tau } \left(\frac{U^{i}}{\Gamma + 1} \right) \,, \nonumber \\ \omega^{\hat{i}}_{\hat{j}\hat{0}} = \omega^{\hat{0}}_{\hat{j}\hat{i}} & = & \left(\Gamma + 1 \right) \, \tilde{\nabla }_{j} \left(\frac{U^{i}}{\Gamma + 1} \right) \,, \nonumber \\ \omega^{\hat{j}}_{\hat{0}\hat{i}} & = & \frac{1}{\Gamma + 1} \, \left(U^{i} \, \frac{{\rm d} U^{j}}{{\rm d} \tau } - U^{j} \, \frac{{\rm d} U^{i}}{{\rm d} \tau } \right) = (\Gamma - 1) \left(\hat{\ell }^{i} \, \frac{{\rm d} \hat{\ell }^{j}}{{\rm d} \tau } - \hat{\ell }^{j} \, \frac{{\rm d} \hat{\ell }^{i}}{{\rm d} \tau } \right)\nonumber \\ \omega^{\hat{j}}_{\hat{k}\hat{i}} & = & \frac{1}{\Gamma + 1} \, \left(U^{i} \, \tilde{\nabla }_{k} U^{j} - U^{j} \, \tilde{\nabla }_{k} U^{i} \right) = (\Gamma - 1) \left(\hat{\ell }^{i} \, \tilde{\nabla }_{k} \hat{\ell }^{j} - \hat{\ell }^{j} \, \tilde{\nabla }_{k} \hat{\ell }^{i} \right) . \end{eqnarray*} (84) Here we write the components of the unit three-vector |$\boldsymbol{\hat{\ell}} = {\boldsymbol{U}}/U$| as |$\hat{\ell }^{i}$| with indices |$(i\,, j)$| and |$(\hat{i}, \hat{j})$| run from 1 to 3. Using the definitions of |${\rm d}/{\rm d} \tau = \Gamma \, ({\rm d}/{\rm d}t)$|⁠, |${\boldsymbol {\tilde{{\nabla}}}}$| (see equations 33 and 34) and |$U^{i} = \Gamma \, V^{i}$| one can show that these expressions are equivalent to Webb’s results, equation (6.9) in Webb (1985) which incidentally does not list all coefficients. We provide two examples: Webb lists |$\omega^{\hat{j}}_{\hat{0}\hat{i}}$| (our notation) by the equivalent expression that follows from |$U^{i} = \Gamma \, V^{i}$| and |${\rm d}/{\rm d}\tau = \Gamma \, ({\rm d}/{\rm d}t)$|⁠: \begin{eqnarray*} \omega^{\hat{j}}_{\hat{0}\hat{i}} = \frac{\Gamma^3}{(\Gamma +1) c^3} \left(V^{i} \, \frac{{\rm d} V^{j}}{{\rm d} t} - V^{j} \, \frac{{\rm d} V^{i}}{{\rm d} t} \right) \,, \end{eqnarray*} (85) where we reinstated c for the comparison. Webb also uses the relation (with the Einstein summation convention for double indices) \begin{eqnarray*} \omega^{\hat{i}}_{\hat{i}\hat{0}} = (\Gamma + 1) \left[ \, {\boldsymbol {\tilde{{\nabla}}}}{\boldsymbol{\cdot}}\left(\frac{{\boldsymbol{U}}}{\Gamma +1} \right) \, \right] = \frac{1}{c} \left(\frac{\partial \Gamma }{\partial t} + {\nabla}{\boldsymbol{\cdot}}(\Gamma {\boldsymbol{V}}) \right) . \end{eqnarray*} (86) 8 THE GENERAL CASE: INCLUDING THE GYRO-PHASE DEPENDENCE In order to fully describe such effects as cosmic ray (or photon) diffusion and cosmic ray viscosity we must calculate the full deviations of |$f(\tilde{{\boldsymbol Z }})$| from isotropy in |$\bar{{\boldsymbol{p}}}$|-space, now including the dependence on gyro-phase that was neglected in earlier calculations and in Sections 5 and 6. For photons the concept of gyro-averaging makes no sense to begin with. To that end one must solve the Vlasov equation \begin{eqnarray*} \dot{\tilde{{\boldsymbol Z }}} {\boldsymbol{\cdot}}\frac{\partial f(\tilde{{\boldsymbol Z }})}{\partial \tilde{{\boldsymbol Z }}} = \dot{t} \, \frac{\partial f}{\partial t} + \dot{{\boldsymbol x }} {\boldsymbol{\cdot}}{\nabla}f + {\boldsymbol {\bar{K}} } {\boldsymbol{\cdot}}\frac{\partial f}{\partial \bar{{\boldsymbol{p}}}} - \bar{\Omega } \frac{\partial f}{\partial \bar{\phi }} = - \frac{f - f_{0}}{\bar{\tau }_{\rm r}} \,, \end{eqnarray*} (87) making an expansion for small anisotropies in the fluid rest frame. Here |${\boldsymbol {\bar{K}} } = {\boldsymbol {\bar{K}} }_{\rm f} + \bar{{\boldsymbol{p}}}{\boldsymbol{\times}}{\boldsymbol {\bar{\omega}}}$| (see equation 44). This is a natural assumption in a relativistic flow of a scattering medium. Generally speaking, this assumption does not mean that the anisotropy in the distribution is small in any other frame. If the fluid moves at relativistic speeds, this is true in particular in the laboratory frame, due to the effect of relativistic beaming. 8.1 The diffusion approximation and the choice for |$f_{0}({\boldsymbol x }\,, t\,, \bar{p})$| We assume that the gradient scale L and the time-scale T of the variations in |$f(\tilde{{\boldsymbol Z }})$| and |${\boldsymbol{U}} ({\boldsymbol x }\,, t)$| satisfy \begin{eqnarray*} \frac{V \bar{\tau }_{\rm r}}{L}\,, \quad \frac{\bar{v}\bar{\tau }_{\rm r}}{L}\,, \quad \frac{V}{\bar{\Omega }L}\,, \quad \frac{\bar{v}}{\bar{\Omega } L}\,, \quad \frac{\bar{\tau }_{\rm r}}{T}\,, \quad \frac{1}{\bar{\Omega } T} = {\cal O}(\epsilon) \ll 1 . \end{eqnarray*} (88) This can be formalized by introducing ε as an ordering parameter (to be put equal to unity at the end of the calculation) and by making the replacements \begin{eqnarray*} \bar{\tau }_{\rm r} \, \Longrightarrow \, \varepsilon \, \bar{\tau }_{\rm r}\,, \quad \bar{\Omega } \, \Longrightarrow \, \bar{\Omega }/\varepsilon . \end{eqnarray*} (89) Note that this leaves the dimensionless parameter |$\bar{\Omega } \, \bar{\tau }_{\rm r}$| unchanged so that this parameter can take any value. Then equation (87) takes the form \begin{eqnarray*} \varepsilon \, \left(\dot{t} \, \frac{\partial f}{\partial t} + \dot{{\boldsymbol x }} {\boldsymbol{\cdot}}{\nabla}f + {\boldsymbol {\bar{K}} } {\boldsymbol{\cdot}}\frac{\partial f}{\partial \bar{{\boldsymbol{p}}}} \right) = \bar{\Omega } \frac{\partial f}{\partial \bar{\phi }} + \overline{{\rm St}}(f) = \bar{\Omega } \frac{\partial f}{\partial \bar{\phi }} - \frac{f - f_{0}}{\bar{\tau }_{\rm r}} . \end{eqnarray*} (90) One expands the Vlasov distribution |$f(\tilde{{\boldsymbol Z }})$| as \begin{eqnarray*} f(\tilde{{\boldsymbol Z }}) = f_{0}(t\,, \quad {\boldsymbol x }\,, \quad |\bar{{\boldsymbol{p}}}|) + \varepsilon \, f_{1}(\tilde{{\boldsymbol Z }}) + \varepsilon^2 \, f_{2}(\tilde{{\boldsymbol Z }}) + \ldots . \end{eqnarray*} (91) The leading term f0 is isotropic in |$\bar{{\boldsymbol{p}}}$|-space while f1, f2, etc. contain the effect of anisotropies. In the diffusion approximation (as is employed here) one breaks off expansion (equation 91) after the first-order term and never calculates |$f_{2}(\tilde{{\boldsymbol Z }})$|⁠. Equation (90) implies that |$f_{0}(t\,, {\boldsymbol x }\,, |\bar{{\boldsymbol{p}}}|)$| should be chosen as the average of the Vlasov distribution over solid angle in |$\bar{p}$|-space to all orders in ε. It is the most general ‘equilibrium distribution’ that (i) annihilates the combined gyration/scattering operator on the right-hand side equation (90) and (ii) is consistent with the exact Vlasov equation (90). This implies \begin{eqnarray*} f_{0}({\boldsymbol x }\,, t\,, \bar{{\boldsymbol{p}}}) = \frac{1}{4 {\pi} } \int_{-1}^{+1} {\rm d} \bar{\mu } \, \int_{0}^{2 {\pi} } {\rm d} \bar{\phi } \, f({\boldsymbol x }\,, t\,, \bar{{\boldsymbol{p}}}) \equiv \left\langle f(\tilde{{\boldsymbol Z }}) \right\rangle \,, \end{eqnarray*} (92) where we introduce the notation ⟨···⟩ for the average over momentum solid angle. In the case of pure scattering, which is the case we consider here, relation (92) also follows from the fact that scattering merely redistributes particles over the solid angle |$\bar{\Omega }_{p}$| in momentum space, so that the scattering operator |$\overline{{\rm St}}(f)$| (or |${\rm St}(f) = \bar{\gamma } \, \overline{{\rm St}}(f)$|⁠) must satisfy \begin{eqnarray*} \int_{4 {\pi} } {\rm d} \bar{\Omega }_{p} \, \overline{{\rm St}}(f) = 0 \,, \end{eqnarray*} (93)regardless of the precise form of the relaxation term |$\overline{{\rm St}}(f)$|⁠. We note in passing that, because of condition (93), the use of the algebraic Anderson–Witting term for the scattering operator does not suffer from the problem that it violates particle conservation. This problem does occur in relativistic theories such as this whenever one wants to use such a term to force the distribution into a certain distribution in |$\bar{{\cal E}}= \sqrt{\bar{p}^2 + m^2}$| or |$\bar{p}$|⁠, such as a relativistic Maxwellian. In the terminology of perturbation theory the evolution equation for f0 (derived below) should contain all ‘secular’ effects. The higher-order terms f1, f2, etc. are ‘oscillating terms’ that contain all anisotropies in momentum space. Because of equation (92) these terms satisfy |$\left\langle f_{\rm n}(\tilde{{\boldsymbol Z }}) \right\rangle = 0$| for n ≥ 1: Their average over solid angle in |$\bar{{\boldsymbol{p}}}$| space vanishes. In particular one has |$\langle f_{1}(\tilde{{\boldsymbol Z }})\rangle = 0$|⁠, which will prove important in our discussion of particle diffusion (this paper) and viscosity (Paper 2) and even more so in the discussion of viscous stresses due to radiation in Paper 3. As a result |$f_{0} = \langle f (\tilde{{\boldsymbol Z }}) \rangle $| will in principle include |${\cal O}(\varepsilon)$| and higher-order corrections that can only be calculated by solving the self-consistent transport equation for f0. The subscript ‘0’ should therefore not be interpreted to mean that f0 is of order |${\cal O}(\varepsilon^{0})$|⁠: It can be expanded as |$f_{0} = f_{00} + \varepsilon \, f_{01} + \ldots$|⁠. In practice this expansion is never needed. One version of this transport equation for the average distribution (to first order in ε) is obtained by averaging equation (90) over solid angle: \begin{eqnarray*} \Gamma \, \left[ \, \frac{\partial f_{0}}{\partial t} + ({\boldsymbol{V}} {\boldsymbol{\cdot}}{\nabla}) f_{0} \, \right] + \left\langle {\boldsymbol {\bar{K}} }_{\rm f} {\boldsymbol{\cdot}}\mbox{${\boldsymbol {\hat{n}}}$}\right\rangle \, \frac{\partial f_{0}}{\partial \bar{p}}= - \varepsilon \left\lbrace \left\langle \bar{{\boldsymbol{V}}}{\boldsymbol{\cdot}}{\boldsymbol {\tilde{{\nabla}}}}f_{1} \right\rangle + \left\langle {\boldsymbol {\tilde{K}}} {\boldsymbol{\cdot}}\frac{\partial f_{1}}{\partial \bar{{\boldsymbol{p}}}} \right\rangle \right\rbrace + {\cal O}(\varepsilon^2) . \end{eqnarray*} (94) Here we have transported all explicit order ε terms to the right-hand side, and defined |${\boldsymbol {\tilde{K}}} \equiv {\boldsymbol {\bar{K}} } - \left\langle {\boldsymbol {\bar{K}} } \right\rangle$|⁠, the oscillating part of the fictitious force. Also |$\mbox{${\boldsymbol {\hat{n}}}$}= \bar{{\boldsymbol{p}}}/\bar{p}$| is the unit vector along the FRF momentum and we have used |$\partial f_{0}/\partial \bar{{\boldsymbol{p}}}= \mbox{${\boldsymbol {\hat{n}}}$}\, (\partial f_{0}/\partial \bar{p})$|⁠. The second term on the left-hand side of equation (94) is the average change of |$\bar{p}$| due to scattering by the flow. This term describes adiabatic losses/gains and was also found by Webb (1985). Using equation (49) it can be written as \begin{eqnarray*} \left\langle {\boldsymbol {\bar{K}} }_{\rm f} {\boldsymbol{\cdot}}\mbox{${\boldsymbol {\hat{n}}}$}\right\rangle = - \frac{(\Gamma + 1) \, \bar{p}}{3} \, \left[ \, {\boldsymbol {\tilde{{\nabla}}}}{\boldsymbol{\cdot}}\left(\frac{{\boldsymbol{U}}}{\Gamma +1} \right) \, \right] = -\frac{\bar{p}}{3} \, \left(\frac{\partial \Gamma }{\partial t} + {\nabla}{\boldsymbol{\cdot}}(\Gamma {\boldsymbol{V}}) \right) = - \frac{\bar{p}}{3} \, (\nabla {\boldsymbol{\cdot}}U) . \end{eqnarray*} (95) The two terms |${\cal O}(\varepsilon)$| terms on the right-hand side of equation (94) respectively lead to spatial diffusion plus a drift velocity (see Section 9) and to irreversible (viscous) heating of the scattered particles, as shown in Paper 2. 8.2 Solution for |$f_{1}(\tilde{{\boldsymbol Z }})$| To first order in ε Vlasov equation (90) yields the following equation for f1: \begin{eqnarray*} \bar{\Omega } \frac{\partial f_{1}}{\partial \bar{\phi }} - \frac{f_{1}}{\bar{\tau }_{\rm r}} = \frac{{\rm d} f_{0}}{{\rm d} \tau } + \bar{{\boldsymbol{V}}}{\boldsymbol{\cdot}}{\boldsymbol {\tilde{{\nabla}}}}f_{0} + \left({\boldsymbol {\bar{K}} }_{\rm f} {\boldsymbol{\cdot}}\mbox{${\boldsymbol {\hat{n}}}$}\right) \, \frac{\partial f_{0}}{\partial \bar{p}} . \end{eqnarray*} (96) Using relation (94) to eliminate df0/dτ, neglecting the |${\cal O}(\varepsilon)$| term so that \begin{eqnarray*} \frac{{\rm d} f_{0}}{{\rm d} \tau } = \frac{(\Gamma + 1) \, \bar{p}}{3} \, \left[ \, {\boldsymbol {\tilde{{\nabla}}}}{\boldsymbol{\cdot}}\left(\frac{{\boldsymbol{U}}}{\Gamma +1} \right) \, \right]\, \bar{p}\, \frac{\partial f_{0}}{\partial \bar{p}} + {\cal O}(\varepsilon) \,, \end{eqnarray*} (97) one finds \begin{eqnarray*} \bar{\Omega } \frac{\partial f_{1}}{\partial \bar{\phi }} - \frac{f_{1}}{\bar{\tau }_{\rm r}} = \bar{{\boldsymbol{V}}}{\boldsymbol{\cdot}}{\boldsymbol {\tilde{{\nabla}}}}f_{0} + \tilde{K}(\tilde{{\boldsymbol Z }}) \, \frac{\partial f_{0}}{\partial \bar{p}} + {\cal O}(\varepsilon) . \end{eqnarray*} (98) Here |$\tilde{K}(\tilde{{\boldsymbol Z }})$| is the ‘oscillating’ part of the rate of change of |$\bar{p}$|⁠, formally defined as |$\tilde{K}(\tilde{{\boldsymbol Z }}) = \dot{\bar{p}} - \left\langle \dot{\bar{p}} \right\rangle$|⁠: \begin{eqnarray*} \tilde{K}(\tilde{{\boldsymbol Z }}) \equiv {\boldsymbol {\bar{K}} }_{\rm f} {\boldsymbol{\cdot}}\mbox{${\boldsymbol {\hat{n}}}$}- \left\langle {\boldsymbol {\bar{K}} }_{\rm f} {\boldsymbol{\cdot}}\mbox{${\boldsymbol {\hat{n}}}$}\right\rangle = \frac{\bar{p}}{3} \, \left[ \left(\Gamma + 1 \right) {\boldsymbol {\tilde{{\nabla}}}}{\boldsymbol{\cdot}}\left(\frac{{\boldsymbol{U}}}{\Gamma +1} \right)\right] - (\Gamma +1) \, m \, \left(\mbox{${\boldsymbol {\hat{n}}}$}{\boldsymbol{\cdot}}\frac{{\rm d}}{{\rm d}\lambda } \left(\frac{{\boldsymbol{U}}}{\Gamma +1} \right)\right) . \end{eqnarray*} (99) When written out in terms of the local Cartesian xyz coordinates introduced in Section 5, equation (98) reads \begin{eqnarray*} \bar{\Omega } \, \frac{\partial f_{1}}{\partial \bar{\phi }} - \frac{f_{1}}{\bar{\tau }_{\rm r}} & = & \bar{v}\, \left\lbrace \sqrt{1 - \bar{\mu }^2}\, \left[ \cos \bar{\phi } \, {{\boldsymbol {\hat{x}}}}+ \sin \bar{\phi } \, {{\boldsymbol {\hat{y}}}}\right] + \bar{\mu } \, \mbox{$\boldsymbol{\hat{b}}$}\right\rbrace {\boldsymbol{\cdot}}{\boldsymbol {\tilde{{\nabla}}}}f_{0} \nonumber \\ && -\, \bar{\gamma } m \, (\Gamma + 1) \, \left\lbrace \left(\sqrt{1 - \bar{\mu }^2}\, \left[ \cos \bar{\phi } \, {{\boldsymbol {\hat{x}}}}+ \sin \bar{\phi } \, {{\boldsymbol {\hat{y}}}}\right] + \bar{\mu } \, \mbox{$\boldsymbol{\hat{b}}$}\right) {\boldsymbol{\cdot}}\frac{{\rm d}}{{\rm d} \tau } \left(\frac{{\boldsymbol{U}}}{\Gamma +1} \right) \right\rbrace \, \frac{\partial f_{0}}{\partial \bar{p}} \nonumber \\ && -\, (\Gamma + 1) \, \left\lbrace \left[ \, \left(1 - \bar{\mu }^2 \right) \, \mbox{$\boldsymbol {\rm M} $}_{1}+ \bar{\mu } \, \sqrt{1 - \bar{\mu }^2}\, \mbox{$\boldsymbol {\rm M} $}_{2}\right] {{\boldsymbol \, : \,}}{\boldsymbol {\tilde{{\nabla}}}}\left(\frac{{\boldsymbol{U}}}{\Gamma + 1} \right) \right\rbrace \, \bar{p}\, \frac{\partial f_{0}}{\partial \bar{p}} \nonumber \\ && -\, (\Gamma + 1) \, \left\lbrace \frac{3 \bar{\mu }^2 - 1}{2} \, \left[ \, \mbox{$\boldsymbol{\hat{b}}$}\, \mbox{$\boldsymbol{\hat{b}}$}{{\boldsymbol \, : \,}}{\boldsymbol {\tilde{{\nabla}}}}\left(\frac{{\boldsymbol{U}}}{\Gamma +1} \right) - \frac{1}{3}\, {\boldsymbol {\tilde{{\nabla}}}}{\boldsymbol{\cdot}}\left(\frac{{\boldsymbol{U}}}{\Gamma +1} \right) \, \right] \right\rbrace \, \bar{p}\, \frac{\partial f_{0}}{\partial \bar{p}} . \end{eqnarray*} (100) Here we have put |${\boldsymbol {\bar{E}} }_{\parallel } = 0$| for simplicity, used |$\bar{{\boldsymbol{V}}}= \bar{v}\, \left(\sqrt{1 - \bar{\mu }^2}\, \left[ \cos \bar{\phi } \, {{\boldsymbol {\hat{x}}}}+ \sin \bar{\phi } \, {{\boldsymbol {\hat{y}}}}\right] + \bar{\mu } \, \mbox{$\boldsymbol{\hat{b}}$}\right)$| with |$\mbox{$\boldsymbol{\hat{b}}$}= {{\boldsymbol {\hat{z}}}}$|⁠, together with the fact that f0 is isotropic. The Coriolis term in the fictitious force does not contribute at this order as |$| {\boldsymbol {\bar{\omega}}}|/\bar{\Omega }$| and |$| {\boldsymbol {\bar{\omega}}}| \, \bar{\tau }_{\rm r}$| are small (first-order in ε) quantities under the assumptions of this derivation, and since f0 is isotropic. The two 3 × 3 matrices |$\mbox{$\boldsymbol {\rm M} $}_{1}$| and |$\mbox{$\boldsymbol {\rm M} $}_{2}$| are formally defined as \begin{eqnarray*} \mbox{$\boldsymbol {\rm M} $}_{1}= \frac{1}{2}\, \left(\begin{array}{@{}c@{\quad}c@{\quad}c@{}}2 \cos^2\bar{\phi } - 1 & 2 \cos {\bar{\phi }}\, \sin {\bar{\phi }}& 0 \\ 2 \cos {\bar{\phi }}\sin {\bar{\phi }}& 2 \sin^2 \bar{\phi } - 1 & 0 \\ 0 & 0 & 0 \\ \end{array} \right) = \frac{1}{2}\, \left(\begin{array}{@{}c@{\quad}c@{\quad}c@{}}\cos 2 \bar{\phi } & \sin 2 \bar{\phi } & 0 \\ \sin 2 \bar{\phi } & - \cos 2 \bar{\phi } & 0 \\ 0 & 0 & 0 \\ \end{array} \right), \mbox{$\boldsymbol {\rm M} $}_{2}= \left(\begin{array}{@{}c@{\quad}c@{\quad}c@{}}0 & 0 & \cos {\bar{\phi }}\\ 0 & 0 & \sin {\bar{\phi }}\\ \cos {\bar{\phi }}& \sin {\bar{\phi }}& 0 \\ \end{array} \right) . \end{eqnarray*} (101) The double contractions |$\mbox{$\boldsymbol {\rm M} $}_{1,2} {{\boldsymbol \, : \,}}{\boldsymbol {\tilde{{\nabla}}}}({\boldsymbol{U}}/\Gamma +1)$| stand for |$M_{1,2}^{ij} \tilde{\nabla }_{j} (U^{i}/\Gamma +1)$|⁠. The gyro-average of |$\mbox{$\boldsymbol {\rm M} $}_{1,2}$| vanishes. This means that all terms on the right-hand side of equation (100) vanish when one averages over momentum solid angle in the FRF, which implies ⟨f1⟩ =0, as required. Integration of equation (100) for |$f_{1}({\boldsymbol x }\,, t\,, \bar{{\boldsymbol{p}}})$| effectively corresponds to the replacement \begin{eqnarray*} \cos (m \bar{\phi }) \, \Longrightarrow \, {\displaystyle - \frac{\bar{\tau }_{\rm r}}{1 + m^2 \bar{\Omega }^2 \bar{\tau }_{\rm r}^2} \left[ \, \cos (m \bar{\phi }) - m \bar{\Omega } \bar{\tau }_{\rm r} \, \sin (m \bar{\phi }) \, \right] }\,, \quad \sin (m \bar{\phi }) \, \Longrightarrow \, {\displaystyle - \frac{\bar{\tau }_{\rm r}}{1 + m^2 \bar{\Omega }^2 \bar{\tau }_{\rm r}^2} \left[ \, \sin (m \bar{\phi }) + m \bar{\Omega } \bar{\tau }_{\rm r} \, \cos (m \bar{\phi }) \, \right] } \,, \end{eqnarray*} (102) for |$m = 0, 1, 2$|⁠. Furthermore, the matrices |$\mbox{$\boldsymbol {\rm M} $}_{1}$| and |$\mbox{$\boldsymbol {\rm M} $}_{2}$| are both traceless and symmetric. This implies that in the double contraction with |${\boldsymbol {\tilde{{\nabla}}}}[{\boldsymbol{U}}/(\Gamma +1)]$| we can make the replacement \begin{eqnarray*} {\boldsymbol {\tilde{{\nabla}}}}\left(\frac{{\boldsymbol{U}}}{\Gamma +1} \right) \, \Longrightarrow \frac{1}{2}\, \left\lbrace {\boldsymbol {\tilde{{\nabla}}}}\left(\frac{{\boldsymbol{U}}}{\Gamma +1} \right) + \left[ {\boldsymbol {\tilde{{\nabla}}}}\left(\frac{{\boldsymbol{U}}}{\Gamma +1} \right) \right]^{\rm tr} - \frac{2}{3} \, \left[ {\boldsymbol {\tilde{{\nabla}}}}{\boldsymbol{\cdot}}\left(\frac{{\boldsymbol{U}}}{\Gamma +1} \right) \right]\, {\mathbf{ I}} \right\rbrace . \end{eqnarray*} (103) Here the superscript ‘tr’ stands for ‘transpose’ and |${ \mathbf{ I}} \equiv {\rm diag}(1\,, 1\,, 1)$| is the unit 3 × 3 tensor in three-space. To facilitate notation we will define the symmetrized and traceless rate-of-strain tensor |$\mbox{$\boldsymbol {\rm \mathit{ W}} $}$|⁠: \begin{eqnarray*} \mbox{$\boldsymbol {\rm \mathit{ W}} $} \equiv (\Gamma + 1) \, \left\lbrace {\boldsymbol {\tilde{{\nabla}}}}\left(\frac{{\boldsymbol{U}}}{\Gamma +1} \right) + \left[ {\boldsymbol {\tilde{{\nabla}}}}\left(\frac{{\boldsymbol{U}}}{\Gamma +1} \right) \right]^{\rm tr} - \frac{2}{3} \, \left[ {\boldsymbol {\tilde{{\nabla}}}}{\boldsymbol{\cdot}}\left(\frac{{\boldsymbol{U}}}{\Gamma +1} \right) \right]\, {\mathbf{I}} \right\rbrace . \end{eqnarray*} (104) For a non-relativistic flow (Γ ≃ 1, |${\boldsymbol{U}} \simeq {\boldsymbol{V}}$| with |$|{\boldsymbol{V}}| \ll 1$|⁠, |${\boldsymbol {\tilde{{\nabla}}}}\simeq {\nabla}$|⁠) this reduces to the rate of strain tensor as defined by Landau & Lifshitz (1959) and Braginskii (1965): \begin{eqnarray*} \mbox{$\boldsymbol {\rm \mathit{ W}} $}_{\rm nr} = {\nabla}{\boldsymbol{V}} + \left({\nabla}{\boldsymbol{V}} \right)^{\rm tr} - \frac{2}{3} \, ({\nabla}{\boldsymbol{\cdot}}{\boldsymbol{V}}) \, {\mathbf{I}}. \end{eqnarray*} (105) With definition (104) one has: \begin{eqnarray*} \mbox{$\boldsymbol {\rm M} $}_{1,2} {{\boldsymbol \, : \,}}{\boldsymbol {\tilde{{\nabla}}}}\left(\frac{{\boldsymbol{U}}}{\Gamma +1} \right) = \frac{1}{2}\,\mbox{$\boldsymbol {\rm M} $}_{1,2} {{\boldsymbol \, : \,}}\mbox{$\boldsymbol {\rm { W}} $}, \mbox{$\boldsymbol{\hat{b}}$}\, \mbox{$\boldsymbol{\hat{b}}$}{{\boldsymbol \, : \,}}{\boldsymbol {\tilde{{\nabla}}}}\left(\frac{{\boldsymbol{U}}}{\Gamma +1} \right) - \frac{1}{3}\, {\boldsymbol {\tilde{{\nabla}}}}{\boldsymbol{\cdot}}\left(\frac{{\boldsymbol{U}}}{\Gamma +1} \right) = \frac{1}{2}\, \mbox{$\boldsymbol{\hat{b}}$}\mbox{$\boldsymbol{\hat{b}}$}{{\boldsymbol \, : \,}}\mbox{$\boldsymbol {\rm { W}} $} . \end{eqnarray*} (106) It is also useful to define \begin{eqnarray*} {\boldsymbol {\bar{a}} } \equiv (\Gamma +1) \, \frac{{\rm d}}{{\rm d} \tau } \left(\frac{{\boldsymbol{U}}}{\Gamma +1} \right) \,, \end{eqnarray*} (107) which corresponds to the non-vanishing spatial part of the four-acceleration of the fluid in the FRF: |$a^{\hat{a}} = (0\ ,, {\boldsymbol {\bar{a}} })$| (see Appendix C). One finds that \begin{eqnarray*} f_{1}({\boldsymbol x }\,, t\,, \bar{{\boldsymbol{p}}}) = - \bar{\tau }_{\rm r} \, \left(\bar{v}\, {\boldsymbol {\bar{\Sigma}}} {\boldsymbol{\cdot}}{\boldsymbol {\tilde{{\nabla}}}}\right) f_{0} + \bar{\tau }_{\rm r} \, \left({\boldsymbol {\bar{\Sigma}}} {\boldsymbol{\cdot}}{\boldsymbol {\bar{a}} } \right) \, \bar{\gamma } m \, \frac{\partial f_{0}}{\partial \bar{p}} + \bar{\tau }_{\rm r} \, \left(\mbox{$\boldsymbol {\rm \hat{M}} $} {{\boldsymbol \, : \,}}\mbox{$\boldsymbol {\rm W} $} \right) \, \bar{p}\, \frac{\partial f_{0}}{\partial \bar{p}} . \end{eqnarray*} (108) The FRF vector |${\boldsymbol {\bar{\Sigma}}}$| and FRF three-tensor |$\mbox{$\boldsymbol {\rm \hat{M}} $}$| contain all the magnetic effects. When expressed in terms of components in the local Cartesian xyz frame in the FRF the three-vector |${\boldsymbol {\bar{\Sigma}}}$| equals \begin{eqnarray*} {\boldsymbol {\bar{\Sigma}}} = \sqrt{1 - \bar{\mu }^2}\, \left\lbrace \frac{\displaystyle \cos {\bar{\phi }}\, {{\boldsymbol {\hat{x}}}}+ \sin {\bar{\phi }}\, {{\boldsymbol {\hat{y}}}}}{\displaystyle {1 + \bar{\Omega }^2 \bar{\tau }_{\rm r}^2}} + \bar{\Omega } \bar{\tau }_{\rm r} \, \left(\frac{\displaystyle \cos {\bar{\phi }}\, {{\boldsymbol {\hat{y}}}}- \sin {\bar{\phi }}\, {{\boldsymbol {\hat{x}}}}}{\displaystyle {1 + \bar{\Omega }^2 \bar{\tau }_{\rm r}^2}} \right) \right\rbrace + \bar{\mu } \, \mbox{$\boldsymbol{\hat{b}}$}. \end{eqnarray*} (109) If we write |$\bar{{\boldsymbol{p}}}= \bar{{\boldsymbol{p}}}_{\perp } + \bar{p}_{\parallel }\, \mbox{$\boldsymbol{\hat{b}}$}$|⁠, it is easily seen that equation (109) is equivalent with \begin{eqnarray*} {\boldsymbol {\bar{\Sigma}}} = \frac{1}{\bar{p}} \left[ \, \frac{\bar{{\boldsymbol{p}}}_{\perp }}{\displaystyle {1 + \bar{\Omega }^2 \bar{\tau }_{\rm r}^2}} + \frac{\bar{\Omega } \bar{\tau }_{\rm r}}{\displaystyle {1 + \bar{\Omega }^2 \bar{\tau }_{\rm r}^2}} \, \left(\mbox{$\boldsymbol{\hat{b}}$}{\boldsymbol{\times}}\bar{{\boldsymbol{p}}}_{\perp } \right) + \bar{p}_{\parallel }\, \mbox{$\boldsymbol{\hat{b}}$}\, \right]. \end{eqnarray*} (110) The symmetric and traceless 3 × 3 matrix |$\mbox{$\boldsymbol {\rm \hat{M}} $}$| equals \begin{eqnarray*} \mbox{$\boldsymbol {\rm \hat{M}} $} \equiv \frac{1 - \bar{\mu }^2}{2} \, \frac{\mbox{$\boldsymbol {\rm M} $}_{1}- 2 \bar{\Omega } \, \bar{\tau }_{\rm r} \, \mbox{$\boldsymbol {\bar{\rm M}} $}_{1}}{\displaystyle {1 + 4\bar{\Omega }^2 \bar{\tau }_{\rm r}^2}} + \frac{\bar{\mu } \, \sqrt{1 - \bar{\mu }^2}}{2} \, \frac{\mbox{$\boldsymbol {\rm M} $}_{2}- \bar{\Omega } \, \bar{\tau }_{\rm r} \, \mbox{$\boldsymbol {\bar{\rm M}} $}_{2}}{\displaystyle {1 + \bar{\Omega }^2 \bar{\tau }_{\rm r}^2}} + \frac{3 \bar{\mu }^2 - 1}{4} \, \mbox{$\boldsymbol{\hat{b}}$}\, \mbox{$\boldsymbol{\hat{b}}$}. \end{eqnarray*} (111) The matrices |$\mbox{$\boldsymbol {\bar{\rm M}} $}_{1}$| and |$\mbox{$\boldsymbol {\bar{\rm M}} $}_{2}$| are defined as \begin{eqnarray*} \mbox{$\boldsymbol {\bar{\rm M}} $}_{1}= - \frac{{\rm d} \mbox{$\boldsymbol {\rm M} $}_{1}}{{\rm d} \bar{\phi }} = \frac{1}{2}\, \left(\begin{array}{@{}c@{\quad}c@{\quad}c@{}}+ \sin 2 \bar{\phi } & - \cos 2 \bar{\phi } & 0 \\ - \cos 2 \bar{\phi } & - \sin 2 \bar{\phi } & 0 \\ 0 & 0 & 0 \\ \end{array} \right)\,, \quad \mbox{$\boldsymbol {\bar{\rm M}} $}_{2}= - \frac{{\rm d} \mbox{$\boldsymbol {\rm M} $}_{2}}{{\rm d} \bar{\phi }} = \left(\begin{array}{@{}c@{\quad}c@{\quad}c@{}}0 & 0 & \sin {\bar{\phi }}\\ 0 & 0 & - \cos {\bar{\phi }}\\ \sin {\bar{\phi }}& - \cos {\bar{\phi }}& 0 \\ \end{array} \right) . \end{eqnarray*} (112) 8.3 The unmagnetized case In the absence of a magnetic field (or for uncharged particles or photons scattered by the flow) the expression for |$f_{1}({\boldsymbol x }\,, t\,, \bar{{\boldsymbol{p}}})$| simplifies considerably. Equation (94) for the isotropic part of the distribution remains valid in this case. If we write the particle three-momentum in the FRF as |$\bar{{\boldsymbol{p}}}= \bar{p}\, \mbox{${\boldsymbol {\hat{n}}}$}$|⁠, |$f_{1}({\boldsymbol x }\,, t\,, \bar{{\boldsymbol{p}}})$| is given by equation (108) with \begin{eqnarray*} {\boldsymbol {\bar{\Sigma}}} = \frac{\bar{{\boldsymbol{p}}}}{\bar{p}} = \mbox{${\boldsymbol {\hat{n}}}$}, \mbox{$\boldsymbol {\rm \hat{M}} $} = \frac{1}{2}\, \left(\mbox{${\boldsymbol {\hat{n}}}$}\mbox{${\boldsymbol {\hat{n}}}$}- \langle \mbox{${\boldsymbol {\hat{n}}}$}\mbox{${\boldsymbol {\hat{n}}}$}\rangle \right) = \frac{1}{2}\, \left(\mbox{${\boldsymbol {\hat{n}}}$}\, \mbox{${\boldsymbol {\hat{n}}}$}- \frac{1}{3}\, {\mathbf{I}} \right) \,. \end{eqnarray*} (113) This yields, using |$\bar{{\boldsymbol{V}}}= \bar{v}\, \mbox{${\boldsymbol {\hat{n}}}$}$| , \begin{eqnarray*} f_{1}({\boldsymbol x }\,, t\,, \bar{{\boldsymbol{p}}}) = - \bar{\tau }_{\rm r} \bar{v}(\mbox{${\boldsymbol {\hat{n}}}$}{\boldsymbol{\cdot}}{\boldsymbol {\tilde{{\nabla}}}}) f_{0} + \bar{\tau }_{\rm r} \, \left(\mbox{${\boldsymbol {\hat{n}}}$}{\boldsymbol{\cdot}}{\boldsymbol {\bar{a}} } \right) \, \bar{\gamma } m \, \frac{\partial f_{0}}{\partial \bar{p}} + \frac{\bar{\tau }_{\rm r}}{2} \, \left(\mbox{${\boldsymbol {\hat{n}}}$}\, \mbox{${\boldsymbol {\hat{n}}}$}{{\boldsymbol \, : \,}}\mbox{$\boldsymbol {\rm W} $} \right) \, \bar{p}\, \frac{\partial f_{0}}{\partial \bar{p}} . \end{eqnarray*} (114) In the last term we have used that |$\mbox{$\boldsymbol {\rm W} $}$| is traceless, which implies |${\mathbf{I}} {{\boldsymbol \, : \,}}\mbox{$\boldsymbol {\rm W} $} = {\rm trace}(\mbox{$\boldsymbol {\rm W} $}) = 0$|⁠. 9 APPLICATION: DIFFUSION AND DRIFT As a first application of the theory presented here we consider spatial diffusion and particle drift in a relativistic flow. The calculation of the viscosity and its effects is rather more involved and needs a more careful discussion of the assumptions when we compare with the existing results from radiation hydrodynamics for a non-relativistic flow. For this reason the calculation of the stress tensor is relegated to Paper 2. Since we are using FRF momentum variables, it is easiest to first calculate the number flux four-vector in the FRF. With c = 1 and |${\boldsymbol {\bar{v}} } \equiv \bar{{\boldsymbol{p}}}/\bar{\gamma }m$| the FRF particle velocity the number flux four-vector |$N^{\hat{a}}$| is given by: \begin{eqnarray*} N^{\hat{a}}(x^{\mu }) = \int \frac{{\rm d}^3 \bar{{\boldsymbol{p}}}}{\bar{{\cal E}}} \, f(\tilde{{\boldsymbol Z }}) \, p^{\hat{a}} = \int {\rm d}^3 \bar{{\boldsymbol{p}}}\, f(\tilde{{\boldsymbol Z }}) \, (1\,, \quad {\boldsymbol {\bar{v}} }) . \end{eqnarray*} (115) The corresponding four-vector in the LF follows from the Lorentz transformation: |$N^{\mu } = {\rm E}^{\mu }_{\hat{a}} \, N^{\hat{a}}$|⁠. The quantity |${\rm d}^3 {\boldsymbol{p}}/{\cal E}$| is a Lorentz invariant (see for instance Landau & Lifshitz 1975, ch. 2). Making the diffusion approximation, |$f \simeq f_{0} + \varepsilon \, f_{1}$|⁠, one finds (in obvious notation) |$N^{\hat{a}} = N^{\hat{a}}_{0} + \varepsilon \, N^{\hat{a}}_{1}$|⁠, with \begin{eqnarray*} N^{\hat{a}}_{0} = (\bar{n}\,, 0\,, 0\,, 0)\,, \bar{n}({\boldsymbol x }\,, t)\equiv \int {\rm d}^3 \bar{{\boldsymbol{p}}}\, f_{0} = \int_{0}^{\infty } {\rm d} \bar{p}\, 4 {\pi} \bar{p}^2 \, f_{0}({\boldsymbol x }\,, t\,, \bar{p}). \end{eqnarray*} (116) The first-order particle flux four-vector |$N^{\hat{a}}_{1}$| contains the diffusive flux as well as a drift term. Since the angular average of |$f_{1}(\tilde{{\boldsymbol Z }})$| satisfies ⟨f1⟩ = 0 we have |$N^{\hat{0}}_{1} = 0$|⁠, |$N_{1}^{\mu } U_{\mu } = 0$| in a general Lorentz frame. So |$N^{\hat{a}}_{1}$| can be written as \begin{eqnarray*} N^{\hat{a}}_{1} ({\boldsymbol x }\,, t)= \int_{0}^{\infty } {\rm d} \bar{p}\, \left(0\,, {\boldsymbol {\bar{S}} }({\boldsymbol x }\,, t\,, \bar{p}) \right) \,, \end{eqnarray*} (117) with \begin{eqnarray*} {\boldsymbol {\bar{S}} }({\boldsymbol x }\,, t\,, \bar{p}) \equiv \bar{p}^2 \, \int_{0}^{2 {\pi} } {\rm d} \bar{\phi } \, \int_{-1}^{1} {\rm d} \bar{\mu }\, {\boldsymbol {\bar{v}} } \, f_{1}({\boldsymbol x }\,, t\,, \bar{{\boldsymbol{p}}}) = 4 {\pi} \bar{p}^2 \, \left\langle {\boldsymbol {\bar{v}} } \, f_{1}({\boldsymbol x }\,, t\,, \bar{{\boldsymbol{p}}}) \right\rangle . \end{eqnarray*} (118) Using |${\boldsymbol {\bar{v}} } = \bar{v}\sqrt{1 - \bar{\mu }^2}\, (\cos {\bar{\phi }}\, {{\boldsymbol {\hat{x}}}}+ \sin {\bar{\phi }}\, {{\boldsymbol {\hat{y}}}}) + \bar{v}\bar{\mu }\, \mbox{$\boldsymbol{\hat{b}}$}$| in the local Cartesian coordinate system, substituting expression (108) for |$f_{1}({\boldsymbol x }\,, t\,, \bar{{\boldsymbol{p}}})$| and performing the integration over momentum solid angle one finds that only the first two terms in the expression for |$f_{1}({\boldsymbol x }\,, t\,, \bar{{\boldsymbol{p}}})$| contribute to the flux: The term involving |$\mbox{$\boldsymbol {\rm \hat{M}} $} {{\boldsymbol \, : \,}}\mbox{$\boldsymbol {\rm W} $}$| has an angular dependence such that it integrates to zero when performing the integration over |$\bar{\mu }$| and |$\bar{\phi }$|⁠. That term will become important when we discuss viscosity in Paper 2. 9.1 Diffusion The first term (⁠|$\propto {\boldsymbol {\tilde{{\nabla}}}}f_{0}$|⁠) in expression (108) gives the diffusive contribution to the flux. It takes a form equivalent to Fick’s law: \begin{eqnarray*} {\boldsymbol {\bar{S}} }_{\rm diff}({\boldsymbol x }\,, t\,, \bar{p}) = - 4 {\pi} \bar{p}^2 \, \left[ \, \mbox{$\boldsymbol {\bar{\rm { \mathit{ D}}}} $} {\boldsymbol{\cdot}}{\boldsymbol {\tilde{{\nabla}}}}f_{0}(\tilde{{\boldsymbol Z }}) \, \right]\,, \end{eqnarray*} (119) with the three-tensor |$\mbox{$\boldsymbol {\bar{\rm \mathit{ D}}} $}$| (diffusion tensor) defined in the FRF as: \begin{eqnarray*} \mbox{$\boldsymbol {\bar{\rm { \mathit{ D}}}} $} = \frac{\bar{v}^2 \bar{\tau }_{\rm r}}{3} \, \left(\begin{array}{@{}c@{\quad}c@{\quad}c@{}}\displaystyle \frac{1}{1 + \bar{\Omega }^2 \bar{\tau }_{\rm r}^2} & {\displaystyle - \frac{\bar{\Omega } \bar{\tau }_{\rm r}}{1 + \bar{\Omega }^2 \bar{\tau }_{\rm r}^2}} & 0 \\ & & \\ {\displaystyle \frac{\bar{\Omega } \bar{\tau }_{\rm r}}{1 + \bar{\Omega }^2 \bar{\tau }_{\rm r}^2}} & {\displaystyle \frac{1}{1 + \bar{\Omega }^2 \bar{\tau }_{\rm r}^2}} & 0 \\ & & \\ 0 & 0 & 1 \\ \end{array} \right) . \end{eqnarray*} (120) As in the non-relativistic case (e.g. Miyamoto 1980, ch. 7.3), diffusion in a magnetic field is characterized entirely by three coefficients for parallel, perpendicular, and cross diffusion (Hall diffusion) with respect to the plasma rest frame magnetic field. They are respectively given by \begin{eqnarray*} \bar{D}_{\parallel } = \frac{\bar{v}^2 \bar{\tau }_{\rm r}}{3}\,, \bar{D}_{\perp } = \frac{\bar{v}^2}{3} \, \frac{\bar{\tau }_{\rm r}}{1 + \bar{\Omega }^2 \bar{\tau }_{\rm r}^2}\,, \quad \bar{D}_{\times } = \frac{\bar{v}^2}{3} \, \frac{\bar{\Omega } \, \bar{\tau }_{\rm r}^2 }{1 + \bar{\Omega }^2 \bar{\tau }_{\rm r}^2} . \end{eqnarray*} (121) Note that these three diffusion coefficients are all defined entirely in terms of quantities in the FRF: |$\bar{v}$|⁠, |$\bar{\tau }_{\rm r}$|⁠, and |$\bar{\Omega }$|⁠. The relaxation time |$\bar{\tau }_{\rm r}$| appearing here was introduced in the Anderson–Witting collision operator (equation 35). While algebraically convenient that collision operator is not what one expects when charged particles are scattered by MHD turbulence in the FRF. In Section 11 below we briefly consider the case of a Fokker–Planck type relaxation operator describing diffusive elastic scattering in the FRF that resembles what one finds in the case of MHD turbulence. There we show that the results for diffusion and drift are formally the same as obtained here. In the viscosity calculation the results differ slightly, as explained in Paper 2. 9.2 Drift velocity The term |$\propto {\boldsymbol {\bar{\Sigma}}} {\boldsymbol{\cdot}}{\boldsymbol {\bar{a}} }$| in the expression for |$f_{1}({\boldsymbol x }\,, t\,, \bar{{\boldsymbol{p}}})$|⁠, due to the acceleration felt in the non-inertial FRF, leads to a drift term of the form \begin{eqnarray*} {\boldsymbol {\bar{S}} }_{\rm drift}({\boldsymbol x }\,, t\,, \bar{p}) = 4 {\pi} \bar{p}^2 \, {\boldsymbol {\bar{V}} }_{\rm D}({\boldsymbol x }\,, t\,, \bar{p}) \, f_{0}({\boldsymbol x }\,, t\,, \bar{p})\,, \end{eqnarray*} (122) with \begin{eqnarray*} {\boldsymbol {\bar{V}} }_{\rm D}({\boldsymbol x }\,, t\,, \bar{p}) = - \frac{1}{3 \bar{p}^2} \, \frac{\partial }{\partial \bar{p}} \left\lbrace \, \bar{p}^3 \bar{\tau }_{\rm r} \, \left[ \, \frac{{\boldsymbol {\bar{a}} }_{\perp }}{\displaystyle {1 + \bar{\Omega }^2 \bar{\tau }_{\rm r}^2}} + \frac{\bar{\Omega } \bar{\tau }_{\rm r}}{\displaystyle {1 + \bar{\Omega }^2 \bar{\tau }_{\rm r}^2}} \, \left({\boldsymbol {\bar{a}} } {\boldsymbol{\times}}\mbox{$\boldsymbol{\hat{b}}$}\right)+ (\mbox{$\boldsymbol{\hat{b}}$}{\boldsymbol{\cdot}}{\boldsymbol {\bar{a}} }) \, \mbox{$\boldsymbol{\hat{b}}$}\, \right]\, \right\rbrace . \end{eqnarray*} (123) Here |${\boldsymbol {\bar{a}} }_{\perp } \equiv {\boldsymbol {\bar{a}} } - (\mbox{$\boldsymbol{\hat{b}}$}{\boldsymbol{\cdot}}{\boldsymbol {\bar{a}} }) \, \mbox{$\boldsymbol{\hat{b}}$}$| is the component of the FRF three-vector |${\boldsymbol {\bar{a}} }$| perpendicular to the FRF magnetic field. This result is most easily obtained by (i) using expression (110) for |${\boldsymbol {\bar{\Sigma}}}$|⁠, (ii) employing that the integration over angles gives 4π times the average over solid angle, with |$<\bar{p}_{x}^2>= <\bar{p}_{y}^2>= <\bar{p}_{z}^{2}>= \bar{p}^2/3$|⁠, and (iii) followed by a partial integration over |$\bar{p}$|⁠. The drift term is small compared to the diffusion term in the non-relativistic limit, and is therefore absent in the results of Kulsrud (1983) and Skilling (1975): Formally it is of the order |$V^2 \, \bar{\tau }_{\rm r}/L$| with L the gradient scale of the flow. The diffusion term is typically of the order |$\bar{v}^2 \, \bar{\tau }_{\rm r}/L$|⁠, that is a factor |$(\bar{v}/V)^2 \gg 1$| larger. In the limit |$\bar{\Omega } \bar{\tau }_{\rm r}\gg 1$| the term involving |${\boldsymbol {\bar{a}} } {\boldsymbol{\times}}\mbox{$\boldsymbol{\hat{b}}$}$| has a simple physical interpretation: It is independent of the collisional relaxation time |$\bar{\tau }_{\rm r}$| and equals (reinstating c for clarity) \begin{eqnarray*} {\boldsymbol {\bar{V}} }_{\times } =- \frac{1}{3 \bar{p}^2} \, \frac{\partial }{\partial \bar{p}} \left[ \, \bar{p}^3 \left(\frac{c}{q B} (\bar{\gamma } m \, {\boldsymbol {\bar{a}} } {\boldsymbol{\times}}\mbox{$\boldsymbol{\hat{b}}$}) \right) \, \right]. \end{eqnarray*} (124) This is the effect of the well-known |${\boldsymbol F } {\boldsymbol{\times}}{\boldsymbol B }$| drift that particles experience in a strong magnetic field, where the force here is the average fictitious force |$- \bar{\gamma } m \, {\boldsymbol {\bar{a}} }$| due to the three-acceleration |${\boldsymbol {\bar{a}} }$| of the FRF. If one defines the average drift velocity as \begin{eqnarray*} {\boldsymbol {\bar{V}} }_{\rm drift} ({\boldsymbol x }\,, t)\equiv \frac{1}{\bar{n}} \, \int_{0}^{\infty } {\rm d} \bar{p}\, 4 {\pi} \bar{p}^2 \, f_{0}(\tilde{{\boldsymbol Z }}) \, {\boldsymbol {\bar{V}} }_{\rm D}({\boldsymbol x }\,, t\,, \bar{p}) \end{eqnarray*} (125) the FRF number flux four-vector is (putting ε = 1) \begin{eqnarray*} N^{\hat{a}} ({\boldsymbol x }\,, t)= \bar{n} \, \left(1\,, {\boldsymbol {\bar{V}} }_{\rm drift} \right) - \int_{0}^{\infty } {\rm d} \bar{p}\, 4 {\pi} \bar{p}^2 \, \left[ \, \mbox{$\boldsymbol {\bar{\rm \mathit{ D}}} $} {\boldsymbol{\cdot}}{\boldsymbol {\tilde{{\nabla}}}}f_{0}(\tilde{{\boldsymbol Z }}) \, \right]. \end{eqnarray*} (126) 10 COVARIANCE OF THE RESULTS One can always find the flux |$N^{\mu }({\boldsymbol x }\ ,, t\ ,, \bar{p})$| in the LF by using the Lorentz transformation from the FRF: |$N^{\mu } = {\rm E}^{\mu }_{\, \hat{a}} \, N^{\hat{a}}$|⁠. There is however a quicker and more general way that uses the fact that physical relations can always be expressed covariantly. The covariance requirement also serves as a check on the algebra of the previous sections that was done in a 3 + 1 split. In what follows we put the ordering parameter ε equal to unity. It is useful to introduce the two symmetric projection four-tensors hμν and Hμν, as well as the anti-symmetric tensor |$\hat{f}^{\mu \nu }$|⁠: \begin{eqnarray*} h^{\mu \nu } = \eta^{\mu \nu } - U^{\mu } U^{\nu }\,, H^{\mu \nu } = h^{\mu \nu } + b^{\mu } b^{\nu }\,, \quad \hat{f}^{\mu \nu } = e^{\mu \nu \sigma \tau }\, U_{\sigma } b_{\tau } = F^{\mu \nu }/{{\rm B}}. \end{eqnarray*} (127) Here eμνστ is the totally antisymmetric Levi (pseudo-)tensor. In the FRF, using the local xyz spatial coordinates introduced in Section 5, the first two tensors take the form |$h^{\hat{a}\hat{b}} = {\rm diag}(0\,, -1\,, -1 \,, -1)$|⁠, |$H^{\hat{a}\hat{b}} = {\rm diag}(0\,, -1\,, -1\,, 0)$|⁠. The only non-vanishing components of the anti-symmetric tensor |$\hat{f}^{\hat{a}\hat{b}}$| in the FRF are |$\hat{f}^{12} = - \hat{f}^{21} = - 1$|⁠, which implies that in the FRF \begin{eqnarray*} \hat{f}^{\hat{a}\hat{b}} \, A_{\hat{b}} = (0\,, \quad {\boldsymbol {\bar{A}} } {\boldsymbol{\times}}\mbox{$\boldsymbol{\hat{b}}$}) . \end{eqnarray*} (128) These tensors satisfy |$h^{\mu \nu } \, U_{\nu } = 0$|⁠, |$H^{\mu \nu } \, U_{\nu } = H^{\mu \nu } \, b_{\nu } = 0$| and |$\hat{f}^{\mu \nu } \, U_{\nu } = \hat{f}^{\mu \nu }\, b_{\nu } = 0$|⁠. It is possible to write the terms involving the vector |${\boldsymbol {\bar{\Sigma}}}$| in expression (108) covariantly. To that end we define the two projection tensors \begin{eqnarray*} G^{\mu \nu }_{\pm } \equiv \frac{H^{\mu \nu }}{\displaystyle {1 + \bar{\Omega }^2 \bar{\tau }_{\rm r}^2}} \pm \frac{ \bar{\Omega } \bar{\tau }_{\rm r}}{\displaystyle {1 + \bar{\Omega }^2 \bar{\tau }_{\rm r}^2}} \, \hat{f}^{\mu \nu } - b^{\mu } \, b^{\nu } . \end{eqnarray*} (129) The antisymmetry of fμν and the symmetry of the two other tensors in this definition imply that |$G_{-}^{\mu \nu } = G_{+}^{\nu \mu }$|⁠. The vector |${\boldsymbol {\bar{\Sigma}}}$| was defined in (110) using FRF momentum variables. We introduce the FRF four-vector |$\Sigma^{\hat{a}} \equiv (0\,, {\boldsymbol {\bar{\Sigma}}})$|⁠. The same vector is given as \begin{eqnarray*} \Sigma^{\mu } \equiv \frac{1}{\bar{p}} \, G_{-}^{\mu \nu } \, p_{\nu } \end{eqnarray*} (130) in an arbitrary Lorentz frame. One can express |$\bar{p}= | \bar{{\boldsymbol{p}}}|$| as the Lorentz-invariant quantity \begin{eqnarray*} \bar{p}= \sqrt{(p {\boldsymbol{\cdot}}U)^2 - (p {\boldsymbol{\cdot}}p)} . \end{eqnarray*} (131) Here we use the notation |$A {\boldsymbol{\cdot}}B \equiv A_{\mu } B^{\mu } = A^{0}B^{0} - {\boldsymbol A } {\boldsymbol{\cdot}}{\boldsymbol B }$| for the scalar product of two four-vectors. One has |$G_{\pm }^{\mu \nu } \, U_{\nu } = U_{\nu } \, G_{\pm }^{\mu \nu } = 0$| and |$G_{\pm }^{\mu \nu } \, b_{\nu } = b^{\mu }$|⁠, |$b_{\mu } \, G_{\pm }^{\mu \nu } = b^{\nu }$|⁠. When there is no magnetic field (or for uncharged particles), the relation |$G_{\pm }^{\mu \nu } = h^{\mu \nu }$| applies. The first two terms in expression (108) for |$f_{1}({\boldsymbol x }\,, t\,, \bar{{\boldsymbol{p}}})$|⁠, which respectively lead to the diffusive flux and the drift, can be written quasi-covariantly in terms of Σμ and the four-acceleration aμ = dUμ/dτ: \begin{eqnarray*} f_{1}^{\rm diff}({\boldsymbol x }\,, t\,, \bar{{\boldsymbol{p}}}) = - \bar{v}\bar{\tau }_{\rm r} \, \Sigma^{\mu } \nabla_{\mu } f_{0}\,, \quad f_{1}^{\rm drift}({\boldsymbol x }\,, t\,, \bar{{\boldsymbol{p}}}) = - \bar{\tau }_{\rm r} \, \left(\Sigma^{\mu } a_{\mu } \right) \, \bar{\gamma } m \, \frac{\partial f_{0}}{\partial \bar{p}}. \end{eqnarray*} (132) Using equation (34) it is easily seen that one can write the FRF diffusive flux as \begin{eqnarray*} N^{\hat{a}}_{\rm diff} = - 4 {\pi} \bar{p}^{2} \, \left[ \, - \bar{D}_{\perp } \, H^{\hat{a}\hat{b}} \, \frac{\partial f_{0}}{\partial x^{\hat{b}}} + \bar{D}_{\parallel } \, b^{\hat{a}} b^{\hat{b}} \, \frac{\partial f_{0}}{\partial x^{\hat{b}}} + \bar{D}_{\times } \, \hat{f}^{\hat{a}\hat{b}} \, \frac{\partial f_{0}}{\partial x^{\hat{b}}} \, \right] \equiv - 4 {\pi} \bar{p}^{2} \, D^{\hat{a}\hat{b}} \, \frac{\partial f_{0}}{\partial x^{\hat{b}}} . \end{eqnarray*} (133) This defines a diffusion tensor |$D^{\hat{a}\hat{b}}\equiv - \bar{D}_{\perp } \, H^{\hat{a}\hat{b}} + \bar{D}_{\parallel } \, b^{\hat{a}} b^{\hat{b}} + \bar{D}_{\times } \, \hat{f}^{\hat{a}\hat{b}} = - \bar{D} \, G_{-}^{\hat{a}\hat{b}}$|⁠, using definition (129) for |$G_{-}^{\mu \nu }$|⁠. The scalar diffusion coefficient |$\bar{D}$| that appears in this definition equals |$\bar{D} \equiv \bar{v}^2 \, \bar{\tau }_{\rm r}/3$|⁠. This tensor has the FRF components |$D^{\hat{0}\hat{0}} = D^{\hat{0} \hat{i}} = D^{\hat{i} \hat{0}} = 0$|⁠, |$D^{\hat{i} \hat{j}} = (\bar{\mbox{$\boldsymbol {\rm D} $}})_{ij}$|⁠. Since relation (133) is formally covariant, the corresponding relation \begin{eqnarray*} N^{\mu }_{\rm diff} = - 4 {\pi} \bar{p}^2 \, D^{\mu \nu } \, \frac{\partial f_{0}({\boldsymbol x }\,, t\,, \bar{p})}{\partial x^{\nu }} \end{eqnarray*} (134) holds in any frame, with Dμν defined as \begin{eqnarray*} D^{\mu \nu } \equiv - \bar{D}_{\perp } \, H^{\mu \nu } + \bar{D}_{\parallel } \, b^{\mu } b^{\nu } + \bar{D}_{\times } \, \hat{f}^{\mu \nu } = - \bar{D} \, G_{-}^{\mu \nu } . \end{eqnarray*} (135) It should be noted that we continue to use FRF momentum variables. The drift term can similarly be cast into a covariant form. Relation (123) in the FRF corresponds to \begin{eqnarray*} V^{\hat{a}}_{\rm drift} = - \frac{1}{3 \bar{p}^2} \, \frac{\partial }{\partial \bar{p}} \left[ \, \bar{p}^3 \bar{\tau }_{\rm r} \, \left\lbrace \frac{H^{\hat{a}\hat{b}} \, a_{\hat{b}}}{\displaystyle {1 + \bar{\Omega }^2 \bar{\tau }_{\rm r}^2}} + \frac{\bar{\Omega } \bar{\tau }_{\rm r}}{\displaystyle {1 + \bar{\Omega }^2 \bar{\tau }_{\rm r}^2}} \, \left(\hat{f}^{\hat{a}\hat{b}} \, a_{\hat{b}} \right) - (b \cdot a) \, b^{\hat{a}} \right\rbrace \, \right]. \end{eqnarray*} (136) The FRF fluid four-acceleration is |$a_{\hat{a}} = (0\ ,, - {\boldsymbol {\bar{a}} })$| (see Appendix C, equation 107). We also use |$b {\boldsymbol{\cdot}}a \equiv b^{\hat{a}} a_{\hat{a}} = - \mbox{$\boldsymbol{\hat{b}}$}{\boldsymbol{\cdot}}{\boldsymbol {\bar{a}} }$|⁠. This implies that one can write the drift velocity in any frame covariantly in terms of the four-acceleration aμ = dUμ/dτ and the projection tensor |$G_{+}^{\mu \nu }$| as \begin{eqnarray*} V^{\mu }_{\rm drift} = - \frac{1}{3 \bar{p}^2} \, \frac{\partial }{\partial \bar{p}} \left(\bar{p}^3 \bar{\tau }_{\rm r} \, G_{+}^{\mu \nu } a_{\nu } \right) . \end{eqnarray*} (137) Note that |$V^{\mu } \, U_{\mu } = 0$|⁠. The fact that |$G_{+}^{\mu \nu }$| appears here while |$G_{-}^{\mu \nu }$| appears in the definition of Σμ and Dμν is the result of isotropy of f0 so that |$\langle p^{\hat{i}} \, p^{\hat{j}}\rangle = \bar{p}^2 \, \delta_{\hat{i}\hat{j}}/3$| and the properties of the cross-product in |$f^{\hat{i}\hat{j}} \, p_{\hat{j}} = \bar{{\boldsymbol{p}}}{\boldsymbol{\times}}\mbox{$\boldsymbol{\hat{b}}$}$|⁠. This leads to the following FRF relation, needed for the solid-angle average of the Hall term |$\propto \bar{\Omega } \bar{\tau }_{\rm r}/(\displaystyle {1 + \bar{\Omega }^2 \bar{\tau }_{\rm r}^2})$| in |$V^{\mu }_{\rm drift}$|⁠: \begin{eqnarray*} \left\langle \left[ (\mbox{$\boldsymbol{\hat{b}}$}{\boldsymbol{\times}}\bar{{\boldsymbol{p}}}) {\boldsymbol{\cdot}}{\boldsymbol {\bar{a}} } \, \right] \bar{{\boldsymbol{p}}}\right\rangle^{\hat{i}} = \epsilon_{\hat{k}\hat{l} \hat{m}} \, b^{\hat{l}} \, \bar{a}^{\hat{k}} \, \left\langle p^{\hat{m}} \, p^{\hat{i}} \right\rangle = \frac{\bar{p}^2}{3} \, \epsilon_{\hat{k} \, \hat{l} \, \hat{i}} \, b^{\hat{l}} \, \bar{a}^{\hat{k}} = \frac{\bar{p}^2}{3} \, ({\boldsymbol {\bar{a}} } {\boldsymbol{\times}}\mbox{$\boldsymbol{\hat{b}}$})^{\hat{i}} . \end{eqnarray*} (138) 10.1 The unmagnetized case If there is no magnetic field or for uncharged particles (in particular photons) scattered by the medium one has |$\bar{D}_{\perp } = \bar{D}_{\parallel } = \bar{v}^2 \, \bar{\tau }_{\rm r}/3 = \bar{D}$|⁠, |$\bar{D}_{\times } = 0$|⁠, and |$G_{\pm }^{\mu \nu } = h^{\mu \nu }$|⁠. Then |$D^{\mu \nu }= - \bar{D} \, h^{\mu \nu }$| and relation (134) for the diffusive part of the particle flux then simplifies to the covariant version of Fick’s law: \begin{eqnarray*} N^{\mu }_{\rm diff} = 4 {\pi} \bar{p}^2 \, \bar{D} \, h^{\mu \nu } \, \frac{\partial f_{0}({\boldsymbol x }\,, t\,, \bar{p})}{\partial x^{\nu }} \, \mbox{(for ${{\rm B}}= 0$ or $q = 0$)} . \end{eqnarray*} (139) The drift term also simplifies as the drift velocity now equals, with Uμaμ= 0, so that |$h^{\mu \nu } \, a_{\mu } = a^{\mu }$|⁠: \begin{eqnarray*} V^{\mu }_{\rm drift} = - \frac{1}{3 \bar{p}^2} \, \frac{\partial }{\partial \bar{p}} \left(\bar{p}^3 \bar{\tau }_{\rm r} \, a^{\mu } \right) \, \mbox{(for ${{\rm B}}= 0$ or $q = 0$)} . \end{eqnarray*} (140) 10.2 Conservation form of the transport equation The number flux four-vector per unit momentum (including the contribution from the zero-order term |$N^{\hat{a}}_{0}$|⁠, equation 116) can be written in terms of |${\cal N}({\boldsymbol x }\,, t\,, \bar{p}) \equiv 4 {\pi} \bar{p}^2 \, f_{0}({\boldsymbol x }\,, t\,, \bar{p}) = {\rm d}{\bar{n}}/{{\rm d} \bar{p}}$| as \begin{eqnarray*} N^{\mu }({\boldsymbol x }\,, t\,, \bar{p}) = \left(U^{\mu } + V^{\mu }_{\rm drift} \right) \, {\cal N} - D^{\mu \nu } \, \frac{\partial {\cal N}}{\partial x^{\nu }} . \end{eqnarray*} (141) Because |$U_{\mu } \, V^{\mu }_{\rm drift} = 0$| and |$U_{\mu } \, {\cal D}^{\mu \nu } = 0$| one has \begin{eqnarray*} N^{\mu } U_{\mu } ={\cal N}({\boldsymbol x }\,, t\,, \bar{p}) = \frac{{\rm d} \bar{n}}{{\rm d} \bar{p}} \,, \end{eqnarray*} (142) so the drift and diffusion terms in |$N^{\mu }({\boldsymbol x }\,, t\,, \bar{p})$| do not change the ‘proper density per unit momentum’ of the particles. We have shown in Section 4 that Liouville’s theorem holds for mixed variables. One can average the Vlasov equation in conservation form over solid angle. With |$f_{0} \equiv \langle f(\tilde{{\boldsymbol Z }})\rangle$| this leads to \begin{eqnarray*} \left\langle (\partial /\partial \tilde{{\boldsymbol Z }}) {\boldsymbol{\cdot}}[\dot{\tilde{{\boldsymbol Z }}} \, f(\tilde{{\boldsymbol Z }})] \right\rangle = \left\langle {\rm St}(f) \right\rangle \equiv - \frac{\langle f\rangle - \langle f\rangle}{\bar{\tau }_{\rm r}} = 0 . \end{eqnarray*} (143) After multiplication of this relation by |$4 {\pi} \bar{p}^2$| one finds an equation for |${\cal N}({\boldsymbol x }\,, t\,, \bar{p})$|⁠: \begin{eqnarray*} \nabla_{\mu } \left\lbrace \left(U^{\mu } + V^{\mu }_{\rm drift} \right) \, {\cal N} - D^{\mu \nu } \, \nabla_{\nu } {\cal N} \right\rbrace - \frac{\partial }{\partial \bar{p}} \left(\frac{\bar{p}}{3} \, (\nabla {\boldsymbol{\cdot}}U) \, {\cal N} \right) + \frac{\partial }{\partial \bar{p}} \left(4 {\pi} \bar{p}^2 \, \left\langle \tilde{K}(\tilde{{\boldsymbol Z }}) \, f_{1}(\tilde{{\boldsymbol Z }}) \right\rangle \right) = 0 . \end{eqnarray*} (144) This is the transport equation in conservation form in the diffusion approximation. Here ∇μ = ∂/∂xμ with FRF momentum fixed, and we have used |${\rm d}/{\rm d} \tau = U {\boldsymbol{\cdot}}\nabla$| as well as relations (94) and (95). We have put ε = 1 in the last term on the left-hand side. The quantity |$\tilde{K}(\tilde{{\boldsymbol Z }})$| was defined in equation (99). The last term on the left-hand side of equation (144) represents viscous heating. This term is usually absent in calculations for non-relativistic flows since it is (by assumption) small in these calculations. The works of Jokipii & Williams (1992) and Williams & Jokipii (1993) are the exception: They calculate this term for a non-relativistic flow. In Paper 2 we evaluate this term explicitly. There we show that it takes the form of momentum diffusion in |$\bar{p}$| and that it indeed represents the irreversible heating of the scattered particles in a thermodynamic sense. As already mentioned the second term (⁠|$\propto \nabla {\boldsymbol{\cdot}}U$|⁠) in equation (144) corresponds to adiabatic momentum gains/losses (see equation 95). Equation (144) contains all terms formally of order ε: spatial drift and spatial diffusion as well as the momentum diffusion (heating) term |$\propto \left\langle \tilde{K}(\tilde{{\boldsymbol Z }}) \, f_{1}(\tilde{{\boldsymbol Z }}) \right\rangle$|⁠, consistent with the diffusion approximation. It shows for the first time that the transport equation for particles scattered in a relativistic flow is an almost straightforward generalization of the non-relativistic result in the case of simple elastic scattering in the FRF. The main differences are [1] the appearance of the drift term |$\propto V^{\mu }_{\rm drift}$|⁠, and [2] the consistent use of the FRF momentum rather than the LF momentum so that the spatial derivatives are taken with the FRF four-momentum fixed. In a non-relativistic flow, with |$| {\boldsymbol{U}}| \ll 1$| and |$V/\bar{v}\ll 1$|⁠, one can neglect the drift term |$\propto V^{\mu }_{\rm drift}$| and (as long as one works to first order in |$V/\bar{v}$|⁠) one can replace the FRF momentum |$\bar{p}$| by the LF momentum p. Then equation (144) is equivalent to the results of Skilling (1975), without the momentum diffusion terms due to acceleration by Alfvén waves that are treated in that paper but that do not appear here since we assumed elastic scattering in the FRF. As discussed we do calculate the additional effects of the diffusion across the large-scale magnetic field and (in Paper 2) of irreversible heating. The simple form (equation 144) of the transport equation is implicit in the results of Webb (1985), but there it is disguised by his 3 + 1 split in the notation in terms of |${\boldsymbol x }$|⁠, t (and the corresponding derivatives), Γ, and |${\boldsymbol{V}}$|⁠, and the different (more complicated) choice for the scattering operator. Webb also uses gyro-averaging in most of his results, which eliminates the terms associated with cross-field diffusion and possible viscous effects from consideration, as argued above. 11 THE EFFECT OF AN ALTERNATIVE CHOICE FOR THE SCATTERING OPERATOR In the calculation of |$f_{1}({\boldsymbol x }\ ,, t\ ,, \bar{{\boldsymbol{p}}})$|⁠, which involves solving equation (100), the precise form of the relaxation operator St(f) does influence results. So far we have chosen to employ the algebraic Anderson–Witting form, the relativistic generalization of the Krook collision term. This operator damps all anisotropies (dipole, quadrupole, etc.) at the same rate determined by the relaxation time |$\bar{\tau }_{\rm r}$|⁠. In this section we briefly discuss the effect of employing an alternative choice: a Fokker–Plank type collision operator. In terms of the pitch angle |$\bar{\theta }= \cos^{-1} \bar{\mu }$| and the gyro-phase angle |$\bar{\phi }$| one has in terms of |$\overline{\rm St}(f) = \bar{\gamma }^{-1} \, {\rm St}(f)$| \begin{eqnarray*} \overline{{\rm St}}_{\rm FP}(f) \equiv \frac{1}{ 2 \bar{\tau }_{\rm r}}\, \left[ \, \frac{1}{\sin \bar{\theta }}\, \frac{\partial }{\partial \bar{\theta }} \left(\sin \bar{\theta }\, \frac{\partial f}{\partial \bar{\theta }} \right) + \frac{1}{\sin^2 \bar{\theta }} \frac{\partial^2 f}{\partial \bar{\phi }^2} \, \right], \end{eqnarray*} (145) which corresponds to isotropic diffusion of the momentum direction in the fluid frame. The factor |$\bar{\gamma } = (p {\boldsymbol{\cdot}}U)/m$| takes care of time dilation and the factor |${1}/{2}$| is a convention so that the relaxation time for dipole anisotropies equals |$\bar{\tau }_{\rm r}$|⁠. The equivalent of equation (100) becomes \begin{eqnarray*} \bar{\Omega } \, \frac{\partial f_{1}}{\partial \bar{\phi }} + \frac{1}{ 2 \bar{\tau }_{\rm r}}\, \left[ \, \frac{1}{\sin \bar{\theta }}\, \frac{\partial }{\partial \bar{\theta }} \left(\sin \bar{\theta }\, \frac{\partial f}{\partial \bar{\theta }} \right) + \frac{1}{\sin^2 \bar{\theta }} \frac{\partial^2 f}{\partial \bar{\phi }^2} \, \right] = \mbox{source term}, \end{eqnarray*} (146) where the source term is the entire term on the right-hand side of equation (100). The obvious approach in this case is to use spherical harmonics |$Y_{\ell }^{m}(\bar{\theta }\ ,, \, \bar{\phi })$|⁠, the eigenfunctions of the diffusive relaxation operator |$\overline{{\rm St}}_{\rm FP}(f)$| that satisfy \begin{eqnarray*} \overline{{\rm St}}_{\rm FP}\left(Y_{\ell }^{m}(\bar{\theta }, \bar{\phi }) \right) = - \frac{\ell (\ell +1)}{2 \bar{\tau }_{\rm r}} \, Y_{\ell }^{m}(\bar{\theta }, \bar{\phi }) . \end{eqnarray*} (147) Here |$\ell = 0\,, 1\,, 2 \,, \ldots$| and −ℓ ≤ m ≤ +ℓ. In what follows we employ the form of |$Y_{\ell }^{m}(\bar{\theta }\ ,, \, \bar{\phi })$| that includes the Condon–Shortley phase factor (c.f. Arfken & Weber 2005, ch.12.6): \begin{eqnarray*} Y_{\ell }^{m}(\bar{\theta }, \bar{\phi }) = (-1)^m \, \sqrt{\displaystyle \frac{2 \ell +1}{4 {\pi} } \, \frac{(\ell - m)!}{(\ell + m)!}} \, P_{\ell }^{\, m}(\cos \bar{\theta }) \, e^{im \bar{\phi }} . \end{eqnarray*} (148) The |$P_{\ell }^{\, m}(\cos \bar{\theta })$| are the associated Legendre polynomials. We expand |$f_{1}({\boldsymbol x }\,, t\,, \bar{{\boldsymbol{p}}})$| in spherical harmonics, \begin{eqnarray*} f_{1}({\boldsymbol x }\,, t\,, \bar{{\boldsymbol{p}}}) = \sum_{\ell = 1}^{\infty } \, \sum_{m = - \ell }^{m = + \ell } a_{\ell }^{\, m}({\boldsymbol x }\,, t\,, \bar{p}) \, Y_{\ell }^{m}(\bar{\theta }, \bar{\phi })\,, \end{eqnarray*} (149) which implies that (as before) |$\left\langle f_{1}({\boldsymbol x }\,, t\,, \bar{{\boldsymbol{p}}}) \right\rangle = 0$|⁠. In a similar fashion we expand the (known) source term on the right-hand side of equation (146), which by definition has a vanishing momentum space solid-angle average: \begin{eqnarray*} \mbox{source term} = \sum_{\ell = 1}^{\infty } \, \sum_{m = - \ell }^{m = + \ell } s_{\ell }^{\, m}({\boldsymbol x }\,, t\,, \bar{p}) \, Y_{\ell }^{m}(\bar{\theta }, \bar{\phi }) . \end{eqnarray*} (150) After substitution of these two expansions into equation (146) one finds the following relation that determines the coefficients |$a_{\ell }^{\, m}$| in terms of the known source-term coefficients |$s_{\ell }^{\, m}$|⁠: \begin{eqnarray*} \left(im \, \bar{\Omega } - \frac{\ell (\ell +1)}{2 \bar{\tau }_{\rm r}} \right) \, a_{\ell }^{\, m} = s_{\ell }^{\, m} . \end{eqnarray*} (151) This follows from the fact that the spherical harmonics are orthogonal eigenfunctions of the Fokker–Planck relaxation operator. Using the identities \begin{eqnarray*} \left(\begin{array}{c}\bar{\mu }\\ \\ \sqrt{1 - \bar{\mu }^2}\, \cos \bar{\phi }\\ \\ \sqrt{1 - \bar{\mu }^2}\, \sin \bar{\phi }\\ \end{array} \right) = \sqrt{\frac{4 {\pi} }{3}} \, \left(\begin{array}{c}Y_{1}^{\, 0}(\bar{\theta }, \bar{\phi }) \\ \\ \frac{1}{\sqrt{2}} \, \left(Y_{1}^{\, -1}(\bar{\theta }, \bar{\phi })- Y_{1}^{\, 1}(\bar{\theta }, \bar{\phi })\right) \\ \\ \frac{i}{\sqrt{2}} \, \left(Y_{1}^{\, -1}(\bar{\theta }, \bar{\phi })+ Y_{1}^{\, 1}(\bar{\theta }, \bar{\phi })\right) \\ \end{array} \right) \end{eqnarray*} (152) it is easily established that the terms in equation (100) that lead to the first two terms in expression (108) for |$f_{1}({\boldsymbol x }\,, t\,, \bar{{\boldsymbol{p}}})$|⁠, the terms that involve the three-vector |${\boldsymbol {\bar{\Sigma}}}$| defined in (109), all come exclusively from the ℓ = 1 terms (dipole terms) in expansion (150). Solving (151) for |$a_{1}^{0}$| and |$a_{1}^{\pm 1}$| one finds that the resulting contributions to |$f_{1}({\boldsymbol x }\,, t\,, \bar{{\boldsymbol{p}}})$| in the case of Fokker–Planck type scattering are unchanged with respect to the result found using the simpler Anderson–Witting relaxation term, i.e. |${\boldsymbol {\bar{\Sigma}}}_{\rm FP} = {\boldsymbol {\bar{\Sigma}}}$|⁠. This result is essentially due to the fact that the Fokker–Planck relaxation operator has the eigenvalue |$-1/\bar{\tau }_{\rm r}$| for the spherical harmonics with ℓ = 1, just as in the Anderson–Witting case. The remaining terms in the source term in equation (100) are the terms that involve the strain tensor |$\mbox{$\boldsymbol {\rm W} $}$| via the replacement rule (equation 103). These terms consist exclusively of terms in the expansion from spherical harmonics with ℓ = 2: They are quadrupole terms. This follows from the identities \begin{eqnarray*} \left(\begin{array}{c}\displaystyle \frac{3 \bar{\mu }^2 -1}{2} \\ \\ \bar{\mu }\sqrt{1 - \bar{\mu }^2}\, \cos \bar{\phi }\\ \\ \bar{\mu }\sqrt{1 - \bar{\mu }^2}\, \sin \bar{\phi }\\ \end{array} \right) = \sqrt{\frac{4 {\pi} }{5}} \, \left(\begin{array}{c}Y_{2}^{\, 0}(\bar{\theta }, \bar{\phi }) \\ \\ \frac{1}{\sqrt{6}} \, \left(Y_{2}^{\, - 1}(\bar{\theta }, \bar{\phi }) - Y_{2}^{\, 1}(\bar{\theta }, \bar{\phi })\right) \\ \\ \frac{i}{\sqrt{6}} \, \left(Y_{2}^{\, -1}(\bar{\theta }, \bar{\phi })+ Y_{2}^{\, - 1}(\bar{\theta }, \bar{\phi })\right) \\ \end{array} \right) \end{eqnarray*} (153) and \begin{eqnarray*} \left(1 - \bar{\mu }^2 \right) \, \left(\begin{array}{c}\cos 2 \bar{\phi }\\ \\ \sin 2 \bar{\phi }\\ \end{array} \right) = \sqrt{\frac{8 {\pi} }{3}} \, \left(\begin{array}{c}Y_{2}^{\, 2}(\bar{\theta }, \bar{\phi }) + Y_{2}^{\, - 2}(\bar{\theta }, \bar{\phi }) \\ \\ i \, \left(Y_{2}^{\, 2}(\bar{\theta }, \bar{\phi })- Y_{2}^{\, - 2}(\bar{\theta }, \bar{\phi })\right) \\ \end{array} \right) . \end{eqnarray*} (154) The eigenvalue for the ℓ = 2 under diffusive scattering is |$-3/\bar{\tau }_{\rm r}$|⁠, corresponding to a relaxation time |$\bar{\tau }_{\rm r}/3$|⁠: The Fokker–Planck relaxation term damps the ℓ = 2 spherical harmonics (quadrupole terms) three times faster than the Anderson–Witting relaxation term, where all terms with ℓ ≠ 0 have the same relaxation time, equal to |$\bar{\tau }_{\rm r}$|⁠. As a result, when dealing with Fokker–Planck (diffusive) scattering one must replace |$\bar{\tau }_{\rm r}$| by |$\bar{\tau }_{\rm r}/3$| in the terms in quadrupole terms in |$f_{1}({\boldsymbol x }\,, t\,, \bar{{\boldsymbol{p}}})$|⁠. The last term in expression (108), the term that governs viscous effects as we will show in Paper 2, becomes \begin{eqnarray*} \frac{\bar{\tau }_{\rm r}}{3} \left(\mbox{$\boldsymbol {\rm \hat{M}} $}_{\rm FP} {{\boldsymbol \, : \,}}\mbox{$\boldsymbol {\rm W} $} \right) \, \bar{p}\, \frac{\partial f_{0}}{\partial \bar{p}} . \end{eqnarray*} (155) Here the 3 × 3 matrix |$\mbox{$\boldsymbol {\rm \hat{M}} $}_{\rm FP}$| equals (analogous to definition 111 for |$\mbox{$\boldsymbol {\rm \hat{M}} $}$| in the Anderson–Witting relaxation) \begin{eqnarray*} \mbox{$\boldsymbol {\rm \hat{M}} $}_{\rm FP} \equiv \frac{1 - \bar{\mu }^2}{2} \, \frac{\mbox{$\boldsymbol {\rm M} $}_{1}- \frac{2}{3}\, \bar{\Omega } \, \bar{\tau }_{\rm r} \, \mbox{$\boldsymbol {\bar{\rm M}} $}_{1}}{\displaystyle {1 + \frac{4}{9} \, \bar{\Omega }^2 \bar{\tau }_{\rm r}^2}} + \frac{\bar{\mu } \, \sqrt{1 - \bar{\mu }^2}}{2} \, \frac{\mbox{$\boldsymbol {\rm M} $}_{2}- \frac{1}{3}\, \bar{\Omega } \, \bar{\tau }_{\rm r} \, \mbox{$\boldsymbol {\bar{\rm M}} $}_{2}}{\displaystyle {1 + \frac{1}{9} \,\bar{\Omega }^2 \bar{\tau }_{\rm r}^2}} + \frac{3 \bar{\mu }^2 - 1}{4} \, \mbox{$\boldsymbol{\hat{b}}$}\, \mbox{$\boldsymbol{\hat{b}}$}. \end{eqnarray*} (156) Otherwise, the functional form of |$f_{1}({\boldsymbol x }\,, t\,, \bar{{\boldsymbol{p}}})$| is unchanged. We therefore find in the Fokker–Planck case \begin{eqnarray*} f_{1 \rm FP}({\boldsymbol x }\,, t\,, \bar{{\boldsymbol{p}}}) = - \bar{\tau }_{\rm r} \, \left(\bar{v}\, {\boldsymbol {\bar{\Sigma}}} {\boldsymbol{\cdot}}{\boldsymbol {\tilde{{\nabla}}}}\right) f_{0} + \bar{\tau }_{\rm r} \, \left({\boldsymbol {\bar{\Sigma}}} {\boldsymbol{\cdot}}{\boldsymbol {\bar{a}} } \right) \, \bar{\gamma } m \, \frac{\partial f_{0}}{\partial \bar{p}} + \frac{\bar{\tau }_{\rm r}}{3} \, \left(\mbox{$\boldsymbol {\rm \hat{M}} $}_{\rm FP} {{\boldsymbol \, : \,}}\mbox{$\boldsymbol {\rm W} $} \right) \, \bar{p}\, \frac{\partial f_{0}}{\partial \bar{p}} \,, \end{eqnarray*} (157) with the vector |${\boldsymbol {\bar{\Sigma}}}$| still given by equation (109). All results of the previous sections are now easily transcribed to this case. In particular, the results obtained for diffusion and drift in Section 9 are unchanged. 12 DISCUSSION In this paper we have considered the transport of particles that are being scattered by a relativistic flow. We rederive and extend the results of Webb (1985). In Section 3 we gave a more intuitive derivation of the fictitious four-force that occurs if the fluid rest frame is not an inertial (Lorentz) frame. This fictitious force is the source of all terms associated with collisional transport: diffusion, drift, and viscosity. In Section 7 we explicitly showed the correspondence with Webb’s results. In Section 4 we demonstrated that it is possible to write all final equations in a quasi-covariant form and proved that Liouville’s theorem, i.e. incompressible particle flow in seven-dimensional phase space (⁠|$x^{\mu }\,, \bar{{\boldsymbol{p}}}$|⁠), still holds if one uses mixed laboratory frame/fluid rest frame variables. We generalized the approaches of Skilling (1975) and Kulsrud (1983) to a relativistic flow, showing in the process (Sections 5 and 6) that these two approaches are largely equivalent. In Section 7 we extended existing results by allowing a dependence on gyro-phase, which opens the possibility of discussing the effect of the large-scale magnetic field on particle transport properties. We derived in Section 8 that the perturbation expansion for |$f({\boldsymbol x }\,, t\,, \bar{{\boldsymbol{p}}}) = f(\tilde{{\boldsymbol Z }})$| that one employs in the diffusion approximation with the Anderson–Witting relaxation term |${\rm St}(f) = - \bar{\gamma }(f - f_{0})/\bar{\tau }_{\rm r}$| should have f0 = ⟨f⟩ with |$\langle f\rangle({\boldsymbol x }\,, t\,, \bar{p})$| the angular average of the exact Vlasov distribution. The perturbation expansion therefore takes the form \begin{eqnarray*} f(\tilde{{\boldsymbol Z }}) = \langle f (\tilde{{\boldsymbol Z }}) \rangle + \varepsilon \, f_{1}({\boldsymbol x }\,, t\,, \bar{{\boldsymbol{p}}}) + \ldots . \end{eqnarray*} (158) In particular, this is true for the case of pure scattering. This implies |$\langle f_{n}(\tilde{{\boldsymbol Z }})\rangle = 0$| for n ≥ 1. It is perhaps good to point out once again that this particular conclusion is independent of the form of the relaxation term St(f), as long as this term corresponds to simple elastic scattering in the FRF. Our use of the algebraic Anderson–Witting scattering term in Sections 4–10 yields the conditions f0 = ⟨f⟩ and |$\langle f_{1}(\tilde{{\boldsymbol Z }})\rangle = 0$| in an almost trivial manner. In Section 11 we discussed an alternative choice for the relaxation term: diffusive scattering in the form of the Fokker–Planck scattering operator. Then the same conclusion holds because of the properties (in particular |$\langle Y_{\ell }^{\, m}\rangle = 0$| for ℓ ≠ 0) and the orthogonality of the spherical harmonics used in the perturbation expansion that leads to the final expression (equation 157) for |$f_{1}({\boldsymbol x }\,, t\,, \bar{{\boldsymbol{p}}})$|⁠. In Section 9 we derived – as a first application – the cross-field diffusion of the scattered particles by explicitly calculating the gyro-phase-dependent terms in the distribution function |$f({\boldsymbol x }\,, t\,, \bar{{\boldsymbol{p}}})$| in a perturbation scheme valid for small anisotropies in the fluid rest frame (i.e. for frequent scattering). For relativistic flow speeds there is an extra drift term that does not appear in the non-relativistic treatment. In the follow-up paper (Paper 2) we show that this argument has consequences for the calculation of the energy-momentum tensor in the diffusion approximation. In particular it implies that the viscous terms in the energy-momentum tensor correspond to the effect of shear viscosity, and that there is no contribution from a bulk viscosity. We will show in Paper 2, by explicitly calculating the term on the right-hand side of equation (144), that the back-reaction of the flow on the scattered particles (or photons) leads to an irreversible heating term in the transport equation. In a second follow-up paper (Paper 3) we explore radiation transport in the diffusion approximation and the radiation stress tensor in a relativistic flow. Published results (e.g. Weinberg 1971; Mihalas & Weibel Mihalas 1984; and Castor 2004), derived for the case |$| {\boldsymbol{U}}| \ll 1$| (⁠|$| {\boldsymbol{V}}| \ll c$|⁠) give diffusive corrections to the comoving number – and energy density of the photons. In Paper 3 we argue that the condition |$\langle f_{n}(\tilde{{\boldsymbol Z }})\rangle = 0$| for n ≥ 1 implies that these corrections should not be included: The leading-order solution of the self-consistent radiation transfer equation in the diffusion approximation shows that these extra terms will vanish to (at least) second order in |$\varepsilon \sim \bar{\lambda }/L \ll 1$| with |$\bar{\lambda } = \bar{v}\bar{\tau }_{\rm r}$| the FRF scattering mean free path and L the gradient scale in the plasma. Our approach, which identifies |$f_{0}({\boldsymbol x }\,, t\,, \bar{p})$| with |$\langle f(\tilde{{\boldsymbol Z }})\rangle$| by including the leading-order transport equation in the expansion for |$f(\tilde{{\boldsymbol Z }})$|⁠, avoids this issue entirely, and is the mathematically correct procedure. These considerations also mean that our result for the viscous stresses from the scattered particles, when applied to photons as discussed in Paper 3, are different from the recent results of Coughlin & Begelman (2014). We demonstrate that for a gas of photons only shear viscosity (and therefore no bulk viscosity) is present, in agreement with the results of Weinberg (1971). ACKNOWLEDGEMENTS We thank: M.Begelman, A. Bell, R. Blandford, J. Kirk, the late D. Lynden-Bell, M. Rees and G. Peletier for stimulating discussions. This work was done in Nijmegen, Oxford, and Cambridge. We thank P. Groot, G. Nelemans, S. Balbus, and A. Fabian for the hospitality and support. REFERENCES Achterberg A. , 1983 , Phys. Rev. A , 28 , 2449 https://doi.org/10.1103/PhysRevA.28.2449 Crossref Search ADS Anderson J. L. , Witting H. R. , 1974 , Physica , 74 , 466 https://doi.org/10.1016/0031-8914(74)90355-3 Crossref Search ADS Arfken G. B. , Weber H. J. , 2005 , Mathematical Methods for Physicists , 6th ed. Elsevier , Amsterdam, Boston Google Preview WorldCat COPAC Bekenstein J. D. , Oron E. , 1978 , Phys. Rev. D , 18 , 1809 https://doi.org/10.1103/PhysRevD.18.1809 Crossref Search ADS Bhatnagar P. L. , Gross E. P. , Krook M. , 1954 , Phys. Rev. , 94 , 511 https://doi.org/10.1103/PhysRev.94.511 Crossref Search ADS Braginskii S. I. , 1965 , Reviews of Plasma Physics , 1 , 205 Carroll S. M. , 2004 , Spacetime and Geometry. An Introduction to General Relativity . Addison-Wesley , San Francisco Google Preview WorldCat COPAC Cary J. R. , Littlejohn R. G. , 1983 , Ann. Phys. NY , 151 , 1 https://doi.org/10.1016/0003-4916(83)90313-5 Crossref Search ADS Castor J. I. , 2014 , PubMed PubMed Radiation Hydrodynamics . Cambridge Univ. Press , Cambridge PubMed PubMed Google Preview WorldCat COPAC Coughlin E. R. , Begelman M. C. , 2014 , ApJ , 797 , 103 https://doi.org/10.1088/0004-637X/797/2/103 Crossref Search ADS Debbasch F. , van Leeuwen W. A. , 2009a , Physica A , 388 , 1079 https://doi.org/10.1016/j.physa.2008.12.023 Crossref Search ADS Debbasch F. , van Leeuwen W. A. , 2009b , Physica A , 388 , 1818 https://doi.org/10.1016/j.physa.2009.01.009 Crossref Search ADS Dolginov A. Z. , Toptygin I. M. , 1967 , Geomagn. Aeronomy , 7 , 785 Gleeson L. J. , Axford W. I. , 1967 , ApJ , 149 , L115 https://doi.org/10.1086/180070 Crossref Search ADS Goldstein H. , Poole C. , Safko J. , 2002 , Classical Mechanics . Addison-Wesley , San Francisco Google Preview WorldCat COPAC Hakim R. J. , 2011 , Introduction to Relativistic Statistical Mechanics: Classical and Quantum . World Scientific , Singapore https://doi.org/10.1142/7881 Google Scholar Crossref Search ADS Google Preview WorldCat COPAC Jackson J. D. , 1998 , Classical Electrodynamics , 3rd edn. J. Wiley & Sons , New York Google Preview WorldCat COPAC Jokipii J. R. , Williams L. L. , 1992 , ApJ , 394 , 184 https://doi.org/10.1086/171569 Crossref Search ADS Jones F. C. , 1990 , ApJ , 361 , 162 https://doi.org/10.1086/169179 Crossref Search ADS Klimontovich Y. L. , 1960 , Sov. Phys. JETP , 37 , 524 https://doi.org/10.1086/169179 Kulsrud R. M. , 1983 , in Galeev A. A. , Sudan R. N. , eds, Basic Plasma Physics: Selected Chapters, Handbook of Plasma Physics , Volume 1 . North Holland Publ Co , Amsterdam , p. 1 Google Preview WorldCat COPAC Landau L. D. , Lifshitz E. M. , 1959 , Fluid Mechanics, Course of Theoretical Physics , Vol. 6 . Pergamon Press , Oxford WorldCat COPAC Landau L. D. , Lifshitz E. M. , 1975 , The Classical Theory of Fields, Course of Theoretical Physics , Vol. 2. : Pergamon Press , Oxford WorldCat COPAC Lichnerowicz A. , 1967 , Relativistic Hydrodynamics and Magnetohydrodynamics, The Mathematical Physics Monograph Series . Benjamin , New York Google Preview WorldCat COPAC Melrose D. B. , 1989 , Instabilities in Space and Laboratory Plasmas. Cambridge Univ. Press , Cambridge , WorldCat COPAC Mihalas D. , Weibel Mihalas B. , 1984 , Foundations of Radiation Hydrodynamics . Oxford University Press , New York Google Preview WorldCat COPAC Misner C. W. , Thorne K. S. , Wheeler J. A. , 1973 , Gravitation . W.H. Freeman and Co , San Francisco Google Preview WorldCat COPAC Miyamoto K. , 1980 , Plasma Physics for Nuclear Fusion . MIT Press , Cambridge, Mass . Google Preview WorldCat COPAC Parker E. N. , 1965 , Planet. Space Sci. , 13 , 9 https://doi.org/10.1016/0032-0633(65)90131-5 Crossref Search ADS Pomraning G. C. , 1973 , The Equations of Radiation Hydrodynamics, International Series of Monographs in Natural Philosophy . Pergamon Press , Oxford WorldCat COPAC Skilling J. , 1975 , MNRAS , 172 , 557 https://doi.org/10.1093/mnras/172.3.557 Crossref Search ADS Webb G. M. , 1985 , ApJ , 296 , 319 https://doi.org/10.1086/163451 Crossref Search ADS Weinberg S. , 1971 , ApJ , 168 , 175 https://doi.org/10.1086/151073 Crossref Search ADS Williams L. L. , Jokipii J. R. , 1993 , ApJ , 417 , 725 https://doi.org/10.1086/173351 Crossref Search ADS Zee A. , 2013 , Einstein Gravity in a Nutshell . Princeton Univ. Press , Princeton, NJ Google Preview WorldCat COPAC APPENDIX A: EQUATION OF MOTION IN MIXED VARIABLES The equation of motion using mixed variables (⁠|$x^{\mu }\ ,, p_{\hat{a}}$|⁠) can be derived using a variational principle. We use lab frame time and position and fluid rest frame momentum and energy, which are not canonically conjugate variables. In that situation one can employ the approach of Cary & Littlejohn (1983) for non-canonical Hamiltonian mechanics. To that end we quite generally define an action |${\cal A}$| for a freely falling particle in curved space–time with the metric gμν: \begin{eqnarray*} {\cal A} = \int_{\lambda_{1}}^{\lambda_{2}} \left[ \, p_{\mu } \, {\rm d} x^{\mu } - {\cal H}(p\,, x) \, {\rm d} \lambda \, \right] = \int_{\lambda_{1}}^{\lambda_{2}} \left[ \, p_{\mu } \, \frac{{\rm d} x^{\mu }}{{\rm d} \lambda } - {\cal H}(p\,, x) \, \right]\, {\rm d} \lambda . \end{eqnarray*} (A1) Here λ is an affine parameter that parametrizes the four-distance along the particle world line, and |$p_{\mu } \equiv g_{\mu \sigma }(x) \, p^{\sigma }$|⁠. For a free-falling particle the Hamiltonian |${\cal H}(p\ ,, x)$| is \begin{eqnarray*} {\cal H}(p\,, x) = \frac{1}{2m}\, p_{\mu } p^{\mu } = \frac{1}{2m} \, g^{\mu \nu }(x) \, p_{\mu } p_{\nu }. \end{eqnarray*} (A2) The quantity \begin{eqnarray*} {\cal L}(x\,, p) = \, p_{\mu } \, \frac{{\rm d} x^{\mu }}{{\rm d} \lambda } - {\cal H}(p\,, x) \end{eqnarray*} (A3) is known as the phase-space Lagrangian. The principle of stationary action is |$\delta {\cal A} = 0$|⁠, where the variations in pμ and xμ vanish at the end points of the integration over λ: δxμ(λ1, 2) = δpμ(λ1, 2) = 0. We have \begin{eqnarray*} \delta {\cal A} = \int_{\lambda_{1}}^{\lambda_{2}} \delta {\cal L}(x\,, p) \, {\rm d} \lambda = \int_{\lambda_{1}}^{\lambda_{2}} \left\lbrace \delta p_{\mu } \, \frac{{\rm d} x^{\mu }}{{\rm d} \lambda } + p_{\mu } \, \delta \left(\, \frac{{\rm d} x^{\mu }}{{\rm d} \lambda } \right) - \left(\frac{\partial {\cal H}}{\partial p_{\mu }} \, \delta p_{\mu } + \frac{\partial {\cal H}}{\partial x^{\mu }} \, \delta x^{\mu } \right) \right\rbrace {\rm d} \lambda . \end{eqnarray*} (A4) Since |$\delta ({\rm d} x^{\mu }/{\rm d} \lambda) = {\rm d} \, \delta x^{\mu }/{\rm d} \lambda$| we can partially integrate the second term: \begin{eqnarray*} \delta {\cal A} = \int_{\lambda_{1}}^{\lambda_{2}} \left\lbrace \delta p_{\mu } \left(\frac{{\rm d} x^{\mu }}{{\rm d} \lambda } - \frac{\partial {\cal H}}{\partial p_{\mu }} \right) - \delta x^{\mu } \, \left(\frac{{\rm d} p_{\mu }}{{\rm d} \lambda } + \frac{\partial {\cal H}}{\partial x^{\mu }} \right) \right\rbrace \, {\rm d} \lambda = 0 . \end{eqnarray*} (A5) With δpμ and δxμ arbitrary (except at the end points) we find the equations of motion (Euler–Lagrange equations): \begin{eqnarray*} \frac{{\rm d} x^{\mu }}{{\rm d} \lambda } = \frac{\partial {\cal H}}{\partial p_{\mu }} = \frac{p^{\mu }}{m}, \frac{{\rm d} p_{\mu }}{{\rm d} \lambda } = - \frac{\partial {\cal H}}{\partial x^{\mu }} = - \frac{1}{2m} \, \left(\frac{\partial g^{\alpha \beta }}{\partial x^{\mu }} \right) \, p_{\alpha } p_{\beta } . \end{eqnarray*} (A6) These are essentially the Hamiltonian equations and pμ and xμ are canonically conjugate variables. The first equation is simply dxμ/dλ = uμ = pμ/m, the correct definition for the four-velocity. The second equation can be cast in a somewhat more familiar form by employing the relation \begin{eqnarray*} \frac{\partial g^{\sigma \tau }}{\partial x^{\mu }} = - g^{\sigma \alpha } g^{\tau \beta } \, \frac{\partial g_{\alpha \beta }}{\partial x^{\mu }} . \end{eqnarray*} (A7) This follows from differentiation of the relation |$g^{\sigma \tau } = g^{\sigma \alpha } g^{\tau \beta } \, g_{\alpha \beta }$|⁠, the symmetry of the metric tensor, and relations like |$g_{\alpha \beta } g^{\tau \beta } = \delta_{\alpha }^{\tau }$|⁠. This allows us to write the second equation in equation (A6) in terms of uμ = pμ/m and the well-known connection coefficients |$\Gamma^{\alpha }_{\beta \gamma }$| (see for instance Carroll 2004, chapter 3): \begin{eqnarray*} \frac{{\rm d} u_{\mu }}{{\rm d} \lambda } = \frac{1}{2}\, \left(\frac{\partial \, g_{\alpha \beta }}{\partial x^{\mu }}\right) \, u^{\alpha } u^{\beta } = \Gamma^{\tau }_{\mu \beta } \, u^{\beta } u_{\tau } . \end{eqnarray*} (A8) This is the correct form for the parallel transport of the covariant momentum vector of a freely falling particle in curved space–time. We now change momentum variables from the LF four-momentum pμ to the FRF four-momentum |$p_{\hat{a}} = {\rm E}_{\, \hat{a}}^{\mu }(x) \, p_{\mu }$|⁠. We require that the new momentum variables are locally Lorentzian, i.e. the base vectors |${\boldsymbol {\hat{e}}}_{\hat{i}}$| employed in the FRF satisfy |${\boldsymbol {\hat{e}}}_{\hat{a}} {\boldsymbol{\cdot}}{\boldsymbol {\hat{e}}}_{\hat{b}} = \eta_{\hat{a}\hat{b}}$|⁠, where |${\boldsymbol {\hat{e}}}_{\hat{a}} = {\rm E}^{\mu }_{\, \hat{a}}(x) \, {\boldsymbol e }_{\mu }$| with the |${\boldsymbol e }_{\mu }$| the base vectors of a coordinated basis that satisfy |${\boldsymbol e }_{\mu } {\boldsymbol{\cdot}}{\boldsymbol e }_{\nu } = g_{\mu \nu }(x)$|⁠. Then |$p_{\mu } \, p^{\mu } = g_{\mu \nu } \, p^{\mu } p^{\nu } = p_{\hat{a}} \, p^{\hat{a}} = \eta_{\hat{a}\hat{b}} \, p^{\hat{a}} p^{\hat{b}}$|⁠. In addition, suppressing the space–time dependence in |${\rm E}_{\, \hat{a}}^{\mu }(x)$| and gμν, \begin{eqnarray*} g^{\mu \nu } \, {\rm E}_{\mu }^{\, \hat{a}} \, {\rm E}_{\nu }^{\, \hat{b}} = \eta^{\hat{a}\hat{b}}\,, \eta_{\hat{a}\hat{b}} \, {\rm E}^{\, \hat{a}}_{\mu } {\rm E}^{\, \hat{b}}_{\nu } = g_{\mu \nu }\,, {\rm E}_{\mu }^{\, \hat{a}} {\rm E}^{\nu }_{\, \hat{a}} = g^{\nu }_{\mu } = \delta^{\nu }_{\mu }\,, {\rm E}_{\mu }^{\, \hat{a}} {\rm E}^{\mu }_{\, \hat{b}} = \eta^{\hat{a}}_{\hat{b}} = \delta^{\hat{a}}_{\hat{b}} . \end{eqnarray*} (A9) These equations are the generalization to curved space–time of the relations in equation (4). Together they form the basis of the tetrad formalism employed by Webb (1985) in his derivation of the transport equation. A discussion of this formalism can be found in Carroll (2004, appendix J), whose notation we largely follow, and in Zee (2013, ch. IX.7), who employs the term Vielbein rather than tetrad. For our derivation of the equation of motion in mixed variables we will only need the relations given above. At the end of this section we show that the results are equivalent to those obtained using the full tetrad formalism. Using the FRF momentum variables the phase-space Lagrangian is \begin{eqnarray*} {\cal L}(x\ ,, p) = \, p_{\hat{a}} \, {\rm E}^{\, \hat{a}}_{\mu }(x) \, \frac{{\rm d} x^{\mu }}{{\rm d} \lambda } - \frac{p^{\hat{a}} p_{\hat{a}}}{2m} \, \equiv p_{\hat{a}} \, {\rm E}^{\, \hat{a}}_{\mu }(x) \, \frac{{\rm d} x^{\mu }}{{\rm d} \lambda } - \hat{{\cal H}}(p) . \end{eqnarray*} (A10) The ‘new’ Hamiltonian |$\hat{{\cal H}}(p) = p^{\hat{a}} p_{\hat{a}}/2m = \eta^{\hat{a}\hat{b}} \, p_{\hat{a}} p_{\hat{b}}/2m$| is now formally independent of space–time position. The principle of stationary action under the variation (⁠|$\delta x^{\mu }\ ,, \delta p_{\hat{a}}$|⁠) yields the following set of Euler–Lagrange equations: \begin{eqnarray*} {\rm E}_{\mu }^{\, \hat{a}} \, \frac{{\rm d} x^{\mu }}{{\rm d} \lambda } = \frac{\partial \hat{{\cal H}}}{\partial p_{\hat{a}}} = \frac{p^{\hat{a}}}{m}, \frac{{\rm d}}{{\rm d} \lambda } \left({\rm E}_{\mu }^{\, \hat{a}} \, p_{\hat{a}} \right) - p_{\hat{a}} \, \left(\frac{\partial {\rm E}_{\nu }^{\, \hat{a}}}{\partial x^{\mu }} \right) \, \frac{{\rm d} x^{\nu }}{{\rm d} \lambda } = - \frac{\partial \hat{\cal H}}{\partial x^{\mu }} = 0 . \end{eqnarray*} (A11) Here we used |$p^{\hat{a}} = \eta^{\hat{a}\hat{b}} \, p_{\hat{b}}$| and |$\delta {\rm E}^{\, \hat{a}}_{\nu } = (\partial {\rm E}^{\, \hat{a}}_{\nu }/\partial x^{\mu }) \, \delta x^{\mu }$|⁠. With equation (A9) the Euler–Lagrange equations can be written as \begin{eqnarray*} \frac{{\rm d} x^{\nu }}{{\rm d} \lambda } = {\rm E}^{\nu }_{\, \hat{a}} \, \frac{p^{\hat{a}}}{m} = u^{\nu }\,, \frac{{\rm d} \, p_{\hat{c}}}{{\rm d} \lambda } = \omega^{\hat{a}}_{\hat{b}\hat{c}}\, \frac{p_{\hat{a}} p^{\hat{b}}}{m} . \end{eqnarray*} (A12) Defining |$\partial_{\hat{a}} = {\rm E}^{\mu }_{\, \hat{a}} (\partial /\partial x^{\mu })$| the |$\omega^{\hat{a}}_{\hat{b}\hat{c}}$| are given by \begin{eqnarray*} \omega^{\hat{a}}_{\hat{b}\hat{c}}\equiv {\rm E}^{\nu }_{\, \hat{b}} \, {\rm E}^{\mu }_{\, \hat{c}} \, \left(\frac{\partial {\rm E}_{\nu }^{\, \hat{a}}}{\partial x^{\mu }} - \frac{\partial {\rm E}_{\mu }^{\, \hat{a}} \,}{\partial x^{\nu }} \right) = {\rm E}^{\nu }_{\, \hat{b}} \, \partial_{\hat{c}} {\rm E}_{\nu }^{\, \hat{a}} - {\rm E}^{\nu }_{\, \hat{c}} \, \partial_{\hat{b}} {\rm E}_{\nu }^{\, \hat{a}} . \end{eqnarray*} (A13) These equations do not take the canonical form: xμ and |$p_{\hat{a}}$| are not canonically conjugate variables. Result (A12) is in agreement with the results of the full tetrad formalism, as we now show briefly. The |$\omega^{\hat{a}}_{\hat{b}\hat{c}}$| are the commutation coefficients of the orthonormal tetrad base vectors: |${\boldsymbol {\hat{e}}}_{\hat{i}}$|⁠: |$[ \, {\boldsymbol {\hat{e}}}_{\hat{a}}\,, {\boldsymbol {\hat{e}}}_{\hat{b}} \, ] = \omega^{\hat{c}}_{\hat{a}\hat{b}} \, {\boldsymbol {\hat{e}}}_{\hat{c}}$|⁠. The second equation in equation (A12) is equivalent with the standard form in terms of |$u_{\hat{a}} = p_{\hat{a}}/m$|⁠: \begin{eqnarray*} \frac{{\rm d} u_{\hat{a}}}{{\rm d} \lambda } + \Gamma_{\hat{a}\hat{c}\hat{b}} \, u^{\hat{b}} u^{\hat{c}} = 0 . \end{eqnarray*} (A14) Here |$\Gamma_{\hat{a}\hat{c}\hat{b}} \equiv \eta_{\hat{a}\hat{d}} \, \Gamma^{\hat{d}}_{\hat{c}\hat{b}}$| and the tetrad connection coefficients |$\Gamma^{\hat{a}}_{\hat{c}\hat{b}}$| are defined in terms of the covariant derivative of the |${\rm E}^{\mu }_{\, \hat{a}}$| as \begin{eqnarray*} \Gamma^{\hat{a}}_{\hat{c}\hat{b}} \equiv {\rm E}_{\nu }^{\, \hat{a}} \, {\rm E}^{\mu }_{\, \hat{b}} \, \left(\frac{\partial {\rm E}^{\nu }_{\, \hat{c}}}{\partial x^{\mu }} + \Gamma^{\nu }_{\mu \lambda } \, {\rm E}^{\lambda }_{\, \hat{c}} \right) = - \, {\rm E}^{\mu }_{\, \hat{b}} \, {\rm E}^{\nu }_{\, \hat{c}} \, \left(\frac{\partial {\rm E}_{\nu }^{\, \hat{a}}}{\partial x^{\mu }} - \Gamma^{\lambda }_{\mu \nu } {\rm E}_{\lambda }^{\, \hat{a}} \right) . \end{eqnarray*} (A15) Defining |$\omega_{\hat{a}\hat{b}\hat{c}} = \eta_{\hat{a}\hat{d}} \, \omega^{\hat{d}}_{\hat{b}\hat{c}}$| and using that in general relativity space is torsion free with |$\Gamma^{\lambda }_{\mu \nu } = \Gamma^{\lambda }_{\nu \mu }$|⁠, so that in the definition (A13) for the |$\omega^{\hat{a}}_{\hat{b}\hat{c}}$| ordinary derivatives can be replaced by covariant derivatives, and one can show (e.g. Misner et al. 1973, p. 358) \begin{eqnarray*} \Gamma_{\hat{a}\hat{b}\hat{c}} = - \Gamma_{\hat{b}\hat{a}\hat{c}} = \frac{1}{2}\left(\omega_{\hat{b}\hat{a}\hat{c}} + \omega_{\hat{c}\hat{a}\hat{b}} - \omega_{\hat{a}\hat{b}\hat{c}} \right)\,, \omega_{\hat{a}\hat{b}\hat{c}} = - \omega_{\hat{a}\hat{c}\hat{b}} . \end{eqnarray*} (A16) With these relations one can prove |$\Gamma_{\hat{a}\hat{c}\hat{b}} \, u^{\hat{b}} u^{\hat{c}} = \frac{1}{2}(\Gamma_{\hat{a}\hat{b}\hat{c}} + \Gamma_{\hat{a}\hat{c}\hat{b}}) \, u^{\hat{b}} u^{\hat{c}} = \omega_{\hat{b}\hat{a}\hat{c}} \, u^{\hat{b}} u^{\hat{c}} = - \omega_{\hat{b}\hat{c}\hat{a}} u^{\hat{b}} u^{\hat{c}} = - \omega^{\hat{b}}_{\hat{c}\hat{a}} \, u_{\hat{b}} u^{\hat{c}}$|⁠, and the second relation in equation (A12) follows after a relabelling of indices and writing everything in terms of momenta rather than four-velocities. Finally, in flat space–time the term involving |${\rm E}^{\nu }_{\, \hat{b}} \, (\partial {\rm E}_{\nu }^{\, \hat{a}}/\partial x^{\mu })$| in the equation for |$p_{\hat{c}}$| does not contribute. From the above definitions it follows generally that \begin{eqnarray*} {\rm E}^{\nu }_{\, \hat{b}}\, \left(\frac{\partial {\rm E}_{\nu }^{\, \hat{a}}}{\partial x^{\mu }} \right) \, \frac{p_{\hat{a}} \, p^{\hat{b}}}{m} = \frac{1}{2m} \left(\frac{\partial \, g_{\alpha \beta }}{\partial x^{\mu }} \right) \, p^{\alpha } p^{\beta } \,, \end{eqnarray*} (A17) which vanishes in flat space–time where gαβ = ηαβ. In that case the second equation in equation (A12) reduces to equation (9) of the main paper. APPENDIX B: CALCULATION OF THE FICTITIOUS FORCE TERMS B1 Energy change The change in the FRF energy |$\bar{{\cal E}}$| follows from using |$p_{\hat{0}} = \bar{{\cal E}}$| and |${\rm E}_{\hat{0}}^{\mu } = U^{\mu }$| in equation (9) with |$\hat{c}= 0$|⁠: \begin{eqnarray*} \frac{{\rm d} p_{\hat{0}}}{{\rm d} \lambda } = \frac{{\rm d} \bar{{\cal E}}}{{\rm d} \lambda } = - {\rm E}^{\mu }_{\hat{0}} \, \frac{{\rm d} {\rm E}_{\mu }^{\hat{a}}}{{\rm d} \lambda } \, p_{\hat{a}} = - U^{\mu } \, \left(\frac{{\rm d} {\boldsymbol{U}}}{{\rm d} \lambda } {\boldsymbol{\cdot}}\frac{\partial E^{\hat{a}}_{\mu }}{\partial {\boldsymbol{U}}} \right) \, p_{\hat{a}} . \end{eqnarray*} (B1) Here we employ the three-vector notation (with implied summation over double roman indices) \begin{eqnarray*} \frac{{\rm d} {\boldsymbol{U}}}{{\rm d} \lambda } {\boldsymbol{\cdot}}\frac{\partial {\rm E}^{\hat{a}}_{\mu }}{\partial {\boldsymbol{U}}} \equiv \frac{{\rm d} U^{k}}{{\rm d} \lambda } \, \frac{\partial {\rm E}_{\mu }^{\hat{a}}}{\partial U^{k}} . \end{eqnarray*} (B2) Since |$p_{\hat{a}} = (\bar{{\cal E}}, - \bar{{\boldsymbol{p}}})$| one has \begin{eqnarray*} p_{\hat{a}} \, \frac{\partial {\rm E}_{\mu }^{\hat{a}}}{\partial {\boldsymbol{U}}} = \bar{{\cal E}}\, \frac{\partial {\rm E}_{\mu }^{\hat{0}}}{\partial {\boldsymbol{U}}} - \bar{p}^{i} \, \frac{\partial {\rm E}_{\mu }^{\hat{i}}}{\partial {\boldsymbol{U}}} \end{eqnarray*} (B3) and \begin{eqnarray*} \frac{{\rm d} \bar{{\cal E}}}{{\rm d} \lambda } = - U^{\mu } \, \left[ \, \bar{{\cal E}}\, \left(\frac{{\rm d} {\boldsymbol{U}}}{{\rm d} \lambda } {\boldsymbol{\cdot}}\frac{\partial {\rm E}_{\mu }^{\hat{0}}}{\partial {\boldsymbol{U}}} \right) - \bar{p}^{i} \, \left(\frac{{\rm d} {\boldsymbol{U}}}{{\rm d} \lambda } {\boldsymbol{\cdot}}\frac{\partial {\rm E}_{\mu }^{\hat{i}}}{\partial {\boldsymbol{U}}} \right) \, \right]. \end{eqnarray*} (B4) It is convenient to look at the contribution from μ = 0 and from |$\mu = 1\,, 2\,, 3$| separately. Using the relations (13) from the main paper in the form \begin{eqnarray*} \frac{\partial {\rm E}_{0}^{\hat{0}}}{\partial {\boldsymbol{U}}} = {\boldsymbol{U}}/\Gamma\,, \bar{p}^{i} \, \frac{\partial {\rm E}_{0}^{\hat{i}}}{\partial {\boldsymbol{U}}} = - \bar{{\boldsymbol{p}}}\,, \end{eqnarray*} (B5) one sees that the μ = 0 term gives a contribution \begin{eqnarray*} \left(\frac{{\rm d} \bar{{\cal E}}}{{\rm d} \lambda } \right)_{\mu = 0} = - \left(\bar{{\cal E}}{\boldsymbol{U}} + \Gamma \, \bar{{\boldsymbol{p}}}\right) {\boldsymbol{\cdot}}\frac{{\rm d} {\boldsymbol{U}}}{{\rm d} \lambda } . \end{eqnarray*} (B6) The terms with |$\mu = 1\,, 2\,, 3$| follow from the second and third relation in equation (13) and add up to \begin{eqnarray*} \left(\frac{{\rm d} \bar{{\cal E}}}{{\rm d} \lambda } \right)_{\mu = 1,2,3} =- U^{j} \left[ \, \bar{{\cal E}}\, (- \delta_{jk}) - \frac{\bar{p}^{i}}{\Gamma + 1} \left(\delta_{jk} U^{i} + \delta_{ik} U^{j} - \frac{U^{i} U^{j} U^{k}}{\Gamma (\Gamma + 1)}\right) \, \right]\, \frac{{\rm d} U^{k}}{{\rm d} \lambda } \,. \end{eqnarray*} (B7) In three-vector notation \begin{eqnarray*} \left(\frac{{\rm d} \bar{{\cal E}}}{{\rm d} \lambda } \right)_{\mu = 1,2,3} = \left[ \, \left(\bar{{\cal E}}+ \frac{\bar{{\boldsymbol{p}}}{\boldsymbol{\cdot}}{\boldsymbol{U}}}{\Gamma + 1} - \frac{U^2}{\Gamma (\Gamma + 1)^2} \, (\bar{{\boldsymbol{p}}}{\boldsymbol{\cdot}}{\boldsymbol{U}}) \right) \, {\boldsymbol{U}} + \left(\frac{U^2}{\Gamma + 1} \right) \, \bar{{\boldsymbol{p}}}\, \right] {\boldsymbol{\cdot}}\frac{{\rm d} {\boldsymbol{U}}}{{\rm d} \lambda } . \end{eqnarray*} (B8) Adding contributions (B6) and (B8) one finds that the terms |$\propto \bar{{\cal E}}$| cancel. Then by using U2 = Γ2 − 1 the remaining terms can be written as \begin{eqnarray*} \frac{{\rm d} \bar{{\cal E}}}{{\rm d} \lambda } = - \left(\bar{{\boldsymbol{p}}}- \frac{(\bar{{\boldsymbol{p}}}{\boldsymbol{\cdot}}{\boldsymbol{U}})}{\Gamma (\Gamma + 1)} \, {\boldsymbol{U}} \right) {\boldsymbol{\cdot}}\frac{{\rm d} {\boldsymbol{U}}}{{\rm d} \lambda } = - (\Gamma + 1) \, \left[ \bar{{\boldsymbol{p}}}{\boldsymbol{\cdot}}\frac{{\rm d}}{{\rm d} \lambda } \left(\frac{{\boldsymbol{U}}}{\Gamma + 1} \right) \right]. \end{eqnarray*} (B9) The second equality follows from |${\boldsymbol{U}} {\boldsymbol{\cdot}}({\rm d} {\boldsymbol{U}}/{\rm d} \lambda) = \Gamma ({\rm d} \Gamma /{\rm d}\lambda)$|⁠. B2 Momentum change The change in FRF three-momentum |$\bar{{\boldsymbol{p}}}$| is calculated from equation (9) using |$p_{\hat{i}} = - p^{\hat{i}}$| so that \begin{eqnarray*} \frac{{\rm d} \bar{p}^{\hat{i}}}{{\rm d} \lambda } = - \frac{{\rm d} \bar{p}_{\hat{i}}}{{\rm d} \lambda } = {\rm E}^{\mu }_{\, \hat{i}} \, \frac{{\rm d} {\rm E}^{\hat{a}}_{\mu }}{{\rm d} \lambda } \, p_{\hat{a}} = {\rm E}^{\mu }_{\, \hat{i}}\, \left(\bar{{\cal E}}\, \frac{{\rm d} {\rm E}^{\hat{0}}_{\mu }}{{\rm d} \lambda } - \bar{p}^{j} \, \frac{{\rm d} {\rm E}^{\hat{j}}_{\mu }}{{\rm d} \lambda } \right) = {\rm E}^{\mu }_{\, \hat{i}}\, \left(\bar{{\cal E}}\, \frac{\partial {\rm E}^{\hat{0}}_{\mu }}{\partial {\boldsymbol{U}}} - \bar{p}^{j} \, \frac{\partial {\rm E}^{\hat{j}}_{\mu }}{\partial {\boldsymbol{U}}} \right) {\boldsymbol{\cdot}}\frac{{\rm d} {\boldsymbol{U}}}{{\rm d} \lambda } . \end{eqnarray*} (B10) Again using the first two relations in equation (13) (or equivalently equation B2) it follows that the contribution of the μ = 0 term is (in three-vector notation) \begin{eqnarray*} \left(\frac{{\rm d} \bar{{\boldsymbol{p}}}}{{\rm d} \lambda } \right)_{\mu = 0} = \left\lbrace \left(\bar{{\cal E}}\, {\boldsymbol{V}} + \bar{{\boldsymbol{p}}}\right) {\boldsymbol{\cdot}}\frac{{\rm d} {\boldsymbol{U}}}{{\rm d} \lambda } \right\rbrace \, {\boldsymbol{U}} . \end{eqnarray*} (B11) The second and third relation in equation (13) yields for the terms from |$\mu = 1\,, 2\,, 3$| \begin{eqnarray*} \left(\frac{{\rm d} \bar{p}^i}{{\rm d} \lambda } \right)_{\mu = 1,2,3} = \left(\delta_{si} + \frac{U^{s} U^{i}}{\Gamma + 1} \right) \left\lbrace \bar{{\cal E}}\, \left(- \delta_{sk} \right) - \frac{\bar{p}^{j}}{\Gamma + 1} \, \left(\delta_{sk} U^{j} + \delta_{jk} U^{s} - \frac{U^{s} U^{j} U^{k}}{\Gamma (\Gamma + 1)} \right) \right\rbrace \, \frac{{\rm d} U^{k}}{{\rm d} \lambda } . \end{eqnarray*} (B12) The terms involving the Kronecker delta function δsi can be written in three-vector notation as \begin{eqnarray*} \left(\frac{{\rm d} \bar{{\boldsymbol{p}}}}{{\rm d} \lambda } \right)_{\mu = 1,2,3}^{(1)} = - \left(\bar{{\cal E}}+ \frac{\bar{{\boldsymbol{p}}}{\boldsymbol{\cdot}}{\boldsymbol{U}}}{\Gamma + 1} \right) \, \frac{{\rm d} {\boldsymbol{U}}}{{\rm d} \lambda } - \frac{{\boldsymbol{U}}}{\Gamma + 1} \, \left\lbrace \left(\bar{{\boldsymbol{p}}}{\boldsymbol{\cdot}}\frac{{\rm d} {\boldsymbol{U}}}{{\rm d} \lambda } \right) - \frac{\bar{{\boldsymbol{p}}}{\boldsymbol{\cdot}}{\boldsymbol{U}}}{\Gamma (\Gamma +1)} \, \left({\boldsymbol{U}} {\boldsymbol{\cdot}}\frac{{\rm d} {\boldsymbol{U}}}{{\rm d} \lambda } \right) \right\rbrace . \end{eqnarray*} (B13) The terms involving the dyadic product UsUi yield: \begin{eqnarray*} \left(\frac{{\rm d} \bar{p}^i}{{\rm d} \lambda } \right)_{\mu = 1,2,3}^{(2)} = - \frac{U^i}{\Gamma + 1} \, \left(\bar{{\cal E}}+ \frac{\bar{{\boldsymbol{p}}}{\boldsymbol{\cdot}}{\boldsymbol{U}}}{\Gamma + 1} \right) \, \left({\boldsymbol{U}} {\boldsymbol{\cdot}}\frac{{\rm d} {\boldsymbol{U}}}{{\rm d} \lambda } \right) - \frac{U^2 \, U^{i}}{(\Gamma + 1)^2} \, \left\lbrace \left(\bar{{\boldsymbol{p}}}{\boldsymbol{\cdot}}\frac{{\rm d} {\boldsymbol{U}}}{{\rm d} \lambda } \right) - \frac{\bar{{\boldsymbol{p}}}{\boldsymbol{\cdot}}{\boldsymbol{U}}}{\Gamma (\Gamma +1)} \, \left({\boldsymbol{U}} {\boldsymbol{\cdot}}\frac{{\rm d} {\boldsymbol{U}}}{{\rm d} \lambda } \right) \right\rbrace . \end{eqnarray*} (B14) Using U2 = Γ2 − 1 = (Γ − 1)(Γ + 1) and changing to three-vector notation this becomes: \begin{eqnarray*} \left(\frac{{\rm d} \bar{{\boldsymbol{p}}}}{{\rm d} \lambda } \right)_{\mu = 1,2,3}^{(2)} = - \frac{{\boldsymbol{U}}}{\Gamma + 1} \, \left\lbrace \, \left(\bar{{\cal E}}+ \frac{\bar{{\boldsymbol{p}}}{\boldsymbol{\cdot}}{\boldsymbol{U}}}{\Gamma (\Gamma + 1)} \right) \, \left({\boldsymbol{U}} {\boldsymbol{\cdot}}\frac{{\rm d} {\boldsymbol{U}}}{{\rm d} \lambda } \right) + (\Gamma - 1) \, \left(\bar{{\boldsymbol{p}}}{\boldsymbol{\cdot}}\frac{{\rm d} {\boldsymbol{U}}}{{\rm d} \lambda } \right) \right\rbrace . \end{eqnarray*} (B15) Adding the two terms (B13) and (B15) leads to a number of cancellations and one finds for the term from |$\mu = 1\ ,, 2\ ,, 3$|⁠: \begin{eqnarray*} \left(\frac{{\rm d} \bar{{\boldsymbol{p}}}}{{\rm d} \lambda } \right)_{\mu = 1,2,3} = - \left(\bar{{\cal E}}+ \frac{\bar{{\boldsymbol{p}}}{\boldsymbol{\cdot}}{\boldsymbol{U}}}{\Gamma + 1} \right) \, \frac{{\rm d} {\boldsymbol{U}}}{{\rm d} \lambda } - \frac{{\boldsymbol{U}}}{\Gamma +1} \, \left\lbrace \Gamma \, \left(\bar{{\boldsymbol{p}}}{\boldsymbol{\cdot}}\frac{{\rm d} {\boldsymbol{U}}}{{\rm d} \lambda } \right) + \bar{{\cal E}}\, \left({\boldsymbol{U}} {\boldsymbol{\cdot}}\frac{{\rm d} {\boldsymbol{U}}}{{\rm d} \lambda } \right) \right\rbrace . \end{eqnarray*} (B16) Finally, adding equation (B11) and (B16) and using |${\boldsymbol{U}} = \Gamma \, {\boldsymbol{V}}$| one finds the total FRF three-momentum change due to the fictitious force: \begin{eqnarray*} \frac{{\rm d} \bar{{\boldsymbol{p}}}}{{\rm d} \lambda } = - \left(\bar{{\cal E}}+ \frac{\bar{{\boldsymbol{p}}}{\boldsymbol{\cdot}}{\boldsymbol{U}}}{\Gamma + 1} \right) \, \frac{{\rm d} {\boldsymbol{U}}}{{\rm d} \lambda } + \left((\bar{{\boldsymbol{p}}}+ \bar{{\cal E}}\, {\boldsymbol{V}}) {\boldsymbol{\cdot}}\frac{{\rm d} {\boldsymbol{U}}}{{\rm d} \lambda } \right) \frac{{\boldsymbol{U}}}{\Gamma + 1} . \end{eqnarray*} (B17) The relation |${\boldsymbol{U}} {\boldsymbol{\cdot}}({\rm d} {\boldsymbol{U}}/{\rm d} \lambda) = \Gamma ({\rm d} \Gamma /{\rm d}\lambda)$| allows one to cast this into the form \begin{eqnarray*} \frac{{\rm d} \bar{{\boldsymbol{p}}}}{{\rm d} \lambda } = - (\Gamma + 1) \, \bar{{\cal E}}\, \frac{{\rm d}}{{\rm d} \lambda } \left(\frac{{\boldsymbol{U}}}{\Gamma + 1} \right) + \frac{1}{\Gamma +1} \, \left[ \, \left(\bar{{\boldsymbol{p}}}{\boldsymbol{\cdot}}\frac{{\rm d} {\boldsymbol{U}}}{{\rm d} \lambda } \right) \, {\boldsymbol{U}} - \left(\bar{{\boldsymbol{p}}}{\boldsymbol{\cdot}}{\boldsymbol{U}} \right) \, \frac{{\rm d} {\boldsymbol{U}}}{{\rm d} \lambda } \, \right]. \end{eqnarray*} (B18) For our purposes this is most conveniently written as \begin{eqnarray*} \frac{{\rm d} \bar{{\boldsymbol{p}}}}{{\rm d} \lambda } = - (\Gamma + 1) \, \bar{{\cal E}}\, \frac{{\rm d}}{{\rm d} \lambda } \left(\frac{{\boldsymbol{U}}}{\Gamma + 1} \right) + \bar{{\boldsymbol{p}}}{\boldsymbol{\times}}\left(\frac{{\boldsymbol{U}} {\boldsymbol{\times}}\, {\displaystyle \frac{{\rm d} {\boldsymbol{U}}}{{\rm d} \lambda }}}{\Gamma +1} \right) \equiv \bar{\gamma } \, \left({\boldsymbol K }_{\rm f} + \bar{{\boldsymbol{p}}}{\boldsymbol{\times}}{\boldsymbol {\bar{\omega}}} \right) . \end{eqnarray*} (B19) In this expression we introduce |${\boldsymbol K }_{\rm f}$| and |${\boldsymbol {\bar{\omega}}}$| that are defined as \begin{eqnarray*} {\boldsymbol K }_{\rm f} \equiv - (\Gamma + 1) \, \frac{{\rm d}}{{\rm d} \lambda } \left(\frac{m {\boldsymbol{U}}}{\Gamma + 1} \right)\,, {\boldsymbol {\bar{\omega}}} \equiv \frac{{\boldsymbol{U}} {\boldsymbol{\times}}\, {\displaystyle \frac{{\rm d} {\boldsymbol{U}}}{{\rm d} \lambda }}}{\bar{\gamma }(\Gamma +1)} . \end{eqnarray*} (B20) Relations (B9) and (B19/B20) are used in the main paper. It is easily checked using (B9) and (B19) that the fictitious four-force respects the Lorentz-invariant |$p^{\mu } p_{\mu } = p^{\hat{a}}p_{\hat{a}} = \bar{{\cal E}}^2 - | \bar{{\boldsymbol{p}}}|^2 = m^2$|⁠, as |$\bar{{\cal E}}({\rm d} \bar{{\cal E}}/{\rm d} \lambda) - \bar{{\boldsymbol{p}}}{\boldsymbol{\cdot}}({\rm d} \bar{{\boldsymbol{p}}}/ {\rm d} \lambda) = 0$|⁠. This serves as an independent check on the calculation presented here. APPENDIX C: RELATIVISTIC MHD We briefly derive a number of relations from relativistic MHD that are needed in the main paper. Additional details can be found in the book by Lichnerowicz (1967), the seminal paper of Bekenstein & Oron (1978), and in Achterberg (1983). Quite generally, one can define electric and magnetic four-vectors using the four-velocity Uμ of the fluid, Faraday tensor Fμν and its dual *Fμν: \begin{eqnarray*} {\cal E}^{\mu } \equiv F^{\mu \nu } \, U_{\nu }\,, {\cal B}^{\mu } ={}^{\ast } F^{\mu \nu } \, U_{\nu } . \end{eqnarray*} (C1) Here Fμν is the Faraday tensor and |${}^{\ast } F^{\mu \nu } = \frac{1}{2}\epsilon^{\mu \nu \sigma \tau } \, F_{\sigma \tau }$| is its Hodge dual. In the FRF one has |${\cal E}^{\mu } = (0\,, {\boldsymbol {\bar{E}} })$| and |${\cal B}^{\mu } = (0\,, {\boldsymbol {\bar{B}} })$| with |$({\boldsymbol {\bar{E}} }\,, {\boldsymbol {\bar{B}} })$| the electromagnetic fields in that frame. The antisymmetry of Fμν and its Hodge dual *Fμν implies that |${\cal E}^{\mu } U_{\mu } = {\cal B}^{\mu } U_{\mu } = 0$|⁠. Ideal MHD holds for the plasma at large when \begin{eqnarray*} {\cal E}^{\mu } = 0 \,, \end{eqnarray*} (C2) which is equivalent with the well-known condition |${\boldsymbol E } + {\boldsymbol{V}} {\boldsymbol{\times}}{\boldsymbol B } = 0$|⁠, valid in any frame. In that case the Faraday tensor and its dual are given by \begin{eqnarray*} F^{\mu \nu } = \epsilon^{\mu \nu \sigma \tau } U_{\sigma } B_{\tau } , {}^{\ast }F^{\mu \nu } = U^{\nu } {\cal B}^{\mu } - U^{\mu } {\cal B}^{\nu } . \end{eqnarray*} (C3) The induction equation is |$\nabla_{\mu } \left({}^{\ast }F^{\mu \nu } \right) = \nabla_{\mu } \left(U^{\nu } {\cal B}^{\mu } - U^{\mu } {\cal B}^{\nu } \right) = 0$|⁠. The resulting set of equations for the magnetic field can be written as \begin{eqnarray*} (U {\boldsymbol{\cdot}}\nabla)\left(\frac{{\cal B}^{\mu }}{n} \right) \equiv \frac{{\rm d}}{{\rm d} \tau } \left(\frac{{\cal B}^{\mu }}{n} \right) = \left(\frac{{\cal B}}{n} {\boldsymbol{\cdot}}\nabla \right) U^{\mu } - \left(a {\boldsymbol{\cdot}}\frac{{\cal B}}{n} \right) U^{\mu } \,, \, \nabla {\boldsymbol{\cdot}}{\cal B}= - a {\boldsymbol{\cdot}}{\cal B}. \end{eqnarray*} (C4) Here n is the proper density of the fluid that satisfies the particle conservation law \begin{eqnarray*} {\rm d}n/{\rm d} \tau = - n \, (\nabla {\boldsymbol{\cdot}}U) \,, \end{eqnarray*} (C5) and aμ ≡ dUμ/dτ is the fluid four-acceleration that satisfies |$U {\boldsymbol{\cdot}}a = 0$|⁠. We employ the notation |$A {\boldsymbol{\cdot}}B \equiv A_{\mu } B^{\mu }$| throughout. The perhaps unfamiliar terms in these two equations (in a direct comparison with the corresponding equations in non-relativistic MHD) involve the four-acceleration aμ. They are needed to maintain the orthogonality between |${\cal B}^{\mu }$| and Uμ, (i.e. |${\cal B}{\boldsymbol{\cdot}}U = 0$|⁠). These terms can be neglected in a non-relativistic flow. Contracting the first equation in equation (C4) with |${\cal B}_{\mu }/n$| it follows with the help of equation (C5) that \begin{eqnarray*} \frac{1}{{{\rm B}}} \, \frac{{\rm d} {{\rm B}}}{{\rm d} \tau } = - \left(b \, b {{\boldsymbol \, : \,}}\nabla U + \nabla {\boldsymbol{\cdot}}U \right) . \end{eqnarray*} (C6) Here |${{\rm B}}\equiv \sqrt{- {\cal B}^{\mu } {\cal B}_{\mu }}$| is the magnetic field strength in the FRF, where |${\cal B}^{\mu } = (0\ ,, {\boldsymbol {\bar{B}} })$| so that |${\cal B}^{\mu } {\cal B}_{\mu } = - | {\boldsymbol {\bar{B}} } |^2$|⁠. The vector bμ is defined as |$b^{\mu } \equiv {\cal B}^{\mu }/{{\rm B}}$| so that |$b^{\mu } \, b_{\mu } = -1$|⁠. It takes the form |$b^{\mu } = (0\,, \mbox{$\boldsymbol{\hat{b}}$})$| in the FRF. We employ the notation |$b \, b {{\boldsymbol \, : \,}}\nabla U \equiv b_{\mu } b^{\nu } \, \nabla_{\nu } U^{\mu }$|⁠. The second relation in equation (C4) leads, after substituting |${\cal B}^{\mu } = {{\rm B}}\, b^{\mu }$| to \begin{eqnarray*} \nabla {\boldsymbol{\cdot}}b = - \frac{(b {\boldsymbol{\cdot}}\nabla) {{\rm B}}}{{{\rm B}}} - a {\boldsymbol{\cdot}}b = - \frac{(b {\boldsymbol{\cdot}}\nabla) {{\rm B}}}{{{\rm B}}} - b \, U {{\boldsymbol \, : \,}}\nabla U . \end{eqnarray*} (C7) We will convert the above results to a mixed form employing FRF magnetic field variables B and |$\mbox{$\boldsymbol{\hat{b}}$}$| and LF time t and position |${\boldsymbol x }$|⁠. In the FRF one has locally|$U^{\hat{a}} = (1\,, 0\,, 0 \,, 0)$| and the covariant gradient operator is (in column notation) \begin{eqnarray*} \nabla_{\hat{a}} = {\rm E}^{\mu }_{\hat{a}} \, \nabla_{\mu } = \left(\begin{array}{c}\displaystyle \Gamma \, \frac{\partial }{\partial t} + ({\boldsymbol{U}} {\boldsymbol{\cdot}}{\nabla}) \\ \\ {\displaystyle {\nabla}+ {\boldsymbol{U}} \, \left(\frac{\partial }{\partial t} + \frac{{\boldsymbol{U}} {\boldsymbol{\cdot}}{\nabla}}{\Gamma + 1} \right)} \\ \end{array} \right) = \left(\begin{array}{c}\displaystyle \frac{{\rm d}}{{\rm d} \tau } \\ \\ \\ {\boldsymbol {\tilde{{\nabla}}}}\\ \end{array} \right) . \end{eqnarray*} (C8) The velocity gradient tensor evaluated in the FRF is \begin{eqnarray*} \nabla_{\hat{a}} U^{\hat{b}} = {\rm E}^{\mu }_{\hat{a}} \, {\rm E}_{\nu }^{\, \hat{b}} \, \left(\nabla_{\mu } U^{\nu } \right) = {\rm E}_{\nu }^{\, \hat{b}} \, \nabla_{\hat{a}} U^{\nu } = \left(\begin{array}{cc}0 & {\displaystyle (\Gamma + 1) \, \frac{{\rm d}}{{\rm d} \tau } \left(\frac{{\boldsymbol{U}}}{\Gamma +1} \right)} \\ & \\ {\boldsymbol 0 } & {\displaystyle (\Gamma + 1) \, {\boldsymbol {\tilde{{\nabla}}}}\left(\frac{{\boldsymbol{U}} }{\Gamma +1} \right)} \\ \end{array} \right) \,, \end{eqnarray*} (C9) using a 4 × 4 matrix notation in the last equality. This implies for the four-divergence of Uμ: \begin{eqnarray*} \nabla {\boldsymbol{\cdot}}U = \frac{\partial \Gamma }{\partial t} + {\nabla}{\boldsymbol{\cdot}}{\boldsymbol{U}} = \nabla_{\hat{a}}U^{\hat{a}} = (\Gamma + 1) \, {\boldsymbol {\tilde{{\nabla}}}}{\boldsymbol{\cdot}}\left(\frac{{\boldsymbol{U}} }{\Gamma +1} \right) \ ,, \end{eqnarray*} (C10) an identity that is easily checked using the definition of |${\boldsymbol {\tilde{{\nabla}}}}$| (equation 34). The FRF four-acceleration is |$a^{\hat{a}} = E^{\hat{a}}_{\, \mu } \, ({\rm d} U^{\mu }/{\rm d} \tau) = \left(0\,, {\boldsymbol {\bar{a}} } \right)$|⁠, with \begin{eqnarray*} {\boldsymbol {\bar{a}} } = (\Gamma + 1) \, \frac{{\rm d}}{{\rm d}\tau } \left(\frac{{\boldsymbol{U}}}{\Gamma +1} \right) . \end{eqnarray*} (C11) Using these relations equation (C6) becomes \begin{eqnarray*} \frac{1}{{{\rm { B}}}} \, \frac{{\rm d} {{\rm { B}}}}{{\rm d} \tau } = (\Gamma +1) \, \left\lbrace \mbox{$\boldsymbol{\hat{b}}$}\, \mbox{$\boldsymbol{\hat{b}}$}{{\boldsymbol \, : \,}}{\boldsymbol {\tilde{{\nabla}}}}\left(\frac{{\boldsymbol{U}} }{\Gamma +1} \right) - {\boldsymbol {\tilde{{\nabla}}}}{\boldsymbol{\cdot}}\left(\frac{{\boldsymbol{U}} }{\Gamma +1} \right) \right\rbrace \,, \end{eqnarray*} (C12) and relation (C7) can be written as \begin{eqnarray*} {\boldsymbol {\tilde{{\nabla}}}}{\boldsymbol{\cdot}}\mbox{$\boldsymbol{\hat{b}}$}= - \frac{(\mbox{$\boldsymbol{\hat{b}}$}{\boldsymbol{\cdot}}{\boldsymbol {\tilde{{\nabla}}}}) {{\rm { B}}}}{{{\rm { B}}}} + (\Gamma + 1) \left\lbrace \mbox{$\boldsymbol{\hat{b}}$}{\boldsymbol{\cdot}}\frac{{\rm d}}{{\rm d}\tau } \left(\frac{{\boldsymbol{U}}}{\Gamma +1} \right) \right\rbrace . \end{eqnarray*} (C13) The second term on the right-hand side is a relativistic correction that occurs when the comoving frame is not an inertial frame. It can be neglected in a non-relativistic flow with |$|{\boldsymbol{U}}| \ll 1$|⁠. In that same limit equations (107) through (72) reduce to known results. Defining |${\boldsymbol B } = {{\rm B}}\, \mbox{$\boldsymbol{\hat{b}}$}$| and |${\rm d}/{\rm d}t = \partial /\partial t + ({\boldsymbol{V}} {\boldsymbol{\cdot}}{\nabla})$|⁠, the comoving derivative \begin{eqnarray*} {\boldsymbol {\bar{a}} } = \frac{{\rm d} {\boldsymbol{V}}}{{\rm d}t}\,, \frac{{\rm d} {{\rm B}}}{{\rm d} t} = {{\rm B}}\, \left(\mbox{$\boldsymbol{\hat{b}}$}\, \mbox{$\boldsymbol{\hat{b}}$}{{\boldsymbol \, : \,}}{\nabla}{\boldsymbol{V}} - {\nabla}{\boldsymbol{\cdot}}{\boldsymbol{V}} \right)\,, {\nabla}{\boldsymbol{\cdot}}\mbox{$\boldsymbol{\hat{b}}$}= - \frac{(\mbox{$\boldsymbol{\hat{b}}$}{\boldsymbol{\cdot}}{\nabla}) {{\rm B}}}{{{\rm B}}} \, \Longleftrightarrow \, {\nabla}{\boldsymbol{\cdot}}{\boldsymbol B } = 0 . \end{eqnarray*} (C14) Relations (C10), (C12), and (C13) are used repeatedly in the main paper. © 2018 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices)
TI - Relativistic theory of particles in a scattering flow – I. Basic equations, diffusion, and drift
JO - Monthly Notices of the Royal Astronomical Society
DO - 10.1093/mnras/sty1449
DA - 2018-09-11
UR - https://www.deepdyve.com/lp/oxford-university-press/relativistic-theory-of-particles-in-a-scattering-flow-i-basic-FPSTcFRiG4
SP - 1747
VL - 479
IS - 2
DP - DeepDyve
ER -