API Reference

All public classes and functions are documented here. Entries are generated automatically from source docstrings.

Quick reference

Class / function

Module

Purpose

MilkyWayAxionHalo

MilkyWayAxionHalo

SHM axion-wind model

axion_lineshape()

MilkyWayAxionHalo

Analytical PSD lineshape (static)

getAmpSpectra()

MilkyWayAxionHalo

Stochastic amplitude spectra

getRabiFreq()

MilkyWayAxionHalo

rms Rabi frequency

GravBoundAxionHalo

GravBoundAxionHalo

Generic TISE solver

EarthBoundAxionHalo

EarthBoundAxionHalo

Earth-specific TISE solver

Station

Station

Geographic lab location

Sample

Sample

NMR sample

Magnet

Apparatus

Static polarising magnet

SimuParams

SimuTypes

Simulation parameter dict (TypedDict)

SimuEntry

SimuTypes

Simulation + parameters pair

MagField

SimuTools

(Pseudo)magnetic field builder

Simulations

SimuTools

Multi-run simulation manager

Simulation

SimuTools

Single Bloch-equation simulation

axionbloch.MilkyWayAxionHalo

Models the axion dark-matter field from the Milky Way galactic halo using the Standard Halo Model (SHM) velocity distribution. The main class is MilkyWayAxionHalo, which provides:

  • Axion field properties: Compton frequency, quality factor, coherence time.

  • Analytical PSD lineshapes for three coupling geometries via the static method axion_lineshape().

  • Stochastic amplitude spectra for time-domain simulations via getAmpSpectra().

  • rms Rabi frequency via getRabiFreq().

.. py:module:: axionbloch.MilkyWayAxionHalo

.. py:class:: MilkyWayAxionHalo(name=’Milky Way Axion Halo’, nu_a: ~astropy.units.quantity.Quantity = None, m_a: ~astropy.units.quantity.Quantity = None, g_aNN: ~astropy.units.quantity.Quantity = None, Qa: ~astropy.units.quantity.Quantity = None, v_0: ~astropy.units.quantity.Quantity = <Quantity 220. km / s>, v_lab: ~astropy.units.quantity.Quantity = <Quantity 233. km / s>, windAngle: ~astropy.units.quantity.Quantity = None, rho_E_DM: ~astropy.units.quantity.Quantity = <Quantity 0.3 GeV / cm3>, verbose: bool = False) :module: axionbloch.MilkyWayAxionHalo

Bases: :py:class:object

Axion dark-matter field (axion wind) from the Milky Way halo.

Models the pseudomagnetic field experienced by nuclear spins due to the coherent oscillation of axion dark matter. Uses the Standard Halo Model (SHM) velocity distribution.

The class provides:

  • Axion field properties (Compton frequency, quality factor, coherence time).

  • Analytical lineshapes for gradient and non-gradient coupling cases (Gramolin et al.).

  • Stochastic amplitude spectra for time-domain simulations.

.. py:attribute:: MilkyWayAxionHalo.Qa :module: axionbloch.MilkyWayAxionHalo :type: ~astropy.units.quantity.Annotated[~astropy.units.quantity.Quantity, Unit(dimensionless)]

.. py:method:: MilkyWayAxionHalo.axion_lineshape(v_0: ~astropy.units.quantity.Annotated[~astropy.units.quantity.Quantity, Unit(“m / s”)], v_lab: ~astropy.units.quantity.Annotated[~astropy.units.quantity.Quantity, Unit(“m / s”)], nu_a: ~astropy.units.quantity.Annotated[~astropy.units.quantity.Quantity, Unit(“Hz”)], nu: ~astropy.units.quantity.Quantity, case: str = ‘non-grad’, alpha: ~astropy.units.quantity.Quantity = <Quantity 0. rad>, verbose: bool = False) :module: axionbloch.MilkyWayAxionHalo :staticmethod:

  Calculate the analytical axion power-spectral-density lineshape.

  Implements Eqs. (12), (19), (20) of the Gramolin lineshape paper for
  the non-gradient, parallel-gradient, and perpendicular-gradient coupling
  cases respectively.

  .. warning::
      ``nu`` should not be too far above ``nu_a``.  The ratio
      ``nu / nu_a`` must remain below ~1.03 to avoid numerical overflow in
      the ``sinh(beta)`` term.

  :param v_0: Local circular-rotation speed of the Milky Way.
  :type v_0: Quantity [m/s]
  :param v_lab: Speed of the laboratory in the galactic rest frame.
  :type v_lab: Quantity [m/s]
  :param nu_a: Axion Compton frequency.
  :type nu_a: Quantity [Hz]
  :param nu: Frequency array at which to evaluate the lineshape.  Must be
             uniformly spaced.
  :type nu: Quantity array [Hz]
  :param case: Coupling geometry: ``'non-grad'``, ``'grad_par'``, or
               ``'grad_perp'``.
  :type case: str
  :param alpha: Angle between the sensitive axis and the axion wind velocity.
  :type alpha: Quantity [rad]
  :param verbose: Print diagnostic information.
  :type verbose: bool

  :returns: Power spectral density lineshape, normalized so that
            ``integral(lineshape * dnu) = 1``.
  :rtype: Quantity array, shape ``(len(nu),)`` [Hz⁻²]

  .. rubric:: References

  A. Gramolin et al., https://github.com/gramolin/lineshape

.. py:attribute:: MilkyWayAxionHalo.g_aNN :module: axionbloch.MilkyWayAxionHalo :type: ~astropy.units.quantity.Annotated[~astropy.units.quantity.Quantity, Unit(“1 / GeV”)]

.. py:method:: MilkyWayAxionHalo.getAmpSpectra(frequencies: ~numpy.ndarray, case: str = ‘grad_perp’, numSpectra: int = 1, rand_seed: int = None, use_stoch: bool = True, verbose: bool = False) -> ~numpy.ndarray :module: axionbloch.MilkyWayAxionHalo

  Return complex amplitude spectra for the axion pseudo-magnetic field.

  Draws ``numSpectra`` realizations of the stochastic axion field in the
  frequency domain using the SHM lineshape as the power spectral density.
  Each realization has random phases uniformly distributed in ``[0, 2π)``.
  When ``use_stoch=True``, amplitudes are additionally drawn from an
  exponential distribution, matching the expected statistics of a
  narrowband stochastic signal.

  :param frequencies: Absolute frequencies at which to evaluate the spectrum.
  :type frequencies: Quantity array [Hz]
  :param case: Lineshape coupling case: ``'non-grad'``, ``'grad_par'``, or
               ``'grad_perp'``.
  :type case: str
  :param numSpectra: Number of independent field realizations to generate.
  :type numSpectra: int
  :param rand_seed: Seed for the random-number generator (for reproducibility).
  :type rand_seed: int, optional
  :param use_stoch: If ``True``, draw stochastic amplitudes; if ``False``, use the
                    deterministic ``sqrt(PSD)`` amplitude.
  :type use_stoch: bool
  :param verbose: Unused; reserved for future diagnostic output.
  :type verbose: bool

  :returns: Complex amplitude spectra.
  :rtype: ndarray, shape ``(numSpectra, len(frequencies))``

.. py:method:: MilkyWayAxionHalo.getRabiFreq(g_aNN: ~astropy.units.quantity.Annotated[~astropy.units.quantity.Quantity, Unit(“1 / GeV”)] | None = None, case=’grad_perp’, verbose=False) -> ~astropy.units.quantity.Quantity :module: axionbloch.MilkyWayAxionHalo

  get the Rabi frequency of the pseudomagnetic field amplitude in [Hz] for the specified case
  case: "non-grad", "grad_par" or "grad_perp", determines the lineshape function to use

.. py:attribute:: MilkyWayAxionHalo.m_a :module: axionbloch.MilkyWayAxionHalo :type: ~astropy.units.quantity.Annotated[~astropy.units.quantity.Quantity, Unit(“kg”)]

.. py:attribute:: MilkyWayAxionHalo.nu_a :module: axionbloch.MilkyWayAxionHalo :type: ~astropy.units.quantity.Annotated[~astropy.units.quantity.Quantity, Unit(“Hz”)]

.. py:attribute:: MilkyWayAxionHalo.nu_a_eff :module: axionbloch.MilkyWayAxionHalo :type: ~astropy.units.quantity.Annotated[~astropy.units.quantity.Quantity, Unit(“Hz”)]

.. py:attribute:: MilkyWayAxionHalo.rho_E_DM :module: axionbloch.MilkyWayAxionHalo :type: ~astropy.units.quantity.Annotated[~astropy.units.quantity.Quantity, Unit(“GeV / cm3”)] :value: <Quantity 0.3 GeV / cm3>

.. py:attribute:: MilkyWayAxionHalo.tau_a_est :module: axionbloch.MilkyWayAxionHalo :type: ~astropy.units.quantity.Annotated[~astropy.units.quantity.Quantity, Unit(“s”)]

.. py:attribute:: MilkyWayAxionHalo.v_0 :module: axionbloch.MilkyWayAxionHalo :type: ~astropy.units.quantity.Annotated[~astropy.units.quantity.Quantity, Unit(“km / s”)] :value: <Quantity 220. km / s>

.. py:attribute:: MilkyWayAxionHalo.v_lab :module: axionbloch.MilkyWayAxionHalo :type: ~astropy.units.quantity.Annotated[~astropy.units.quantity.Quantity, Unit(“km / s”)] :value: <Quantity 233. km / s>

axionbloch.GravBoundAxionHalo

Provides the general TISE solver GravBoundAxionHalo for axions gravitationally bound to a compact body. The solver:

  • Builds a 1-D radial finite-difference Hamiltonian \(H = T + V_\mathrm{eff}\).

  • Diagonalizes it with scipy.linalg.eigh for each angular-momentum channel l.

  • Stores eigen-energies and normalized reduced radial wavefunctions \(u_{n_r l}(r)\).

  • Computes the 3-D gradient of the total wavefunction toward a geographic Station via findGradientsAtDirection(), and its time evolution via findGradientsOverTime().

  • Computes the axion-nucleon coupling frequency (Omega_a) from gradients via findOmega_aOverTime(), with RMS value computation via findRmsOmega_aOverTime().

Sub-class EarthBoundAxionHalo is pre-configured with the Earth’s gravitational potential.

.. py:module:: axionbloch.GravBoundAxionHalo

Solver for gravitationally-bound axion halos around bodies.

:class:GravBoundAxionHalo solves the time-independent Schrödinger equation (TISE) on a 1-D radial grid using a finite-difference Hamiltonian, returning bound-state wavefunctions and energies for each orbital quantum number l.

:class:EarthBoundAxionHalo (in :mod:axionbloch.EarthBoundAxionHalo) is a subclass pre-configured with the Earth’s gravitational potential.

.. py:class:: GravBoundAxionHalo(name=’Gravitationally Bound Axion Halo’, nu_a: ~astropy.units.quantity.Quantity | None = None, N: int = 4096, extent: ~astropy.units.quantity.Quantity = <Quantity 128. earthRad>, getPot=None, mass_enclosed=None, g_aNN=None, a_0=None, verbose: bool = False) :module: axionbloch.GravBoundAxionHalo

Bases: :py:class:object

Solve the TISE for axions gravitationally bound to compact bodies.

Constructs a 1-D radial finite-difference Hamiltonian H = T + V on a uniform grid spanning [-extent/2, extent/2], diagonalizes it for each requested angular-momentum channel l, and stores the resulting wavefunctions and energy expectation values in :attr:states.

Note: input / output units are in SI, while internal calculations are in atomic units.

.. py:attribute:: GravBoundAxionHalo.a_0 :module: axionbloch.GravBoundAxionHalo :type: ~astropy.units.quantity.Quantity | None

.. py:method:: GravBoundAxionHalo.compareGradientsOverStates(station: ~axionbloch.Station.Station, meas_time: ~astropy.time.core.Time, stateNamesDict: dict, truncRadius: ~astropy.units.quantity.Quantity | None = None, showPlot: bool = True, verbose: bool = False) -> dict :module: axionbloch.GravBoundAxionHalo

  Compare gradient components across different eigenstate combinations.

  Calls :meth:`findGradientsOverStates` and extracts gradient values at
  the station's location. Plots the three gradient components as scatter
  points for each state combination.

  :param station: Geographic location defining the direction and evaluation radius.
  :type station: :class:`~axionbloch.Station.Station`
  :param meas_time: Measurement epoch.
  :type meas_time: :class:`astropy.time.Time`
  :param stateNamesDict: Dictionary mapping labels to state name lists. Example:
                         ``{"2p": ["2p"], "2p and 3p": ["2p", "3p"]}``.
  :type stateNamesDict: dict
  :param truncRadius: Radial truncation passed through to :meth:`findGradientsOverStates`.
  :type truncRadius: Quantity, optional
  :param showPlot: If True, display the comparison plot (default: True).
  :type showPlot: bool
  :param verbose: Print per-step progress.
  :type verbose: bool

  :returns: * *dict with keys*
            * \* ``'state_labels'`` — list of state combination labels
            * \* ``'grad_r'``       — Quantity array of grad_r values at station radius
            * \* ``'grad_theta'``   — Quantity array of grad_theta values at station radius
            * \* ``'grad_phi'``     — Quantity array of grad_phi values at station radius
            * \* ``'r_eval'``       — station radius used for evaluation

.. py:method:: GravBoundAxionHalo.findGradientsAtDirection(stateNames: list[str], station: ~axionbloch.Station.Station | None = None, meas_time: ~astropy.time.core.Time | None = None, truncRadius: ~astropy.units.quantity.Quantity | None = None, showPlot: bool = True, verbose: bool = False) :module: axionbloch.GravBoundAxionHalo

  Compute and plot the 3-D gradient of the total wavefunction at a direction.

  Superimposes the requested eigenstates (with equal weight), computes the
  spherical-coordinate gradient (∂_r, ∂_θ/r, ∂_φ/(r sinθ)), interpolates
  each component onto a radial line pointing toward the requested direction,
  and plots all four quantities (wavefunction + three gradient components).

  The direction should be specified by a :class:`~axionbloch.Station.Station`
  object.

  :param stateNames: Labels of eigenstates to include (e.g. ``['1s', '2p']``).
                     Defaults to the lowest-energy state.
  :type stateNames: list of str
  :param station: Geographic station whose latitude, longitude, and elevation define
                  the direction.
  :type station: Station, optional
  :param truncRadius: Truncation radius for reducing computation and plotting time.
  :type truncRadius: Quantity

.. py:method:: GravBoundAxionHalo.findGradientsOverStates(station: ~axionbloch.Station.Station, meas_time: ~astropy.time.core.Time, stateNamesDict: dict, truncRadius: ~astropy.units.quantity.Quantity | None = None, showPlot: bool = True, verbose: bool = False) -> dict :module: axionbloch.GravBoundAxionHalo

  Compute and plot gradient components for multiple state combinations.

  Calls :meth:`findGradientsAtDirection` for each state combination in
  ``stateNamesDict``, computes the three gradient components at a given
  direction, and optionally plots them in a combined figure.

  :param station: Geographic location defining the direction.
  :type station: :class:`~axionbloch.Station.Station`
  :param meas_time: Measurement epoch.
  :type meas_time: :class:`astropy.time.Time`
  :param stateNamesDict: Dictionary mapping labels to state name lists. Example:
                         ``{"2p": ["2p"], "2p and 3p": ["2p", "3p"]}``.
  :type stateNamesDict: dict
  :param truncRadius: Radial truncation passed through to :meth:`findGradientsAtDirection`.
  :type truncRadius: Quantity, optional
  :param showPlot: If True, display the plot of gradient components (default: True).
  :type showPlot: bool
  :param verbose: Print per-step progress.
  :type verbose: bool

  :returns: * *dict with keys*
            * \* ``'state_labels'`` — list of state combination labels
            * \* ``'grad_r'``     — dict of Quantity arrays for each state combination
            * \* ``'grad_theta'`` — dict of Quantity arrays for each state combination
            * \* ``'grad_phi'``   — dict of Quantity arrays for each state combination
            * \* ``'r_line'``     — common radial grid for all combinations

.. py:method:: GravBoundAxionHalo.findGradientsOverTime(stateNames: list[str], station: ~axionbloch.Station.Station, meas_times: list[~astropy.time.core.Time], truncRadius: ~astropy.units.quantity.Quantity | None = None, verbose: bool = False) -> dict :module: axionbloch.GravBoundAxionHalo

  Gradient components at a station evaluated over a list of epochs.

  Calls :meth:`findGradientsAtDirection` for each time and collects the
  three gradient values at Earth's surface (r = 1 R_earth).

  :param stateNames: Eigenstate labels to superimpose.
  :type stateNames: list of str
  :param station: Geographic location.
  :type station: :class:`~axionbloch.Station.Station`
  :param times: Epochs at which to evaluate the gradients.
  :type times: iterable of :class:`astropy.time.Time`
  :param truncRadius: Radial truncation passed through to :meth:`findGradientsAtDirection`.
  :type truncRadius: Quantity, optional
  :param verbose: Print per-step progress.
  :type verbose: bool

  :returns: * *dict with keys*
            * \* ``'times'``      — input time list
            * \* ``'grad_r'``     — Quantity array, shape ``(N_times,)``
            * \* ``'grad_theta'`` — Quantity array, shape ``(N_times,)``
            * \* ``'grad_phi'``   — Quantity array, shape ``(N_times,)``

.. py:method:: GravBoundAxionHalo.findHighProbStates(radius_range=[<Quantity 0.9 earthRad>, <Quantity 1.1 earthRad>], threshold=0.01) :module: axionbloch.GravBoundAxionHalo

  Print and plot eigenstates with significant probability in a radial shell.

  Iterates over all solved states, computes the integrated probability
  inside ``radius_range``, and reports the norm and partial integral for
  each state.

  :param radius_range: ``[r_min, r_max]`` defining the radial shell of interest.
  :type radius_range: list of Quantity
  :param threshold: Minimum probability fraction to flag a state (currently unused,
                    reserved for future filtering).
  :type threshold: float

.. py:method:: GravBoundAxionHalo.findOmega_aOverTime(stateNames: list[str], station: ~axionbloch.Station.Station, meas_times, truncRadius: ~astropy.units.quantity.Quantity | None = None, g_aNN: ~astropy.units.quantity.Annotated[~astropy.units.quantity.Quantity, Unit(“1 / GeV”)] = <Quantity 1.e-09 1 / GeV>, verbose: bool = False) -> dict :module: axionbloch.GravBoundAxionHalo

  Axion-nucleon coupling frequency (Omega_a) evaluated over a list of epochs.

  Calls :meth:`findGradientsOverTime` and computes Omega_a = g_aNN * gradient
  (with gradient in units of m^-1, ignoring radian units).

  :param stateNames: Eigenstate labels to superimpose.
  :type stateNames: list of str
  :param station: Geographic location.
  :type station: :class:`~axionbloch.Station.Station`
  :param meas_times: Epochs at which to evaluate Omega_a.
  :type meas_times: iterable of :class:`astropy.time.Time`
  :param truncRadius: Radial truncation passed through to :meth:`findGradientsOverTime`.
  :type truncRadius: Quantity, optional
  :param verbose: Print per-step progress.
  :type verbose: bool

  :returns: * *dict with keys*
            * \* ``'times'``      — input time list
            * \* ``'Omega_a_r'``     — Quantity array, shape ``(N_times,)``, Omega_a from radial gradient
            * \* ``'Omega_a_theta'`` — Quantity array, shape ``(N_times,)``, Omega_a from theta gradient
            * \* ``'Omega_a_phi'``   — Quantity array, shape ``(N_times,)``, Omega_a from phi gradient

.. py:method:: GravBoundAxionHalo.findrmsOmega_aOverTime(stateNames: list[str], station: ~axionbloch.Station.Station, meas_times, truncRadius: ~astropy.units.quantity.Quantity | None = None, g_aNN: ~astropy.units.quantity.Annotated[~astropy.units.quantity.Quantity, Unit(“1 / GeV”)] = <Quantity 1.e-09 1 / GeV>, verbose: bool = False) -> dict :module: axionbloch.GravBoundAxionHalo

  RMS Omega_a over a list of epochs.

  Calls :meth:`findOmega_aOverTime` and computes the root-mean-square
  (RMS) value for each gradient component.

  :param stateNames: Eigenstate labels to superimpose.
  :type stateNames: list of str
  :param station: Geographic location.
  :type station: :class:`~axionbloch.Station.Station`
  :param meas_times: Epochs at which to evaluate Omega_a.
  :type meas_times: iterable of :class:`astropy.time.Time`
  :param truncRadius: Radial truncation passed through to :meth:`findOmega_aOverTime`.
  :type truncRadius: Quantity, optional
  :param g_aNN: Axion-nucleon coupling constant (default: 1e-9 GeV⁻¹).
  :type g_aNN: Quantity, optional
  :param verbose: Print per-step progress.
  :type verbose: bool

  :returns: * *dict with keys*
            * \* ``'rms_Omega_a_r'``     — RMS of Omega_a from radial gradient
            * \* ``'rms_Omega_a_theta'`` — RMS of Omega_a from theta gradient
            * \* ``'rms_Omega_a_phi'``   — RMS of Omega_a from phi gradient

.. py:attribute:: GravBoundAxionHalo.g_aNN :module: axionbloch.GravBoundAxionHalo :type: ~astropy.units.quantity.Quantity | None

.. py:method:: GravBoundAxionHalo.getStateEnergies() :module: axionbloch.GravBoundAxionHalo

  Return a list of eigen-energies in the units stored in ``self.states``.

.. py:method:: GravBoundAxionHalo.getStateNames() :module: axionbloch.GravBoundAxionHalo

  Return a list of state label strings (e.g. ``['1s', '2s', '2p']``).

.. py:method:: GravBoundAxionHalo.listEigenStates() :module: axionbloch.GravBoundAxionHalo

  Print a formatted table of all solved eigenstates.

  Columns: radial quantum number, *l*, principal quantum number, label,
  eigen-energy (eV), kinetic energy (eV), and mean axion speed (m/s).
  States are listed in ascending eigen-energy order.

.. py:attribute:: GravBoundAxionHalo.mass_enclosed :module: axionbloch.GravBoundAxionHalo :type: ~astropy.units.quantity.Quantity | None

.. py:attribute:: GravBoundAxionHalo.orbitalLabels :module: axionbloch.GravBoundAxionHalo :value: [‘s’, ‘p’, ‘d’, ‘f’, ‘g’, ‘h’, ‘i’, ‘k’, ‘l’, ‘m’]

.. py:method:: GravBoundAxionHalo.plotEigenEnergiesInPot() :module: axionbloch.GravBoundAxionHalo

  Plot eigen-energies superimposed on the gravitational potential.

  Draws a horizontal line for each bound state at its eigen-energy,
  spanning the classical turning points where the potential crosses that
  energy level.

.. py:method:: GravBoundAxionHalo.plotEigenStates(numStates: int = 8, startState: int = 0, truncRadius: ~astropy.units.quantity.Quantity | None = None, xlim=(-0.3, 5.3), ylim=None, showPlot=True) :module: axionbloch.GravBoundAxionHalo

  Plot a vertical stack of reduced radial wavefunctions.

  Eigenstates are plotted from lowest energy (bottom) to highest (top),
  sharing a common y-scale so amplitudes can be compared visually.

  :param numStates: How many consecutive states (starting from ``startState``) to show.
  :type numStates: int
  :param startState: Index into the energy-sorted state list at which to begin.
  :type startState: int
  :param xlim: x-axis limits in units of Earth radii.
  :type xlim: tuple or None
  :param ylim: Shared y-axis limits.  Auto-computed from peak amplitude if None.
  :type ylim: tuple or None

.. py:method:: GravBoundAxionHalo.plotEigenstate(n_r: int, l: int) :module: axionbloch.GravBoundAxionHalo

  Plot the reduced radial wavefunction for a single eigenstate.

  If the requested state has not been solved yet, calls
  :meth:`solve_TISE_3D_l` automatically.

  :param n_r: Radial quantum number (number of radial nodes).
  :type n_r: int
  :param l: Orbital angular-momentum quantum number.
  :type l: int

.. py:method:: GravBoundAxionHalo.plotGradients(stateNames: str, r, R_r, r_line, grad_r_line, grad_theta_line, grad_phi_line, station: ~axionbloch.Station.Station | None = None, label: str | None = None) :module: axionbloch.GravBoundAxionHalo

  Plotting helper for :meth:`findGradients`.

  :param station: If provided, ``station.name`` is used as the plot title.
  :type station: Station, optional
  :param label: Fallback title when ``station`` is ``None``.
  :type label: str, optional

.. py:attribute:: GravBoundAxionHalo.pot :module: axionbloch.GravBoundAxionHalo :type: ~astropy.units.quantity.Quantity

.. py:method:: GravBoundAxionHalo.showValueAndUnits() :module: axionbloch.GravBoundAxionHalo

  Print values and units of the key physical quantities.

  Ideally the numerical values should be close to 1 and units should be
  identical for quantities compared or added together, making unit
  mismatches easy to spot during development.

.. py:method:: GravBoundAxionHalo.solve_TISE_3D(l_vals=[3], max_n_r: int = 64, showPlot: bool = False, verbose: bool = False) :module: axionbloch.GravBoundAxionHalo

  Solve the TISE for all requested angular-momentum channels.

  Iterates over ``l_vals``, calling :meth:`solve_TISE_3D_l` for each,
  then sorts all accumulated states by eigen-energy.

  :param l_vals: Orbital angular-momentum quantum numbers to solve.
  :type l_vals: list of int
  :param max_n_r: Number of radial eigenstates to retain per channel.
  :type max_n_r: int
  :param showPlot: Forwarded to :meth:`solve_TISE_3D_l`.
  :type showPlot: bool
  :param verbose: Forwarded to :meth:`solve_TISE_3D_l`.
  :type verbose: bool

.. py:method:: GravBoundAxionHalo.solve_TISE_3D_l(l: int, showPlot: bool = False, max_n_r: int = 10, verbose: bool = False) :module: axionbloch.GravBoundAxionHalo

  Solve the TISE for a single angular-momentum channel *l*.

  Builds the finite-difference Hamiltonian ``H = T + V_eff`` (with the
  centrifugal barrier included in ``V_eff``), diagonalizes it with
  ``scipy.linalg.eigh``, normalizes the lowest ``max_n_r`` eigenstates,
  computes expectation values of T, V, and V_eff, and stores the results
  in :attr:`states`.

  :param l: Orbital angular-momentum quantum number.
  :type l: int
  :param showPlot: If True, plot the reduced radial wavefunctions after solving.
  :type showPlot: bool
  :param max_n_r: Number of radial eigenstates to retain (counted from the ground
                  state of channel *l*).
  :type max_n_r: int
  :param verbose: Print timing and eigen-energy tables.
  :type verbose: bool

.. py:method:: GravBoundAxionHalo.sortByEigenE() :module: axionbloch.GravBoundAxionHalo

  Sort ``self.states`` in-place by ascending eigen-energy.

.. py:method:: GravBoundAxionHalo.stackEigenStates(numStates: int = 8, startState: int = 0, xlim=(-0.3, 5.3), ylim=None) :module: axionbloch.GravBoundAxionHalo

  Plot a vertical stack of reduced radial wavefunctions.

  Eigenstates are plotted from lowest energy (bottom) to highest (top),
  sharing a common y-scale so amplitudes can be compared visually.

  :param numStates: How many consecutive states (starting from ``startState``) to show.
  :type numStates: int
  :param startState: Index into the energy-sorted state list at which to begin.
  :type startState: int
  :param xlim: x-axis limits in units of Earth radii.
  :type xlim: tuple or None
  :param ylim: Shared y-axis limits.  Auto-computed from peak amplitude if None.
  :type ylim: tuple or None

axionbloch.EarthBoundAxionHalo

Sub-class of GravBoundAxionHalo pre-configured with Earth’s gravitational potential from the Preliminary Earth Model (PEM) density profile. Adds Earth-specific gradient helpers: findGradients() (takes a Station), findGradientsAtEarthSurface(), findGradientsWithMilkyWay(), and the pseudomagnetic-field conversion getBfield(). Also exports helper functions used to build the potential:

Function

Description

loadPEMdata

Load PEM density-profile data

PREM_density

PREM piecewise density at a given radius

earth_grav_potential_earth_center

Gravitational potential referenced to Earth’s centre

earth_grav_potential_infty

Gravitational potential referenced to infinity

get_CumulativeMass

Cumulative enclosed mass as a function of radius

plot_earth_grav_potential

Plot density, mass, and potential profiles

.. py:module:: axionbloch.EarthBoundAxionHalo

.. py:class:: EarthBoundAxionHalo(name=’Earth-Bound Axion Halo’, nu_a: ~astropy.units.quantity.Annotated[~astropy.units.quantity.Quantity, Unit(“Hz”)] | None = None, N: int = 4096, extent: ~astropy.units.quantity.Annotated[~astropy.units.quantity.Quantity, Unit(“m”)] = <Quantity 128. earthRad>, getPot=, a_0: ~astropy.units.quantity.Annotated[~astropy.units.quantity.Quantity, Unit(“eV”)] | None = None, totalMassEnclosed: ~astropy.units.quantity.Annotated[~astropy.units.quantity.Quantity, Unit(“kg”)] | None = <Quantity 4.e-09 earthMass>, g_aNN: ~astropy.units.quantity.Annotated[~astropy.units.quantity.Quantity, Unit(“1 / GeV”)] = <Quantity 1.e-09 1 / GeV>, verbose: bool = False) :module: axionbloch.EarthBoundAxionHalo

Bases: :py:class:~axionbloch.GravBoundAxionHalo.GravBoundAxionHalo

.. py:attribute:: EarthBoundAxionHalo.N :module: axionbloch.EarthBoundAxionHalo :type: int :value: 4096

.. py:attribute:: EarthBoundAxionHalo.a_0 :module: axionbloch.EarthBoundAxionHalo :type: ~astropy.units.quantity.Annotated[~astropy.units.quantity.Quantity, Unit(“eV”)] | None :value: None

.. py:method:: EarthBoundAxionHalo.coh_time_g1() :module: axionbloch.EarthBoundAxionHalo

  x : complex-valued time series
  dt: sampling interval
  method: "1e" or "integral"

.. py:attribute:: EarthBoundAxionHalo.extent :module: axionbloch.EarthBoundAxionHalo :type: ~astropy.units.quantity.Annotated[~astropy.units.quantity.Quantity, Unit(“m”)] :value: <Quantity 128. earthRad>

.. py:method:: EarthBoundAxionHalo.findGradients(stateNames: list[str] | None = None, station: ~axionbloch.Station.Station | None = None, truncRadius: ~astropy.units.quantity.Quantity | None = <Quantity 3. earthRad>, showPlot: bool = False, verbose: bool = False) -> tuple :module: axionbloch.EarthBoundAxionHalo

  Compute the wavefunction gradient toward a geographic direction.

  The direction is specified by a :class:`~axionbloch.Station.Station`;
  raises :exc:`ValueError` when none is provided.

  :param stateNames:
  :type stateNames: list of str, optional
  :param station:
  :type station: Station, optional
  :param truncRadius:
  :type truncRadius: Quantity, optional
  :param showPlot:
  :type showPlot: bool
  :param verbose:
  :type verbose: bool

  :rtype: r, R_r, r_line, grad_r, grad_theta, grad_phi

.. py:method:: EarthBoundAxionHalo.findGradientsAtEarthSurface(station: ~axionbloch.Station.Station, stateNames: list[str] | None = None, truncRadius: ~astropy.units.quantity.Quantity | None = <Quantity 3. earthRad>, showPlot: bool = False, verbose: bool = False) -> tuple :module: axionbloch.EarthBoundAxionHalo

  Compute the wavefunction gradient toward a station at the Earth's surface.

  A convenience wrapper around :meth:`findGradientsAtDirection` with a
  default truncation radius suited to ground-based experiments.

  :param station: Geographic station whose latitude, longitude, and elevation define
                  the direction.
  :type station: Station
  :param stateNames: Eigenstates to superimpose (e.g. ``['2p']``).  Defaults to the
                     lowest-energy solved state.
  :type stateNames: list of str, optional
  :param truncRadius: Radial cutoff for the interpolation grid.
  :type truncRadius: Quantity, optional
  :param showPlot: Plot the wavefunction and gradient profiles.
  :type showPlot: bool
  :param verbose: Print timing and diagnostic information.
  :type verbose: bool

  :returns: Same six arrays as :meth:`findGradients`.
  :rtype: r, R_r, r_line, grad_r, grad_theta, grad_phi

.. py:method:: EarthBoundAxionHalo.findGradientsWithMilkyWay(mw, stateNames: list[str] | None = None, truncRadius: ~astropy.units.quantity.Quantity | None = None, showPlot: bool = False, verbose: bool = False) -> dict :module: axionbloch.EarthBoundAxionHalo

  Compute the gradient at a station with time-dependent galactic context.

  Uses the station embedded in the :class:`~axionbloch.MilkyWay.MilkyWay`
  instance for the geographic direction, then enriches the result with
  galactic kinematics derived from that same object:

  - Lab velocity :math:`\mathbf{v}_\mathrm{lab}` (magnitude and direction).
  - Wind angle between :math:`\mathbf{v}_\mathrm{lab}` and the
    station's sensitive axis.
  - Gradient in Cartesian ITRS coordinates (useful for projecting onto
    a non-vertical :math:`\mathbf{B}_0`).
  - Projection of the gradient onto the sensitive axis (radial direction).

  :param mw: Must have :attr:`~axionbloch.MilkyWay.MilkyWay.station` set.
  :type mw: :class:`~axionbloch.MilkyWay.MilkyWay`
  :param stateNames: Eigenstates to include (e.g. ``['2p']``).
  :type stateNames: list of str, optional
  :param truncRadius: Radial cutoff.
  :type truncRadius: Quantity, optional
  :param showPlot: Plot the wavefunction and gradient profiles.
  :type showPlot: bool
  :param verbose: Print diagnostics.
  :type verbose: bool

  :returns: ``r``, ``R_r``, ``r_line`` — radial grid and wavefunction.

            ``grad_r``, ``grad_theta``, ``grad_phi`` — spherical gradient
            components along the full r_line.

            ``grad_r_surface``, ``grad_theta_surface``, ``grad_phi_surface`` —
            values at Earth's surface (:math:`r = R_\oplus`).

            ``grad_cartesian_itrs`` — Cartesian gradient vector in ITRS
            (x = prime meridian, z = north pole) at Earth's surface.

            ``nvec_gcrs`` — station's unit normal in GCRS (equatorial inertial),
            changes with time as Earth rotates.

            ``v_lab``, ``v_lab_magnitude`` — lab velocity in galactic frame.

            ``wind_angle`` — angle between v_lab and the sensitive axis [rad].
  :rtype: dict with keys

  .. rubric:: Examples

  >>> from astropy.time import Time
  >>> from axionbloch.MilkyWay import MilkyWay
  >>> from axionbloch.Station import Baltimore
  >>> mw = MilkyWay(time=Time('2024-06-21T14:00:00'), station=Baltimore)
  >>> result = halo.findGradientsWithMilkyWay(
  ...     mw, stateNames=['2p'], truncRadius=2 * unit.R_earth
  ... )
  >>> print(result['wind_angle'].to('deg'))
  >>> print(result['grad_r_surface'])

.. py:attribute:: EarthBoundAxionHalo.g_aNN :module: axionbloch.EarthBoundAxionHalo :type: ~astropy.units.quantity.Annotated[~astropy.units.quantity.Quantity, Unit(“1 / GeV”)] :value: <Quantity 1.e-09 1 / GeV>

.. py:method:: EarthBoundAxionHalo.getBfield(rate_Hz: float, timeLen: int, rand_seed: int, numFields: int = 1, verbose: bool = False) :module: axionbloch.EarthBoundAxionHalo

.. py:attribute:: EarthBoundAxionHalo.nu_a :module: axionbloch.EarthBoundAxionHalo :type: ~astropy.units.quantity.Annotated[~astropy.units.quantity.Quantity, Unit(“Hz”)] | None :value: None

.. py:method:: EarthBoundAxionHalo.plotEigenStates(numStates: int = 8, startState: int = 0, truncRadius: ~astropy.units.quantity.Quantity | None = <Quantity 3. earthRad>, xlim=(-0.3, 5.3), ylim=None, showPlot=True, savefig=False) :module: axionbloch.EarthBoundAxionHalo

  Plot a vertical stack of reduced radial wavefunctions.

  Eigenstates are plotted from lowest energy (bottom) to highest (top),
  sharing a common y-scale so amplitudes can be compared visually.

  :param numStates: How many consecutive states (starting from ``startState``) to show.
  :type numStates: int
  :param startState: Index into the energy-sorted state list at which to begin.
  :type startState: int
  :param xlim: x-axis limits in units of Earth radii.
  :type xlim: tuple or None
  :param ylim: Shared y-axis limits.  Auto-computed from peak amplitude if None.
  :type ylim: tuple or None

.. py:attribute:: EarthBoundAxionHalo.totalMassEnclosed :module: axionbloch.EarthBoundAxionHalo :type: ~astropy.units.quantity.Annotated[~astropy.units.quantity.Quantity, Unit(“kg”)] | None :value: <Quantity 4.e-09 earthMass>

.. py:function:: PREM_density(radius_km) :module: axionbloch.EarthBoundAxionHalo

.. py:function:: earth_grav_potential_earth_center() :module: axionbloch.EarthBoundAxionHalo

Returns a function Phi(r), valid both inside and outside Earth. Uses PREM-like model for interior, point-mass approximation for exterior.

.. py:function:: earth_grav_potential_infty() :module: axionbloch.EarthBoundAxionHalo

Returns a function Phi(r[m]) [J/kg], valid both inside and outside Earth. Uses PREM-like model for interior, point-mass approximation for exterior.

.. py:function:: get_CumulativeMass() :module: axionbloch.EarthBoundAxionHalo

Returns the cumulative mass as a function of radius. Uses PREM-like model for interior, point-mass approximation for exterior.

.. py:function:: loadPEMdata(filepath=’axionbloch/data/Earth_Models/PREM_data.txt’) :module: axionbloch.EarthBoundAxionHalo

.. py:function:: plot_earth_grav_potential(showplot=True) :module: axionbloch.EarthBoundAxionHalo

axionbloch.Station

Defines Station, which converts geographic coordinates (latitude, longitude, elevation) to spherical coordinates \((\theta, \phi)\) and the outward unit normal \(\hat{n}\). in_solarZ_frame() expresses the station position in a solar-oriented frame at a given measurement time.

Pre-built station instances exported from this module:

Mainz, Geneva, Baltimore, Tokyo, Mumbai, Sanya, Sydney, CapeTown, BuenosAires.

.. py:module:: axionbloch.Station

Geographic station definitions for axion-wind calculations.

A :class:Station stores the location of an experimental site on Earth and exposes the corresponding unit vector and distance from the Earth’s centre.

.. py:class:: Station(name: str, *, NSsemisphere: str, EWsemisphere: str, latitude: ~astropy.units.quantity.Quantity, longitude: ~astropy.units.quantity.Quantity, elevation: ~astropy.units.quantity.Quantity, verbose: bool = False) :module: axionbloch.Station

Bases: :py:class:object

A geographic location on Earth’s surface.

Converts latitude / longitude / elevation to spherical coordinates (theta, phi) with the following convention:

  • theta = 0 at the north pole, increases towards the south pole.

  • phi = 0 on the prime meridian (Greenwich), increases eastward.

  • nvec is the outward unit normal at the station’s location.

  • R is the distance from Earth’s centre to the station.

.. attribute:: name

  :type: str

.. attribute:: NSsemisphere

  ``'N'`` or ``'S'``.

  :type: str

.. attribute:: EWsemisphere

  ``'E'`` or ``'W'``.

  :type: str

.. attribute:: latitude, longitude

  Absolute values; hemisphere is tracked via ``NSsemisphere`` /
  ``EWsemisphere``.

  :type: Quantity

.. attribute:: elevation

  Height above sea level.

  :type: Quantity

.. attribute:: theta, phi

  Spherical coordinates (rad).

  :type: Quantity

.. attribute:: nvec

  Outward unit normal vector in Cartesian coordinates.

  :type: ndarray, shape (3,)

.. attribute:: R

  Distance from Earth's centre.

  :type: Quantity

.. py:method:: Station.direction_to_spherical(direction) -> tuple[~astropy.units.quantity.Quantity, ~astropy.units.quantity.Quantity] :module: axionbloch.Station :staticmethod:

  Spherical coordinates (theta, phi) of a Cartesian direction in the Earth frame.

  Uses the same convention as :class:`Station`:
  theta = 0 at north pole, phi = 0 on prime meridian increasing eastward.

  :param direction: Direction vector ``(x, y, z)`` in ECEF Cartesian coordinates.
                    Need not be normalized.
  :type direction: array-like, shape (3,)

  :returns: * **theta** (*Quantity [rad]*) -- Polar angle from north pole in ``[0, π]``.
            * **phi** (*Quantity [rad]*) -- Azimuthal angle in ``(-π, π]``.

.. py:method:: Station.in_solarZ_frame(meas_time, output: str = ‘spherical’) :module: axionbloch.Station

  Station position expressed in the Sun–Earth-velocity frame.

  Builds a right-handed orthonormal frame at the given epoch:

  * **ẑ** — Earth → Sun direction.
  * **x̂** — Earth's heliocentric velocity, orthogonalized against ẑ.
  * **ŷ = ẑ × x̂** — completes the right-hand triad.

  The station's geocentric Cartesian position (GCRS) is projected onto
  this triad.

  :param time: Observation epoch.
  :type time: :class:`astropy.time.Time`
  :param output: ``'cartesian'`` → ``(x, y, z)`` each a Quantity [m].
                 ``'spherical'`` → ``(r, theta, phi)`` where *theta* ∈ [0, π] is
                 the polar angle from ẑ and *phi* ∈ (−π, π] is the azimuthal
                 angle from x̂.
  :type output: str

  :rtype: tuple of Quantity

axionbloch.Sample

Defines Sample, which encapsulates all NMR sample properties needed for a simulation:

  • Gyromagnetic ratio \(\gamma\), magnetic dipole moment \(\mu\).

  • Mass density, molar mass, spin count per molecule → spin number density.

  • Relaxation times \(T_1\), \(T_2\).

  • Volume, temperature, and initial polarization.

Key methods:

Method

Description

getThermalPol

Exact thermal polarization at a given field and temperature

getM0

Magnetisation \(M_0\) from a given polarization

getM0_SPN

Magnetisation from polarization and spin number density

getM0eqb

Equilibrium magnetization at a given field and temperature

.. py:module:: axionbloch.Sample

Physical sample descriptor for NMR/spin-precession experiments.

:class:Sample holds the material and geometric properties of a single-phase spin sample and provides helpers for computing spin density, thermal polarization, and equilibrium magnetization.

.. py:class:: Sample(name: str | None = None, gamma: ~astropy.units.quantity.Annotated[~astropy.units.quantity.Quantity, Unit(“Hz rad / T”)] | None = None, massDensity: ~astropy.units.quantity.Quantity | None = None, molarMass: ~astropy.units.quantity.Quantity | None = None, numOfSpinsPerMolecule: ~astropy.units.quantity.Quantity | None = None, T2: ~astropy.units.quantity.Quantity | None = None, T1: ~astropy.units.quantity.Quantity | None = None, vol: ~astropy.units.quantity.Quantity | None = None, mu: ~astropy.units.quantity.Quantity | None = None, temp: ~astropy.units.quantity.Quantity | None = None, pol: ~astropy.units.quantity.Quantity | None = None, verbose: bool = False) :module: axionbloch.Sample

Bases: :py:class:object

Single-phase spin sample for NMR/spin-precession experiments.

Stores material properties and geometry, and derives spin number density and total spin count on construction. Density inputs should be evaluated at STP (273.15 K, 10⁵ Pa per IUPAC 1982).

.. py:method:: Sample.getM0(pol: ~astropy.units.quantity.Quantity | None = None, verbose: bool = False) :module: axionbloch.Sample

  Return the bulk magnetization M0 = μ × pol × n_s [A/m].

  :param pol: Spin polarization; falls back to ``self.pol`` when omitted.
  :type pol: Quantity [dimensionless], optional

  :returns: Magnetization density.
  :rtype: Quantity [A/m]

.. py:method:: Sample.getM0_SPN(verbose: bool = False) :module: axionbloch.Sample

  Return the spin-projection-noise magnetization M0_SPN = μ × √N / V [A/m].

  Spin projection noise arises from quantum fluctuations: for N uncorrelated
  spin-1/2 particles the rms fluctuation in spin number is √N, giving an
  effective magnetization density μ√N / V.

  :returns: Spin-projection-noise magnetization density.
  :rtype: Quantity [A/m]

.. py:method:: Sample.getM0eqb(B_pol: ~astropy.units.quantity.Quantity, temp: ~astropy.units.quantity.Quantity | None = None, verbose: bool = False) :module: axionbloch.Sample

  Return the thermal-equilibrium magnetization at given B_pol and temperature.

  Convenience wrapper: calls :meth:`getThermalPol` then :meth:`getM0`.

  :param B_pol: Polarizing magnetic field strength.
  :type B_pol: Quantity [T]
  :param temp: Sample temperature; falls back to ``self.temp`` when omitted.
  :type temp: Quantity [K], optional

  :returns: Equilibrium magnetization density.
  :rtype: Quantity [A/m]

.. py:method:: Sample.getThermalPol(B_pol: ~astropy.units.quantity.Quantity, temp: ~astropy.units.quantity.Quantity | None = None, verbose: bool = False) -> ~astropy.units.quantity.Quantity :module: axionbloch.Sample

  Return the thermal spin polarization at a given field and temperature.

  Uses the exact result ``tanh(ℏγB / 2k_BT)``.

  :param B_pol: Polarizing magnetic field strength.
  :type B_pol: Quantity [T]
  :param temp: Sample temperature; falls back to ``self.temp`` when omitted.
  :type temp: Quantity [K], optional

  :returns: Thermal polarization in [0, 1].
  :rtype: Quantity [dimensionless]

axionbloch.Apparatus

Defines Magnet, which models the static DC polarising magnet. Field inhomogeneity is characterised by a Lorentzian distribution with a user-specified fractional FWHM (e.g. 2 ppm), discretised into numPt spin packets. The resulting spread of Larmor frequencies gives the effective \(T_2^*\) dephasing in simulations.

.. py:module:: axionbloch.Apparatus

.. py:class:: Magnet(name=’magnet’, B0: ~astropy.units.quantity.Annotated[~astropy.units.quantity.Quantity, Unit(“T”)] | None = None, direction: list | None = None, FWHM: ~astropy.units.quantity.Annotated[~astropy.units.quantity.Quantity, Unit(“ppm”)] | None = None, numPt: float = 1, nFWHM: float = 10.0, verbose: bool = False) :module: axionbloch.Apparatus

Bases: :py:class:object

DC magnetic field apparatus.

Models a static polarising magnet whose field may be spatially inhomogeneous. Inhomogeneity is characterised by a Lorentzian lineshape with a given FWHM, and the field distribution is discretised into numPt spin packets for simulation.

.. attribute:: name

  :type: str

.. attribute:: B0

  Nominal field strength.

  :type: astropy Quantity

.. attribute:: direction

  Unit vector along B0 (laboratory frame).

  :type: list or None

.. attribute:: FWHM

  Fractional field inhomogeneity (dimensionless, e.g. ppm).

  :type: astropy Quantity

.. attribute:: B_spread

  Sampled B0 values for each spin packet.

  :type: ndarray

.. attribute:: ratios

  Fractional weight of each spin packet (sums to 1).

  :type: ndarray

.. py:attribute:: Magnet.B0 :module: axionbloch.Apparatus :type: ~astropy.units.quantity.Quantity | None :value: None

.. py:attribute:: Magnet.FWHM :module: axionbloch.Apparatus :type: ~astropy.units.quantity.Quantity | None :value: None

.. py:attribute:: Magnet.direction :module: axionbloch.Apparatus :type: list | None :value: None

.. py:attribute:: Magnet.nFWHM :module: axionbloch.Apparatus :type: float :value: 10.0

.. py:attribute:: Magnet.name :module: axionbloch.Apparatus :type: str

.. py:attribute:: Magnet.numPt :module: axionbloch.Apparatus :type: float :value: 1

.. py:method:: Magnet.setHomogeneity(numPt: int | float = None, showPlot: bool = False, verbose: bool = False) :module: axionbloch.Apparatus

  Set the spin-packet sampling of the field inhomogeneity.

  Discretises a Lorentzian field distribution into ``numPt`` spin packets
  using a non-uniform sampling scheme (density proportional to
  ``|x|^(1/2)`` in normalized coordinates) to oversample the peak region.
  Each packet is assigned a weight ``ratio`` equal to its fraction of the
  total Lorentzian area so that weighted averages reproduce the continuous
  lineshape.

  For ``numPt == 1`` or zero FWHM the field is treated as homogeneous and
  a single packet at ``B0`` with weight 1 is used.

  :param numPt: Override the stored ``self.numPt``.
  :type numPt: int or float, optional
  :param showPlot: Display diagnostic histograms of the sampled B values and weights.
  :type showPlot: bool
  :param verbose: Print diagnostic information.
  :type verbose: bool

axionbloch.SimuTypes

Shared type definitions for the simulation layer.

  • SimuParams — TypedDict describing the full parameter set for one simulation run. Passed as all_params=[params] to Simulations.

  • SimuEntry — dataclass pairing a completed Simulation with its SimuParams. Stored in pool after a run.

.. py:module:: axionbloch.SimuTypes

Type definitions shared across the simulation layer.

:class:SimuParams (TypedDict) describes the full set of parameters passed to a Bloch-equation simulation run. :class:SimuEntry (dataclass) pairs a live :class:~axionbloch.SimuTools.Simulation object with its :class:SimuParams for bookkeeping across multiple runs.

.. py:class:: SimuEntry(simu: Simulation, params: SimuParams) :module: axionbloch.SimuTypes

Bases: :py:class:object

Pair of a completed :class:~axionbloch.SimuTools.Simulation and its parameters.

Stored in :attr:~axionbloch.SimuTools.Simulations.pool after a run.

.. attribute:: simu

  The executed simulation instance with result trajectories.

  :type: Simulation

.. attribute:: params

  The parameter dictionary used to configure and run ``simu``.

  :type: SimuParams

.. py:attribute:: SimuEntry.params :module: axionbloch.SimuTypes :type: SimuParams

.. py:attribute:: SimuEntry.simu :module: axionbloch.SimuTypes :type: Simulation

.. py:class:: SimuParams :module: axionbloch.SimuTypes

Bases: :py:class:~typing.TypedDict

Typed dictionary describing all parameters for a single simulation run.

Pass an instance of this dict to :class:~axionbloch.SimuTools.Simulations as one element of all_params.

.. attribute:: key_info

  Arbitrary metadata (e.g. ``{"nu_a": axion.nu_a}``), printed during
  verbose runs and stored alongside results.

  :type: dict

.. attribute:: axion

  Axion field model object.

  :type: MilkyWayAxionHalo

.. attribute:: sample

  NMR sample object.

  :type: Sample

.. attribute:: magnet

  Static bias-field object.

  :type: Magnet

.. attribute:: excField

  Excitation / pseudomagnetic field object.

  :type: MagField

.. attribute:: B_a_rms

  RMS axion-induced pseudomagnetic field amplitude (T).

  :type: Quantity, optional

.. attribute:: numFields

  Number of independent stochastic field realizations.

  :type: int

.. attribute:: rand_seed

  Random seed for reproducible field realizations.

  :type: int

.. attribute:: init_M

  Initial magnetization magnitude (dimensionless, normalized by M0).

  :type: Quantity

.. attribute:: init_M_theta

  Initial polar angle of the magnetization vector.

  :type: Quantity [rad]

.. attribute:: init_M_phi

  Initial azimuthal angle of the magnetization vector.

  :type: Quantity [rad]

.. attribute:: rate

  Output sampling rate of the simulation.

  :type: Quantity [Hz]

.. attribute:: duration

  Total duration of the simulation.

  :type: Quantity [s]

.. py:attribute:: SimuParams.B_a_rms :module: axionbloch.SimuTypes :type: Quantity | None

.. py:attribute:: SimuParams.axion :module: axionbloch.SimuTypes :type: MilkyWayAxionHalo

.. py:attribute:: SimuParams.duration :module: axionbloch.SimuTypes :type: Quantity | None

.. py:attribute:: SimuParams.excField :module: axionbloch.SimuTypes :type: MagField

.. py:attribute:: SimuParams.init_M :module: axionbloch.SimuTypes :type: Quantity | None

.. py:attribute:: SimuParams.init_M_phi :module: axionbloch.SimuTypes :type: Quantity

.. py:attribute:: SimuParams.init_M_theta :module: axionbloch.SimuTypes :type: Quantity

.. py:attribute:: SimuParams.key_info :module: axionbloch.SimuTypes :type: object | None :value: {}

.. py:attribute:: SimuParams.magnet :module: axionbloch.SimuTypes :type: Magnet

.. py:attribute:: SimuParams.numFields :module: axionbloch.SimuTypes :type: int

.. py:attribute:: SimuParams.rand_seed :module: axionbloch.SimuTypes :type: int

.. py:attribute:: SimuParams.rate :module: axionbloch.SimuTypes :type: Quantity | None

.. py:attribute:: SimuParams.sample :module: axionbloch.SimuTypes :type: Sample

axionbloch.SimuTools

The main simulation engine. Three classes work together:

MagField builds the time-domain effective magnetic field array. Call one of its set* methods to fill B_vec:

Method

Field type

setXYPulse

Continuous circularly-polarised excitation

set90DegPulse

Calibrated 90° hard pulse

setCPMGPulseTrain

CPMG spin-echo pulse train

setNPulsesArbDelay

N pulses with arbitrary delays

setAxionFields

Stochastic axion pseudomagnetic fields

Simulation wraps the C++ blochsimulation RK4 kernel for a single run: generates field realizations, loops over spin packets, integrates the Bloch equations, and stores results in trjries.

Simulations manages a collection of Simulation runs. Call run() to execute all runs and saveToPkl() to persist results.

axionbloch.constants

Physical constants and unit conversion factors used throughout the package. All values are astropy.units.Quantity objects so unit arithmetic is fully supported.

Symbol

Description

gamma_p

Proton gyromagnetic ratio (rad Hz T⁻¹)

mu_p

Proton magnetic dipole moment

gamma_Xe129

¹²⁹Xe gyromagnetic ratio (rad Hz T⁻¹)

mu_Xe129

¹²⁹Xe magnetic dipole moment

.. py:module:: axionbloch.constants

Physical constants and unit conversion factors used throughout the axionbloch package.

Provides:

  • Nuclear magneton helper mu_N

  • Proton and Xe-129 gyromagnetic ratios / magnetic moments

All constants with units are astropy Quantity objects.

.. py:function:: mu_N(m) :module: axionbloch.constants

Return the nuclear magneton ehbar / (2m) for a nucleus of mass m.

axionbloch.utils

General-purpose helpers shared across the package.

Item

Description

check(arg)

Debug-print a variable with its name, value, and unit

PhysicalObject

Base class with useCommonUnits() for consistent unit display

Lorentzian, Gaussian, ExpCos

Common NMR lineshape functions

estimateLorzfit, estimateGaussfit

Initial-parameter estimators for curve fits

getDateAndTime()

Compact timestamp string YYYYMMDD_HHMMSS

sci_fmt(x, pos)

Scientific-notation tick formatter for matplotlib

high_contrast_extended

Colour palette for multi-line plots

.. py:module:: axionbloch.utils

General-purpose utilities for the axionbloch package.

Contents

  • Debugging helpers: :func:check

  • Simple polynomial and lineshape functions: :func:poly1, :func:poly2, :func:Lorentzian, :func:dualLorentzian, :func:tribLorentzian

  • Curve-fit estimators: :func:estimateLorzfit, :func:estimatedualLorzfit

  • Unit-aware base class: :class:PhysicalObject

  • Misc: :func:giveDateAndTime, :func:sci_fmt, :func:save_phys_quantity

.. py:function:: Add_vector(ax, start=None, end=None, mutation_scale=10, lw=1.6, color=’k’, alpha=1, zorder=5, linestyle=’-’, verbose=False) :module: axionbloch.utils

.. py:class:: Arrow3D(xs, ys, zs, *args, **kwargs) :module: axionbloch.utils

Bases: :py:class:~matplotlib.patches.FancyArrowPatch

.. py:method:: Arrow3D.do_3d_projection(renderer=None) :module: axionbloch.utils

.. py:method:: Arrow3D.draw(renderer) :module: axionbloch.utils

  Draw the Artist (and its children) using the given renderer.

  This has no effect if the artist is not visible (`.Artist.get_visible`
  returns False).

  :param renderer:
  :type renderer: `~matplotlib.backend_bases.RendererBase` subclass.

  .. rubric:: Notes

  This method is overridden in the Artist subclasses.

.. py:method:: Arrow3D.set(*, agg_filter=, alpha=, animated=, antialiased=, arrowstyle=, capstyle=, clip_box=, clip_on=, clip_path=, color=, connectionstyle=, edgecolor=, facecolor=, fill=, gid=, hatch=, hatch_linewidth=, in_layout=, joinstyle=, label=, linestyle=, linewidth=, mouseover=, mutation_aspect=, mutation_scale=, patchA=, patchB=, path_effects=, picker=, positions=, rasterized=, sketch_params=, snap=, transform=, url=, visible=, zorder=) :module: axionbloch.utils

  Set multiple properties at once.

  Supported properties are

  Properties:
      agg_filter: a filter function, which takes a (m, n, 3) float array and a dpi value, and returns a (m, n, 3) array and two offsets from the bottom left corner of the image
      alpha: float or None
      animated: bool
      antialiased or aa: bool or None
      arrowstyle: [ '-' | '<-' | '->' | '<->' | '<|-' | '-|>' | '<|-|>' | ']-' | '-[' | ']-[' | '|-|' | ']->' | '<-[' | 'simple' | 'fancy' | 'wedge' ]
      capstyle: `.CapStyle` or {'butt', 'projecting', 'round'}
      clip_box: `~matplotlib.transforms.BboxBase` or None
      clip_on: bool
      clip_path: Patch or (Path, Transform) or None
      color: :mpltype:`color`
      connectionstyle: [ 'arc3' | 'angle3' | 'angle' | 'arc' | 'bar' ]
      edgecolor or ec: :mpltype:`color` or None
      facecolor or fc: :mpltype:`color` or None
      figure: `~matplotlib.figure.Figure` or `~matplotlib.figure.SubFigure`
      fill: bool
      gid: str
      hatch: {'/', '\\', '|', '-', '+', 'x', 'o', 'O', '.', '*'}
      hatch_linewidth: unknown
      in_layout: bool
      joinstyle: `.JoinStyle` or {'miter', 'round', 'bevel'}
      label: object
      linestyle or ls: {'-', '--', '-.', ':', '', (offset, on-off-seq), ...}
      linewidth or lw: float or None
      mouseover: bool
      mutation_aspect: float
      mutation_scale: float
      patchA: `.patches.Patch`
      patchB: `.patches.Patch`
      path_effects: list of `.AbstractPathEffect`
      picker: None or bool or float or callable
      positions: unknown
      rasterized: bool
      sketch_params: (scale: float, length: float, randomness: float)
      snap: bool or None
      transform: `~matplotlib.transforms.Transform`
      url: str
      visible: bool
      zorder: float

.. py:function:: DTRC_filter(signal, samprate: float, TC: float, order: int) :module: axionbloch.utils

Discrete-Time RC Filter

.. rubric:: References

[1] Zurich Instruments, MFIJ User Manual 500 kHz / 5 MHz Impedance Analyzer P185 6.4. Discrete-Time Filters https://docs.zhinst.com/pdf/ziMFIA_UserManual.pdf

.. py:function:: ExpCos(t=None, Amp=None, T2=None, nu=None, phi0=None, offset=None) :module: axionbloch.utils

Exponentially decay cos wave s=Aexp(-t/T2)sin(2pinu*t+phi0)+offset

.. py:function:: ExpCosiSin(t=None, Amp=None, T2=None, nu=None, phi0=None, offsetx=None, offsety=None) :module: axionbloch.utils

Exponentially decay cos wave s=Aexp(-t/T2)sin(2pinu*t+phi0)+offset

.. py:function:: ExpCosiSinResidual(params, t, s) :module: axionbloch.utils

Exponentially decay cos wave s=Aexp(-t/T2)sin(2pinu*t+phi0)+offset

.. py:function:: Gaussian(x, center, sigma, area, offset) :module: axionbloch.utils

Return the value of the Gaussian function

                      area                 1  (x-center)^2
   offset + ─────────────────────── exp(- ─── ──────────────)
               sigma * sqrt(2 Pi)          2     sigma^2

:param x: argument of the Lorentzian function :type x: scalar or array_like :param center: the position of the Lorentzian peak :type center: scalar :param sigma: variance of Gaussian function. FWHM = 2.35482 sigma, FWTM = 4.29193 sigma :type sigma: scalar :param area: area under the Lorentzian curve (without taking offset into consideration) :type area: scalar :param offset: offset for the curve :type offset: scalar

:returns: the value of the Lorentzian function :rtype: ndarray or scalar

.. rubric:: Examples

.. rubric:: References

https://en.wikipedia.org/wiki/Gaussian_function

.. py:function:: Init_0090sphere(ax, verbose=False) :module: axionbloch.utils

.. py:function:: Init_3020sphere(ax, verbose=False) :module: axionbloch.utils

.. py:function:: LIAFilterFunction(x=None, tau=None, order=None) :module: axionbloch.utils

.. py:function:: LIAFilterHomega(datax=None, datay=None, frequency=None, taun=None, order=None) :module: axionbloch.utils

Return the complex array H(ω) for digital filter correction H(ω) = 1 / (1 + 2 * np.pi * 1j * frequency * taun) ** order

:param frequency: Dmodualtor frequency of lock-in amplifier. :type frequency: scalar :param taun: taun equals TC (Time constant of the exponential running average filter). :type taun: scalar :param order: Number of cascaded digital filters. :type order: scalar :param verbose: choose True to display processing information :type verbose: bool

:returns: Complex array H(ω) :rtype: ndarray

.. rubric:: Examples

.. rubric:: References

Zurich Instruments, MFIA User Manual, Page 275, 6.4.1. Discrete-Time RC Filter https://docs.zhinst.com/pdf/ziMFIA_UserManual.pdf

.. py:function:: LIAFilterHomegaSquared(datax=None, datay=None, frequency=None, taun=None, order=None) :module: axionbloch.utils

.. py:function:: LIAFilterHomegaSquared1(frequency=None, taun=None, dmodfreq=None, a=None) :module: axionbloch.utils

.. py:function:: LIAFilterHomegaSquared2(frequency=None, taun=None, dmodfreq=None, a=None) :module: axionbloch.utils

.. py:function:: LIAFilterPSD(frequency=None, taun=None, order=None, verbose=False) :module: axionbloch.utils

Return the absolute value of complex array H(ω) for digital filter correction np.abs(H(ω)) = np.abs( 1 / (1 + 2 * np.pi * 1j * frequency * taun) ** order )

:param frequency: Dmodualtor frequency of lock-in amplifier. :type frequency: scalar :param taun: taun equals TC (Time constant of the exponential running average filter). :type taun: scalar :param order: Number of cascaded digital filters. :type order: scalar :param verbose: choose True to display processing information :type verbose: bool

:returns: Absolute value of complex array H(ω) :rtype: ndarray

.. rubric:: Examples

.. rubric:: References

Zurich Instruments, MFIJ User Manual, Page 275, 6.4.1. Discrete-Time RC Filter https://docs.zhinst.com/pdf/ziMFIA_UserManual.pdf

.. py:function:: Lorentzian(x, center, FWHM, area: float = 1.0, offset: float = 0.0) :module: axionbloch.utils

Return the value of the Lorentzian function offset + 0.5FWHMarea / (np.pi * ( (x-center)**2 + (0.5*FWHM)**2 ) )

                      FWHM A
   offset + ───────────────────────
             2π ((x-c)^2+(FWHM/2)^2 )

:param x: argument of the Lorentzian function :type x: scalar or array_like :param center: the position of the Lorentzian peak :type center: scalar :param FWHM: full width of half maximum (FWHM) / linewidth of the Lorentzian peak :type FWHM: scalar :param area: area under the Lorentzian curve (without taking offset into consideration) :type area: scalar :param offset: offset for the curve :type offset: scalar

:returns: the value of the Lorentzian function :rtype: ndarray or scalar

.. rubric:: Examples

.. rubric:: References

Null

.. py:function:: Lorentzian_0edge(x, center, FWHM, area: float = 1.0, offset: float = 0.0) :module: axionbloch.utils

Return the value of the Lorentzian function with values = 0 at edges offset + 0.5FWHMarea / (np.pi * ( (x-center)**2 + (0.5*FWHM)**2 ) )

                      FWHM A
   offset + ───────────────────────
             2π ((x-c)^2+(FWHM/2)^2 )

:param x: argument of the Lorentzian function :type x: scalar or array_like :param center: the position of the Lorentzian peak :type center: scalar :param FWHM: full width of half maximum (FWHM) / linewidth of the Lorentzian peak :type FWHM: scalar :param area: area under the Lorentzian curve (without taking offset into consideration) :type area: scalar :param offset: offset for the curve :type offset: scalar

:returns: the value of the Lorentzian function :rtype: ndarray or scalar

.. rubric:: Examples

.. rubric:: References

Null

.. py:function:: MethanolCS2temp(CSval=None, CSunit=’ppm’, referfreq=1000000.0, tempunit=’K’) :module: axionbloch.utils

.. py:function:: MovAvgByStep(xstamp=None, rawsignal=None, weights=None, step_len=1, verbose=False) :module: axionbloch.utils

A moving average with tunable step length, especially designed for axion signal search.

:param rawsignal: raw signal :type rawsignal: array :param weights: The weights for doing the averaging :type weights: array :param step_len: The step length for doing the moving average. Default to 1. :type step_len: int :param verbose: It is here for no reason. :type verbose: bool

:returns: np.array(prcdsiganl) – The processed signal. :rtype: array

.. py:function:: Npole2station(theta_e=None, phi_e=None, theta_s=None, phi_s=None, verbose=False) :module: axionbloch.utils

return in cartesian coordinates

.. py:class:: PhysicalObject() :module: axionbloch.utils

Bases: :py:class:object

Base class for physical objects with Quantity attributes. Automatically converts units and saves quantities to HDF5.

.. py:method:: PhysicalObject.loadFromH5group(group) :module: axionbloch.utils

  Load all attributes listed in self.quantities and self.generalQuantities
  from an HDF5 group.

.. py:method:: PhysicalObject.loadFromPkl(path: str, verbose: bool = False) :module: axionbloch.utils

  Load an instance of this class from a pickle file.

.. py:method:: PhysicalObject.saveToH5group(group: ~h5py._hl.group.Group, verbose: bool = False) :module: axionbloch.utils

  Save all Quantity attributes to the HDF5 group.

.. py:method:: PhysicalObject.saveToPkl(fileDir: str = ‘’, fileName: str | None = None, overwrite: bool = False, verbose: bool = False) :module: axionbloch.utils

  Save this instance to a pickle file.

.. py:method:: PhysicalObject.useCommonUnits(verbose: bool = False) :module: axionbloch.utils

  Convert all Quantity attributes to their common units.
  Subclasses should define a dict `quantities` mapping attribute names
  to desired units.

.. py:function:: PolyEven(x, C0, C2, C4, C6, C8, C10, center, verbose=False) :module: axionbloch.utils

Return the value of the polynomial C0 + C2 * (x-center)^2 + C4 * (x-center)^4 + C6 * (x-center)^6 + C8 * (x-center)^8 + C10 * (x-center)^10

:param x: argument of the polynomial :type x: scalar or array_like

C0, C2, C4, C6, C8 : scalar

verbose : bool the option for displaying assistive information

:returns: the value of the polynomial :rtype: scalar or array_like

.. rubric:: Examples

.. rubric:: References

Null

.. py:function:: angle_between(v1, v2) :module: axionbloch.utils

Returns the angle in radians between vectors ‘v1’ and ‘v2’::

angle_between((1, 0, 0), (0, 1, 0)) 1.5707963267948966 angle_between((1, 0, 0), (1, 0, 0)) 0.0 angle_between((1, 0, 0), (-1, 0, 0)) 3.141592653589793

.. py:function:: axisfmt_C2K(x, y, format_string) :module: axionbloch.utils

.. py:function:: axisfmt_K2C(x, y, format_string) :module: axionbloch.utils

.. py:function:: boltzmann_probabilities(energies: ~typing.Sequence[~astropy.units.quantity.Quantity], T: ~astropy.units.quantity.Quantity) -> ~numpy.ndarray :module: axionbloch.utils

Compute Boltzmann probabilities for a set of energy eigenstates.

:param energies: energies with units of energy :type energies: sequence of Quantity :param T: Temperature (must be > 0) :type T: Quantity

:returns: Dimensionless probabilities :rtype: np.ndarray

.. py:function:: calculate_fwhm(x, y, peak=True) :module: axionbloch.utils

Calculate the Full Width at Half Maximum (FWHM) of a curve. Works for both peaks and dips.

:param x: The x-values of the curve. :type x: array-like :param y: The y-values of the curve (can be positive or negative). :type y: array-like :param peak: If True, calculate FWHM for a peak (maximum). If False, calculate for a dip (minimum). :type peak: bool

:returns: The FWHM of the curve. :rtype: float

.. py:function:: check(arg) :module: axionbloch.utils

Print information of input arg

.. rubric:: Example

import numpy as np

a = np.zeros((2, 4))

check(a)

a+=1

check(a)

check(len(a))

TERMINAL OUTPUT:

casper-gradient-code\testofcheckpoint.py @45 a : ndarray(array([[0., 0., 0., 0.], [0., 0., 0., 0.]])) [shape=(2, 4)]

casper-gradient-code\testofcheckpoint.py @47 a : ndarray(array([[1., 1., 1., 1.], [1., 1., 1., 1.]])) [shape=(2, 4)]

casper-gradient-code\testofcheckpoint.py @48 len(a) : int(2)

casper-gradient-code\testofcheckpoint.py @49 a.shape : tuple((2, 4)) [len=2]