physicskit.particle#

physicskit.particle hub · History · Examples · API reference

physicskit.particle: relativistic kinematics, decays, scattering, and nuclear physics.

Works throughout in natural units with \(c=1\) (nuclear specifically uses MeV, the standard convention for its formulas).

  • physicskit.particle.kinematics – four-vectors, Lorentz boosts, invariant mass, rapidity, and center-of-momentum frames.

  • physicskit.particle.decays – two-body decay kinematics, radioactive decay chains (the Bateman equations), and the muon’s three-body Michel-spectrum decay.

  • physicskit.particle.scattering – Mandelstam kinematics and Rutherford scattering.

  • physicskit.particle.nuclear – the semi-empirical mass formula and nuclear reaction Q-values.

  • physicskit.particle.collider – toy branching cascades (particle showers, parton showers) and a schematic detector geometry.

  • physicskit.particle.confinement – toy quark confinement and QCD string breaking.

  • physicskit.particle.electroweak – toy classical symmetry breaking (Higgs mechanism), leading-order QED annihilation, and neutral-meson CP violation.

  • physicskit.particle.neutrinos – two-flavor vacuum neutrino oscillations.

  • physicskit.particle.visualizers – decay-chain and differential-cross-section plots, plus animations of all of the above.

class physicskit.particle.FourVector(E, px, py, pz)[source]#

Bases: object

An energy-momentum four-vector \(p^\mu=(E,p_x,p_y,p_z)\).

Parameters:

  • E (float) – Energy.

  • px (float) – Momentum components.

  • py (float) – Momentum components.

  • pz (float) – Momentum components.

Examples

>>> p = FourVector(2.0, 0.0, 0.0, 1.0)
>>> round(p.mass, 6)
1.732051
property beta#

Speed, \(\beta=|\vec p|/E\).

property gamma#

Lorentz factor, \(\gamma=E/m\).

Raises:

ZeroDivisionError – If the four-vector is massless (no rest frame).

property mass#

Invariant mass, \(m=\sqrt{E^2-|\vec p|^2}\).

The radicand is clipped at 0 before the square root, so a four-vector that is numerically spacelike only by floating-point noise (e.g. the sum of several on-shell four-vectors) returns 0.0 rather than raising on a tiny negative argument.

property p_mag#

Momentum magnitude, \(|\vec p|\).

property p_vec#

the momentum 3-vector.

Type:

ndarray of shape (3,)

class physicskit.particle.ShowerParticle(four_vector, generation, flavor='', parent=None, children=<factory>)[source]#

Bases: object

One node of a branching-cascade tree.

Variables:

  • four_vector (FourVector) – This particle’s lab-frame four-momentum.

  • generation (int) – Branching depth (0 for the initial/primary particle).

  • flavor (str, default="") – Parton flavor label ("q", "qbar", "g"); unused (left empty) by simple_shower().

  • parent (ShowerParticle, optional) – The particle this one was produced from (None for the root).

  • children (list of ShowerParticle) – The two daughters, if this particle branched further (empty for a final-state / leaf particle).

Parameters:
children: list#
flavor: str = ''#
four_vector: FourVector#
generation: int#
property is_leaf#

Whether this particle is final-state (did not branch further).

parent: ShowerParticle | None = None#
physicskit.particle.activity(N, decay_constant_)[source]#

Decay rate (activity), \(A=\lambda N\).

Examples

>>> activity(100.0, 0.01)
1.0
physicskit.particle.animate_cp_asymmetry(t, delta_m, gamma_s, gamma_l, epsilon, interval=60, ax=None)[source]#

Animate the neutral-meson decay-rate CP asymmetry \(A(t)\) building up, from physicskit.particle.electroweak.cp_asymmetry().

Parameters:

  • t (ndarray) – Times.

  • delta_m – See physicskit.particle.electroweak.meson_decay_rates_cp_eigenstate().

  • gamma_s – See physicskit.particle.electroweak.meson_decay_rates_cp_eigenstate().

  • gamma_l – See physicskit.particle.electroweak.meson_decay_rates_cp_eigenstate().

  • epsilon – See physicskit.particle.electroweak.meson_decay_rates_cp_eigenstate().

  • interval (int, default=60) – Delay between frames, in milliseconds.

  • ax (matplotlib.axes.Axes, optional) – Axes to draw into; a new figure is created if omitted.

Returns:

matplotlib.animation.FuncAnimation

Examples

>>> import numpy as np
>>> t = np.linspace(0, 20, 100)
>>> anim = animate_cp_asymmetry(t, delta_m=0.5, gamma_s=1.0, gamma_l=0.1, epsilon=0.002 + 0.001j)
>>> anim.__class__.__name__
'FuncAnimation'
physicskit.particle.animate_decay_chain_bars(t, populations, labels=None, interval=150, ax=None)[source]#

Animate each species’ population evolving over time in a decay chain.

Parameters:

  • t (ndarray of shape (n_t,)) – Times.

  • populations (ndarray of shape (n_species, n_t)) – Population of each species over time, e.g. from physicskit.particle.decays.bateman_decay_chain().

  • labels (sequence of str, optional) – Bar labels; defaults to "species 1", "species 2", …

  • interval (int, default=150) – Delay between frames, in milliseconds.

  • ax (matplotlib.axes.Axes, optional) – Axes to draw into; a new figure is created if omitted.

Returns:

matplotlib.animation.FuncAnimation

Examples

>>> import numpy as np
>>> from physicskit.particle.decays import bateman_decay_chain
>>> t = np.linspace(0, 10, 30)
>>> N = bateman_decay_chain(1000.0, [0.5, 0.2], t)
>>> anim = animate_decay_chain_bars(t, N)
>>> anim.__class__.__name__
'FuncAnimation'
physicskit.particle.animate_detector_event(four_vectors, charges, B, detector_radii=(3.0, 6.0), n_frames=60, interval=100, ax=None)[source]#

Animate final-state tracks growing outward from the vertex and landing in a schematic calorimeter.

Charged particles curve (radius \(r=p_T/(qB)\), see physicskit.particle.collider.charged_track_points()); neutral particles fly in straight lines. Each frame reveals one more step along every track; a marker is drawn at each track’s current endpoint, standing in for a calorimeter energy deposit once the track reaches the outer radius.

Parameters:

  • four_vectors (sequence of physicskit.particle.kinematics.FourVector) – Final-state particles.

  • charges (sequence of float) – Charge of each particle, in units of the elementary charge.

  • B (float) – Magnetic field strength (see charged_track_points()).

  • detector_radii ((float, float), default=(3.0, 6.0)) – Radii of the schematic tracker and calorimeter circles.

  • n_frames (int, default=60) – Number of animation frames (and points sampled along each track).

  • interval (int, default=100) – Delay between frames, in milliseconds.

  • ax (matplotlib.axes.Axes, optional) – Axes to draw into; a new figure is created if omitted.

Returns:

matplotlib.animation.FuncAnimation

Examples

>>> from physicskit.particle.kinematics import FourVector
>>> ps = [FourVector(5.0, 3.0, 0.0, 3.0), FourVector(4.0, 0.0, 3.0, 2.0)]
>>> anim = animate_detector_event(ps, charges=[1.0, -1.0], B=1.0)
>>> anim.__class__.__name__
'FuncAnimation'
physicskit.particle.animate_higgs_rollover(phi_t, a, b, interval=30, ax=None)[source]#

Animate a scalar field rolling from the unstable symmetric point into a broken-symmetry vacuum, from physicskit.particle.electroweak.higgs_field_rollover().

The static curve is physicskit.particle.electroweak.higgs_potential(); a point (with a fading trail) traces \((\phi(t), V(\phi(t)))\) as the field rolls down one side of the double well.

Parameters:

  • phi_t (ndarray) – Field trajectory \(\phi(t)\), e.g. the first output of higgs_field_rollover().

  • a (float) – Potential parameters, see higgs_potential().

  • b (float) – Potential parameters, see higgs_potential().

  • interval (int, default=30) – Delay between frames, in milliseconds.

  • ax (matplotlib.axes.Axes, optional) – Axes to draw into; a new figure is created if omitted.

Returns:

matplotlib.animation.FuncAnimation

Examples

>>> import numpy as np
>>> from physicskit.particle.electroweak import higgs_field_rollover
>>> t = np.linspace(0, 40, 100)
>>> phi, _ = higgs_field_rollover(1e-3, 0.0, a=1.0, b=1.0, t_eval=t, damping=0.05)
>>> anim = animate_higgs_rollover(phi, a=1.0, b=1.0)
>>> anim.__class__.__name__
'FuncAnimation'
physicskit.particle.animate_michel_histogram(x_samples, batch_size=20, n_bins=30, interval=100, ax=None)[source]#

Animate a histogram of sampled Michel-spectrum electron energies building up event-by-event, converging to physicskit.particle.decays.michel_spectrum().

Parameters:

  • x_samples (ndarray) – Sampled scaled electron energies, e.g. from physicskit.particle.decays.sample_michel_electron_energies().

  • batch_size (int, default=20) – Number of additional samples revealed per frame.

  • n_bins (int, default=30) – Number of histogram bins.

  • interval (int, default=100) – Delay between frames, in milliseconds.

  • ax (matplotlib.axes.Axes, optional) – Axes to draw into; a new figure is created if omitted.

Returns:

matplotlib.animation.FuncAnimation

Examples

>>> import numpy as np
>>> from physicskit.particle.decays import sample_michel_electron_energies
>>> xs = sample_michel_electron_energies(500, rng=np.random.default_rng(0))
>>> anim = animate_michel_histogram(xs, batch_size=100)
>>> anim.__class__.__name__
'FuncAnimation'
physicskit.particle.animate_neutrino_oscillation(L, E, theta, delta_m2, interval=60, ax=None)[source]#

Animate \(P(\nu_e)\) and \(P(\nu_\mu)\) growing as a function of distance traveled, from physicskit.particle.neutrinos.oscillation_probability().

Parameters:

  • L (ndarray) – Baselines (distance traveled), in km, in increasing order.

  • E (float) – Neutrino energy, in GeV.

  • theta (float) – Mixing angle, radians.

  • delta_m2 (float) – Mass-squared splitting, in \({\rm eV}^2\).

  • interval (int, default=60) – Delay between frames, in milliseconds.

  • ax (matplotlib.axes.Axes, optional) – Axes to draw into; a new figure is created if omitted.

Returns:

matplotlib.animation.FuncAnimation

Examples

>>> import numpy as np
>>> L = np.linspace(0, 1000, 100)
>>> anim = animate_neutrino_oscillation(L, E=1.0, theta=0.6, delta_m2=2.5e-3)
>>> anim.__class__.__name__
'FuncAnimation'
physicskit.particle.animate_particle_cascade(root, detector_radii=(3.0, 6.0), track_length=0.6, interval=150, ax=None)[source]#

Animate a toy shower cascade (physicskit.particle.collider.simple_shower()) propagating outward from the collision vertex through a schematic detector.

Concentric circles mark a “tracker” (inner) and “calorimeter” (outer) region; each branch of the shower tree appears as one more track segment per frame, in order of production (generation).

Parameters:

  • root (physicskit.particle.collider.ShowerParticle) – Root of the shower tree, e.g. from physicskit.particle.collider.simple_shower().

  • detector_radii ((float, float), default=(3.0, 6.0)) – Radii of the schematic tracker and calorimeter circles.

  • track_length (float, default=0.6) – Schematic (not physically derived) drawn length of each track segment – purely a layout choice, not a decay-length calculation.

  • interval (int, default=150) – Delay between frames, in milliseconds.

  • ax (matplotlib.axes.Axes, optional) – Axes to draw into; a new figure is created if omitted.

Returns:

matplotlib.animation.FuncAnimation

Examples

>>> import numpy as np
>>> from physicskit.particle.collider import simple_shower
>>> root = simple_shower(50.0, E_threshold=5.0, rng=np.random.default_rng(0))
>>> anim = animate_particle_cascade(root)
>>> anim.__class__.__name__
'FuncAnimation'
physicskit.particle.animate_parton_shower(root, detector_radii=(3.0, 6.0), track_length=0.6, n_jets=2, interval=150, ax=None)[source]#

Animate a toy parton shower (physicskit.particle.collider.parton_shower()) developing into jets.

Reuses the tree-reveal machinery of animate_particle_cascade(), color-coding branches by parton flavor (quark/antiquark vs. gluon) and overlaying the final jet-axis directions found by physicskit.particle.collider.cluster_into_jets().

Parameters:
Returns:

matplotlib.animation.FuncAnimation

Examples

>>> import numpy as np
>>> from physicskit.particle.collider import parton_shower
>>> root = parton_shower(50.0, E_threshold=5.0, rng=np.random.default_rng(0))
>>> anim = animate_parton_shower(root)
>>> anim.__class__.__name__
'FuncAnimation'
physicskit.particle.animate_qed_angular_distribution(sqrt_s_values, alpha=0.0072973525692838015, n_theta=200, interval=200, ax=None)[source]#

Animate the \(e^+e^-\to\mu^+\mu^-\) angular distribution’s shape as \(\sqrt s\) is swept, from physicskit.particle.electroweak.qed_dsigma_domega_mumu().

Parameters:

  • sqrt_s_values (ndarray) – Center-of-mass energies to sweep through (one per frame).

  • alpha (float, default=ALPHA_FS) – Fine-structure constant.

  • n_theta (int, default=200) – Number of polar-angle samples per frame.

  • interval (int, default=200) – Delay between frames, in milliseconds.

  • ax (matplotlib.axes.Axes, optional) – Polar axes to draw into; new polar axes are created if omitted.

Returns:

matplotlib.animation.FuncAnimation

Examples

>>> import numpy as np
>>> sqrt_s = np.linspace(5.0, 20.0, 20)
>>> anim = animate_qed_angular_distribution(sqrt_s)
>>> anim.__class__.__name__
'FuncAnimation'
physicskit.particle.animate_string_breaking(sim, interval=120, ax=None)[source]#

Animate a receding quark pair’s confining string stretching and breaking, from physicskit.particle.confinement.string_break_chain().

A thick line (whose width and color track the stored tension, energy_per_segment) stretches between the two receding quark markers; each time the simulation records a break, it snaps into one more, shorter segment, with a new marker appearing at each newly created quark-antiquark endpoint. (The endpoints of the n_segments(t) current strings are drawn evenly spaced across the current total extent – a schematic simplification, see string_break_chain()’s docstring, rather than tracking each new pair’s exact, individually frozen creation position.)

Parameters:

  • sim (dict) – Output of physicskit.particle.confinement.string_break_chain().

  • interval (int, default=120) – Delay between frames, in milliseconds.

  • ax (matplotlib.axes.Axes, optional) – Axes to draw into; a new figure is created if omitted.

Returns:

matplotlib.animation.FuncAnimation

Examples

>>> import numpy as np
>>> from physicskit.particle.confinement import string_break_chain
>>> sim = string_break_chain(np.linspace(0, 10, 40), v=0.3, kappa=1.0, m_q=1.0, n_breaks=2)
>>> anim = animate_string_breaking(sim)
>>> anim.__class__.__name__
'FuncAnimation'
physicskit.particle.bateman_decay_chain(N0, decay_constants, t)[source]#

Population of every species in a linear radioactive decay chain.

Solves the Bateman equations for a chain \(1\to2\to\cdots\to n\) with decay constants decay_constants (species 1 starts with population N0, all daughters start at zero):

\[N_k(t) = N_0\Big(\prod_{i=1}^{k-1}\lambda_i\Big) \sum_{i=1}^{k}\frac{e^{-\lambda_i t}} {\prod_{j=1,j\neq i}^{k}(\lambda_j-\lambda_i)}\]

with \(N_1(t)=N_0e^{-\lambda_1 t}\). A zero decay constant for the final species represents a stable end product.

Parameters:

  • N0 (float) – Initial population of species 1.

  • decay_constants (sequence of float) – Decay constant of each of the n species in the chain. Must be pairwise distinct among the species actually populated (the closed-form sum divides by lambda_j - lambda_i).

  • t (array_like) – Times at which to evaluate the populations.

Returns:

ndarray of shape (n, len(t)) – Population of each species over time.

Examples

>>> t = np.array([0.0, 1.0])
>>> N = bateman_decay_chain(100.0, [1.0, 2.0], t)
>>> N.shape
(2, 2)
>>> round(float(N[0, 0]), 6)
100.0
physicskit.particle.binding_energy_per_nucleon(Z, A)[source]#

Binding energy per nucleon, semf_binding_energy(Z, A) / A, in MeV.

Examples

>>> bpn = binding_energy_per_nucleon(26, 56)
>>> 8.0 < bpn < 9.0
True
physicskit.particle.boost(four_vector, beta, axis='z')[source]#

Lorentz-boost a FourVector along a Cartesian axis.

Parameters:

  • four_vector (FourVector) – The four-vector to boost.

  • beta (float) – Boost velocity, \(-1 < \beta < 1\), in units of \(c=1\).

  • axis ({"x", "y", "z"}, default="z") – Cartesian axis of the boost.

Returns:

FourVector – The boosted four-vector.

Examples

>>> p = FourVector(1.0, 0.0, 0.0, 0.0)  # a particle at rest, mass 1
>>> pb = boost(p, 0.6, axis="z")
>>> round(pb.E, 6), round(pb.pz, 6)
(1.25, 0.75)
physicskit.particle.boost_generic(four_vector, beta_vec)[source]#

Lorentz-boost a FourVector along an arbitrary 3-velocity.

Generalizes boost() to a boost direction that need not lie along a Cartesian axis – the transformation used, e.g., to take a daughter four-momentum computed in a parent’s rest frame back to the lab frame when the parent’s own flight direction is not one of the axes (as happens at every vertex of a branching cascade; see physicskit.particle.collider).

\[E' = \gamma(E + \vec\beta\cdot\vec p), \qquad \vec p' = \vec p + \left[\frac{\gamma-1}{\beta^2}(\vec\beta\cdot\vec p) + \gamma E\right]\vec\beta\]

which reduces to boost() when \(\vec\beta\) is aligned with a single axis.

Parameters:

  • four_vector (FourVector) – The four-vector to boost.

  • beta_vec (array_like of shape (3,)) – Boost velocity, \(|\vec\beta| < 1\), in units of \(c=1\).

Returns:

FourVector – The boosted four-vector.

Examples

>>> p = FourVector(1.0, 0.0, 0.0, 0.0)  # a particle at rest, mass 1
>>> pb = boost_generic(p, [0.0, 0.0, 0.6])
>>> round(pb.E, 6), round(pb.pz, 6)
(1.25, 0.75)
physicskit.particle.boost_to_com(four_vectors)[source]#

The 3-velocity of the center-of-momentum frame of a system.

Parameters:

four_vectors (iterable of FourVector) – The particles making up the system.

Returns:

ndarray of shape (3,) – The velocity (units of \(c=1\)) of the system’s total momentum relative to the current frame, \(\vec\beta_{\rm com} = \vec p_{\rm tot}/E_{\rm tot}\).

Examples

>>> p1 = FourVector(2.0, 0.0, 0.0, 1.0)
>>> p2 = FourVector(1.0, 0.0, 0.0, 0.0)
>>> beta_com = boost_to_com([p1, p2])
>>> [round(float(b), 6) for b in beta_com]
[0.0, 0.0, 0.333333]
physicskit.particle.charged_track_points(four_vector, charge, B, vertex=(0.0, 0.0), n_points=100, path_length=1.0)[source]#

2D transverse-plane trajectory of a particle in a uniform axial magnetic field.

A charged particle with transverse momentum \(p_T\) follows a circular arc of radius

\[r = \frac{p_T}{qB}\]

the standard relation for the curvature of a charged track in a uniform axial field (see e.g. the Particle Data Group’s “Passage of particles through matter” review); here in schematic natural units where \(q\) is in units of the elementary charge and \(B\) is scaled so that \(r\) comes out directly in the detector’s length unit (the familiar practical-unit form is \(r[{\rm m}] =p_T[{\rm GeV}/c]/(0.3\,B[{\rm T}]\,q[e])\)). A neutral particle (charge=0) instead follows a straight line.

Parameters:

  • four_vector (FourVector) – The particle’s four-momentum.

  • charge (float) – Charge, in units of the elementary charge (0 for neutral).

  • B (float) – Magnetic field along \(+z\), out of the transverse plane. The Lorentz force \(q\,\mathbf{v}\times\mathbf{B}\) then bends positive charges clockwise (for B > 0) and negative charges counterclockwise, as seen looking down the \(z\) axis.

  • vertex (array_like of shape (2,), default=(0, 0)) – Starting point of the track.

  • n_points (int, default=100) – Number of points to sample along the trajectory.

  • path_length (float, default=1.0) – For a neutral track, the length of the straight segment drawn; for a charged track, the arc length traced out (capped implicitly by how many points are requested along it).

Returns:

ndarray of shape (n_points, 2) – Points along the trajectory, starting at vertex and moving off in four_vector’s initial transverse direction.

Examples

>>> from physicskit.particle.kinematics import FourVector
>>> p = FourVector(10.0, 3.0, 0.0, 5.0)
>>> pts = charged_track_points(p, charge=0.0, B=1.0, path_length=2.0)
>>> np.allclose(pts[0], [0.0, 0.0]) and np.allclose(pts[-1], [2.0, 0.0])
True
physicskit.particle.cluster_into_jets(leaves, n_jets=2, n_iter=25, seed=0)[source]#

Group final-state shower particles into jets by momentum direction.

A minimal, dependency-free k-means clustering on the unit momentum directions (i.e. by angular proximity) – a schematic stand-in for a real sequential jet algorithm (\(k_T\), anti-\(k_T\), …), adequate for visually grouping a toy shower’s final state into a couple of back-to-back jet cones.

Parameters:

  • leaves (sequence of ShowerParticle) – Final-state particles, e.g. from shower_leaves().

  • n_jets (int, default=2) – Number of jets to form (capped at len(leaves)).

  • n_iter (int, default=25) – Number of k-means (Lloyd) iterations.

  • seed (int, default=0) – Seed for the deterministic random initialization.

Returns:

list of list of ShowerParticle – One list of member particles per jet.

Examples

>>> import numpy as np
>>> root = parton_shower(100.0, E_threshold=5.0, rng=np.random.default_rng(0))
>>> jets = cluster_into_jets(shower_leaves(root), n_jets=2)
>>> len(jets) == 2
True
physicskit.particle.cp_asymmetry(t, delta_m, gamma_s, gamma_l, epsilon)[source]#

Decay-rate CP asymmetry \(A(t)\), from meson_decay_rates_cp_eigenstate().

\[A(t) = \frac{\Gamma(\bar P^0(t)\to f) - \Gamma(P^0(t)\to f)} {\Gamma(\bar P^0(t)\to f) + \Gamma(P^0(t)\to f)}\]

which oscillates at frequency \(\Delta m\) inside a decaying envelope, vanishing identically for \(\varepsilon=0\) (no CP violation).

Parameters:
Returns:

ndarray or float

Examples

>>> import numpy as np
>>> t = np.linspace(0, 5, 6)
>>> A = cp_asymmetry(t, delta_m=1.0, gamma_s=5.0, gamma_l=0.5, epsilon=0.0)
>>> np.allclose(A, 0.0)
True
physicskit.particle.decay_constant(half_life_)[source]#

Decay constant from half-life, \(\lambda=\ln2/t_{1/2}\).

Examples

>>> round(decay_constant(np.log(2)), 6)
1.0
physicskit.particle.flatten_shower(root)[source]#

All nodes of a shower tree (root included), in no particular order.

Parameters:

root (ShowerParticle)

Returns:

list of ShowerParticle

physicskit.particle.half_life(decay_constant_)[source]#

Half-life from decay constant, \(t_{1/2}=\ln2/\lambda\).

Examples

>>> round(half_life(1.0), 6)
0.693147
physicskit.particle.higgs_field_rollover(phi0, phidot0, a, b, t_eval, damping=0.0)[source]#

Classically evolve a homogeneous scalar field rolling in higgs_potential().

Integrates the (0+1 dimensional – no spatial gradient/kink structure, a single point standing in for a spatially uniform field configuration) classical field equation

\[\ddot\phi = -\frac{dV}{d\phi} - \gamma\dot\phi = 2a\phi - 4b\phi^3 - \gamma\dot\phi\]

Started near the unstable symmetric point \(\phi\approx0\) with a tiny perturbation, the field is repelled from it and rolls down into one of the two true vacua \(\phi=\pm v\) (higgs_vev()) – whichever side the initial perturbation pushes it toward, a toy realization of spontaneous symmetry breaking. A small damping \(\gamma>0\) (representing energy radiated into other field modes, not modeled explicitly) lets the field settle into the chosen vacuum instead of oscillating in it forever.

Parameters:

  • phi0 (float) – Initial field value and velocity.

  • phidot0 (float) – Initial field value and velocity.

  • a (float) – Potential parameters, see higgs_potential().

  • b (float) – Potential parameters, see higgs_potential().

  • t_eval (array_like) – Times at which to report the solution.

  • damping (float, default=0.0) – Damping coefficient \(\gamma\ge0\).

Returns:

phi, phidot (ndarray) – Field value and velocity at each time in t_eval.

Examples

>>> import numpy as np
>>> t = np.linspace(0, 60, 600)
>>> phi, phidot = higgs_field_rollover(1e-3, 0.0, a=1.0, b=1.0, t_eval=t, damping=0.08)
>>> v = higgs_vev(1.0, 1.0)
>>> bool(abs(abs(phi[-1]) - v) < 0.05)
True
physicskit.particle.higgs_potential(phi, a, b)[source]#

The classical \(\phi^4\) double-well potential, \(V(\phi)=-a\phi^2+b\phi^4\).

For \(a,b>0\) this has an unstable extremum at \(\phi=0\) (the symmetric, “false vacuum” point) and two degenerate true minima at \(\phi=\pm v\), \(v=\sqrt{a/(2b)}\) – the field must “choose” one of them, spontaneously breaking the \(\phi\to-\phi\) symmetry. This is the same quartic potential that gives the Higgs field its nonzero vacuum expectation value.

Parameters:

  • phi (array_like) – Field value.

  • a (float) – Potential parameters, \(a,b>0\).

  • b (float) – Potential parameters, \(a,b>0\).

Returns:

ndarray or float

Examples

>>> float(higgs_potential(0.0, a=1.0, b=1.0))
0.0
physicskit.particle.higgs_vev(a, b)[source]#

The true-vacuum field value, \(v=\sqrt{a/(2b)}\), minimizing higgs_potential().

Examples

>>> round(higgs_vev(a=2.0, b=1.0), 6)
1.0
physicskit.particle.impact_parameter(theta, Z1, Z2, E_kin, alpha=0.0072973525692838015)[source]#

Classical Rutherford impact parameter producing scattering angle theta.

\[b(\theta) = \frac{Z_1Z_2\alpha}{2E_{\rm kin}}\cot(\theta/2)\]
Parameters:

  • theta (array_like) – Scattering angle (radians).

  • Z1 (float) – Projectile and target charge numbers.

  • Z2 (float) – Projectile and target charge numbers.

  • E_kin (float) – Projectile kinetic energy.

  • alpha (float, default=ALPHA_FS) – Coupling constant.

Returns:

ndarray or float

Examples

>>> round(float(impact_parameter(np.pi, Z1=1, Z2=79, E_kin=5.0)), 6)
0.0
physicskit.particle.invariant_mass(four_vectors)[source]#

The invariant mass of a system of four-vectors.

Parameters:

four_vectors (iterable of FourVector) – The particles making up the system.

Returns:

float – The mass of the sum of the four-vectors – e.g. the reconstructed mass of a resonance from its decay products.

Examples

>>> p1 = FourVector(1.0, 0.0, 0.0, 0.6)
>>> p2 = FourVector(1.0, 0.0, 0.0, -0.6)
>>> round(invariant_mass([p1, p2]), 6)
2.0
physicskit.particle.mandelstam_s(p1, p2)[source]#

Mandelstam \(s=(p_1+p_2)^2\), the total invariant mass squared.

Parameters:
Returns:

float

Examples

>>> from physicskit.particle.kinematics import FourVector
>>> p1 = FourVector(1.0, 0.0, 0.0, 0.6)
>>> p2 = FourVector(1.0, 0.0, 0.0, -0.6)
>>> round(mandelstam_s(p1, p2), 6)
4.0
physicskit.particle.mandelstam_t(p1, p3)[source]#

Mandelstam \(t=(p_1-p_3)^2\), the momentum-transfer invariant.

Parameters:
Returns:

float

physicskit.particle.mandelstam_u(p1, p4)[source]#

Mandelstam \(u=(p_1-p_4)^2\), the other momentum-transfer invariant.

Parameters:
Returns:

float

Notes

For a 2-to-2 process \(1+2\to3+4\), \(s+t+u=\sum_i m_i^2\).

physicskit.particle.meson_decay_rates_cp_eigenstate(t, delta_m, gamma_s, gamma_l, epsilon)[source]#

Time-dependent \(P^0,\bar P^0\to f_{CP}\) decay rates, Wigner-Weisskopf toy model.

A neutral meson \(P^0\) (modeled on the \(K^0\)-\(\bar K^0\) system) mixes through the effective non-Hermitian Hamiltonian \(H=M-i\Gamma/2\); its mass eigenstates are the short/long-lived combinations \(P_{S,L}\approx P^0\pm\bar P^0\) (up to the small CP-violating admixture parametrized by the complex parameter \(\varepsilon\), \(|\varepsilon|\ll1\)). Decaying to a common CP eigenstate final state \(f\) (e.g. \(\pi^+\pi^-\)), the two flavor-tagged rates interfere as

\[\begin{split}\Gamma(P^0(t)\to f) &= e^{-\Gamma_S t} + |\varepsilon|^2e^{-\Gamma_L t} + 2|\varepsilon|e^{-\Gamma t}\cos(\Delta m\,t-\phi_\varepsilon) \\ \Gamma(\bar P^0(t)\to f) &= e^{-\Gamma_S t} + |\varepsilon|^2e^{-\Gamma_L t} - 2|\varepsilon|e^{-\Gamma t}\cos(\Delta m\,t-\phi_\varepsilon)\end{split}\]

with \(\Gamma=(\Gamma_S+\Gamma_L)/2\); this is the standard interference pattern used to measure \(\Delta m\) and the \(K_S\) lifetime from \(K^0/\bar K^0\to\pi\pi\) decay-rate oscillations (see the Particle Data Group’s “CP violation in \(K_L\) decays” review, or Griffiths, Introduction to Elementary Particles, Ch. 8).

Parameters:

  • t (array_like) – Proper time since production (tagged as a pure \(P^0\) or \(\bar P^0\) flavor eigenstate at \(t=0\)).

  • delta_m (float) – Mass splitting between the two mass eigenstates, \(\Delta m=m_L-m_S\).

  • gamma_s (float) – Decay widths of the short- and long-lived mass eigenstates.

  • gamma_l (float) – Decay widths of the short- and long-lived mass eigenstates.

  • epsilon (complex) – CP-violation parameter, \(|\varepsilon|\ll1\).

Returns:

gamma_meson, gamma_mesonbar (ndarray)

Examples

>>> import numpy as np
>>> t = np.linspace(0, 5, 6)
>>> g, gbar = meson_decay_rates_cp_eigenstate(t, delta_m=1.0, gamma_s=5.0, gamma_l=0.5, epsilon=0.0)
>>> np.allclose(g, gbar)
True
physicskit.particle.michel_spectrum(x)[source]#

The (unpolarized, massless-electron) Michel spectrum shape for muon decay.

\[\frac{d\Gamma}{dx} = 2x^2(3-2x), \qquad x=\frac{2E_e}{m_\mu}\in[0,1]\]

the standard leading-order result for \(\mu^-\to e^-\bar\nu_\mu\nu_e\) from the V-A four-fermion weak interaction (Michel 1950; see e.g. Commins & Bucksbaum, Weak Interactions of Leptons and Quarks, Ch. 4, or the muon-decay chapter of any standard particle physics text), in the limit \(m_e/m_\mu\to0\) and for an unpolarized muon (no Michel parameter \(\rho,\eta,\xi,\delta\) dependence – those parametrize corrections from polarization and a nonzero electron mass/anomalous couplings, all set to their V-A/Standard-Model, unpolarized values here). Normalized so that \(\int_0^1 (d\Gamma/dx)\,dx = 1\).

Parameters:

x (array_like) – Scaled electron energy, \(x=2E_e/m_\mu\in[0,1]\).

Returns:

ndarray or float

Examples

>>> import numpy as np
>>> x = np.linspace(0, 1, 100001)
>>> round(float(np.trapezoid(michel_spectrum(x), x)), 4)
1.0
physicskit.particle.muon_decay_event(m_mu, rng=None)[source]#

Sample one \(\mu^-\to e^-+\bar\nu_\mu+\nu_e\) decay in the muon rest frame.

Toy-but-exact three-body kinematics built from two chained exact two-body decays (reusing two_body_decay(), so energy-momentum is conserved to machine precision at both vertices):

  1. The electron’s energy is drawn from the Michel spectrum (massless electron approximation) via sample_michel_electron_energies(), fixing the invariant mass of the recoiling neutrino pair through \(m_{\nu\nu}^2=m_\mu^2-2m_\mu E_e\) (exact for \(m_e=0\)). two_body_decay() then gives the electron and the neutrino-pair “system” four-momenta, back to back in the muon rest frame, with an isotropically sampled decay angle (the standard unpolarized-muon assumption).

  2. The neutrino-pair system is split into its two (massless) neutrinos isotropically in its own rest frame via a second call to two_body_decay(), then boosted back to the muon rest frame with boost_generic().

Both neutrino masses are set to zero – an excellent approximation given their sub-eV masses compared to \(m_\mu\sim100\) MeV.

Parameters:

  • m_mu (float) – Muon mass (in whatever energy unit the caller uses consistently).

  • rng (numpy.random.Generator, optional) – Random number generator; a fresh default one is used if omitted.

Returns:

p_e, p_numu_bar, p_nue (FourVector) – The three decay products, in the muon rest frame.

Examples

>>> from physicskit.particle.kinematics import invariant_mass
>>> rng = np.random.default_rng(1)
>>> p_e, p_numu_bar, p_nue = muon_decay_event(105.658, rng=rng)
>>> round(invariant_mass([p_e, p_numu_bar, p_nue]), 6)
105.658
physicskit.particle.oscillation_probability(L, E, theta, delta_m2)[source]#

Two-flavor vacuum oscillation probability \(P(\nu_e\to\nu_\mu)\).

\[P(\nu_e\to\nu_\mu) = \sin^2(2\theta)\, \sin^2\!\left(1.267\,\frac{\Delta m^2 L}{E}\right)\]

the standard PDG form in practical units: \(L\) in km, \(E\) in GeV, \(\Delta m^2\) in \({\rm eV}^2\) (the numerical constant 1.267 absorbs \(\hbar,c\) and the unit conversion; see the Particle Data Group’s “Neutrino Mixing” review).

Parameters:

  • L (array_like) – Baseline (distance traveled), in km.

  • E (array_like) – Neutrino energy, in GeV.

  • theta (float) – Mixing angle, in radians.

  • delta_m2 (float) – Mass-squared splitting, in \({\rm eV}^2\).

Returns:

ndarray or float

Examples

>>> round(float(oscillation_probability(0.0, 1.0, 0.5, 2.5e-3)), 10)
0.0
>>> theta = np.pi / 4  # maximal mixing
>>> L_E = np.pi / 2 / (1.267 * 2.5e-3)  # argument = pi/2 -> sin^2 = 1
>>> round(float(oscillation_probability(L_E, 1.0, theta, 2.5e-3)), 6)
1.0
physicskit.particle.parton_shower(E0, flavor0='q', m0_fraction=0.3, daughter_fraction=0.42, E_threshold=1.0, gluon_splitting_prob=0.2, max_generations=10, rng=None)[source]#

A toy DGLAP-flavored parton shower: quarks radiate gluons, gluons split into gluon pairs or quark-antiquark pairs.

Uses the same exact-two-body-split-and-boost construction as simple_shower() (see the module docstring), with a simplified flavor rule standing in for the DGLAP splitting functions \(P_{qq},P_{gg},P_{qg}\):

  • a quark or antiquark always radiates a gluon (q -> q g, qbar -> qbar g), the only flavor-conserving 1-to-2 splitting available to it;

  • a gluon splits into two gluons with probability 1-gluon_splitting_prob (g -> g g) or into a quark-antiquark pair with probability gluon_splitting_prob (g -> q qbar).

Parameters:

  • E0 (float) – Energy of the initial parton.

  • flavor0 (str, default="q") – Flavor of the initial parton ("q", "qbar", or "g").

  • m0_fraction – See simple_shower().

  • daughter_fraction – See simple_shower().

  • E_threshold – See simple_shower().

  • max_generations – See simple_shower().

  • rng – See simple_shower().

  • gluon_splitting_prob (float, default=0.2) – Probability that a gluon branching produces a quark-antiquark pair rather than two gluons.

Returns:

ShowerParticle – The root of the shower tree.

Examples

>>> import numpy as np
>>> root = parton_shower(100.0, E_threshold=5.0, rng=np.random.default_rng(0))
>>> leaves = shower_leaves(root)
>>> from physicskit.particle.kinematics import invariant_mass
>>> bool(abs(invariant_mass([leaf.four_vector for leaf in leaves]) - root.four_vector.mass) < 1e-6)
True
physicskit.particle.plot_decay_chain(t, populations, labels=None, ax=None)[source]#

Plot each species’ population over time in a radioactive decay chain.

Parameters:

  • t (ndarray of shape (n_t,)) – Times.

  • populations (ndarray of shape (n_species, n_t)) – Population of each species over time, e.g. from physicskit.particle.decays.bateman_decay_chain().

  • labels (sequence of str, optional) – Label for each species; defaults to "species 1", "species 2", …

  • 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.particle.decays import bateman_decay_chain
>>> t = np.linspace(0, 10, 50)
>>> N = bateman_decay_chain(1000.0, [0.5, 0.2], t)
>>> fig, ax = plot_decay_chain(t, N, labels=["Parent", "Daughter"])
>>> isinstance(fig, plt.Figure)
True
physicskit.particle.plot_differential_cross_section(theta, dsigma_domega, ax=None, log_scale=True)[source]#

Plot a differential cross section \(d\sigma/d\Omega\) vs. scattering angle.

Parameters:

  • theta (ndarray of shape (n,)) – Scattering angles (radians).

  • dsigma_domega (ndarray of shape (n,)) – Differential cross section values, e.g. from physicskit.particle.scattering.rutherford_dsigma_domega().

  • ax (matplotlib.axes.Axes, optional) – Axes to draw into; a new figure is created if omitted.

  • log_scale (bool, default=True) – Plot the cross section on a log y-axis, appropriate for the Rutherford formula’s steep small-angle divergence.

Returns:

  • fig (matplotlib.figure.Figure)

  • ax (matplotlib.axes.Axes)

Examples

>>> import numpy as np
>>> from physicskit.particle.scattering import rutherford_dsigma_domega
>>> theta = np.linspace(0.1, np.pi - 0.1, 50)
>>> vals = rutherford_dsigma_domega(theta, Z1=2, Z2=79, E_kin=5.0)
>>> fig, ax = plot_differential_cross_section(theta, vals)
>>> isinstance(fig, plt.Figure)
True
physicskit.particle.q_value(reactant_masses, product_masses)[source]#

Reaction Q-value, \(Q=\sum m_{\rm reactants}-\sum m_{\rm products}\).

A positive Q means the reaction is exothermic (releases energy); a negative Q means it is endothermic and requires that much additional kinetic energy to proceed. Masses (or mass-energies) may be in any consistent unit.

Parameters:

  • reactant_masses (iterable of float) – Masses of the reactants.

  • product_masses (iterable of float) – Masses of the products.

Returns:

float

Examples

>>> round(q_value([2.014102, 3.016049], [4.002602, 1.008665]), 6)  # D + T -> He4 + n
0.018884
physicskit.particle.qed_dsigma_domega_mumu(cos_theta, sqrt_s, alpha=0.0072973525692838015)[source]#

Leading-order QED differential cross section for \(e^+e^-\to\mu^+\mu^-\).

\[\frac{d\sigma}{d\Omega} = \frac{\alpha^2}{4s}\left(1+\cos^2\theta\right)\]

the standard tree-level, single-photon-exchange, unpolarized result in the ultrarelativistic (massless-fermion) limit – e.g. Peskin & Schroeder, An Introduction to Quantum Field Theory, eq. (5.15); or Halzen & Martin, Quarks and Leptons, eq. (6.23). \(\theta\) is the muon’s polar angle relative to the beam in the center-of-mass frame, and \(s=(\sqrt s)^2\) is the Mandelstam invariant (see physicskit.particle.scattering.mandelstam_s()).

Simplification: electron and muon masses are neglected (valid for \(\sqrt s\gg2m_\mu\approx0.2\) GeV); higher-order QED corrections and any \(Z\)-boson contribution (relevant near and above the \(Z\) pole, \(\sqrt s\sim91\) GeV) are not included.

Parameters:

  • cos_theta (array_like) – Cosine of the muon’s center-of-mass polar angle.

  • sqrt_s (float) – Center-of-mass energy.

  • alpha (float, default=ALPHA_FS) – Fine-structure constant.

Returns:

ndarray or float

Examples

>>> round(float(qed_dsigma_domega_mumu(0.0, sqrt_s=10.0) / qed_dsigma_domega_mumu(1.0, sqrt_s=10.0)), 6)
0.5
physicskit.particle.qed_total_cross_section_mumu(sqrt_s, alpha=0.0072973525692838015)[source]#

Total leading-order QED cross section for \(e^+e^-\to\mu^+\mu^-\).

\[\sigma = \frac{4\pi\alpha^2}{3s}\]

obtained by integrating qed_dsigma_domega_mumu() over the full solid angle; the same simplifications apply (massless fermions, no \(Z\) contribution).

Parameters:

  • sqrt_s (float) – Center-of-mass energy.

  • alpha (float, default=ALPHA_FS) – Fine-structure constant.

Returns:

float

Examples

>>> round(qed_total_cross_section_mumu(10.0), 10) == round(4 * np.pi * ALPHA_FS**2 / (3 * 100.0), 10)
True
physicskit.particle.radioactive_decay_number(N0, decay_constant_, t)[source]#

Surviving population, \(N(t)=N_0e^{-\lambda t}\).

Parameters:

  • N0 (float) – Initial population.

  • decay_constant_ – Decay constant \(\lambda\).

  • t (array_like) – Time(s), same units as \(1/\lambda\).

Returns:

ndarray or float

Examples

>>> round(float(radioactive_decay_number(100.0, 1.0, 0.693147)), 4)
50.0
physicskit.particle.rapidity(four_vector, axis='z')[source]#

Rapidity along a Cartesian axis, \(y=\tfrac12\ln\!\big[(E+p_i)/(E-p_i)\big]\).

Rapidity is additive under boosts along the same axis: boosting a particle by velocity \(\beta\) shifts its rapidity by exactly \(\operatorname{artanh}\beta\).

Parameters:

  • four_vector (FourVector) – The particle.

  • axis ({"x", "y", "z"}, default="z") – Cartesian axis.

Returns:

float

Examples

>>> p = FourVector(2.0, 0.0, 0.0, 1.0)
>>> round(rapidity(p), 6)
0.549306
physicskit.particle.rutherford_dsigma_domega(theta, Z1, Z2, E_kin, alpha=0.0072973525692838015)[source]#

Rutherford differential cross section, in natural units.

\[\frac{d\sigma}{d\Omega} = \left(\frac{Z_1Z_2\alpha}{4E_{\rm kin}}\right)^2 \frac{1}{\sin^4(\theta/2)}\]
Parameters:

  • theta (array_like) – Scattering angle (radians). theta=0 gives inf – the genuine, integrable-but-divergent small-angle Rutherford behavior.

  • Z1 (float) – Projectile and target charge numbers.

  • Z2 (float) – Projectile and target charge numbers.

  • E_kin (float) – Projectile kinetic energy.

  • alpha (float, default=ALPHA_FS) – Coupling constant (the fine-structure constant, by default).

Returns:

ndarray or float

Examples

>>> import numpy as np
>>> theta = np.array([np.pi / 4, np.pi / 2])
>>> vals = rutherford_dsigma_domega(theta, Z1=1, Z2=79, E_kin=5.0)
>>> ratio = (vals[0] * np.sin(theta[0] / 2) ** 4) / (vals[1] * np.sin(theta[1] / 2) ** 4)
>>> round(float(ratio), 6)
1.0
physicskit.particle.sample_michel_electron_energies(n, rng=None)[source]#

Sample n scaled electron energies \(x=2E_e/m_\mu\) from the Michel spectrum.

Rejection sampling against the shape’s bound: \(x^2(3-2x)\le1\) on \([0,1]\) (attained at \(x=1\)), so a uniform proposal on \([0,1]\times[0,1]\) accepted below the (unnormalized) curve \(x^2(3-2x)\) reproduces michel_spectrum().

Parameters:

  • n (int) – Number of samples.

  • rng (numpy.random.Generator, optional) – Random number generator; a fresh default one is used if omitted.

Returns:

ndarray of shape (n,)

Examples

>>> rng = np.random.default_rng(0)
>>> xs = sample_michel_electron_energies(20000, rng=rng)
>>> bool(0.68 < xs.mean() < 0.72)  # theoretical mean <x> = 7/10
True
physicskit.particle.semf_binding_energy(Z, A)[source]#

Semi-empirical (Weizsacker) total nuclear binding energy, in MeV.

\[B(Z,A) = a_V A - a_S A^{2/3} - a_C\frac{Z(Z-1)}{A^{1/3}} - a_A\frac{(A-2Z)^2}{A} + \delta(A,Z)\]

with the standard coefficients \(a_V=15.75\), \(a_S=17.8\), \(a_C=0.711\), \(a_A=23.7\) MeV, and the pairing term \(\delta=+a_PA^{-1/2}\) for even-even nuclei, \(\delta=-a_PA^{-1/2}\) for odd-odd nuclei, and \(\delta=0\) for odd-\(A\) nuclei, with \(a_P\approx11.18\) MeV.

Parameters:

  • Z (int) – Atomic (proton) number.

  • A (int) – Mass number.

Returns:

float – Total binding energy, in MeV.

Examples

>>> B = semf_binding_energy(26, 56)  # Iron-56
>>> 8.0 < B / 56 < 9.0
True
physicskit.particle.shower_leaves(root)[source]#

The final-state (non-branching) particles of a shower tree.

Parameters:

root (ShowerParticle)

Returns:

list of ShowerParticle

physicskit.particle.simple_shower(E0, m0_fraction=0.3, daughter_fraction=0.42, E_threshold=1.0, max_generations=8, rng=None)[source]#

A toy branching cascade: a single high-energy particle splitting repeatedly until its descendants fall below an energy threshold.

Each split is an exact 1-to-2 branching (see the module docstring for the shared construction); the initial particle is given a toy “virtuality” \(m_0=\) m0_fraction \({}\times E_0\) representing how far off-shell it is when it enters the cascade, and each generation’s daughters inherit a virtuality shrunk by daughter_fraction (< 0.5, so it always halves-or-better, guaranteeing the branching stays kinematically allowed).

Parameters:

  • E0 (float) – Energy of the initial incident particle (along the beam/z axis).

  • m0_fraction (float, default=0.3) – Initial virtuality as a fraction of E0.

  • daughter_fraction (float, default=0.42) – Fraction of the parent’s virtuality inherited by each daughter.

  • E_threshold (float, default=1.0) – Energy below which a particle is treated as final-state (stops branching).

  • max_generations (int, default=8) – Hard cap on branching depth (safety net against numerical edge cases).

  • rng (numpy.random.Generator, optional) – Random number generator; a fresh default one is used if omitted.

Returns:

ShowerParticle – The root of the cascade tree (the initial particle; its descendants are reached via .children).

Examples

>>> import numpy as np
>>> root = simple_shower(100.0, E_threshold=5.0, rng=np.random.default_rng(0))
>>> leaves = shower_leaves(root)
>>> from physicskit.particle.kinematics import invariant_mass
>>> bool(abs(invariant_mass([leaf.four_vector for leaf in leaves]) - root.four_vector.mass) < 1e-6)
True
physicskit.particle.string_break_chain(t, v, kappa, m_q, r0=0.0, n_breaks=3)[source]#

Sequence of string-breaking events as a receding quark pair separates.

The original quark pair separates at the prescribed constant speed v, \(r(t)=r_0+vt\). The k-th break occurs when the total stored string energy \(\kappa r(t)\) first reaches \(2m_qk\) – the energy needed for each of the k string segments existing just before that break to independently pair-produce a light \(q\bar q\) pair at threshold (each of the k segments is assumed, for this bookkeeping, to carry an equal share \(\kappa r(t)/k\) of the total stored energy). After the break there are k+1 segments.

Parameters:

  • t (array_like) – Times at which to evaluate the state of the system.

  • v (float) – Constant recession speed of the original quark pair, \(0<v<1\).

  • kappa (float) – String tension.

  • m_q (float) – Constituent mass of the light quark pair created at each break (the pair-production threshold for one string segment is \(2m_q\)).

  • r0 (float, default=0.0) – Initial quark-antiquark separation at \(t=0\).

  • n_breaks (int, default=3) – Maximum number of sequential breaks to compute.

Returns:

dict – t, r (total extent), energy_total (\(\kappa r\)), n_segments (number of string segments at each time), energy_per_segment, break_times, r_break_unit (\(2m_q/\kappa\), the extent of a single fresh segment at the instant it breaks), kappa, m_q.

Examples

>>> import numpy as np
>>> sim = string_break_chain(np.linspace(0, 10, 5), v=0.5, kappa=1.0, m_q=1.0, n_breaks=2)
>>> sim["break_times"]
array([4., 8.])
physicskit.particle.string_tension_energy(r, kappa)[source]#

Energy stored in a confining linear-potential string, \(V(r)=\kappa r\).

Parameters:

  • r (array_like) – Quark-antiquark separation.

  • kappa (float) – String tension.

Returns:

ndarray or float

Examples

>>> float(string_tension_energy(2.0, 0.9))
1.8
physicskit.particle.survival_probability(L, E, theta, delta_m2)[source]#

Two-flavor vacuum survival probability, \(P(\nu_e\to\nu_e)=1-P(\nu_e\to\nu_\mu)\).

Parameters:
Returns:

ndarray or float

Examples

>>> round(float(survival_probability(0.0, 1.0, 0.5, 2.5e-3)), 10)
1.0
physicskit.particle.two_body_decay(M, m1, m2, cos_theta, phi)[source]#

The two daughters’ four-momenta in the parent’s rest frame.

Parameters:

  • M (float) – Parent mass.

  • m1 (float) – Daughter masses.

  • m2 (float) – Daughter masses.

  • cos_theta (float) – Cosine of daughter 1’s polar angle from the z-axis.

  • phi (float) – Daughter 1’s azimuthal angle.

Returns:

p1, p2 (FourVector) – The two daughters, back-to-back in the parent’s rest frame.

Examples

>>> p1, p2 = two_body_decay(1.0, 0.4, 0.4, cos_theta=1.0, phi=0.0)
>>> round(p1.E, 6), round(p1.pz, 6)
(0.5, 0.3)
physicskit.particle.two_body_decay_momentum(M, m1, m2)[source]#

Daughter momentum magnitude in a two-body decay, in the parent’s rest frame.

\[p^* = \frac{\sqrt{\lambda(M^2,m_1^2,m_2^2)}}{2M}, \qquad \lambda(a,b,c) = a^2+b^2+c^2-2ab-2bc-2ca\]
Parameters:

  • M (float) – Parent mass.

  • m1 (float) – Daughter masses.

  • m2 (float) – Daughter masses.

Returns:

float

Raises:

ValueError – If M < m1 + m2 (energetically forbidden decay).

Examples

>>> round(two_body_decay_momentum(1.0, 0.4, 0.4), 6)
0.3

Relativistic four-vector kinematics.

Works throughout in natural units with \(c=1\): energy, momentum, and mass all share the same energy-like unit (GeV for particle physics, MeV for nuclear physics – the caller’s choice, so long as it is used consistently), the same convention physicskit.relativity uses for its own geometrized units (\(G=c=1\)).

  • FourVector – an energy-momentum four-vector \(p^\mu = (E, p_x, p_y, p_z)\), with its invariant mass, speed, and Lorentz factor as derived properties.

  • boost() – a Lorentz boost along a Cartesian axis.

  • boost_generic() – a Lorentz boost along an arbitrary 3-velocity.

  • invariant_mass() – the invariant mass of a system of particles.

  • rapidity() – the additive (under boosts) rapidity variable.

  • boost_to_com() – the velocity of a system’s center-of-momentum frame.

class physicskit.particle.kinematics.FourVector(E, px, py, pz)[source]#

Bases: object

An energy-momentum four-vector \(p^\mu=(E,p_x,p_y,p_z)\).

Parameters:

  • E (float) – Energy.

  • px (float) – Momentum components.

  • py (float) – Momentum components.

  • pz (float) – Momentum components.

Examples

>>> p = FourVector(2.0, 0.0, 0.0, 1.0)
>>> round(p.mass, 6)
1.732051
property beta#

Speed, \(\beta=|\vec p|/E\).

property gamma#

Lorentz factor, \(\gamma=E/m\).

Raises:

ZeroDivisionError – If the four-vector is massless (no rest frame).

property mass#

Invariant mass, \(m=\sqrt{E^2-|\vec p|^2}\).

The radicand is clipped at 0 before the square root, so a four-vector that is numerically spacelike only by floating-point noise (e.g. the sum of several on-shell four-vectors) returns 0.0 rather than raising on a tiny negative argument.

property p_mag#

Momentum magnitude, \(|\vec p|\).

property p_vec#

the momentum 3-vector.

Type:

ndarray of shape (3,)

physicskit.particle.kinematics.boost(four_vector, beta, axis='z')[source]#

Lorentz-boost a FourVector along a Cartesian axis.

Parameters:

  • four_vector (FourVector) – The four-vector to boost.

  • beta (float) – Boost velocity, \(-1 < \beta < 1\), in units of \(c=1\).

  • axis ({"x", "y", "z"}, default="z") – Cartesian axis of the boost.

Returns:

FourVector – The boosted four-vector.

Examples

>>> p = FourVector(1.0, 0.0, 0.0, 0.0)  # a particle at rest, mass 1
>>> pb = boost(p, 0.6, axis="z")
>>> round(pb.E, 6), round(pb.pz, 6)
(1.25, 0.75)
physicskit.particle.kinematics.boost_generic(four_vector, beta_vec)[source]#

Lorentz-boost a FourVector along an arbitrary 3-velocity.

Generalizes boost() to a boost direction that need not lie along a Cartesian axis – the transformation used, e.g., to take a daughter four-momentum computed in a parent’s rest frame back to the lab frame when the parent’s own flight direction is not one of the axes (as happens at every vertex of a branching cascade; see physicskit.particle.collider).

\[E' = \gamma(E + \vec\beta\cdot\vec p), \qquad \vec p' = \vec p + \left[\frac{\gamma-1}{\beta^2}(\vec\beta\cdot\vec p) + \gamma E\right]\vec\beta\]

which reduces to boost() when \(\vec\beta\) is aligned with a single axis.

Parameters:

  • four_vector (FourVector) – The four-vector to boost.

  • beta_vec (array_like of shape (3,)) – Boost velocity, \(|\vec\beta| < 1\), in units of \(c=1\).

Returns:

FourVector – The boosted four-vector.

Examples

>>> p = FourVector(1.0, 0.0, 0.0, 0.0)  # a particle at rest, mass 1
>>> pb = boost_generic(p, [0.0, 0.0, 0.6])
>>> round(pb.E, 6), round(pb.pz, 6)
(1.25, 0.75)
physicskit.particle.kinematics.boost_to_com(four_vectors)[source]#

The 3-velocity of the center-of-momentum frame of a system.

Parameters:

four_vectors (iterable of FourVector) – The particles making up the system.

Returns:

ndarray of shape (3,) – The velocity (units of \(c=1\)) of the system’s total momentum relative to the current frame, \(\vec\beta_{\rm com} = \vec p_{\rm tot}/E_{\rm tot}\).

Examples

>>> p1 = FourVector(2.0, 0.0, 0.0, 1.0)
>>> p2 = FourVector(1.0, 0.0, 0.0, 0.0)
>>> beta_com = boost_to_com([p1, p2])
>>> [round(float(b), 6) for b in beta_com]
[0.0, 0.0, 0.333333]
physicskit.particle.kinematics.invariant_mass(four_vectors)[source]#

The invariant mass of a system of four-vectors.

Parameters:

four_vectors (iterable of FourVector) – The particles making up the system.

Returns:

float – The mass of the sum of the four-vectors – e.g. the reconstructed mass of a resonance from its decay products.

Examples

>>> p1 = FourVector(1.0, 0.0, 0.0, 0.6)
>>> p2 = FourVector(1.0, 0.0, 0.0, -0.6)
>>> round(invariant_mass([p1, p2]), 6)
2.0
physicskit.particle.kinematics.rapidity(four_vector, axis='z')[source]#

Rapidity along a Cartesian axis, \(y=\tfrac12\ln\!\big[(E+p_i)/(E-p_i)\big]\).

Rapidity is additive under boosts along the same axis: boosting a particle by velocity \(\beta\) shifts its rapidity by exactly \(\operatorname{artanh}\beta\).

Parameters:

  • four_vector (FourVector) – The particle.

  • axis ({"x", "y", "z"}, default="z") – Cartesian axis.

Returns:

float

Examples

>>> p = FourVector(2.0, 0.0, 0.0, 1.0)
>>> round(rapidity(p), 6)
0.549306

Two-body decay kinematics and radioactive decay chains.

Uses the same natural-unit convention (\(c=1\)) as physicskit.particle.kinematics.

physicskit.particle.decays.activity(N, decay_constant_)[source]#

Decay rate (activity), \(A=\lambda N\).

Examples

>>> activity(100.0, 0.01)
1.0
physicskit.particle.decays.bateman_decay_chain(N0, decay_constants, t)[source]#

Population of every species in a linear radioactive decay chain.

Solves the Bateman equations for a chain \(1\to2\to\cdots\to n\) with decay constants decay_constants (species 1 starts with population N0, all daughters start at zero):

\[N_k(t) = N_0\Big(\prod_{i=1}^{k-1}\lambda_i\Big) \sum_{i=1}^{k}\frac{e^{-\lambda_i t}} {\prod_{j=1,j\neq i}^{k}(\lambda_j-\lambda_i)}\]

with \(N_1(t)=N_0e^{-\lambda_1 t}\). A zero decay constant for the final species represents a stable end product.

Parameters:

  • N0 (float) – Initial population of species 1.

  • decay_constants (sequence of float) – Decay constant of each of the n species in the chain. Must be pairwise distinct among the species actually populated (the closed-form sum divides by lambda_j - lambda_i).

  • t (array_like) – Times at which to evaluate the populations.

Returns:

ndarray of shape (n, len(t)) – Population of each species over time.

Examples

>>> t = np.array([0.0, 1.0])
>>> N = bateman_decay_chain(100.0, [1.0, 2.0], t)
>>> N.shape
(2, 2)
>>> round(float(N[0, 0]), 6)
100.0
physicskit.particle.decays.decay_constant(half_life_)[source]#

Decay constant from half-life, \(\lambda=\ln2/t_{1/2}\).

Examples

>>> round(decay_constant(np.log(2)), 6)
1.0
physicskit.particle.decays.half_life(decay_constant_)[source]#

Half-life from decay constant, \(t_{1/2}=\ln2/\lambda\).

Examples

>>> round(half_life(1.0), 6)
0.693147
physicskit.particle.decays.michel_spectrum(x)[source]#

The (unpolarized, massless-electron) Michel spectrum shape for muon decay.

\[\frac{d\Gamma}{dx} = 2x^2(3-2x), \qquad x=\frac{2E_e}{m_\mu}\in[0,1]\]

the standard leading-order result for \(\mu^-\to e^-\bar\nu_\mu\nu_e\) from the V-A four-fermion weak interaction (Michel 1950; see e.g. Commins & Bucksbaum, Weak Interactions of Leptons and Quarks, Ch. 4, or the muon-decay chapter of any standard particle physics text), in the limit \(m_e/m_\mu\to0\) and for an unpolarized muon (no Michel parameter \(\rho,\eta,\xi,\delta\) dependence – those parametrize corrections from polarization and a nonzero electron mass/anomalous couplings, all set to their V-A/Standard-Model, unpolarized values here). Normalized so that \(\int_0^1 (d\Gamma/dx)\,dx = 1\).

Parameters:

x (array_like) – Scaled electron energy, \(x=2E_e/m_\mu\in[0,1]\).

Returns:

ndarray or float

Examples

>>> import numpy as np
>>> x = np.linspace(0, 1, 100001)
>>> round(float(np.trapezoid(michel_spectrum(x), x)), 4)
1.0
physicskit.particle.decays.muon_decay_event(m_mu, rng=None)[source]#

Sample one \(\mu^-\to e^-+\bar\nu_\mu+\nu_e\) decay in the muon rest frame.

Toy-but-exact three-body kinematics built from two chained exact two-body decays (reusing two_body_decay(), so energy-momentum is conserved to machine precision at both vertices):

  1. The electron’s energy is drawn from the Michel spectrum (massless electron approximation) via sample_michel_electron_energies(), fixing the invariant mass of the recoiling neutrino pair through \(m_{\nu\nu}^2=m_\mu^2-2m_\mu E_e\) (exact for \(m_e=0\)). two_body_decay() then gives the electron and the neutrino-pair “system” four-momenta, back to back in the muon rest frame, with an isotropically sampled decay angle (the standard unpolarized-muon assumption).

  2. The neutrino-pair system is split into its two (massless) neutrinos isotropically in its own rest frame via a second call to two_body_decay(), then boosted back to the muon rest frame with boost_generic().

Both neutrino masses are set to zero – an excellent approximation given their sub-eV masses compared to \(m_\mu\sim100\) MeV.

Parameters:

  • m_mu (float) – Muon mass (in whatever energy unit the caller uses consistently).

  • rng (numpy.random.Generator, optional) – Random number generator; a fresh default one is used if omitted.

Returns:

p_e, p_numu_bar, p_nue (FourVector) – The three decay products, in the muon rest frame.

Examples

>>> from physicskit.particle.kinematics import invariant_mass
>>> rng = np.random.default_rng(1)
>>> p_e, p_numu_bar, p_nue = muon_decay_event(105.658, rng=rng)
>>> round(invariant_mass([p_e, p_numu_bar, p_nue]), 6)
105.658
physicskit.particle.decays.radioactive_decay_number(N0, decay_constant_, t)[source]#

Surviving population, \(N(t)=N_0e^{-\lambda t}\).

Parameters:

  • N0 (float) – Initial population.

  • decay_constant_ – Decay constant \(\lambda\).

  • t (array_like) – Time(s), same units as \(1/\lambda\).

Returns:

ndarray or float

Examples

>>> round(float(radioactive_decay_number(100.0, 1.0, 0.693147)), 4)
50.0
physicskit.particle.decays.sample_michel_electron_energies(n, rng=None)[source]#

Sample n scaled electron energies \(x=2E_e/m_\mu\) from the Michel spectrum.

Rejection sampling against the shape’s bound: \(x^2(3-2x)\le1\) on \([0,1]\) (attained at \(x=1\)), so a uniform proposal on \([0,1]\times[0,1]\) accepted below the (unnormalized) curve \(x^2(3-2x)\) reproduces michel_spectrum().

Parameters:

  • n (int) – Number of samples.

  • rng (numpy.random.Generator, optional) – Random number generator; a fresh default one is used if omitted.

Returns:

ndarray of shape (n,)

Examples

>>> rng = np.random.default_rng(0)
>>> xs = sample_michel_electron_energies(20000, rng=rng)
>>> bool(0.68 < xs.mean() < 0.72)  # theoretical mean <x> = 7/10
True
physicskit.particle.decays.two_body_decay(M, m1, m2, cos_theta, phi)[source]#

The two daughters’ four-momenta in the parent’s rest frame.

Parameters:

  • M (float) – Parent mass.

  • m1 (float) – Daughter masses.

  • m2 (float) – Daughter masses.

  • cos_theta (float) – Cosine of daughter 1’s polar angle from the z-axis.

  • phi (float) – Daughter 1’s azimuthal angle.

Returns:

p1, p2 (FourVector) – The two daughters, back-to-back in the parent’s rest frame.

Examples

>>> p1, p2 = two_body_decay(1.0, 0.4, 0.4, cos_theta=1.0, phi=0.0)
>>> round(p1.E, 6), round(p1.pz, 6)
(0.5, 0.3)
physicskit.particle.decays.two_body_decay_momentum(M, m1, m2)[source]#

Daughter momentum magnitude in a two-body decay, in the parent’s rest frame.

\[p^* = \frac{\sqrt{\lambda(M^2,m_1^2,m_2^2)}}{2M}, \qquad \lambda(a,b,c) = a^2+b^2+c^2-2ab-2bc-2ca\]
Parameters:

  • M (float) – Parent mass.

  • m1 (float) – Daughter masses.

  • m2 (float) – Daughter masses.

Returns:

float

Raises:

ValueError – If M < m1 + m2 (energetically forbidden decay).

Examples

>>> round(two_body_decay_momentum(1.0, 0.4, 0.4), 6)
0.3

Cross-sections and Mandelstam kinematics.

Uses the same natural-unit convention (\(\hbar=c=1\)) as physicskit.particle.kinematics.

physicskit.particle.scattering.ALPHA_FS = 0.0072973525692838015#

The fine-structure constant, \(\alpha\approx1/137.036\).

physicskit.particle.scattering.impact_parameter(theta, Z1, Z2, E_kin, alpha=0.0072973525692838015)[source]#

Classical Rutherford impact parameter producing scattering angle theta.

\[b(\theta) = \frac{Z_1Z_2\alpha}{2E_{\rm kin}}\cot(\theta/2)\]
Parameters:

  • theta (array_like) – Scattering angle (radians).

  • Z1 (float) – Projectile and target charge numbers.

  • Z2 (float) – Projectile and target charge numbers.

  • E_kin (float) – Projectile kinetic energy.

  • alpha (float, default=ALPHA_FS) – Coupling constant.

Returns:

ndarray or float

Examples

>>> round(float(impact_parameter(np.pi, Z1=1, Z2=79, E_kin=5.0)), 6)
0.0
physicskit.particle.scattering.mandelstam_s(p1, p2)[source]#

Mandelstam \(s=(p_1+p_2)^2\), the total invariant mass squared.

Parameters:
Returns:

float

Examples

>>> from physicskit.particle.kinematics import FourVector
>>> p1 = FourVector(1.0, 0.0, 0.0, 0.6)
>>> p2 = FourVector(1.0, 0.0, 0.0, -0.6)
>>> round(mandelstam_s(p1, p2), 6)
4.0
physicskit.particle.scattering.mandelstam_t(p1, p3)[source]#

Mandelstam \(t=(p_1-p_3)^2\), the momentum-transfer invariant.

Parameters:
Returns:

float

physicskit.particle.scattering.mandelstam_u(p1, p4)[source]#

Mandelstam \(u=(p_1-p_4)^2\), the other momentum-transfer invariant.

Parameters:
Returns:

float

Notes

For a 2-to-2 process \(1+2\to3+4\), \(s+t+u=\sum_i m_i^2\).

physicskit.particle.scattering.rutherford_dsigma_domega(theta, Z1, Z2, E_kin, alpha=0.0072973525692838015)[source]#

Rutherford differential cross section, in natural units.

\[\frac{d\sigma}{d\Omega} = \left(\frac{Z_1Z_2\alpha}{4E_{\rm kin}}\right)^2 \frac{1}{\sin^4(\theta/2)}\]
Parameters:

  • theta (array_like) – Scattering angle (radians). theta=0 gives inf – the genuine, integrable-but-divergent small-angle Rutherford behavior.

  • Z1 (float) – Projectile and target charge numbers.

  • Z2 (float) – Projectile and target charge numbers.

  • E_kin (float) – Projectile kinetic energy.

  • alpha (float, default=ALPHA_FS) – Coupling constant (the fine-structure constant, by default).

Returns:

ndarray or float

Examples

>>> import numpy as np
>>> theta = np.array([np.pi / 4, np.pi / 2])
>>> vals = rutherford_dsigma_domega(theta, Z1=1, Z2=79, E_kin=5.0)
>>> ratio = (vals[0] * np.sin(theta[0] / 2) ** 4) / (vals[1] * np.sin(theta[1] / 2) ** 4)
>>> round(float(ratio), 6)
1.0

Nuclear binding energy and reaction Q-values.

Masses and energies here are in MeV, the universal convention for the semi-empirical mass formula’s fitted coefficients.

physicskit.particle.nuclear.binding_energy_per_nucleon(Z, A)[source]#

Binding energy per nucleon, semf_binding_energy(Z, A) / A, in MeV.

Examples

>>> bpn = binding_energy_per_nucleon(26, 56)
>>> 8.0 < bpn < 9.0
True
physicskit.particle.nuclear.q_value(reactant_masses, product_masses)[source]#

Reaction Q-value, \(Q=\sum m_{\rm reactants}-\sum m_{\rm products}\).

A positive Q means the reaction is exothermic (releases energy); a negative Q means it is endothermic and requires that much additional kinetic energy to proceed. Masses (or mass-energies) may be in any consistent unit.

Parameters:

  • reactant_masses (iterable of float) – Masses of the reactants.

  • product_masses (iterable of float) – Masses of the products.

Returns:

float

Examples

>>> round(q_value([2.014102, 3.016049], [4.002602, 1.008665]), 6)  # D + T -> He4 + n
0.018884
physicskit.particle.nuclear.semf_binding_energy(Z, A)[source]#

Semi-empirical (Weizsacker) total nuclear binding energy, in MeV.

\[B(Z,A) = a_V A - a_S A^{2/3} - a_C\frac{Z(Z-1)}{A^{1/3}} - a_A\frac{(A-2Z)^2}{A} + \delta(A,Z)\]

with the standard coefficients \(a_V=15.75\), \(a_S=17.8\), \(a_C=0.711\), \(a_A=23.7\) MeV, and the pairing term \(\delta=+a_PA^{-1/2}\) for even-even nuclei, \(\delta=-a_PA^{-1/2}\) for odd-odd nuclei, and \(\delta=0\) for odd-\(A\) nuclei, with \(a_P\approx11.18\) MeV.

Parameters:

  • Z (int) – Atomic (proton) number.

  • A (int) – Mass number.

Returns:

float – Total binding energy, in MeV.

Examples

>>> B = semf_binding_energy(26, 56)  # Iron-56
>>> 8.0 < B / 56 < 9.0
True

Plotting and animation helpers for physicskit.particle.

  • physicskit.particle.visualizers.static – decay-chain and differential-cross-section line plots.

  • physicskit.particle.visualizers.animations – animations of the toy shower/confinement/electroweak/neutrino/decay demonstrations built on physicskit.particle’s collider, confinement, electroweak, and neutrinos modules.

physicskit.particle.visualizers.animate_cp_asymmetry(t, delta_m, gamma_s, gamma_l, epsilon, interval=60, ax=None)[source]#

Animate the neutral-meson decay-rate CP asymmetry \(A(t)\) building up, from physicskit.particle.electroweak.cp_asymmetry().

Parameters:

  • t (ndarray) – Times.

  • delta_m – See physicskit.particle.electroweak.meson_decay_rates_cp_eigenstate().

  • gamma_s – See physicskit.particle.electroweak.meson_decay_rates_cp_eigenstate().

  • gamma_l – See physicskit.particle.electroweak.meson_decay_rates_cp_eigenstate().

  • epsilon – See physicskit.particle.electroweak.meson_decay_rates_cp_eigenstate().

  • interval (int, default=60) – Delay between frames, in milliseconds.

  • ax (matplotlib.axes.Axes, optional) – Axes to draw into; a new figure is created if omitted.

Returns:

matplotlib.animation.FuncAnimation

Examples

>>> import numpy as np
>>> t = np.linspace(0, 20, 100)
>>> anim = animate_cp_asymmetry(t, delta_m=0.5, gamma_s=1.0, gamma_l=0.1, epsilon=0.002 + 0.001j)
>>> anim.__class__.__name__
'FuncAnimation'
physicskit.particle.visualizers.animate_decay_chain_bars(t, populations, labels=None, interval=150, ax=None)[source]#

Animate each species’ population evolving over time in a decay chain.

Parameters:

  • t (ndarray of shape (n_t,)) – Times.

  • populations (ndarray of shape (n_species, n_t)) – Population of each species over time, e.g. from physicskit.particle.decays.bateman_decay_chain().

  • labels (sequence of str, optional) – Bar labels; defaults to "species 1", "species 2", …

  • interval (int, default=150) – Delay between frames, in milliseconds.

  • ax (matplotlib.axes.Axes, optional) – Axes to draw into; a new figure is created if omitted.

Returns:

matplotlib.animation.FuncAnimation

Examples

>>> import numpy as np
>>> from physicskit.particle.decays import bateman_decay_chain
>>> t = np.linspace(0, 10, 30)
>>> N = bateman_decay_chain(1000.0, [0.5, 0.2], t)
>>> anim = animate_decay_chain_bars(t, N)
>>> anim.__class__.__name__
'FuncAnimation'
physicskit.particle.visualizers.animate_detector_event(four_vectors, charges, B, detector_radii=(3.0, 6.0), n_frames=60, interval=100, ax=None)[source]#

Animate final-state tracks growing outward from the vertex and landing in a schematic calorimeter.

Charged particles curve (radius \(r=p_T/(qB)\), see physicskit.particle.collider.charged_track_points()); neutral particles fly in straight lines. Each frame reveals one more step along every track; a marker is drawn at each track’s current endpoint, standing in for a calorimeter energy deposit once the track reaches the outer radius.

Parameters:

  • four_vectors (sequence of physicskit.particle.kinematics.FourVector) – Final-state particles.

  • charges (sequence of float) – Charge of each particle, in units of the elementary charge.

  • B (float) – Magnetic field strength (see charged_track_points()).

  • detector_radii ((float, float), default=(3.0, 6.0)) – Radii of the schematic tracker and calorimeter circles.

  • n_frames (int, default=60) – Number of animation frames (and points sampled along each track).

  • interval (int, default=100) – Delay between frames, in milliseconds.

  • ax (matplotlib.axes.Axes, optional) – Axes to draw into; a new figure is created if omitted.

Returns:

matplotlib.animation.FuncAnimation

Examples

>>> from physicskit.particle.kinematics import FourVector
>>> ps = [FourVector(5.0, 3.0, 0.0, 3.0), FourVector(4.0, 0.0, 3.0, 2.0)]
>>> anim = animate_detector_event(ps, charges=[1.0, -1.0], B=1.0)
>>> anim.__class__.__name__
'FuncAnimation'
physicskit.particle.visualizers.animate_higgs_rollover(phi_t, a, b, interval=30, ax=None)[source]#

Animate a scalar field rolling from the unstable symmetric point into a broken-symmetry vacuum, from physicskit.particle.electroweak.higgs_field_rollover().

The static curve is physicskit.particle.electroweak.higgs_potential(); a point (with a fading trail) traces \((\phi(t), V(\phi(t)))\) as the field rolls down one side of the double well.

Parameters:

  • phi_t (ndarray) – Field trajectory \(\phi(t)\), e.g. the first output of higgs_field_rollover().

  • a (float) – Potential parameters, see higgs_potential().

  • b (float) – Potential parameters, see higgs_potential().

  • interval (int, default=30) – Delay between frames, in milliseconds.

  • ax (matplotlib.axes.Axes, optional) – Axes to draw into; a new figure is created if omitted.

Returns:

matplotlib.animation.FuncAnimation

Examples

>>> import numpy as np
>>> from physicskit.particle.electroweak import higgs_field_rollover
>>> t = np.linspace(0, 40, 100)
>>> phi, _ = higgs_field_rollover(1e-3, 0.0, a=1.0, b=1.0, t_eval=t, damping=0.05)
>>> anim = animate_higgs_rollover(phi, a=1.0, b=1.0)
>>> anim.__class__.__name__
'FuncAnimation'
physicskit.particle.visualizers.animate_michel_histogram(x_samples, batch_size=20, n_bins=30, interval=100, ax=None)[source]#

Animate a histogram of sampled Michel-spectrum electron energies building up event-by-event, converging to physicskit.particle.decays.michel_spectrum().

Parameters:

  • x_samples (ndarray) – Sampled scaled electron energies, e.g. from physicskit.particle.decays.sample_michel_electron_energies().

  • batch_size (int, default=20) – Number of additional samples revealed per frame.

  • n_bins (int, default=30) – Number of histogram bins.

  • interval (int, default=100) – Delay between frames, in milliseconds.

  • ax (matplotlib.axes.Axes, optional) – Axes to draw into; a new figure is created if omitted.

Returns:

matplotlib.animation.FuncAnimation

Examples

>>> import numpy as np
>>> from physicskit.particle.decays import sample_michel_electron_energies
>>> xs = sample_michel_electron_energies(500, rng=np.random.default_rng(0))
>>> anim = animate_michel_histogram(xs, batch_size=100)
>>> anim.__class__.__name__
'FuncAnimation'
physicskit.particle.visualizers.animate_neutrino_oscillation(L, E, theta, delta_m2, interval=60, ax=None)[source]#

Animate \(P(\nu_e)\) and \(P(\nu_\mu)\) growing as a function of distance traveled, from physicskit.particle.neutrinos.oscillation_probability().

Parameters:

  • L (ndarray) – Baselines (distance traveled), in km, in increasing order.

  • E (float) – Neutrino energy, in GeV.

  • theta (float) – Mixing angle, radians.

  • delta_m2 (float) – Mass-squared splitting, in \({\rm eV}^2\).

  • interval (int, default=60) – Delay between frames, in milliseconds.

  • ax (matplotlib.axes.Axes, optional) – Axes to draw into; a new figure is created if omitted.

Returns:

matplotlib.animation.FuncAnimation

Examples

>>> import numpy as np
>>> L = np.linspace(0, 1000, 100)
>>> anim = animate_neutrino_oscillation(L, E=1.0, theta=0.6, delta_m2=2.5e-3)
>>> anim.__class__.__name__
'FuncAnimation'
physicskit.particle.visualizers.animate_particle_cascade(root, detector_radii=(3.0, 6.0), track_length=0.6, interval=150, ax=None)[source]#

Animate a toy shower cascade (physicskit.particle.collider.simple_shower()) propagating outward from the collision vertex through a schematic detector.

Concentric circles mark a “tracker” (inner) and “calorimeter” (outer) region; each branch of the shower tree appears as one more track segment per frame, in order of production (generation).

Parameters:

  • root (physicskit.particle.collider.ShowerParticle) – Root of the shower tree, e.g. from physicskit.particle.collider.simple_shower().

  • detector_radii ((float, float), default=(3.0, 6.0)) – Radii of the schematic tracker and calorimeter circles.

  • track_length (float, default=0.6) – Schematic (not physically derived) drawn length of each track segment – purely a layout choice, not a decay-length calculation.

  • interval (int, default=150) – Delay between frames, in milliseconds.

  • ax (matplotlib.axes.Axes, optional) – Axes to draw into; a new figure is created if omitted.

Returns:

matplotlib.animation.FuncAnimation

Examples

>>> import numpy as np
>>> from physicskit.particle.collider import simple_shower
>>> root = simple_shower(50.0, E_threshold=5.0, rng=np.random.default_rng(0))
>>> anim = animate_particle_cascade(root)
>>> anim.__class__.__name__
'FuncAnimation'
physicskit.particle.visualizers.animate_parton_shower(root, detector_radii=(3.0, 6.0), track_length=0.6, n_jets=2, interval=150, ax=None)[source]#

Animate a toy parton shower (physicskit.particle.collider.parton_shower()) developing into jets.

Reuses the tree-reveal machinery of animate_particle_cascade(), color-coding branches by parton flavor (quark/antiquark vs. gluon) and overlaying the final jet-axis directions found by physicskit.particle.collider.cluster_into_jets().

Parameters:
Returns:

matplotlib.animation.FuncAnimation

Examples

>>> import numpy as np
>>> from physicskit.particle.collider import parton_shower
>>> root = parton_shower(50.0, E_threshold=5.0, rng=np.random.default_rng(0))
>>> anim = animate_parton_shower(root)
>>> anim.__class__.__name__
'FuncAnimation'
physicskit.particle.visualizers.animate_qed_angular_distribution(sqrt_s_values, alpha=0.0072973525692838015, n_theta=200, interval=200, ax=None)[source]#

Animate the \(e^+e^-\to\mu^+\mu^-\) angular distribution’s shape as \(\sqrt s\) is swept, from physicskit.particle.electroweak.qed_dsigma_domega_mumu().

Parameters:

  • sqrt_s_values (ndarray) – Center-of-mass energies to sweep through (one per frame).

  • alpha (float, default=ALPHA_FS) – Fine-structure constant.

  • n_theta (int, default=200) – Number of polar-angle samples per frame.

  • interval (int, default=200) – Delay between frames, in milliseconds.

  • ax (matplotlib.axes.Axes, optional) – Polar axes to draw into; new polar axes are created if omitted.

Returns:

matplotlib.animation.FuncAnimation

Examples

>>> import numpy as np
>>> sqrt_s = np.linspace(5.0, 20.0, 20)
>>> anim = animate_qed_angular_distribution(sqrt_s)
>>> anim.__class__.__name__
'FuncAnimation'
physicskit.particle.visualizers.animate_string_breaking(sim, interval=120, ax=None)[source]#

Animate a receding quark pair’s confining string stretching and breaking, from physicskit.particle.confinement.string_break_chain().

A thick line (whose width and color track the stored tension, energy_per_segment) stretches between the two receding quark markers; each time the simulation records a break, it snaps into one more, shorter segment, with a new marker appearing at each newly created quark-antiquark endpoint. (The endpoints of the n_segments(t) current strings are drawn evenly spaced across the current total extent – a schematic simplification, see string_break_chain()’s docstring, rather than tracking each new pair’s exact, individually frozen creation position.)

Parameters:

  • sim (dict) – Output of physicskit.particle.confinement.string_break_chain().

  • interval (int, default=120) – Delay between frames, in milliseconds.

  • ax (matplotlib.axes.Axes, optional) – Axes to draw into; a new figure is created if omitted.

Returns:

matplotlib.animation.FuncAnimation

Examples

>>> import numpy as np
>>> from physicskit.particle.confinement import string_break_chain
>>> sim = string_break_chain(np.linspace(0, 10, 40), v=0.3, kappa=1.0, m_q=1.0, n_breaks=2)
>>> anim = animate_string_breaking(sim)
>>> anim.__class__.__name__
'FuncAnimation'
physicskit.particle.visualizers.plot_decay_chain(t, populations, labels=None, ax=None)[source]#

Plot each species’ population over time in a radioactive decay chain.

Parameters:

  • t (ndarray of shape (n_t,)) – Times.

  • populations (ndarray of shape (n_species, n_t)) – Population of each species over time, e.g. from physicskit.particle.decays.bateman_decay_chain().

  • labels (sequence of str, optional) – Label for each species; defaults to "species 1", "species 2", …

  • 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.particle.decays import bateman_decay_chain
>>> t = np.linspace(0, 10, 50)
>>> N = bateman_decay_chain(1000.0, [0.5, 0.2], t)
>>> fig, ax = plot_decay_chain(t, N, labels=["Parent", "Daughter"])
>>> isinstance(fig, plt.Figure)
True
physicskit.particle.visualizers.plot_differential_cross_section(theta, dsigma_domega, ax=None, log_scale=True)[source]#

Plot a differential cross section \(d\sigma/d\Omega\) vs. scattering angle.

Parameters:

  • theta (ndarray of shape (n,)) – Scattering angles (radians).

  • dsigma_domega (ndarray of shape (n,)) – Differential cross section values, e.g. from physicskit.particle.scattering.rutherford_dsigma_domega().

  • ax (matplotlib.axes.Axes, optional) – Axes to draw into; a new figure is created if omitted.

  • log_scale (bool, default=True) – Plot the cross section on a log y-axis, appropriate for the Rutherford formula’s steep small-angle divergence.

Returns:

  • fig (matplotlib.figure.Figure)

  • ax (matplotlib.axes.Axes)

Examples

>>> import numpy as np
>>> from physicskit.particle.scattering import rutherford_dsigma_domega
>>> theta = np.linspace(0.1, np.pi - 0.1, 50)
>>> vals = rutherford_dsigma_domega(theta, Z1=2, Z2=79, E_kin=5.0)
>>> fig, ax = plot_differential_cross_section(theta, vals)
>>> isinstance(fig, plt.Figure)
True
On this page
Show Source