physicskit.plasma #

Returns:

float – Alfven speed in m/s.

See also

magnetosonic_speeds: The compressive (fast/slow) counterparts of this purely magnetic wave.

Examples

>>> round(float(alfven_speed(B=1.0, rho=1e-6)), 2)
892062.06

physicskit.plasma.alfven_wave_pulse_ic(x, x0, width, amplitude)[source]#

Initial condition for a transverse Alfven-wave pulse launched from rest.

A Gaussian transverse-field pulse \(B_y(x,0)=B_1 e^{-[(x-x_0)/w]^2}\) with the transverse velocity perturbation initially zero – a “plucked string” initial condition. Because the linearized ideal-MHD Alfven-wave equations (simulate_alfven_wave()) are the same non-dispersive wave equation a string obeys, a disturbance released from rest splits exactly in half and propagates as two identical, oppositely directed pulses at \(\pm v_A\) – exactly the field-line-plucking picture Alfven’s original 1942 analogy describes.

Parameters:

x (ndarray) – Spatial grid (periodic).
x0 (float) – Pulse center.
width (float) – Gaussian pulse width.
amplitude (float) – Peak transverse field perturbation \(B_1\).

Return type:

Returns:

By0, vy0 (ndarray) – Initial transverse magnetic field and velocity perturbations, same shape as x (vy0 identically zero).

See also

simulate_alfven_wave: Evolves this initial condition forward in time.

Examples

>>> import numpy as np
>>> x = np.linspace(-10, 10, 64, endpoint=False)
>>> By0, vy0 = alfven_wave_pulse_ic(x, x0=0.0, width=1.0, amplitude=0.1)
>>> bool(np.all(vy0 == 0.0))
True

physicskit.plasma.animate_alfven_wave(By0, vy0, x, dt, steps_per_frame, n_frames, B0, rho0, mu0=1.0, interval=60)[source]#

Animate a transverse Alfven-wave pulse propagating (and splitting) along the background field.

Repeatedly calls physicskit.plasma.waves.simulate_alfven_wave() for steps_per_frame steps at a time and animates the resulting transverse-field snapshots.

Parameters:

By0 (ndarray) – Initial transverse field and velocity perturbations, e.g. from physicskit.plasma.waves.alfven_wave_pulse_ic().
vy0 (ndarray) – Initial transverse field and velocity perturbations, e.g. from physicskit.plasma.waves.alfven_wave_pulse_ic().
x (ndarray) – Uniformly spaced periodic spatial grid.
dt (float) – Time step per sub-step.
steps_per_frame (int) – Number of RK4 steps advanced between animation frames.
n_frames (int) – Number of animation frames.
B0 (float) – Background field magnitude.
rho0 (float) – Background mass density.
mu0 (float) – Vacuum permeability (normalized units by default; see physicskit.plasma.waves.simulate_alfven_wave()).
interval (int) – Delay between frames in milliseconds.

Return type:

FuncAnimation

Returns:

matplotlib.animation.FuncAnimation

Examples

>>> import numpy as np
>>> from physicskit.plasma.waves import alfven_wave_pulse_ic
>>> x = np.linspace(-20, 20, 256, endpoint=False)
>>> By0, vy0 = alfven_wave_pulse_ic(x, x0=0.0, width=1.0, amplitude=0.1)
>>> anim = animate_alfven_wave(By0, vy0, x, dt=0.002, steps_per_frame=200, n_frames=4, B0=1.0, rho0=1.0)
>>> len(list(anim.new_frame_seq()))
4

physicskit.plasma.animate_drift_wave_turbulence(phi0, dt, steps_per_frame, n_frames, length, nu=0.03, interval=60)[source]#

Animate a Hasegawa-Mima potential field developing turbulent structure from small-amplitude noise.

Repeatedly calls physicskit.plasma.turbulence.simulate_hasegawa_mima() for steps_per_frame steps at a time and animates the resulting potential-field snapshots.

Parameters:

phi0 (ndarray) – Initial potential field, e.g. from physicskit.plasma.turbulence.drift_wave_noise_ic().
dt (float) – Time step per sub-step.
steps_per_frame (int) – Number of RK4 steps advanced between animation frames.
n_frames (int) – Number of animation frames.
length (float) – Physical domain size.
nu (float) – Dissipation coefficient (see physicskit.plasma.turbulence.simulate_hasegawa_mima()).
interval (int) – Delay between frames in milliseconds.

Return type:

FuncAnimation

Returns:

matplotlib.animation.FuncAnimation

Examples

>>> import numpy as np
>>> from physicskit.plasma.turbulence import drift_wave_noise_ic
>>> phi0 = drift_wave_noise_ic(48, 2 * np.pi, amplitude=0.05, seed=0)
>>> anim = animate_drift_wave_turbulence(phi0, dt=0.02, steps_per_frame=10, n_frames=4, length=2 * np.pi)
>>> len(list(anim.new_frame_seq()))
4

physicskit.plasma.animate_ion_acoustic_soliton(u0, x, dt, steps_per_frame, n_frames, interval=60)[source]#

Animate an ion-acoustic soliton propagating without change of shape.

Repeatedly calls physicskit.plasma.waves.ion_acoustic_soliton_evolve() for steps_per_frame steps at a time and animates the resulting density-pulse snapshots.

Parameters:

u0 (ndarray) – Initial soliton profile, e.g. from physicskit.plasma.waves.ion_acoustic_soliton_profile().
x (ndarray) – Uniformly spaced periodic spatial grid.
dt (float) – Time step per sub-step.
steps_per_frame (int) – Number of steps advanced between animation frames.
n_frames (int) – Number of animation frames.
interval (int) – Delay between frames in milliseconds.

Return type:

FuncAnimation

Returns:

matplotlib.animation.FuncAnimation

Examples

>>> import numpy as np
>>> from physicskit.plasma.waves import ion_acoustic_soliton_profile
>>> x = np.linspace(-30, 30, 256, endpoint=False)
>>> u0 = ion_acoustic_soliton_profile(x, speed=4.0, x0=-15.0)
>>> anim = animate_ion_acoustic_soliton(u0, x, dt=0.0005, steps_per_frame=200, n_frames=4)
>>> len(list(anim.new_frame_seq()))
4

physicskit.plasma.animate_langmuir_wave(x0, v0, L, ng, dt, steps_per_frame, n_frames, interval=60)[source]#

Animate a Langmuir wave’s electron density oscillating in place at (approximately) the plasma frequency.

Repeatedly calls physicskit.plasma.kinetic.pic_simulate() for steps_per_frame steps at a time and, each frame, deposits the current particle positions onto the grid with physicskit.plasma.kinetic.deposit_number_density() – reusing the exact charge-assignment kernel the PIC field solve itself uses, so the density shown is precisely what physicskit.plasma.kinetic.pic_step() sees when it solves Poisson’s equation for the field driving the next sub-step.

Parameters:

x0 (ndarray) – Initial particle positions and velocities, e.g. from physicskit.plasma.kinetic.langmuir_wave_ic().
v0 (ndarray) – Initial particle positions and velocities, e.g. from physicskit.plasma.kinetic.langmuir_wave_ic().
L (float) – Domain length.
ng (int) – Number of grid points.
dt (float) – Time step per leapfrog sub-step.
steps_per_frame (int) – Number of PIC steps advanced between animation frames.
n_frames (int) – Number of animation frames.
interval (int) – Delay between frames in milliseconds.

Return type:

FuncAnimation

Returns:

matplotlib.animation.FuncAnimation

Examples

>>> import numpy as np
>>> from physicskit.plasma.kinetic import langmuir_wave_ic
>>> k = 2 * np.pi / 4.0
>>> x0, v0 = langmuir_wave_ic(4000, L=4.0, k_mode=k, alpha=0.05, v_th=0.05, seed=0)
>>> anim = animate_langmuir_wave(x0, v0, L=4.0, ng=32, dt=0.05, steps_per_frame=4, n_frames=4)
>>> len(list(anim.new_frame_seq()))
4

physicskit.plasma.animate_reconnection(psi0, eta, v0, dt, steps_per_frame, n_frames, Lx, Ly, interval=60)[source]#

Animate the flux function \(\psi(x,y,t)\) reconnecting at the X-point of a resistive current sheet.

Repeatedly calls physicskit.plasma.instabilities.simulate_reconnection() for steps_per_frame steps at a time, using each call’s final state as the next call’s initial condition, and shows the accumulated flux snapshots as an imshow animation – the antiparallel field lines above and below the sheet visibly merge into a single reconnected topology at the X-point, and the squeezed-out reconnected flux forms the outflow “jets” along the sheet.

Each snapshot has its \(x\)-mean subtracted, \(\psi(x,y,t) - \langle\psi\rangle_x(y,t)\), before display. The unperturbed Harris profile \(-L\ln\cosh(y/L)\) grows without bound away from the sheet, so on a fixed color scale it dwarfs the localized island/X-point structure that is actually reconnecting – left in, the animation reads as visually static even though the underlying field is evolving. Removing the (x-independent) background isolates exactly the x-varying perturbation that breaks and reconnects, which is what makes the merging visible frame to frame.

Parameters:

psi0 (ndarray) – Initial flux function, e.g. from physicskit.plasma.instabilities.reconnection_harris_ic().
eta (float) – Resistivity.
v0 (float) – Inflow speed.
dt (float) – Time step per simulate_reconnection() sub-step.
steps_per_frame (int) – Number of sub-steps advanced between animation frames.
n_frames (int) – Number of animation frames.
Lx (float) – Domain size.
Ly (float) – Domain size.
interval (int) – Delay between frames in milliseconds.

Return type:

FuncAnimation

Returns:

matplotlib.animation.FuncAnimation

Examples

>>> from physicskit.plasma.instabilities import reconnection_harris_ic
>>> psi0 = reconnection_harris_ic(32, 32, Lx=20.0, Ly=20.0, sheet_width=1.0, perturbation_amplitude=0.2)
>>> anim = animate_reconnection(psi0, eta=0.02, v0=0.05, dt=0.02, steps_per_frame=5, n_frames=4, Lx=20.0, Ly=20.0)
>>> len(list(anim.new_frame_seq()))
4

physicskit.plasma.animate_two_stream_phase_space(x0, v0, L, ng, dt, steps_per_frame, n_frames, interval=60)[source]#

Animate two-stream-instability phase space \((x, v)\) developing its characteristic vortex.

Repeatedly calls physicskit.plasma.kinetic.pic_simulate() for steps_per_frame leapfrog steps at a time, redrawing a scatter of every particle’s position and velocity each frame – the two initially separate beams of physicskit.plasma.kinetic.two_stream_ic() visibly wrap around each other into a single phase-space “hole” as the instability saturates.

Parameters:

x0 (ndarray) – Initial particle positions and velocities, e.g. from physicskit.plasma.kinetic.two_stream_ic().
v0 (ndarray) – Initial particle positions and velocities, e.g. from physicskit.plasma.kinetic.two_stream_ic().
L (float) – Domain length.
ng (int) – Number of grid points for the PIC field solve.
dt (float) – Time step per leapfrog sub-step.
steps_per_frame (int) – Number of PIC steps advanced between animation frames.
n_frames (int) – Number of animation frames.
interval (int) – Delay between frames in milliseconds.

Return type:

FuncAnimation

Returns:

matplotlib.animation.FuncAnimation

Examples

>>> from physicskit.plasma.kinetic import two_stream_ic
>>> x0, v0 = two_stream_ic(2000, L=10.0, v_drift=3.0, v_th=0.5, seed=0)
>>> anim = animate_two_stream_phase_space(x0, v0, L=10.0, ng=32, dt=0.05, steps_per_frame=5, n_frames=4)
>>> len(list(anim.new_frame_seq()))
4

physicskit.plasma.animate_wakefield_acceleration(x0, v0, q, m, E0, k, v_phase, dt, steps, frame_stride=20, interval=60)[source]#

Animate a test charge surfing a prescribed traveling wakefield, with an energy-gain trace inset.

Runs physicskit.plasma.acceleration.simulate_wakefield_acceleration() once for the full trajectory, then animates a snapshot of the wakefield \(E_z(x, t)\) with the particle’s position marked on it, alongside a running plot of its kinetic energy – showing both the spatial “surfing” picture and the resulting energy gain simultaneously.

Parameters:

x0 (float) – Initial particle position and velocity.
v0 (float) – Initial particle position and velocity.
q (float) – Particle charge and mass.
m (float) – Particle charge and mass.
E0 (float) – Wakefield amplitude, wavenumber, and phase velocity, as in physicskit.plasma.acceleration.wakefield_e_field().
k (float) – Wakefield amplitude, wavenumber, and phase velocity, as in physicskit.plasma.acceleration.wakefield_e_field().
v_phase (float) – Wakefield amplitude, wavenumber, and phase velocity, as in physicskit.plasma.acceleration.wakefield_e_field().
dt (float) – Time step.
steps (int) – Total number of Boris-pusher steps to integrate.
frame_stride (int) – Number of integration steps between animation frames.
interval (int) – Delay between frames in milliseconds.

Return type:

FuncAnimation

Returns:

matplotlib.animation.FuncAnimation

Examples

>>> anim = animate_wakefield_acceleration(
...     x0=0.0, v0=0.9, q=1.0, m=1.0, E0=0.05, k=1.0, v_phase=1.0, dt=0.01, steps=400, frame_stride=40
... )
>>> len(list(anim.new_frame_seq()))
11

physicskit.plasma.animate_weibel_filamentation(x, t, wpe, temperature_anisotropy, n_modes=12, seed=0, interval=60)[source]#

Animate the transverse current filaments growing under the reduced quasi-linear Weibel model.

Parameters:

x (ndarray) – Spatial grid.
t (ndarray) – Animation frame times.
wpe (float) – Electron plasma frequency.
temperature_anisotropy (float) – The ratio \(T_\perp/T_\parallel\); see physicskit.plasma.instabilities.weibel_growth_rate().
n_modes (int) – Number of seeded Fourier modes.
seed (int) – Random seed.
interval (int) – Delay between frames in milliseconds.

Return type:

FuncAnimation

Returns:

matplotlib.animation.FuncAnimation

See also

physicskit.plasma.instabilities.simulate_weibel_filamentation: Supplies the current-density snapshots animated here.

Examples

>>> import numpy as np
>>> x = np.linspace(0, 20.0, 128, endpoint=False)
>>> t = np.linspace(0, 5.0, 6)
>>> anim = animate_weibel_filamentation(x, t, wpe=1.0, temperature_anisotropy=4.0, n_modes=6)
>>> len(list(anim.new_frame_seq()))
6

physicskit.plasma.boris_integrate(pos0, vel0, q, m, E, B, dt, steps)[source]#

Integrate a charged particle’s trajectory in uniform E and B fields with the Boris pusher.

Repeatedly applies boris_push() inside a single Numba-compiled loop, avoiding Python-level overhead per step – the throughput this buys is what makes the Boris scheme practical for particle-in-cell codes tracking millions of particles.

Parameters:

pos0 (ndarray) – Initial position in meters.
vel0 (ndarray) – Initial velocity in m/s.
q (float) – Particle charge in Coulombs.
m (float) – Particle mass in kilograms.
E (ndarray) – Uniform electric field in V/m.
B (ndarray) – Uniform magnetic field in Tesla.
dt (float) – Time step in seconds.
steps (int) – Number of steps to advance.

Return type:

Returns:

pos_hist (ndarray, shape (steps + 1, 3)) – Position at every step, including the initial condition.
vel_hist (ndarray, shape (steps + 1, 3)) – Velocity at every step, including the initial condition.

See also

boris_push: The single-step update repeated here.
exb_drift: The closed-form drift velocity this trajectory averages to when a non-zero E is present.

Examples

Kinetic energy is conserved to machine precision over many gyrations in a pure magnetic field:

>>> import numpy as np
>>> B = np.array([0.0, 0.0, 1.0])
>>> E = np.zeros(3)
>>> omega_c = cyclotron_frequency(QE, MP, 1.0)
>>> dt = (2 * np.pi / omega_c) / 200
>>> pos_hist, vel_hist = boris_integrate(
...     np.zeros(3), np.array([1e5, 0.0, 0.0]), QE, MP, E, B, dt, steps=2000
... )
>>> speeds = np.linalg.norm(vel_hist, axis=1)
>>> bool(np.max(np.abs(speeds - 1e5)) < 1e-6)
True

physicskit.plasma.boris_push(pos, vel, q, m, E, B, dt)[source]#

Advance a charged particle one leapfrog step with the Boris integrator.

The Boris (1970) scheme splits each step into an electric half-acceleration, an exact rotation about \(\mathbf{B}\) (via the Boris “E-cross-B” trick, avoiding any explicit trigonometry), and a second electric half-acceleration. The rotation step preserves speed exactly in a pure magnetic field, which is why the scheme conserves kinetic energy over arbitrarily many gyro-orbits where a naive forward-Euler or even RK4 update would spiral outward.

Parameters:

pos (ndarray) – Particle position \([x, y, z]\) in meters.
vel (ndarray) – Particle velocity \([v_x, v_y, v_z]\) in m/s.
q (float) – Particle charge in Coulombs.
m (float) – Particle mass in kilograms.
E (ndarray) – Electric field \([E_x, E_y, E_z]\) in V/m.
B (ndarray) – Magnetic field \([B_x, B_y, B_z]\) in Tesla.
dt (float) – Time step in seconds.

Return type:

Returns:

pos_new, vel_new (ndarray, shape (3,)) – Position and velocity after one step of size dt.

See also

boris_integrate: Repeated application of this step over many time steps.

Examples

A proton launched perpendicular to a uniform field gyrates without gaining or losing speed:

>>> import numpy as np
>>> r0 = np.array([0.0, 0.0, 0.0])
>>> v0 = np.array([1e5, 0.0, 0.0])
>>> E = np.zeros(3)
>>> B = np.array([0.0, 0.0, 1.0])
>>> r1, v1 = boris_push(r0, v0, q=QE, m=MP, E=E, B=B, dt=1e-9)
>>> round(float(np.linalg.norm(v1)), 6)
100000.0

physicskit.plasma.cma_coordinates(omega, wpe, wce)[source]#

Dimensionless \((X, Y)\) coordinates of the Clemmow-Mullaly-Allis (CMA) diagram.

The CMA diagram partitions the plane \(X=\omega_{pe}^2/\omega^2\) (density axis) vs. \(Y=\omega_{ce}/\omega\) (field-strength axis) into regions of distinct wave topology – cutoffs, resonances, and the number and polarization of propagating modes – giving a single map of every cold-plasma wave regime from ordinary light waves (\(X, Y \to 0\)) to the Alfven wave and whistler branches (\(Y \gg 1\)).

Parameters:

omega (float) – Wave angular frequency in rad/s.
wpe (float) – Electron plasma frequency in rad/s.
wce (float) – Electron cyclotron frequency in rad/s.

Return type:

Returns:

X, Y (float) – CMA diagram coordinates.

Examples

>>> cma_coordinates(omega=1.0, wpe=2.0, wce=3.0)
(4.0, 3.0)

physicskit.plasma.cold_plasma_dispersion(theta, S, D, P)[source]#

Solve the cold-plasma dispersion relation for the squared refractive index \(n^2\).

Substituting a plane wave into \(\mathbf{n}\times(\mathbf{n}\times\mathbf{E})+\mathbf{K}\cdot\mathbf{E}=0\) for propagation at angle \(\theta\) to \(\mathbf{B}_0\) gives the Appleton-Hartree biquadratic \(An^4 - Bn^2 + C = 0\) with

\[A = S\sin^2\theta + P\cos^2\theta, \quad B = RL\sin^2\theta + PS(1+\cos^2\theta), \quad C = PRL,\]

which reduces at \(\theta=0\) to the decoupled R- and L-waves (\(n^2=R\) or \(L\)) and at \(\theta=\pi/2\) to the O-mode (\(n^2=P\)) and X-mode (\(n^2=RL/S\)).

Parameters:

theta (float) – Angle between the wavevector and \(\mathbf{B}_0\), in radians.
S (float) – Stix parameters from stix_parameters().
D (float) – Stix parameters from stix_parameters().
P (float) – Stix parameters from stix_parameters().

Return type:

Returns:

n_sq_plus, n_sq_minus (float) – The two roots of the biquadratic (the two cold-plasma wave branches at this angle and frequency). A negative root means that branch is evanescent rather than propagating.

See also

stix_parameters: Supplies S, D, P.

Examples

Parallel propagation recovers the pure R- and L-wave refractive indices:

>>> import numpy as np
>>> S, D, P = -2.314262935091991, 90.69535621969021, -7959.851640342367
>>> n2_plus, n2_minus = cold_plasma_dispersion(theta=0.0, S=S, D=D, P=P)
>>> R, L = rl_parameters(S, D)
>>> bool(np.isclose(sorted([n2_plus, n2_minus]), sorted([R, L])).all())
True

physicskit.plasma.curvature_drift(v_par, q, m, B, R_c)[source]#

The curvature drift \(\mathbf{v}_R = \dfrac{m v_\parallel^2}{q}\dfrac{\mathbf{R}_c\times\mathbf{B}}{R_c^2 B^2}\).

A particle streaming along a curved field line feels a centrifugal force in the guiding-center frame, which – crossed with \(\mathbf{B}\) – produces a drift perpendicular to the plane of curvature. In a low-beta toroidal equilibrium this combines with grad_b_drift() (curvature and field-strength gradients point the same way when \(\nabla\times\mathbf{B}=0\) in vacuum) into the single \(\nabla B\) + curvature drift responsible for charge separation and the resulting E-cross-B rotation in a tokamak.

Parameters:

v_par (float) – Speed parallel to the field, in m/s.
q (float) – Particle charge in Coulombs (signed).
m (float) – Particle mass in kilograms.
B (ndarray) – Local magnetic field vector in Tesla.
R_c (ndarray) – Radius-of-curvature vector, pointing from the field line’s local center of curvature to the particle, magnitude \(R_c\) in meters.

Return type:

Returns:

ndarray, shape (3,) – Drift velocity in m/s.

Examples

>>> import numpy as np
>>> B = np.array([0.0, 0.0, 1.0])
>>> R_c = np.array([1.0, 0.0, 0.0])
>>> drift = curvature_drift(1e5, QE, MP, B, R_c)
>>> round(float(drift[1]), 2)
-104.4

physicskit.plasma.cyclotron_frequency(q, m, B)[source]#

Angular cyclotron (gyration) frequency \(\omega_c = qB/m\).

Parameters:

q (float) – Particle charge in Coulombs (signed).
m (float) – Particle mass in kilograms.
B (float) – Magnetic field magnitude in Tesla.

Return type:

Returns:

float – Signed angular gyrofrequency in rad/s (negative for negatively charged particles, reflecting the opposite sense of rotation).

See also

larmor_radius: The orbit radius set by this frequency.

Examples

A proton in a 1 T field gyrates about 15.3 MHz:

>>> round(cyclotron_frequency(QE, MP, 1.0) / (2 * 3.141592653589793) / 1e6, 2)
15.25

physicskit.plasma.deposit_number_density(x, L, ng, n0=1.0)[source]#

Deposit particle positions onto a grid as a number density, via cloud-in-cell (CIC) weighting.

Each particle represents a “cloud” of physical charge spanning one grid cell, split linearly between its two nearest grid points – the standard first-order PIC weighting scheme, chosen because it is exact for a uniform density and (unlike nearest-grid-point deposit) produces a smooth, differentiable force with no self-force discontinuities as particles cross cell boundaries.

Parameters:

x (ndarray) – Particle positions in \([0, L)\).
L (float) – Domain length (periodic).
ng (int) – Number of grid points.
n0 (float) – Equilibrium number density (sets each particle’s statistical weight, \(n_0 L / n_{particles}\)).

Return type:

Returns:

ndarray, shape (ng,) – Number density on the grid.

See also

interpolate_field: The companion gather operation.

Examples

Total deposited charge exactly equals the physical charge represented, regardless of how the particles are distributed:

>>> import numpy as np
>>> x = np.array([0.1, 2.4, 4.9, 7.7])
>>> rho = deposit_number_density(x, L=10.0, ng=20, n0=2.0)
>>> dx = 10.0 / 20
>>> round(float(np.sum(rho) * dx), 8)
20.0

physicskit.plasma.drift_wave_noise_ic(n, length, amplitude, seed=0)[source]#

Small-amplitude random-phase noise initial condition for Hasegawa-Mima drift-wave turbulence.

A featureless, isotropic random field with the specified root-mean-square amplitude – the standard “let it find its own structure” initial condition for turbulence simulations, in place of a specific unstable linear eigenmode, since real drift-wave turbulence in a fusion device is continuously driven by many linearly unstable modes simultaneously rather than growing from one clean seed.

Parameters:

n (int) – Number of grid points along each axis.
length (float) – Physical domain size. Unused, since the noise is independent per grid point (white) and so has no length scale; kept for signature symmetry with simulate_hasegawa_mima().
amplitude (float) – Root-mean-square amplitude of the seeded potential noise.
seed (int) – Random seed.

Return type:

Returns:

ndarray, shape (n, n) – Initial potential \(\phi(x, y)\).

See also

simulate_hasegawa_mima: Evolves this initial condition forward in time.

Examples

>>> phi0 = drift_wave_noise_ic(64, 2 * 3.141592653589793, amplitude=0.01, seed=0)
>>> phi0.shape
(64, 64)

physicskit.plasma.exb_drift(E, B)[source]#

The \(\mathbf{E}\times\mathbf{B}\) drift velocity \(\mathbf{v}_E = (\mathbf{E}\times\mathbf{B})/B^2\).

Unlike every other guiding-center drift, this one is independent of charge, mass, and energy: electrons and ions drift together at the same velocity, so \(\mathbf{E}\times\mathbf{B}\) drift carries no net current and instead advects the whole plasma as a fluid – the drift underlying tokamak radial-electric-field rotation and magnetospheric convection alike.

Parameters:

E (ndarray) – Electric field in V/m.
B (ndarray) – Magnetic field in Tesla.

Return type:

Returns:

ndarray, shape (3,) – Drift velocity in m/s.

Examples

>>> import numpy as np
>>> E = np.array([0.0, 1e3, 0.0])
>>> B = np.array([0.0, 0.0, 1.0])
>>> exb_drift(E, B)
array([1000.,    0.,    0.])

physicskit.plasma.grad_b_drift(v_perp, q, m, B, grad_B)[source]#

The grad-B drift \(\mathbf{v}_{\nabla B} = \dfrac{m v_\perp^2}{2qB^3}(\mathbf{B}\times\nabla B)\).

A particle gyrating in a field whose magnitude varies across the orbit sees a tighter turn (smaller Larmor radius) on the strong-field side than the weak-field side, so the orbit fails to close and the guiding center creeps sideways. Because the drift is inversely proportional to charge, electrons and ions drift in opposite directions, producing a net current – the origin of the ring current in planetary magnetospheres.

Parameters:

v_perp (float) – Speed perpendicular to the field, in m/s.
q (float) – Particle charge in Coulombs (signed).
m (float) – Particle mass in kilograms.
B (ndarray) – Local magnetic field vector in Tesla.
grad_B (ndarray) – Gradient of the field magnitude, \(\nabla |\mathbf{B}|\), in Tesla/meter.

Return type:

Returns:

ndarray, shape (3,) – Drift velocity in m/s.

See also

curvature_drift: The companion drift from field-line curvature, which combines with this one in any real toroidal field.

Examples

>>> import numpy as np
>>> B = np.array([0.0, 0.0, 1.0])
>>> grad_B = np.array([0.1, 0.0, 0.0])
>>> drift = grad_b_drift(1e5, QE, MP, B, grad_B)
>>> round(float(drift[1]), 4)
5.2198

physicskit.plasma.grad_shafranov_rhs_solovev(R, c1, c2)[source]#

Source term \(\Delta^*\psi = c_1 R^2 + c_2\) of the linear Solov’ev equilibrium.

Equivalent to \(-\mu_0 R^2 p'(\psi) - FF'(\psi)\) in the Grad-Shafranov equation for the special case of constant \(p'(\psi)\) and \(FF'(\psi)\).

Parameters:

R (ndarray) – Major-radius coordinate, in meters.
c1 (float) – Source coefficients (related to \(p'\) and \(FF'\)).
c2 (float) – Source coefficients (related to \(p'\) and \(FF'\)).

Return type:

Returns:

ndarray – The right-hand side \(\Delta^*\psi\), same shape as R.

Examples

>>> import numpy as np
>>> grad_shafranov_rhs_solovev(np.array([1.0, 2.0]), c1=1.0, c2=-2.0)
array([-1.,  2.])

physicskit.plasma.hasegawa_mima_rhs(q, KX, KY, K2)[source]#

Non-dissipative right-hand side of the Hasegawa-Mima potential-vorticity equation, evaluated pseudo-spectrally.

Evolves the potential vorticity \(q=\nabla^2\phi-\phi\) (from which \(\phi\) is recovered via \(\hat{\phi}=-\hat{q}/(1+k^2)\), the elliptic inversion analogous to the streamfunction Poisson solve of 2D Navier-Stokes) under \(\mathbf{E}\times\mathbf{B}\) advection by the drift velocity \((u,v)=(-\partial_y\phi,\partial_x\phi)\) and the linear drift-wave term \(-\partial_y\phi\). Excludes the dissipative regularization simulate_hasegawa_mima() adds on top, which – being stiffer at grid scale than this advective term – is integrated separately via an exact integrating factor rather than folded into this explicit right-hand side.

Parameters:

q (ndarray) – Potential vorticity \(q = \nabla^2\phi - \phi\).
KX (ndarray) – Wavenumber grids from _spectral_grid().
KY (ndarray) – Wavenumber grids from _spectral_grid().
K2 (ndarray) – Wavenumber grids from _spectral_grid().

Return type:

Returns:

ndarray, shape (n, n) – \(dq/dt\), excluding dissipation.

physicskit.plasma.interpolate_field(x, field_grid, L)[source]#

Interpolate a grid-defined field to particle positions, via cloud-in-cell (CIC) weighting.

The gather step dual to deposit_number_density(): using the same linear weights for both deposit and gather is what makes the PIC method momentum-conserving (no self-force on an isolated particle).

Parameters:

x (ndarray) – Particle positions in \([0, L)\).
field_grid (ndarray) – Field values on the grid, e.g. from solve_poisson_1d().
L (float) – Domain length (periodic).

Return type:

Returns:

ndarray, shape (n_particles,) – Field value at each particle’s position.

Examples

A particle sitting exactly on a grid node picks up that node’s value:

>>> import numpy as np
>>> ng = 8
>>> field_grid = np.array([1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0])
>>> x = np.array([2 * (10.0 / ng)])
>>> interpolate_field(x, field_grid, L=10.0)
array([3.])

physicskit.plasma.ion_acoustic_soliton_evolve(u0, x, dt, steps)[source]#

Evolve an ion-acoustic KdV initial condition forward in time on a periodic domain.

A pseudo-spectral Strang-split scheme: the stiff linear dispersion \(\partial_\xi^3\) is advanced exactly in Fourier space, and the non-stiff nonlinear advection \(6uu_\xi\) with RK4 in between – structurally the standard splitting for any KdV-type equation, applied here directly to the ion-acoustic reduction rather than by importing a generic KdV solver, since the physical field this equation governs (a Debye-length-normalized density/potential perturbation moving at order the ion-sound speed) is specific to this module.

Parameters:

u0 (ndarray) – Initial field, sampled on the periodic grid x.
x (ndarray) – Uniformly spaced periodic spatial grid.
dt (float) – Time step.
steps (int) – Number of steps to advance.

Return type:

Returns:

ndarray – Field after steps * dt time units.

See also

ion_acoustic_soliton_profile: Exact traveling-wave solution this reproduces.

Examples

A soliton launched at speed \(c=4\) has advanced by very close to \(c\cdot(\text{steps}\cdot dt)\) and kept its amplitude, since KdV solitons propagate without changing shape:

>>> import numpy as np
>>> N, L = 512, 60.0
>>> x = np.linspace(-L / 2, L / 2, N, endpoint=False)
>>> u0 = ion_acoustic_soliton_profile(x, speed=4.0, x0=-15.0)
>>> u = ion_acoustic_soliton_evolve(u0, x, dt=0.0005, steps=4000)
>>> shift = x[np.argmax(u)] - x[np.argmax(u0)]
>>> bool(abs(shift - 4.0 * 4000 * 0.0005) < 0.5)
True
>>> bool(abs(u.max() - u0.max()) < 0.05)
True

physicskit.plasma.ion_acoustic_soliton_profile(x, speed, x0=0.0)[source]#

Exact single-soliton solution of the ion-acoustic Korteweg-de Vries reduction.

The reductive-perturbation (Washimi & Taniuti, 1966) expansion of the cold-ion-fluid/Boltzmann-electron equations in the weakly nonlinear, weakly dispersive limit reduces the ion-acoustic wave problem to the KdV equation for the normalized density (or potential) perturbation \(u\) in a frame moving at the ion-sound speed; after the standard rescaling of the stretched coordinates it takes the canonical form \(u_t + 6uu_\xi + u_{\xi\xi\xi} = 0\) – the same normal form every weakly-dispersive weakly-nonlinear wave problem reduces to, with the ion-acoustic-specific physics (electron Boltzmann response supplying the nonlinearity, ion inertia and Debye-length dispersion supplying the \(\partial_\xi^3\) term) fixing only the physical unit conversions between \(u,\xi,t\) and density, position, and time in the ion-sound-speed frame. Its exact traveling-wave solution is a single soliton of speed \(c\) (in the stretched frame) and amplitude \(c/2\), propagating without change of shape – consistent with a real ion-acoustic soliton, whose speed always exceeds the linear sound speed by an amount set by its amplitude.

Parameters:

x (ndarray) – Spatial grid, in the stretched (ion-sound-speed) frame.
speed (float) – Soliton speed \(c\) in the stretched frame (amplitude \(=c/2\)); must be positive.
x0 (float) – Initial center position.

Return type:

Returns:

ndarray – \(u(x, 0) = \tfrac{c}{2}\,\mathrm{sech}^2\!\big(\tfrac{\sqrt{c}}{2}(x-x_0)\big)\).

See also

ion_acoustic_soliton_evolve: Propagate this (or any) initial condition forward in time.

Examples

>>> round(float(ion_acoustic_soliton_profile(0.0, speed=4.0)), 6)
2.0

physicskit.plasma.landau_damping_ic(n_particles, L, k_mode, alpha, v_th, seed=0)[source]#

Quiet-start initial condition seeding a single-mode density perturbation for a Landau damping test.

Displaces an otherwise uniform particle load by \(\delta x = (\alpha/k)\sin(kx_0)\), which – since particle number is conserved through the displacement’s Jacobian – produces exactly the density perturbation \(n(x) \approx n_0(1 - \alpha\cos(kx))\) used in the classic linear Landau damping test problem, without the sampling noise a random density draw would add on top of the intended signal.

Parameters:

n_particles (int) – Number of particles.
L (float) – Domain length (periodic), in normalized length units. Choosing \(L = 2\pi/k_{mode}\) fits exactly one wavelength.
k_mode (float) – Wavenumber of the seeded perturbation.
alpha (float) – Perturbation amplitude (\(\alpha \ll 1\) for the linear regime).
v_th (float) – Thermal speed of the background Maxwellian.
seed (int) – Random seed for the velocity sampling.

Return type:

Returns:

x, v (ndarray, shape (n_particles,)) – Particle positions (in \([0, L)\)) and velocities.

See also

pic_simulate: Evolve this initial condition forward in time.
landau_damping_rate: The analytic decay rate this test is checked against.

Examples

>>> import numpy as np
>>> x, v = landau_damping_ic(4, L=4.0, k_mode=2 * 3.141592653589793 / 4.0, alpha=0.1, v_th=1.0, seed=0)
>>> bool(np.all((x >= 0) & (x < 4.0)))
True
>>> x.shape, v.shape
((4,), (4,))

physicskit.plasma.landau_damping_rate(k, v_th, omega_pe=1.0)[source]#

Analytic linear Landau damping rate for a Maxwellian electron plasma.

The classic weak-damping result (Landau, 1946) for a wave of wavenumber \(k\) on a Maxwellian of thermal speed \(v_{th}\), with Debye length \(\lambda_D = v_{th}/\omega_{pe}\):

\[\gamma = -\omega_{pe}\sqrt{\frac{\pi}{8}}\, \frac{1}{(k\lambda_D)^3}\, \exp\!\left(-\frac{1}{2(k\lambda_D)^2} - \frac{3}{2}\right).\]

Electrons resonant with the wave’s phase velocity (\(v = \omega/k\)) surf it, extracting energy from the field on net because the Maxwellian has slightly more slower particles being accelerated than faster particles being decelerated – a purely collisionless damping mechanism with no dissipation at the particle level, reproduced here by pic_simulate() without any explicit damping term in the equations of motion.

Parameters:

k (float) – Wavenumber, in units of inverse Debye length times \(k\lambda_D\) conventions – concretely, pass the physical wavenumber and set v_th/omega_pe consistently.
v_th (float) – Electron thermal speed.
omega_pe (float) – Electron plasma frequency (1.0 in the normalized units used throughout this module).

Return type:

Returns:

float – Damping rate \(\gamma\) (negative, since the wave decays). The formula is only accurate for weak damping, \(k\lambda_D \lesssim 0.5\).

See also

pic_simulate: Numerically reproduces this decay from first principles.

Examples

The standard textbook benchmark case, \(k\lambda_D = 0.5\):

>>> round(float(landau_damping_rate(k=0.5, v_th=1.0)), 4)
-0.1514

physicskit.plasma.langmuir_wave_ic(n_particles, L, k_mode, alpha, v_th, seed=0)[source]#

Small-amplitude single-mode initial condition for a Langmuir (electron plasma) wave.

Uses the same quiet-start displacement as landau_damping_ic(), \(\delta x=(\alpha/k)\sin(kx_0)\) (density perturbation \(n_0(1-\alpha\cos(kx))\)), and additionally gives each particle the coherent velocity \(\delta v = (\alpha\,\omega_{pe}/k)\sin(kx_0)\) (\(\omega_{pe}=1\) in these normalized units). In cold-fluid linear theory the displacement then evolves as \(\xi(x_0,t) = (\alpha/k)\sin(kx_0)\,[\cos\omega_{pe}t + \sin\omega_{pe}t]\): a standing wave ringing in place at \(\omega_{pe}\), with \(\sqrt2\) the amplitude (twice the field energy) of the displacement-only start and a \(\pi/4\) phase shift.

How long it rings is set by Landau damping, i.e. by \(k\lambda_D\) (see landau_damping_rate()), not by the velocity kick: for \(k_{mode}\,v_{th} \ll \omega_{pe}\) this and landau_damping_ic() both oscillate essentially undamped, while for \(k\lambda_D \gtrsim 0.3\) both damp.

Parameters:

n_particles (int) – Number of particles.
L (float) – Domain length (periodic); \(L=2\pi/k_{mode}\) fits one wavelength.
k_mode (float) – Wavenumber of the seeded standing wave.
alpha (float) – Perturbation amplitude (\(\alpha \ll 1\)).
v_th (float) – Thermal speed of the background Maxwellian; keep \(k_{mode}\,v_{th} \ll \omega_{pe}=1\) for the wave to be only weakly Landau-damped.
seed (int) – Random seed for the thermal velocity sampling.

Return type:

Returns:

x, v (ndarray, shape (n_particles,)) – Particle positions and velocities.

See also

landau_damping_ic: The companion density-only perturbation (no coherent velocity kick); it damps or rings for the same \(k\lambda_D\) as this one.
pic_simulate: Evolve this initial condition forward in time.

Examples

>>> import numpy as np
>>> k = 2 * np.pi / 4.0
>>> x, v = langmuir_wave_ic(4, L=4.0, k_mode=k, alpha=0.05, v_th=0.05, seed=0)
>>> bool(np.all((x >= 0) & (x < 4.0)))
True
>>> x.shape, v.shape
((4,), (4,))

physicskit.plasma.larmor_radius(v_perp, q, m, B)[source]#

Larmor (gyro) radius \(r_L = m v_\perp / (|q| B)\).

Parameters:

v_perp (float) – Speed perpendicular to the magnetic field, in m/s.
q (float) – Particle charge in Coulombs (sign is ignored).
m (float) – Particle mass in kilograms.
B (float) – Magnetic field magnitude in Tesla.

Return type:

Returns:

float – Gyroradius in meters.

Examples

>>> round(larmor_radius(1e5, QE, MP, 1.0), 6)
0.001044

physicskit.plasma.lundquist_number(L, vA, eta)[source]#

Lundquist number \(S = L v_A/\eta\), the magnetic Reynolds number built from the Alfven speed.

Parameters:

L (float) – Characteristic length scale (e.g. current-sheet length), in meters.
vA (float) – Alfven speed in m/s, from alfven_speed().
eta (float) – Magnetic diffusivity in m^2/s.

Return type:

Returns:

float – Lundquist number (dimensionless). Fusion and astrophysical plasmas typically have \(S \sim 10^{6}\) – \(10^{14}\).

See also

sweet_parker_rate: Reconnection rate scaling as \(S^{-1/2}\).

Examples

>>> lundquist_number(L=1.0, vA=1e6, eta=1.0)
1000000.0

physicskit.plasma.magnetic_mirror_bounce(z0, v_par0, v_perp0, m, B_func, dz=1e-06, dt=1e-10, steps=20000)[source]#

Simulate 1D guiding-center bounce motion between the throats of a magnetic mirror.

Integrates \(m\dot{v}_\parallel = -\mu\, dB/dz\) with \(\mu\) fixed at its initial value (the adiabatic invariant), using a symmetric leapfrog step. A particle with too little pitch angle to reflect before reaching the mirror throat’s peak field instead falls into the loss cone and would be lost from confinement in a real device; mirror_force() supplies the underlying force law.

Parameters:

z0 (float) – Initial position along the field line, in meters.
v_par0 (float) – Initial parallel velocity, in m/s.
v_perp0 (float) – Initial perpendicular velocity, in m/s (sets \(\mu\) via magnetic_moment()).
m (float) – Particle mass in kilograms.
B_func (callable) – Field-strength profile B_func(z) -> float along the field line, in Tesla.
dz (float) – Finite-difference step used to evaluate \(dB/dz\), in meters.
dt (float) – Time step in seconds.
steps (int) – Number of leapfrog steps to advance.

Return type:

Returns:

z_hist (ndarray, shape (steps + 1,)) – Position along the field line at every step.
v_par_hist (ndarray, shape (steps + 1,)) – Parallel velocity at every step.

Examples

A particle launched from the mirror midplane with enough perpendicular energy reflects before reaching the throat, reversing the sign of its parallel velocity:

>>> import numpy as np
>>> B_func = lambda z: 1.0 + 4.0 * (z / 0.05) ** 2
>>> z_hist, v_par_hist = magnetic_mirror_bounce(
...     z0=0.0, v_par0=2e4, v_perp0=8e4, m=MP, B_func=B_func, steps=6000
... )
>>> bool(v_par_hist[0] > 0 and v_par_hist[-1] < 0)
True

physicskit.plasma.magnetic_moment(v_perp, m, B)[source]#

The first adiabatic invariant \(\mu = m v_\perp^2 / (2B)\).

Conserved for a charged particle whose gyration is fast compared to any change in the field it sees, \(\mu\) acts as a magnetic “potential energy per unit field”: as the particle drifts into stronger \(B\), \(v_\perp\) must grow to keep \(\mu\) fixed, converting parallel kinetic energy into perpendicular kinetic energy. That conversion is the mechanism behind mirror_force().

Parameters:

v_perp (float) – Speed perpendicular to the magnetic field, in m/s.
m (float) – Particle mass in kilograms.
B (float) – Magnetic field magnitude in Tesla.

Return type:

Returns:

float – Magnetic moment in Joules/Tesla.

See also

mirror_force: The parallel force derived from this invariant.

Examples

>>> round(magnetic_moment(1e5, MP, 1.0) * 1e18, 4)
8.3631

physicskit.plasma.magnetosonic_speeds(vA, cs, theta)[source]#

Fast and slow magnetosonic phase speeds at propagation angle \(\theta\) to \(\mathbf{B}\).

The two compressive MHD normal modes solve \(v_{f,s}^2 = \tfrac{1}{2}\left[(v_A^2+c_s^2) \pm \sqrt{(v_A^2+c_s^2)^2 - 4v_A^2c_s^2\cos^2\theta}\right]\). At \(\theta=0\) (propagation along \(\mathbf{B}\)) they reduce to \(\max(v_A, c_s)\) and \(\min(v_A, c_s)\); at \(\theta=\pi/2\) the fast mode becomes the purely compressive \(\sqrt{v_A^2+c_s^2}\) and the slow mode vanishes, since a perpendicular perturbation cannot bend field lines that are already perpendicular to its wavevector.

Parameters:

vA (float) – Alfven speed in m/s, from alfven_speed().
cs (float) – Sound speed in m/s, from sound_speed().
theta (float) – Angle between the wavevector and \(\mathbf{B}\), in radians.

Return type:

Returns:

v_fast, v_slow (float) – Fast and slow magnetosonic phase speeds in m/s.

Examples

>>> import numpy as np
>>> vf, vs = magnetosonic_speeds(vA=892062.06, cs=1e5, theta=np.pi / 2)
>>> round(vf, 2)
897649.55
>>> round(vs, 2)
0.0

physicskit.plasma.maxwellian_velocities(n_particles, v_th, seed=0)[source]#

Sample particle velocities from a Maxwellian (Gaussian) distribution.

Parameters:

n_particles (int) – Number of particles to sample.
v_th (float) – Thermal speed (standard deviation of the Gaussian), in normalized velocity units.
seed (int) – Seed for the pseudo-random number generator, for reproducibility.

Return type:

Returns:

ndarray, shape (n_particles,) – Sampled velocities.

Examples

>>> import numpy as np
>>> v = maxwellian_velocities(4, v_th=1.0, seed=0)
>>> bool(np.allclose(v, [0.12573022, -0.13210486, 0.64042265, 0.10490012]))
True

physicskit.plasma.mirror_force(mu, grad_B_parallel)[source]#

The mirror force \(F_\parallel = -\mu\, \partial B/\partial \ell\) along a field line.

As a particle’s guiding center moves into a region of stronger field (larger \(\partial B/\partial\ell\)), the conservation of magnetic_moment() forces \(v_\perp\) to grow at the expense of \(v_\parallel\); this is the reaction force decelerating the parallel motion. If the field is strong enough, \(v_\parallel\) reaches zero before the particle passes the throat and it reflects – magnetic mirror confinement, simulated end to end in magnetic_mirror_bounce().

Parameters:

mu (float) – Magnetic moment in J/T, from magnetic_moment().
grad_B_parallel (float) – Gradient of the field magnitude along the field line, in Tesla/meter.

Return type:

Returns:

float – Force in Newtons, directed to push the particle toward weaker field.

Examples

>>> mirror_force(mu=1e-17, grad_B_parallel=2.0)
-2e-17

physicskit.plasma.petschek_rate(S)[source]#

Petschek reconnection rate \(v_{in}/v_A \approx \pi/(8\ln S)\).

Petschek (1964) showed that if the diffusion region shrinks to a small X-point rather than the full Sweet-Parker sheet length, four standing slow-mode shocks can carry most of the inflowing flux and energy conversion, giving a reconnection rate that falls only logarithmically with \(S\) instead of as \(S^{-1/2}\) – fast enough to plausibly explain solar flare and magnetospheric substorm timescales.

Parameters:: S (float) – Lundquist number, from lundquist_number().
Return type:: float
Returns:: float – Dimensionless reconnection rate \(v_{in}/v_A\).

See also

sweet_parker_rate: The slower, steady-sheet reconnection rate this improves on.

Examples

>>> round(float(petschek_rate(S=1e6)), 4)
0.0284

physicskit.plasma.pic_simulate(x0, v0, L, ng, dt, steps, qm=-1.0, n0=1.0)[source]#

Run an electrostatic PIC simulation forward in time, recording the field-energy history.

Correctly initializes the leapfrog velocity offset (staggering v0 back by half a step using the field at x0) before repeatedly applying pic_step().

Parameters:

x0 (ndarray) – Initial particle positions, e.g. from landau_damping_ic() or two_stream_ic().
v0 (ndarray) – Initial particle velocities.
L (float) – Domain length (periodic).
ng (int) – Number of grid points.
dt (float) – Time step.
steps (int) – Number of steps to advance.
qm (float) – Charge-to-mass ratio in normalized units.
n0 (float) – Equilibrium number density.

Return type:

Returns:

dict – {"t": ndarray of shape (steps,), "field_energy": ndarray of shape (steps,), "x": final positions, "v": final velocities}.

See also

pic_step: The single-step update repeated here.
landau_damping_rate: Analytic decay rate to compare the field-energy history against.

Examples

>>> import numpy as np
>>> x0, v0 = landau_damping_ic(2000, L=4 * np.pi, k_mode=0.5, alpha=0.01, v_th=1.0, seed=0)
>>> result = pic_simulate(x0, v0, L=4 * np.pi, ng=32, dt=0.1, steps=20)
>>> result["field_energy"].shape
(20,)
>>> bool(np.all(np.isfinite(result["field_energy"])))
True

physicskit.plasma.pic_step(x, v, L, ng, dt, qm=-1.0, n0=1.0)[source]#

Advance the electrostatic PIC system one leapfrog step.

Deposits the electron density, solves for the self-consistent field against a uniform neutralizing ion background, gathers the field back onto the particles, and kicks/drifts them – one full cycle of the deposit-solve-gather-push loop at the heart of every PIC code. Velocities are staggered a half step behind positions (standard leapfrog); see pic_simulate() for a driver that initializes that offset correctly.

Parameters:

x (ndarray) – Particle positions in \([0, L)\).
v (ndarray) – Particle velocities, staggered a half step behind x.
L (float) – Domain length (periodic).
ng (int) – Number of grid points.
dt (float) – Time step.
qm (float) – Charge-to-mass ratio in normalized units (-1.0 for electrons with a fixed, uniform ion background of density n0).
n0 (float) – Equilibrium number density.

Return type:

Returns:

x_new, v_new (ndarray) – Updated positions and velocities.
field_energy (float) – \(\int E^2/2\,dx\), evaluated at the field used for this step’s kick.

See also

pic_simulate: Repeated application of this step with correct leapfrog initialization.

Examples

>>> import numpy as np
>>> x = np.linspace(0, 10, 50, endpoint=False)
>>> v = np.zeros(50)
>>> x_new, v_new, fe = pic_step(x, v, L=10.0, ng=32, dt=0.1)
>>> x_new.shape, v_new.shape
((50,), (50,))
>>> bool(np.isfinite(fe))
True

physicskit.plasma.plasma_frequency(n, q=1.602176634e-19, m=9.1093837139e-31)[source]#

Species plasma frequency \(\omega_p = \sqrt{nq^2/(\varepsilon_0 m)}\).

The natural oscillation frequency of a species displaced from quasineutrality: the restoring electric field it builds up is proportional to the displacement, making every unmagnetized plasma a harmonic oscillator at this frequency – the very phenomenon Langmuir identified in 1928.

Parameters:

n (float) – Number density in m^-3.
q (float) – Species charge magnitude in Coulombs.
m (float) – Species mass in kilograms.

Return type:

Returns:

float – Angular plasma frequency in rad/s.

Examples

>>> round(float(plasma_frequency(1e19)) / 1e9, 3)
178.399

physicskit.plasma.plot_cma_diagram(X, Y, ax=None)[source]#

Scatter a set of plasma states on log-log Clemmow-Mullaly-Allis (CMA) diagram axes.

Parameters:

X (ndarray) – CMA coordinates, as from physicskit.plasma.waves.cma_coordinates().
Y (ndarray) – CMA coordinates, as from physicskit.plasma.waves.cma_coordinates().
ax (matplotlib.axes.Axes, optional) – Axes to draw into; a new figure is created if omitted.

Returns:

fig (matplotlib.figure.Figure)
ax (matplotlib.axes.Axes)

Examples

>>> import numpy as np
>>> from physicskit.plasma.waves import cma_coordinates
>>> omega = np.linspace(0.1, 5.0, 30)
>>> X, Y = cma_coordinates(omega, wpe=2.0, wce=3.0)
>>> fig, ax = plot_cma_diagram(X, Y)
>>> isinstance(fig, plt.Figure)
True

physicskit.plasma.plot_drift_trajectory(pos_hist, ax=None)[source]#

Plot the guiding-center drift path as seen from above (the x-y plane).

Parameters:

pos_hist (ndarray) – Position history; only the first two components are used.
ax (matplotlib.axes.Axes, optional) – Axes to draw into; a new figure is created if omitted.

Returns:

fig (matplotlib.figure.Figure)
ax (matplotlib.axes.Axes)

See also

plot_particle_orbit_3d: The full 3D gyro-orbit this drift path averages over.

Examples

>>> import numpy as np
>>> from physicskit.plasma.single_particle import boris_integrate, QE, MP
>>> E = np.array([0.0, 1e3, 0.0])
>>> B = np.array([0.0, 0.0, 1.0])
>>> pos_hist, vel_hist = boris_integrate(np.zeros(3), np.zeros(3), QE, MP, E, B, 1e-10, steps=500)
>>> fig, ax = plot_drift_trajectory(pos_hist)
>>> isinstance(fig, plt.Figure)
True

physicskit.plasma.plot_field_energy_history(t, field_energy, ax=None)[source]#

Semilog plot of electrostatic field energy vs. time, showing Landau damping (or two-stream growth).

Parameters:

t (ndarray) – Time array, as from physicskit.plasma.kinetic.pic_simulate().
field_energy (ndarray) – Field energy at each time, same shape as t.
ax (matplotlib.axes.Axes, optional) – Axes to draw into; a new figure is created if omitted.

Returns:

fig (matplotlib.figure.Figure)
ax (matplotlib.axes.Axes)

Examples

>>> import numpy as np
>>> from physicskit.plasma.kinetic import landau_damping_ic, pic_simulate
>>> x0, v0 = landau_damping_ic(2000, L=4 * np.pi, k_mode=0.5, alpha=0.05, v_th=1.0, seed=0)
>>> result = pic_simulate(x0, v0, L=4 * np.pi, ng=32, dt=0.1, steps=30)
>>> fig, ax = plot_field_energy_history(result["t"], result["field_energy"])
>>> isinstance(fig, plt.Figure)
True

physicskit.plasma.plot_flux_surfaces(R, Z, psi, ax=None, levels=20)[source]#

Contour-plot poloidal flux surfaces \(\psi(R, Z)\) of a toroidal equilibrium.

Parameters:

R (ndarray) – Cylindrical coordinate grids.
Z (ndarray) – Cylindrical coordinate grids.
psi (ndarray) – Poloidal flux, as from physicskit.plasma.mhd.solve_grad_shafranov().
ax (matplotlib.axes.Axes, optional) – Axes to draw into; a new figure is created if omitted.
levels (int) – Number of contour levels.

Returns:

fig (matplotlib.figure.Figure)
ax (matplotlib.axes.Axes)

Examples

>>> import numpy as np
>>> from physicskit.plasma.mhd import solve_grad_shafranov
>>> R = np.linspace(0.5, 1.5, 41)
>>> Z = np.linspace(-0.5, 0.5, 41)
>>> psi = solve_grad_shafranov(R, Z, c1=1.0, c2=-2.0)
>>> fig, ax = plot_flux_surfaces(R, Z, psi)
>>> isinstance(fig, plt.Figure)
True

physicskit.plasma.plot_particle_orbit_3d(pos_hist, ax=None)[source]#

Plot a charged particle’s 3D trajectory (e.g. Boris-pusher gyro-orbit).

Parameters:

pos_hist (ndarray) – Position history, as from physicskit.plasma.single_particle.boris_integrate().
ax (matplotlib.axes.Axes3D, optional) – 3D axes to draw into; a new figure is created if omitted.

Returns:

fig (matplotlib.figure.Figure)
ax (matplotlib.axes.Axes3D)

Examples

>>> import numpy as np
>>> from physicskit.plasma.single_particle import boris_integrate, QE, MP
>>> pos_hist, vel_hist = boris_integrate(
...     np.zeros(3), np.array([1e5, 0.0, 0.0]), QE, MP, np.zeros(3), np.array([0.0, 0.0, 1.0]), 1e-10, steps=200
... )
>>> fig, ax = plot_particle_orbit_3d(pos_hist)
>>> isinstance(fig, plt.Figure)
True

physicskit.plasma.plot_phase_space(x, v, ax=None, bins=64)[source]#

Heatmap the particle-in-cell phase-space density \(f(x, v)\) from a particle snapshot.

Parameters:

x (ndarray) – Particle positions and velocities, e.g. from physicskit.plasma.kinetic.pic_simulate().
v (ndarray) – Particle positions and velocities, e.g. from physicskit.plasma.kinetic.pic_simulate().
ax (matplotlib.axes.Axes, optional) – Axes to draw into; a new figure is created if omitted.
bins (int) – Number of histogram bins along each axis.

Returns:

fig (matplotlib.figure.Figure)
ax (matplotlib.axes.Axes)

See also

plot_phase_space_interactive: An interactive Plotly scatter of the same data.

Examples

>>> import numpy as np
>>> from physicskit.plasma.kinetic import landau_damping_ic
>>> x, v = landau_damping_ic(2000, L=4 * np.pi, k_mode=0.5, alpha=0.1, v_th=1.0, seed=0)
>>> fig, ax = plot_phase_space(x, v)
>>> isinstance(fig, plt.Figure)
True

physicskit.plasma.plot_phase_space_interactive(x, v)[source]#

Interactive Plotly scatter of a particle-in-cell phase-space snapshot.

Unlike plot_phase_space(), points remain individually identifiable under pan and zoom – useful for inspecting fine structure like phase-space vortices in a developed two-stream instability.

Parameters:

x (ndarray) – Particle positions and velocities.
v (ndarray) – Particle positions and velocities.

Returns:

plotly.graph_objects.Figure

See also

plot_phase_space: The static Matplotlib heatmap equivalent.

Examples

>>> import numpy as np
>>> from physicskit.plasma.kinetic import two_stream_ic
>>> x, v = two_stream_ic(500, L=10.0, v_drift=3.0, v_th=0.5, seed=0)
>>> fig = plot_phase_space_interactive(x, v)
>>> isinstance(fig, go.Figure)
True

physicskit.plasma.plot_q_profile(r, q, ax=None)[source]#

Plot the tokamak safety factor \(q(r)\) against minor radius.

Parameters:

r (ndarray) – Minor-radius coordinate.
q (ndarray) – Safety factor, e.g. from physicskit.plasma.mhd.safety_factor_large_aspect_ratio() evaluated at each r.
ax (matplotlib.axes.Axes, optional) – Axes to draw into; a new figure is created if omitted.

Returns:

fig (matplotlib.figure.Figure)
ax (matplotlib.axes.Axes)

Examples

>>> import numpy as np
>>> from physicskit.plasma.mhd import safety_factor_large_aspect_ratio
>>> r = np.linspace(0.05, 0.5, 20)
>>> q = np.array([safety_factor_large_aspect_ratio(ri, R0=1.0, Bt=2.0, Bp=0.2) for ri in r])
>>> fig, ax = plot_q_profile(r, q)
>>> isinstance(fig, plt.Figure)
True

physicskit.plasma.reconnection_field_from_flux(psi, dx, dy)[source]#

Recover the in-plane magnetic field \((B_x, B_y) = (\partial_y\psi, -\partial_x\psi)\) from the flux function.

Parameters:

psi (ndarray) – Flux function, periodic in both directions.
dx (float) – Grid spacing along x and y.
dy (float) – Grid spacing along x and y.

Return type:

Returns:

Bx, By (ndarray, shape (nx, ny)) – In-plane magnetic field components.

Examples

>>> import numpy as np
>>> psi = reconnection_harris_ic(32, 32, Lx=20.0, Ly=20.0, sheet_width=1.0, perturbation_amplitude=0.0)
>>> Bx, By = reconnection_field_from_flux(psi, dx=20.0 / 32, dy=20.0 / 32)
>>> bool(np.all(Bx[:, 16] * Bx[:, 15] <= 0))
True

physicskit.plasma.reconnection_harris_ic(nx, ny, Lx, Ly, sheet_width, perturbation_amplitude, k_modes=1)[source]#

Flux-function initial condition: a Harris current sheet with an X-point-seeding ripple.

The unperturbed Harris sheet \(\psi_0(y) = -B_0 L\ln\cosh(y/L)\) gives an antiparallel reconnecting field \(B_x = \partial_y\psi_0 = -B_0\tanh(y/L)\) (with the reconnection_field_from_flux() convention \(B_x=\partial_y\psi\)) that reverses sign across \(y=0\) – the classic current-sheet configuration in which reconnection is normally studied. Adding a small ripple \(\psi_1 = \epsilon\cos(k x)\,\mathrm{sech}^2(y/L)\), localized to the sheet and periodic in \(x\), is the standard tearing-mode-like seed: it breaks the translational symmetry along the sheet just enough to create one X-point/O-point pair per wavelength, without which the sheet would stay a stationary equilibrium forever regardless of resistivity.

Parameters:

nx (int) – Grid points along x (periodic, along the sheet) and y (across the sheet).
ny (int) – Grid points along x (periodic, along the sheet) and y (across the sheet).
Lx (float) – Domain size in x and y; the domain spans [0, Lx) x [-Ly/2, Ly/2).
Ly (float) – Domain size in x and y; the domain spans [0, Lx) x [-Ly/2, Ly/2).
sheet_width (float) – Current-sheet half-thickness \(L\) (with \(B_0=1\) in these normalized units).
perturbation_amplitude (float) – Amplitude \(\epsilon\) of the seeding ripple.
k_modes (int) – Number of full wavelengths of the ripple fit into Lx.

Return type:

Returns:

ndarray, shape (nx, ny) – Flux function \(\psi(x, y)\).

See also

simulate_reconnection: Evolves this initial condition forward in time.

Examples

>>> psi0 = reconnection_harris_ic(64, 64, Lx=20.0, Ly=20.0, sheet_width=1.0, perturbation_amplitude=0.1)
>>> psi0.shape
(64, 64)

physicskit.plasma.reconnection_inflow_velocity(X, Y, v0, Lx, Ly)[source]#

Prescribed stagnation-point inflow/outflow velocity field around the reconnection X-point.

A simple incompressible stagnation flow, \(v_x=+v_0 x'/L_x\), \(v_y=-v_0 y'/L_y\) (with \(x', y'\) measured from the nearest X-point at the domain center), that converges onto the X-point from above and below and squirts back out sideways as reconnection outflow “jets” – exactly the inflow/outflow pattern Sweet-Parker and Petschek reconnection assume, here imposed directly rather than derived from a momentum equation (see the module docstring).

Parameters:

X (ndarray) – Coordinate grids, as built by reconnection_harris_ic().
Y (ndarray) – Coordinate grids, as built by reconnection_harris_ic().
v0 (float) – Characteristic inflow speed.
Lx (float) – Domain size in x and y, used to center the stagnation point.
Ly (float) – Domain size in x and y, used to center the stagnation point.

Return type:

Returns:

vx, vy (ndarray) – Velocity components, same shape as X.

Examples

>>> import numpy as np
>>> X, Y = np.meshgrid(np.linspace(0, 20, 8, endpoint=False), np.linspace(-10, 10, 8, endpoint=False), indexing="ij")
>>> vx, vy = reconnection_inflow_velocity(X, Y, v0=0.1, Lx=20.0, Ly=20.0)
>>> vx.shape
(8, 8)

physicskit.plasma.resistive_diffusion_time(L, eta)[source]#

Resistive diffusion time \(\tau_\eta = L^2/\eta\) for magnetic field to decay through a length \(L\).

Parameters:

L (float) – Length scale in meters.
eta (float) – Magnetic diffusivity in m^2/s.

Return type:

Returns:

float – Diffusion time in seconds.

Examples

>>> resistive_diffusion_time(L=1.0, eta=1.0)
1.0

physicskit.plasma.rl_parameters(S, D)[source]#

Right- and left-hand Stix parameters \(R = S+D\), \(L = S-D\).

\(R\) and \(L\) are the dielectric response seen by a purely right- or left-hand circularly polarized wave propagating exactly along \(\mathbf{B}_0\); the R-wave resonates at the electron cyclotron frequency and the L-wave at the ion cyclotron frequency.

Parameters:

S (float) – Stix parameters from stix_parameters().
D (float) – Stix parameters from stix_parameters().

Return type:

Returns:

R, L (float) – Right- and left-hand dielectric parameters.

Examples

>>> rl_parameters(S=1.0, D=0.5)
(1.5, 0.5)

physicskit.plasma.safety_factor_large_aspect_ratio(r, R0, Bt, Bp)[source]#

Tokamak safety factor \(q \approx rB_t/(R_0 B_p)\) in the large-aspect-ratio approximation.

Counts how many times a field line winds the long way (toroidally) around the torus for each time it winds the short way (poloidally). Field lines with rational \(q = m/n\) close on themselves after \(n\) toroidal transits and are resonant surfaces for magnetic-island-forming instabilities; \(q=1\) in particular marks the sawtooth-unstable region at a tokamak’s core.

Parameters:

r (float) – Minor-radius coordinate of the flux surface, in meters.
R0 (float) – Major radius of the torus, in meters.
Bt (float) – Toroidal field strength at the flux surface, in Tesla.
Bp (float) – Poloidal field strength at the flux surface, in Tesla.

Return type:

Returns:

float – Safety factor (dimensionless).

Examples

>>> round(safety_factor_large_aspect_ratio(r=0.3, R0=1.0, Bt=2.0, Bp=0.2), 10)
3.0

physicskit.plasma.simulate_alfven_wave(By0, vy0, x, dt, steps, B0, rho0, mu0=1.0)[source]#

Time-step the linearized 1D ideal-MHD Alfven-wave equations with a pseudo-spectral RK4 scheme.

Advances the coupled transverse induction and momentum equations

\[\partial_t B_y = B_0\,\partial_x v_y, \qquad \rho_0\,\partial_t v_y = \frac{B_0}{\mu_0}\,\partial_x B_y,\]

which combine into the non-dispersive wave equation \(\partial_t^2 B_y = v_A^2\,\partial_x^2 B_y\) with \(v_A=B_0/\sqrt{\mu_0\rho_0}\) (physicskit.plasma.mhd.alfven_speed()). Spatial derivatives are evaluated exactly via FFT (as elsewhere in this package’s periodic-domain solvers) and advanced in time with classical RK4; a uniform background field and density are assumed throughout (linear, ideal, cold-background MHD – no thermal pressure term enters a purely transverse, incompressible perturbation like this one).

Unlike physicskit.plasma.mhd.alfven_speed() (SI units), this solver defaults to normalized units (mu0=1.0, matching order-unity B0/rho0) since with the true SI \(\mu_0\approx1.26\times10^{-6}\) a Tesla-scale field gives an Alfven speed of order \(10^5\)-\(10^6\) m/s, requiring correspondingly tiny time steps to satisfy the CFL bound below; pass mu0=physicskit.plasma.mhd.MU0 explicitly for SI-consistent parameters and scale dt accordingly.

Parameters:

By0 (ndarray) – Initial transverse field and velocity perturbations, e.g. from alfven_wave_pulse_ic().
vy0 (ndarray) – Initial transverse field and velocity perturbations, e.g. from alfven_wave_pulse_ic().
x (ndarray) – Uniformly spaced periodic spatial grid.
dt (float) – Time step; the CFL condition \(v_A\,dt \le dx\) should be respected for the explicit RK4 stepping to remain stable.
steps (int) – Number of RK4 steps to advance.
B0 (float) – Background field magnitude.
rho0 (float) – Background mass density.
mu0 (float) – Vacuum permeability (normalized units by default; see above).

Return type:

Returns:

dict – {"By": final transverse field, "vy": final transverse velocity}.

See also

alfven_wave_pulse_ic: Builds the initial condition consumed here.
physicskit.plasma.mhd.alfven_speed: The propagation speed this recovers.

Examples

>>> import numpy as np
>>> x = np.linspace(-20, 20, 256, endpoint=False)
>>> By0, vy0 = alfven_wave_pulse_ic(x, x0=0.0, width=1.0, amplitude=0.1)
>>> result = simulate_alfven_wave(By0, vy0, x, dt=0.002, steps=2000, B0=1.0, rho0=1.0)
>>> result["By"].shape
(256,)
>>> bool(np.isfinite(result["By"]).all())
True

physicskit.plasma.simulate_hasegawa_mima(phi0, dt, steps, length, nu=0.03)[source]#

Time-step the Hasegawa-Mima drift-wave-turbulence equation with pseudo-spectral RK4.

Starting from a potential field \(\phi_0\) (e.g. small-amplitude noise from drift_wave_noise_ic()), converts to the potential vorticity \(q=\nabla^2\phi-\phi\) and advances it one step at a time: the non-stiff advection and linear drift-wave terms (hasegawa_mima_rhs()) with explicit RK4, Strang-split around an exact diffusive decay \(\hat{q}\mathrel{*}=e^{-\nu k^2 dt/2}\) applied before and after – the same splitting idea physicskit.fields.solitons’s KdV solver uses for its stiff dispersive term, applied here because a small Laplacian-type dissipation \(-\nu\nabla^2 q\) (not part of the ideal Hasegawa-Mima equation) is needed to drain enstrophy piling up at the grid scale once the flow turns turbulent, and integrating it exactly avoids the explicit-RK4 stability restriction that treating it as an ordinary right-hand-side term would impose. Reproduces the characteristic Hasegawa-Mima phenomenology of small-scale drift waves nonlinearly steepening and merging into a field of long-lived coherent vortices (“turbulent blobs”) that then dominates the cross-field transport, in analogy to two-dimensional Navier-Stokes turbulence’s inverse energy cascade.

Parameters:

phi0 (ndarray) – Initial electrostatic potential on a doubly periodic [0, length)^2 domain.
dt (float) – Time step.
steps (int) – Number of RK4 steps to advance.
length (float) – Physical domain size.
nu (float) – Dissipation coefficient (see above); purely a numerical stabilizer, not part of the ideal physics. Too small a value lets the (undealiased) pseudo-spectral nonlinear term blow up once the flow turns turbulent; the default was chosen to remain stable for the noise amplitudes and grid resolutions used throughout this module’s examples and tests.

Return type:

Returns:

dict – {"phi": final potential, "q": final potential vorticity}.

See also

drift_wave_noise_ic: Builds a typical initial condition consumed here.
hasegawa_mima_rhs: The single-step right-hand side repeated here.

Examples

>>> import numpy as np
>>> phi0 = drift_wave_noise_ic(48, 2 * np.pi, amplitude=0.01, seed=0)
>>> result = simulate_hasegawa_mima(phi0, dt=0.02, steps=20, length=2 * np.pi)
>>> result["phi"].shape
(48, 48)
>>> bool(np.all(np.isfinite(result["phi"])))
True

physicskit.plasma.simulate_reconnection(psi0, eta, v0, dt, steps, Lx, Ly)[source]#

Time-step the kinematic resistive induction equation \(\partial_t\psi = \eta\nabla^2\psi - \mathbf{v}\cdot\nabla\psi\).

A doubly periodic, second-order central-difference, explicit (forward Euler in time) finite-difference solve. With a prescribed velocity field rather than one obtained from solving the momentum equation (see the module docstring), this reduces reconnection to pure flux transport-and-diffusion: the inflow (reconnection_inflow_velocity()) advects oppositely directed flux into the X-point, where resistive diffusion (the \(\eta\nabla^2\psi\) term) is the only mechanism that can actually break and reconnect field lines – exactly Faraday’s law with an Ohmic (rather than ideal) Ohm’s law, restricted to the kinematic (fixed-flow) limit.

Explicit forward-Euler time-stepping is only conditionally stable: the diffusive term requires \(\eta\,dt \lesssim \tfrac{1}{4}\min(dx,dy)^2\) and the advective term requires the Courant condition \(v_0\,dt \lesssim \min(dx, dy)\). Both are the caller’s responsibility to satisfy by choosing dt appropriately; unstable combinations manifest as exponentially growing grid-scale noise.

Parameters:

psi0 (ndarray) – Initial flux function, e.g. from reconnection_harris_ic().
eta (float) – Resistivity (magnetic diffusivity).
v0 (float) – Characteristic inflow speed for reconnection_inflow_velocity().
dt (float) – Time step.
steps (int) – Number of forward-Euler steps to advance.
Lx (float) – Domain size in x and y.
Ly (float) – Domain size in x and y.

Return type:

Returns:

dict – {"psi": final flux function, "Bx": ..., "By": ...}.

See also

reconnection_harris_ic: Builds the initial condition consumed here.
physicskit.plasma.mhd.sweet_parker_rate: The steady-state reconnection rate this kinematic model is a simplified, time-dependent analogue of.

Examples

>>> import numpy as np
>>> psi0 = reconnection_harris_ic(48, 48, Lx=20.0, Ly=20.0, sheet_width=1.0, perturbation_amplitude=0.2)
>>> result = simulate_reconnection(psi0, eta=0.02, v0=0.05, dt=0.02, steps=50, Lx=20.0, Ly=20.0)
>>> result["psi"].shape
(48, 48)
>>> bool(np.all(np.isfinite(result["psi"])))
True

physicskit.plasma.simulate_wakefield_acceleration(x0, v0, q, m, E0, k, v_phase, dt, steps)[source]#

Integrate a test charge’s 1D motion in the prescribed traveling wakefield, with the Boris pusher.

Calls physicskit.plasma.single_particle.boris_push() once per step with \(\mathbf{B}=0\) and \(\mathbf{E}=(E_z(x,t),0,0)\) re-evaluated at the particle’s current position and time each step – unlike physicskit.plasma.single_particle.boris_integrate(), which assumes a field uniform in space and time, this wakefield varies in both, so the field must be recomputed every step rather than baked into a single Numba-compiled loop. The pusher itself remains exactly energy-conserving in a pure magnetic field and exact-leapfrog in a pure electric field; it is otherwise an ordinary (non-relativistic) classical integrator, so results should be read qualitatively once the particle’s speed approaches a meaningful fraction of v_phase (real wakefield acceleration is an intrinsically relativistic problem).

Parameters:

x0 (float) – Initial position along the wave.
v0 (float) – Initial longitudinal velocity.
q (float) – Particle charge.
m (float) – Particle mass.
E0 (float) – Peak wakefield amplitude, from wakefield_e_field().
k (float) – Wakefield wavenumber.
v_phase (float) – Wakefield phase velocity.
dt (float) – Time step.
steps (int) – Number of steps to advance.

Return type:

Returns:

dict – {"t": ndarray shape (steps+1,), "x": position history, "v": longitudinal velocity history, "kinetic_energy": :math:`\tfrac12 m v^2` history}.

See also

wakefield_e_field: The prescribed field this particle rides.

Examples

A particle launched exactly at the wave’s phase velocity, sitting on the accelerating part of the field, gains kinetic energy:

>>> import numpy as np
>>> result = simulate_wakefield_acceleration(
...     x0=0.0, v0=0.9, q=1.0, m=1.0, E0=0.05, k=1.0, v_phase=1.0, dt=0.01, steps=2000
... )
>>> bool(result["kinetic_energy"][-1] > result["kinetic_energy"][0])
True

physicskit.plasma.simulate_weibel_filamentation(x, t, wpe, temperature_anisotropy, n_modes=12, k_max_factor=0.9, seed=0)[source]#

Reduced quasi-linear Weibel filamentation model: a spectrum of independently growing current modes.

Seeds n_modes transverse-current Fourier modes with small random amplitudes and phases, spanning wavenumbers up to k_max_factor times the linear cutoff \(k_{max}\) of weibel_growth_rate(), and evolves each mode’s amplitude independently as \(a_k(t) = a_k(0)\,e^{\gamma(k)t}\) – exact linear theory, valid until the fastest-growing filaments approach nonlinear saturation (order-unity current perturbation), at which point real filamentation departs from this superposition (mode coupling, filament merging) that this reduced model does not capture. This stands in for a full 2D-in-velocity electromagnetic PIC simulation; see the module docstring for why.

Parameters:

x (ndarray) – Spatial grid (periodic).
t (ndarray) – Times at which to evaluate the current pattern.
wpe (float) – Electron plasma frequency.
temperature_anisotropy (float) – The ratio \(T_\perp/T_\parallel\); must exceed 1.
n_modes (int) – Number of Fourier modes to seed.
k_max_factor (float) – Fraction of the linear cutoff wavenumber up to which modes are seeded (kept below 1 so every seeded mode is genuinely unstable).
seed (int) – Random seed for the initial mode amplitudes and phases.

Return type:

Returns:

ndarray, shape (nt, nx) – Transverse current density \(J_y(x, t)\), in units where the initial per-mode amplitude is order unity (i.e. a normalized, not absolute, current).

See also

weibel_growth_rate: The per-mode growth rate used here.

Examples

>>> import numpy as np
>>> x = np.linspace(0, 20.0, 128, endpoint=False)
>>> t = np.linspace(0, 5.0, 10)
>>> J = simulate_weibel_filamentation(x, t, wpe=1.0, temperature_anisotropy=4.0, n_modes=6, seed=0)
>>> J.shape
(10, 128)
>>> bool(np.std(J[-1]) > np.std(J[0]))
True

physicskit.plasma.solovev_particular_solution(R, Z, c1, c2)[source]#

Closed-form particular solution of the linear (Solov’ev) Grad-Shafranov equation.

When the source term is linear in \(R^2\) – i.e. \(p'(\psi) = \text{const}\) and \(FF'(\psi) = \text{const}\), so \(\Delta^*\psi = c_1 R^2 + c_2\) – the Grad-Shafranov operator \(\Delta^*\psi = \partial_R^2\psi - R^{-1}\partial_R\psi + \partial_Z^2\psi\) admits the exact polynomial solution \(\psi_p = \tfrac{c_1}{8}R^4 + \tfrac{c_2}{2}Z^2\) (Solov’ev, 1968); adding any solution of the homogeneous equation \(\Delta^*\psi=0\) shapes the boundary into a D-shaped or elongated cross-section without affecting the pressure and current profile. Used here to validate solve_grad_shafranov() against an exact answer, and to supply consistent Dirichlet boundary data for it.

Parameters:

R (ndarray) – Cylindrical coordinates (broadcastable), in meters.
Z (ndarray) – Cylindrical coordinates (broadcastable), in meters.
c1 (float) – Coefficients of the linear source term \(\Delta^*\psi = c_1 R^2 + c_2\).
c2 (float) – Coefficients of the linear source term \(\Delta^*\psi = c_1 R^2 + c_2\).

Return type:

Returns:

ndarray – Poloidal flux \(\psi(R, Z)\), same shape as R/Z.

See also

grad_shafranov_rhs_solovev: The corresponding source term.
solve_grad_shafranov: Numerical solver validated against this solution.

Examples

>>> import numpy as np
>>> solovev_particular_solution(np.array([1.0]), np.array([0.0]), c1=1.0, c2=-2.0)
array([0.125])

physicskit.plasma.solve_grad_shafranov(R, Z, c1, c2, omega=1.8, max_iter=4000)[source]#

Solve the axisymmetric Grad-Shafranov equation by successive over-relaxation (SOR).

Finite-differences the elliptic operator \(\Delta^*\psi = \partial_R^2\psi - R^{-1}\partial_R\psi + \partial_Z^2\psi\) on a rectangular \((R, Z)\) grid and relaxes it toward the linear (Solov’ev) source \(c_1 R^2 + c_2\), using Dirichlet boundary data taken from the exact solovev_particular_solution() – so the interior solution this converges to is known analytically and can be checked directly, rather than only visually.

Parameters:

R (ndarray) – Major-radius grid points, in meters; must be strictly positive (the operator is singular at \(R=0\)).
Z (ndarray) – Vertical grid points, in meters.
c1 (float) – Coefficients of the linear source term, as in grad_shafranov_rhs_solovev().
c2 (float) – Coefficients of the linear source term, as in grad_shafranov_rhs_solovev().
omega (float) – SOR relaxation parameter, \(1 < \omega < 2\).
max_iter (int) – Number of relaxation sweeps.

Return type:

Returns:

ndarray, shape (nr, nz) – Poloidal flux \(\psi(R, Z)\) on the grid.

See also

solovev_particular_solution: The exact solution this converges to.
safety_factor_large_aspect_ratio: A downstream equilibrium diagnostic.

Examples

>>> import numpy as np
>>> R = np.linspace(0.5, 1.5, 41)
>>> Z = np.linspace(-0.5, 0.5, 41)
>>> c1, c2 = 1.0, -2.0
>>> psi = solve_grad_shafranov(R, Z, c1, c2)
>>> RR, ZZ = np.meshgrid(R, Z, indexing="ij")
>>> psi_exact = solovev_particular_solution(RR, ZZ, c1, c2)
>>> bool(np.max(np.abs(psi - psi_exact)) < 1e-3)
True

physicskit.plasma.solve_poisson_1d(rho, L)[source]#

Solve the 1D periodic Poisson equation \(dE/dx = \rho\) (normalized \(\varepsilon_0=1\)) via FFT.

Parameters:

rho (ndarray) – Net charge density on the grid (e.g. ion background minus deposited electron density). Its mean is discarded, since a uniform charge density has no periodic solution and physically should integrate to zero net charge in the box.
L (float) – Domain length (periodic).

Return type:

Returns:

ndarray, shape (ng,) – Electric field on the grid.

See also

deposit_number_density: Supplies the density this solves for.
interpolate_field: Gathers this field back onto the particles.

Examples

>>> import numpy as np
>>> ng = 64
>>> L = 2 * np.pi
>>> x_grid = np.linspace(0, L, ng, endpoint=False)
>>> rho = np.sin(x_grid)
>>> E = solve_poisson_1d(rho, L)
>>> E_exact = -np.cos(x_grid)
>>> bool(np.max(np.abs(E - E_exact)) < 1e-10)
True

physicskit.plasma.sound_speed(gamma, p, rho)[source]#

Ordinary adiabatic sound speed \(c_s = \sqrt{\gamma p/\rho}\).

Parameters:

gamma (float) – Adiabatic index (5/3 for an ideal monatomic gas).
p (float) – Pressure in Pa.
rho (float) – Mass density in kg/m^3.

Return type:

Returns:

float – Sound speed in m/s.

Examples

>>> round(float(sound_speed(gamma=5 / 3, p=1.0, rho=1e-6)), 2)
1290.99

physicskit.plasma.stix_parameters(omega, B, species)[source]#

Stix cold-plasma dielectric tensor components \(S\), \(D\), \(P\).

Summing each species’ contribution to the linearized fluid response gives the Hermitian dielectric tensor

\[\begin{split}\mathbf{K} = \begin{pmatrix} S & -iD & 0 \\ iD & S & 0 \\ 0 & 0 & P \end{pmatrix}, \qquad S = 1 - \sum_s \frac{\omega_{ps}^2}{\omega^2-\omega_{cs}^2}, \quad D = \sum_s \frac{\omega_{cs}}{\omega}\frac{\omega_{ps}^2}{\omega^2-\omega_{cs}^2}, \quad P = 1 - \sum_s \frac{\omega_{ps}^2}{\omega^2},\end{split}\]

with each species’ signed cyclotron frequency \(\omega_{cs}=q_sB/m_s\) entering \(D\) with its own sign – the source of the circular-polarization asymmetry between the R- and L-waves.

Parameters:

omega (float) – Wave angular frequency in rad/s.
B (float) – Background magnetic field magnitude in Tesla.
species (sequence of (float, float, float)) – (n, q, m) for each plasma species: number density in m^-3, signed charge in Coulombs, mass in kilograms.

Return type:

Returns:

S, D, P (float) – Stix dielectric tensor components (dimensionless).

See also

cold_plasma_dispersion: Solves the dispersion relation built from these.
rl_parameters: The right/left-hand combinations \(R=S+D\), \(L=S-D\).

Examples

>>> electrons = (1e19, -QE, ME)
>>> ions = (1e19, QE, MP)
>>> S, D, P = stix_parameters(omega=2e9, B=1.0, species=[electrons, ions])
>>> round(S, 4), round(D, 4), round(P, 4)
(-2.3143, 90.6954, -7959.8516)

physicskit.plasma.sweet_parker_layer_width(L, S)[source]#

Sweet-Parker current-sheet thickness \(\delta = L/\sqrt{S}\).

Parameters:

L (float) – Current-sheet length in meters.
S (float) – Lundquist number, from lundquist_number().

Return type:

Returns:

float – Current-sheet thickness in meters.

Examples

>>> round(float(sweet_parker_layer_width(L=1e7, S=1e6)), 4)
10000.0

physicskit.plasma.sweet_parker_rate(S)[source]#

Sweet-Parker reconnection rate \(v_{in}/v_A = S^{-1/2}\).

Sweet and Parker’s 1957/1958 model treats reconnection as steady inflow through a long, thin resistive current sheet of aspect ratio \(\delta/L \sim S^{-1/2}\); mass conservation through that narrow sheet throttles the inflow (and hence the whole reconnection process) to the same \(S^{-1/2}\) scaling. For solar-flare-scale Lundquist numbers (\(S\sim10^{12}\)) this predicts reconnection millions of times too slow to explain observed flare energy-release times – the puzzle petschek_rate() was proposed to resolve.

Parameters:: S (float) – Lundquist number, from lundquist_number().
Return type:: float
Returns:: float – Dimensionless reconnection rate \(v_{in}/v_A\).

See also

sweet_parker_layer_width: The current-sheet thickness behind this rate.
petschek_rate: The faster, X-point reconnection alternative.

Examples

>>> round(float(sweet_parker_rate(S=1e6)), 6)
0.001

physicskit.plasma.two_stream_ic(n_particles, L, v_drift, v_th, seed=0)[source]#

Initial condition for the two-stream instability: two counter-streaming Maxwellian beams.

Splits the particles into two equal populations drifting at \(\pm v_{drift}\), each with thermal spread \(v_{th}\), and seeds the fastest-growing long-wavelength mode with a small density ripple. When \(v_{drift}\) exceeds the thermal spread by enough to make the combined velocity distribution doubly-peaked, the positive-slope region between the two peaks violates the (kinetic) Penrose stability criterion and the ripple grows exponentially, eventually rolling the two beams up into a single phase-space vortex.

Parameters:

n_particles (int) – Number of particles (split evenly between the two beams).
L (float) – Domain length (periodic).
v_drift (float) – Drift speed of each beam (beams move at \(+v_{drift}\) and \(-v_{drift}\)).
v_th (float) – Thermal spread of each beam.
seed (int) – Random seed.

Return type:

Returns:

x, v (ndarray, shape (n_particles,)) – Particle positions and velocities.

See also

landau_damping_ic: The companion (stable) single-beam initial condition.

Examples

>>> import numpy as np
>>> x, v = two_stream_ic(1000, L=10.0, v_drift=3.0, v_th=0.5, seed=0)
>>> x.shape, v.shape
((1000,), (1000,))
>>> bool(np.mean(v) < 0.5)
True

physicskit.plasma.wakefield_e_field(x, t, E0, k, v_phase)[source]#

Prescribed traveling-wave longitudinal wakefield \(E_z(x,t) = E_0\cos[k(x - v_{ph}t)]\).

A rigid sinusoidal accelerating structure moving at the fixed phase velocity v_phase, standing in for the self-consistent plasma wakefield a real driver beam or laser wake would generate (see the module docstring for the simplification this represents).

Parameters:

x (ndarray) – Position(s) at which to evaluate the field.
t (float) – Time.
E0 (float) – Peak field amplitude.
k (float) – Wakefield wavenumber.
v_phase (float) – Phase velocity of the traveling wave.

Return type:

Returns:

ndarray or float – \(E_z(x, t)\), same shape as x.

See also

simulate_wakefield_acceleration: Injects a test particle into this field.

Examples

>>> round(float(wakefield_e_field(0.0, t=0.0, E0=1.0, k=1.0, v_phase=1.0)), 6)
1.0

physicskit.plasma.weibel_fastest_growing_mode(wpe, temperature_anisotropy, c=299792458.0)[source]#

Wavenumber and growth rate of the fastest-growing Weibel mode.

For the growth-rate law of weibel_growth_rate(), \(\gamma(k)^2\) is maximized at \(k=0\) (long wavelength, where the stabilizing field-tension term vanishes), giving \(\gamma_{max}=\omega_{pe}\sqrt{T_\perp/T_\parallel-1}\) – the growth rate used to set the overall filament-growth timescale in simulate_weibel_filamentation(), even though the dominant filament spacing in a real (finite-seed) simulation is set by which finite-\(k\) modes happen to be seeded with the largest initial amplitude, since every \(k\) below cutoff grows, just more slowly farther from \(k=0\).

Parameters:

wpe (float) – Electron plasma frequency.
temperature_anisotropy (float) – The ratio \(T_\perp/T_\parallel\). At or below 1 the plasma is Weibel-stable and (0.0, 0.0) is returned.
c (float) – Speed of light. Unused, since the maximum sits at \(k=0\) where the \(k^2c^2\) term vanishes; kept for signature symmetry with weibel_growth_rate().

Return type:

Returns:

k_max_growth (float) – Wavenumber of maximum growth (always 0.0).
gamma_max (float) – The corresponding growth rate \(\gamma_{max}\).

Examples

>>> k0, gamma_max = weibel_fastest_growing_mode(wpe=1.0, temperature_anisotropy=4.0)
>>> k0
0.0
>>> round(gamma_max, 6)
1.732051

physicskit.plasma.weibel_growth_rate(k, wpe, temperature_anisotropy, c=299792458.0)[source]#

Linear growth rate of the (non-relativistic, cold-parallel-limit) Weibel/filamentation instability.

For an electron distribution anisotropic between the direction of wave propagation (temperature \(T_\parallel\), taken cold here) and the transverse direction (temperature \(T_\perp\)), the electromagnetic mode with wavevector \(\mathbf{k}\) transverse to the anisotropy axis is purely growing (Weibel, 1959) with

\[\gamma(k)^2 = \omega_{pe}^2\left(\frac{T_\perp}{T_\parallel}-1\right) - k^2 c^2,\]

unstable for \(T_\perp > T_\parallel\) and only up to the cutoff \(k_{max} = (\omega_{pe}/c)\sqrt{T_\perp/T_\parallel - 1}\) – beyond which the field’s magnetic tension (the \(k^2c^2\) term) overcomes the free energy the anisotropy supplies. Growing transverse magnetic perturbations pinch counter-streaming electron sub-populations into spatially separated current channels: the filaments this instability is named for.

Parameters:

k (ndarray) – Wavenumber(s), transverse to the anisotropy axis.
wpe (float) – Electron plasma frequency.
temperature_anisotropy (float) – The ratio \(T_\perp/T_\parallel\); must exceed 1 for any instability.
c (float) – Speed of light, in the same units k and wpe combine to give (SI by default).

Return type:

Returns:

ndarray or float – Growth rate \(\gamma(k)\), clipped to zero for stable (real or imaginary-\(\gamma^2<0\)) wavenumbers.

See also

weibel_fastest_growing_mode: The wavenumber maximizing this growth rate.
simulate_weibel_filamentation: Superposes modes growing at this rate.

Examples

>>> import numpy as np
>>> round(float(weibel_growth_rate(k=0.0, wpe=1.0, temperature_anisotropy=4.0)), 6)
1.732051
>>> bool(weibel_growth_rate(k=1e10, wpe=1.0, temperature_anisotropy=4.0) == 0.0)
True

physicskit.plasma.whistler_dispersion(omega, wpe, wce)[source]#

Electron-only, parallel-propagation whistler-mode refractive index \(n^2 = \omega_{pe}^2/[\omega(\omega_{ce}-\omega)]\).

The low-frequency (\(\omega \ll \omega_{ce}\)), right-hand circularly polarized branch of the R-wave, with ion motion neglected. Because \(n^2\) grows without bound as \(\omega\to\omega_{ce}\), whistlers slow down sharply as they approach the electron cyclotron frequency – the falling-tone “whistle” heard in radio receivers after a lightning stroke launches a broadband pulse that disperses along Earth’s field lines, with the highest frequencies (closest to \(\omega_{ce}\)) arriving last.

Parameters:

omega (float) – Wave angular frequency in rad/s, with \(\omega \ll \omega_{ce}\).
wpe (float) – Electron plasma frequency in rad/s, from plasma_frequency().
wce (float) – Electron cyclotron frequency in rad/s (positive).

Return type:

Returns:

float – Squared refractive index \(n^2\).

See also

cold_plasma_dispersion: The general dispersion relation this is a limit of.

Examples

Well below the electron cyclotron frequency, it agrees closely with the exact electron-only R-wave root of cold_plasma_dispersion():

>>> wpe, wce = 1.784e11, 1.759e11
>>> omega = 0.001 * wce
>>> S, D, P = stix_parameters(omega, B=1.0, species=[(wpe ** 2 * EPS0 * ME / QE ** 2, -QE, ME)])
>>> R, _ = rl_parameters(S, D)
>>> n2_approx = whistler_dispersion(omega, wpe, wce)
>>> round(abs(n2_approx - R) / R, 2)
0.0

Single-particle motion in electromagnetic fields: the Boris pusher and guiding-center drifts.

Two complementary pictures of charged-particle motion in a magnetized plasma live here. The full-orbit picture integrates the exact Lorentz force \(m\dot{\mathbf{v}} = q(\mathbf{E} + \mathbf{v}\times\mathbf{B})\) with boris_push() / boris_integrate(), Jay Boris’s 1970 energy-conserving leapfrog scheme – still the workhorse integrator of every particle-in-cell (PIC) plasma code. The guiding-center picture instead averages over the fast gyration and tracks only the slow drift of the orbit’s center: exb_drift(), grad_b_drift(), and curvature_drift() give that drift velocity directly, while magnetic_moment() and mirror_force() capture the adiabatic invariant \(\mu = m v_\perp^2 / 2B\) responsible for magnetic-mirror confinement.

The Boris kernels are compiled with Numba for particle-in-cell-scale performance; every other function here is a closed-form NumPy expression.

physicskit.plasma.single_particle.ME = 9.1093837139e-31#: Electron mass, in kilograms.

physicskit.plasma.single_particle.MP = 1.67262192595e-27#: Proton mass, in kilograms.

physicskit.plasma.single_particle.QE = 1.602176634e-19#: Elementary charge, in Coulombs.

physicskit.plasma.single_particle.boris_integrate(pos0, vel0, q, m, E, B, dt, steps)[source]#

Integrate a charged particle’s trajectory in uniform E and B fields with the Boris pusher.

Repeatedly applies boris_push() inside a single Numba-compiled loop, avoiding Python-level overhead per step – the throughput this buys is what makes the Boris scheme practical for particle-in-cell codes tracking millions of particles.

Parameters:

pos0 (ndarray) – Initial position in meters.
vel0 (ndarray) – Initial velocity in m/s.
q (float) – Particle charge in Coulombs.
m (float) – Particle mass in kilograms.
E (ndarray) – Uniform electric field in V/m.
B (ndarray) – Uniform magnetic field in Tesla.
dt (float) – Time step in seconds.
steps (int) – Number of steps to advance.

Return type:

Returns:

pos_hist (ndarray, shape (steps + 1, 3)) – Position at every step, including the initial condition.
vel_hist (ndarray, shape (steps + 1, 3)) – Velocity at every step, including the initial condition.

See also

boris_push: The single-step update repeated here.
exb_drift: The closed-form drift velocity this trajectory averages to when a non-zero E is present.

Examples

Kinetic energy is conserved to machine precision over many gyrations in a pure magnetic field:

>>> import numpy as np
>>> B = np.array([0.0, 0.0, 1.0])
>>> E = np.zeros(3)
>>> omega_c = cyclotron_frequency(QE, MP, 1.0)
>>> dt = (2 * np.pi / omega_c) / 200
>>> pos_hist, vel_hist = boris_integrate(
...     np.zeros(3), np.array([1e5, 0.0, 0.0]), QE, MP, E, B, dt, steps=2000
... )
>>> speeds = np.linalg.norm(vel_hist, axis=1)
>>> bool(np.max(np.abs(speeds - 1e5)) < 1e-6)
True

physicskit.plasma.single_particle.boris_push(pos, vel, q, m, E, B, dt)[source]#

Advance a charged particle one leapfrog step with the Boris integrator.

The Boris (1970) scheme splits each step into an electric half-acceleration, an exact rotation about \(\mathbf{B}\) (via the Boris “E-cross-B” trick, avoiding any explicit trigonometry), and a second electric half-acceleration. The rotation step preserves speed exactly in a pure magnetic field, which is why the scheme conserves kinetic energy over arbitrarily many gyro-orbits where a naive forward-Euler or even RK4 update would spiral outward.

Parameters:

pos (ndarray) – Particle position \([x, y, z]\) in meters.
vel (ndarray) – Particle velocity \([v_x, v_y, v_z]\) in m/s.
q (float) – Particle charge in Coulombs.
m (float) – Particle mass in kilograms.
E (ndarray) – Electric field \([E_x, E_y, E_z]\) in V/m.
B (ndarray) – Magnetic field \([B_x, B_y, B_z]\) in Tesla.
dt (float) – Time step in seconds.

Return type:

Returns:

pos_new, vel_new (ndarray, shape (3,)) – Position and velocity after one step of size dt.

See also

boris_integrate: Repeated application of this step over many time steps.

Examples

A proton launched perpendicular to a uniform field gyrates without gaining or losing speed:

>>> import numpy as np
>>> r0 = np.array([0.0, 0.0, 0.0])
>>> v0 = np.array([1e5, 0.0, 0.0])
>>> E = np.zeros(3)
>>> B = np.array([0.0, 0.0, 1.0])
>>> r1, v1 = boris_push(r0, v0, q=QE, m=MP, E=E, B=B, dt=1e-9)
>>> round(float(np.linalg.norm(v1)), 6)
100000.0

physicskit.plasma.single_particle.curvature_drift(v_par, q, m, B, R_c)[source]#

The curvature drift \(\mathbf{v}_R = \dfrac{m v_\parallel^2}{q}\dfrac{\mathbf{R}_c\times\mathbf{B}}{R_c^2 B^2}\).

A particle streaming along a curved field line feels a centrifugal force in the guiding-center frame, which – crossed with \(\mathbf{B}\) – produces a drift perpendicular to the plane of curvature. In a low-beta toroidal equilibrium this combines with grad_b_drift() (curvature and field-strength gradients point the same way when \(\nabla\times\mathbf{B}=0\) in vacuum) into the single \(\nabla B\) + curvature drift responsible for charge separation and the resulting E-cross-B rotation in a tokamak.

Parameters:

v_par (float) – Speed parallel to the field, in m/s.
q (float) – Particle charge in Coulombs (signed).
m (float) – Particle mass in kilograms.
B (ndarray) – Local magnetic field vector in Tesla.
R_c (ndarray) – Radius-of-curvature vector, pointing from the field line’s local center of curvature to the particle, magnitude \(R_c\) in meters.

Return type:

Returns:

ndarray, shape (3,) – Drift velocity in m/s.

Examples

>>> import numpy as np
>>> B = np.array([0.0, 0.0, 1.0])
>>> R_c = np.array([1.0, 0.0, 0.0])
>>> drift = curvature_drift(1e5, QE, MP, B, R_c)
>>> round(float(drift[1]), 2)
-104.4

physicskit.plasma.single_particle.cyclotron_frequency(q, m, B)[source]#

Angular cyclotron (gyration) frequency \(\omega_c = qB/m\).

Parameters:

q (float) – Particle charge in Coulombs (signed).
m (float) – Particle mass in kilograms.
B (float) – Magnetic field magnitude in Tesla.

Return type:

Returns:

float – Signed angular gyrofrequency in rad/s (negative for negatively charged particles, reflecting the opposite sense of rotation).

See also

larmor_radius: The orbit radius set by this frequency.

Examples

A proton in a 1 T field gyrates about 15.3 MHz:

>>> round(cyclotron_frequency(QE, MP, 1.0) / (2 * 3.141592653589793) / 1e6, 2)
15.25

physicskit.plasma.single_particle.exb_drift(E, B)[source]#

The \(\mathbf{E}\times\mathbf{B}\) drift velocity \(\mathbf{v}_E = (\mathbf{E}\times\mathbf{B})/B^2\).

Unlike every other guiding-center drift, this one is independent of charge, mass, and energy: electrons and ions drift together at the same velocity, so \(\mathbf{E}\times\mathbf{B}\) drift carries no net current and instead advects the whole plasma as a fluid – the drift underlying tokamak radial-electric-field rotation and magnetospheric convection alike.

Parameters:

E (ndarray) – Electric field in V/m.
B (ndarray) – Magnetic field in Tesla.

Return type:

Returns:

ndarray, shape (3,) – Drift velocity in m/s.

Examples

>>> import numpy as np
>>> E = np.array([0.0, 1e3, 0.0])
>>> B = np.array([0.0, 0.0, 1.0])
>>> exb_drift(E, B)
array([1000.,    0.,    0.])

physicskit.plasma.single_particle.grad_b_drift(v_perp, q, m, B, grad_B)[source]#

The grad-B drift \(\mathbf{v}_{\nabla B} = \dfrac{m v_\perp^2}{2qB^3}(\mathbf{B}\times\nabla B)\).

A particle gyrating in a field whose magnitude varies across the orbit sees a tighter turn (smaller Larmor radius) on the strong-field side than the weak-field side, so the orbit fails to close and the guiding center creeps sideways. Because the drift is inversely proportional to charge, electrons and ions drift in opposite directions, producing a net current – the origin of the ring current in planetary magnetospheres.

Parameters:

v_perp (float) – Speed perpendicular to the field, in m/s.
q (float) – Particle charge in Coulombs (signed).
m (float) – Particle mass in kilograms.
B (ndarray) – Local magnetic field vector in Tesla.
grad_B (ndarray) – Gradient of the field magnitude, \(\nabla |\mathbf{B}|\), in Tesla/meter.

Return type:

Returns:

ndarray, shape (3,) – Drift velocity in m/s.

See also

curvature_drift: The companion drift from field-line curvature, which combines with this one in any real toroidal field.

Examples

>>> import numpy as np
>>> B = np.array([0.0, 0.0, 1.0])
>>> grad_B = np.array([0.1, 0.0, 0.0])
>>> drift = grad_b_drift(1e5, QE, MP, B, grad_B)
>>> round(float(drift[1]), 4)
5.2198

physicskit.plasma.single_particle.larmor_radius(v_perp, q, m, B)[source]#

Larmor (gyro) radius \(r_L = m v_\perp / (|q| B)\).

Parameters:

v_perp (float) – Speed perpendicular to the magnetic field, in m/s.
q (float) – Particle charge in Coulombs (sign is ignored).
m (float) – Particle mass in kilograms.
B (float) – Magnetic field magnitude in Tesla.

Return type:

Returns:

float – Gyroradius in meters.

Examples

>>> round(larmor_radius(1e5, QE, MP, 1.0), 6)
0.001044

physicskit.plasma.single_particle.magnetic_mirror_bounce(z0, v_par0, v_perp0, m, B_func, dz=1e-06, dt=1e-10, steps=20000)[source]#

Simulate 1D guiding-center bounce motion between the throats of a magnetic mirror.

Integrates \(m\dot{v}_\parallel = -\mu\, dB/dz\) with \(\mu\) fixed at its initial value (the adiabatic invariant), using a symmetric leapfrog step. A particle with too little pitch angle to reflect before reaching the mirror throat’s peak field instead falls into the loss cone and would be lost from confinement in a real device; mirror_force() supplies the underlying force law.

Parameters:

z0 (float) – Initial position along the field line, in meters.
v_par0 (float) – Initial parallel velocity, in m/s.
v_perp0 (float) – Initial perpendicular velocity, in m/s (sets \(\mu\) via magnetic_moment()).
m (float) – Particle mass in kilograms.
B_func (callable) – Field-strength profile B_func(z) -> float along the field line, in Tesla.
dz (float) – Finite-difference step used to evaluate \(dB/dz\), in meters.
dt (float) – Time step in seconds.
steps (int) – Number of leapfrog steps to advance.

Return type:

Returns:

z_hist (ndarray, shape (steps + 1,)) – Position along the field line at every step.
v_par_hist (ndarray, shape (steps + 1,)) – Parallel velocity at every step.

Examples

A particle launched from the mirror midplane with enough perpendicular energy reflects before reaching the throat, reversing the sign of its parallel velocity:

>>> import numpy as np
>>> B_func = lambda z: 1.0 + 4.0 * (z / 0.05) ** 2
>>> z_hist, v_par_hist = magnetic_mirror_bounce(
...     z0=0.0, v_par0=2e4, v_perp0=8e4, m=MP, B_func=B_func, steps=6000
... )
>>> bool(v_par_hist[0] > 0 and v_par_hist[-1] < 0)
True

physicskit.plasma.single_particle.magnetic_moment(v_perp, m, B)[source]#

The first adiabatic invariant \(\mu = m v_\perp^2 / (2B)\).

Conserved for a charged particle whose gyration is fast compared to any change in the field it sees, \(\mu\) acts as a magnetic “potential energy per unit field”: as the particle drifts into stronger \(B\), \(v_\perp\) must grow to keep \(\mu\) fixed, converting parallel kinetic energy into perpendicular kinetic energy. That conversion is the mechanism behind mirror_force().

Parameters:

v_perp (float) – Speed perpendicular to the magnetic field, in m/s.
m (float) – Particle mass in kilograms.
B (float) – Magnetic field magnitude in Tesla.

Return type:

Returns:

float – Magnetic moment in Joules/Tesla.

See also

mirror_force: The parallel force derived from this invariant.

Examples

>>> round(magnetic_moment(1e5, MP, 1.0) * 1e18, 4)
8.3631

physicskit.plasma.single_particle.mirror_force(mu, grad_B_parallel)[source]#

The mirror force \(F_\parallel = -\mu\, \partial B/\partial \ell\) along a field line.

As a particle’s guiding center moves into a region of stronger field (larger \(\partial B/\partial\ell\)), the conservation of magnetic_moment() forces \(v_\perp\) to grow at the expense of \(v_\parallel\); this is the reaction force decelerating the parallel motion. If the field is strong enough, \(v_\parallel\) reaches zero before the particle passes the throat and it reflects – magnetic mirror confinement, simulated end to end in magnetic_mirror_bounce().

Parameters:

mu (float) – Magnetic moment in J/T, from magnetic_moment().
grad_B_parallel (float) – Gradient of the field magnitude along the field line, in Tesla/meter.

Return type:

Returns:

float – Force in Newtons, directed to push the particle toward weaker field.

Examples

>>> mirror_force(mu=1e-17, grad_B_parallel=2.0)
-2e-17

Magnetohydrodynamics: wave speeds, toroidal equilibrium, and magnetic reconnection.

Treats the plasma as a single conducting fluid rather than a collection of orbiting particles. Three regimes are covered: the characteristic wave speeds of ideal MHD (alfven_speed(), sound_speed(), magnetosonic_speeds()); static toroidal equilibrium, obtained by solving the Grad-Shafranov equation \(\Delta^*\psi = -\mu_0 R^2 p'(\psi) - FF'(\psi)\) for the poloidal flux function \(\psi(R, Z)\) that balances pressure against the magnetic force in an axisymmetric device (solve_grad_shafranov(), safety_factor_large_aspect_ratio()); and resistive magnetic reconnection, where a thin current sheet lets field lines of opposite polarity break and reconnect, converting magnetic energy to heat and flow (sweet_parker_rate(), petschek_rate()).

physicskit.plasma.mhd.MU0 = 1.2566370614359173e-06#: Vacuum permeability, in H/m.

physicskit.plasma.mhd.alfven_speed(B, rho)[source]#

Alfven speed \(v_A = B/\sqrt{\mu_0 \rho}\), at which a perturbation propagates along tensioned field lines.

Hannes Alfven’s 1942 discovery that a magnetized, perfectly conducting fluid supports a transverse wave – field lines behaving like strings under tension \(B^2/\mu_0\), plucked by the inertia of the frozen-in plasma – founded MHD as a distinct discipline and earned the 1970 Nobel Prize.

Parameters:

B (float) – Magnetic field magnitude in Tesla.
rho (float) – Mass density in kg/m^3.

Return type:

Returns:

float – Alfven speed in m/s.

See also

magnetosonic_speeds: The compressive (fast/slow) counterparts of this purely magnetic wave.

Examples

>>> round(float(alfven_speed(B=1.0, rho=1e-6)), 2)
892062.06

physicskit.plasma.mhd.grad_shafranov_rhs_solovev(R, c1, c2)[source]#

Source term \(\Delta^*\psi = c_1 R^2 + c_2\) of the linear Solov’ev equilibrium.

Equivalent to \(-\mu_0 R^2 p'(\psi) - FF'(\psi)\) in the Grad-Shafranov equation for the special case of constant \(p'(\psi)\) and \(FF'(\psi)\).

Parameters:

R (ndarray) – Major-radius coordinate, in meters.
c1 (float) – Source coefficients (related to \(p'\) and \(FF'\)).
c2 (float) – Source coefficients (related to \(p'\) and \(FF'\)).

Return type:

Returns:

ndarray – The right-hand side \(\Delta^*\psi\), same shape as R.

Examples

>>> import numpy as np
>>> grad_shafranov_rhs_solovev(np.array([1.0, 2.0]), c1=1.0, c2=-2.0)
array([-1.,  2.])

physicskit.plasma.mhd.lundquist_number(L, vA, eta)[source]#

Lundquist number \(S = L v_A/\eta\), the magnetic Reynolds number built from the Alfven speed.

Parameters:

L (float) – Characteristic length scale (e.g. current-sheet length), in meters.
vA (float) – Alfven speed in m/s, from alfven_speed().
eta (float) – Magnetic diffusivity in m^2/s.

Return type:

Returns:

float – Lundquist number (dimensionless). Fusion and astrophysical plasmas typically have \(S \sim 10^{6}\) – \(10^{14}\).

See also

sweet_parker_rate: Reconnection rate scaling as \(S^{-1/2}\).

Examples

>>> lundquist_number(L=1.0, vA=1e6, eta=1.0)
1000000.0

physicskit.plasma.mhd.magnetosonic_speeds(vA, cs, theta)[source]#

Fast and slow magnetosonic phase speeds at propagation angle \(\theta\) to \(\mathbf{B}\).

The two compressive MHD normal modes solve \(v_{f,s}^2 = \tfrac{1}{2}\left[(v_A^2+c_s^2) \pm \sqrt{(v_A^2+c_s^2)^2 - 4v_A^2c_s^2\cos^2\theta}\right]\). At \(\theta=0\) (propagation along \(\mathbf{B}\)) they reduce to \(\max(v_A, c_s)\) and \(\min(v_A, c_s)\); at \(\theta=\pi/2\) the fast mode becomes the purely compressive \(\sqrt{v_A^2+c_s^2}\) and the slow mode vanishes, since a perpendicular perturbation cannot bend field lines that are already perpendicular to its wavevector.

Parameters:

vA (float) – Alfven speed in m/s, from alfven_speed().
cs (float) – Sound speed in m/s, from sound_speed().
theta (float) – Angle between the wavevector and \(\mathbf{B}\), in radians.

Return type:

Returns:

v_fast, v_slow (float) – Fast and slow magnetosonic phase speeds in m/s.

Examples

>>> import numpy as np
>>> vf, vs = magnetosonic_speeds(vA=892062.06, cs=1e5, theta=np.pi / 2)
>>> round(vf, 2)
897649.55
>>> round(vs, 2)
0.0

physicskit.plasma.mhd.petschek_rate(S)[source]#

Petschek reconnection rate \(v_{in}/v_A \approx \pi/(8\ln S)\).

Petschek (1964) showed that if the diffusion region shrinks to a small X-point rather than the full Sweet-Parker sheet length, four standing slow-mode shocks can carry most of the inflowing flux and energy conversion, giving a reconnection rate that falls only logarithmically with \(S\) instead of as \(S^{-1/2}\) – fast enough to plausibly explain solar flare and magnetospheric substorm timescales.

Parameters:: S (float) – Lundquist number, from lundquist_number().
Return type:: float
Returns:: float – Dimensionless reconnection rate \(v_{in}/v_A\).

See also

sweet_parker_rate: The slower, steady-sheet reconnection rate this improves on.

Examples

>>> round(float(petschek_rate(S=1e6)), 4)
0.0284

physicskit.plasma.mhd.resistive_diffusion_time(L, eta)[source]#

Resistive diffusion time \(\tau_\eta = L^2/\eta\) for magnetic field to decay through a length \(L\).

Parameters:

L (float) – Length scale in meters.
eta (float) – Magnetic diffusivity in m^2/s.

Return type:

Returns:

float – Diffusion time in seconds.

Examples

>>> resistive_diffusion_time(L=1.0, eta=1.0)
1.0

physicskit.plasma.mhd.safety_factor_large_aspect_ratio(r, R0, Bt, Bp)[source]#

Tokamak safety factor \(q \approx rB_t/(R_0 B_p)\) in the large-aspect-ratio approximation.

Counts how many times a field line winds the long way (toroidally) around the torus for each time it winds the short way (poloidally). Field lines with rational \(q = m/n\) close on themselves after \(n\) toroidal transits and are resonant surfaces for magnetic-island-forming instabilities; \(q=1\) in particular marks the sawtooth-unstable region at a tokamak’s core.

Parameters:

r (float) – Minor-radius coordinate of the flux surface, in meters.
R0 (float) – Major radius of the torus, in meters.
Bt (float) – Toroidal field strength at the flux surface, in Tesla.
Bp (float) – Poloidal field strength at the flux surface, in Tesla.

Return type:

Returns:

float – Safety factor (dimensionless).

Examples

>>> round(safety_factor_large_aspect_ratio(r=0.3, R0=1.0, Bt=2.0, Bp=0.2), 10)
3.0

physicskit.plasma.mhd.solovev_particular_solution(R, Z, c1, c2)[source]#

Closed-form particular solution of the linear (Solov’ev) Grad-Shafranov equation.

When the source term is linear in \(R^2\) – i.e. \(p'(\psi) = \text{const}\) and \(FF'(\psi) = \text{const}\), so \(\Delta^*\psi = c_1 R^2 + c_2\) – the Grad-Shafranov operator \(\Delta^*\psi = \partial_R^2\psi - R^{-1}\partial_R\psi + \partial_Z^2\psi\) admits the exact polynomial solution \(\psi_p = \tfrac{c_1}{8}R^4 + \tfrac{c_2}{2}Z^2\) (Solov’ev, 1968); adding any solution of the homogeneous equation \(\Delta^*\psi=0\) shapes the boundary into a D-shaped or elongated cross-section without affecting the pressure and current profile. Used here to validate solve_grad_shafranov() against an exact answer, and to supply consistent Dirichlet boundary data for it.

Parameters:

R (ndarray) – Cylindrical coordinates (broadcastable), in meters.
Z (ndarray) – Cylindrical coordinates (broadcastable), in meters.
c1 (float) – Coefficients of the linear source term \(\Delta^*\psi = c_1 R^2 + c_2\).
c2 (float) – Coefficients of the linear source term \(\Delta^*\psi = c_1 R^2 + c_2\).

Return type:

Returns:

ndarray – Poloidal flux \(\psi(R, Z)\), same shape as R/Z.

See also

grad_shafranov_rhs_solovev: The corresponding source term.
solve_grad_shafranov: Numerical solver validated against this solution.

Examples

>>> import numpy as np
>>> solovev_particular_solution(np.array([1.0]), np.array([0.0]), c1=1.0, c2=-2.0)
array([0.125])

physicskit.plasma.mhd.solve_grad_shafranov(R, Z, c1, c2, omega=1.8, max_iter=4000)[source]#

Solve the axisymmetric Grad-Shafranov equation by successive over-relaxation (SOR).

Finite-differences the elliptic operator \(\Delta^*\psi = \partial_R^2\psi - R^{-1}\partial_R\psi + \partial_Z^2\psi\) on a rectangular \((R, Z)\) grid and relaxes it toward the linear (Solov’ev) source \(c_1 R^2 + c_2\), using Dirichlet boundary data taken from the exact solovev_particular_solution() – so the interior solution this converges to is known analytically and can be checked directly, rather than only visually.

Parameters:

R (ndarray) – Major-radius grid points, in meters; must be strictly positive (the operator is singular at \(R=0\)).
Z (ndarray) – Vertical grid points, in meters.
c1 (float) – Coefficients of the linear source term, as in grad_shafranov_rhs_solovev().
c2 (float) – Coefficients of the linear source term, as in grad_shafranov_rhs_solovev().
omega (float) – SOR relaxation parameter, \(1 < \omega < 2\).
max_iter (int) – Number of relaxation sweeps.

Return type:

Returns:

ndarray, shape (nr, nz) – Poloidal flux \(\psi(R, Z)\) on the grid.

See also

solovev_particular_solution: The exact solution this converges to.
safety_factor_large_aspect_ratio: A downstream equilibrium diagnostic.

Examples

>>> import numpy as np
>>> R = np.linspace(0.5, 1.5, 41)
>>> Z = np.linspace(-0.5, 0.5, 41)
>>> c1, c2 = 1.0, -2.0
>>> psi = solve_grad_shafranov(R, Z, c1, c2)
>>> RR, ZZ = np.meshgrid(R, Z, indexing="ij")
>>> psi_exact = solovev_particular_solution(RR, ZZ, c1, c2)
>>> bool(np.max(np.abs(psi - psi_exact)) < 1e-3)
True

physicskit.plasma.mhd.sound_speed(gamma, p, rho)[source]#

Ordinary adiabatic sound speed \(c_s = \sqrt{\gamma p/\rho}\).

Parameters:

gamma (float) – Adiabatic index (5/3 for an ideal monatomic gas).
p (float) – Pressure in Pa.
rho (float) – Mass density in kg/m^3.

Return type:

Returns:

float – Sound speed in m/s.

Examples

>>> round(float(sound_speed(gamma=5 / 3, p=1.0, rho=1e-6)), 2)
1290.99

physicskit.plasma.mhd.sweet_parker_layer_width(L, S)[source]#

Sweet-Parker current-sheet thickness \(\delta = L/\sqrt{S}\).

Parameters:

L (float) – Current-sheet length in meters.
S (float) – Lundquist number, from lundquist_number().

Return type:

Returns:

float – Current-sheet thickness in meters.

Examples

>>> round(float(sweet_parker_layer_width(L=1e7, S=1e6)), 4)
10000.0

physicskit.plasma.mhd.sweet_parker_rate(S)[source]#

Sweet-Parker reconnection rate \(v_{in}/v_A = S^{-1/2}\).

Sweet and Parker’s 1957/1958 model treats reconnection as steady inflow through a long, thin resistive current sheet of aspect ratio \(\delta/L \sim S^{-1/2}\); mass conservation through that narrow sheet throttles the inflow (and hence the whole reconnection process) to the same \(S^{-1/2}\) scaling. For solar-flare-scale Lundquist numbers (\(S\sim10^{12}\)) this predicts reconnection millions of times too slow to explain observed flare energy-release times – the puzzle petschek_rate() was proposed to resolve.

Parameters:: S (float) – Lundquist number, from lundquist_number().
Return type:: float
Returns:: float – Dimensionless reconnection rate \(v_{in}/v_A\).

See also

sweet_parker_layer_width: The current-sheet thickness behind this rate.
petschek_rate: The faster, X-point reconnection alternative.

Examples

>>> round(float(sweet_parker_rate(S=1e6)), 6)
0.001

Cold plasma waves: the Stix dielectric tensor, dispersion relations, and CMA mapping.

Linearizes the multi-fluid cold-plasma equations around a uniform background threaded by a uniform \(\mathbf{B}_0\) to get the dielectric tensor components \(S\), \(D\), \(P\) (Stix, 1962) in stix_parameters(), then solves the general dispersion relation \(\mathbf{n}\times(\mathbf{n}\times\mathbf{E}) + \mathbf{K}\cdot\mathbf{E}=0\) for the refractive index \(n=ck/\omega\) at any propagation angle \(\theta\) to \(\mathbf{B}_0\) with cold_plasma_dispersion(). Every named mode – the parallel-propagating R- and L-waves, the perpendicular O- and X-modes, and the low-frequency whistler branch – is a special case of that one quartic in \(n^2\). cma_coordinates() maps a plasma state onto the two dimensionless axes of the Clemmow-Mullaly-Allis diagram that organizes all of them.

physicskit.plasma.waves.EPS0 = 8.8541878188e-12#: Vacuum permittivity, in F/m.

physicskit.plasma.waves.ME = 9.1093837139e-31#: Electron mass, in kilograms.

physicskit.plasma.waves.MP = 1.67262192595e-27#: Proton mass, in kilograms.

physicskit.plasma.waves.QE = 1.602176634e-19#: Elementary charge, in Coulombs.

physicskit.plasma.waves.alfven_wave_pulse_ic(x, x0, width, amplitude)[source]#

Initial condition for a transverse Alfven-wave pulse launched from rest.

A Gaussian transverse-field pulse \(B_y(x,0)=B_1 e^{-[(x-x_0)/w]^2}\) with the transverse velocity perturbation initially zero – a “plucked string” initial condition. Because the linearized ideal-MHD Alfven-wave equations (simulate_alfven_wave()) are the same non-dispersive wave equation a string obeys, a disturbance released from rest splits exactly in half and propagates as two identical, oppositely directed pulses at \(\pm v_A\) – exactly the field-line-plucking picture Alfven’s original 1942 analogy describes.

Parameters:

x (ndarray) – Spatial grid (periodic).
x0 (float) – Pulse center.
width (float) – Gaussian pulse width.
amplitude (float) – Peak transverse field perturbation \(B_1\).

Return type:

Returns:

By0, vy0 (ndarray) – Initial transverse magnetic field and velocity perturbations, same shape as x (vy0 identically zero).

See also

simulate_alfven_wave: Evolves this initial condition forward in time.

Examples

>>> import numpy as np
>>> x = np.linspace(-10, 10, 64, endpoint=False)
>>> By0, vy0 = alfven_wave_pulse_ic(x, x0=0.0, width=1.0, amplitude=0.1)
>>> bool(np.all(vy0 == 0.0))
True

physicskit.plasma.waves.cma_coordinates(omega, wpe, wce)[source]#

Dimensionless \((X, Y)\) coordinates of the Clemmow-Mullaly-Allis (CMA) diagram.

The CMA diagram partitions the plane \(X=\omega_{pe}^2/\omega^2\) (density axis) vs. \(Y=\omega_{ce}/\omega\) (field-strength axis) into regions of distinct wave topology – cutoffs, resonances, and the number and polarization of propagating modes – giving a single map of every cold-plasma wave regime from ordinary light waves (\(X, Y \to 0\)) to the Alfven wave and whistler branches (\(Y \gg 1\)).

Parameters:

omega (float) – Wave angular frequency in rad/s.
wpe (float) – Electron plasma frequency in rad/s.
wce (float) – Electron cyclotron frequency in rad/s.

Return type:

Returns:

X, Y (float) – CMA diagram coordinates.

Examples

>>> cma_coordinates(omega=1.0, wpe=2.0, wce=3.0)
(4.0, 3.0)

physicskit.plasma.waves.cold_plasma_dispersion(theta, S, D, P)[source]#

Solve the cold-plasma dispersion relation for the squared refractive index \(n^2\).

Substituting a plane wave into \(\mathbf{n}\times(\mathbf{n}\times\mathbf{E})+\mathbf{K}\cdot\mathbf{E}=0\) for propagation at angle \(\theta\) to \(\mathbf{B}_0\) gives the Appleton-Hartree biquadratic \(An^4 - Bn^2 + C = 0\) with

\[A = S\sin^2\theta + P\cos^2\theta, \quad B = RL\sin^2\theta + PS(1+\cos^2\theta), \quad C = PRL,\]

which reduces at \(\theta=0\) to the decoupled R- and L-waves (\(n^2=R\) or \(L\)) and at \(\theta=\pi/2\) to the O-mode (\(n^2=P\)) and X-mode (\(n^2=RL/S\)).

Parameters:

theta (float) – Angle between the wavevector and \(\mathbf{B}_0\), in radians.
S (float) – Stix parameters from stix_parameters().
D (float) – Stix parameters from stix_parameters().
P (float) – Stix parameters from stix_parameters().

Return type:

Returns:

n_sq_plus, n_sq_minus (float) – The two roots of the biquadratic (the two cold-plasma wave branches at this angle and frequency). A negative root means that branch is evanescent rather than propagating.

See also

stix_parameters: Supplies S, D, P.

Examples

Parallel propagation recovers the pure R- and L-wave refractive indices:

>>> import numpy as np
>>> S, D, P = -2.314262935091991, 90.69535621969021, -7959.851640342367
>>> n2_plus, n2_minus = cold_plasma_dispersion(theta=0.0, S=S, D=D, P=P)
>>> R, L = rl_parameters(S, D)
>>> bool(np.isclose(sorted([n2_plus, n2_minus]), sorted([R, L])).all())
True

physicskit.plasma.waves.ion_acoustic_soliton_evolve(u0, x, dt, steps)[source]#

Evolve an ion-acoustic KdV initial condition forward in time on a periodic domain.

A pseudo-spectral Strang-split scheme: the stiff linear dispersion \(\partial_\xi^3\) is advanced exactly in Fourier space, and the non-stiff nonlinear advection \(6uu_\xi\) with RK4 in between – structurally the standard splitting for any KdV-type equation, applied here directly to the ion-acoustic reduction rather than by importing a generic KdV solver, since the physical field this equation governs (a Debye-length-normalized density/potential perturbation moving at order the ion-sound speed) is specific to this module.

Parameters:

u0 (ndarray) – Initial field, sampled on the periodic grid x.
x (ndarray) – Uniformly spaced periodic spatial grid.
dt (float) – Time step.
steps (int) – Number of steps to advance.

Return type:

Returns:

ndarray – Field after steps * dt time units.

See also

ion_acoustic_soliton_profile: Exact traveling-wave solution this reproduces.

Examples

A soliton launched at speed \(c=4\) has advanced by very close to \(c\cdot(\text{steps}\cdot dt)\) and kept its amplitude, since KdV solitons propagate without changing shape:

>>> import numpy as np
>>> N, L = 512, 60.0
>>> x = np.linspace(-L / 2, L / 2, N, endpoint=False)
>>> u0 = ion_acoustic_soliton_profile(x, speed=4.0, x0=-15.0)
>>> u = ion_acoustic_soliton_evolve(u0, x, dt=0.0005, steps=4000)
>>> shift = x[np.argmax(u)] - x[np.argmax(u0)]
>>> bool(abs(shift - 4.0 * 4000 * 0.0005) < 0.5)
True
>>> bool(abs(u.max() - u0.max()) < 0.05)
True

physicskit.plasma.waves.ion_acoustic_soliton_profile(x, speed, x0=0.0)[source]#

Exact single-soliton solution of the ion-acoustic Korteweg-de Vries reduction.

The reductive-perturbation (Washimi & Taniuti, 1966) expansion of the cold-ion-fluid/Boltzmann-electron equations in the weakly nonlinear, weakly dispersive limit reduces the ion-acoustic wave problem to the KdV equation for the normalized density (or potential) perturbation \(u\) in a frame moving at the ion-sound speed; after the standard rescaling of the stretched coordinates it takes the canonical form \(u_t + 6uu_\xi + u_{\xi\xi\xi} = 0\) – the same normal form every weakly-dispersive weakly-nonlinear wave problem reduces to, with the ion-acoustic-specific physics (electron Boltzmann response supplying the nonlinearity, ion inertia and Debye-length dispersion supplying the \(\partial_\xi^3\) term) fixing only the physical unit conversions between \(u,\xi,t\) and density, position, and time in the ion-sound-speed frame. Its exact traveling-wave solution is a single soliton of speed \(c\) (in the stretched frame) and amplitude \(c/2\), propagating without change of shape – consistent with a real ion-acoustic soliton, whose speed always exceeds the linear sound speed by an amount set by its amplitude.

Parameters:

x (ndarray) – Spatial grid, in the stretched (ion-sound-speed) frame.
speed (float) – Soliton speed \(c\) in the stretched frame (amplitude \(=c/2\)); must be positive.
x0 (float) – Initial center position.

Return type:

Returns:

ndarray – \(u(x, 0) = \tfrac{c}{2}\,\mathrm{sech}^2\!\big(\tfrac{\sqrt{c}}{2}(x-x_0)\big)\).

See also

ion_acoustic_soliton_evolve: Propagate this (or any) initial condition forward in time.

Examples

>>> round(float(ion_acoustic_soliton_profile(0.0, speed=4.0)), 6)
2.0

physicskit.plasma.waves.plasma_frequency(n, q=1.602176634e-19, m=9.1093837139e-31)[source]#

Species plasma frequency \(\omega_p = \sqrt{nq^2/(\varepsilon_0 m)}\).

The natural oscillation frequency of a species displaced from quasineutrality: the restoring electric field it builds up is proportional to the displacement, making every unmagnetized plasma a harmonic oscillator at this frequency – the very phenomenon Langmuir identified in 1928.

Parameters:

n (float) – Number density in m^-3.
q (float) – Species charge magnitude in Coulombs.
m (float) – Species mass in kilograms.

Return type:

Returns:

float – Angular plasma frequency in rad/s.

Examples

>>> round(float(plasma_frequency(1e19)) / 1e9, 3)
178.399

physicskit.plasma.waves.rl_parameters(S, D)[source]#

Right- and left-hand Stix parameters \(R = S+D\), \(L = S-D\).

\(R\) and \(L\) are the dielectric response seen by a purely right- or left-hand circularly polarized wave propagating exactly along \(\mathbf{B}_0\); the R-wave resonates at the electron cyclotron frequency and the L-wave at the ion cyclotron frequency.

Parameters:

S (float) – Stix parameters from stix_parameters().
D (float) – Stix parameters from stix_parameters().

Return type:

Returns:

R, L (float) – Right- and left-hand dielectric parameters.

Examples

>>> rl_parameters(S=1.0, D=0.5)
(1.5, 0.5)

physicskit.plasma.waves.simulate_alfven_wave(By0, vy0, x, dt, steps, B0, rho0, mu0=1.0)[source]#

Time-step the linearized 1D ideal-MHD Alfven-wave equations with a pseudo-spectral RK4 scheme.

Advances the coupled transverse induction and momentum equations

\[\partial_t B_y = B_0\,\partial_x v_y, \qquad \rho_0\,\partial_t v_y = \frac{B_0}{\mu_0}\,\partial_x B_y,\]

which combine into the non-dispersive wave equation \(\partial_t^2 B_y = v_A^2\,\partial_x^2 B_y\) with \(v_A=B_0/\sqrt{\mu_0\rho_0}\) (physicskit.plasma.mhd.alfven_speed()). Spatial derivatives are evaluated exactly via FFT (as elsewhere in this package’s periodic-domain solvers) and advanced in time with classical RK4; a uniform background field and density are assumed throughout (linear, ideal, cold-background MHD – no thermal pressure term enters a purely transverse, incompressible perturbation like this one).

Unlike physicskit.plasma.mhd.alfven_speed() (SI units), this solver defaults to normalized units (mu0=1.0, matching order-unity B0/rho0) since with the true SI \(\mu_0\approx1.26\times10^{-6}\) a Tesla-scale field gives an Alfven speed of order \(10^5\)-\(10^6\) m/s, requiring correspondingly tiny time steps to satisfy the CFL bound below; pass mu0=physicskit.plasma.mhd.MU0 explicitly for SI-consistent parameters and scale dt accordingly.

Parameters:

By0 (ndarray) – Initial transverse field and velocity perturbations, e.g. from alfven_wave_pulse_ic().
vy0 (ndarray) – Initial transverse field and velocity perturbations, e.g. from alfven_wave_pulse_ic().
x (ndarray) – Uniformly spaced periodic spatial grid.
dt (float) – Time step; the CFL condition \(v_A\,dt \le dx\) should be respected for the explicit RK4 stepping to remain stable.
steps (int) – Number of RK4 steps to advance.
B0 (float) – Background field magnitude.
rho0 (float) – Background mass density.
mu0 (float) – Vacuum permeability (normalized units by default; see above).

Return type:

Returns:

dict – {"By": final transverse field, "vy": final transverse velocity}.

See also

alfven_wave_pulse_ic: Builds the initial condition consumed here.
physicskit.plasma.mhd.alfven_speed: The propagation speed this recovers.

Examples

>>> import numpy as np
>>> x = np.linspace(-20, 20, 256, endpoint=False)
>>> By0, vy0 = alfven_wave_pulse_ic(x, x0=0.0, width=1.0, amplitude=0.1)
>>> result = simulate_alfven_wave(By0, vy0, x, dt=0.002, steps=2000, B0=1.0, rho0=1.0)
>>> result["By"].shape
(256,)
>>> bool(np.isfinite(result["By"]).all())
True

physicskit.plasma.waves.stix_parameters(omega, B, species)[source]#

Stix cold-plasma dielectric tensor components \(S\), \(D\), \(P\).

Summing each species’ contribution to the linearized fluid response gives the Hermitian dielectric tensor

\[\begin{split}\mathbf{K} = \begin{pmatrix} S & -iD & 0 \\ iD & S & 0 \\ 0 & 0 & P \end{pmatrix}, \qquad S = 1 - \sum_s \frac{\omega_{ps}^2}{\omega^2-\omega_{cs}^2}, \quad D = \sum_s \frac{\omega_{cs}}{\omega}\frac{\omega_{ps}^2}{\omega^2-\omega_{cs}^2}, \quad P = 1 - \sum_s \frac{\omega_{ps}^2}{\omega^2},\end{split}\]

with each species’ signed cyclotron frequency \(\omega_{cs}=q_sB/m_s\) entering \(D\) with its own sign – the source of the circular-polarization asymmetry between the R- and L-waves.

Parameters:

omega (float) – Wave angular frequency in rad/s.
B (float) – Background magnetic field magnitude in Tesla.
species (sequence of (float, float, float)) – (n, q, m) for each plasma species: number density in m^-3, signed charge in Coulombs, mass in kilograms.

Return type:

Returns:

S, D, P (float) – Stix dielectric tensor components (dimensionless).

See also

cold_plasma_dispersion: Solves the dispersion relation built from these.
rl_parameters: The right/left-hand combinations \(R=S+D\), \(L=S-D\).

Examples

>>> electrons = (1e19, -QE, ME)
>>> ions = (1e19, QE, MP)
>>> S, D, P = stix_parameters(omega=2e9, B=1.0, species=[electrons, ions])
>>> round(S, 4), round(D, 4), round(P, 4)
(-2.3143, 90.6954, -7959.8516)

physicskit.plasma.waves.whistler_dispersion(omega, wpe, wce)[source]#

Electron-only, parallel-propagation whistler-mode refractive index \(n^2 = \omega_{pe}^2/[\omega(\omega_{ce}-\omega)]\).

The low-frequency (\(\omega \ll \omega_{ce}\)), right-hand circularly polarized branch of the R-wave, with ion motion neglected. Because \(n^2\) grows without bound as \(\omega\to\omega_{ce}\), whistlers slow down sharply as they approach the electron cyclotron frequency – the falling-tone “whistle” heard in radio receivers after a lightning stroke launches a broadband pulse that disperses along Earth’s field lines, with the highest frequencies (closest to \(\omega_{ce}\)) arriving last.

Parameters:

omega (float) – Wave angular frequency in rad/s, with \(\omega \ll \omega_{ce}\).
wpe (float) – Electron plasma frequency in rad/s, from plasma_frequency().
wce (float) – Electron cyclotron frequency in rad/s (positive).

Return type:

Returns:

float – Squared refractive index \(n^2\).

See also

cold_plasma_dispersion: The general dispersion relation this is a limit of.

Examples

Well below the electron cyclotron frequency, it agrees closely with the exact electron-only R-wave root of cold_plasma_dispersion():

>>> wpe, wce = 1.784e11, 1.759e11
>>> omega = 0.001 * wce
>>> S, D, P = stix_parameters(omega, B=1.0, species=[(wpe ** 2 * EPS0 * ME / QE ** 2, -QE, ME)])
>>> R, _ = rl_parameters(S, D)
>>> n2_approx = whistler_dispersion(omega, wpe, wce)
>>> round(abs(n2_approx - R) / R, 2)
0.0

Electrostatic particle-in-cell (PIC) solution of the 1D1V Vlasov-Poisson system.

Rather than discretizing the distribution function \(f(x, v, t)\) on a phase-space grid, a PIC code samples it with a finite set of computational “super-particles” and lets them stream along exact single-particle orbits, recovering the self-consistent field by depositing their charge onto a spatial grid and solving Poisson’s equation there each step. This is Vlasov’s collisionless equation \(\partial_t f + v\,\partial_x f - (e/m)E\,\partial_v f = 0\) solved by the method of characteristics: each particle is one characteristic. The whole pipeline – deposit_number_density(), solve_poisson_1d(), interpolate_field(), pic_step() – reproduces collisionless (Landau) damping and the two-stream instability without ever assuming a collision operator, exactly as the underlying kinetic theory predicts.

Units throughout are the standard PIC-normalized units: electron charge \(e=1\), mass \(m_e=1\), vacuum permittivity \(\varepsilon_0=1\), and equilibrium density \(n_0=1\), so that the electron plasma frequency \(\omega_{pe}=1\) and velocities are in units of the thermal speed \(v_{th}\). Ions are a fixed, uniform, charge-neutralizing background (infinite mass limit).

The charge-deposit and field-interpolation kernels – called once per particle per step – are compiled with Numba, since they dominate the cost of every PIC time step.

physicskit.plasma.kinetic.deposit_number_density(x, L, ng, n0=1.0)[source]#

Deposit particle positions onto a grid as a number density, via cloud-in-cell (CIC) weighting.

Each particle represents a “cloud” of physical charge spanning one grid cell, split linearly between its two nearest grid points – the standard first-order PIC weighting scheme, chosen because it is exact for a uniform density and (unlike nearest-grid-point deposit) produces a smooth, differentiable force with no self-force discontinuities as particles cross cell boundaries.

Parameters:

x (ndarray) – Particle positions in \([0, L)\).
L (float) – Domain length (periodic).
ng (int) – Number of grid points.
n0 (float) – Equilibrium number density (sets each particle’s statistical weight, \(n_0 L / n_{particles}\)).

Return type:

Returns:

ndarray, shape (ng,) – Number density on the grid.

See also

interpolate_field: The companion gather operation.

Examples

Total deposited charge exactly equals the physical charge represented, regardless of how the particles are distributed:

>>> import numpy as np
>>> x = np.array([0.1, 2.4, 4.9, 7.7])
>>> rho = deposit_number_density(x, L=10.0, ng=20, n0=2.0)
>>> dx = 10.0 / 20
>>> round(float(np.sum(rho) * dx), 8)
20.0

physicskit.plasma.kinetic.interpolate_field(x, field_grid, L)[source]#

Interpolate a grid-defined field to particle positions, via cloud-in-cell (CIC) weighting.

The gather step dual to deposit_number_density(): using the same linear weights for both deposit and gather is what makes the PIC method momentum-conserving (no self-force on an isolated particle).

Parameters:

x (ndarray) – Particle positions in \([0, L)\).
field_grid (ndarray) – Field values on the grid, e.g. from solve_poisson_1d().
L (float) – Domain length (periodic).

Return type:

Returns:

ndarray, shape (n_particles,) – Field value at each particle’s position.

Examples

A particle sitting exactly on a grid node picks up that node’s value:

>>> import numpy as np
>>> ng = 8
>>> field_grid = np.array([1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0])
>>> x = np.array([2 * (10.0 / ng)])
>>> interpolate_field(x, field_grid, L=10.0)
array([3.])

physicskit.plasma.kinetic.landau_damping_ic(n_particles, L, k_mode, alpha, v_th, seed=0)[source]#

Quiet-start initial condition seeding a single-mode density perturbation for a Landau damping test.

Displaces an otherwise uniform particle load by \(\delta x = (\alpha/k)\sin(kx_0)\), which – since particle number is conserved through the displacement’s Jacobian – produces exactly the density perturbation \(n(x) \approx n_0(1 - \alpha\cos(kx))\) used in the classic linear Landau damping test problem, without the sampling noise a random density draw would add on top of the intended signal.

Parameters:

n_particles (int) – Number of particles.
L (float) – Domain length (periodic), in normalized length units. Choosing \(L = 2\pi/k_{mode}\) fits exactly one wavelength.
k_mode (float) – Wavenumber of the seeded perturbation.
alpha (float) – Perturbation amplitude (\(\alpha \ll 1\) for the linear regime).
v_th (float) – Thermal speed of the background Maxwellian.
seed (int) – Random seed for the velocity sampling.

Return type:

Returns:

x, v (ndarray, shape (n_particles,)) – Particle positions (in \([0, L)\)) and velocities.

See also

pic_simulate: Evolve this initial condition forward in time.
landau_damping_rate: The analytic decay rate this test is checked against.

Examples

>>> import numpy as np
>>> x, v = landau_damping_ic(4, L=4.0, k_mode=2 * 3.141592653589793 / 4.0, alpha=0.1, v_th=1.0, seed=0)
>>> bool(np.all((x >= 0) & (x < 4.0)))
True
>>> x.shape, v.shape
((4,), (4,))

physicskit.plasma.kinetic.landau_damping_rate(k, v_th, omega_pe=1.0)[source]#

Analytic linear Landau damping rate for a Maxwellian electron plasma.

The classic weak-damping result (Landau, 1946) for a wave of wavenumber \(k\) on a Maxwellian of thermal speed \(v_{th}\), with Debye length \(\lambda_D = v_{th}/\omega_{pe}\):

\[\gamma = -\omega_{pe}\sqrt{\frac{\pi}{8}}\, \frac{1}{(k\lambda_D)^3}\, \exp\!\left(-\frac{1}{2(k\lambda_D)^2} - \frac{3}{2}\right).\]

Electrons resonant with the wave’s phase velocity (\(v = \omega/k\)) surf it, extracting energy from the field on net because the Maxwellian has slightly more slower particles being accelerated than faster particles being decelerated – a purely collisionless damping mechanism with no dissipation at the particle level, reproduced here by pic_simulate() without any explicit damping term in the equations of motion.

Parameters:

k (float) – Wavenumber, in units of inverse Debye length times \(k\lambda_D\) conventions – concretely, pass the physical wavenumber and set v_th/omega_pe consistently.
v_th (float) – Electron thermal speed.
omega_pe (float) – Electron plasma frequency (1.0 in the normalized units used throughout this module).

Return type:

Returns:

float – Damping rate \(\gamma\) (negative, since the wave decays). The formula is only accurate for weak damping, \(k\lambda_D \lesssim 0.5\).

See also

pic_simulate: Numerically reproduces this decay from first principles.

Examples

The standard textbook benchmark case, \(k\lambda_D = 0.5\):

>>> round(float(landau_damping_rate(k=0.5, v_th=1.0)), 4)
-0.1514

physicskit.plasma.kinetic.langmuir_wave_ic(n_particles, L, k_mode, alpha, v_th, seed=0)[source]#

Small-amplitude single-mode initial condition for a Langmuir (electron plasma) wave.

Uses the same quiet-start displacement as landau_damping_ic(), \(\delta x=(\alpha/k)\sin(kx_0)\) (density perturbation \(n_0(1-\alpha\cos(kx))\)), and additionally gives each particle the coherent velocity \(\delta v = (\alpha\,\omega_{pe}/k)\sin(kx_0)\) (\(\omega_{pe}=1\) in these normalized units). In cold-fluid linear theory the displacement then evolves as \(\xi(x_0,t) = (\alpha/k)\sin(kx_0)\,[\cos\omega_{pe}t + \sin\omega_{pe}t]\): a standing wave ringing in place at \(\omega_{pe}\), with \(\sqrt2\) the amplitude (twice the field energy) of the displacement-only start and a \(\pi/4\) phase shift.

How long it rings is set by Landau damping, i.e. by \(k\lambda_D\) (see landau_damping_rate()), not by the velocity kick: for \(k_{mode}\,v_{th} \ll \omega_{pe}\) this and landau_damping_ic() both oscillate essentially undamped, while for \(k\lambda_D \gtrsim 0.3\) both damp.

Parameters:

n_particles (int) – Number of particles.
L (float) – Domain length (periodic); \(L=2\pi/k_{mode}\) fits one wavelength.
k_mode (float) – Wavenumber of the seeded standing wave.
alpha (float) – Perturbation amplitude (\(\alpha \ll 1\)).
v_th (float) – Thermal speed of the background Maxwellian; keep \(k_{mode}\,v_{th} \ll \omega_{pe}=1\) for the wave to be only weakly Landau-damped.
seed (int) – Random seed for the thermal velocity sampling.

Return type:

Returns:

x, v (ndarray, shape (n_particles,)) – Particle positions and velocities.

See also

landau_damping_ic: The companion density-only perturbation (no coherent velocity kick); it damps or rings for the same \(k\lambda_D\) as this one.
pic_simulate: Evolve this initial condition forward in time.

Examples

>>> import numpy as np
>>> k = 2 * np.pi / 4.0
>>> x, v = langmuir_wave_ic(4, L=4.0, k_mode=k, alpha=0.05, v_th=0.05, seed=0)
>>> bool(np.all((x >= 0) & (x < 4.0)))
True
>>> x.shape, v.shape
((4,), (4,))

physicskit.plasma.kinetic.maxwellian_velocities(n_particles, v_th, seed=0)[source]#

Sample particle velocities from a Maxwellian (Gaussian) distribution.

Parameters:

n_particles (int) – Number of particles to sample.
v_th (float) – Thermal speed (standard deviation of the Gaussian), in normalized velocity units.
seed (int) – Seed for the pseudo-random number generator, for reproducibility.

Return type:

Returns:

ndarray, shape (n_particles,) – Sampled velocities.

Examples

>>> import numpy as np
>>> v = maxwellian_velocities(4, v_th=1.0, seed=0)
>>> bool(np.allclose(v, [0.12573022, -0.13210486, 0.64042265, 0.10490012]))
True

physicskit.plasma.kinetic.pic_simulate(x0, v0, L, ng, dt, steps, qm=-1.0, n0=1.0)[source]#

Run an electrostatic PIC simulation forward in time, recording the field-energy history.

Correctly initializes the leapfrog velocity offset (staggering v0 back by half a step using the field at x0) before repeatedly applying pic_step().

Parameters:

x0 (ndarray) – Initial particle positions, e.g. from landau_damping_ic() or two_stream_ic().
v0 (ndarray) – Initial particle velocities.
L (float) – Domain length (periodic).
ng (int) – Number of grid points.
dt (float) – Time step.
steps (int) – Number of steps to advance.
qm (float) – Charge-to-mass ratio in normalized units.
n0 (float) – Equilibrium number density.

Return type:

Returns:

dict – {"t": ndarray of shape (steps,), "field_energy": ndarray of shape (steps,), "x": final positions, "v": final velocities}.

See also

pic_step: The single-step update repeated here.
landau_damping_rate: Analytic decay rate to compare the field-energy history against.

Examples

>>> import numpy as np
>>> x0, v0 = landau_damping_ic(2000, L=4 * np.pi, k_mode=0.5, alpha=0.01, v_th=1.0, seed=0)
>>> result = pic_simulate(x0, v0, L=4 * np.pi, ng=32, dt=0.1, steps=20)
>>> result["field_energy"].shape
(20,)
>>> bool(np.all(np.isfinite(result["field_energy"])))
True

physicskit.plasma.kinetic.pic_step(x, v, L, ng, dt, qm=-1.0, n0=1.0)[source]#

Advance the electrostatic PIC system one leapfrog step.

Deposits the electron density, solves for the self-consistent field against a uniform neutralizing ion background, gathers the field back onto the particles, and kicks/drifts them – one full cycle of the deposit-solve-gather-push loop at the heart of every PIC code. Velocities are staggered a half step behind positions (standard leapfrog); see pic_simulate() for a driver that initializes that offset correctly.

Parameters:

x (ndarray) – Particle positions in \([0, L)\).
v (ndarray) – Particle velocities, staggered a half step behind x.
L (float) – Domain length (periodic).
ng (int) – Number of grid points.
dt (float) – Time step.
qm (float) – Charge-to-mass ratio in normalized units (-1.0 for electrons with a fixed, uniform ion background of density n0).
n0 (float) – Equilibrium number density.

Return type:

Returns:

x_new, v_new (ndarray) – Updated positions and velocities.
field_energy (float) – \(\int E^2/2\,dx\), evaluated at the field used for this step’s kick.

See also

pic_simulate: Repeated application of this step with correct leapfrog initialization.

Examples

>>> import numpy as np
>>> x = np.linspace(0, 10, 50, endpoint=False)
>>> v = np.zeros(50)
>>> x_new, v_new, fe = pic_step(x, v, L=10.0, ng=32, dt=0.1)
>>> x_new.shape, v_new.shape
((50,), (50,))
>>> bool(np.isfinite(fe))
True

physicskit.plasma.kinetic.solve_poisson_1d(rho, L)[source]#

Solve the 1D periodic Poisson equation \(dE/dx = \rho\) (normalized \(\varepsilon_0=1\)) via FFT.

Parameters:

rho (ndarray) – Net charge density on the grid (e.g. ion background minus deposited electron density). Its mean is discarded, since a uniform charge density has no periodic solution and physically should integrate to zero net charge in the box.
L (float) – Domain length (periodic).

Return type:

Returns:

ndarray, shape (ng,) – Electric field on the grid.

See also

deposit_number_density: Supplies the density this solves for.
interpolate_field: Gathers this field back onto the particles.

Examples

>>> import numpy as np
>>> ng = 64
>>> L = 2 * np.pi
>>> x_grid = np.linspace(0, L, ng, endpoint=False)
>>> rho = np.sin(x_grid)
>>> E = solve_poisson_1d(rho, L)
>>> E_exact = -np.cos(x_grid)
>>> bool(np.max(np.abs(E - E_exact)) < 1e-10)
True

physicskit.plasma.kinetic.two_stream_ic(n_particles, L, v_drift, v_th, seed=0)[source]#

Initial condition for the two-stream instability: two counter-streaming Maxwellian beams.

Splits the particles into two equal populations drifting at \(\pm v_{drift}\), each with thermal spread \(v_{th}\), and seeds the fastest-growing long-wavelength mode with a small density ripple. When \(v_{drift}\) exceeds the thermal spread by enough to make the combined velocity distribution doubly-peaked, the positive-slope region between the two peaks violates the (kinetic) Penrose stability criterion and the ripple grows exponentially, eventually rolling the two beams up into a single phase-space vortex.

Parameters:

n_particles (int) – Number of particles (split evenly between the two beams).
L (float) – Domain length (periodic).
v_drift (float) – Drift speed of each beam (beams move at \(+v_{drift}\) and \(-v_{drift}\)).
v_th (float) – Thermal spread of each beam.
seed (int) – Random seed.

Return type: