Foundation | Prerequisites

Mathematical Prerequisites for Solar Physics

The essential mathematics, physics, and computational tools you need before diving into the solar physics course

1. Vector Calculus Review

Solar physics is fundamentally a theory of vector fields — magnetic fields, velocity fields, and force densities. Fluency in vector calculus, especially in curvilinear coordinates, is non-negotiable.

1.1 Gradient, Divergence, and Curl

Gradient of a scalar field $\phi$ in Cartesian coordinates:

$$\nabla\phi = \frac{\partial \phi}{\partial x}\hat{\mathbf{x}} + \frac{\partial \phi}{\partial y}\hat{\mathbf{y}} + \frac{\partial \phi}{\partial z}\hat{\mathbf{z}}$$

In spherical coordinates $(r, \theta, \phi)$:

$$\nabla\phi = \frac{\partial \phi}{\partial r}\hat{\mathbf{r}} + \frac{1}{r}\frac{\partial \phi}{\partial \theta}\hat{\boldsymbol{\theta}} + \frac{1}{r\sin\theta}\frac{\partial \phi}{\partial \varphi}\hat{\boldsymbol{\varphi}}$$

Divergence in Cartesian: $\nabla \cdot \mathbf{F} = \partial F_x/\partial x + \partial F_y/\partial y + \partial F_z/\partial z$. In spherical coordinates:

$$\nabla \cdot \mathbf{F} = \frac{1}{r^2}\frac{\partial(r^2 F_r)}{\partial r} + \frac{1}{r\sin\theta}\frac{\partial(\sin\theta\, F_\theta)}{\partial \theta} + \frac{1}{r\sin\theta}\frac{\partial F_\varphi}{\partial \varphi}$$

Curl in Cartesian: $(\nabla\times\mathbf{F})_i = \epsilon_{ijk}\partial_j F_k$. In spherical coordinates:

$$(\nabla\times\mathbf{F})_r = \frac{1}{r\sin\theta}\left[\frac{\partial(\sin\theta\, F_\varphi)}{\partial \theta} - \frac{\partial F_\theta}{\partial \varphi}\right]$$

$$(\nabla\times\mathbf{F})_\theta = \frac{1}{r}\left[\frac{1}{\sin\theta}\frac{\partial F_r}{\partial \varphi} - \frac{\partial(r F_\varphi)}{\partial r}\right]$$

$$(\nabla\times\mathbf{F})_\varphi = \frac{1}{r}\left[\frac{\partial(r F_\theta)}{\partial r} - \frac{\partial F_r}{\partial \theta}\right]$$

1.2 Derivation: Divergence-Free Magnetic Field in Spherical Coordinates

Maxwell's equations require that the magnetic field is solenoidal: $\nabla \cdot \mathbf{B} = 0$. In spherical coordinates, this constraint reads:

Step 1. Write out the spherical divergence applied to $\mathbf{B} = B_r\hat{\mathbf{r}} + B_\theta\hat{\boldsymbol{\theta}} + B_\varphi\hat{\boldsymbol{\varphi}}$:

$$\boxed{\nabla \cdot \mathbf{B} = \frac{1}{r^2}\frac{\partial(r^2 B_r)}{\partial r} + \frac{1}{r\sin\theta}\frac{\partial(\sin\theta\, B_\theta)}{\partial \theta} + \frac{1}{r\sin\theta}\frac{\partial B_\varphi}{\partial \varphi} = 0}$$

Step 2. For a purely radial field (e.g., a monopole), $B_\theta = B_\varphi = 0$, so the constraint reduces to $\partial(r^2 B_r)/\partial r = 0$, giving$B_r \propto 1/r^2$. This is why magnetic monopoles, if they existed, would produce inverse-square fields.

Step 3. For a dipole field, we write $B_r = B_0(R_\odot/r)^3 \cdot 2\cos\theta$and $B_\theta = B_0(R_\odot/r)^3 \cdot \sin\theta$. One can verify:

$$\frac{1}{r^2}\frac{\partial}{\partial r}\left(r^2 \cdot \frac{2B_0 R_\odot^3\cos\theta}{r^3}\right) = \frac{-6B_0 R_\odot^3\cos\theta}{r^4}$$

$$\frac{1}{r\sin\theta}\frac{\partial}{\partial\theta}\left(\sin\theta \cdot \frac{B_0 R_\odot^3\sin\theta}{r^3}\right) = \frac{2B_0 R_\odot^3\cos\theta \cdot 1}{r^4} \cdot \frac{r}{r\sin\theta}\cdot\sin\theta$$

Computing carefully: $\partial(\sin^2\theta)/\partial\theta = 2\sin\theta\cos\theta$, so the $\theta$-term gives $+2B_0 R_\odot^3\cos\theta/(r^4) \times 3 = +6B_0 R_\odot^3\cos\theta/r^4$. The two terms cancel exactly, confirming $\nabla\cdot\mathbf{B}=0$.

1.3 Integral Theorems

Gauss's (Divergence) Theorem: For a volume $V$ bounded by surface $S$:

$$\oint_S \mathbf{F} \cdot d\mathbf{A} = \int_V \nabla \cdot \mathbf{F} \, dV$$

Applied to $\nabla\cdot\mathbf{B}=0$: the total magnetic flux through any closed surface vanishes. This means every field line that enters a volume must also leave it — there are no magnetic charges.

Stokes's Theorem: For a surface $S$ bounded by curve $C$:

$$\oint_C \mathbf{F} \cdot d\boldsymbol{\ell} = \int_S (\nabla \times \mathbf{F}) \cdot d\mathbf{A}$$

Applied to Ampere's law: $\oint \mathbf{B}\cdot d\boldsymbol{\ell} = \mu_0 I_{\mathrm{enc}}$. This is how we compute the current flowing through a coronal loop from the magnetic field measured at its boundary.

1.4 Key Vector Identity for MHD

The induction equation requires the identity for the curl of a cross product. Starting from the general vector identity:

$$\boxed{\nabla\times(\mathbf{v}\times\mathbf{B}) = (\mathbf{B}\cdot\nabla)\mathbf{v} - (\mathbf{v}\cdot\nabla)\mathbf{B} + \mathbf{v}(\nabla\cdot\mathbf{B}) - \mathbf{B}(\nabla\cdot\mathbf{v})}$$

Derivation: Use the BAC-CAB rule: $\mathbf{A}\times(\mathbf{B}\times\mathbf{C}) = \mathbf{B}(\mathbf{A}\cdot\mathbf{C}) - \mathbf{C}(\mathbf{A}\cdot\mathbf{B})$combined with the product rule for the $\nabla$ operator. Since $\nabla\cdot\mathbf{B}=0$in MHD, the third term vanishes, giving:

$$\nabla\times(\mathbf{v}\times\mathbf{B}) = (\mathbf{B}\cdot\nabla)\mathbf{v} - (\mathbf{v}\cdot\nabla)\mathbf{B} - \mathbf{B}(\nabla\cdot\mathbf{v})$$

The three surviving terms represent: (1) stretching of field lines by velocity gradients along $\mathbf{B}$, (2) advection of $\mathbf{B}$ by the flow, and (3) compression/expansion of field lines.

Where this appears in the course:

Ch. 9 — Magnetic Fields of the Sun: spherical harmonic decomposition of the photospheric field
Ch. 10 — Solar Dynamo: the induction equation uses this identity throughout
Ch. 12 — Solar Flares: curl of B determines current sheets and reconnection sites

2. Differential Equations

Nearly every quantitative result in solar physics comes from solving a differential equation. Here we review the types that appear most frequently.

2.1 First-Order ODEs

Separation of variables: If $dy/dx = f(x)g(y)$, then$\int dy/g(y) = \int f(x)\,dx$. Example: radioactive decay$dN/dt = -\lambda N$ gives $N(t) = N_0 e^{-\lambda t}$.

Integrating factor: For $dy/dx + P(x)y = Q(x)$, multiply by $\mu(x) = e^{\int P\,dx}$:

$$\frac{d}{dx}\left[\mu(x) y\right] = \mu(x) Q(x) \implies y = \frac{1}{\mu(x)}\int \mu(x) Q(x)\,dx$$

2.2 Second-Order ODEs

Harmonic oscillator: The equation $d^2y/dx^2 + \omega^2 y = 0$has solutions $y = A\cos(\omega x) + B\sin(\omega x)$. This governs all wave phenomena in the Sun — acoustic waves, MHD waves, and gravity waves.

Bessel's equation: In cylindrical or spherical geometries, separation of variables leads to:

$$x^2 \frac{d^2y}{dx^2} + x\frac{dy}{dx} + (x^2 - n^2)y = 0$$

with solutions $J_n(x)$ (Bessel functions of the first kind). In helioseismology, the radial eigenfunctions of solar oscillations involve Bessel functions. The zeros of$J_n$ determine the allowed eigenfrequencies.

2.3 The Parker Wind Equation

A central result of solar wind theory is Parker's (1958) equation for a steady, spherically symmetric, isothermal wind. Starting from the momentum equation and continuity:

Step 1. Mass conservation: $\rho v r^2 = \text{const}$, so$d(\rho v r^2)/dr = 0$, giving:

$$\frac{1}{\rho}\frac{d\rho}{dr} = -\frac{1}{v}\frac{dv}{dr} - \frac{2}{r}$$

Step 2. The momentum equation for isothermal flow ($P = \rho c_s^2$):

$$v\frac{dv}{dr} = -\frac{c_s^2}{\rho}\frac{d\rho}{dr} - \frac{GM_\odot}{r^2}$$

Step 3. Substitute the density gradient from Step 1:

$$\boxed{\left(v - \frac{c_s^2}{v}\right)\frac{dv}{dr} = \frac{2c_s^2}{r} - \frac{GM_\odot}{r^2}}$$

Critical point analysis: The right-hand side vanishes at the critical radius$r_c = GM_\odot/(2c_s^2)$. The left-hand side vanishes when $v = c_s$(the sonic point). The unique transonic solution that passes smoothly through $(r_c, c_s)$is the Parker solar wind solution. All other solution branches are either subsonic everywhere (breezes) or unphysical.

2.4 The Diffusion Equation

Heat conduction in the solar corona is governed by the diffusion equation:

$$\boxed{\frac{\partial T}{\partial t} = \kappa \nabla^2 T}$$

where $\kappa$ is the thermal diffusivity. For Spitzer conductivity in a fully ionized plasma, $\kappa \propto T^{5/2}$, making this equation nonlinear. The characteristic diffusion time across a length $L$ is $\tau_{\text{diff}} \sim L^2/\kappa$.

Similarly, magnetic diffusion obeys $\partial\mathbf{B}/\partial t = \eta\nabla^2\mathbf{B}$where $\eta = 1/(\mu_0\sigma)$ is the magnetic diffusivity. The magnetic diffusion time $\tau_\eta = L^2/\eta$ is enormous for the Sun (~10$^{10}$ years), explaining why solar magnetic fields are effectively frozen into the plasma.

Where this appears in the course:

Ch. 3 — Helioseismology: Bessel functions describe radial oscillation eigenmodes
Ch. 8 — Solar Wind Physics: Parker wind equation and critical point analysis
Ch. 7 — Corona: heat conduction determines coronal loop temperature profiles

3. Electromagnetism Essentials

The Sun is a magnetic star. Its activity, from sunspots to flares to the solar wind, is driven by electromagnetic forces. We need Maxwell's equations and their consequences.

3.1 Maxwell's Equations in SI

$$\nabla \cdot \mathbf{E} = \frac{\rho_q}{\varepsilon_0}, \qquad \nabla \cdot \mathbf{B} = 0$$

$$\nabla \times \mathbf{E} = -\frac{\partial \mathbf{B}}{\partial t}, \qquad \nabla \times \mathbf{B} = \mu_0\mathbf{J} + \mu_0\varepsilon_0\frac{\partial \mathbf{E}}{\partial t}$$

3.2 Derivation: The Electromagnetic Wave Equation

Step 1. Take the curl of Faraday's law:$\nabla\times(\nabla\times\mathbf{E}) = -\partial(\nabla\times\mathbf{B})/\partial t$.

Step 2. Use the vector identity $\nabla\times(\nabla\times\mathbf{E}) = \nabla(\nabla\cdot\mathbf{E}) - \nabla^2\mathbf{E}$. In vacuum, $\nabla\cdot\mathbf{E} = 0$.

Step 3. Substitute Ampere-Maxwell in vacuum ($\mathbf{J}=0$):$\nabla\times\mathbf{B} = \mu_0\varepsilon_0\partial\mathbf{E}/\partial t$:

$$\boxed{\nabla^2\mathbf{E} = \mu_0\varepsilon_0\frac{\partial^2\mathbf{E}}{\partial t^2} = \frac{1}{c^2}\frac{\partial^2\mathbf{E}}{\partial t^2}}$$

where $c = 1/\sqrt{\mu_0\varepsilon_0}$. This demonstrates that electromagnetic disturbances propagate as waves at the speed of light.

3.3 Magnetic Pressure and Tension

The magnetic force per unit volume can be decomposed using the identity:

$$\frac{1}{\mu_0}(\nabla\times\mathbf{B})\times\mathbf{B} = -\nabla\left(\frac{B^2}{2\mu_0}\right) + \frac{1}{\mu_0}(\mathbf{B}\cdot\nabla)\mathbf{B}$$

The first term is the magnetic pressure gradient: $P_B = B^2/(2\mu_0)$acts isotropically, like gas pressure. For a 100 G sunspot field,$P_B \approx 4 \times 10^3$ Pa, comparable to the photospheric gas pressure.

The second term is the magnetic tension: $(\mathbf{B}\cdot\nabla)\mathbf{B}/\mu_0$acts like a restoring force on bent field lines, analogous to tension on a rubber band. This is the force that drives Alfven waves.

3.4 Magnetic Reynolds Number and Frozen-in Flux

The induction equation combines Faraday's law with Ohm's law ($\mathbf{J} = \sigma(\mathbf{E} + \mathbf{v}\times\mathbf{B})$):

$$\frac{\partial \mathbf{B}}{\partial t} = \nabla\times(\mathbf{v}\times\mathbf{B}) + \eta\nabla^2\mathbf{B}$$

The ratio of the advection term to the diffusion term defines the magnetic Reynolds number:

$$\boxed{R_m = \frac{vL}{\eta}}$$

For $R_m \gg 1$ (typical in the solar corona: $R_m \sim 10^{10}$), diffusion is negligible and the field is "frozen" into the plasma. Flux conservation means $d\Phi_B/dt = 0$ for any surface comoving with the fluid. This is Alfven's frozen-in theorem.

3.5 The Lorentz Force

A charged particle with charge $q$ and velocity $\mathbf{v}$ in electric and magnetic fields experiences:

$$\boxed{\mathbf{F} = q(\mathbf{E} + \mathbf{v}\times\mathbf{B})}$$

Since $\mathbf{v}\times\mathbf{B}$ is perpendicular to $\mathbf{v}$, the magnetic force does no work. It deflects particles into helical orbits around field lines, with gyrofrequency $\omega_c = |q|B/m$ and gyroradius $r_L = mv_\perp/(|q|B)$.

Where this appears in the course:

Ch. 9 — Magnetic Fields: magnetic pressure balances gas pressure in sunspots
Ch. 12 — Solar Flares: reconnection occurs where frozen-in breaks down ($R_m \sim 1$)
Ch. 14 — Energetic Particles: Lorentz force governs particle acceleration and transport

4. Fluid Mechanics

The solar interior, atmosphere, and wind are treated as fluids. Understanding the equations of fluid dynamics is essential before adding magnetic fields to get MHD.

4.1 The Euler Equation

For an inviscid fluid, Newton's second law per unit volume gives the Euler equation:

$$\boxed{\rho\left(\frac{\partial \mathbf{v}}{\partial t} + (\mathbf{v}\cdot\nabla)\mathbf{v}\right) = -\nabla P + \rho\mathbf{g}}$$

The left-hand side is $\rho$ times the material derivative $D\mathbf{v}/Dt$. The term $(\mathbf{v}\cdot\nabla)\mathbf{v}$ is the advective acceleration, responsible for nonlinear effects like shock formation in the solar wind.

4.2 Navier-Stokes with Viscosity

Adding viscous stress to Euler's equation gives the Navier-Stokes equation:

$$\rho\frac{D\mathbf{v}}{Dt} = -\nabla P + \rho\mathbf{g} + \mu\nabla^2\mathbf{v} + \left(\zeta + \frac{\mu}{3}\right)\nabla(\nabla\cdot\mathbf{v})$$

where $\mu$ is the dynamic viscosity and $\zeta$ is the bulk viscosity. For an incompressible fluid ($\nabla\cdot\mathbf{v}=0$), the last term vanishes. In solar physics, viscous effects are usually small compared to magnetic forces, but they matter in the photospheric boundary layer and in certain models of coronal heating.

4.3 Bernoulli's Equation

For steady, incompressible, inviscid flow along a streamline:

$$\boxed{\frac{1}{2}v^2 + \frac{P}{\rho} + gz = \text{const along a streamline}}$$

In solar physics, a generalized Bernoulli equation (including enthalpy and magnetic terms) is used to analyze solar wind acceleration. The Parker wind solution is essentially the Bernoulli integral for an isothermal, spherically symmetric flow with gravity.

4.4 Reynolds Number and Turbulence

The Reynolds number compares inertial to viscous forces:

$$Re = \frac{\rho v L}{\mu} = \frac{vL}{\nu}$$

where $\nu = \mu/\rho$ is the kinematic viscosity. Turbulence onset typically occurs for $Re \gtrsim 10^3$. In the solar convection zone, $Re \sim 10^{12}$, meaning the flow is wildly turbulent. This is why granulation and supergranulation are observed.

Where this appears in the course:

Ch. 5 — Photosphere: granulation modeled as convective fluid flow
Ch. 8 — Solar Wind: Bernoulli principle extended for accelerating wind
Ch. 1 — Solar Interior: convection zone dynamics driven by high Reynolds number turbulence

5. Thermodynamics & Statistical Mechanics

The Sun spans an enormous range of temperatures (6000 K surface to 15 million K core) and densities. Thermodynamic and statistical methods are essential for understanding its radiation, ionization state, and energy transport.

5.1 Ideal Gas Law for Plasmas

The equation of state for an ideal gas:

$$\boxed{P = nkT = \frac{\rho k T}{\mu m_H}}$$

where $n$ is the number density, $k$ is Boltzmann's constant,$\mu$ is the mean molecular weight (which depends on ionization state!), and$m_H$ is the hydrogen mass. For a fully ionized hydrogen-helium plasma with solar abundances ($X = 0.7, Y = 0.28$), $\mu \approx 0.6$.

5.2 The Saha Equation

The ionization equilibrium between stage $r$ and $r+1$ of an element is given by the Saha equation:

$$\boxed{\frac{n_{r+1} n_e}{n_r} = \frac{2U_{r+1}}{U_r}\left(\frac{2\pi m_e k T}{h^2}\right)^{3/2} \exp\left(-\frac{\chi_r}{kT}\right)}$$

where $U_r$ is the partition function and $\chi_r$ is the ionization potential. At the solar photosphere ($T \approx 5800$ K), hydrogen is mostly neutral, but metals like Na and Ca have their first ionization stages populated, giving rise to prominent absorption lines.

5.3 Boltzmann Distribution

The population ratio between energy levels $i$ and $0$ (ground state) is:

$$\boxed{\frac{n_i}{n_0} = \frac{g_i}{g_0}\exp\left(-\frac{E_i}{kT}\right)}$$

where $g_i$ is the statistical weight (degeneracy) of level $i$. This determines the strength of absorption and emission lines. For the H$\alpha$line at $T = 10{,}000$ K, the ratio of $n=2$ to $n=1$population is $\sim 10^{-4}$, yet this is enough to produce the strong chromospheric emission seen in solar observations.

5.4 Radiation: Planck Function and Stefan-Boltzmann Law

The spectral radiance of a blackbody at temperature $T$:

$$B_\nu(T) = \frac{2h\nu^3}{c^2}\frac{1}{e^{h\nu/(kT)} - 1}$$

The peak shifts according to Wien's law: $\lambda_{\max} T = 2.898 \times 10^{-3}$ m K. For $T_\odot \approx 5778$ K, $\lambda_{\max} \approx 502$ nm (green-yellow).

Integrating over all frequencies gives the Stefan-Boltzmann law for total luminosity:

$$\boxed{L = 4\pi R^2 \sigma T^4}$$

For the Sun: $L_\odot = 4\pi(6.96\times 10^8)^2 \times 5.67\times 10^{-8} \times 5778^4 \approx 3.85 \times 10^{26}$ W.

Where this appears in the course:

Ch. 5 — Photosphere: Saha and Boltzmann determine spectral line formation
Ch. 6 — Chromosphere: departure from LTE requires non-equilibrium treatment
Ch. 1 — Solar Interior: Stefan-Boltzmann determines the photospheric boundary condition

6. Plasma Physics Basics

The Sun is made of plasma from its core to the outer heliosphere. Plasma has collective behavior that sets it apart from neutral gases. Here are the key scales and parameters.

6.1 Debye Shielding

A test charge in a plasma is screened by a cloud of oppositely charged particles. The characteristic screening length is the Debye length:

$$\boxed{\lambda_D = \sqrt{\frac{\varepsilon_0 k T}{n_e e^2}}}$$

Derivation: Poisson's equation $\nabla^2\phi = -\rho_q/\varepsilon_0$combined with Boltzmann-distributed electrons $n_e = n_0 e^{e\phi/(kT)}$. Linearizing for $e\phi \ll kT$: $\nabla^2\phi = \phi/\lambda_D^2$, giving exponentially decaying potential $\phi \propto e^{-r/\lambda_D}/r$.

In the solar corona ($T = 10^6$ K, $n_e = 10^{14}$ m$^{-3}$):$\lambda_D \approx 0.07$ m. The number of particles in a Debye sphere,$N_D = (4/3)\pi n_e \lambda_D^3 \sim 10^9 \gg 1$, confirming that collective plasma behavior dominates over individual particle interactions.

6.2 Plasma Frequency

The natural oscillation frequency of electron density perturbations:

$$\boxed{\omega_{pe} = \sqrt{\frac{n_e e^2}{\varepsilon_0 m_e}}}$$

Derivation: Displace electrons by $\delta x$. The restoring electric field is $E = n_e e \delta x / \varepsilon_0$. Newton's law for an electron:$m_e \ddot{\delta x} = -eE = -n_e e^2 \delta x/\varepsilon_0$. This is a harmonic oscillator with $\omega^2 = n_e e^2/(\varepsilon_0 m_e)$.

Electromagnetic waves can only propagate in the plasma if $\omega > \omega_{pe}$. This is why radio waves below a critical frequency cannot escape from the denser solar atmosphere — and why solar radio bursts reveal information about coronal density.

6.3 Gyrofrequency and Gyroradius

A charged particle in a uniform magnetic field orbits with:

$$\omega_c = \frac{|q|B}{m}, \qquad r_L = \frac{mv_\perp}{|q|B}$$

The gyrofrequency $\omega_c$ (or cyclotron frequency) sets the fastest timescale for particle dynamics in a magnetic field. The gyroradius $r_L$ (or Larmor radius) sets the smallest spatial scale. For a thermal proton in the corona ($T = 10^6$ K,$B = 10$ G): $r_L \approx 0.1$ m, $\omega_c \approx 10^6$ rad/s.

6.4 Alfven Speed

The characteristic speed for magnetic disturbances propagating along field lines:

$$\boxed{v_A = \frac{B}{\sqrt{\mu_0 \rho}}}$$

Derivation: Consider a small transverse perturbation $\delta\mathbf{v}$of a uniform plasma with density $\rho_0$ and field $\mathbf{B}_0$. The linearized momentum equation gives $\rho_0\partial\delta\mathbf{v}/\partial t = (\mathbf{B}_0\cdot\nabla)\delta\mathbf{B}/\mu_0$and the induction equation gives $\partial\delta\mathbf{B}/\partial t = (\mathbf{B}_0\cdot\nabla)\delta\mathbf{v}$. Combining: $\partial^2\delta\mathbf{v}/\partial t^2 = v_A^2 (\hat{b}\cdot\nabla)^2\delta\mathbf{v}$, a wave equation with phase speed $v_A$ along $\mathbf{B}_0$.

In the corona ($B \sim 10$ G, $n \sim 10^{14}$ m$^{-3}$):$v_A \sim 2000$ km/s, comparable to the coronal sound speed.

6.5 The MHD Equations

The full set of ideal MHD equations combines fluid dynamics with Maxwell's equations:

Continuity:

$$\frac{\partial \rho}{\partial t} + \nabla\cdot(\rho\mathbf{v}) = 0$$

Momentum:

$$\rho\frac{D\mathbf{v}}{Dt} = -\nabla P + \frac{1}{\mu_0}(\nabla\times\mathbf{B})\times\mathbf{B} + \rho\mathbf{g}$$

Induction:

$$\frac{\partial\mathbf{B}}{\partial t} = \nabla\times(\mathbf{v}\times\mathbf{B})$$

Energy:

$$\frac{d}{dt}\left(\frac{P}{\rho^\gamma}\right) = 0 \quad \text{(adiabatic)}$$

Constraint:

$$\nabla\cdot\mathbf{B} = 0$$

These eight equations (in conservative form) plus the solenoidal constraint form the foundation of essentially all theoretical and numerical solar physics.

Where this appears in the course:

Ch. 9 — Magnetic Fields: Alfven speed determines coronal loop oscillation periods
Ch. 7 — Corona: plasma beta ($\beta = P/(B^2/2\mu_0)$) classification of coronal regions
Ch. 12 — Solar Flares: full MHD simulations of reconnection and eruptions
Ch. 13 — CMEs: MHD models of flux rope eruption

7. Nuclear Physics Basics

The Sun's energy source is thermonuclear fusion. Understanding the relevant nuclear physics requires concepts from quantum mechanics that govern the reaction rates.

7.1 Binding Energy and Mass Defect

The mass of a nucleus is less than the sum of its constituent nucleon masses. The difference, called the mass defect $\Delta m$, corresponds to the binding energy via Einstein's relation:

$$\boxed{E_B = \Delta m \cdot c^2 = \left[Z m_p + (A-Z) m_n - M_{\text{nucleus}}\right] c^2}$$

For the pp-chain: $4p \to {}^4\text{He} + 2e^+ + 2\nu_e$. The mass defect is$\Delta m = 4(1.00728) - 4.00260 = 0.02652$ u, giving$E = 0.02652 \times 931.5 = 24.7$ MeV per helium nucleus produced. After subtracting the neutrino losses (~2%), this corresponds to $\epsilon \approx 0.7\%$mass-to-energy conversion efficiency.

7.2 Coulomb Barrier and Quantum Tunneling

Two nuclei with charges $Z_1, Z_2$ must overcome their mutual Coulomb repulsion to get close enough for the strong force to act. The Coulomb barrier height is:

$$E_C = \frac{Z_1 Z_2 e^2}{4\pi\varepsilon_0 r_0} \approx \frac{Z_1 Z_2 \times 1.44 \text{ MeV fm}}{r_0}$$

For pp fusion, $E_C \sim 550$ keV, yet the core temperature gives thermal energy$kT \approx 1.3$ keV. Classically, fusion is impossible! The resolution is quantum tunneling through the Coulomb barrier. The tunneling probability is characterized by the Gamow factor:

$$\boxed{T \propto \exp(-2\pi\eta), \qquad \eta = \frac{Z_1 Z_2 e^2}{\hbar v} = \frac{Z_1 Z_2 e^2}{4\pi\varepsilon_0 \hbar v}}$$

where $\eta$ is the Sommerfeld parameter and $v$ is the relative velocity. The exponential sensitivity to velocity means only particles in the high-energy tail of the Maxwell-Boltzmann distribution contribute significantly to the reaction rate.

7.3 Cross-Section and Reaction Rate

The nuclear cross-section is conventionally written:

$$\sigma(E) = \frac{S(E)}{E}\exp\left(-\sqrt{\frac{E_G}{E}}\right)$$

where $S(E)$ is the astrophysical S-factor (a slowly varying function of energy) and $E_G = 2\mu(Z_1 Z_2 e^2 \pi/\hbar)^2$ is the Gamow energy. The thermally averaged reaction rate is:

$$\boxed{\langle\sigma v\rangle = \sqrt{\frac{8}{\pi \mu (kT)^3}}\int_0^\infty S(E)\exp\left(-\frac{E}{kT} - \sqrt{\frac{E_G}{E}}\right)dE}$$

The integrand peaks at the Gamow peak energy $E_0 = (E_G (kT)^2 / 4)^{1/3}$. For pp at $T = 1.5\times 10^7$ K: $E_0 \approx 6$ keV, with a width of about 6 keV — the Gamow window.

Where this appears in the course:

Ch. 2 — Nuclear Fusion: pp-chain and CNO cycle reaction rates depend on tunneling
Ch. 4 — Solar Neutrinos: neutrino flux is proportional to reaction rate
Ch. 1 — Solar Interior: energy generation rate $\epsilon \propto T^4$ for pp-chain

8. Dimensional Analysis

Dimensional analysis is a physicist's most powerful tool for rapid estimation. Before solving any equation exactly, you should always be able to estimate the answer from dimensions alone.

8.1 Worked Example: The Jeans Length

Problem: Find the critical length scale below which a self-gravitating gas cloud is supported by thermal pressure against gravitational collapse.

Relevant parameters: gravitational constant $G$ [m$^3$ kg$^{-1}$ s$^{-2}$], density $\rho$ [kg m$^{-3}$], sound speed $c_s$ [m s$^{-1}$].

Dimensional argument: We need a length. From $G\rho$ we get [s$^{-2}$] (this is the free-fall rate squared). The sound speed gives [m s$^{-1}$]. So the only combination with dimensions of length is:

$$\boxed{\lambda_J \sim \frac{c_s}{\sqrt{G\rho}}}$$

The exact result from a stability analysis includes a factor of $\sqrt{\pi}$:$\lambda_J = c_s\sqrt{\pi/(G\rho)}$. Dimensional analysis got us the right answer up to an order-unity constant.

8.2 Worked Example: Alfven Crossing Time

Problem: How long does an Alfven wave take to cross a coronal loop of length $L$?

Parameters: loop length $L \sim 10^8$ m, magnetic field $B \sim 10^{-3}$ T, density $\rho \sim 10^{-12}$ kg/m$^3$.

$$v_A = \frac{B}{\sqrt{\mu_0\rho}} = \frac{10^{-3}}{\sqrt{4\pi\times 10^{-7} \times 10^{-12}}} \approx 9 \times 10^5 \text{ m/s}$$

$$\boxed{\tau_A = \frac{L}{v_A} \approx \frac{10^8}{9\times 10^5} \approx 110 \text{ s} \approx 2 \text{ minutes}}$$

This matches the observed periods of transverse coronal loop oscillations (2–10 minutes), as discovered by TRACE and SDO.

8.3 Worked Example: Magnetic Diffusion Time

Problem: How long does it take for a magnetic field to diffuse across a sunspot?

Parameters: sunspot radius $L \sim 10^7$ m, magnetic diffusivity$\eta = 1/(\mu_0\sigma) \sim 1$ m$^2$/s (photospheric plasma).

$$\boxed{\tau_\eta = \frac{L^2}{\eta} \approx \frac{(10^7)^2}{1} = 10^{14} \text{ s} \approx 3 \text{ million years}}$$

Yet sunspots live only days to weeks. This means sunspot decay is NOT governed by simple Ohmic diffusion. Turbulent diffusion (enhanced by convective motions) with an effective$\eta_{\text{turb}} \sim 10^8$ m$^2$/s gives a more realistic timescale of $\sim 10^6$ s $\sim$ 10 days.

Where this appears in the course:

Ch. 11 — Sunspots: magnetic diffusion time sets sunspot lifetimes
Ch. 10 — Solar Dynamo: diffusion time constrains the 11-year cycle period
Ch. 7 — Corona: Alfven crossing time sets coronal loop oscillation periods

Interactive Simulation: Mathematical Prerequisites in Action

This four-panel simulation visualizes key mathematical objects from the prerequisites. Panel 1 shows a magnetic dipole vector field (vector calculus). Panel 2 plots Bessel functions relevant to helioseismology modes. Panel 3 maps plasma parameters as functions of temperature and density. Panel 4 displays MHD wave dispersion relations.

Mathematical Prerequisites: Dipole Field, Bessel Functions, Plasma Parameters, MHD Dispersion

Python

script.py280 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from scipy.special import jn as Jn

# ============================================================
# Mathematical Prerequisites for Solar Physics: 4-Panel Simulation
# ============================================================

fig, axes = plt.subplots(2, 2, figsize=(14, 12))
fig.suptitle('Mathematical Prerequisites for Solar Physics', fontsize=16,
             fontweight='bold', color='white')
fig.patch.set_facecolor('#0a0a0a')

for ax in axes.flat:
    ax.set_facecolor('#0a0a0a')
    ax.tick_params(colors='white', labelsize=9)
    ax.xaxis.label.set_color('white')
    ax.yaxis.label.set_color('white')
    ax.title.set_color('white')
    for spine in ax.spines.values():
        spine.set_color('#444')

# ── Panel 1: Magnetic Dipole Vector Field ──
ax1 = axes[0, 0]

# Create grid in x-z plane (y=0 slice)
x = np.linspace(-3, 3, 20)
z = np.linspace(-3, 3, 20)
X, Z = np.meshgrid(x, z)
r = np.sqrt(X**2 + Z**2)
r = np.where(r < 0.5, np.nan, r)  # mask interior

# Dipole field: B_r = 2m cos(theta)/r^3, B_theta = m sin(theta)/r^3
cos_theta = Z / r
sin_theta = X / r

# In Cartesian: dipole with moment along z
Bx = 3.0 * X * Z / r**5
Bz = (3.0 * Z**2 / r**5) - 1.0 / r**3
Bmag = np.sqrt(Bx**2 + Bz**2)

# Normalize arrows for visibility
Bx_norm = Bx / (Bmag + 1e-10)
Bz_norm = Bz / (Bmag + 1e-10)

# Color by field strength
log_B = np.log10(Bmag + 1e-10)

quiv = ax1.quiver(X, Z, Bx_norm, Bz_norm, log_B,
                   cmap='hot', scale=25, width=0.004, alpha=0.9)

# Draw the "Sun" as a circle
theta_circle = np.linspace(0, 2*np.pi, 100)
ax1.fill(0.5*np.cos(theta_circle), 0.5*np.sin(theta_circle),
         color='#ff6b35', alpha=0.8, zorder=5)
ax1.plot(0.5*np.cos(theta_circle), 0.5*np.sin(theta_circle),
         color='#fbbf24', linewidth=1.5, zorder=6)

# Draw some field lines
for theta0 in np.linspace(0.15, np.pi/2, 6):
    r_line = np.linspace(0.5, 3.0, 200)
    # Dipole field line: r = r_eq * sin^2(theta)
    r_eq_val = 0.5 / np.sin(theta0)**2
    th_line = np.linspace(theta0, np.pi - theta0, 200)
    r_fl = r_eq_val * np.sin(th_line)**2
    mask = r_fl >= 0.5
    x_fl = r_fl[mask] * np.sin(th_line[mask])
    z_fl = r_fl[mask] * np.cos(th_line[mask])
    valid = r_fl[mask] < 3.0
    if np.any(valid):
        ax1.plot(x_fl[valid], z_fl[valid], color='#fbbf24', linewidth=0.8, alpha=0.5)

ax1.set_xlim(-3, 3)
ax1.set_ylim(-3, 3)
ax1.set_xlabel('x / R')
ax1.set_ylabel('z / R')
ax1.set_title('Magnetic Dipole Field (Vector Calculus)', fontsize=11)
ax1.set_aspect('equal')

# ── Panel 2: Bessel Functions (Helioseismology Modes) ──
ax2 = axes[0, 1]

x_bessel = np.linspace(0, 20, 500)
colors_bessel = ['#f43f5e', '#fb923c', '#fbbf24', '#34d399', '#60a5fa']
labels_bessel = ['J_0', 'J_1', 'J_2', 'J_3', 'J_4']

for n_order in range(5):
    y_bessel = Jn(n_order, x_bessel)
    ax2.plot(x_bessel, y_bessel, color=colors_bessel[n_order],
             linewidth=2, label=labels_bessel[n_order], alpha=0.9)
    # Mark zeros
    from scipy.special import jn_zeros
    zeros = jn_zeros(n_order, 5)
    zeros_in_range = zeros[zeros <= 20]
    ax2.scatter(zeros_in_range, np.zeros_like(zeros_in_range),
                color=colors_bessel[n_order], s=25, zorder=5, marker='o')

ax2.axhline(y=0, color='#555', linewidth=0.5, linestyle='--')
ax2.set_xlabel('x')
ax2.set_ylabel('J_n(x)')
ax2.set_title('Bessel Functions (Helioseismology Modes)', fontsize=11)
ax2.legend(fontsize=9, facecolor='#0a0a0a', edgecolor='#444', labelcolor='white',
           loc='upper right')
ax2.set_ylim(-0.5, 1.1)

# Annotate: "Zeros = eigenfrequencies"
ax2.annotate('Zeros determine\neigenfrequencies',
             xy=(2.4, 0), xytext=(5, 0.7),
             fontsize=8, color='#fbbf24',
             arrowprops=dict(arrowstyle='->', color='#fbbf24', lw=1.2))

# ── Panel 3: Plasma Parameters vs Temperature/Density ──
ax3 = axes[1, 0]

# Constants (SI)
eps0 = 8.854e-12
k_B = 1.381e-23
e_charge = 1.602e-19
m_e = 9.109e-31
m_p = 1.673e-27
mu0 = 4 * np.pi * 1e-7

T_range = np.logspace(4, 8, 200)  # 10^4 to 10^8 K
n_e_fixed = 1e15  # m^-3 (typical corona)

# Debye length
lambda_D = np.sqrt(eps0 * k_B * T_range / (n_e_fixed * e_charge**2))

# Plasma frequency (constant with T at fixed n_e)
omega_pe = np.sqrt(n_e_fixed * e_charge**2 / (eps0 * m_e))
omega_pe_arr = np.full_like(T_range, omega_pe)

# Thermal speed
v_th = np.sqrt(2 * k_B * T_range / m_p)

# Gyroradius (for B=10 G = 1e-3 T, proton)
B_field = 1e-3
r_L = m_p * v_th / (e_charge * B_field)

ax3_twin = ax3.twinx()

l1, = ax3.plot(T_range, lambda_D, color='#f43f5e', linewidth=2, label='Debye length')
l2, = ax3.plot(T_range, r_L, color='#fb923c', linewidth=2, label='Gyroradius (proton)')
l3, = ax3_twin.plot(T_range, v_th / 1e3, color='#60a5fa', linewidth=2, linestyle='--',
                     label='Thermal speed (km/s)')

ax3.set_xscale('log')
ax3.set_yscale('log')
ax3_twin.set_yscale('log')
ax3.set_xlabel('Temperature [K]')
ax3.set_ylabel('Length scale [m]', color='#f43f5e')
ax3_twin.set_ylabel('Speed [km/s]', color='#60a5fa')
ax3_twin.yaxis.label.set_color('#60a5fa')
ax3_twin.tick_params(colors='#60a5fa')
ax3.set_title('Plasma Parameters vs Temperature (n_e = 10^15 m^-3, B = 10 G)', fontsize=10)

# Shade solar regions
ax3.axvspan(5000, 7000, alpha=0.15, color='#fbbf24', label='Photosphere')
ax3.axvspan(1e6, 3e6, alpha=0.15, color='#f43f5e', label='Corona')

lines = [l1, l2, l3]
labels = [l.get_label() for l in lines]
ax3.legend(lines, labels, fontsize=8, facecolor='#0a0a0a', edgecolor='#444',
           labelcolor='white', loc='upper left')

# ── Panel 4: MHD Wave Dispersion Relations ──
ax4 = axes[1, 1]

# Parameters
v_A_val = 1000e3    # 1000 km/s (Alfven speed)
c_s_val = 200e3     # 200 km/s (sound speed)

theta_prop = np.linspace(0, np.pi, 300)

# Phase speeds for the three MHD modes
# v_fast^2 = 0.5*(vA^2 + cs^2) + 0.5*sqrt((vA^2+cs^2)^2 - 4*vA^2*cs^2*cos^2(theta))
# v_slow^2 = 0.5*(vA^2 + cs^2) - 0.5*sqrt(...)
# v_alfven = vA * |cos(theta)|

cos2 = np.cos(theta_prop)**2
sum_sq = v_A_val**2 + c_s_val**2
disc = np.sqrt(sum_sq**2 - 4 * v_A_val**2 * c_s_val**2 * cos2)

v_fast = np.sqrt(0.5 * (sum_sq + disc)) / 1e3  # km/s
v_slow = np.sqrt(np.maximum(0.5 * (sum_sq - disc), 0)) / 1e3
v_alfven = v_A_val * np.abs(np.cos(theta_prop)) / 1e3

# Polar plot of phase speeds
# Convert to Cartesian for polar diagram
vf_x = v_fast * np.sin(theta_prop)
vf_z = v_fast * np.cos(theta_prop)
vs_x = v_slow * np.sin(theta_prop)
vs_z = v_slow * np.cos(theta_prop)
va_x = v_alfven * np.sin(theta_prop)
va_z = v_alfven * np.cos(theta_prop)

# Also plot the negative theta half
ax4.plot(vf_x, vf_z, color='#f43f5e', linewidth=2, label='Fast magnetosonic')
ax4.plot(-vf_x, vf_z, color='#f43f5e', linewidth=2)
ax4.plot(vs_x, vs_z, color='#fbbf24', linewidth=2, label='Slow magnetosonic')
ax4.plot(-vs_x, vs_z, color='#fbbf24', linewidth=2)
ax4.plot(va_x, va_z, color='#60a5fa', linewidth=2, linestyle='--', label='Alfven')
ax4.plot(-va_x, va_z, color='#60a5fa', linewidth=2, linestyle='--')

# Draw B0 direction arrow
ax4.annotate('', xy=(0, 1100), xytext=(0, 900),
             arrowprops=dict(arrowstyle='->', color='white', lw=2))
ax4.text(50, 1050, 'B_0', color='white', fontsize=12, fontweight='bold')

ax4.set_xlabel('v_perp [km/s]')
ax4.set_ylabel('v_parallel [km/s]')
ax4.set_title('MHD Wave Phase Speed Diagram', fontsize=11)
ax4.set_aspect('equal')
ax4.legend(fontsize=9, facecolor='#0a0a0a', edgecolor='#444', labelcolor='white',
           loc='lower right')
ax4.axhline(y=0, color='#333', linewidth=0.5)
ax4.axvline(x=0, color='#333', linewidth=0.5)

plt.tight_layout(rect=[0, 0, 1, 0.95])
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0a0a0a')
plt.close()
print("Plot saved to output.png")

# ── Numerical output ──
print("\n" + "=" * 65)
print("MATHEMATICAL PREREQUISITES: KEY NUMERICAL VALUES")
print("=" * 65)

print("\n--- Plasma Parameters (Corona: T=10^6 K, n_e=10^15 m^-3) ---")
T_cor = 1e6
n_cor = 1e15
lD = np.sqrt(eps0 * k_B * T_cor / (n_cor * e_charge**2))
print("  Debye length = " + str(round(lD, 4)) + " m")
wpe = np.sqrt(n_cor * e_charge**2 / (eps0 * m_e))
print("  Plasma frequency = " + str(round(wpe/1e6, 1)) + " MHz (rad/s)")
fpe = wpe / (2*np.pi)
print("  Plasma frequency = " + str(round(fpe/1e6, 1)) + " MHz (Hz)")
vth_p = np.sqrt(2*k_B*T_cor / m_p)
print("  Proton thermal speed = " + str(round(vth_p/1e3, 1)) + " km/s")
B_cor = 1e-3
rL_p = m_p * vth_p / (e_charge * B_cor)
print("  Proton gyroradius (B=10 G) = " + str(round(rL_p, 3)) + " m")
omega_cp = e_charge * B_cor / m_p
print("  Proton gyrofrequency = " + str(round(omega_cp, 0)) + " rad/s")
vA = B_cor / np.sqrt(mu0 * n_cor * m_p)
print("  Alfven speed = " + str(round(vA/1e3, 0)) + " km/s")

print("\n--- Dimensional Analysis Results ---")
L_loop = 1e8
tau_A = L_loop / vA
print("  Alfven crossing time (L=100 Mm) = " + str(round(tau_A, 1)) + " s")
eta_phot = 1.0
L_spot = 1e7
tau_diff = L_spot**2 / eta_phot
print("  Ohmic diffusion time (sunspot) = " + str(round(tau_diff/3.15e7, 0)) + " years")
eta_turb = 1e8
tau_turb = L_spot**2 / eta_turb
print("  Turbulent diffusion time = " + str(round(tau_turb/86400, 1)) + " days")

print("\n--- Nuclear Physics ---")
delta_m = 4*1.00728 - 4.00260  # in atomic mass units
E_pp = delta_m * 931.5  # MeV
print("  pp-chain mass defect = " + str(round(delta_m, 5)) + " u")
print("  pp-chain energy release = " + str(round(E_pp, 1)) + " MeV")
print("  Mass-to-energy efficiency = " + str(round(delta_m/(4*1.00728)*100, 2)) + "%")
E_G_pp = 493.0  # keV, Gamow energy for pp
kT_core = 1.38e-23 * 1.5e7 / 1.602e-16  # keV
E0_gamow = (E_G_pp * kT_core**2 / 4)**(1.0/3.0)
print("  Core kT = " + str(round(kT_core, 2)) + " keV")
print("  Gamow peak energy (pp) = " + str(round(E0_gamow, 1)) + " keV")

print("\n--- MHD Wave Speeds ---")
print("  Alfven speed v_A = " + str(round(v_A_val/1e3, 0)) + " km/s")
print("  Sound speed c_s = " + str(round(c_s_val/1e3, 0)) + " km/s")
v_fast_par = np.sqrt(v_A_val**2 + c_s_val**2)
print("  Fast mode (perp to B) = " + str(round(v_fast_par/1e3, 0)) + " km/s")
print("  Slow mode (along B) = " + str(round(c_s_val/1e3, 0)) + " km/s")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Quick Reference: Key Equations

Topic	Equation	Where Used
Solenoidal condition	$\nabla\cdot\mathbf{B} = 0$	Everywhere (MHD constraint)
Parker wind	$(v - c_s^2/v) dv/dr = 2c_s^2/r - GM/r^2$	Ch. 8: Solar Wind
Mag. Reynolds #	$R_m = vL/\eta$	Ch. 9, 10, 12: Frozen flux
Euler equation	$\rho D\mathbf{v}/Dt = -\nabla P + \rho\mathbf{g}$	Ch. 1, 5, 8: Fluid dynamics
Ideal gas law	$P = \rho kT/(\mu m_H)$	Ch. 1, 5, 7, 8: Equation of state
Saha equation	$n_{r+1}n_e/n_r \propto T^{3/2}\exp(-\chi/kT)$	Ch. 5, 6: Ionization balance
Debye length	$\lambda_D = \sqrt{\varepsilon_0 kT/(n_e e^2)}$	Ch. 7, 14: Plasma scales
Alfven speed	$v_A = B/\sqrt{\mu_0\rho}$	Ch. 7, 9, 12, 13: MHD waves
Gamow tunneling	$T \propto \exp(-2\pi\eta)$	Ch. 2, 4: Nuclear reactions
Stefan-Boltzmann	$L = 4\pi R^2\sigma T^4$	Ch. 1, 5: Luminosity

← Course Overview Begin: Solar Interior →

Topic	Equation	Where Used
Solenoidal condition	\(\nabla\cdot\mathbf{B} = 0\)	Everywhere (MHD constraint)
Parker wind	\((v - c_s^2/v) dv/dr = 2c_s^2/r - GM/r^2\)	Ch. 8: Solar Wind
Mag. Reynolds #	\(R_m = vL/\eta\)	Ch. 9, 10, 12: Frozen flux
Euler equation	\(\rho D\mathbf{v}/Dt = -\nabla P + \rho\mathbf{g}\)	Ch. 1, 5, 8: Fluid dynamics
Ideal gas law	\(P = \rho kT/(\mu m_H)\)	Ch. 1, 5, 7, 8: Equation of state
Saha equation	\(n_{r+1}n_e/n_r \propto T^{3/2}\exp(-\chi/kT)\)	Ch. 5, 6: Ionization balance
Debye length	\(\lambda_D = \sqrt{\varepsilon_0 kT/(n_e e^2)}\)	Ch. 7, 14: Plasma scales
Alfven speed	\(v_A = B/\sqrt{\mu_0\rho}\)	Ch. 7, 9, 12, 13: MHD waves
Gamow tunneling	\(T \propto \exp(-2\pi\eta)\)	Ch. 2, 4: Nuclear reactions
Stefan-Boltzmann	\(L = 4\pi R^2\sigma T^4\)	Ch. 1, 5: Luminosity