Categories: Physics, Plasma physics.

Two-fluid equations

The two-fluid model describes a plasma as two separate but overlapping fluids, one for ions and one for electrons. Instead of tracking individual particles, it gives the dynamics of fluid elements dV\dd{V} (i.e. small “blobs”). These blobs are assumed to be much larger than the Debye length, such that electromagnetic interactions between nearby blobs can be ignored.

From Newton’s second law, we know that the velocity v\vb{v} of a particle with mass mm and charge qq is as follows, when subjected only to the Lorentz force:

mdvdt=q(E+v×B)\begin{aligned} m \dv{\vb{v}}{t} = q (\vb{E} + \vb{v} \cross \vb{B}) \end{aligned}

From here, the derivation is similar to that of the Navier-Stokes equations. We replace d/dt\idv{}{t} with a material derivative D/Dt\mathrm{D}/\mathrm{D}t, and define u\vb{u} as the blob’s center-of-mass velocity:

mnDuDt=qn(E+u×B)\begin{aligned} m n \frac{\mathrm{D} \vb{u}}{\mathrm{D} t} = q n (\vb{E} + \vb{u} \cross \vb{B}) \end{aligned}

Where we have multiplied by the number density nn of the particles. Due to particle collisions in the fluid, stresses become important. Therefore, we include the Cauchy stress tensor P^\hat{P}, leading to the following two equations:

miniDuiDt=qini(E+ui×B)+P^imeneDueDt=qene(E+ue×B)+P^e\begin{aligned} m_i n_i \frac{\mathrm{D} \vb{u}_i}{\mathrm{D} t} &= q_i n_i (\vb{E} + \vb{u}_i \cross \vb{B}) + \nabla \cdot \hat{P}_i{}^\top \\ m_e n_e \frac{\mathrm{D} \vb{u}_e}{\mathrm{D} t} &= q_e n_e (\vb{E} + \vb{u}_e \cross \vb{B}) + \nabla \cdot \hat{P}_e{}^\top \end{aligned}

Where the subscripts ii and ee refer to ions and electrons, respectively. Finally, we also account for momentum transfer between ions and electrons due to Rutherford scattering, leading to these two-fluid momentum equations:

miniDuiDt=qini(E+ui×B)+P^ifiemini(uiue)meneDueDt=qene(E+ue×B)+P^efeimene(ueui)\begin{aligned} \boxed{ \begin{aligned} m_i n_i \frac{\mathrm{D} \vb{u}_i}{\mathrm{D} t} &= q_i n_i (\vb{E} + \vb{u}_i \cross \vb{B}) + \nabla \cdot \hat{P}_i{}^\top - f_{ie} m_i n_i (\vb{u}_i - \vb{u}_e) \\ m_e n_e \frac{\mathrm{D} \vb{u}_e}{\mathrm{D} t} &= q_e n_e (\vb{E} + \vb{u}_e \cross \vb{B}) + \nabla \cdot \hat{P}_e{}^\top - f_{ei} m_e n_e (\vb{u}_e - \vb{u}_i) \end{aligned} } \end{aligned}

Where fief_{ie} is the mean frequency at which an ion collides with electrons, and vice versa for feif_{ei}. For simplicity, we assume that the plasma is isotropic and that shear stresses are negligible, in which case the stress term can be replaced by the gradient p- \nabla p of a scalar pressure pp:

miniDuiDt=qini(E+ui×B)pifiemini(uiue)meneDueDt=qene(E+ue×B)pefeimene(ueui)\begin{aligned} m_i n_i \frac{\mathrm{D} \vb{u}_i}{\mathrm{D} t} &= q_i n_i (\vb{E} + \vb{u}_i \cross \vb{B}) - \nabla p_i - f_{ie} m_i n_i (\vb{u}_i - \vb{u}_e) \\ m_e n_e \frac{\mathrm{D} \vb{u}_e}{\mathrm{D} t} &= q_e n_e (\vb{E} + \vb{u}_e \cross \vb{B}) - \nabla p_e - f_{ei} m_e n_e (\vb{u}_e - \vb{u}_i) \end{aligned}

Next, we demand that matter is conserved. In other words, the rate at which particles enter/leave a volume VV must be equal to the flux through the enclosing surface SS:

0=tVndV+SnudS=V(nt+(nu))dV\begin{aligned} 0 &= \pdv{}{t}\int_V n \dd{V} + \oint_S n \vb{u} \cdot \dd{\vb{S}} = \int_V \Big( \pdv{n}{t} + \nabla \cdot (n \vb{u}) \Big) \dd{V} \end{aligned}

Where we have used the divergence theorem. Since VV is arbitrary, we can remove the integrals, leading to the following continuity equations:

nit+(niui)=0net+(neue)=0\begin{aligned} \boxed{ \pdv{n_i}{t} + \nabla \cdot (n_i \vb{u}_i) = 0 \qquad \quad \pdv{n_e}{t} + \nabla \cdot (n_e \vb{u}_e) = 0 } \end{aligned}

These are 8 equations (2 scalar continuity, 2 vector momentum), but 16 unknowns ui\vb{u}_i, ue\vb{u}_e, E\vb{E}, B\vb{B}, nin_i, nen_e, pip_i and pep_e. We would like to close this system, so we need 8 more. An obvious choice is Maxwell’s equations, in particular Faraday’s and Ampère’s law (since Gauss’ laws are redundant; see the article on Maxwell’s equations):

×E=Bt×B=μ0(niqiui+neqeue+ε0Et)\begin{aligned} \boxed{ \nabla \cross \vb{E} = - \pdv{\vb{B}}{t} \qquad \quad \nabla \cross \vb{B} = \mu_0 \Big( n_i q_i \vb{u}_i + n_e q_e \vb{u}_e + \varepsilon_0 \pdv{\vb{E}}{t} \Big) } \end{aligned}

Now we have 14 equations, so we need 2 more, for the pressures pip_i and pep_e. This turns out to be the thermodynamic equation of state: for quasistatic, reversible, adiabatic compression of a gas with constant heat capacity (i.e. a calorically perfect gas), it turns out that:

DDt(pVγ)=0γCPCV=N+2N\begin{aligned} \frac{\mathrm{D}}{\mathrm{D} t} \big( p V^\gamma \big) = 0 \qquad \quad \gamma \equiv \frac{C_P}{C_V} = \frac{N + 2}{N} \end{aligned}

Where γ\gamma is the heat capacity ratio, and can be calculated from the number of degrees of freedom NN of each particle in the gas. In a fully ionized plasma, N=3N = 3.

The density n1/Vn \propto 1/V, so since pVγp V^\gamma is constant in time, for some constant CC:

DDt(pnγ)=0    p=Cnγ\begin{aligned} \frac{\mathrm{D}}{\mathrm{D} t} \Big( \frac{p}{n^\gamma} \Big) = 0 \quad \implies \quad p = C n^\gamma \end{aligned}

In the two-fluid model, we thus have the following two equations of state, giving us a set of 16 equations for 16 unknowns:

DDt(piniγ)=0DDt(peneγ)=0\begin{aligned} \boxed{ \frac{\mathrm{D}}{\mathrm{D} t} \Big( \frac{p_i}{n_i^\gamma} \Big) = 0 \qquad \quad \frac{\mathrm{D}}{\mathrm{D} t} \Big( \frac{p_e}{n_e^\gamma} \Big) = 0 } \end{aligned}

Note that from the relation p=Cnγp = C n^\gamma, we can calculate the p\nabla p term in the momentum equation, using simple differentiation and the ideal gas law:

p=Cnγ    p=γCnγnn=γpnn=γkBTn\begin{aligned} p = C n^\gamma \quad \implies \quad \nabla p = \gamma \frac{C n^{\gamma}}{n} \nabla n = \gamma p \frac{\nabla n}{n} = \gamma k_B T \nabla n \end{aligned}

Note that the ideal gas law was not used immediately, to allow for γ1\gamma \neq 1.

Fluid drifts

The momentum equations reduce to the following if we assume the flow is steady u/t=0\ipdv{\vb{u}}{t} = 0, and neglect electron-ion momentum transfer on the right:

mini(ui)uiqini(E+ui×B)pimene(ue)ueqene(E+ue×B)pe\begin{aligned} m_i n_i (\vb{u}_i \cdot \nabla) \vb{u}_i &\approx q_i n_i (\vb{E} + \vb{u}_i \cross \vb{B}) - \nabla p_i \\ m_e n_e (\vb{u}_e \cdot \nabla) \vb{u}_e &\approx q_e n_e (\vb{E} + \vb{u}_e \cross \vb{B}) - \nabla p_e \end{aligned}

We take the cross product with B\vb{B}, which leaves only the component u\vb{u}_\perp of u\vb{u} perpendicular to B\vb{B} in the Lorentz term:

0=qn(E+u×B)×Bp×Bmn((u)u)×B=qn(E×BuB2)p×Bmn((u)u)×B\begin{aligned} 0 &= q n (\vb{E} + \vb{u}_\perp \cross \vb{B}) \cross \vb{B} - \nabla p \cross \vb{B} - m n \big( (\vb{u} \cdot \nabla) \vb{u} \big) \cross \vb{B} \\ &= q n (\vb{E} \cross \vb{B} - \vb{u}_\perp B^2) - \nabla p \cross \vb{B} - m n \big( (\vb{u} \cdot \nabla) \vb{u} \big) \cross \vb{B} \end{aligned}

Isolating for u\vb{u}_\perp tells us that the fluids drifts perpendicularly to B\vb{B}, with velocity u\vb{u}_\perp:

u=E×BB2p×BqnB2m((u)u)×BqB2\begin{aligned} \vb{u}_\perp = \frac{\vb{E} \cross \vb{B}}{B^2} - \frac{\nabla p \cross \vb{B}}{q n B^2} - \frac{m \big( (\vb{u} \cdot \nabla) \vb{u} \big) \cross \vb{B}}{q B^2} \end{aligned}

The last term is often neglected, which turns out to be a valid approximation if E=0\vb{E} = 0, or if E\vb{E} is parallel to p\nabla p. The first term is the familiar E×B\vb{E} \cross \vb{B} drift vE\vb{v}_E from guiding center theory, and the second term is called the diamagnetic drift vD\vb{v}_D:

vE=E×BB2vD=p×BqnB2\begin{aligned} \boxed{ \vb{v}_E = \frac{\vb{E} \cross \vb{B}}{B^2} } \qquad \quad \boxed{ \vb{v}_D = - \frac{\nabla p \cross \vb{B}}{q n B^2} } \end{aligned}

It is called diamagnetic because it creates a current that induces a magnetic field opposite to the original B\vb{B}. In a quasi-neutral plasma qene=qiniq_e n_e = - q_i n_i, the current density J\vb{J} is given by:

J=qene(vDevDi)=qene(pi×BqiniB2pe×BqeneB2)=B×(pi+pe)B2\begin{aligned} \vb{J} = q_e n_e (\vb{v}_{De} - \vb{v}_{Di}) = q_e n_e \Big( \frac{\nabla p_i \cross \vb{B}}{q_i n_i B^2} - \frac{\nabla p_e \cross \vb{B}}{q_e n_e B^2} \Big) = \frac{\vb{B} \cross \nabla (p_i + p_e)}{B^2} \end{aligned}

Using the ideal gas law p=kBTnp = k_B T n, this can be rewritten as follows:

J=kBB×(Tini+Tene)B2\begin{aligned} \vb{J} = k_B \frac{\vb{B} \cross \nabla (T_i n_i + T_e n_e)}{B^2} \end{aligned}

Curiously, vD\vb{v}_D does not involve any net movement of particles, because a pressure gradient does not necessarily cause particles to move. Instead, there is a higher density of gyration paths in the high-pressure region, so that the particle flux through a reference plane is higher. This causes the fluid elements to drift, but not the guiding centers.


  1. F.F. Chen, Introduction to plasma physics and controlled fusion, 3rd edition, Springer.
  2. M. Salewski, A.H. Nielsen, Plasma physics: lecture notes, 2021, unpublished.