Resistive MHD: deriving the equations¶
Advanced level
This page builds on MHD and frozen-in flux, Vector and scalar potentials and Reconnection rate and Ohm's law. The equations derived here are the ones the Simulation Lab code solves.
Resistive magnetohydrodynamics (MHD) describes a plasma as a single electrically conducting fluid with a small but finite resistivity. It is the simplest model in which reconnection can happen at all: ideal MHD freezes field lines into the plasma (frozen-in flux), and resistivity is the one ingredient that lets them slip. The Sweet–Parker model, the tearing instability and the plasmoid instability are all resistive-MHD results.
This page derives the equations step by step, says what is thrown away at each step and why, and ends with the dimensionless form that a computer solves.
The plan¶
| Step | What we do | What we assume |
|---|---|---|
| 1 | Describe each species (ions, electrons) by fluid moments | Many particles per Debye sphere |
| 2 | Add the species into one fluid | Quasi-neutrality \(n_i \approx n_e\) |
| 3 | Get Ohm's law from the electron equation | Scales much larger than the ion inertial length \(d_i\) |
| 4 | Simplify Maxwell's equations | Speeds much less than \(c\) |
| 5 | Combine into the induction equation | — |
| 6 | Close the system with an energy equation | Isotropic pressure, ideal gas |
| 7 | Check energy conservation | — |
| 8 | Rewrite with the vector potential and in dimensionless form | — |
Step 1: from particles to fluids¶
Each species \(s\) (ions \(i\), electrons \(e\)) has a distribution function \(f_s(\mathbf{x},\mathbf{v},t)\): the number of particles per unit volume of position and velocity space. Its velocity moments are the fluid quantities:
These are the number density, the mean (bulk) velocity and the pressure tensor (the momentum flux carried by random thermal motion). Multiplying the kinetic equation for \(f_s\) by \(1\) and by \(m_s\mathbf{v}\) and integrating over velocity gives the two fluid equations of species \(s\):
The first says particles are neither created nor destroyed. The second is Newton's law for a fluid element: Lorentz force, minus the pressure force, plus \(\mathbf{R}_s\), the momentum species \(s\) gains from collisions with the other species. Collisions conserve total momentum, so \(\mathbf{R}_i = -\mathbf{R}_e\).
Why the pressure tensor becomes a scalar
When collisions (or wave scattering) keep the velocity distribution close to isotropic, \(\mathsf{P}_s \approx p_s\mathsf{I}\) and \(\nabla\cdot\mathsf{P}_s = \nabla p_s\). The small anisotropic part that remains is viscosity, which we keep as a term with a kinematic viscosity \(\nu\).
Step 2: one fluid¶
Define the single-fluid variables
Quasi-neutrality. On scales much larger than the Debye length the plasma shields any charge imbalance (see What is a plasma?), so \(n_i = n_e = n\) to excellent accuracy. The net charge density \(e(n_i-n_e)\) is tiny, although the current is not.
Mass. Multiply each continuity equation by \(m_s\) and add:
Momentum. Add the two momentum equations. Three things happen:
- The electric forces cancel: \(e n_i\mathbf{E} - e n_e\mathbf{E} = 0\) by quasi-neutrality.
- The magnetic forces combine into \(e(n_i\mathbf{u}_i - n_e\mathbf{u}_e)\times\mathbf{B} = \mathbf{J}\times\mathbf{B}\).
- The collision terms cancel, \(\mathbf{R}_i + \mathbf{R}_e = 0\).
With the electron inertia neglected (\(m_e \ll m_i\)) and gravity \(\mathbf{g}\) added,
where \(\boldsymbol{\tau}\) is the viscous stress. For a constant kinematic viscosity \(\nu\) its force per unit mass is approximately \(\nu\left[\nabla^2\mathbf{v} + \tfrac13\nabla(\nabla\cdot\mathbf{v})\right]\).
Answer: where did the electric field go?
The electric force on a fluid element is the charge density times \(\mathbf{E}\). Quasi-neutrality makes the charge density essentially zero, so the bulk fluid feels no net electric force; it is pushed only through \(\mathbf{J}\times\mathbf{B}\). The electric field has not disappeared. It is determined by Ohm's law (next step), not by Gauss's law.
The \(\mathbf{J}\times\mathbf{B}\) force splits into magnetic pressure and tension (Bending flux tubes):
Step 3: Ohm's law¶
The one-fluid momentum equation does not contain \(\mathbf{E}\), so a separate equation must determine it. It comes from the electron momentum equation. Write the electron velocity in one-fluid variables, \(\mathbf{u}_e = \mathbf{v} - \mathbf{J}/(ne)\), and model the collisional drag as friction proportional to the relative velocity of electrons and ions:
where \(\nu_{ei}\) is the electron–ion collision frequency. Solving the electron equation for \(\mathbf{E}\) gives the generalized Ohm's law:
\(\eta_{\text{res}}\) is the electrical resistivity (units Ω·m). For a fully ionized plasma it is the Spitzer value, \(\eta_{\text{res}} \approx 5.2\times10^{-5}\,\ln\Lambda\,T_{\text{eV}}^{-3/2}\) Ω·m, which is independent of density and falls steeply with temperature.
Which terms survive? Compare each term with \(\mathbf{v}\times\mathbf{B}\) for a structure of size \(L\) and flows near the Alfvén speed:
| Term | Size relative to \(\mathbf{v}\times\mathbf{B}\) | Small when |
|---|---|---|
| Hall | \(d_i/L\) | \(L \gg d_i = c/\omega_{pi}\) (ion inertial length) |
| Electron pressure | \(\sim d_i/L\) (for \(\beta\sim1\)) | \(L \gg d_i\) |
| Electron inertia | \((d_e/L)^2\) | \(L \gg d_e = c/\omega_{pe}\) |
| Resistive | \(1/R_m\) | \(R_m \gg 1\) |
Resistive MHD keeps only the first term on the right:
From here on \(\eta\) is the magnetic diffusivity, the same \(\eta\) as on the frozen-in flux page.
The price of this step
Dropping the Hall and electron-pressure terms is valid only when every scale, including the diffusion region, is much larger than \(d_i\). In the magnetosphere the diffusion region is about \(d_i\) thick, so collisionless reconnection needs the extra terms (Reconnection rate and Ohm's law). Resistive MHD is still the right starting point: it is the large-scale limit, and it gives the Sweet–Parker and tearing results.
Step 4: low-frequency Maxwell equations¶
Why drop the displacement current? With \(E\sim vB\) (from Ohm's law) and a time scale \(L/v\),
The result is Ampère's law \(\mu_0\mathbf{J} = \nabla\times\mathbf{B}\). It is now an equation for \(\mathbf{J}\), not for \(\mathbf{E}\), and it implies \(\nabla\cdot\mathbf{J} = 0\): currents flow in closed loops. Gauss's law \(\nabla\cdot\mathbf{E} = \rho_q/\epsilon_0\) is not needed; it only fixes the tiny charge density that quasi-neutrality ignores.
Step 5: the induction equation¶
Solve Ohm's law for \(\mathbf{E}\), use Ampère's law, and substitute into Faraday's law:
For uniform \(\eta\), the identity \(\nabla\times\nabla\times\mathbf{B} = \nabla(\nabla\cdot\mathbf{B}) - \nabla^2\mathbf{B}\) with \(\nabla\cdot\mathbf{B}=0\) gives
The first term freezes the field into the plasma; the second lets it diffuse through the plasma at a rate set by \(\eta\). Taking the divergence shows \(\partial_t(\nabla\cdot\mathbf{B}) = 0\): if \(\nabla\cdot\mathbf{B} = 0\) initially, it stays zero.
Step 6: closing the system with an energy equation¶
The momentum equation contains \(p\), so we need an equation for it. For an ideal gas with ratio of specific heats \(\gamma\) (5/3 for a monatomic plasma), the internal energy per unit volume is \(p/(\gamma-1)\). The first law of thermodynamics for a moving fluid element is
where \(Q\) is the heating rate per unit volume. With no heating this is \(d(p\rho^{-\gamma})/dt = 0\): adiabatic compression. Using the continuity equation, \(d\rho/dt = -\rho\nabla\cdot\mathbf{v}\), so
- Ohmic heating \(\mu_0\eta J^2 = \mathbf{J}\cdot(\mathbf{E}+\mathbf{v}\times\mathbf{B})\) is the work done by the electric field in the frame moving with the plasma. It is where reconnection turns magnetic energy into heat.
- Viscous heating turns the kinetic energy of sheared and compressed flows into heat. In full Cartesian tensor notation, \(\Phi_\nu = 2\rho\nu\left[S_{ij}S_{ij} - \tfrac13(\nabla\cdot\mathbf{v})^2\right]\), where \(S_{ij} = \tfrac12(\partial_i v_j + \partial_j v_i)\) is the strain rate tensor. In general curvilinear coordinates without full Christoffel symbols, this is often represented in positive-definite scalar form as \(\Phi_\nu \approx \rho\nu\left[|\nabla\times\mathbf{v}|^2 + \tfrac43(\nabla\cdot\mathbf{v})^2\right]\).
- Conduction \(\mathbf{q} = -\kappa\nabla T\) carries heat from hot to cold.
The complete set¶
These are 8 equations (1 + 3 + 1 + 3) for 8 unknowns: \(\rho\), the three components of \(\mathbf{v}\), \(p\) and the three components of \(\mathbf{B}\). \(\mathbf{J}\) and \(\mathbf{E}\) are not separate unknowns; they follow from \(\mathbf{B}\) and \(\mathbf{v}\).
Step 7: where the energy goes¶
Dot the momentum equation with \(\mathbf{v}\), and dot Faraday's law with \(\mathbf{B}/\mu_0\) (using Poynting's theorem). With \(\mathbf{J}\cdot\mathbf{E} = \mathbf{v}\cdot(\mathbf{J}\times\mathbf{B}) + \mu_0\eta J^2\) (Problem 5), the three energy equations are
Every source on the right appears twice with opposite signs, so the total energy \(W = \rho v^2/2 + p/(\gamma-1) + B^2/2\mu_0\) (plus gravitational energy) is conserved; it only moves around. The exchange terms are the whole story of reconnection energetics:
| Exchange term | From | To | In reconnection |
|---|---|---|---|
| \(\mathbf{v}\cdot(\mathbf{J}\times\mathbf{B})\) | magnetic | kinetic | Tension of reconnected field lines drives the outflow jets |
| \(\mu_0\eta J^2\) | magnetic | thermal | Ohmic heating in the current sheet |
| \(p\,\nabla\cdot\mathbf{v}\) | thermal | kinetic | Expanding plasma does work; compression (\(\nabla\cdot\mathbf{v}<0\)) heats it instead |
| \(\Phi_\nu\) | kinetic | thermal | Viscous damping of the jets |
Answer: why does reconnection need resistivity to release energy?
The only direct sink of magnetic energy that is not reversible is \(\mu_0\eta J^2\). The \(\mathbf{v}\cdot(\mathbf{J}\times\mathbf{B})\) exchange works in both directions, as in an Alfvén wave. In ideal MHD (\(\eta = 0\)) field lines cannot change connectivity, so the stressed field cannot relax to a lower-energy topology. A tiny \(\eta\) acting in a thin sheet unlocks the topology, and then most of the energy is released through the outflow jets (\(\mathbf{v}\cdot\mathbf{J}\times\mathbf{B}\)), not by ohmic heating.
Step 8a: the vector-potential form¶
Write \(\mathbf{B} = \nabla\times\mathbf{A}\) and \(\mathbf{E} = -\partial_t\mathbf{A} - \nabla\phi\) (Vector and scalar potentials). Then \(\nabla\cdot\mathbf{B} = 0\) holds automatically, and Faraday's law is satisfied identically. We are free to choose the gauge; the Weyl gauge \(\phi = 0\) turns Ohm's law directly into an evolution equation for \(\mathbf{A}\):
This is the form the Simulation Lab code uses: it stores \(\mathbf{A}\), computes \(\mathbf{B} = \nabla\times\mathbf{A}\) whenever it needs the field, and so never creates magnetic monopoles. It also works for non-uniform \(\eta\) (for example, a resistivity enhanced near the X-line), because \(\mathbf{E} = -\mathbf{v}\times\mathbf{B} + \mu_0\eta\mathbf{J}\) holds point by point.
In 2D (nothing depends on \(z\)), the in-plane field is \(\mathbf{B}_\perp = \nabla A_z\times\hat{\mathbf{z}}\) and the \(z\) component of the equation for \(\mathbf{A}\) becomes (Problem 8)
\(A_z\) is the flux function ψ: with \(L = x\), \(N = y\) and \(M = -z\), \(\psi = -A_M = A_z\). Its contours are field lines, and the change of \(A_z\) at the X-line is the reconnected flux.
Step 8b: dimensionless form¶
Pick a length \(L_0\), a field \(B_0\) and a density \(\rho_0\). Then the natural speed is the Alfvén speed and the natural time is the Alfvén crossing time:
Measuring every quantity in these units (\(\tilde{\mathbf{x}} = \mathbf{x}/L_0\), \(\tilde t = t/\tau_A\), \(\tilde{\mathbf{B}} = \mathbf{B}/B_0\), …) removes \(\mu_0\) and leaves dimensionless numbers in place of \(\eta\) and \(\nu\):
So in code units, \(\eta\) is \(1/S\) and \(\nu\) is \(1/Re\). Two more numbers set the physics: \(\gamma\) and the plasma beta \(\beta = 2\mu_0 p/B^2\). Any two systems with the same \(S\), \(Re\), \(\beta\) and \(\gamma\) behave identically after rescaling: a lab experiment and a solar flare differ "only" in \(S\).
From a real system to code parameters
Solar corona: \(S\sim10^{13}\). No computer can resolve a current sheet with \(\delta/L = S^{-1/2}\sim3\times10^{-7}\), so simulations use \(S\sim10^3\)–\(10^5\) (\(\eta = 10^{-3}\)–\(10^{-5}\) in code units) and study how results scale with \(S\). Above \(S\approx10^4\) the Sweet–Parker sheet itself goes unstable to plasmoids, and the rate becomes nearly independent of \(S\) (about 0.01).
Practice problems¶
Try each problem before opening its solution. Click an option to check a concept question;
type a number and press Check for a numerical one (forms like 6.2e7 or 6.2x10^7 both work).
Concept questions¶
Problem 1
Why does the electric force \(\rho_q\mathbf{E}\) not appear in the one-fluid momentum equation?
- Because \(\mathbf{E} = 0\) in a plasma
- Because \(\mathbf{E}\) is parallel to \(\mathbf{B}\)
- Because quasi-neutrality makes the charge density negligible
- Because the displacement current was dropped
The force on the electrons and on the ions are \(\mp en\mathbf{E}\); with \(n_i = n_e\) they cancel. \(\mathbf{E}\) itself is not zero: it is set by Ohm's law.
Problem 2
Which approximation turns Ampère–Maxwell's law into \(\mu_0\mathbf{J} = \nabla\times\mathbf{B}\)?
- Quasi-neutrality
- Flow speeds much less than the speed of light
- Large Lundquist number
- Isotropic pressure
The displacement current is smaller than \(\nabla\times\mathbf{B}\) by \(v^2/c^2\).
Problem 3
A simulation is run with \(\eta = 2\times10^{-3}\) in code units. What Lundquist number does it represent?
- 500
- 0.002
- 2000
- \(\sqrt{500}\approx 22\)
In code units \(\eta = 1/S\), so \(S = 1/0.002 = 500\) (based on the unit length and the unit Alfvén speed).
Problem 4
Which energy-exchange term mainly powers the outflow jets of reconnection?
- Ohmic heating \(\mu_0\eta J^2\)
- \(\mathbf{v}\cdot(\mathbf{J}\times\mathbf{B})\), work done by the magnetic force
- \(p\nabla\cdot\mathbf{v}\), work done by expansion
- Viscous heating
The tension of the newly reconnected, sharply bent field lines does work on the plasma, \(\mathbf{v}\cdot(\mathbf{J}\times\mathbf{B}) > 0\), accelerating it to about the Alfvén speed.
Show that…¶
Problem 5: the energy exchange identity
Using \(\mathbf{E} = -\mathbf{v}\times\mathbf{B} + \mu_0\eta\mathbf{J}\), show that \(\mathbf{J}\cdot\mathbf{E} = \mathbf{v}\cdot(\mathbf{J}\times\mathbf{B}) + \mu_0\eta J^2\).
Hint
Use the cyclic property of the triple product: \(\mathbf{a}\cdot(\mathbf{b}\times\mathbf{c}) = \mathbf{b}\cdot(\mathbf{c}\times\mathbf{a})\).
Proof
\(\mathbf{J}\cdot\mathbf{E} = -\mathbf{J}\cdot(\mathbf{v}\times\mathbf{B}) + \mu_0\eta J^2\). By the cyclic rule \(\mathbf{J}\cdot(\mathbf{v}\times\mathbf{B}) = \mathbf{v}\cdot(\mathbf{B}\times\mathbf{J}) = -\mathbf{v}\cdot(\mathbf{J}\times\mathbf{B})\). So \(\mathbf{J}\cdot\mathbf{E} = \mathbf{v}\cdot(\mathbf{J}\times\mathbf{B}) + \mu_0\eta J^2\): the field gives energy to the flow (first term) and to heat (second term).
Problem 6: the induction equation keeps ∇·B = 0
Show that \(\partial_t\mathbf{B} = \nabla\times(\mathbf{v}\times\mathbf{B} - \mu_0\eta\mathbf{J})\) implies \(\partial_t(\nabla\cdot\mathbf{B}) = 0\), even for non-uniform \(\eta\).
Hint
The divergence of any curl is zero.
Proof
\(\partial_t(\nabla\cdot\mathbf{B}) = \nabla\cdot\partial_t\mathbf{B} = \nabla\cdot[\nabla\times(\ldots)] = 0\) for any vector inside the curl. Writing the induction equation as a curl, rather than expanding it into \(\eta\nabla^2\mathbf{B}\), keeps this property even when \(\eta\) varies in space.
Problem 7: adiabatic pressure equation
Starting from \(d(p\rho^{-\gamma})/dt = 0\) and the continuity equation, show that \(\partial_t p + \mathbf{v}\cdot\nabla p + \gamma p\nabla\cdot\mathbf{v} = 0\).
Hint
Expand the derivative of the product, then use \(d\rho/dt = -\rho\nabla\cdot\mathbf{v}\).
Proof
\(\rho^{-\gamma}\,dp/dt - \gamma p\rho^{-\gamma-1}\,d\rho/dt = 0\), so \(dp/dt = (\gamma p/\rho)\,d\rho/dt = -\gamma p\nabla\cdot\mathbf{v}\). With \(dp/dt = \partial_t p + \mathbf{v}\cdot\nabla p\) this is the result.
Problem 8: the 2D flux function is carried and diffused
In 2D with \(\mathbf{B} = \nabla\times(A_z\hat{\mathbf{z}})\) and uniform \(\eta\), show that the \(z\) component of \(\partial_t\mathbf{A} = \mathbf{v}\times\mathbf{B} - \mu_0\eta\mathbf{J}\) gives \(\partial_t A_z + \mathbf{v}\cdot\nabla A_z = \eta\nabla^2 A_z\).
Hint
\(B_x = \partial_y A_z\), \(B_y = -\partial_x A_z\), and \(\mu_0 J_z = \partial_x B_y - \partial_y B_x\).
Proof
\((\mathbf{v}\times\mathbf{B})_z = v_xB_y - v_yB_x = -v_x\partial_xA_z - v_y\partial_yA_z = -\mathbf{v}\cdot\nabla A_z\). \(\mu_0J_z = -\partial_x^2A_z - \partial_y^2A_z = -\nabla^2A_z\), so \(-\mu_0\eta J_z = \eta\nabla^2 A_z\).
Numbers¶
Problem 9
What is the Spitzer resistivity of a 100 eV plasma with \(\ln\Lambda = 20\), in Ω·m?
Solution
\(\eta_{\text{res}} \approx 5.2\times10^{-5}\times20\times100^{-3/2} = 1.04\times10^{-3}/1000 = 1.04\times10^{-6}\) Ω·m.
Problem 10
Convert that resistivity into a magnetic diffusivity \(\eta = \eta_{\text{res}}/\mu_0\), in m²/s.
Solution
\(\eta = 1.04\times10^{-6}/(1.257\times10^{-6}) \approx 0.83\) m²/s.
Problem 11
The ion inertial length is \(d_i \approx 228\ \text{km}/\sqrt{n\,[\text{cm}^{-3}]}\). In the magnetotail \(n = 0.3\) cm\(^{-3}\). What is \(d_i\), in km?
Solution
\(d_i = 228/\sqrt{0.3} \approx 416\) km. A tail current sheet is often only a few \(d_i\) thick, so the Hall term matters there and resistive MHD can describe only the large-scale flows.
Problem 12
A code-unit current sheet has \(J = 3\) and \(\eta = 5\times10^{-3}\). What is the ohmic heating rate \(\eta J^2\) in code units?
Solution
\(\eta J^2 = 5\times10^{-3}\times9 = 0.045\) (in units of \(B_0^2/(\mu_0\tau_A)\)). In code units \(\mu_0 = 1\).