8  Plasmas as Fluid

In a plasma the situation is much more complicated than that single particle orbits; the \(\mathbf{E}\) and \(\mathbf{B}\) fields are not prescribed but are determined by the positions and motions of the charges themselves. One must solve a self-consistent problem; that is, find a set of particle trajectories and field patterns such that the particles will generate the fields as they move along their orbits and the fields will cause the particles to move in those exact orbits. And this must be done in a time-varying situation. It sounds very hard, but it is not.

We have seen that a typical plasma density might be \(10^{18}\) ion–electron pairs per \(\text{m}^3\). If each of these particles follows a complicated trajectory and it is necessary to follow each of these, predicting the plasma’s behavior would be a hopeless task. Fortunately, this is not usually necessary because, surprisingly, the majority-perhaps as much as 80%-of plasma phenomena observed in real experiments can be explained by a rather crude model. This model is that used in fluid mechanics, in which the identity of the individual particle is neglected, and only the motion of fluid elements is taken into account. Of course, in the case of plasmas, the fluid contains electrical charges. In an ordinary fluid, frequent collisions between particles keep the particles in a fluid element moving together. It is surprising that such a model works for plasmas, which generally have infrequent collisions. But we shall see that there is a reason for this.

A more refined treatment-the kinetic theory of plasmas-requires more mathematical calculation. An introduction to kinetic theory is given in Chapter 12.

In some plasma problems, neither fluid theory nor kinetic theory is sufficient to describe the plasma’s behavior. Then one has to fall back on the tedious process of following the individual trajectories. Brute-force computer simulation can play an important role in filling the gap between theory and experiment in those instances where even kinetic theory cannot come close to explaining what is observed.

8.1 Definitions

Variable	Description
\(m_s\)	mass
\(q_s\)	charge
\(n_s\)	number density
\(\rho_s\)	mass density
\(\sigma_s, \rho^\ast\)	charge density (\(=n_s q_s\))
\(T_s\)	temperature
\(p_s\)	scalar pressure
\(\mathbf{u}_s\)	flow velocity
\(\mathbf{J}_s\)	current density (\(=\rho_s^\ast \mathbf{u}_s\))
\(e_s\)	internal energy
\(\phi_s\)	potential energy
\(\epsilon_s\)	total energy (\(=e_s+\frac{\mathbf{u_s}^2}{2}+\phi_s\))

where \(s\) denotes the species (e.g. \(\text{H}^+, \text{O}^+\)). Do not confuse \(\sigma\) here with conductivity. Then the total quantities without subscripts can be written as

\[ \begin{aligned} &n=\sum_s{n_s} \\ &\rho=\sum_{s}\rho_s=\sum_{s}{m_s n_s} \\ &p=\sum_s{p_s} \\ &T=\sum_s{\frac{n_s}{n}T_s} \\ &\mathbf{u}=\frac{1}{\rho}\sum_s \rho_s \mathbf{u}_s \\ &\mathbf{v}_s=\mathbf{u}_s -\mathbf{u},\quad \text{relative velocity of the } s^{th} \text{ species} \\ &\sigma = \sum_s{\sigma_s}=\sum_s{q_s n_s} \\ &\mathbf{J}=\sum_s{\mathbf{J}_s}=\sum_s{\sigma_s \mathbf{u}_s}=\sum_s{\sigma_s}\mathbf{v}_s+\mathbf{u}\sum_s{\sigma_s}=\sum_s{\sigma_s \mathbf{v}_s}+\sigma\mathbf{u}=\mathbf{J}_{cd}+ \mathbf{J}_{cv} \\ &\quad \text{where } \mathbf{J}_{cd}=\sum_s{\sigma_s}\mathbf{v}_s \text{ is the conduction current density} \\ &\quad \quad \quad \quad\mathbf{J}_{cv}=\sigma \mathbf{u} \text{ is the convection current density} \\ &\epsilon=\frac{1}{\rho}\sum_s{\rho_s\epsilon_s}=e+\frac{u^2}{2}+\phi \text{ (internal + kinetic + potential)} \end{aligned} \]

It can be easily verified that

\[ \sum_s{\rho_s\mathbf{v}_s}=0 \]

In this note we tried to use \(\mathbf{v}\) for all species/particle-based velocities, and \(\mathbf{u}\) for all bulk velocities.

We have 5 independent unknown quantities for each species: \(n_j, \mathbf{u}_j, p_j\). The pressure can be replaced by temperature \(T_j\). Additionally, we have the 6 EM field quantities: \(\mathbf{B}, \mathbf{E}\). For a 2-fluid despcription with ions and electrons, altogether we have \(5*2+6=16\) unknowns, so we need 16 equations to determine the system.

8.2 Relation of Plasma to Ordinary Electromagnetics

8.2.1 Maxwell’s Equations

In vacuum:

\[ \begin{aligned} \epsilon_0 \nabla\cdot\mathbf{E} &= \sigma \\ \nabla\times\mathbf{E} &= -\dot{\mathbf{B}} \\ \nabla\cdot\mathbf{B} &= 0 \\ \nabla\times\mathbf{B} &= \mu_0 (\mathbf{J}+\epsilon_0\dot{\mathbf{E}}) \end{aligned} \tag{8.1}\]

In a medium:

\[ \begin{aligned} \nabla\cdot\mathbf{D} &= \sigma \\ \nabla\times\mathbf{E} &= -\dot{\mathbf{B}} \\ \nabla\cdot\mathbf{B} &= 0 \\ \nabla\times\mathbf{H} &= \mathbf{J}+\dot{\mathbf{D}} \\ \mathbf{D} &= \epsilon\mathbf{E} \\ \mathbf{B} &= \mu\mathbf{H} \end{aligned} \tag{8.2}\]

\(\sigma\) and \(\mathbf{J}\) stand for the “free” charge and current densities. The “bound” charge and current densities arising from polarization and magnetization of the medium are included in the definition of the quantities \(\mathbf{D}\) and \(\mathbf{H}\) in terms of \(\epsilon\) and \(\mu\). In a plasma, the ions and electrons comprising the plasma are the equivalent of the “bound” charges and currents. Since these charges move in a complicated way, it is impractical to try to lump their effects into two constants \(\epsilon\) and \(\mu\). Consequently, in plasma physics, one generally works with the vacuum equations, in which \(\sigma\) and \(\mathbf{J}\) include all the charges and currents, both external and internal.

Note that we have used \(\mathbf{E}\) and \(\mathbf{B}\) in the vacuum equations rather than their counterparts \(\mathbf{D}\) and \(\mathbf{H}\), which are related by the constants \(\epsilon_0\) and \(\mu_0\). This is because the forces \(q\mathbf{E}\) and \(\mathbf{J}\times\mathbf{B}\) depend on \(\mathbf{E}\) and \(\mathbf{B}\) rather than \(\mathbf{D}\) and \(\mathbf{H}\), and it is not necessary to introduce the latter quantities as long as one is dealing with the vacuum equations.

8.2.2 Classical Treatment of Magnetic Materials

Since each gyrating particle has a magnetic moment, it would seem that the logical thing to do would be to consider a plasma as a magnetic material with a permeability \(\mu_m\).¹ To see why this is not done in practice, let us review the way magnetic materials are usually treated.

The ferromagnetic domains, say, of a piece of iron have magnetic moments \(\mu_i\), giving rise to a bulk magnetization \[ \mathbf{M} = \frac{1}{V}\sum_i\mu_i \] per unit volume. This has the same effect as a bound current density equal to \[ \mathbf{J}_b = \nabla\times\mathbf{M} \]

In the vacuum Ampère’s law, we must include in \(\mathbf{J}\) both this current and the “free”, or externally applied, current \(\mathbf{J}_f\):

\[ \mu_0^{-1}\nabla\times\mathbf{B} = \mathbf{J}_f + \mathbf{J}_b + \epsilon_0 \dot{\mathbf{E}} \]

We wish to write this equation in the simple form \[ \nabla\times\mathbf{H} = \mathbf{J}_f + \epsilon_0\dot{\mathbf{E}} \] by including \(\mathbf{J}_b\) in the definition of \(\mathbf{H}\). This can be done if we let \[ \mathbf{H} = \mu_0^{-1}\mathbf{B} - \mathbf{M} \]

To get a simple relation between \(\mathbf{B}\) and \(\mathbf{H}\), we assume \(\mathbf{M}\) to be proportional to \(\mathbf{B}\) or \(\mathbf{H}\): \[ \mathbf{M} = \chi_m \mathbf{H} \]

The constant \(\chi_m\) is the magnetic susceptibility. We now have \[ \mathbf{B} = \mu_0(1+\chi_m)\mathbf{H} \equiv \mu_m \mathbf{H} \]

This simple relation between \(\mathbf{B}\) and \(\mathbf{H}\) is possible because of the linear relation between \(\mathbf{M}\) and \(\mathbf{H}\).

In a plasma with a magnetic field, each particle has a magnetic moment \(\mu_\alpha\), and the quantity \(\mathbf{M}\) is the sum of all these \(\mu_\alpha\)’s in 1 \(\text{m}^3\). But now we have \[ \mu_\alpha = \frac{mv_{\perp\alpha}^2}{2B}\propto \frac{1}{B}\quad \mathbf{M}\propto \frac{1}{B} \]

The relation between \(\mathbf{M}\) and \(\mathbf{H}\) (or \(\mathbf{B}\)) is no longer linear, and we cannot write \(\mathbf{B} = \mu_m\mathbf{H}\) with \(\mu_m\) constant. It is therefore not useful to consider a plasma as a magnetic medium.

8.2.3 Classical Treatment of Dielectrics

The polarization \(\mathbf{P}\) per unit volume is the sum over all the individual moments \(\mathbf{p}_i\) of the electric dipoles. This gives rise to a bound charge density

\[ \sigma_b = -\nabla\cdot\mathbf{P} \tag{8.3}\]

In the vacuum equation, we must include both the bound charge and the free charge:

\[ \epsilon_0\nabla\cdot\mathbf{E} = \sigma_f + \sigma_b \]

We wish to write this in the simple form

\[ \nabla\cdot\mathbf{D} = \sigma_f \]

by including \(\sigma_b\) in the definition of \(\mathbf{D}\). This can be done by letting

\[ \mathbf{D} = \epsilon_0\mathbf{E} + \mathbf{P} \equiv \epsilon \mathbf{E} \]

If \(\mathbf{P}\) is linearly proportional to \(\mathbf{E}\),

\[ \mathbf{P} = \epsilon_0 \chi_e\mathbf{E} \]

then \(\epsilon\) is a constant given by

\[ \epsilon = (1+\chi_e)\epsilon_0 \]

There is no a priori reason why a relation like the above cannot be valid in a plasma, so we may proceed to try to get an expression for \(\epsilon\) in a plasma.

8.2.4 The Dielectric Constant of a Plasma

We have seen in Section 7.4 that a fluctuating \(\mathbf{E}\) field gives rise to a polarization current \(\mathbf{J}_p\). This leads, in turn, to a polarization charge given by the equation of continuity:

\[ \frac{\partial \sigma_p}{\partial t} + \nabla\cdot\mathbf{J}_p = 0 \]

This is the equivalent of Equation 8.3, except that, as we noted before, a polarization effect does not arise in a plasma unless the electric field is time varying. Since we have an explicit expression for \(\mathbf{J}_p\) but not for \(\sigma_p\), it is easier to work with the Ampère’s law:

\[ \nabla\times\mathbf{B} = \mu_0(\mathbf{J}_f +\mathbf{J}_p + \epsilon\dot{\mathbf{E}}) \]

We wish to write this in the form

\[ \nabla\times\mathbf{B} = \mu_0(\mathbf{J}_f + \epsilon\dot{\mathbf{E}}) \]

This can be done if we let

\[ \epsilon = \epsilon_0 + \frac{j_p}{\dot{E}} \]

From Equation 7.12 for \(\mathbf{J}_p\), we have

\[ \epsilon = \epsilon_0 + \frac{\rho}{B^2}\quad \text{or}\quad \epsilon_R\equiv \frac{\epsilon}{\epsilon_0} = 1+\frac{\mu_0\rho c^2}{B^2} \tag{8.4}\]

This is the low-frequency plasma dielectric constant for transverse motions. The qualifications are necessary because our expression for \(\mathbf{J}_p\) is valid only for \(\omega^2\ll\omega_c^2\) and for \(\mathbf{E}\) perpendicular to \(\mathbf{B}\). The general expression for \(\epsilon\), of course, is very complicated and hardly fits on one page.

Note that as \(\rho\rightarrow 0\), \(\epsilon_R\) approaches its vacuum value, unity, as it should. As \(B\rightarrow\infty\), \(\epsilon_R\) also approaches unity. This is because the polarization drift \(\mathbf{v}_p\) then vanishes, and the particles do not move in response to the transverse electric field. In a usual laboratory plasma, the second term in Equation 8.4 is large compared with unity. For instance, if \(n=10^{16}\,\text{m}^{-3}\) and \(B=0.1\,\text{T}\) we have (for hydrogen) \[ \frac{\mu_0\rho c^2}{B^2} = \frac{(4\pi\times 10^{-7})(10^{16})(1.67\times 10^{-27})(3\times 10^8)^2}{(0.1)^2} = 189 \]

This means that the electric fields due to the particles in the plasma greatly alter the fields applied externally. A plasma with large \(\epsilon\) shields out alternating fields, just as a plasma with small \(\lambda_D\) shields out dc fields.

8.3 Fluid Equations

In (Chen 2016), this is introduced in the diffusion chapter 5.7, which is a bit weird. In (Bellan 2008), this is introduced in chapter 2 where Vlasov equation is first derived and then follows the simplifications which lead to 2-fluid and MHD.

8.3.1 Equation of Continuity

The integral form of mass conservation for each species is \[ \frac{\mathrm{d}}{\mathrm{d}t}\int_V \rho_s dx^3=0 \]

The conservation of matter requires that the total number of particles \(N_s\) in a volume V can change only if there is a net flux of particles across the surface S bounding that volume. Since the particle flux density is \(n_s\mathbf{u}_s\), we have, by the divergence theorem, \[ \frac{\partial N_s}{\partial t} = \int_V\frac{\partial n_s}{\partial t}dV = -\oint n\mathbf{u}_s\cdot \mathrm{d}\mathbf{S} = - \int_V \nabla\cdot(n_s\mathbf{u}_s)dV \]

Since this must hold for any volume V, the integrands must be equal: \[ \frac{\partial n_s}{\partial t} + \nabla\cdot(n_s\mathbf{u}_s) = 0 \tag{8.5}\]

There is one such equation of continuity for each species. Any sources or sinks of particles are to be added to the right-hand side.

8.3.2 Momentum Equation

Maxwell’s equations tell us what \(\mathbf{E}\) and \(\mathbf{B}\) are for a given state of the plasma. To solve the self-consistent problem, we must also have an equation giving the plasma’s response to given \(\mathbf{E}\) and \(\mathbf{B}\). In the fluid approximation, we consider the plasma to be composed of two or more interpenetrating fluids, one for each species. In the simplest case, when there is only one species of ion, we shall need two equations of motion, one for the positively charged ion fluid and one for the negatively charged electron fluid. In a partially ionized gas, we shall also need an equation for the fluid of neutral atoms. The neutral fluid will interact with the ions and electrons only through collisions. The ion and electron fluids will interact with each other even in the absence of collisions, because of the \(\mathbf{E}\) and \(\mathbf{B}\) fields they generate.

The equation of motion for a single particle is \[ m\frac{\mathrm{d}\mathbf{v}}{\mathrm{d}t} = q(\mathbf{E}+\mathbf{v}\times\mathbf{B}) \tag{8.6}\]

Assume first that there are no collisions and no thermal motions. Then all the particles in a fluid element move together, and the average velocity \(\mathbf{u}\) of the particles in the element is the same as the individual particle velocity \(\mathbf{v}\). The fluid equation is obtained simply by multiplying Equation 8.6 by the density \(n\): \[ mn\frac{\mathrm{d}\mathbf{u}}{\mathrm{d}t} = qn(\mathbf{E}+\mathbf{u}\times\mathbf{B}) \tag{8.7}\]

This is, however, not a convenient form to use. In Equation 8.6, the time derivative is to be taken at the position of the particles. On the other hand, we wish to have an equation for fluid elements fixed in space, because it would be impractical to do otherwise. Consider a drop of cream in a cup of coffee as a fluid element. As the coffee is stirred, the drop distorts into a filament and finally disperses all over the cup, losing its identity. A fluid element at a fixed spot in the cup, however, retains its identity although particles continually go in and out of it.

To make the transformation to variables in a fixed frame, consider \(\mathbf{G}(x,t)\) to be any property of a fluid in one-dimensional x space. The change of \(\mathbf{G}\) with time in a frame moving with the fluid is the sum of two terms: \[ \frac{\mathrm{d}\mathbf{G}(x,t)}{\mathrm{d}t} = \frac{\partial \mathbf{G}}{\partial t} + \frac{\partial \mathbf{G}}{\partial x}\frac{\partial x}{\partial t} = \frac{\partial \mathbf{G}}{\partial t} + u_x\frac{\partial \mathbf{G}}{\partial x} \]

The first term on the right represents the change of \(\mathbf{G}\) at a fixed point in space, and the second term represents the change of G as the observer moves with the fluid into a region in which \(\mathbf{G}\) is different. In three dimensions, this generalizes to \[ \frac{\mathrm{d}\mathbf{G}}{\mathrm{d}t} = \frac{\partial \mathbf{G}}{\partial t} + (\mathbf{u}\cdot\nabla)\mathbf{G} \]

This is called the convective derivative (Equation 3.2).

In the case of a plasma, we take \(\mathbf{G}\) to be the fluid velocity \(\mathbf{u}\) and write Equation 8.7 as \[ mn\Big[\frac{\partial\mathbf{u}}{\partial t} + (\mathbf{u}\cdot\nabla)\mathbf{u}\Big] = qn(\mathbf{E}+\mathbf{u}\times\mathbf{B}) \]

where \(\partial\mathbf{u}/\partial t\) is the time derivative in a fixed frame.

Stress Tensor –> scalar pressure

\[ mn\Big[\frac{\partial\mathbf{u}}{\partial t} + (\mathbf{u}\cdot\nabla)\mathbf{u}\Big] = qn(\mathbf{E}+\mathbf{u}\times\mathbf{B}) - \nabla p \tag{8.8}\]

What we have shown here is only a special case: the transfer of \(x\) momentum by motion in the \(x\) direction; and we have assumed that the fluid is isotropic, so that the same result holds in the y and z directions. But it is also possible to transfer \(y\) momentum by motion in the \(x\) direction, for instance. This shear stress cannot be represented by a scalar \(p\) but must be given by a tensor \(\mathbf{P}\), the stress tensor, whose components \(P_{ij}= mn\overline{v_i v_j}\) specify both the direction of motion and the component of momentum involved. In the general case the term \(\nabla p\) is replaced by \(-\nabla\cdot\mathbf{P}\).

When the distribution function is an isotropic Maxwellian, \(\mathbf{P}\) is written \[ \mathbf{P} = \begin{pmatrix}p & 0 & 0 \\ 0 & p & 0 \\ 0 & 0 & p\end{pmatrix} \]

\(-\nabla\cdot\mathbf{P} = \nabla p\). A plasma could have two temperatures \(T_\perp\) and \(T_\parallel\) in the presence of a magnetic field. In that case, there would be two pressures \(T_\perp\) and \(T_\parallel\) in the presence of a magnetic field. In that case, there would be two pressure \(p_\perp = nk_B T_\perp\) and \(p_\parallel = nk_B T_\parallel\). The stress tensor is then \[ \mathbf{P} = \begin{pmatrix}p_\perp & 0 & 0 \\ 0 & p_\perp & 0 \\ 0 & 0 & p_\parallel\end{pmatrix} \]

where the coordinate of the third row or column is the direction of \(\mathbf{B}\). This is still diagonal and shows isotropy in a plane perpendicular to \(\mathbf{B}\).

In an ordinary fluid, the off-diagonal elements of \(\mathbf{P}\) are usually associated with viscosity. When particles make collisions, they come off with an average velocity in the direction of the fluid velocity \(\mathbf{u}\) at the point where they made their last collision. This momentum is transferred to another fluid element upon the next collision. This tends to equalize \(\mathbf{u}\) at different points, and the resulting resistance to shear flow is what we intuitively think of as viscosity. The longer the mean free path, the farther momentum is carried, and the larger is the viscosity. In a plasma there is a similar effect which occurs even in the absence of collisions. The Larmor gyration of particles (particularly ions) brings them into different parts of the plasma and tends to equalize the fluid velocities there. The Larmor radius rather than the mean free path sets the scale of this kind of collisionless viscosity. It is a finite-Larmor-radius effect which occurs in addition to collisional viscosity and is closely related to the \(\mathbf{v}_E\) drift in a nonuniform \(\mathbf{E}\) field (Equation 7.11).

8.3.2.1 Collisions

If there is a neutral gas, the charged fluid will exchange momentum with it through collisions. The momentum lost per collision will be proportional to the relative velocity \(\mathbf{u}-\mathbf{u}_0\), where \(\mathbf{u}_0\) is the velocity of the neutral fluid. If \(\tau\), the mean free time between collisions, is approximately constant, the resulting force term can be roughly written as \(-mn(\mathbf{u}-\mathbf{u}_0)/\tau\). The equation of motion can be generalized to include anisotropic pressure and neutral collisions as follows: \[ mn\Big[\frac{\partial\mathbf{u}}{\partial t} + (\mathbf{u}\cdot\nabla)\mathbf{u}\Big] = qn(\mathbf{E}+\mathbf{u}\times\mathbf{B}) - \nabla\cdot\overleftrightarrow{P} - \frac{mn(\mathbf{u}-\mathbf{u}_0)}{\tau} \tag{8.9}\]

This can also be written as (including the pressure term and other forces) \[ \rho\frac{\mathrm{d}\mathbf{u}}{\mathrm{d}t}=(\rho^\ast \mathbf{E}+\mathbf{J}\times\mathbf{B})-\nabla\cdot\overleftrightarrow{P}+\mathbf{f}_n \]

or, in an equivalent conservative form, \[ \frac{\partial (\rho\mathbf{u})}{\partial t}+\nabla\cdot(\rho\mathbf{u}\mathbf{u})=(\rho^\ast \mathbf{E}+\mathbf{J}\times\mathbf{B})-\nabla\cdot\overleftrightarrow{P}+\mathbf{f}_n \]

Collisions between charged particles have not been included; these will be discussed in sec-collision ADD IT!.

8.3.2.2 Comparison with Ordinary Hydrodynamics

Ordinary fluids obey the Navier–Stokes equation

\[ \rho\Big[ \frac{\partial\mathbf{u}}{\partial t} + (\mathbf{u}\cdot\nabla)\mathbf{u} \Big] = -\nabla p + \rho\nu\nabla^2\mathbf{u} \tag{8.10}\]

This is the same as Equation 8.9 except for the absence of electromagnetic forces and collisions between species (there being only one species). The viscosity term \(\rho\nu\nabla^2\mathbf{u}\), where \(\nu\) is the kinematic viscosity coefficient, is just the collisional part of \(\nabla\cdot\overleftrightarrow{P} - \nabla p\) in the absence of magnetic fields. Equation 8.10 describes a fluid in which there are frequent collisions between particles. Equation 8.9, on the other hand, was derived without any explicit statement of the collision rate. Since the two equations are identical except for the \(\mathbf{E}\) and \(\mathbf{B}\) terms, can Equation 8.9 really describe a plasma species? The answer is a guarded yes, and the reasons for this will tell us the limitations of the fluid theory. This is extremely important to clarify — there are still quite many people think that collision is assumed for the fluid theory, therefore do not believe that MHD can be used to describe plasmas.

In the derivation of Equation 8.9, we did actually assume implicitly that there were many collisions. This assumption came in the derivation of the pressure tensor (Section 12.4) when we took the velocity distribution to be Maxwellian. Such a distribution generally comes about as the result of frequent collisions. However, this assumption was used only to take the average of \(v^2\). Any other distribution with the same average would give us the same answer. The fluid theory, therefore, is not very sensitive to deviations from the Maxwellian distribution, but there are instances in which these deviations are important. Kinetic theory must then be used.

There is also an empirical observation by Irving Langmuir which helps the fluid theory. In working with the electrostatic probes which bear his name, Langmuir discovered that the electron distribution function was far more nearly Maxwellian than could be accounted for by the collision rate. This phenomenon, called Langmuir’s paradox, has been attributed at times to high-frequency oscillations. There has been no satisfactory resolution of the paradox, but this seems to be one of the few instances in plasma physics where nature works in our favor.

Another reason the fluid model works for plasmas is that the magnetic field, when there is one, can play the role of collisions in a certain sense. When a particle is accelerated, say by an \(\mathbf{E}\) field, it would continuously increase in velocity if it were allowed to free-stream. When there are frequent collisions, the particle comes to a limiting velocity proportional to \(\mathbf{E}\). The electrons in a copper wire, for instance, drift together with a velocity \(\mathbf{v}=\mu\mathbf{E}\), where \(\mu\) is the mobility. A magnetic field also limits free-streaming by forcing particles to gyrate in Larmor orbits. The electrons in a plasma also drift together with a velocity proportional to \(\mathbf{E}\), namely, \(\mathbf{v}_E=\mathbf{E}\times\mathbf{B}/B^2\). In this sense, a collisionless plasma behaves like a collisional fluid. Of course, particles do free-stream along the magnetic field, and the fluid picture is not particularly suitable for motions in that direction. For motions perpendicular to \(\mathbf{B}\), the fluid theory is a good approximation.

8.3.3 Equation of State

One more relation is needed to close the system of equations. A complete description requires the equation of energy. However, the simplest way is to use the thermodynamic equation of state relating \(p\) to \(n\): \[ p = C\rho^\gamma \tag{8.11}\]

where \(C\) is a constant and \(\gamma\) is the ratio of specific heats \(C_p/C_\nu\). The term \(\nabla p\) is therefore given by \[ \frac{\nabla p}{p} = \gamma \frac{\nabla n}{n} \tag{8.12}\]

For isothermal compression, we have \[ \nabla p = \nabla(nk_B T) = k_B T\nabla n \]

so that, clearly, \(\gamma=1\). For adiabatic compression, \(k_B T\) will also change, giving \(\gamma\) a value larger than one. If \(N\) is the number of degrees of freedom, \(\gamma\) is given by \[ \gamma = (2+N)/N \]

The validity of the equation of state requires that heat flow be negligible; that is, that thermal conductivity be low. Again, this is more likely to be true in directions perpendicular to \(\mathbf{B}\) than parallel to it. Fortunately, most basic phenomena can be described adequately by the crude assumption of the equation of state.

8.4 Two-Fluid Model

Besides the Vlasov theory (Chapter 12), we can apply the simpler but equally powerful 2-fluid model, in which electrons and ions are treated as two different fluids. Depending on the form of the pressure term, we have the 5/6/10-moment equations.

8.4.1 Five-moment

Five-moment refers to \((\rho, \mathbf{u}, p)\). Under either adiabatic or isothermal assumptions, we can get rid of the pressure/energy equation and replace it with an equation of state. For simplicity, let the plasma have only two species: ions and electrons; extension to more species is trivial. The charge and current densities are then given by \[ \begin{aligned} \sigma &= n_i q_i + n_e q_e \\ \mathbf{J} &= n_i q_i \mathbf{v}_i + n_e q_e \mathbf{v}_e \end{aligned} \]

We shall neglect collisions and viscosity. The complete equation set is: \[ \begin{aligned} \epsilon_0\nabla\cdot\mathbf{E} = n_i q_i + n_e q_e \\ \nabla\times\mathbf{E} = -\dot{\mathbf{B}} \\ \nabla\cdot\mathbf{B} = 0 \\ \mu_0^{-1}\nabla\times\mathbf{B} = n_i q_i \mathbf{u}_i + n_e q_e \mathbf{u}_e + \epsilon_0\dot{\mathbf{E}} \\ m_jn_j\Big[ \frac{\partial\mathbf{u}_j}{\partial t}+(\mathbf{u}_j\cdot\nabla)\mathbf{u}_j\Big] = q_j n_j (\mathbf{E}+\mathbf{u}_j\times\mathbf{B}) - \nabla p_j\quad j=i,e \\ \frac{\partial n_j}{\partial t} + \nabla\cdot(n_j\mathbf{u}_j) = 0\quad j=i,e \\ p_j = C_j n_j^\gamma\quad j=i,e \end{aligned} \tag{8.13}\] where \(C_j\) is a constant and \(\gamma\) is the adiabatic index with the approximations in Table 8.1.

Table 8.1: Equation of state approximations.

Regime	Equation of state	Name
\(v_\mathrm{ph} \gg v_{ts}\)	\(p_s\sim n_s^\gamma\)	adiabatic
\(v_\mathrm{ph} \ll v_{ts}\)	\(p_s=n_s k_B T_s\), \(T_s=\mathrm{const.}\)	isothermal

There are 16 scalar unknowns: \(n_i, n_e, p_i, p_e, \mathbf{u}_i, \mathbf{u}_e, \mathbf{E}\), and \(\mathbf{B}\). There are apparently 18 scalar equations if we count each vector equation as three scalar equations. However, two of Maxwell’s equations are superfluous, since the two of the equations can be recovered from the divergences of the other two. The simultaneous solution of this set of 16 equations in 16 unknowns gives a self-consistent set of fields and motions in the fluid approximation.

8.4.2 Six-moment

For monatomic gases, the six-moment equations for all charged fluids (indexed by \(s\)) can be written as: \[ \begin{aligned} \frac{\partial \rho_s}{\partial t} + \nabla\cdot(\rho_s\mathbf{u}_s) = 0 \\ \frac{\partial (\rho_s\mathbf{u}_s)}{\partial t} + \nabla\cdot\left[ \rho_s\mathbf{u}_s\mathbf{u}_s + p_{s\perp}\pmb{I} + (p_{s\parallel} - p_{s\perp})\hat{b}\hat{b} \right] = \frac{q_s}{m_s}\rho_s(\mathbf{E} + \mathbf{u}_s\times\mathbf{B}) \\ \frac{\partial p_{s\parallel}}{\partial t} + \nabla\cdot(p_{s\parallel}\mathbf{u}_s) = -2p_{s\parallel}\hat{b}\cdot(\hat{b}\cdot\nabla)\mathbf{u}_s \\ \frac{\partial p_{s\perp}}{\partial t} + \nabla\cdot(p_{s\perp}\mathbf{u}_s) = -p_{s\perp}(\nabla\cdot\mathbf{u}_s) + p_{s\perp}\hat{b}\cdot(\hat{b}\cdot\nabla)\mathbf{u}_s \end{aligned} \tag{8.14}\]

The two pressure equations in Equation 8.14 can be combined to give the equation for the average pressure \(p=(2p_\perp + p_\parallel)/3\): \[ \frac{\partial p_s}{\partial t} + \nabla\cdot(p_s\mathbf{u}_s) = (p_s - p_{s\parallel})\hat{b}\cdot(\hat{b}\cdot\nabla)\mathbf{u}_s - \left( p_s - \frac{p_{s\parallel}}{3} \right)\nabla\cdot\mathbf{u}_s \tag{8.15}\]

Alternatively, we can solve for the hydrodynamic energy density \(e = \frac{n\mathbf{u}^2}{2} + \frac{3}{2}p\) for each species: \[ \frac{\partial e_s}{\partial t} + \nabla\cdot\left[ \mathbf{u}_s(e_s+p_s) + \mathbf{u}_s\cdot(p_{s\parallel} - p_{s\perp})\hat{b}\hat{b} \right] = \frac{q_s}{m_s}\rho_s \mathbf{u}_s\cdot\mathbf{E} \]

which is more of a conservative form and can be beneficial to get better jump conditions across shock waves. Note, however, that the parallel pressure equation is still solved with the adiabatic assumption, so non-adiabatic heating is not properly captured. In addition, the magnetic energy is not included into the energy density, so the jump conditions are only approximate (I think this is more of a numerical consideration). In general, there can be many more source terms on the right hand sides of the above equations corresponding to gravity, charge exchange, chemical reactions, collisions, etc.

The electric field \(\mathbf{E}\) and magnetic field \(\mathbf{B}\) are obtained from Maxwell’s equations: \[ \begin{aligned} \frac{\partial \mathbf{B}}{\partial t} + \nabla\times\mathbf{E} &= 0 \\ \frac{\partial \mathbf{E}}{\partial t} - c^2\nabla\times\mathbf{B} &= -c^2\mu_0\mathbf{J} \\ \nabla\cdot\mathbf{E} &= \frac{\rho_c}{\epsilon_0} \\ \nabla\cdot\mathbf{B} &= 0 \end{aligned} \] where \(\rho_c = \sum_s (q_s/m_s)\rho_s\) is the total charge density and \(\mathbf{J}=\sum_s(q_s/m_s)\rho_s\mathbf{u}_s\) is the current density.

The last two equations are constraints on the initial conditions; these are not guaranteed to hold numerically. Classical tricks involves using hyperbolic/parabolic cleaning or a facially-collocated Yee-type mesh.

8.4.3 Five/Ten-moment

(Wang et al. 2020)

• Continuity equation for each species \[ \frac{\partial(m_s n_s)}{\partial t} + \frac{\partial(m_sn_su_{j,s})}{\partial x_j} = 0 \]

• Momentum equation for each species \[ n_s m_s\frac{\mathrm{d}\mathbf{u}_s}{\mathrm{d}t} = n_s q_s(\mathbf{E}+\mathbf{u}\times\mathbf{B}) - \nabla\cdot\overleftrightarrow{P}_s - \mathbf{R}_{s\alpha} \] or in Einstein’s notation, \[ \begin{aligned} \frac{\partial(m_s n_s u_{i,s})}{\partial t} = n_s q_s (E_i + \epsilon_{ijk}u_{j,s}B_k) -\frac{\partial \mathcal{P}_{ij,s}}{\partial x_j} - R_{i,s\alpha} \end{aligned} \] where \(\epsilon_{ijk}\) is the Levi-Civita symbol. The moments are defined as \[ \begin{aligned} n_s(\mathbf{x}) &\equiv \int f_s \mathrm{d}\mathbf{v} \\ u_{i,s}(\mathbf{x}) &\equiv \frac{1}{n_s(\mathbf{x})}\int v_i f_s \mathrm{d}\mathbf{v} \\ \mathcal{P}_{ij,s}(\mathbf{x}) &\equiv m_s \int v_i v_j f_s \mathrm{d}\mathbf{v} \end{aligned} \] with \(f_s(\mathbf{x},\mathbf{v},t)\) being the phase space distribution function. We will neglect the subscript \(s\) hereinafter for convenience. For completeness, \(\mathcal{P}_{ij}\) relates to the more familiar thermal pressure tensor \[ P_{ij} \equiv m\int(v_i - u_i)(v_i - u_i)f\mathrm{d}\mathbf{v} \] by \[ \mathcal{P}_{ij} = P_{ij} + nm u_i u_j \]

For simplicity, non-ideal effects like viscous dissipation are neglected. The electric and magnetic fields \(\mathbf{E}\) and \(\mathbf{B}\) are evolved using Maxwell equations \[ \begin{aligned} \frac{\partial\mathbf{B}}{\partial t} + \nabla\times\mathbf{E} = 0 \\ \frac{\partial\mathbf{E}}{\partial t} - c^2\nabla\times\mathbf{B} = -\frac{1}{\epsilon_0}\sum_s q_s n_s \mathbf{u}_s \\ \nabla\cdot\mathbf{B} = 0 \\ \nabla\cdot\mathbf{E} = \epsilon_0^{-1}\sum_s n_s q_s \end{aligned} \] with \(c=1/\sqrt(\mu_0\epsilon_0)\) being the speed of light.

To close the system, the second order moment \(\mathcal{P}_{ij}\) or \(P_{ij}\) must be specified. For example, a cold fluid closure simply sets \(P_{ij}=0\), while an isothermal equation of state (EOS) assumes that the temperature is constant. Or, assuming zero heat flux and that the pressure tensor is isotropic, we can write an adiabatic EOS for \(P_{ij} = p I_{ij}\), \[ \frac{\partial\mathcal{\epsilon}}{\partial t} + \nabla\cdot\left[ (p+\mathcal{\epsilon})\mathbf{u} \right] = nq\mathbf{u}\cdot\mathbf{E} \] where \[ \mathcal{\epsilon} \equiv \frac{p}{\gamma-1} + \frac{1}{2}mn\mathbf{u}^2 \] is the total fluid (thermal + kinetic) energy and \(\gamma\) is the adiabatic index, setting to \(5/3\) for a fully ionized plasma. For a plasma with \(S\) species (\(s=1,...,S\)) this system is closed and has a total of \(5S+6\) equations, and are here referred to as the five-moment model. More general models can be obtained by retaining the evolution equations for all six components of the pressure tensor in the so-called ten-moment model \[ \frac{\mathcal{P}_{ij,s}}{\partial t} + \frac{\partial \mathcal{Q}_{ijm,s}}{\partial x_m} = n_sq_su_{i,s}E_j + \frac{q_s}{m_s}\epsilon_{iml}\mathcal{P}_{mj,s}B_l \tag{8.16}\] where the third moment \[ \mathcal{Q}_{ijm,s}(\mathbf{x}) \equiv m_s\int v_i v_j v_m f_s\mathrm{d}\mathbf{v} \] relates to the heat flux tensor defined in the fluid frame \[ Q_{ijm} \equiv m\int (v_i-u_i)(v_j-u_j)(v_m-u_m)\mathrm{d}\mathbf{v} \] by \[ \mathcal{Q}_{ijm} = Q_{ijm} + u_i\mathcal{P}_{jm} - 2nm u_i u_j u_m \]

Again, the equations here must be closed by some approximation for the heat-flux tensor. Another option is to include evolution equations for even higher order moments, e.g., the ten independent components of the heat-flux tensor.

Theoretically, the multifluid-Maxwell equations approach the Hall magnetohydrodynamics (MHD) under asymptotic limits of vanishing electron mass (\(m_e \rightarrow 0\)) and infinite speed of light (\(c\rightarrow\infty\)). All waves and effects within the two-fluid picture are retained, for example, the light wave, electron and ion inertial effects like the ion cyclotron wave and whistler wave. Particularly, through properly devised heat-flux closures, the ten-moment model could partially capture nonlocal kinetic effects like Landau damping, in a manner similar to the gyrokinetic models.

Although multifluid-Maxwell models provide a more complete description of the plasma than reduced, asymptotic models like MHD, they are less frequently used. The reason for this is the fast kinetic scales involved. Retaining the electron inertia adds plasma-frequency and cyclotron time-scale, while non-neutrality adds Debye length spatial-scales. Further, inclusion of the displacement currents means that EM waves must be resolved when using an explicit scheme. Fortunately, the restrictions due to kinetic scales are introduced only through the non-hyperbolic source terms of the momentum equation, the Ampère’s law, and the pressure equation. Therefore we may eliminate these restrictions by updating the source term separately either exactly or using an implicit algorithm (Note: BATSRUS applies the point-implicit scheme.). This allows larger time steps and leads to significant speedup, especially with realistic electron/ion mass ratios. The speed of light constraint still exists, however, can be greatly relaxed, using reduced values for the speed of light and/or sub-cycling Maxwell equations. Of course, an implicit Maxwell solver, or a reduced set of electromagnetic equations like the Darwin approximation², can also relax the time-step restrictions.

8.4.4 Characteristic wave speeds

The fastest wave speed in the 5/6/10-moment equations is the speed of light \(c\). Theoretically, the characteristic wave speeds shall be consistent and fallback to MHD in the isotropic 5-moment case.

Numerically, however, using \(c\) in the numerical fluxes makes the scheme rather diffusive (?). The discretization time step is limited by the Courant-Friedrichs-Lewy (CFL) condition based on the speed of light. In practice, we can reduce the speed of light to a value that is a factor of 2-3 faster than the fastest flow and fast wave speed to speed up the simulation.

8.5 Generalized Ohm’s Law

The generalized Ohm’s law can be derived from two-fluid equations. Landau once pointed out that the validity of describing plasma behaviors with two-fluid equations is based on the fact that the time scale of ions and electrons to reach an equilibrium Maxwellian distribution is much smaller than that of the heat exchange between species. The Ohm’s law is established upon the collisions between electrons, ions and neutral species, and is used to describe the relation between \(\mathbf{J}\) and \(\mathbf{E},\,\mathbf{B},\,\mathbf{u}\).

8.5.1 Basic definitions and assumptions

Assume plasma consists of electrons and ionized \(\textrm{H}^+\). We use \(e,\, i,\, a\) as subscripts for electron, ion and neutral species respectively. \[ \begin{aligned} q_e=-q_i=-e \\ m_e\ll m_i \approx m_a \end{aligned} \]

With the definition of ionization degree in Equation 2.1, \[ \alpha=\frac{n_i}{n_i+n_a} \tag{8.17}\] we have \[ \begin{aligned} \frac{n_i}{n_a}&=\frac{\alpha}{1-\alpha}\\ \frac{n_i}{2n_i+n_a}&=\frac{\alpha}{1+\alpha}\\ \frac{n_a}{2n_i+n_a}&=\frac{1-\alpha}{1+\alpha} \end{aligned} \]

The quasi-neutrality assumption gives \[ \begin{aligned} \rho_e^\ast+\rho_i^\ast =0\ \Rightarrow\ n_e=n_i \\ \rho_e=m_e n_e \ll \rho_i=m_i n_i \end{aligned} \]

Assume a system of ideal gas in thermal equilibrium, \[ \begin{aligned} T_e=T_i=T_a \\ p=p_e+p_i+p_a \end{aligned} \]

Using the ionization degree defined in Equation 8.17, we have \[ \begin{aligned} p_e = p_i = \frac{n_i}{2n_i+n_a} p = \frac{\alpha}{1+\alpha} p \\ p_a = \frac{n_a}{2n_i+n_a}p = \frac{1-\alpha}{1+\alpha}p \end{aligned} \] Also, we have \[ v_{e,th}\gg v_{i,th} \]

Assume that the flow velocity is much smaller than the thermal speed, \[ |\mathbf{v}_s|\ll |\mathbf{v}_{s,th}|,\quad s=e,\, i,\, a. \]

Typical time scale must satisfy \(\tau_0\gg \tau_{\text{collision}}\).

Finally, we ignore viscosity and non-EM force in the derivation.

The average velocity of protons with respect to bulk velocity of the plasma, the velocity difference between electron and ion, and the average velocity of neutral particles are \[ \begin{aligned} \mathbf{v}_i&=\mathbf{u}_i-\mathbf{u} \\ \mathbf{v}_{ei}&=\mathbf{u}_e-\mathbf{u}_i \\ \mathbf{u}_a&=\mathbf{u}-\frac{n_i}{n_a}\left(\mathbf{v}_i+\frac{m_e}{m_i}\mathbf{v}_{ei}\right) \end{aligned} \]

The average momenta of electrons relative to protons and neutral particles are \[ \begin{aligned} \mathbf{I}_{ei}&=m_e\mathbf{v}_{ei} \\ \mathbf{I}_{ea}&=m_e\left[\left( 1+\frac{n_im_e}{n_a m_i}\right)\mathbf{v}_{ei}+\left( 1+\frac{n_i}{n_a}\right)\mathbf{v}_i\right]\approx m_e\left[ \mathbf{v}_{ei}+\left(\frac{1}{1-\alpha}\right)\mathbf{v}_i \right] \\ \mathbf{I}_{ia}&=m_i(\mathbf{u}_i-\mathbf{u}_a)=m_i (\mathbf{v}_i+\mathbf{u}-\mathbf{u}_a)=\frac{1}{1-\alpha}m_i\mathbf{v}_i+\frac{\alpha}{1-\alpha}m_e\mathbf{v}_{ei} \end{aligned} \]

Note that from Newton’s third law, \(\Delta \mathbf{I}_{sk}=-\Delta \mathbf{I}_{ks}\).

8.5.2 Equation of motion for each species

8.5.2.1 For electrons

\[ \begin{aligned} m_en_e\frac{d\mathbf{u}_e}{dt}&=m_en_i\frac{d}{dt}\left(\mathbf{v}_{ei}+\mathbf{v}_i+\mathbf{u}\right) \\ &=-\nabla p_e - n_i e \left[ \mathbf{E}+(\mathbf{v}_{ei}+\mathbf{v}_i+\mathbf{u})\times\mathbf{B} \right]+\mathbf{f}_e^c \end{aligned} \tag{8.18}\] where the collisional force \[ \mathbf{f}_e^c=-n_i m_e\mathbf{v}_{ei}\left(\tau_{ei}^{-1}+\tau_{ea}^{-1}\right)-\frac{1}{1-\alpha}n_i m_e\mathbf{v}_i\tau_{ea}^{-1} \tag{8.19}\]

Proof:

\(\mathbf{f}_{sk}^c\cdot\tau_{sk}=n_s\Delta \mathbf{I}_{sk}\), where \(\tau_{sk}\) is the average collision interval. \[ \Delta \mathbf{I}_{sk}=-\Delta \mathbf{I}_{ks} \]

Assuming elastic collision, we have \[ \begin{aligned} \Delta\mathbf{I}_{sk}&=-\frac{m_sm_k}{m_s+m_k}\mathbf{v}_{sk} \\ \Rightarrow \Delta\mathbf{I}_{ek}&=-\frac{m_k}{m_e+m_k}m_e\mathbf{v}_{ek}\approx -\mathbf{I}_{ek}\text{(ignore electron mass)}\\ \Delta\mathbf{I}_{ia}&=-\frac{m_a}{m_a+m_i}m_i\mathbf{v}_{ia}\approx -\frac{1}{2}\mathbf{I}_{ia}\text{(ion mass=neutron mass)} \end{aligned} \]

Therefore the collision force of electron is \[ \mathbf{f}_e^c=-n_i\mathbf{I}_{ei}\tau_{ei}^{-1}-n_i\mathbf{I}_{ea}\tau_{ea}^{-1}=-n_i m_e\mathbf{v}_{ei}\left(\tau_{ei}^{-1}+\tau_{ea}^{-1}\right)-\frac{1}{1-\alpha}n_i m_e\mathbf{v}_i\tau_{ea}^{-1} \] The collision force between ion and neutral is \[ \mathbf{f}_{ia}^c=-\frac{1}{2}n_i\mathbf{I}_{ia}\tau_{ia}^{-1} \] □

8.5.2.2 For protons

\[ m_i n_i\frac{d\mathbf{u}_i}{dt}=m_i n_i\frac{d}{dt}(\mathbf{v}_i+\mathbf{u})=-\nabla p_i+n_i e\left[ \mathbf{E}+\left(\mathbf{v}_i+\mathbf{u}\times\mathbf{B}\right) \right]+\mathbf{f}_i^c \tag{8.20}\] where the collision force is \[ \mathbf{f}_{i}^c=\mathbf{f}_{ia}^c+\mathbf{f}_{ie}^c=n_im_e\mathbf{v}_{ei}\tau_{ei}^{-1}-\frac{\alpha}{2(1-\alpha)}n_im_e\mathbf{v}_{ei}\tau_{ia}^{-1}-\frac{1}{2(1-\alpha)}n_i m_i\mathbf{v}_i\tau_{ia}^{-1} \tag{8.21}\]

8.5.2.3 For the whole plasma

If we do not consider the effect of other forces on plasma and ignore electron mass (\(m_en_e\ll m_in_i\)), \[ m_i(n_i+n_a)\frac{d\mathbf{u}}{dt}=-\nabla p+\mathbf{J}\times\mathbf{B}=-\nabla p-n_i e\mathbf{v}_{ei}\times\mathbf{B} \tag{8.22}\] where the total current density is \[ \mathbf{J}=\sum_s{n_s q_s\mathbf{u}_s}=-n_i e(\mathbf{v}_{ei}+\mathbf{v}_i+\mathbf{u})+n_i e(\mathbf{v}_i+\mathbf{u})=-n_i e\mathbf{v}_{ei} \] Physically we usually think of current in plasma as the velocity difference between electrons and ions.

8.5.3 Preconditions

If an external B-field is present and the conductor is not at rest but moving at velocity \(\mathbf{u}\), then an extra term must be added to account for the current induced by the Lorentz force on the charge carriers: \[ \mathbf{J}=\sigma(\mathbf{E}+\mathbf{u}\times\mathbf{B}) \] where \(\sigma\) is the electric conductivity.

In the rest frame of the moving conductor this term drops out because \(\mathbf{u}=0\). There is no contradiction because the electric field in the rest frame differs from the E-field in the lab frame: \(\mathbf{E}^\prime = \mathbf{E} + \mathbf{u}\times\mathbf{B}\). Electric and magnetic fields are relative and related by the Lorentz transformation.

If the current \(\mathbf{J}\) is alternating because the applied voltage or E-field varies in time, then reactance must be added to resistance to account for self-inductance (see electrical impedance). The reactance may be strong if the frequency is high or the conductor is coiled. (hyzhou: This is related to the Ponderamotive force. I have some questions here because the time derivative of velocities are all neglected here, which in turn makes the derivation only applicable to a steady system. The obvious problem then is that the generalized Ohm’s law may not hold in a time-dependent system!)

8.5.4 Derivation

Negligible terms:

• \(\frac{d\mathbf{v}_e}{dt}\) in Equation 8.18 and \(\frac{d\mathbf{v}_i}{dt}\) in Equation 8.20, because \(|\mathbf{v}|\ll |\mathbf{v}_\mathrm{th}|\)

• \(\frac{d\mathbf{u}}{dt}\) in Equation 8.18. From Equation 8.22 and \(m_e \ll m_i\), \[ m_en_i\frac{d\mathbf{u}}{dt}=\frac{m_e}{m_i}\frac{n_i}{n_i+n_a}\left(-\nabla p-n_i e\mathbf{v}_{ei}\times\mathbf{B}\right)\ll 1 \]

From Equation 8.18, \[ m_e n_i\frac{d\mathbf{u}}{dt}=-\nabla p_e -n_i e\mathbf{v}_{ei}\times\mathbf{B}-n_ie [\mathbf{E}+(\mathbf{v}_i+\mathbf{u})\times\mathbf{B}]+\mathbf{f}_e^c \]

The above two equations give \[ \begin{aligned} 0=&-\nabla p_e -n_i e\mathbf{v}_{ei}\times\mathbf{B}-n_ie [\mathbf{E}+(\mathbf{v}_i+\mathbf{u})\times\mathbf{B}] \\ & -n_i m_e\mathbf{v}_{ei}\left(\tau_{ei}^{-1}+\tau_{ea}^{-1}\right)-\frac{1}{1-\alpha}n_i m_e\mathbf{v}_i\tau_{ea}^{-1} \end{aligned} \tag{8.23}\]

Equation 8.20, Equation 8.22 and Equation 8.23 are the three main equations we need to derive the generalized Ohm’s law.

If we define the ratios of gyro-period and collision time as \[ \begin{aligned} \kappa_{ei} &\equiv \frac{1}{\Omega_e \tau_{ei}}=\frac{m_e}{eB}\tau_{ei}^{-1} \\ \kappa_{ea} &\equiv \frac{1}{\Omega_e \tau_{ea}}=\frac{m_e}{eB}\tau_{ea}^{-1} \\ \kappa_{ia}&\equiv \frac{1}{2\Omega_i \tau_{ia}}=\frac{m_i}{2eB}\tau_{ia}^{-1} \end{aligned} \] and ion current density as \[ \mathbf{J}_i=n_i e\mathbf{v}_i \] we can write the three main equations above as \[ 0=-\nabla p_e -n_i e(\mathbf{E}+\mathbf{u}\times\mathbf{B}) + \mathbf{J}\times\mathbf{B}-\mathbf{J}_i\times\mathbf{B}+(\kappa_{ei}+\kappa_{ea})B\mathbf{J}-\frac{1}{1-\alpha}\kappa_{ea}B\mathbf{J}_i \tag{8.24}\]

\[ \begin{aligned} m_in_i\frac{d\mathbf{u}}{dt}=&-\nabla p_i+n_i e(\mathbf{E}+\mathbf{u}\times\mathbf{B})+\mathbf{J}_i\times\mathbf{B} \\ &-(\kappa_{ei}-\frac{\alpha}{1-\alpha}\frac{m_e}{m_i}\kappa_{ia})B\mathbf{J}-\frac{1}{1-\alpha}\kappa_{ia}B\mathbf{J}_i \end{aligned} \tag{8.25}\]

\[ \frac{1}{\alpha}m_i n_i\frac{d\mathbf{u}}{dt}=-\nabla p+\mathbf{J}\times\mathbf{B} \tag{8.26}\]

Subtracting Equation 8.25 from Equation 8.26 and Equation 8.24, we get the ion current density as \[ \begin{aligned} \mathbf{J}_i&=\frac{1-\alpha}{(\kappa_{ea}+\kappa_{ia})B}\left[ \alpha\nabla p -\nabla(p_e+p_i)+(1-\alpha)\mathbf{J}\times\mathbf{B}+(\kappa_{ea}+\frac{\alpha}{1-\alpha}\frac{m_e}{m_i}\kappa_{ia})B\mathbf{J}\right] \\ &=\frac{1-\alpha}{(\kappa_{ea}+\kappa_{ia})B}\left[ \frac{\alpha(\alpha-1)}{\alpha+1}\nabla p+(1-\alpha)\mathbf{J}\times\mathbf{B}+(\kappa_{ea}+\frac{\alpha}{1-\alpha}\frac{m_e}{m_i}\kappa_{ia})B\mathbf{J}\right] \end{aligned} \]

Then we can eliminate \(\frac{d\mathbf{u}}{dt}\) and \(\mathbf{J}_i\) from the three main equations and get the generalized Ohm’s law \[ \begin{aligned} \nabla p_e +\frac{\alpha(\alpha-1)}{\alpha+1}r_e \nabla p=&-n_i e(\mathbf{E}+\mathbf{u}\times\mathbf{B}) +\left[1-2r_e(1-\alpha)-\alpha\frac{m_e}{m_i}r_i\right]\mathbf{J}\times\mathbf{B} \\ &+\left[\kappa_{ei}-(1-r_e)\kappa-\frac{\alpha}{1-\alpha}\frac{m_e}{m_i}r_e\kappa_{ia}\right]B\mathbf{J} \\ &+\frac{(1-\alpha)^2}{(\kappa_{ea}+\kappa_{ia})B}\left[ \frac{\alpha}{\alpha+1}\nabla p\times\mathbf{B}-\mathbf{J}\times\mathbf{B}\times\mathbf{B} \right] \end{aligned} \] where \[ r_e=\frac{\kappa_{ea}}{\kappa_{ea}+\kappa_{ia}},\quad r_i=\frac{\kappa_{ia}}{\kappa_{ea}+\kappa_{ia}} \]

8.5.5 Simplification of Generalized Ohm’s Law

1. If the mean free path of protons is much smaller than that of electrons, i.e. the cross section of proton is much smaller than that of electron, \[ \begin{aligned} &\frac{\tau_{ia}}{\tau_{ea}}=\frac{l_{ia}\cdot v_{e,th}}{l_{ea}\cdot v_{i,th}}\propto \Big(\frac{m_i}{m_e}\Big)^{1/2}\cdot\frac{l_{ia}}{l_ea}\ll \Big(\frac{m_i}{m_e}\Big)^{1/2} \\ \Rightarrow &\frac{\kappa_{ea}}{\kappa_{ia}}=2\frac{\Omega_i \tau_{ia}}{\Omega_e\tau_{ea}}\propto 2\frac{m_e}{m_i}\cdot\frac{\tau_{ia}}{\tau_{ea}}\propto2\Big( \frac{m_e}{m_i}\Big)^{1/2}\cdot\frac{l_{ia}}{l_{ea}}\ll 1 \\ \Rightarrow &r_e=\frac{\kappa_{ea}}{\kappa_{ea}+\kappa_{ia}}\ll1\quad r_i=\frac{\kappa_{ia}}{\kappa_{ea}+\kappa_{ia}}\approx 1 \end{aligned} \]

Therefore we have \[ \begin{aligned} \nabla p_e=&-n_i e(\mathbf{E}+\mathbf{u}\times\mathbf{B})+\mathbf{J}\times\mathbf{B}+(\kappa_{ei}+\kappa_{ea})B\mathbf{J}+ \frac{(1-\alpha)^2}{\kappa_{ia}B}\left[ \frac{\alpha}{\alpha+1}\nabla p\times\mathbf{B}-\mathbf{J}\times\mathbf{B}\times\mathbf{B} \right] \\ \mathbf{J}=&\frac{n_i e}{(\kappa_{ei}+\kappa_{ea})B}\left(\mathbf{E}+\mathbf{u}\times\mathbf{B}\right)+ \frac{1}{(\kappa_{ei}+\kappa_{ea})B}\left( \nabla p_e-\mathbf{J}\times\mathbf{B}\right) \\ &-\frac{(1-\alpha)^2}{\kappa_{ia}(\kappa_{ei}+\kappa_{ea})B}\left[ \frac{\alpha}{\alpha+1}\nabla p\times\mathbf{B}-\mathbf{J}\times\mathbf{B}\times\mathbf{B} \right] \end{aligned} \]

2. If EM field is dominant, i.e. EM force \(\gg\) thermal pressure gradient, \(|\mathbf{J}\times\mathbf{B}|\gg |\nabla p|\sim |\nabla p_e|\), then \[ \mathbf{J}=\sigma(\mathbf{E}+\mathbf{u}\times\mathbf{B})-\frac{\sigma}{n_i e}\mathbf{J}\times\mathbf{B}+ \frac{\sigma}{n_i e}\frac{2(1-\alpha)^2 \Omega_i\tau_{ia}}{B}\mathbf{J}\times\mathbf{B}\times\mathbf{B} \] where \(\sigma\) is the electrical conductivity, \[ \sigma=\frac{n_i e^2}{m_e(\tau_{ei}^{-1}+\tau_{ea}^{-1})} \]

3. For a fully ionized plasma, \(\alpha=1\ \Rightarrow\ \kappa_{ea}=\kappa_{ia}=0\), we have \[ \mathbf{J}=\sigma(\mathbf{E}+\mathbf{u}\times\mathbf{B})+\frac{\sigma}{n_i e}\left(\nabla p_e -\mathbf{J}\times\mathbf{B}\right) \] where \[ \sigma=\frac{n_i e^2\tau_{ei}}{m_e} \]

Note that the electron pressure gradient term here shows an amazing effect: it can generate EM field out of nothing! This is one theory in astrophysics which describes the origin of EM field.

4. For a fully ionized plasma, if EM field is dominant, we have \[ \begin{aligned} \mathbf{J}=\sigma(\mathbf{E}+\mathbf{u}\times\mathbf{B}) -\frac{\sigma}{n_i e}\mathbf{J}\times\mathbf{B}= \sigma(\mathbf{E}+\mathbf{u}\times\mathbf{B}) -\Omega_e \tau_{ei}\mathbf{J}\times\mathbf{b} \end{aligned} \]

5. If there is weak magnetic field and dense plasma, i.e. \(n_i \gg 1,\ |\mathbf{B}|\ll 1\), \[ \mathbf{J}=\sigma(\mathbf{E}+\mathbf{u}\times\mathbf{B}) \]

8.5.6 Discussions

In Section 8.6.1, we will see that a relatively simple form can be obtained by assuming that the electron inertia is much smaller than the Lorentz force, which is a good approximation when considering only the perpendicular velocity to the magnetic field.

\[ \mathbf{E} = -\mathbf{u}\times\mathbf{B} + \eta\mathbf{J} + \frac{1}{en}\mathbf{J}\times\mathbf{B}-\frac{1}{en}\nabla\cdot\overleftrightarrow{P}_e +\frac{m_e}{ne^2}\Big[ \frac{\partial \mathbf{J}}{\partial t}+\nabla\cdot(\mathbf{J}\mathbf{u}+\mathbf{u}\mathbf{J}) \Big] \tag{8.27}\]

where the first term on the right-hand side is the convection term, the second term is the resistivity term (conductive term), the third term is called the Hall term, the fourth term is the electron pressure term, and the fifth term is called the electron inertia term, since it is proportional to the mass of electrons.

Note that both \(\mathbf{u}\) and \(\mathbf{J}\) are the first-order moments, with \(\mathbf{u}\) being the (weighted) sum of the first-order moment of electrons and ions while \(\mathbf{J}\) being the difference between them. The generalized Ohm’s law is actually the difference between the electrons’ and ions’ first-order moment equations. The generalized Ohm’s law is an equation that governs the time evolution of \(\mathbf{J}\). Also note that Ampère’s law, with the displacement current retained, is an equation governing the time evolution of \(\mathbf{E}\). However, in the approximation of the resistive MHD, the time derivative terms \(\partial \mathbf{E}/\partial t\) and \(\partial \mathbf{J} /\partial t\) are ignored in Ampère’s law and Ohm’s law, respectively. In this approximation, Ohm’s law is directly solved to determine \(\mathbf{E}\) and Ampère’s law is directly solved to determine \(\mathbf{J}\). [Introduction to plasma physics: with space and laboratory applications, D. A. Gurnett and A. Bhattacharjee.]

\[ \underbrace{\mathbf{E}}_{a} + \underbrace{\mathbf{u}\times\mathbf{B}}_{b} - \underbrace{\frac{1}{\sigma}\mathbf{J}}_{c} - \underbrace{\frac{\mathbf{J}\times\mathbf{B}}{en_e}}_{d} + \underbrace{\frac{1}{n_e e}\nabla(n_e k_B T_e)}_{e} = 0 \]

Denote each term above with a to e. From the single-fluid MHD momentum equation, let \(\omega\) be the oscillation frequency of the perturbed velocity, and variable in scalar form be the characteristic magnitude of that quantity, \[ \begin{aligned} \rho\big( \frac{\partial \mathbf{u}}{\partial t}+\mathbf{u}\cdot\nabla\mathbf{u}\big) = \mathbf{J}\times\mathbf{B} - \nabla p +\rho\mathbf{g} \\ \Rightarrow \omega\rho U\sim J B \end{aligned} \tag{8.28}\]

Using and the assumptions of MHD, the relations between each term in generalized Ohm’s law are \[ \begin{aligned} \frac{b}{d} &= \frac{UB}{JB/en_e} = \frac{en_eUB}{\omega \rho U} = \frac{en_e B}{\omega n_e m_i} = \frac{\Omega_i}{\omega} \gg 1 \\ \frac{c}{d} &= \frac{J/\sigma}{JB/en_e} = \frac{en_e}{B\sigma} = \frac{en_e}{Be^2 n_e/\nu_{ei} m_e} = \frac{\nu_{ei}}{\Omega_e} \ll 1 \\ \frac{b}{d} &= \frac{U}{J/en_e} \ll 1,\ \text{if } \mathbf{J} \text{ is carried by } \mathbf{u}_e \end{aligned} \] which generates direct contradiction with the MHD assumption that \(\mathbf{u}\ll \mathbf{u}_e\) if currents are mostly carried by electrons. (WHAT ABOUT THE PRESSURE TERM SCALING? The gradient implies that it is related to the system size; the pressure implies that it is also related to thermal motion?)

There is an interesting point about the electron pressure term. If we assume an isotropic plasma with \(n_i=n_e\) and an adiabatic process \(P/n^\gamma=C\) where \(C\) is a constant and \(\gamma\) is the adiabatic index, we have \[ \begin{aligned} \mathbf{E}_{\nabla P_e} &= \frac{\nabla P_e}{n} = \frac{\nabla P_i}{n} = \frac{\nabla(Cn^\gamma)}{n} \\ &= C\gamma n^{\gamma-2}\nabla n = C\frac{\gamma}{\gamma-1}\nabla n^{\gamma-1} \end{aligned} \tag{8.29}\]

Equation 8.29 indicates that the associated electric field is a potential field that only relates to density.

8.5.7 Ohm’s Law in a Multi-Ion System

\[ \begin{aligned} \mathbf{E} &= -\sum_i\mathbf{u}_i\times\mathbf{B} + \eta\mathbf{J} + \frac{1}{en_e}\mathbf{J}\times\mathbf{B}-\frac{1}{en_e}\nabla\cdot\overleftrightarrow{P}_e \\ &+\frac{m_e}{n_e e^2}\Big[ \frac{\partial \mathbf{J}}{\partial t}+\nabla\cdot(\mathbf{J}\mathbf{u}+\mathbf{u}\mathbf{J}) \Big] \end{aligned} \tag{8.30}\]

where the quasi-neutrality holds \[ en_e = \sum_i q_i n_i \]

8.6 Magnetohydrodynamics

Particle motion in the two-fluid system was described by the individual species’ mean velocities \(\mathbf{u}_e, \mathbf{u}_i\) and by the pressures \(\overleftrightarrow{P}_e, \overleftrightarrow{P}_i\), which provide an accounting for the mean square of the random deviation of the velocity from its average value. Magnetohydrodynamics is an alternate description of the plasma where, instead of using \(\mathbf{u}_e,\mathbf{u}_i\) to describe mean motion, two new velocity variables that are a linear combination of \(\mathbf{u}_e,\mathbf{u}_i\) are used. As will be seen below, this means a slightly different definition for pressure must also be used.

The new velocity-like variables are

1. the current density \[ \mathbf{J} = \sum_s n_s q_s\mathbf{u}_s \] which is essentially the relative velocity between ions and electrons, and

2. the center-of-mass velocity \[ \mathbf{u} = \frac{1}{\rho}\sum_s m_s n_s \mathbf{u}_s \] where \[ \rho = \sum_s m_s n_s \] is the total mass density.

Magnetohydrodynamics is primarily concerned with low-frequency, long-wavelength, magnetic behavior of the plasma. There are three basic assumptions in MHD: \[ \begin{aligned} \epsilon\omega/4\pi\sigma& \ll 1, \\ (v/c)^2&\ll 1,\\ \lambda/ L& \ll 1, \end{aligned} \] where \(\omega\) is the plasma frequency, \(v\) is the plasma bulk speed, \(\lambda\) is the average free distance, and \(L\) is the system characteristic length.

Standard orderings of ideal MHD can also be written as \[ \epsilon \sim \omega/\Omega_c \sim k\rho \] where the plasma varies on frequency scales \(\omega\) small compared to the gyrofrequency \(\Omega_c\), and varies on spatial scales \(1/k\) long compared to the gyroradius \(\rho\). (Here, \(\epsilon\) means a “small” value.) Thus it covers phenomenon related to compressional and shear Alfvén waves and instabilities, ion acoustic waves, and ion and electron kinetic effects such as Landau damping. However, it does not include drift-waves or other micro-instabilities because they result from finite Larmor radius (FLR) effects which vanish in the usual MHD ordering.

Single fluid MHD is somehow inconsistent since there is only one velocity. The definition of current using velocity cannot be applied, and the current can only be given by Ampère’s law (without the displacement current), \(\mathbf{J}=\nabla\times\mathbf{B}/\mu_0\).

Multiplying Equation 8.5 by \(m_s\) and summing over species gives the MHD continuity equation \[ \frac{\partial\rho}{\partial t} + \nabla\cdot(\rho\mathbf{u}) = 0 \tag{8.31}\]

To obtain an MHD equation of motion we take the first moment of the Vlasov equation and multiply the single species momentum equation by \(m_s\) and sum over species, \[ \frac{\partial}{\partial t}\sum_s m_s \int\mathbf{v}f_s\mathrm{d}\mathbf{v} + \frac{\partial}{\partial\mathbf{x}}\cdot\sum_s\int m_s \mathbf{v}\mathbf{v}f_s\mathrm{d}\mathbf{v} + \sum_s q_s \int\mathbf{v}\frac{\partial}{\partial\mathbf{v}}\cdot [(\mathbf{E}+\mathbf{v}\times\mathbf{B})f_s] = 0 \tag{8.32}\]

The right-hand side is zero since \(\mathbf{R}_{ei} + \mathbf{R}_{ie}=0\), i.e., the total plasma cannot exert drag on itself. We now define random velocities relative to \(\mathbf{u}\) (rather than to \(\mathbf{u}_s\) as was the case for the two-fluid equations) so that the second term can be written as \[ \sum_s\int m_s \mathbf{v}\mathbf{v}f_s\mathrm{d}\mathbf{v} = \sum_s\int m_s (\mathbf{v}^\prime + \mathbf{u})(\mathbf{v}^\prime + \mathbf{u})f_s\mathrm{d}\mathbf{v} = \sum_s\int m_s \mathbf{v}^\prime\mathbf{v}^\prime f_s\mathrm{d}\mathbf{v} + \rho\mathbf{u}\mathbf{u} \tag{8.33}\] where \(\sum_s\int m_s \mathbf{v}^\prime f\mathrm{d}\mathbf{v} = 0\) has been used to eliminate terms linear in \(\mathbf{v}^\prime\). The MHD pressure tensor is now defined in terms of the random velocities relative to \(\mathbf{u}\) and is given by \[ \overleftrightarrow{P}^{\mathrm{MHD}} = \sum_s \int m_s \mathbf{v}^\prime\mathbf{v}^\prime f_s\mathrm{d}\mathbf{v} \tag{8.34}\]

We insert Equation 8.33 and Equation 8.34 in Equation 8.32, integrate by parts on the acceleration term, and perform the summation over species to obtain the MHD equation of motion \[ \frac{\partial(\rho\mathbf{u})}{\partial t} + \nabla\cdot(\rho\mathbf{u}\mathbf{u}) = \left( \sum_s n_s q_s \right)\mathbf{E} + \mathbf{J}\times\mathbf{B} - \nabla\cdot\overleftrightarrow{P}^{\mathrm{MHD}} \tag{8.35}\]

MHD is typically used to describe phenomena having spatial scales large enough for the plasma to be essentially neutral, i.e., \(\sum_s n_s q_s = 0\) so that the \(\mathbf{E}\) term can be dropped. Just as in the two-fluid situation, the left-hand side of Equation 8.35 contains a factor times the MHD continuity equation, since the left-hand side can be expanded as \[ \begin{aligned} \frac{\partial(\rho\mathbf{u})}{\partial t} + \nabla\cdot(\rho\mathbf{u}\mathbf{u}) &= \left[ \cancel{\frac{\partial\rho}{\partial t} + \nabla\cdot(\rho\mathbf{u})} \right]\mathbf{u} + \rho\frac{\partial\mathbf{u}}{\partial t} + \rho\mathbf{u}\cdot\nabla\mathbf{u} \\ &= \rho \left( \frac{\partial\mathbf{u}}{\partial t} + \mathbf{u}\cdot\nabla\mathbf{u} \right) \end{aligned} \]

This leads to the standard form for the MHD equation of motion, \[ \rho\frac{\mathrm{D}\mathbf{u}}{\mathrm{D}t} = \mathbf{J}\times\mathbf{B} - \nabla\cdot\overleftrightarrow{P}^{\mathrm{MHD}} \tag{8.36}\] where \[ \frac{\mathrm{D}}{\mathrm{D}t} \equiv \frac{\partial}{\partial t} + \mathbf{u}\cdot\nabla \] is the convective derivative defined using the center-of-mass velocity. This notation is used to emphasize the difference from single species treatment. Scalar approximations of the MHD pressure tensor will be postponed until after discussing the MHD Ohm’s law and its implications.

8.6.1 MHD Ohm’s law

Equation 8.36 provides one equation relating \(\mathbf{J}\) and \(\mathbf{u}\); let us now find the other one. In order to do this, consider the two-fluid electron equation of motion, \[ m_e\frac{\mathrm{d}\mathbf{u}_e}{\mathrm{d}t} = -e(\mathbf{E}+\mathbf{u}_e\times\mathbf{B}) - \frac{1}{n_e}\nabla(n_e k_B T_e) - \nu_{ei}m_e(\mathbf{u}_e - \mathbf{u}_i) \tag{8.37}\]

In MHD we are interested in low-frequency phenomena with large spatial scales. If the characteristic time scale of the phenomenon is long compared to the electron cyclotron motion, then the electron inertia term \(m_e \mathrm{d}\mathbf{u}_e/\mathrm{d}t\) can be dropped since it is small compared to the magnetic force term \(-e(\mathbf{u}_e\times\mathbf{B})\). This assumption is reasonable for velocities perpendicular to \(\mathbf{B}\), but can be a poor approximation for the velocity component parallel to \(\mathbf{B}\), since parallel velocities do not provide a magnetic force. Since \(\mathbf{u}_e - \mathbf{u}_i = -\mathbf{J}/n_e e\) and \(\mathbf{u}_i\simeq \mathbf{u}\), Equation 8.37 reduces to the generalized Ohm’s law \[ \mathbf{E} + \mathbf{u}\times\mathbf{B} - \frac{1}{n_e e}\mathbf{J}\times\mathbf{B} + \frac{1}{n_e e}\nabla(n_e k_B T_e) = \eta\mathbf{J} \tag{8.38}\]

(What if we consider ion species with charge larger than 1?) The term \(-\mathbf{J}\times\mathbf{B}/n_e e\) on the left-hand side is called the Hall term and can be neglected in either of the following two cases:

1. The pressure term in the MHD equation of motion, Equation 8.36 is negligible compared to the other two terms, which therefore must balance, giving \[ |\mathbf{J}|\sim \omega\rho|\mathbf{u}|/|\mathbf{B}| \]

Here \(\omega\sim\mathrm{D}/\mathrm{D}t\) is the characteristic frequency of the phenomenon. In this case comparison of the Hall term with the U×B term shows that the Hall term is small by a factor \(\sim \omega/\omega_{ci}\), where \(\omega_{ci}=q_i B/m_i\) is the ion cyclotron frequency. Thus dropping the Hall term is justified for phenomena having characteristic frequencies small compared to \(\omega_{ci}\).

2. The electron-ion collision frequency is large compared to the electron cyclotron frequency \(\omega_{ce}=q_e B /m_e\), in which case the Hall term may be dropped since it is small by a factor \(\omega_{ce}/\nu_{ei}\) compared to the right-hand side resistive term \(\eta\mathbf{J}=(m_e\nu_{ei}/n_e e^2)\mathbf{J}\).

8.6.2 Ideal MHD

If the Hall term is dropped, the system is called ideal MHD. Typically, Equation 8.38 will not be used directly; instead its curl will be used to provide the induction equation \[ -\frac{\partial\mathbf{B}}{\partial t} + \nabla\times(\mathbf{u}\times\mathbf{B}) - \frac{1}{n_e e}\nabla n_e\times\nabla k_B T_e = \nabla\times\left( \frac{\eta}{\mu_0}\nabla\times\mathbf{B} \right) \]

Usually the density gradient is parallel to the temperature gradient so that the thermal electromotive force term \((n_e e)^{-1}\nabla n_e\times\nabla k_B T_e\) can be dropped, or the thermal term is often simply ignored, in which case the induction equation reduces to \[ -\frac{\partial\mathbf{B}}{\partial t} + \nabla\times(\mathbf{u}\times\mathbf{B}) = \nabla\times\left( \frac{\eta}{\mu_0}\nabla\times\mathbf{B} \right) \tag{8.39}\] or written in the form of Ohm’s law \[ \mathbf{E} + \mathbf{u}\times\mathbf{B} = \eta\mathbf{J} \tag{8.40}\]

If the resistive term \(\eta\mathbf{J}\) is so small as to be negligible compared to the other terms in Equation 8.40, then the plasma is said to be ideal or perfectly conducting. From the Lorentz transformation of electromagnetic theory we realize that \[ \mathbf{E}+\mathbf{u}\times\mathbf{B} = \mathbf{E}^\prime \] where \(\mathbf{E}^\prime\) is the electric field observed in the frame moving with velocity \(\mathbf{u}\). This implies that the magnetic flux in ideal plasmas is time-invariant in the frame moving with velocity \(\mathbf{u}\), because otherwise Faraday’s law would imply the existence of an electric field in the moving frame. The frozen-in flux concept is the essential “bed-rock” concept underlying ideal MHD. While this concept is often an excellent approximation, it must be kept in mind that the concept becomes invalid in situations when any one of the electron inertia, electron pressure, or Hall terms becomes important and leads to different, more complex behavior.

The frozen-in flux concept is frequently expressed in a slightly different form as a frozen-in field concept, i.e., it is often said that magnetic field lines are frozen into the plasma so that the plasma and magnetic field lines move together as an ensemble. While this point of view can be quite intuitive and useful, it contains some ambiguity because an ideal plasma can actually move across magnetic field lines in certain situations. These situations are such that magnetic flux is preserved within the plasma even though it is moving across field lines.

A formal proof of the frozen-in flux property will now be established by direct calculation of the rate of change of the magnetic flux through a surface \(S(t)\) bounded by a material line \(C(t)\), i.e., a closed contour that moves with the plasma. This magnetic flux is \[ \Phi(t) = \int_{S(t)} \mathbf{B}(\mathbf{x},t)\cdot\mathrm{d}\mathbf{s} \] and the flux changes with respect to time due to either (i) the explicit time dependence of \(\mathbf{B}(t)\) or (ii) changes in the surface \(S(t)\) resulting from plasma motion. The rate of change of flux is thus \[ \frac{\mathrm{D}\Phi}{\mathrm{D}t} = \lim_{\delta t\rightarrow 0}\left( \frac{ \int_{S(t+\delta t)} \mathbf{B}(\mathbf{x},t+\delta t)\cdot\mathrm{d}\mathbf{s} - \int_{S(t)} \mathbf{B}(\mathbf{x},t)\cdot\mathrm{d}\mathbf{s} }{\delta t} \right) \]

The displacement of a segment \(\mathrm{d}\mathbf{l}\) of the bounding contour \(C\) is \(\mathbf{u}\delta t\), where \(\mathbf{u}\) is the velocity of this segment. The incremental change in surface area due to this displacement is \(\Delta S = \mathbf{u}\delta t\times\mathrm{d}\mathbf{l}\). The rate of change of flux can thus be expressed as \[ \begin{aligned} \frac{\mathrm{D}\Phi}{\mathrm{D}t} &= \lim_{\delta t\rightarrow 0}\left( \frac{ \int_{S(t+\delta t)} \left(\mathbf{B}+\delta t\frac{\partial\mathbf{B}}{\partial t}\right)\cdot\mathrm{d}\mathbf{s} - \int_{S(t)} \mathbf{B}\cdot\mathrm{d}\mathbf{s} }{\delta t} \right) \\ &= \lim_{\delta t\rightarrow 0}\left( \frac{ \int_{S(t)} \left(\mathbf{B}+\delta t\frac{\partial\mathbf{B}}{\partial t}\right)\cdot\mathrm{d}\mathbf{s} + \oint_C\mathbf{B}\cdot\mathbf{u}\delta t\times\mathrm{d}\mathbf{l} - \int_{S(t)} \mathbf{B}\cdot\mathrm{d}\mathbf{s} }{\delta t} \right) \\ &= \int_{S(t)}\frac{\partial\mathbf{B}}{\partial t}\cdot\mathrm{d}\mathbf{s} + \oint_C\mathbf{B}\cdot\mathbf{u}\times\mathrm{d}\mathbf{l} \\ &= \int_{S(t)}\frac{\partial\mathbf{B}}{\partial t}\cdot\mathrm{d}\mathbf{s} + \oint_C\mathrm{d}\mathbf{l}\cdot(\mathbf{B}\cdot\mathbf{u}) \\ &= \int_{S(t)}\left[ \frac{\partial\mathbf{B}}{\partial t} + \nabla\times(\mathbf{B}\times\mathbf{u}) \right]\cdot\mathrm{d}\mathbf{S} \end{aligned} \]

Thus, if \[ \frac{\partial\mathbf{B}}{\partial t} + \nabla\times(\mathbf{B}\times\mathbf{u}) = 0 \tag{8.41}\] then \[ \frac{\mathrm{D}\Phi}{\mathrm{D}t} = 0 \] so that the magnetic flux linked by any closed material line is constant. Therefore, magnetic flux is frozen into an ideal plasma because Equation 8.39 reduces to Equation 8.41 if \(\eta = 0\). Equation 8.41 is called the ideal MHD induction equation.

8.6.2.1 Equation of state

Double adiabatic laws

A procedure analogous to that which led to Equation 8.11 or Equation 12.28 gives the MHD adiabatic relation \[ \frac{p^{\mathrm{MHD}}}{\rho^\gamma} = \mathrm{const}. \] where again \(\gamma=(N+2)/N\) and \(N\) is the number of dimensions of the system. It was shown in the previous section that magnetic flux is conserved in the plasma frame. This means that, as shown in Figure 8.1, a tube of plasma initially occupying the same volume as a magnetic flux tube is constrained to evolve in such a way that \(\int\mathbf{B}\cdot\mathrm{d}\mathbf{s}\) stays constant over the plasma tube cross-section. For a flux tube of infinitesimal cross-section, the magnetic field is approximately uniform over the cross-section and we may write this as \(BA = \mathrm{const}\), where \(A\) is the cross-sectional area.

Figure 8.1: Magnetic flux tube with flux \(\Phi=BA\).

Let us define two temperatures for this magnetized plasma, namely \(T_\perp\), the temperature corresponding to motions perpendicular to the magnetic field, and \(T_\parallel\), the temperature corresponding to motions parallel to the magnetic field. If for some reason (e.g., anisotropic heating or compression) the temperature develops an anisotropy such that \(T_\perp \neq T_\parallel\) and if collisions are infrequent, this anisotropy will persist for a long time, since collisions are the means by which the two temperatures equilibrate. Thus, rather than assuming that the MHD pressure is fully isotropic, we consider the less restrictive situation where the MHD pressure tensor is given by \[ \overleftrightarrow{P}^{\mathrm{MHD}} = \begin{bmatrix} p_\perp & 0 & 0 \\ 0 & p_\perp & 0 \\ 0 & 0 & p_\parallel \end{bmatrix} = p_\perp\mathbf{I} + (p_\parallel - p_\perp)\hat{b}\hat{b} \]

The first two coordinates (\(x, y\)-like) in the above matrix refer to the directions perpendicular to the local magnetic field \(\mathbf{B}\) and the third coordinate (\(z\)-like) refers to the direction parallel to \(\mathbf{B}\). The tensor expression on the right-hand side is equivalent (here \(\mathbf{I}\) is the unit tensor) but allows for arbitrary, curvilinear geometry. We now develop separate adiabatic relations for the perpendicular and parallel directions:

• Parallel direction: here the number of dimensions is \(N = 1\) so that \(\gamma = 3\) and so the adiabatic law gives \[ \frac{p_\parallel^{1D}}{\rho_{1D}^3}=\mathrm{const.} \tag{8.42}\] where \(\rho_{1D}\) is the one-dimensional mass density; i.e., \(\rho_{1D}\sim 1/L\), where \(L\) is the length along the flux tube in Figure 8.1. The three-dimensional mass density, which has been used implicitly until now, has the proportionality \(\rho\sim 1/LA\), where \(A\) is the cross-section of the flux tube; similarly the three-dimensional pressure has the proportionality \(p_\parallel \sim \rho T_\parallel\). However, we must be careful to realize that \(p_\parallel^{1D}\sim\rho_{1D}T_\parallel\) so, using \(BA=\mathrm{const.}\), Equation 8.42 can be recast as \[ \mathrm{const.} = \frac{p_\parallel^{1D}}{\rho_{1D}^3} \sim \frac{\rho_{1D}T_\parallel}{\rho_{1D}^3} \sim T_\parallel \rho_{1D}^2 \sim \left( \frac{1}{LA} \right) T_\parallel (LA)^3 B^2 \] or \[ \frac{p_\parallel B^2}{\rho^3} = \mathrm{const.} \tag{8.43}\]

• Perpendicular direction: here the number of dimensions is \(N = 2\) so that \(\gamma = 2\) and the adiabatic law gives \[ \frac{p_\perp^{2D}}{\rho_{2D}^2}=\mathrm{const.} \tag{8.44}\] where \(p_\perp^{2D}\) is the 2-D perpendicular pressure, and has dimensions of energy per unit area, while \(\rho_{2D}\) is the 2-D mass density and has dimensions of mass per unit area. Thus, \(\rho_{2D}\sim 1/A\) so \(p_\perp^{2D}\sim\rho_{2D}T_\perp \sim T_\perp/A\), in which case Equation 8.44 can be recast as \[ \mathrm{const.} = \frac{p_\parallel B^2}{\rho^3} = T_\perp A \sim \left( \frac{1}{LA} \right) T_\perp \frac{LA}{B} \] or \[ \frac{p_\perp}{\rho B} = \mathrm{const.} \tag{8.45}\]

Equation 8.43 and Equation 8.45 are called the double adiabatic or CGL laws after Chew, Goldberger, and Low (Chew, Goldberger, and Low 1956) who first developed these laws.

Single adiabatic law

If collisions are sufficiently frequent to equilibrate the perpendicular and parallel temperatures, then the pressure tensor becomes fully isotropic and the dimensionality of the system is \(N = 3\) so that \(\gamma = 5/3\). There is now just one pressure and one temperature and the adiabatic relation becomes \[ \frac{p}{\rho^{5/3}} = \mathrm{const.} \tag{8.46}\]

8.6.2.2 MHD approximations for Maxwell’s equations

The various assumptions contained in MHD lead to a simplifying approximation of Maxwell’s equations. In particular, the assumption of charge neutrality in MHD makes Poisson’s equation superfluous because Poisson’s equation prescribes the relationship between non-neutrality and the electrostatic component of the electric field. The assumption of charge neutrality also has implications for the current density. To see this, the two-fluid continuity equation is multiplied by \(q_s\) and then summed over species to obtain the charge conservation equation \[ \frac{\partial}{\partial t}\left( \sum_s n_s q_s \right) + \nabla\cdot\mathbf{J} = 0 \tag{8.47}\]

Thus, charge neutrality implies \[ \nabla\cdot\mathbf{J} = 0 \]

Let us now consider Ampère’s law \[ \nabla\times\mathbf{B} = \mu_0\mathbf{J} + \mu_0\epsilon_0\dot{\mathbf{E}} \]

Taking the divergence gives \[ \nabla\cdot \] which is equivalent to Equation 8.47 if Poisson’s equation is invoked.

Finally, MHD is restricted to phenomena having characteristic velocities \(V_{\mathrm{ph}}\) slow compared to the speed of light in vacuum, \(c=(\mu_0 \epsilon_0)^{-1/2}\). Again, \(t_{\mathrm{char}}\) is assumed to represent the characteristic time scale for a given phenomenon and \(l_{\mathrm{char}}\) is assumed to represent the corresponding characteristic length scale so that \(V_{\mathrm{ph}}\sim l_{\mathrm{char}}/t_{\mathrm{char}}\). Faraday’s equation gives the scaling \[ \nabla\times\mathbf{E} = -\dot{\mathbf{B}} \rightarrow E \sim B l_{\mathrm{char}}/t_{\mathrm{char}} \] On comparing the magnitude of the displacement current term in Ampère’s law to the left-hand side it is seen that \[ \frac{\mu_0\epsilon_0|\dot{\mathbf{E}}|}{|\nabla\times\mathbf{B}|} \sim \frac{c^{-2} E/t_{\mathrm{char}}}{B/l_{\mathrm{char}}} \sim \left( \frac{V_{\mathrm{ph}}}{c} \right)^2 \]

Thus, if \(V_{\mathrm{ph}} \ll c\) the displacement current term can be dropped from Amperè’s law resulting in the so-called “pre-Maxwell” form (i.e. Darwin approximation) \[ \nabla\times\mathbf{B} = \mu_0 \mathbf{J} \tag{8.48}\]

The divergence of Equation 8.48 gives \(\nabla\cdot\mathbf{J}=0\) so it is unnecessary to specify it separately.

8.6.2.3 Complete equation set

In the most common sense, ideal MHD involves equations of compressible, adiabatic and inviscid fluid: \[ \begin{aligned} \frac{\partial \rho}{\partial t} + \nabla\cdot(\rho\mathbf{u}) = 0 \\ \frac{\partial \rho\mathbf{u}}{\partial t} + \nabla\cdot\left( \rho\mathbf{u}\mathbf{u} - \frac{1}{\mu_0}\mathbf{B}\mathbf{B} + p^\ast \right) = 0 \\ \frac{\partial \mathcal{E}}{\partial t} + \nabla\cdot\left[ (\mathcal{E} + p^\ast)\mathbf{u} - \frac{1}{\mu_0}\mathbf{B}(\mathbf{B}\cdot\mathbf{u}) \right] = 0 \\ \frac{\partial \mathbf{B}}{\partial t} - \nabla\times(\mathbf{u}\times\mathbf{B}) = 0 \\ p^\ast = p + \frac{\mathbf{B}\cdot\mathbf{B}}{2\mu_0} \\ \mathcal{E} = \frac{p}{\gamma - 1} + \frac{\rho \mathbf{u}\cdot\mathbf{u}}{2} + \frac{\mathbf{B}\cdot\mathbf{B}}{2\mu_0} \end{aligned} \tag{8.49}\] where \(\rho\) is the mass density, \(\mathbf{u}\) the velocity, \(\mathcal{E}\) the total energy density, \(\mathbf{B}\) the magnetic field, \(p\) the thermal pressure, and \(\gamma\) the adiabatic index (ratio of specific heats). Microscopic dissipations of any kind (viscosity, resistivity, or conduction) are not included in ideal MHD. Note that there is only one constant \(\mu_0\) appeared in Equation 8.49. The introduction of temperature comes with a new constant \(R\equiv 2k_B/m\): \[ p = \rho\, R\, T \]

Many often in simulations, the dimensionless units are applied. Equation 8.49 under the dimensionless units can be written as \[ \begin{aligned} \frac{\partial \rho_\ast}{\partial t_\ast} + \nabla_\ast\cdot(\rho_\ast\mathbf{u}_\ast) = 0 \\ \frac{\partial \rho_\ast\mathbf{u}_\ast}{\partial t_\ast} + \nabla_\ast\cdot\left( \rho_\ast\mathbf{u}_\ast\mathbf{u}_\ast - \mathbf{B}_\ast\mathbf{B}_\ast + p_\ast^\ast \right) = 0 \\ \frac{\partial \mathcal{E}_\ast}{\partial t_\ast} + \nabla_\ast\cdot\left[ (\mathcal{E}_\ast + p_\ast^\ast)\mathbf{u}_\ast - \mathbf{B}_\ast(\mathbf{B}_\ast\cdot\mathbf{u}_\ast) \right] = 0 \\ \frac{\partial \mathbf{B}_\ast}{\partial t_\ast} - \nabla_\ast\times(\mathbf{u}_\ast\times\mathbf{B}_\ast) = 0 \\ p_\ast^\ast = p_\ast + \frac{\mathbf{B}_\ast\cdot\mathbf{B}_\ast}{2} \\ \mathcal{E}_\ast = \frac{p_\ast}{\gamma - 1} + \frac{\rho_\ast \mathbf{u}_\ast\cdot\mathbf{u}_\ast}{2} + \frac{\mathbf{B}_\ast\cdot\mathbf{B}_\ast}{2} \end{aligned} \tag{8.50}\] where the subscripts \(\ast\) denotes the variables in dimensionless units, and the adiabatic equation of state is given as \[ p_\ast = \rho_\ast T_\ast \]

In more theoretical derivations, for simplicity the energy equation is substituted with either single adiabatic law Equation 8.46 or double adiabatic law Equation 8.43 and Equation 8.45.³ If we drop the isotropic assumption and instead use separate pressure equations for the parallel and perpendicular components as in Equation 8.14, we get the anisotropic MHD equations:⁴ \[ \begin{aligned} \frac{\partial \rho_s}{\partial t} + \nabla\cdot(\rho_s\mathbf{u}_s) = 0 \\ \frac{\partial (\rho_s\mathbf{u}_s)}{\partial t} + \nabla\cdot\left[ \rho_s\mathbf{u}_s\mathbf{u}_s + p_{s\perp}\pmb{I} + (p_{s\parallel} - p_{s\perp})\hat{b}\hat{b} - \frac{1}{\mu_0}\left( \mathbf{B}\mathbf{B} - \frac{B^2}{2}\pmb{I} \right) \right] = 0 \\ \frac{\partial\mathbf{B}}{\partial t} - \nabla\times\left( \mathbf{u}_s \times\mathbf{B} \right) = 0 \\ \frac{\partial p_{s\parallel}}{\partial t} + \nabla\cdot(p_{s\parallel}\mathbf{u}_s) + 2p_{s\parallel}\hat{b}\cdot(\hat{b}\cdot\nabla)\mathbf{u}_s = \frac{\delta p_{s\parallel}}{\delta t} \\ \frac{\partial p_{s}}{\partial t} + \nabla\cdot(p_{s}\mathbf{u}_s) + \frac{2}{3}p_{s\perp}(\nabla\cdot\mathbf{u}_s) + \frac{2}{3}(p_\parallel - p_\perp)\hat{b}\cdot(\hat{b}\cdot\nabla)\mathbf{u}_s = 0 \end{aligned} \tag{8.51}\]

For capturing jump conditions across a discontinuity, the conservative energy equation is needed to substitute the total pressure equation: \[ \frac{\partial e}{\partial t} + \nabla\cdot\left[ \mathbf{u}\left( e+p_\perp + \frac{\mathbf{B}^2}{2\mu_0} + \mathbf{u}\cdot\left( (p_\parallel - p_\perp)\hat{b}\hat{b} - \frac{\mathbf{B}\mathbf{B}}{2\mu_0} \right) \right) \right] \tag{8.52}\] with total energy density \[ e = \frac{\rho \mathbf{u}^2}{2} + \frac{\mathbf{B}^2}{2\mu_0} + \frac{3}{2}p \tag{8.53}\]

The plasma instabilities associated with anisotropic pressure limit, i.e. firehose, mirror, and ion cyclotron instability may contain kinetic effects that cannot be fully described by MHD (Chapter 14). Moreover, the grid resolution in normal global MHD simulations may not be fine enough to resolve even hydromagnetic instabilities.

• Criterion for the firehose instability \[ \frac{p_\parallel}{p_\perp} > 1 + \frac{\mathbf{B}^2}{\mu_0 p_\perp} \tag{8.54}\]

• Criterion for the mirror instability \[ \frac{p_\perp}{p_\parallel} > 1 + \frac{\mathbf{B}^2}{2\mu_0 p_\perp} \tag{8.55}\]

• Criterion for the ion cyclotron instability \[ \frac{p_\perp}{p_\parallel} > 1 + C_1\left(\frac{\mathbf{B}^2}{2\mu_0 p_\parallel}\right)^{C_2} \tag{8.56}\] where from observations in the magnetosphere [Anderson+ 1996; Gary+ 1995] we can take the average values \(C_1 = 0.3\) and \(C_2 = 0.5\).

The effect of these instabilities is represented by the collision term of the parallel pressure on the RHS of the fourth equation in Equation 8.51: \[ \frac{\delta p_\parallel}{\delta t} = \frac{\bar{p}_\parallel - p_\parallel}{\tau} \tag{8.57}\] where \(\bar{p}_\parallel\) is the marginally stable parallel pressure obtained from Equation 8.54, Equation 8.55, and Equation 8.56. Equation 8.57 is called the pressure relaxation term due to the instabilities, with the relaxation time \(\tau\) related to their growth rates.

These equations provide a self-consistent description of phenomena that satisfy all the various assumptions we have made, namely:

1. The plasma is charge-neutral since characteristic lengths are much longer than a Debye length;
2. The characteristic velocity of the phenomenon under consideration is slow compared to the speed of light;
3. The pressure and density gradients are parallel, so there is no electrothermal EMF;
4. The time scale is long compared to both the electron and ion cyclotron periods.

Even though these assumptions are self-consistent, they may not accurately portray a real plasma and so MHD models, while intuitively appealing, must be used with caution.

Finally, it is worth mentioning that MHD plasmas can be categorized yet another way, namely according to the relative importance of the magnetic force \(\mathbf{J}\times\mathbf{B}\) compared to the hydrodynamic force \(\nabla\cdot\overleftrightarrow{P}\). If the magnetic force is negligible compared to the hydrodynamic force, then there is not much point in using MHD because in this case the system of equations is simply classical hydrodynamics. Thus, the only non-trivial MHD situations are where the magnetic and hydrodynamic forces are of comparable importance or where the magnetic force is much more important than the hydrodynamic force. Using Equation 8.48, the nominal ratio of the hydrodynamic force to the magnetic force is defined as \[ \beta = \frac{p}{B^2/2\mu_0} \sim \frac{\nabla\cdot\overleftarrow{P}}{\mathbf{J}\times\mathbf{B}} \sim \frac{p/L}{B^2/2\mu_0 L} \]

The characteristic gradient scale length \(L\) is assumed to be comparable for both types of forces and so cancels out in the comparison. Low-\(\beta\) plasmas are those where \(B^2/2\mu_0\) is much larger than \(p\) so the hydrodynamic force is negligible compared to the magnetic force, whereas \(\beta = \mathcal{O}(1)\) plasmas are those where the magnetic and hydrodynamic forces are comparable. The different \(\beta\) regime has important implications on the validity of many approximations seen in later chapters.

8.6.3 Hall MHD

As an extension to the ideal/resistive MHD model, Hall MHD decouples the electron and ion motions by retrieving the Hall term in the generalized Ohm’s law \[ \mathbf{E} = -\mathbf{u}\times\mathbf{B}+\eta\mathbf{J}+\frac{1}{en}\mathbf{J}\times\mathbf{B} - \frac{\nabla p_e}{ne} \tag{8.58}\] where the first term on the right-hand-side is the convective term, the second is the resistive term, the third is the Hall term, and the fourth is the electron pressure gradient term. The electron pressure scalar is simplified from the electron pressure tensor, which is obtained from an independent electron pressure equation.

The Hall MHD equations are \[ \begin{aligned} \frac{\partial\rho}{\partial t}&=-\nabla\cdot \left(\rho \mathbf{u} \right) \\ \frac{\partial (\rho\mathbf{u})}{\partial t}&=-\nabla\cdot \left(\rho\mathbf{u}\mathbf{u}+(p+p_e)\bar{\bar{I}}+\frac{B^2}{2\mu_0}\bar{\bar{I}}-\frac{\mathbf{B}\mathbf{B}}{\mu_0} \right) \\ \frac{\partial e}{\partial t}&=-\nabla\cdot\left[ \left(\epsilon+p\right)\mathbf{u} +(\epsilon_e + p_e)\mathbf{u}_e+\mathbf{u}_e\cdot\left(\frac{\mathbf{B}^2}{\mu_0}\bar{\bar{I}}-\frac{\mathbf{B}\mathbf{B}}{\mu_0}\right) - \mathbf{B}\times\eta\mathbf{J} \right] \\ \frac{\partial \mathbf{B}}{\partial t} &= -\nabla\times\left[ \mathbf{u}_e \times\mathbf{B} + \eta\mathbf{J} - \frac{\nabla p_e}{ne}\right] \\ \frac{\partial p_e}{\partial t} + \nabla\cdot \left(p_e\mathbf{u}_e \right)&=-(\gamma-1)p_e\nabla\cdot\mathbf{u}_e \end{aligned} \tag{8.59}\] where \(\bar{\bar{I}}\) is the identity matrix, \(\rho\) is the mass density, \(\mathbf{u}\) is the plasma bulk velocity, \(\mathbf{B}\) is the magnetic field, \(p_e\) is the electron pressure, \(p\) is the ion thermal pressure, \(\mathbf{J}=\nabla\times\mathbf{B}/\mu_0\) is the current density, \(\epsilon=\rho\mathbf{u}^2/2\) is the hydrodynamic energy density, and \(\epsilon_e=\rho\mathbf{u}_e^2/2\) is the electron hydrodynamic energy density. The Hall velocity is defined as \[ \mathbf{v}_H=-\frac{\mathbf{J}}{ne} \tag{8.60}\] and the electron bulk velocity is given by \[ \mathbf{u}_e = \mathbf{u} + \mathbf{v}_H \tag{8.61}\]

The total energy density is (I may have a coefficient error here for \(p_e\)!) \[ e= \frac{1}{2}\rho \mathbf{u}^2+\frac{1}{\gamma-1}(p+p_e)+\frac{\mathbf{B}^2}{2\mu_0} \] where \(\gamma\) is the adiabatic index. Note that in Hall MHD only \((\rho,\mathbf{u},\mathbf{B},p,p_e)\) are unknowns; all others are derived quantities.

One important thing to recognize is that while ideal MHD has no intrinsic scales, Hall MHD introduces a new characteristic length scale into the system – the ion skin depth \(d_i\). In ideal MHD, when the plasma beta is low, the speed of the shear Alfvén mode and the fast magnetosonic mode become approximately equal. This is known as degeneracy, which means it is difficult to distinguish between the two modes, as they essentially travel at the same speed. In Hall MHD, when the wavelength of the waves becomes comparable to \(d_i\), the distinction between the shear Alfvén and the fast magnetosonic modes becomes clear:

1. Fast magnetosonic waves (whistler waves) are affected by electron dynamics and propagate at higher frequencies within the Hall-MHD framework.
2. Shear Alfvén waves (kinetic Alfvén waves, or Alfvén-ion-cyclotron waves) are predominantly driven by ion dynamics and tend to propagate at lower frequencies.

8.7 Fluid Drifts Perpendicular to B

Since a fluid element is composed of many individual particles, one would expect the fluid to have drifts perpendicular to \(\mathbf{B}\) if the individual guiding centers have such drifts. However, since the \(\nabla p\) term appears only in the fluid equations, there is a drift associated with it which the fluid elements have but the particles do not have. For each species, we have an equation of motion \[ mn\Big[ \underbrace{\frac{\partial\mathbf{u}}{\partial t}}_{1}+\underbrace{(\mathbf{u}\cdot\nabla)\mathbf{u}}_{2}\Big] = q n (\mathbf{E}+\underbrace{\mathbf{u}\times\mathbf{B}}_{3}) - \nabla p \tag{8.62}\]

Consider the ratio of term 1 to term 3: \[ \frac{\text{term }1}{\text{term }3} \approx \bigg\lvert\frac{mni\omega v_\perp}{qnv_\perp B}\bigg\rvert \approx \frac{\omega}{\omega_c} \]

Here we have taken \(\partial/\partial t=i\omega\) and are concerned only with \(\mathbf{v}_\perp\). For drifts slow compared with the time scale associated with \(\omega_c\), we may neglect term 1. We shall also neglect the \((\mathbf{u}\cdot\nabla)\mathbf{u}\) term and show a posteriori that this is all right. Let \(\mathbf{E}\) and \(\mathbf{B}\) be uniform, but let \(n\) and \(p\) have a gradient. This is the usual situation in a magnetically confined plasma column (Fig. 3.4 ADD FIGURE!). Taking the cross product of Equation 8.62 with B, we have (neglecting the left-hand side)

\[ \begin{aligned} 0 &= qn[\mathbf{E}\times\mathbf{B}+(\mathbf{u}_\perp\times\mathbf{B})\times\mathbf{B}] - \nabla p\times\mathbf{B} \\ &= qn[\mathbf{E}\times\mathbf{B}+\mathbf{B}(\cancel{\mathbf{u}_\perp\cdot\mathbf{B}})-\mathbf{u}_\perp\mathbf{B}^2] - \nabla p\times\mathbf{B} \\ &= qn[\mathbf{E}\times\mathbf{B}-\mathbf{u}_\perp\mathbf{B}^2] - \nabla p\times\mathbf{B} \end{aligned} \]

Therefore, \[ \mathbf{u}_\perp = \frac{\mathbf{E}\times\mathbf{B}}{B^2} -\frac{\nabla p \times\mathbf{B}}{qnB^2} \equiv \mathbf{u}_E + \mathbf{u}_D \] where \[ \mathbf{u}_E\equiv\frac{\mathbf{E}\times\mathbf{B}}{B^2}\quad\mathbf{E}\times\mathbf{B}\,\text{drift} \]

\[ \mathbf{u}_D\equiv-\frac{\nabla p\times\mathbf{B}}{qnB^2}\quad\text{Diamagnetic drift} \tag{8.63}\]

The drift \(\mathbf{u}_E\) is the same as for guiding centers, but there is now a new drift \(\mathbf{u}_D\), called the diamagnetic drift. Since \(\mathbf{u}_D\) is perpendicular to the direction of the gradient, our neglect of \((\mathbf{u}\cdot\nabla)\mathbf{u}\) is justified if \(\mathbf{E}=0\). If \(\mathbf{E}=-\nabla\phi\neq 0\) (i.e. a potential field), \((\mathbf{u}\cdot\nabla)\mathbf{u}\) is still zero if \(\nabla\phi\) and \(\nabla p\) are in the same direction; otherwise, there could be a more complicated solution involving \((\mathbf{u}\cdot\nabla)\mathbf{u}\).

With the help of Equation 8.12, we can write the diamagnetic drift as \[ \mathbf{u}_D = \pm \frac{\gamma k_B T}{eB}\frac{\hat{z}\times\nabla n}{n} \tag{8.64}\]

In particular, for an isothermal plasma in the geometry of Fig.3.4 (ADD IT!), in which \(\nabla n = \partial n\partial r\hat{r} = n^\prime\hat{r}\) (\(n^\prime < 0\)), we have the following formulas familiar to experimentalists who have worked with Q-machines \[ \begin{aligned} \mathbf{u}_{Di} &= \frac{k_B T_i}{eB}\frac{n^\prime}{n}\hat{\theta} \\ \mathbf{u}_{De} &= -\frac{k_B T_e}{eB}\frac{n^\prime}{n}\hat{\theta} \end{aligned} \]

The physical reason for this drift can be seen from Figure 8.2. Here we have drawn the orbits of ions gyrating in a magnetic field. There is a density gradient toward the left, as indicated by the density of orbits. Through any fixed volume element there are more ions moving downward than upward, since the downward-moving ions come from a region of higher density. There is, therefore, a fluid drift perpendicular to \(\nabla n\) and \(\mathbf{B}\), even though the guiding centers are stationary. The diamagnetic drift reverses sign with \(q\) because the direction of gyration reverses. The magnitude of \(\mathbf{u}_D\) does not depend on mass because the \(m^{-1/2}\) dependence of the velocity is cancelled by the \(m^{1/2}\) dependence of the Larmor radius—less of the density gradient is sampled during a gyration if the mass is small.

Figure 8.2: Origin of the diamagnetic drift.

Since ions and electrons drift in opposite directions, there is a diamagnetic current. For \(\gamma=Z=1\), this is given by \[ \mathbf{J}_D = ne(\mathbf{u}_{Di} + \mathbf{u}_{De}) = \frac{\mathbf{B}}{B^2}\times\nabla\left( p_i + p_e \right) = (k_B T_i + k_B T_e)\frac{\mathbf{B}\times\nabla n}{B^2} \tag{8.65}\]

In the particle picture, one would not expect to measure a current if the guiding centers do not drift. In the fluid picture, the current \(\mathbf{J}_D\) flows wherever there is a pressure gradient. These two viewpoints can be reconciled if one considers that all experiments must be carried out in a finite-sized plasma. Suppose the plasma were in a rigid box (Fig. 3.6 ADD IT!). If one were to calculate the current from the single-particle picture, one would have to take into account the particles at the edges which have cycloidal paths. Since there are more particles on the left than on the right, there is a net current downward, in agreement with the fluid picture.

The reader may not be satisfied with this explanation because it was necessary to specify reflecting walls. If the walls were absorbing or if they were removed, one would find that electric fields would develop because more of one species—the one with larger Larmor radius—would collected than the other. Then the guiding centers would drift, and the simplicity of the model would be lost. Alternatively, one could imagine trying to measure the diamagnetic current with a current probe (Fig. 3.7 ADD IT!). This is just a transformer with a core of magnetic material. The primary winding is the plasma current threading the core, and the secondary is a multiturn winding all around the core. Let the whole thing be infinitesimally thin, so it does not intercept any particles. It is clear from Fig. 3.7 that a net upward current would be measured, there being higher density on the left than on the right, so that the diamagnetic current is a real current. From this example, one can see that it can be quite tricky to work with the single-particle picture. The fluid theory usually gives the right results when applied straightforwardly, even though it contains “fictitious” drifts like the diamagnetic drift.

What about the grad-B and curvature drifts which appeared in the single-particle picture? The curvature drift also exists in the fluid picture, since the centrifugal force is felt by all the particles in a fluid element as they move around a bend in the magnetic field. A term \(\bar{F}_{cf} = \overline{nmv_\parallel^2}/R_c = nk_B T_\parallel / R_c\) has to be added to the right-hand side of the fluid equation of motion. This is equivalent to a gravitational force \(mng\), with \(g = k_BT_\parallel/m R_c\), and leads to a drift \(\mathbf{u}_g=(m/q)(\mathbf{g}\times\mathbf{B})/B^2\), as in the particle picture (Equation 7.5).

Figure 8.3: In a nonuniform B field the guiding centers drift but the fluid elements do not.

The grad-B drift, however, does not exist for fluids. It can be shown on thermodynamic grounds that a magnetic field does not affect a Maxwellian distribution. This is because the Lorentz force is perpendicular to \(\mathbf{v}\) and cannot change the energy of any particle. The most probable distribution \(f(\mathbf{v})\) in the absence of \(\mathbf{B}\) is also the most probable distribution in the presence of \(\mathbf{B}\). If \(f(\mathbf{v})\) remains Maxwellian in a nonuniform \(\mathbf{B}\) field, and there is no density gradient, then the net momentum carried into any fixed fluid element is zero. There is no fluid drift even though the individual guiding centers have drifts; the particle drift in any fixed fluid element cancel out. To see this pictorially, consider the orbits of two particles moving through a fluid element in a nonuniform \(\mathbf{B}\) field (Figure 8.3). Since there is no \(\mathbf{E}\) field, the Larmor radius changes only because of the gradient in B; there is no acceleration, and the particle energy remains constant during the motion. If the two particles have the same energy, they will have the same velocity and Larmor radius while inside the fluid element. There is thus a perfect cancellation between particle pairs when their velocities are added to give the fluid velocity.

When there is a nonuniform \(\mathbf{E}\) field, it is not easy to reconcile the fluid and particle pictures. Then the finite-Larmor-radius effect of Section 7.3 causes both a guiding center drift and a fluid drift, but these are not the same; in fact, they have opposite signs! The particle drift was calculated in Chapter 7, and the fluid drift can be calculated from the off-diagonal elements of \(\mathbf{P}\). It is extremely difficult to explain how the finite-Larmor-radius effects differ. A simple picture like Fig. 3.6 will not work because one has to take into account subtle points like the following: in the absence of a density gradient, the density of guiding centers is not the same as the density of particles! (???)

(I need to think carefully about these pictures.)

8.8 Fluid Drifts Parallel to B

The \(z\) component of the fluid equation of motions is \[ mn\left[ \frac{\partial u_z}{\partial t} + (\mathbf{u}\cdot\nabla)u_z\right] = qn E_z - \frac{\partial\phi}{\partial z} \tag{8.66}\]

The convective term can often be neglected because it is much smaller than the \(\partial u_z/\partial t\) term. We shall avoid complicated arguments here and simply consider cases in which \(u_z\) is spatially uniform. Using Equation 8.66, we have \[ \frac{\partial u_z}{\partial t} = \frac{q}{m}E_z - \frac{\gamma k_B T}{mn}\frac{\partial n}{\partial z} \tag{8.67}\]

This shows that the fluid is accelerated along \(\mathbf{B}\) under the combined electrostatic and pressure gradient forces. A particularly important result is obtained by applying Equation 8.67 to massless electrons. Taking the limit \(m\rightarrow 0\) and specifying \(q=-e\) and \(\mathbf{E}=-\nabla\phi\), we have \[ qE_z = e\frac{\partial\phi}{\partial z} = \frac{\gamma k_B T_e}{n}\frac{\partial n}{\partial z} \]

Electrons are so mobile that their heat conductivity is almost infinite. We may then assume isothermal electrons and take \(\gamma = 1\). Integrating, we have \[ e\phi = k_B T_e \ln n + C \] or \[ n = n_0 \exp(e\phi/k_B T_e) \tag{8.68}\]

This is just the Boltzmann relation for electrons.

What this means physically is that electrons, being light, are very mobile and would be accelerated to high energies very quickly if there were a net force on them. Since electrons cannot leave a region en masse without leaving behind a large ion charge, the electrostatic and pressure gradient forces on the electrons must be closely in balance. This condition leads to the Boltzmann relation. Note that Equation 8.68 applies to each magnetic field lines separately. Different field lines may be charged to different potentials arbitrarily unless a mechanism is provided for the electrons to move across \(\mathbf{B}\). The conductors on which field lines terminate can provide such a mechanism, and the experimentalist has to take these end effects into account carefully.

Figure 8.4: Physical reason for the Boltzmann relation between density and potential.

Figure 8.4 shows graphically what occurs when there is a local density clump in the plasma. Let the density gradient be toward the center of the diagram, and suppose \(k_B T\) is constant. There is then a pressure gradient toward the center. Since the plasma is quasineutral, the gradient exists for both the electron and ion fluids. Consider the pressure gradient force \(\mathbf{F}_p\) on the electron fluid. It drives the mobile electrons away from the center, leaving the ions behind. The resulting positive charge generates a field \(\mathbf{E}\) whose force \(\mathbf{F}_E\) on the electrons opposes \(\mathbf{F}_p\). Only when \(\mathbf{F}_E = -\mathbf{F}_p\) is a steady state achieved. If \(\mathbf{B}\) is a constant, \(\mathbf{E}\) is an electrostatic field \(\mathbf{E}=-\nabla\phi\), and \(\phi\) must be large at the center, where \(n\) is large. This is just what Equation 8.68 tells us. The deviation from stric neutrality adjusts itself so that there is just enough charge to set up the \(\mathbf{E}\) field required to balance the forces on the electrons.

8.9 The Plasma Approximation

The previous example reveals an important characteristic of plasmas that has wide application. We are used to solving for \(\mathbf{E}\) from Poisson’s equation when we are given the charge density \(\sigma\). In a plasma, the opposite procedure is generally used. \(\mathbf{E}\) is found from the equations of motion, and Poisson’s equation is used only to find \(\sigma\). The reason is that a plasma has an overriding tendency to remain neutral. If the ions move, the electrons will follow. \(\mathbf{E}\) must adjust itself so that the orbits of the electrons and ions preserve neutrality. The charge density is of secondary importance; it will adjust itself so that Poisson’s equation is satisfied. This is true, of course, only for low-frequency motions in which the electron inertia is not a factor.

In a plasma, it is usually possible to assume \(n_i=n_e\) and \(\nabla\cdot\mathbf{E}\neq 0\) at the same time. We shall call this the plasma approximation. It is a fundamental trait of plasmas, one which is difficult for the novice to understand. Do not use Poisson’s equation to obtain \(\mathbf{E}\) unless it is unavoidable! In the set of fluid equations Equation 8.13, we may now eliminate Poisson’s equation and also eliminate one of the unknowns by setting \(n_i = n_e = n\).

The plasma approximation is almost the same as the condition of quasineutrality discussed earlier but has a more exact meaning. Whereas quasineutrality refers to a general tendency for a plasma to be neutral in its state of rest, the plasma approximation is a mathematical shortcut that one can use even for wave motions. As long as these motions are slow enough that both ions and electrons have time to move, it is a good approximation to replace Poisson’s equation by the equation \(n_i=n_e\). Of course, if only one species can move and the other cannot follow, such as in high-frequency electron waves, then the plasma approximation is not valid, and \(\mathbf{E}\) must be found from Maxwell’s Equation 8.1 rather than from the ion and electron equations of motion. Ampère’s law with the displacement current retained gives the time evolution of \(\mathbf{E}\):

\[ \frac{\mathrm{d}\mathbf{E}}{\mathrm{d}t} = \frac{1}{\epsilon_0}\left( \frac{\nabla\times\mathbf{B}}{\mu_0} - \mathbf{J} \right) \]

We shall return to the question of the validity of the plasma approximation when we come to the theory of ion waves. At that time, it will become clear why we had to use Poisson’s equation in the derivation of Debye shielding. For electron waves, we can also show that with a simple “modified” electron test particle model that updates the electric field while keeping the magnetic field constant, we can retain the plasma oscillation and electron Berstein wave.

Bellan, Paul M. 2008. Fundamentals of Plasma Physics. Cambridge University Press.

Chen, Francis F. 2016. Introduction to Plasma Physics and Controlled Fusion. Vol. 1. Springer.

Chew, GF, ML Goldberger, and FE Low. 1956. “The Boltzmann Equation and the One-Fluid Hydromagnetic Equations in the Absence of Particle Collisions.” Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 236 (1204): 112–18. https://doi.org/10.1098/rspa.1956.0116.

Wang, Liang, Ammar H Hakim, Jonathan Ng, Chuanfei Dong, and Kai Germaschewski. 2020. “Exact and Locally Implicit Source Term Solvers for Multifluid-Maxwell Systems.” Journal of Computational Physics 415: 109510. https://doi.org/10.1016/j.jcp.2020.109510.

---

1. We have put a subscript m on the permeability to distinguish it from the adiabatic invariant \(\mu\).

2. the Darwin approximation ignores light waves by neglecting \(\partial\mathbf{E}\partial t\) in Amperè’s law.

3. In the case of discontinuities, adiabatic laws fail to describe the jump, and the energy equation needs to be used due to its conservation property.

4. This is slightly different in some terms compared to Equation 8.14. Meng Xing mentioned in the BATSRUS anistropic MHD paper that the polytropic index is taken to be 5/3. Also note that here the unknown is the total pressure \(p\) instead of perpendicular pressure \(p_\perp\).