plasma_plots.theory.orbits

Charged-particle orbits: gyromotion, guiding-center drifts and trapped particles in a tokamak.

The formulas hold in any consistent units (SI, or Struphy’s normalized units with charge and mass in units of e and m_p) unless a function says otherwise. Vectors are arrays with a trailing axis of length 3, (x, y, z) or any right-handed Cartesian basis; everything broadcasts.

Tokamak model

The trapped-particle functions use a large-aspect-ratio tokamak with circular flux surfaces, B(θ) = B₀ / (1 + ε cos θ) with ε = r/R₀ the inverse aspect ratio of the flux surface and θ the poloidal angle (θ = 0 on the outboard midplane, where B is smallest). A particle of speed v and magnetic moment μ has the pitch parameter λ = μB₀/E = (v⊥²/v²) B₀/B, and

v∥² = v² (1 + ε cos θ − λ) / (1 + ε cos θ) = 2ε v² (κ² − sin²(θ/2)) / (1 + ε cos θ)

with the trapping parameter κ² = (1 + ε − λ)/(2ε): particles with κ² < 1 are trapped (they bounce at sin²(θ/2) = κ²), κ² = 0 is deeply trapped at θ = 0, κ² > 1 passes. The bounce and transit frequencies and banana widths are to leading order in ε (the factor 1 + ε cos θ is dropped, and the field line length is q R₀ dθ).

References

P. Helander and D. J. Sigmar, Collisional Transport in Magnetized Plasmas (Cambridge, 2002), chapter 7.

J. Wesson, Tokamaks, 4th ed. (Oxford, 2011), section 3.12.

F. F. Chen, Introduction to Plasma Physics and Controlled Fusion, 3rd ed. (Springer, 2016), chapter 2.

Examples

The bounce frequency of a barely trapped (κ² = 0.9) particle with v = 1 on the ε = 0.1 surface of a tokamak with q = 2 and R₀ = 3:

>>> print(f"{bounce_frequency(1.0, 0.9, 0.1, 2.0, 3.0):.5f}")
0.02271

Functions

Name	Description
banana_width	Compute the full radial width of a banana orbit, Δr = 2√2 q κ ρ / √ε.
bounce_frequency	Compute the bounce frequency of a trapped particle, ω_b = π v √(2ε) / (4 q R₀ K(κ²)).
exb_drift	Compute the E×B drift velocity v_E = E × B / B².
grad_b_drift	Compute the grad-B drift velocity v_∇B = (m v⊥² / (2q)) B × ∇B / B³.
gyrofrequency	Compute the gyrofrequency Ω = q|B|/m.
gyroradius	Compute the gyroradius ρ = m v⊥ / (|q| B).
magnetic_moment	Compute the magnetic moment μ = m v⊥² / (2B), the adiabatic invariant of the gyromotion.
pitch_parameter	Compute the pitch parameter λ = μB₀/E from the pitch v∥/v at the poloidal angle θ.
transit_frequency	Compute the poloidal transit frequency of a passing particle.
trapped_fraction	Compute the effective fraction of trapped particles on a flux surface of a circular tokamak.
trapping_boundary	Compute the trapped–passing boundary on a flux surface.
trapping_parameter	Compute the trapping parameter κ² = (1 + ε − λ) / (2ε).

banana_widthfunction#

View source

def banana_width(gyroradius, kappa2, epsilon, safety_factor)

Compute the full radial width of a banana orbit, Δr = 2√2 q κ ρ / √ε.

From the conservation of the canonical toroidal momentum, the radial excursion is δr = q δv∥ / (ε Ω); v∥ swings between ±κ v √(2ε) at the outboard midplane. The width is largest near the trapped–passing boundary, 2√2 q ρ / √ε, the “q ρ / √ε” of scaling estimates, and vanishes for deeply trapped particles.

Parameters

Name	Type	Description
gyroradius	float or array_like	Gyroradius ρ = m v / (|q| B₀) with the total speed v (for trapped particles v⊥ ≈ v).
kappa2	float or array_like	Trapping parameter κ², 0 ≤ κ² ≤ 1.
epsilon	float or array_like	Inverse aspect ratio ε of the flux surface.
safety_factor	float or array_like	Safety factor q of the flux surface.

Returns

float or numpy.ndarray

Δr, in the units of gyroradius; nan for passing particles (κ² > 1).

References

Helander and Sigmar, Collisional Transport in Magnetized Plasmas, section 7.3.

Wesson, Tokamaks, section 3.12.

Examples

>>> print(f"{banana_width(0.01, 1.0, 0.1, 2.0):.4f}")
0.1789

bounce_frequencyfunction#

View source

def bounce_frequency(speed, kappa2, epsilon, safety_factor, major_radius)

Compute the bounce frequency of a trapped particle, ω_b = π v √(2ε) / (4 q R₀ K(κ²)).

The period of the banana orbit is τ_b = 2π/ω_b = 8 q R₀ K(κ²) / (v √(2ε)), to leading order in ε. Deeply trapped particles (κ² = 0) bounce at ω_b = v √(ε/2) / (q R₀); ω_b → 0 logarithmically at the trapped–passing boundary κ² → 1.

Parameters

Name	Type	Description
speed	float or array_like	The particle speed v.
kappa2	float or array_like	Trapping parameter κ² (see trapping_parameter()), 0 ≤ κ² < 1.
epsilon	float or array_like	Inverse aspect ratio ε of the flux surface.
safety_factor	float or array_like	Safety factor q of the flux surface.
major_radius	float or array_like	Major radius R₀.

Returns

float or numpy.ndarray

ω_b in rad per unit time; nan for passing particles (κ² > 1).

References

Helander and Sigmar, Collisional Transport in Magnetized Plasmas, section 7.2.

Wesson, Tokamaks, section 3.12 (the deeply trapped limit).

Examples

>>> # √(ε/2) deeply trapped
>>> print(f"{bounce_frequency(1.0, 0.0, 0.02, 1.0, 1.0):.4f}")
0.1000

exb_driftfunction#

View source

def exb_drift(E, B)

Compute the E×B drift velocity v_E = E × B / B².

Parameters

Name	Type	Description
E	array_like	Electric field, shape (…, 3).
B	array_like	Magnetic field, shape (…, 3).

Returns

numpy.ndarray

v_E, shape (…, 3); independent of charge and mass, |v_E| = E⊥/B.

Raises

ValueError

If a vector has no trailing axis of length 3.

References

Chen, Introduction to Plasma Physics, section 2.2.2.

Examples

>>> print(exb_drift([1.0, 0.0, 0.0], [0.0, 0.0, 2.0]))
[ 0.  -0.5  0. ]

grad_b_driftfunction#

View source

def grad_b_drift(perpendicular_speed, B, grad_B, charge=1.0, mass=1.0)

Compute the grad-B drift velocity v_∇B = (m v⊥² / (2q)) B × ∇B / B³.

Parameters

Name	Type	Default	Description
perpendicular_speed	float or array_like	required	The speed v⊥ perpendicular to B.
B	array_like	required	Magnetic field, shape (…, 3).
grad_B	array_like	required	Gradient of the field strength ∇|B|, shape (…, 3).
charge	float or array_like	1.0	Charge q; ions and electrons drift in opposite directions. Default: 1.
mass	float or array_like	1.0	Mass m. Default: 1.

Returns

numpy.ndarray

v_∇B, shape (…, 3).

Raises

ValueError

If a vector has no trailing axis of length 3.

References

Chen, Introduction to Plasma Physics, section 2.3.1.

Examples

>>> # B along z, ∇B along x
>>> print(grad_b_drift(1.0, [0.0, 0.0, 1.0], [0.1, 0.0, 0.0]))
[0.   0.05 0.  ]

gyrofrequencyfunction#

View source

def gyrofrequency(field, charge=1.0, mass=1.0)

Compute the gyrofrequency Ω = q|B|/m.

Parameters

Name	Type	Default	Description
field	float or array_like	required	Magnetic field strength |B| (its sign doesn’t matter).
charge	float or array_like	1.0	Charge q; Ω has its sign (negative for electrons). Default: 1.
mass	float or array_like	1.0	Mass m. Default: 1.

Returns

float or numpy.ndarray

Ω, in rad per unit time (a particle with q > 0 gyrates clockwise when viewed along B).

References

Chen, Introduction to Plasma Physics, section 2.2.

Examples

>>> print(gyrofrequency(2.0, charge=-1.0, mass=0.5))
-4.0

gyroradiusfunction#

View source

def gyroradius(perpendicular_speed, field, charge=1.0, mass=1.0)

Compute the gyroradius ρ = m v⊥ / (|q| B).

Parameters

Name	Type	Default	Description
perpendicular_speed	float or array_like	required	The speed v⊥ perpendicular to B.
field	float or array_like	required	Magnetic field strength |B|.
charge	float or array_like	1.0	Charge q (its sign doesn’t matter). Default: 1.
mass	float or array_like	1.0	Mass m. Default: 1.

Returns

float or numpy.ndarray

ρ (positive).

References

Chen, Introduction to Plasma Physics, section 2.2.

Examples

>>> print(gyroradius(3.0, 2.0))
1.5

magnetic_momentfunction#

View source

def magnetic_moment(perpendicular_speed, field, mass=1.0)

Compute the magnetic moment μ = m v⊥² / (2B), the adiabatic invariant of the gyromotion.

Parameters

Name	Type	Default	Description
perpendicular_speed	float or array_like	required	The speed v⊥ perpendicular to B.
field	float or array_like	required	Magnetic field strength |B|.
mass	float or array_like	1.0	Mass m. Default: 1.

Returns

float or numpy.ndarray

μ.

References

Chen, Introduction to Plasma Physics, section 2.3.

Examples

>>> print(magnetic_moment(2.0, 4.0, mass=2.0))
1.0

pitch_parameterfunction#

View source

def pitch_parameter(pitch, epsilon, theta=0.0)

Compute the pitch parameter λ = μB₀/E from the pitch v∥/v at the poloidal angle θ.

λ = (1 − (v∥/v)²) B₀/B(θ) = (1 − (v∥/v)²)(1 + ε cos θ).

Parameters

Name	Type	Default	Description
pitch	float or array_like	required	v∥/v, between −1 and 1, at the angle theta.
epsilon	float or array_like	required	Inverse aspect ratio ε of the flux surface.
theta	float or array_like	0.0	Poloidal angle θ where the pitch is given (0: outboard midplane). Default: 0.

Returns

float or numpy.ndarray

λ.

References

Helander and Sigmar, Collisional Transport in Magnetized Plasmas, section 7.1.

Examples

>>> print(f"{pitch_parameter(0.3, 0.1):.4f}")
1.0010

transit_frequencyfunction#

View source

def transit_frequency(speed, kappa2, epsilon, safety_factor, major_radius)

Compute the poloidal transit frequency of a passing particle.

ω_t = π κ v √(2ε) / (2 q R₀ K(1/κ²)), to leading order in ε: 2π/ω_t is the time for one poloidal turn. Far from the boundary (κ² ≫ 1) it tends to |v∥| / (q R₀); ω_t → 0 logarithmically at κ² → 1.

Parameters

Name	Type	Description
speed	float or array_like	The particle speed v.
kappa2	float or array_like	Trapping parameter κ² (see trapping_parameter()), κ² > 1.
epsilon	float or array_like	Inverse aspect ratio ε of the flux surface.
safety_factor	float or array_like	Safety factor q of the flux surface.
major_radius	float or array_like	Major radius R₀.

Returns

float or numpy.ndarray

ω_t in rad per unit time; nan for trapped particles (κ² < 1).

References

Helander and Sigmar, Collisional Transport in Magnetized Plasmas, section 7.2.

Examples

>>> # ≈ v∥/(qR₀) = 1
>>> print(f"{transit_frequency(1.0, 50.0, 0.01, 1.0, 1.0):.4f}")
0.9950

trapped_fractionfunction#

View source

def trapped_fraction(epsilon, approximation='exact')

Compute the effective fraction of trapped particles on a flux surface of a circular tokamak.

f_t = 1 − (3/4) ⟨B²⟩ ∫₀^{1/B_max} λ dλ / ⟨√(1 − λB)⟩ with B = B₀/(1 + ε cos θ) and the flux-surface average ⟨A⟩ = ∮ A (1 + ε cos θ) dθ / 2π, the fraction that enters neoclassical transport (bootstrap current, neoclassical resistivity).

Parameters

Name	Type	Default	Description
epsilon	float or array_like	required	Inverse aspect ratio ε = r/R₀ of the flux surface, 0 ≤ ε < 1.
approximation	('exact', 'lin-liu', 'sqrt')	"exact"	"exact": the integral above, by quadrature (relative accuracy 1e-8 or better for ε ≥ 1e-8); "lin-liu": Lin-Liu and Miller’s fit 1 − (1 − ε)² / (√(1 − ε²) (1 + 1.46 √ε)), within a few per cent at any ε; "sqrt": the small-ε limit 1.46 √ε. Default: "exact".

Returns

float or numpy.ndarray

f_t, 0 at ε = 0 and 1 at ε = 1.

Raises

ValueError

If approximation is unknown.

References

Helander and Sigmar, Collisional Transport in Magnetized Plasmas, section 11.2 (f_t ≈ 1.46 √ε).

Y. R. Lin-Liu and R. L. Miller, “Upper and lower bounds of the effective trapped particle fraction in general tokamak equilibria”, Phys. Plasmas 2, 1666 (1995).

Examples

>>> fractions = [
...     trapped_fraction(0.1, a) for a in ("exact", "lin-liu", "sqrt")
... ]
>>> print(", ".join(f"{f:.4f}" for f in fractions))
0.4492, 0.4431, 0.4617

trapping_boundaryfunction#

View source

def trapping_boundary(epsilon, quantity='lambda', theta=0.0)

Compute the trapped–passing boundary on a flux surface.

Particles with λ = μB₀/E > λ_c = B₀/B_max = 1 − ε are trapped, equivalently those with a pitch |v∥/v| < √(1 − (1 − ε)/(1 + ε cos θ)) at the angle θ (√(2ε/(1 + ε)) on the outboard midplane).

Parameters

Name	Type	Default	Description
epsilon	float or array_like	required	Inverse aspect ratio ε of the flux surface.
quantity	('lambda', 'pitch')	"lambda"	"lambda": the critical λ_c; "pitch": the critical |v∥/v| at theta. Default: "lambda".
theta	float or array_like	0.0	Poloidal angle for quantity="pitch". Default: 0.

Returns

float or numpy.ndarray

λ_c or the critical pitch.

Raises

ValueError

If quantity is unknown.

References

Wesson, Tokamaks, section 3.12.

Examples

>>> print(
...     f"{trapping_boundary(0.1):.2f}, {trapping_boundary(0.1, 'pitch'):.4f}"
... )
0.90, 0.4264

trapping_parameterfunction#

View source

def trapping_parameter(lam, epsilon)

Compute the trapping parameter κ² = (1 + ε − λ) / (2ε).

κ² < 1 for trapped particles, which bounce at sin²(θ/2) = κ²; κ² = 0 for particles deeply trapped at θ = 0 (λ = 1 + ε), κ² = 1 on the trapped–passing boundary (λ = 1 − ε) and κ² > 1 for passing ones.

Parameters

Name	Type	Description
lam	float or array_like	Pitch parameter λ = μB₀/E (see pitch_parameter()).
epsilon	float or array_like	Inverse aspect ratio ε > 0 of the flux surface.

Returns

float or numpy.ndarray

κ².

References

Helander and Sigmar, Collisional Transport in Magnetized Plasmas, section 7.2.

Examples

>>> print(f"{trapping_parameter(1.0, 0.1):.2f}")
0.50