Internals

BUILTIN_SWEEPS

Dictionary of built-in sweep profiles for motion instructions.

Supported keys:

• :linear – linear sweep s ↦ s
• :min_jerk – minimum-jerk profile 10s^3 - 15s^4 + 6s^5
• :cosine – smooth cosine-based profile

move(atoms, beams, displacement, duration, sweep, dt) -> (moves, nsteps)

Helper to create a scheduled displacement for beams.

Arguments:

• atoms: atom collection (used by modifiers if needed)
• beams: collection of beam objects to move
• displacement: 3-vector total displacement
• duration: total time for the move (seconds)
• sweep: motion profile, either a built-in symbol (:linear, :min_jerk, :cosine) or a custom function s -> f(s) mapping [0,1] → [0,1]
• dt: time step (seconds)

Returns a tuple (moves, nsteps) where moves is a vector of MoveModifier objects and nsteps is the number of steps for the move segment.

ramp(beams, amplitudes_final, ramp_time, dt) -> (ramps, tspan)

Linearly ramp each beam in beams from its current amplitude to the corresponding value in amplitudes_final over ramp_time seconds.

Returns (ramps, nsteps), where ramps is a vector of AmplitudeModifier and nsteps is the number of time steps.

has_beams(sys::System) -> Bool

Return true if sys has any associated AbstractBeam objects.

is_classical(sys::System) -> Bool

Return true if sys has no explicit initial quantum state.

A classical system is characterized here by initial_state === nothing, indicating that only classical degrees of freedom are relevant.

is_quantum(sys::System) -> Bool

Return true if sys contains any quantum field objects (instances of Dynamiq.AbstractField) in its node outputs.

This is used to distinguish purely classical configurations from fully quantum or semiclassical systems.

is_semiclassical(sys::System) -> Bool

Return true if sys combines quantum fields with optical tweezers.

A semiclassical system is one where is_quantum(sys) is true and beams is non-empty, corresponding to quantum internal dynamics coupled to classical external potentials.

per_atom_indices(atom::AbstractAtom, s::AbstractLevel)

Return a list of (index, coefficient) pairs corresponding to placing atom in the pure level s.

This looks up the basis index of s in atom.level_indices and returns a single entry [(idx, 1.0)]. It is used internally by getqstate when building product states across multiple atoms.

per_atom_indices(atom::AbstractAtom, s::Superposition)

Return a list of (index, coefficient) pairs corresponding to the superposition s for a single atom.

For each level–amplitude pair in s.coeffs, the corresponding basis index is extracted from atom.level_indices. The result is a list that encodes the local state of the atom and is used by getqstate to assemble many-body product states.

Base.copy(sys::System) -> System

Create a shallow copy of sys with copied atoms, beams, basis, nodes, and detector specifications.

Atoms and beams are duplicated via copy per element; basis, nodes, and detector_specs are copied container-wise; and the state reference is copied if present. Nodes are shared (they hold compiled field references), so this is safe for concurrent reads but not for concurrent writes to node state.

Base.getindex(sys::System, idx::Int)

Return the idx-th DAG node in sys.nodes.

Base.push!(sys::System, node::AbstractNode)

Append a DAG node to sys.nodes.

Typically called by add_coupling!, add_detuning!, add_decay!, etc. after constructing and building the node.

_play(job::SimulationJob; savefinalstate::Bool=false) -> NamedTuple

Execute a compiled simulation job for a single quantum trajectory shot.

Returns a NamedTuple with:

• detectors: Dict{String, Array} of detector outputs
• times: Vector{Float64} of time points
• final_state: Copy of final quantum state (only if savefinalstate=true)

get_noise_spectrum(noise::LaserPhaseNoiseModel, tspan) -> (freqs, psd_values)

Compute the one-sided laser frequency-noise power spectral density psd_values (units: Hz²/Hz) at a set of logarithmically spaced frequencies freqs (Hz), given a LaserPhaseNoiseModel and a time span tspan.

• The frequency grid is chosen adaptively from the total duration and sampling interval (fundamental and Nyquist limits).
• The PSD combines the servo bump and power-law background components of the model.

Returns a tuple (freqs, psd_values), which can be used for phase-noise synthesis or for plotting.

update_noisy_field_time_refs!(obj, global_time_ref, base_resolve_target=identity) -> obj

Update time references in NoisyField objects so that their noise generation uses a shared global_time_ref.

• If obj is a NoisyField, returns a new NoisyField with the same noise model and RNG but with:
  • global_time_ref set to the provided reference, and
  • the wrapped coupling mapped through base_resolve_target.
• If obj is not a NoisyField, returns it unchanged.

This is called during parameter resolution inside play/play! to ensure time-correlated noise across segments.

build_choi_matrix(input_vecs, outputs, comp_indices)

Construct the unnormalized Choi matrix for the process, restricted to the computational subspace (indices supplied via comp_indices). Reconstructs the action on the operator basis using outputs for |0⟩, |1⟩, |+⟩, |+i⟩.

Returns a complex d^2 × d^2 matrix in the standard Choi (swap) convention: Choi{ik, jl} = [𝔈(|i⟩⟨j|)]{k,l}

choi_to_ptm(choi)

Convert a Choi matrix choi into the corresponding Pauli transfer matrix (PTM) for a single-qubit channel, using the {I, X, Y, Z} basis.

extract_kraus_operators(choi)

Extract Kraus operators from a Choi matrix choi via its singular-value decomposition, discarding numerically negligible singular values.

input_statevectors(system, kets; density_matrix=true)

Generate initial density matrices or statevectors (per density_matrix) for each input state in kets by calling getqstate(system, [ket]) for each. Returns a Vector of statevectors or density matrices.

pauli_input_states(levels::Vector{<:AbstractLevel})

Generate the minimal set of input states for single-qubit process tomography: |0⟩, |1⟩, |+⟩ = (|0⟩+|1⟩)/√2, and |+i⟩ = (|0⟩+i|1⟩)/√2. Returns a Vector of AbstractLevel and Superposition instances, suitable for use in getqstate.

simulate_process(sys, seq, input_states; density_matrix=nothing, shots=1, kwargs...)

Simulate the quantum process for each input state in input_states by calling play for each, properly initializing the system for process tomography analysis.

• If density_matrix is not specified:

  • Uses deterministic evolution (density_matrix = true) when shots == 1.
  • Uses quantum trajectories (density_matrix = false) when shots > 1.

• If density_matrix = true, returns the final density matrix for each input, averaging over all shots if shots > 1.

• If density_matrix = false, runs quantum trajectory simulations. For shots > 1, returns the average output density matrix constructed from the projectors of each trajectory. For shots == 1, emits a warning and returns the projector of the single statevector.

Returns a vector of estimated output density matrices (one per input state).

All additional keyword arguments are forwarded to play.

update!(obj, ::Val{name}, val)

Apply a resolved parameter value val to obj for the parameter named name.

Define methods of this function for each object type that accepts parameters. Called by compile and recompile! for each registered parameter binding.

tweezers_in_col(ta::TweezerArray, col_idx::Int) -> Vector{GaussianBeam}
tweezers_in_col(ta::TweezerArray, col_idxs::AbstractVector{Int}) -> Vector{GaussianBeam}

Return the beams belonging to one or more columns of a TweezerArray.

• tweezers_in_col(ta, j) returns all beams in column j (varying row index).
• tweezers_in_col(ta, col_idxs) concatenates the beams from each requested column.

The beams are returned in row-major order within each column.

tweezers_in_row(ta::TweezerArray, row_idx::Int) -> Vector{GaussianBeam}
tweezers_in_row(ta::TweezerArray, row_idxs::AbstractVector{Int}) -> Vector{GaussianBeam}

Return the beams belonging to one or more rows of a TweezerArray.

• tweezers_in_row(ta, i) returns all beams in row i (varying column index).
• tweezers_in_row(ta, row_idxs) concatenates the beams from each requested row.

These helpers are convenient for addressing or visualizing 1D subsets of a 2D tweezer grid.

copy(tw::TweezerArray) -> TweezerArray

Create a shallow copy of a TweezerArray, duplicating the frequency and amplitude vectors and copying each GaussianBeam.

Useful when deriving modified arrays from an existing configuration without mutating the original.

Base.getindex(t::TweezerArray, i, j)

Index into the underlying beams vector.

This allows TweezerArray to be used like a 2D collection of GaussianBeam objects, e.g. ta[1,2] returns the beam in the first row and the second column.

Beam at (row=i, col=j) has linear index (i-1)*nrow + j

Base.getindex(t::TweezerArray, i...)

Index into the underlying beams vector.

This allows TweezerArray to be used like a 1D collection of GaussianBeam objects, e.g. ta[1] returns the first beam in the flattened row–major ordering.

Base.iterate(t::TweezerArray, state...)

Iterate over the beams in a TweezerArray.

Enables idioms such as

for beam in ta
# do something with each GaussianBeam
end

Base.length(t::TweezerArray)

Return the number of beams in the array (length(ta.beams)), i.e. length(row_freqs) * length(col_freqs).

_resolve(obj, param_values, cache::IdDict)

Fallback resolver: for unhandled types, return obj unchanged.

This ensures that objects which do not depend on parameters pass through the resolution stage transparently.

_resolve(mc::MoveCol, param_values, cache::IdDict)
_resolve(inst::Pulse, param_values, cache::IdDict)
_resolve(inst::On,   param_values, cache::IdDict)
_resolve(inst::Off,  param_values, cache::IdDict)

Resolve parametric fields and timing parameters appearing in instruction objects (MoveCol, Pulse, On, Off), returning fully concrete copies suitable for simulation.

_resolve(x::Number, param_values, cache::IdDict)

Internal helper: numbers are already concrete and are returned unchanged.

_resolve(p::Parameter, param_values, cache::IdDict)

Resolve a scalar Parameter to a concrete numeric value.

If param_values contains an entry for p.name, that value (or parameter) overrides the default p.default / p.std. A nonzero std leads to a random draw mean + randn()*std, enabling static noise sampling; otherwise the mean value is returned.

_resolve(expr::ParametricExpression, param_values, cache::IdDict)

Evaluate a ParametricExpression by first resolving all its arguments and then applying the encoded operation.

Currently supports the binary operators :* and :+, corresponding to multiplication and addition. An error is thrown if an unknown operator symbol is encountered.

_resolve(seq::Sequence, param_values, cache::IdDict)

Resolve all instructions in a Sequence, preserving its time step dt.

Returns a new Sequence with the same dt and resolved instructions.

_resolve(x::String, param_values, cache::IdDict)
_resolve(x::Symbol, param_values, cache::IdDict)
_resolve(x::Function, param_values, cache::IdDict)
_resolve(x::TweezerArray, param_values, cache::IdDict)

Internal helpers: symbols, functions, and TweezerArray instances are treated as already-resolved primitives and returned unchanged.

_resolve(sys::System, param_values, cache::IdDict)

Resolve a System by resolving its atoms and beams, while passing through initial_state, state, basis, nodes, and detector_specs unchanged.

Returns a new System instance with parametric components concretized according to param_values.

_resolve(x::Vector, param_values, cache::IdDict)
_resolve(x::Tuple,  param_values, cache::IdDict)

Recursively resolve all elements of a vector or tuple, threading the same cache so that shared substructures and deferred objects are handled consistently.

_resolve(b::ParametricBeam{T}, param_values, cache::IdDict) where T

Resolve a ParametricBeam{T} by resolving all parametric arguments and constructing a concrete beam of type T.

BeamRabiFrequency

A value type that computes a Rabi frequency from a beam's E-field at the atom's position. Used as the Ω field of a CouplingNode.

Holds a reference to a BeamNode whose _compiled[] field is populated before this value is evaluated (insertion-order guarantee in sys.nodes). Atom position is read from atom.inner.x, which is updated in-place by initialize! before coupling nodes are compiled.

Future derived-value types (e.g. BeamACStarkShift for light-shift detunings) follow the same pattern: store physics inputs, implement _resolve_node_default and _resolve_node_value.

Efield_spherical(beam, r; q_axis=[0,0,1])

Project the lab-frame E field of beam at position r into the spherical components (E0, Eplus, Eminus) with respect to q_axis. Returns (E0, Eplus, Eminus) as ComplexF64.

_add_coupling!(system, atom, level::Pair{<:AbstractLevel,<:AbstractLevel}, Ω;
               beam = nothing, noise = nothing, active = true)

Internal helper: add a single coupling DAG node between two specified atomic sublevels.

Constructs the appropriate node type:

• PlanarCouplingNode if beam isa PlanarBeam
• NoisyCouplingNode if noise isa AbstractNoiseModel
• CouplingNode otherwise

Ω may be a number, Parameter, or ParametricExpression; it is stored in the node and resolved at play time. The node is built with the default (or zero if inactive) coefficient, pushed to sys.nodes, and its compiled field is returned. Returns nothing if either level is not found.

_add_couplings(system, atom, g_levels, e_levels,
               Ω_π, Ω_σ⁺, Ω_σ⁻, noise, active, tol; beam = nothing)

Internal helper: add all allowed couplings between two sets of sublevels.

For each pair (g, e) in g_levels × e_levels, the corresponding matrix element is computed via compute_coupling_strength. Pairs with abs(Ω) < tol are skipped. Remaining couplings are created via _add_coupling! and collected into a vector of Dynamiq.AbstractField instances.

Returns the vector of constructed coupling fields.

compute_coupling_strength(g::AbstractLevel, e::AbstractLevel,
                          atom, Ω_π, Ω_σ⁺, Ω_σ⁻)

Fallback coupling-strength implementation.

For level types without a more specific method, returns 0.0, meaning no allowed dipole coupling is assumed.

compute_coupling_strength(g::FineLevel, e::FineLevel,
                          atom, Ω_π, Ω_σ⁺, Ω_σ⁻)

Compute the electric-dipole coupling strength between two fine-structure levels g and e (J↔J) for given polarization-resolved Rabi frequencies.

Uses the Clebsch–Gordan coefficient ⟨J_g, m_Jg; 1, Δm | J_e, m_Je⟩ and selects Ω_π, Ω_σ⁺, or Ω_σ⁻ according to Δm = m_Je − m_Jg. Returns 0.0 for forbidden transitions (|Δm| > 1) or evaluation failures.

compute_coupling_strength(g::FineLevel, e::HyperfineLevel,
                          atom, Ω_π, Ω_σ⁺, Ω_σ⁻)

Compute the effective coupling strength from a fine-structure level to a hyperfine level.

This is defined symmetrically via compute_coupling_strength(e, g, atom, Ω_π, Ω_σ⁺, Ω_σ⁻).

compute_coupling_strength(g::HyperfineLevel, e::FineLevel,
                          atom, Ω_π, Ω_σ⁺, Ω_σ⁻)

Compute the effective electric-dipole coupling strength from a hyperfine level g (F,m_F) to a fine-structure level e (J, unresolved hyperfine) by summing over compatible projections.

The algorithm decomposes the hyperfine state into |J_g, m_J; I, m_I⟩ components using a clebschgordan factor, then applies dipole Clebsch–Gordan coefficients for J_g → J_e with Δm = m_Je − m_J. Contributions for Δm = 0, ±1 are weighted by Ω_π, Ω_σ⁺, Ω_σ⁻ respectively and accumulated. Returns the total effective Rabi frequency.

compute_coupling_strength(g::HyperfineLevel, e::HyperfineLevel,
                          atom, Ω_π, Ω_σ⁺, Ω_σ⁻)

Compute the electric-dipole coupling strength between two hyperfine levels g and e (F↔F) for given polarization-resolved Rabi frequencies.

The appropriate Clebsch–Gordan coefficient ⟨F_g, m_Fg; 1, Δm | F_e, m_Fe⟩ is evaluated using clebschgordan, and multiplied by Ω_π, Ω_σ⁺, or Ω_σ⁻ depending on Δm = m_Fe − m_Fg ∈ {0, ±1}. Returns 0.0 if the transition is forbidden or if the CG coefficient cannot be evaluated.

dipole_matrix_element(atom::AbstractAtom, g::HyperfineLevel, e::HyperfineLevel; norm=false)

Return the real-valued electric dipole transition amplitude between two hyperfine states

⟨F′ m′ | d_q | F m⟩ in units of the reduced matrix element ⟨J || er || J'⟩_R

in the Racah phase convention (Condon–Shortley phases).

The implementation includes:

• Zeeman (3j) angular factor • Hyperfine recoupling (6j) factor • The spherical tensor phase factor (-1)^q • A global phase choice such that a stretched σ⁺ transition (if it exists) is positive

The formula used is

⟨F′ m′ | d_q | F m⟩ =
    (-1)^q
    (-1)^(F′ - m′)
    ( F′  1   F
     -m′  q   m )
    (-1)^(F′ + J + 1 + I)
    √[(2F′+1)(2F+1)]
    { J′  F′  I
      F   J   1 }

where

q = m′ - m
|q| ≤ 1

The reduced fine-structure matrix element ⟨J′||d||J⟩ is set to 1, so the function returns relative transition amplitudes.

Keyword

norm=false : If true and a stretched σ⁺ transition exists, the result is normalized so that the stretched transition equals +1.

Returns

A real Float64 amplitude in Racah convention. Forbidden transitions return 0.0.

rabi_frequencies(atom, beam; q_axis=[0,0,1], d_red::Float64)

Polarization-resolved Rabi frequencies for an effective E1 operator, independent of J/F. Returns (Ωπ, Ωσ⁺, Ωσ⁻) computed from the spherical components of the beam E-field and a reduced dipole `dred`.

update!(c::GlobalCoupling, ::Val, val)

Update the Rabi rate of a GlobalCoupling by rescaling its operator entries. Called by compile/recompile! for parameter bindings.

update!(d::Detuning, ::Val, val)

Update the energy shift of a Detuning operator to val (rad/s) in-place. Called by compile_node! and recompile_node! for DetuningNode.

update!(j::Jump, ::Val, val)

Rescale the jump operator to a new decay rate val (rad/s) in-place. Updates J.forward entries by sqrt of the rate ratio and clears cached non-Hermitian Hamiltonian and LdagL diagonal.

precompute!(j::Jump, ::Type{<:AbstractMatrix})

Precompute and cache data required for Lindblad master equation propagation.

The diagonal of ((L^\dagger L)) is stored in j.LdagL_diag and reused in the anticommutator part of the Liouvillian.

Returns

• j::Jump with LdagL_diag::Vector{Float64} containing LdagL_diag[j] = ∑ₘ |Lₘⱼ|².

precompute!(j::Jump, ::Type{<:AbstractVector})

Precompute and cache the non-Hermitian contribution (-\mathrm{i}/2 L^\dagger L) for wavefunction Monte Carlo propagation. The result is stored as an Op in j.Hnh and reused during time evolution.

update!(d::BlockadeCoupling, step)

No-op update for blockade couplings. The blockade effect is encoded in the static operator H.

update!(::Detuning, step)

No-op update for static detuning terms. The coefficient remains fixed.

update!(d::GlobalCoupling, step)

No-op update for global couplings. The coefficient is assumed to be handled externally or remain constant in time.

update!(f::StarkShiftAC, step)

Update the AC Stark shift coefficient from the instantaneous beam intensity at the atomic position. The stored coefficient is (lpha I / \hbar) in angular-frequency units.

update!(d::Interaction, step)

No-op update for static pairwise interactions. The operator is fixed and its coefficient is assumed constant unless modified externally.

update!(drive::PlanarCoupling, step)

Update the complex amplitude of a planar coupling using the current atomic position and beam wavevector. This is typically called by the time integrator.

copy(a::NLevelAtom)

Create a deep copy of an NLevelAtom, including position, velocity, polarizabilities, transition wavelengths, and internal cache fields.

Atom(x, v, m, λs, γs, ls)

Simple atomic model with position, velocity, and scalar polarizability.

Fields:

• x::Vector{Float64}: atomic position in real space.
• v::Vector{Float64}: atomic velocity.
• m::Float64: atomic mass.
• λ::NTuple{N,Float64}: wavelengths of relevant transitions (in meters).
• γ::NTuple{N,Float64}: spontaneous emission rates (in Hz).
• ls::NTuple{N,Int}: line-strength factors for each transition.
• alpha::F: polarizability as a function of probe wavelength.

The constructor builds a dispersion-like polarizability

[ \alpha(\omega) \propto \sumi \frac{\text{ls}i \gammai}{\omega - \omegai} ]

from the transition data, where ω(λ) = 2πc/λ.

BasisPolarization(unit_k, polarization)

Convert a Jones vector polarization into the spherical (q-basis) polarization components for a beam with propagation direction unit_k.

• unit_k: 3D unit vector along the beam propagation direction.
• polarization: transverse Jones vector in the lab frame.

Returns a 3-component complex vector corresponding to (\sigma^+), (\pi), and (\sigma^-) components in the quantization basis aligned with unit_k.

fclassical!(dt, atom, beams)

Single-step classical update of position and velocity for atom under forces from beams, using a simple Euler integrator.

• First updates position via (\dot{\mathbf{x}} = \mathbf{v}),
• then updates velocity via (\dot{\mathbf{v}} = \mathbf{F} / m).

fdipole!(dt, psi, atom, drives, jumps, _dpsi)

Update the atomic velocity due to radiation-pressure forces from a set of classical driving fields, using Monte Carlo wavefunction data.

The total scattering rate is obtained from the jump operators as Rtot = Σj ⟨Lj† Lj⟩. For each planar drive, a driven component |dψb⟩ = Hb |ψ⟩ is computed into _dpsi, and a dimensionless weight wb ∝ |ψg* (dψb)e|^2 is formed to quantify how strongly that beam drives the atom. The mean force

F = ħ R_tot Σ_b (w_b / Σ_b' w_b') k_b

is then used to update atom.v over the time step dt.

force(atom, beams)

Compute the total optical dipole force on atom at position atom.x from a set of beams.

• Beams with the same wavelength are summed coherently in the complex field (interference included).
• Beams with different wavelengths contribute incoherently (no interference).

The dipole force is computed as

[ \mathbf{F} = -2 \alpha(\lambda) P \operatorname{Re}\big( E^{\ast} \nabla E \big), ]

where (\alpha(\lambda)) is the polarizability at wavelength (\lambda), (P) is the population of the level, and (E) is the total complex field.

fquantum!(dt, ρ, Hlist, _ρ1, _ρ2; order = 4)

Evolve a density matrix ρ under unitary dynamics generated by a sum of Hamiltonian terms.

The commutator (-i[H, \rho]) is applied using the Op split representation of each term in Hlist, without constructing full dense matrices.

fquantum!(dt, ρ, Hlist, Jlist, _ρ1, _ρ2; order = 4)

Evolve a density matrix ρ by a time step dt under combined Hamiltonian and Lindblad evolution.

Arguments

• dt::Float64: Time step size.
• ρ::Matrix{ComplexF64}: Density matrix, updated in-place.
• Hlist::Vector{Tuple{Base.RefValue{ComplexF64},Op}}: Hamiltonian terms (H_j), each stored as (coeff, H).
• Jlist::Vector{Tuple{Base.RefValue{ComplexF64},Op,Vector{Float64}}}: Lindblad jump data. Each entry is (coeff, L, LdagL_diag) where:
  • coeff[]::ComplexF64: Jump amplitude (physical rate (\gamma = |\text{coeff}|^2)).
  • L::Op: Jump operator (L); only L.forward is used.
  • LdagL_diag::Vector{Float64}: Diagonal of (L^\dagger L) for the anticommutator part.
• _ρ1, _ρ2::Matrix{ComplexF64}: Workspace matrices.

Physics

Integrates the Lindblad master equation

[ \frac{\mathrm{d}\rho}{\mathrm{d}t} = -i[H, \rho] + \sumj \gammaj \Big( Lj \rho Lj^{\dagger} - \tfrac{1}{2} { Lj^{\dagger} Lj, \rho } \Big) ]

with (\gammaj = |\text{coeff}j|^2) and (H = \sumj \text{coeff}j H_j).

fquantum!(dt, qstate, Hlist, _q1, _q2; order = 4)

Evolve a pure state vector qstate forward by a time step dt under a sum of Hamiltonian terms using a multi-stage Magnus-like scheme.

Arguments

• dt::Float64: Time step size.
• qstate::Vector{ComplexF64}: State vector, updated in-place.
• Hlist::Vector{Tuple{Base.RefValue{ComplexF64},Op}}: Hamiltonian terms. Each entry is (coeff, H) where coeff[]::ComplexF64 is a (possibly time-dependent) scalar prefactor and H::Op stores the operator via H.forward and H.reverse.
• _q1, _q2::Vector{ComplexF64}: Work buffers with the same length as qstate.
• order::Int = 4: Number of internal stages used for the time-stepping scheme.

Details

For each term (coeff, H) in Hlist, the Hermitian contribution is applied as

-i dt ( coeff[] * H.forward + conj(coeff[]) * H.reverse ) * qstate

in a multi-stage expansion, accumulating the effect of all terms into qstate. The split storage of H avoids reconstructing full matrices while keeping the Hermitian structure explicit through the forward/reverse parts.

fquantum!(dt, qstate, Hlist, Hnhlist, _q1, _q2; order = 4)

Evolve a pure state vector with both Hermitian and non-Hermitian contributions, e.g. for effective non-Hermitian Hamiltonians in wavefunction Monte Carlo.

• Hlist: Hermitian terms, applied as coeff[] * H.forward + conj(coeff[]) * H.reverse.
• Hnhlist: Non-Hermitian terms, applied as coeff[] * H.forward only.

jump!(psi, J, _prob, _psi1, _psi2, rng)

Apply a randomly selected jump operator from the list J to the state vector psi and return the chosen Jump.

• For a single jump in J, that jump is always applied.
• For multiple jumps, the probability of each process is estimated from the norm of J[k].J * psi, and one is sampled accordingly.

The state vector is projected and renormalized, with a safeguard to avoid division by zero if the jump has vanishing norm.

qme(rho, L, J, tspan; kwargs...)

Quantum master equation solver for density matrices with frozen atoms.

• L: Hamiltonian terms (coeff, Op).
• J: Jump data (coeff, L_op, LdagL_diag).

The evolution is implemented by repeated calls to fquantum! for the combined Hamiltonian and dissipative contributions.

qme_semiclassical(rho, atoms, L, J, tspan; kwargs...)

Quantum master equation solver with semiclassical atomic motion.

Combines:

• Density-matrix evolution under Hamiltonian and Lindblad terms.
• Classical motion updated in parallel by fclassical!.

recoil!(J, rng)

Apply a random recoil kick to the atom associated with jump process J.

A random direction is sampled isotropically, and a momentum kick of magnitude (\hbar k) is applied, where (k = 2\pi / \lambda) is set by the transition wavelength stored in J.atom.lambda[J.transition].

tdse_semiclassical(psi, atoms, H, tspan; kwargs...)

Time-dependent Schrödinger solver with semiclassical atomic motion.

Quantum dynamics is driven by H while classical trajectories evolve under the optical forces from beams. tspan is assumed equidistant.

updatepop!(atom, psi)

Update the internal level populations atom._P from the many-body state vector psi, using the precomputed index lists atom._pidx.

wfmc(psi, H, Hnh, jumps, tspan; kwargs...)

Time-dependent wavefunction Monte Carlo solver.

• H contains Hermitian Hamiltonian terms.
• Hnh contains effective non-Hermitian contributions.
• jumps is a list of Jump processes applied stochastically.
• tspan is an equidistant time grid.

At each step, the state is propagated by fquantum! and a jump is applied with probability set by the norm loss.

wfmc_semiclassical(psi, atoms, H, Hnh, jumps, tspan; kwargs...)

Wavefunction Monte Carlo solver with semiclassical atomic motion.

Combines:

• Quantum trajectory evolution under H and Hnh.
• Classical motion updated by fclassical!.
• Additional dipole-force-induced velocity changes via fdipole!.

expect(op, psi)

Expectation value ⟨psi|Op|psi⟩ for the full operator (forward + reverse).

issame(b1, b2; except = [])

Check whether two basis elements b1 and b2 match on all sites except those in except.

Used internally when constructing local operators on a many-body basis.

mul!(op, psi, psi_)

In-place application of op to psi using a temporary buffer psi_.

mul!(psi_, coeff, op, psi)

In-place application of op with complex prefactor coeff.

Applies the full operator A = A_forward + A_reverse such that psi_ = coeff * A * psi.

operator1(b, atoms, transition, rate; jump = false, blockade = 0)

Build single-atom operator triplets on basis b.

Returns two vectors:

• forward::Vector{Tuple{Int,Int,ComplexF64}}
• reverse::Vector{Tuple{Int,Int,ComplexF64}}

If jump == true, the operator is a one-way collapse operator and reverse is empty. Otherwise, forward corresponds to the "raising"/"lowering" direction defined by transition, and reverse contains the Hermitian-conjugate entries. The blockade parameter optionally suppresses matrix elements when too many atoms occupy a given level.

operator2(b, atoms, transition1, transition2; jump = false)

Two-atom operator constructed from a pair of single-atom transitions.

Returns (forward, reverse) triplet lists for use in Op(forward, reverse, b.dim).

op * psi

Apply the full operator (forward + reverse) to a state vector psi.

coeff * op

Scale the full operator (forward + reverse) by a scalar coeff.

op1 + op2

Sum of two Op objects, implemented by converting both to sparse matrices and constructing a new Op. This is simple and robust, and can be optimized further if needed.

sparse(op::Op)

Reconstruct the sparse matrix representation from an Op.

Returns the full operator A = A_forward + A_reverse, built from op.forward and op.reverse triplets.

build_detector(specs, tspan, vals, resolve_target, system) -> Vector{AbstractDetector}

Instantiate detectors for a simulation segment.

Since detector specs travel with the system and are deep-copied together, atom references in specs already point to the correct copied atoms. resolve_target is still applied for consistency and for resolving deferred parameters.

Arguments

• spec::Vector{DetectorSpec}: Detector specification.
• tspan::Vector{Float64}: Time vector for this segment.
• vals::AbstractVector{T}
• resolve_target::Function: Mapping from original to resolved objects (e.g. for deferred parameters).
• system: System object providing quantum state, basis, etc.

write!(d::CoherenceDetector, i)

Record the current coherence value at time step i.

For a pure state (\vert\psi\rangle), the coherence (\rho{ij}) is computed as (\psii \psij^{\ast}) summed over all basis elements corresponding to (|i\rangle\langle j|). For a density matrix (\rho), the corresponding matrix elements (\rho{ij}) are summed.

write!(d::FieldDetector{<:AbstractField}, i)

Sample the field amplitude at time step i and store it in d.vals[i].

The value recorded is the scalar coefficient _coeff[] of the underlying field or coupling, accessed via base_coupling(d.obj).

write!(d::MotionDetector{<:AbstractBeam}, i)

Sample the reference position of the specified beam at time step i, writing each requested coordinate into d.vals[i, :].

write!(d::MotionDetector{NLevelAtom}, i)

Sample the position of the specified atom at time step i, writing each requested coordinate into d.vals[i, :].

write!(detectors::Vector{MotionDetector{A}}, i)

Sample all motion detectors in the vector at time step i.

write!(d::PhotoDetector, i)

Register a single detection event at time step i by incrementing the count in d.vals[i].

prep!(d::PopulationDetector)

Update the atom's internal population d.atom._P[d.level] by computing the current population of the monitored level from the system's quantum state.

For a pure state (\vert\psi\rangle), the population is

[ Pi = \sum{k \in \mathcal{S}i} |\psik|^2, ]

where (\mathcal{S}_i) is the support of the projector for level i. For a density matrix (\rho), the population is

[ Pi = \sum{k \in \mathcal{S}i} \rho{kk}. ]

update!(m::AmplitudeModifier, i)

Set the complex amplitude _coeff[] of the underlying field or beam to m.vals[i] at time step i, if i is within bounds.

update!(m::MoveModifier, i)

Increment the beam position r0 at time step i by the stored displacement m.vals[i] on the components listed in m.dims, if i is within bounds.

update!(m::PositionModifier, i)

Overwrite the beam position components r0[d] at time step i using the stored trajectory m.vals[i].

Only the indices listed in m.dims are updated.

GHz

Frequency unit: gigahertz.

• Definition: (1\,\text{GHz} = 10^9\,\text{Hz}).

Hz

Frequency unit: hertz.

• Definition: (1\,\text{Hz} = 1\,\text{s}^{-1}).

J

Energy unit: joule.

• Definition: (1\,\text{J} = 1\,\text{kg}\,\text{m}^2\,\text{s}^{-2}).

K

Temperature unit: kelvin.

• Definition: SI base unit of thermodynamic temperature.

MHz

Frequency unit: megahertz.

• Definition: (1\,\text{MHz} = 10^6\,\text{Hz}).

THz

Frequency unit: terahertz.

• Definition: (1\,\text{THz} = 10^{12}\,\text{Hz}).

W

Power unit: watt.

• Definition: (1\,\text{W} = 1\,\text{J/s}).

a0

Bohr radius (a_0).

• Physical meaning: characteristic length scale of the hydrogen atom.
• SI units: meters (m).

amu

Atomic mass unit.

• Physical meaning: 1 unified atomic mass unit.
• SI units: kilograms (kg).

c

Speed of light in vacuum.

• SI units: meters per second (m/s).

cm

Length unit: centimeter.

• Definition: (1\,\text{cm} = 10^{-2}\,\text{m}).

e

Elementary charge.

• SI units: coulombs (C).

g

Mass unit: gram.

• Definition: (1\,\text{g} = 10^{-3}\,\text{kg}).

h

Planck constant (h = 2\pi\hbar).

• Physical meaning: quantum of action.
• SI units: joule seconds (J·s).

hbar

Reduced Planck constant (\hbar).

• Physical meaning: quantum of action divided by (2\pi).
• SI units: joule seconds (J·s).

kHz

Frequency unit: kilohertz.

• Definition: (1\,\text{kHz} = 10^3\,\text{Hz}).

kb

Boltzmann constant (k_B).

• Physical meaning: conversion between temperature and energy.
• SI units: joules per kelvin (J/K).

kg

Mass unit: kilogram.

• Definition: SI base unit of mass.

m

Length unit: meter.

• Definition: SI base unit of length.

mJ

Energy unit: millijoule.

• Definition: (1\,\text{mJ} = 10^{-3}\,\text{J}).

mK

Temperature unit: millikelvin.

• Definition: (1\,\text{mK} = 10^{-3}\,\text{K}).

mW

Power unit: milliwatt.

• Definition: (1\,\text{mW} = 10^{-3}\,\text{W}).

m_e

Electron mass.

• SI units: kilograms (kg).

m_p

Proton mass.

• SI units: kilograms (kg).

mg

Mass unit: milligram.

• Definition: (1\,\text{mg} = 10^{-6}\,\text{kg}).

mm

Length unit: millimeter.

• Definition: (1\,\text{mm} = 10^{-3}\,\text{m}).

ms

Time unit: millisecond.

• Definition: (1\,\text{ms} = 10^{-3}\,\text{s}).

nK

Temperature unit: nanokelvin.

• Definition: (1\,\text{nK} = 10^{-9}\,\text{K}).

nm

Length unit: nanometer.

• Definition: (1\,\text{nm} = 10^{-9}\,\text{m}).

ns

Time unit: nanosecond.

• Definition: (1\,\text{ns} = 10^{-9}\,\text{s}).

ps

Time unit: picosecond.

• Definition: (1\,\text{ps} = 10^{-12}\,\text{s}).

s

Time unit: second.

• Definition: SI base unit of time.

ε0

Vacuum permittivity (\varepsilon_0).

• SI units: farads per meter (F/m).

μ0

Vacuum permeability (\mu_0).

• SI units: newtons per ampere squared (N/A²).

Bohr magneton (\mu_B).

• Physical meaning: quantum of action.
• SI units: joule / T (J·T⁻¹).

µJ

Energy unit: microjoule.

• Definition: (1\,\text{µJ} = 10^{-6}\,\text{J}).

µK

Temperature unit: microkelvin.

• Definition: (1\,\text{µK} = 10^{-6}\,\text{K}).

µW

Power unit: microwatt.

• Definition: (1\,\text{µW} = 10^{-6}\,\text{W}).

µm

Length unit: micrometer.

• Definition: (1\,\text{µm} = 10^{-6}\,\text{m}).

µs

Time unit: microsecond.

• Definition: (1\,\text{µs} = 10^{-6}\,\text{s}).