API

HilbertSpace

Abstract supertype for Hilbert spaces. Concrete subtypes: FockSpace, NLevelSpace, PauliSpace, SpinSpace, PhaseSpace, and ProductSpace for tensor products. Compose with ⊗ or tensor.

ProductSpace{T} <: HilbertSpace

Composite Hilbert space formed by the tensor product of multiple subspaces. The type parameter T is a concrete Tuple type for type-stable storage. Constructed via ⊗ or tensor, not directly.

Operators on a ProductSpace take a positional index specifying which subspace they act on, e.g. Destroy(h, :a, 1).

Examples

julia> FockSpace(:cavity) ⊗ NLevelSpace(:atom, 2)
ℋ(cavity) ⊗ ℋ(atom)

See also ⊗, tensor.

FockSpace(name::Symbol) <: HilbertSpace

Hilbert space for a bosonic mode. Hosts Destroy and Create operators satisfying $[a, a^\dagger] = 1$.

Examples

julia> FockSpace(:cavity)
ℋ(cavity)

See also Destroy, Create, ⊗.

NLevelSpace(name::Symbol, n::Int)
NLevelSpace(name::Symbol, n::Int, ground_state::Int)
NLevelSpace(name::Symbol, levels::Tuple{Vararg{Symbol}})

Hilbert space for an N-level system (atoms, qubits, qudits). Hosts Transition operators $|i\rangle\langle j| \cdot |k\rangle\langle l| = \delta_{jk} |i\rangle\langle l|$.

The ground_state (default 1) selects the projector $|g\rangle\langle g|$ that the canonical basis omits; arithmetic keeps it atomic and expand_completeness materialises $|g\rangle\langle g| = 1 - \sum_{k \neq g}|k\rangle\langle k|$ on demand. Levels can be integers (default 1:n) or symbolic names.

Examples

julia> NLevelSpace(:atom, 3)
ℋ(atom)

julia> NLevelSpace(:atom, (:g, :e))
ℋ(atom)

See also Transition, expand_completeness.

PauliSpace(name::Symbol) <: HilbertSpace

Hilbert space for a two-level system. Hosts Pauli operators satisfying $\sigma_j \sigma_k = \delta_{jk} I + i \epsilon_{jkl} \sigma_l$.

Examples

julia> PauliSpace(:spin)
ℋ(spin)

See also Pauli, SpinSpace.

SpinSpace(name::Symbol) <: HilbertSpace

Hilbert space for collective spin angular momentum. Hosts Spin operators satisfying $[S_j, S_k] = i \epsilon_{jkl} S_l$. The algebra is independent of the spin size $S$, which enters only via the QuantumOpticsBase.SpinBasis chosen for numeric evaluation.

Examples

julia> SpinSpace(:S)
ℋ(S)

See also Spin, PauliSpace.

PhaseSpace(name::Symbol) <: HilbertSpace

Hilbert space for position and momentum quadratures. Hosts Position and Momentum operators satisfying $[p, x] = -i$; arithmetic canonicalizes products to place position left of momentum.

Examples

julia> PhaseSpace(:osc)
ℋ(osc)

See also Position, Momentum, FockSpace.

⊗(spaces::HilbertSpace...)

Create a ProductSpace from multiple Hilbert spaces. Flattens nested ProductSpace arguments so that (A ⊗ B) ⊗ C == A ⊗ B ⊗ C. Unicode input \otimes<tab>; ASCII alias tensor.

Examples

julia> FockSpace(:a) ⊗ FockSpace(:b) ⊗ NLevelSpace(:atom, 2)
ℋ(a) ⊗ ℋ(b) ⊗ ℋ(atom)

See also ProductSpace, tensor.

tensor(spaces::HilbertSpace...)

ASCII alias for ⊗: create a ProductSpace from multiple Hilbert spaces.

Examples

julia> tensor(FockSpace(:a), FockSpace(:b))
ℋ(a) ⊗ ℋ(b)

See also ⊗, ProductSpace.

QSym <: QField

Abstract type for fundamental (leaf) operators in the expression tree.

Every QSym carries three fields identifying its site:

• name::Symbol — display name
• space_index::Int — position in ProductSpace (1 for single spaces)
• index::Index — symbolic summation index (Index), or NO_INDEX

Concrete subtypes: Destroy, Create, Transition, Pauli, Spin, Position, Momentum.

QAdd <: QField

The sole compound expression type: a sum of eagerly-ordered operator products.

All arithmetic on QSym operators produces a QAdd. Iterating over a QAdd yields Pair{QTerm, CNum} entries; read term.ops for the operator sequence and term.ne for the scoped index constraints.

See also QTerm, prefactor, operators, Σ, constraint_pairs.

QTerm

A single entry of a QAdd sum: an ordered operator product plus the pairwise index inequality constraints that scope that product.

Fields

• ops::Vector{Op} — the ordered operator product
• ne::Vector{NonEqualPair} — pairwise α ≠ β constraints (canonicalized)

QTerm is the dict key in QTermDict; iterating a QAdd yields Pair{QTerm, CNum}, and callers reach term.ops / term.ne directly.

Op <: QSym

The single concrete leaf operator. A kind::OpKind tag selects the physical role (annihilation, creation, transition, Pauli, spin, position, momentum); the remaining fields are shared storage interpreted per kind:

• name::Symbol: display name
• space_index::Int: position in a ProductSpace
• index::Index: symbolic site index, or NO_INDEX
• l1, l2, g, nlev::Int32: packed level/axis data. Transition uses all four (i, j, ground state, number of levels); Pauli/Spin store the axis in l1; Fock/PhaseSpace leave them zero.

Construct via the role-named functions Destroy, Create, Transition, Pauli, Spin, Position, Momentum; test the role via is_destroy and siblings, or read it with optype. Collapsing the former per-type hierarchy into one concrete struct makes the operator vector concrete-eltype, so the per-operator hooks dispatch statically. See docs/src/devdocs.md.

prefactor(s::QAdd) -> Complex{Num}

Return the Complex{Num} prefactor of a single-term QAdd.

Throws ArgumentError if s contains more than one term. For multi-term expressions, iterate over the QAdd directly.

Examples

julia> h = FockSpace(:f);

julia> @qnumbers a::Destroy(h);

julia> prefactor(2 * a' * a)
2

See also operators, sorted_arguments.

operators(eqs::AbstractMeanfieldEquations)

The operator products on the left-hand sides, one per tracked moment. Paired with states by position (states(eqs)[i] == average(operators(eqs)[i])).

operators(s::QAdd) -> Vector{Op}

Return the ordered operator sequence of a single-term QAdd.

Throws ArgumentError if s contains more than one term. For multi-term expressions, iterate over the QAdd directly and read term.ops from each QTerm.

Examples

julia> h = FockSpace(:f);

julia> @qnumbers a::Destroy(h);

julia> length(operators(a' * a))
2

See also prefactor, sorted_arguments.

normal_order(expr::QField) -> QAdd

Route every term of expr through the canonicalization pipeline.

In practice this is the identity on anything built through public arithmetic: *, +, -, ^, commutator, Σ, substitute, and adjoint all canonicalize eagerly, so the result of any such call is already normal-ordered. Reach for normal_order explicitly only when an expression was assembled through low-level internals that bypass the arithmetic, or when interfacing with code that expects a finalizer call. simplify uses it internally before simplifying coefficients.

Examples

julia> h = FockSpace(:f);

julia> @qnumbers a::Destroy(h);

julia> a * a'
1 + a' * a

julia> normal_order(a * a')
1 + a' * a

See also simplify, expand_completeness, normal_to_symmetric, symmetric_to_normal.

qadjoint(x)

Hermitian conjugate that distributes through mixed operator/symbolic expressions.

On a QField returns adjoint(x). On a Number returns adjoint(x). On a SymbolicUtils.BasicSymbolic tree, recurses into arguments so coefficients distribute (e.g. qadjoint(2im * a) becomes -2im * a') rather than producing an opaque conj(...) wrapper.

Aliased as qconj and dagger.

inner_adjoint(x)

Push the adjoint inside ⟨...⟩ averages: rewrites conj(⟨X⟩) as ⟨X†⟩ so the result stays expressed as an average of an operator rather than a conj wrapper around one. Used when building equations of motion where both sides must share the canonical "average-of-operator" form for substitution and hashing.

On non-average sub-expressions, behaves like qadjoint.

symmetric_to_normal(expr) -> QAdd

Convert expr from symmetric (Weyl) to normal ordering. Exact inverse of normal_to_symmetric.

Examples

julia> hf = FockSpace(:f);

julia> @qnumbers a::Destroy(hf);

julia> symmetric_to_normal(normal_to_symmetric(a' * a))
a' * a

julia> hp = PhaseSpace(:osc);

julia> @qnumbers x::Position(hp) p::Momentum(hp);

julia> symmetric_to_normal(normal_to_symmetric(x * p))
x * p

See also normal_to_symmetric, normal_order.

normal_to_symmetric(expr) -> QAdd

Convert expr from normal to symmetric (Weyl) ordering. Exact inverse of symmetric_to_normal.

Examples

julia> hf = FockSpace(:f);

julia> @qnumbers a::Destroy(hf);

julia> normal_to_symmetric(a' * a)
-0.5 + a' * a

julia> hp = PhaseSpace(:osc);

julia> @qnumbers x::Position(hp) p::Momentum(hp);

julia> normal_to_symmetric(x * p)
0.5im + x * p

See also symmetric_to_normal, normal_order.

expand_completeness(q) -> QAdd

Rewrite every ground-state projector $\sigma^{gg}$ in q via the completeness relation $\sigma^{gg} = 1 - \sum_{k \neq g} \sigma^{kk}$.

* keeps $\sigma^{gg}$ atomic; reach for expand_completeness when downstream code needs the projector eliminated, e.g. before converting to a numeric basis where the $\sigma^{kk}$ for $k \neq g$ form the independent degrees of freedom.

Examples

julia> h = NLevelSpace(:atom, 2);

julia> σ11 = Transition(h, :σ, 1, 1);

julia> expand_completeness(σ11)
1 - σ₂₂

See also assume_distinct_index, normal_order.

anticommutator(a, b) -> QAdd

Return the anticommutator $\{a, b\} = a\,b + b\,a$ as a QAdd.

Examples

julia> h = FockSpace(:f);

julia> @qnumbers a::Destroy(h);

julia> anticommutator(a, a')
1 + 2 * a' * a

See also commutator.

is_average(x) -> Bool

Whether x is a symbolic average object created by average.

julia> h = FockSpace(:f);

julia> @qnumbers a::Destroy(h);

julia> is_average(average(a)), is_average(a)
(true, false)

is_indexed_sum(x) -> Bool

Whether x is a moment-layer indexed-sum node created by average on a summed QAdd. The summation scope is recovered with get_sum_indices and get_sum_non_equal.

undo_average(expr) -> QAdd

Recursively strip symbolic averages and return the underlying operator expression. Summation metadata is restored. Also accepts a Symbolics.Equation, returning a Pair{QAdd, QAdd}.

julia> h = FockSpace(:f);

julia> @qnumbers a::Destroy(h);

julia> undo_average(average(a' * a)) == a' * a
true

substitute(expr, d::Dict)

Substitute symbolic parameters and/or operators in expr using dictionary d.

Supports both symbolic coefficient replacement (for example x => 2) and operator replacement. The result is re-canonicalized and returned as a QAdd.

Examples

julia> h = FockSpace(:f);

julia> @qnumbers a::Destroy(h);

julia> @variables x;

julia> substitute(x * a' * a, Dict(x => 2))
2 * a' * a

See also change_index.

simplify(expr::QField) -> QAdd

Normal-order expr, then simplify each coefficient symbolically and drop summation indices that no surviving term depends on.

The operator-level work (commutation, same-site composition, like-term collection) all happens inside normal_order. What simplify adds on top is purely at the symbolic-coefficient layer: Symbolics.expand followed by SymbolicUtils.simplify runs on each surviving prefactor, and any term whose coefficient simplifies to zero is dropped. A final pass removes summation indices that no remaining term references.

That symbolic step is expensive, so reach for simplify as a finalizer when cancellations or accumulated symbolic factors need to be folded; use normal_order for intermediate steps.

Examples

julia> h = FockSpace(:f);

julia> @qnumbers a::Destroy(h);

julia> @variables x y;

julia> expr = (x^2 + 2x*y + y^2) * a' * a - (x + y)^2 * a' * a;

julia> simplify(expr)
0

See also normal_order, expand, expand_completeness.

expand(expr::QField) -> QAdd

Expand the symbolic prefactor of each term via Symbolics.expand.

Examples

julia> h = FockSpace(:f);

julia> @qnumbers a::Destroy(h);

julia> @variables x y;

julia> expand((x + y)^2 * a)
(x^2 + 2x*y + y^2) * a

See also simplify.

@qnumbers ops...

Convenience macro for constructing named quantum operators.

Each argument has the form name::OperatorType(hilbert_space, args...). The macro calls OperatorType(hilbert_space, :name, args...) and binds the result to name in the calling scope. Multiple operators can be declared in one call.

Examples

julia> h = FockSpace(:fock);

julia> @qnumbers a::Destroy(h)
(a,)

See also Destroy, Transition, Pauli, Spin.

Destroy(h::FockSpace, name::Symbol) -> Op

Bosonic annihilation operator $a$ on a FockSpace. Satisfies the canonical commutation relation $[a, a^\dagger] = 1$. The adjoint a' returns the corresponding Create operator. Returns an Op tagged OP_DESTROY.

Examples

julia> h = FockSpace(:cavity);

julia> @qnumbers a::Destroy(h);

julia> a * a'
1 + a' * a

See also Create, FockSpace, @qnumbers.

Create(h::FockSpace, name::Symbol) -> Op

Bosonic creation operator $a^\dagger$ on a FockSpace. The adjoint of Destroy; typically obtained via a' rather than constructed directly. Returns an Op tagged OP_CREATE.

See also Destroy, FockSpace.

Transition(h::NLevelSpace, name::Symbol, i, j) -> Op

Transition operator $|i\rangle\langle j|$ on an NLevelSpace. Satisfies the composition rule $|i\rangle\langle j| \cdot |k\rangle\langle l| = \delta_{jk} |i\rangle\langle l|$, with adjoint $|i\rangle\langle j|^\dagger = |j\rangle\langle i|$. Returns an Op tagged OP_TRANSITION, packing i, j, the ground state, and the number of levels into l1, l2, g, nlev.

The carried ground-state level and level count let the eager arithmetic keep $\sigma^{gg}$ atomic in canonical form (use expand_completeness to materialise $|g\rangle\langle g| = 1 - \sum_{k \neq g}|k\rangle\langle k|$ when needed).

Examples

julia> h = NLevelSpace(:atom, 3);

julia> σ = Transition(h, :σ, 1, 2);

julia> σ * σ'
σ₁₁

julia> σ' * σ
σ₂₂

See also NLevelSpace, expand_completeness, @qnumbers.

Pauli(h::PauliSpace, name::Symbol, axis::Int) -> Op

Pauli operator $\sigma_x, \sigma_y, \sigma_z$ on a PauliSpace. The axis argument selects the component: 1 = x, 2 = y, 3 = z. Hermitian (σ' == σ) and satisfies $\sigma_j \sigma_k = \delta_{jk} I + i \epsilon_{jkl} \sigma_l$. Returns an Op tagged OP_PAULI with the axis stored in l1.

Examples

julia> h = PauliSpace(:s);

julia> σx = Pauli(h, :σ, 1); σy = Pauli(h, :σ, 2);

julia> σx * σy
im * σz

julia> σx * σx
1

See also PauliSpace, Spin.

Spin(h::SpinSpace, name::Symbol, axis::Int) -> Op

Angular momentum operator $S_x, S_y, S_z$ on a SpinSpace. The axis argument selects the component: 1 = x, 2 = y, 3 = z. Hermitian (S' == S) and satisfies $[S_j, S_k] = i \epsilon_{jkl} S_l$ (applied eagerly by *). Returns an Op tagged OP_SPIN with the axis stored in l1.

Examples

julia> h = SpinSpace(:S);

julia> Sx = Spin(h, :S, 1); Sy = Spin(h, :S, 2);

julia> Sy * Sx
-im * Sz + Sx * Sy

See also SpinSpace, Pauli.

Position(h::PhaseSpace, name::Symbol) -> Op

Position (quadrature) operator $x$ on a PhaseSpace. Hermitian (x' == x) and related to Fock operators by $x = (a + a^\dagger)/\sqrt{2}$. Canonical pair with Momentum: $[p, x] = -i$. Returns an Op tagged OP_POSITION.

Examples

julia> h = PhaseSpace(:osc);

julia> @qnumbers x::Position(h) p::Momentum(h);

julia> p * x
-im + x * p

See also Momentum, PhaseSpace.

Momentum(h::PhaseSpace, name::Symbol) -> Op

Momentum (quadrature) operator $p$ on a PhaseSpace. Hermitian (p' == p) and related to Fock operators by $p = i(a^\dagger - a)/\sqrt{2}$. Canonical pair with Position: $[p, x] = -i$. Returns an Op tagged OP_MOMENTUM.

See also Position, PhaseSpace.

optype(o::Op) -> OpKind

Return the OpKind tag of o.

operator_name(op::Op) -> Symbol

The display name of op (recovered from the intern table).

is_destroy(o::Op) -> Bool

True iff o is a Destroy (annihilation) operator. One of the role predicates that replace isa on the collapsed Op type; see also is_create, is_transition, is_pauli, is_spin, is_position, is_momentum, optype.

is_create(o::Op) -> Bool

True iff o is a Create (creation) operator. See is_destroy.

is_transition(o::Op) -> Bool

True iff o is a Transition operator. See is_destroy.

is_pauli(o::Op) -> Bool

True iff o is a Pauli operator. See is_destroy.

is_spin(o::Op) -> Bool

True iff o is a Spin operator. See is_destroy.

is_position(o::Op) -> Bool

True iff o is a Position operator. See is_destroy.

is_momentum(o::Op) -> Bool

True iff o is a Momentum operator. See is_destroy.

meanfield(ops, H, J=QField[]; kwargs...)
meanfield(op, H, J=QField[]; kwargs...)

Equations of motion for the averages of ops under the Hamiltonian H and the collapse operators J. The result is the Heisenberg (quantum Langevin) equation with noise neglected: each jump adds a Lindblad term $\sum_i r_i \left(J_i^\dagger O J_i - \tfrac{1}{2}\{J_i^\dagger J_i, O\}\right)$. Returns a MeanfieldEquations, or a NoiseMeanfieldEquations when efficiencies is supplied (measurement backaction).

The returned right-hand sides are left unsimplified; call simplify! (or SymbolicUtils.simplify) on the equations you want to inspect.

Arguments

• ops: the operator(s) whose equations of motion are derived.
• H: the operator defining the system Hamiltonian.
• J: the collapse (jump) operators of the system; defaults to none.

Keyword arguments

• Jdagger=adjoint.(J): the adjoints of the jump operators.
• rates=ones(length(J)): decay rates for J. A square Matrix instead selects collective dissipation across the jump vector.
• efficiencies=nothing: detector efficiencies per jump. When given, the result is a NoiseMeanfieldEquations carrying a measurement-backaction noise drift.
• direction=Forward(): Forward for ordinary evolution, Backward for retrodiction.
• order=nothing: if set, a cumulant_expansion to this order is applied. An Int caps every subspace; a Vector{Int} caps each Hilbert subspace separately.
• mix_choice=maximum: for a Vector order, how to combine the per-subspace caps on a term that acts on several subspaces.
• iv=ModelingToolkitBase.t_nounits: the independent (time) variable of the system.

commutator(a, b) -> QAdd

Return the commutator $[a, b] = a\,b - b\,a$ as a QAdd.

Examples

julia> h = FockSpace(:f);

julia> @qnumbers a::Destroy(h);

julia> commutator(a, a')
1

See also anticommutator, normal_order.

acts_on(expr) -> Vector{Int}

Sorted unique space_index values that expr acts on. Works on QSym, QAdd, averaged BasicSymbolic expressions, and Numbers (Int[]).

julia> h = FockSpace(:a) ⊗ NLevelSpace(:b, 2);

julia> @qnumbers a::Destroy(h, 1) σ::Transition(h, 1, 2, 2);

julia> acts_on(a' * a), acts_on(a' * σ)
([1], [1, 2])

AbstractMeanfieldEquations

Supertype of MeanfieldEquations and NoiseMeanfieldEquations.

MeanfieldEquations

Concrete equation set produced by meanfield when no measurement backaction is requested. The authoritative cumulant hierarchy is graph; the vector fields are a derived view regenerated from it.

Fields

• graph: the source of truth (moments, drifts, treatments, sys, ctx).
• equations: the averaged differential equations (derived view).
• operator_equations: the same equations at the operator level.
• states: the averages on the left-hand sides.
• operators: the operators on the left-hand sides.
• hamiltonian: the system Hamiltonian.
• jumps, jumps_dagger: the collapse operators and their adjoints.
• rates: the decay rates corresponding to jumps.
• iv: the independent (time) variable.
• order: the cumulant-expansion order, or nothing.
• direction: Forward or Backward evolution.

average(expr) -> BasicSymbolic | Number

Build the symbolic average $\langle \mathrm{expr} \rangle$. Distributes over sums, pulls c-number prefactors out, leaves scalars unchanged. Displayed as ⟨…⟩.

julia> h = FockSpace(:f);

julia> @qnumbers a::Destroy(h);

julia> avg = average(a' * a);

julia> is_average(avg)
true

julia> avg + 1
1 + ⟨a' * a⟩

See also undo_average, numeric_average.

cumulant_expansion(expr, order; mix_choice=maximum)

Expand every average in expr whose order exceeds order into lower-order moments via the joint-cumulant identity, neglecting the joint cumulant above order. order is an Int (one cap for all subspaces) or a Vector{Int} (a cap per Hilbert subspace, combined by mix_choice on mixed terms). The output is raw; apply SymbolicUtils.simplify for a canonical form.

See also: https://en.wikipedia.org/wiki/Cumulant#Joint_cumulants

Examples

julia> h = FockSpace(:a) ⊗ FockSpace(:b);

julia> @qnumbers a::Destroy(h, 1) b::Destroy(h, 2);

julia> cumulant_expansion(average(a * b), 1)
⟨b⟩*⟨a⟩

cumulant_expansion(eqs::MeanfieldEquations, order)

Cumulant-expand every RHS of eqs to order, returning a new MeanfieldEquations of the same shape with that order stored. The system's mix_choice (held in eqs.graph.sys) drives the mixed-order truncation. Noise systems take their order through meanfield, not here.

cumulant(op, n=get_order(op))

The n-th joint cumulant of op (either an operator or an average), the signed sum over partitions of its operator factors into n blocks. The result is raw and unsimplified; apply SymbolicUtils.simplify for a canonical form.

Examples

julia> h = FockSpace(:a) ⊗ FockSpace(:b);

julia> @qnumbers a::Destroy(h, 1) b::Destroy(h, 2);

julia> cumulant(a * b)
⟨a * b⟩ - ⟨b⟩*⟨a⟩

julia> cumulant(a * b, 1)
⟨a * b⟩

julia> cumulant(a * b, 3)
0

get_order(expr)

Highest moment order appearing in expr, the order used to decide whether a term is expanded by cumulant_expansion. Numbers have order 0, a bare operator order 1, and a product the number of its operator factors.

Examples

julia> h = FockSpace(:a) ⊗ FockSpace(:b);

julia> @qnumbers a::Destroy(h, 1) b::Destroy(h, 2);

julia> get_order(a)
1

julia> get_order(a * b)
2

julia> get_order(1)
0

states(eqs::AbstractMeanfieldEquations)

The average variables $⟨…⟩$ evolved by eqs, i.e. the dynamical unknowns, one per tracked moment. Paired with operators by position.

moments(eqs::AbstractMeanfieldEquations)

Ordered operator$↔$moment map: each tracked operator product to its average $⟨op⟩$. The correspondence the cumulant hierarchy is built on.

moment_variable_map(eqs::AbstractMeanfieldEquations)

Ordered map from each moment average $⟨op⟩$ to the time-dependent MTK unknown u(t) that System assigns it. The bridge initial_values, get_solution and parameter_map resolve moments against.

closure_report(eqs::AbstractMeanfieldEquations; filter_func=nothing)

Summarise the closure status of eqs:

• n_states: number of tracked moments.
• by_order: tracked-moment count per cumulant order.
• missing: moments still needed to close the system (find_missing).
• closed: whether nothing is missing.

filter_func mirrors find_missing / complete!.

noise_channels(eqs::AbstractMeanfieldEquations)

The monitored measurement channels: one entry per jump with nonzero efficiency, each an independent Brownian in the SDE System builds. Each entry is (; index, jump, rate, efficiency); a deterministic system has none.

CorrelationFunction

Two-time correlation function ⟨op1(t+τ) op2(t)⟩ of a solved mean-field system, constructed via the quantum regression theorem. Build one with the CorrelationFunction(op1, op2, eqs0) constructor, then feed it to ModelingToolkitBase.System together with correlation_u0 and correlation_p0 to solve the τ-evolution, or to Spectrum for the power spectrum.

Fields

• op1, op2: the operators whose correlation is computed.
• op2_ancilla: op2 re-embedded on the ancilla subspace aon_ancilla (internal).
• aon_ancilla: Hilbert-subspace index of the ancilla carrying op2 (internal).
• eqs: the closed equations of motion in the delay τ.
• eqs0: the steady-state system the correlation is built on.
• τ: the delay independent variable.
• steady_state: if true, ambient averages are held constant; if false, they are evolved alongside the τ-states.
• ambient: the τ-system's ambient averages (steady-state coefficients, not τ-states), computed once from the closed eqs so the u0/System/spectrum paths need not re-walk it.

ModelingToolkitBase.System(eqs::MeanfieldEquations; name)
ModelingToolkitBase.System(eqs::NoiseMeanfieldEquations; name)

Build an MTK System from the equation set: one u(t) variable per state, the drift resolved to those variables. A NoiseMeanfieldEquations becomes an SDE with one independent Brownian column per monitored channel, each carrying that channel's noise drift (sign +1 Forward, -1 Backward). RHS leaves not matching a state become parameters.

ModelingToolkitBase.System(c::CorrelationFunction; name)

MTK System for the τ-evolution. Ambient steady-state averages become parameters (values supplied via correlation_p0).

Spectrum(c::CorrelationFunction, ps)

The Laplace-transformed correlation as a callable: S(ω, u_end, p0) returns the symmetric power spectrum 2·Re{∫₀^∞ g(τ)e^{-iωτ}dτ}. Solves (iω·I - A)X̃ = x̃(0) where A is the τ-system Jacobian at steady state and x̃ = x - x_∞ is centred so the one-sided transform converges. The linear/anti-linear split (the conjugate columns from get_adjoints=false) is handled by the 2n augmentation.

Examples

julia> h = FockSpace(:cavity);

julia> @qnumbers a::Destroy(h);

julia> eqs = meanfield([a' * a], a' * a, [a]; rates = [1.0], order = 2);

julia> c = CorrelationFunction(a', a, eqs);

julia> Spectrum(c)
ℱ(⟨a' * a⟩)(ω)

correlation_u0(c, u_end)

Initial values for the τ-evolution of the CorrelationFunction c. Each ancilla state ⟨X(τ) op2ancilla⟩ starts from the steady-state value ⟨X op2⟩ (with op2 on its original subspace), looked up in `uend`.

See also: CorrelationFunction, correlation_p0.

correlation_p0(c, u_end, ps_p0)

Parameter dictionary for the τ-ODE of the CorrelationFunction c: the user parameters ps_p0 together with the numeric values of the ambient steady-state averages, looked up in u_end.

See also: CorrelationFunction, correlation_u0.

Index(h::HilbertSpace, name::Symbol, range, space)
Index(h::HilbertSpace, name::Symbol, range, space_index::Int)

Symbolic summation index for site- or mode-resolved operator expressions.

Represents labels such as i, j, or k that attach to operators and parameters. Pass a subspace type or integer position to specify which factor of a ProductSpace the index ranges over. Use with IndexedOperator to build objects like $a_i$ and with Σ to build sums $\sum_i a_i^\dagger a_i$.

Examples

julia> h = FockSpace(:site) ⊗ FockSpace(:cavity);

julia> i = Index(h, :i, 10, FockSpace(:site));

julia> j = Index(h, :j, 10, 1); # same subspace via integer position

julia> i == j
false

See also has_index, IndexedOperator, Σ, change_index.

IndexedOperator(op::Op, i::Index) -> Op

Return the indexed version of an operator.

IndexedOperator keeps the operator kind and quantum numbers, and only changes the symbolic index label. Use it to build objects such as a_i, σ_j₁₂, or x_k before forming sums with Σ.

Use IndexedOperator(op, NO_INDEX) to remove an index.

Examples

julia> h = FockSpace(:cavity);

julia> @qnumbers a::Destroy(h);

julia> i = Index(h, :i, 5, h);

julia> IndexedOperator(a, i)
a_i

See also Index, Σ, change_index.

IndexedVariable(name::Symbol, i::Index) -> Num

Create a symbolic single-indexed parameter $\mathrm{name}(i)$.

Use this for site-dependent quantities in Hamiltonians, for example detunings or local couplings.

Examples

julia> h = FockSpace(:site) ⊗ FockSpace(:cavity);

julia> i = Index(h, :i, 10, 1);

julia> ω = IndexedVariable(:ω, i)
ω(i)

See also DoubleIndexedVariable.

DoubleIndexedVariable(name::Symbol, i::Index, j::Index; identical::Bool=true) -> Num

Create a symbolic two-index parameter $\mathrm{name}(i, j)$.

Use this for pairwise interactions such as hopping amplitudes or coupling matrices.

Keyword arguments

• identical::Bool = true — if false, the variable evaluates to zero when i == j, enforcing that the parameter is only defined for distinct sites.

Examples

julia> h = FockSpace(:site) ⊗ FockSpace(:cavity);

julia> i = Index(h, :i, 10, 1);

julia> j = Index(h, :j, 10, 1);

julia> J = DoubleIndexedVariable(:J, i, j; identical = false)
J(i, j)

See also IndexedVariable.

Σ(expr, i::Index, non_equal::Vector{Index} = Index[])
Σ(expr, i::Index, j::Index, rest::Index...)
∑(expr, i::Index, ...)

Build symbolic sums over indices, for example $\sum_i expr$.

Returns a QAdd carrying i in its indices field. If expr does not depend on i, the sum is evaluated eagerly as i.range * expr.

The optional non_equal records pairwise inequality constraints i ≠ j per term. Diagonal splitting is performed automatically: if expr carries another free index j on the same Hilbert subspace as i, the contribution at i = j is emitted as a separate diagonal term and the off-diagonal term gains the constraint (i, j).

Multiple positional indices create nested sums: Σ(expr, i, j) is equivalent to Σ(Σ(expr, i), j). The Unicode alias ∑ is also exported.

Examples

julia> h = FockSpace(:site);

julia> @qnumbers a::Destroy(h);

julia> i = Index(h, :i, 3, h);

julia> Σ(IndexedOperator(a', i) * IndexedOperator(a, i), i)
Σ(i=1:3) a_i' * a_i

See also Index, IndexedOperator, constraint_pairs.

change_index(expr, from::Index, to::Index)

Rename an index in an expression by replacing from with to everywhere.

Operator indices are swapped directly. Symbolic prefactors (e.g. IndexedVariable, DoubleIndexedVariable) are substituted via Symbolics.substitute using the sym fields of the indices. DoubleIndexedVariable nodes with identical=false automatically evaluate to zero if the substitution makes both arguments equal.

Examples

julia> h = FockSpace(:site);

julia> @qnumbers a::Destroy(h);

julia> i = Index(h, :i, 5, h); j = Index(h, :j, 5, h);

julia> expr = IndexedVariable(:ω, i) * IndexedOperator(a, i);

julia> change_index(expr, i, j)
ω(j) * a_j

See also get_indices.

change_index(expr, pairs::AbstractDict{Index, Index})

Simultaneous (batched) variant of change_index: substitute every key in pairs by its mapped value in one pass.

Unlike a sequence of single-pair change_index calls, this never produces an intermediate state where a renamed-into index temporarily collides with another encountered index. Use it when the rename describes a permutation or any map whose image overlaps its domain, e.g. swapping two indices:

julia> change_index(expr, Dict(i => j, j => i))   # swap, not fuse

A pair (idx => idx) is a no-op; change_index(expr, Dict()) returns expr unchanged.

get_indices(expr) -> Vector{Index}

Return all symbolic indices that appear in expr.

Returns a Vector{Index} of unique indices found in operator fields and summation metadata. Excludes the sentinel NO_INDEX.

Examples

julia> h = FockSpace(:site);

julia> @qnumbers a::Destroy(h);

julia> i = Index(h, :i, 5, h);

julia> get_indices(IndexedOperator(a, i))
1-element Vector{Index}:
 i

has_index(idx::Index) -> Bool

Return true when idx is an actual symbolic index.

Returns false for the sentinel NO_INDEX, used for operators without an attached index.

assume_distinct_index(q::QAdd, pairs::Vector{Tuple{Index, Index}}) -> QAdd

Re-canonicalize q under the declared pairwise ≠ constraints on free indices, then run expand_completeness so any ground-state projectors that emerge from same-site composition are expanded.

Use this when two free indices semantically range over distinct atoms or modes but no Σ supplies the constraint. SQA cannot infer "different symbol implies different site" automatically: two operators carrying different free indices on the same Hilbert subspace remain in their physical order, and no same-site collapse fires between them. Declaring the pair distinct here lets the canonical sort resolve their relationship and triggers any composition or completeness rewriting it unlocks.

Examples

julia> h = NLevelSpace(:atom, 2);

julia> j = Index(h, :j, 5, h); k = Index(h, :k, 5, h);

julia> σ(i, m, idx) = IndexedOperator(Transition(h, :σ, i, m), idx);

julia> q = σ(2, 1, k) * σ(1, 2, j);

julia> assume_distinct_index(q, [(j, k)])
σ_j₁₂ * σ_k₂₁

See also expand_completeness.

index_slot(x) -> Union{Int, Nothing}
index_slot(idx::Index) -> Union{Int, Nothing}

Concrete slot integer of a per-slot index, or nothing for an abstract index. On an Index reads the slot field directly; on a symbol reads the IndexSlot metadata minted by (i::Index)(k), so consumers can recover the position k without parsing the symbol's name.

index_range(idx::Index) -> Num

The summation range of idx (the user's Num, recovered from the intern table). Returns Num(0) for the sentinel NO_INDEX.

index_name(idx::Index) -> Symbol

The display name of idx (recovered from the intern table).

index_sym(idx::Index) -> Num

The symbolic variable for idx, reconstructed from its interned name (plus the IndexSlot metadata for a per-slot index). SymbolicUtils hashconsing makes the reconstruction identical (===) to the originally minted symbol, so substitution and get_variables are unaffected. An anonymous concrete site (name_id == 0, slot == k) reconstructs to the integer Num(k), matching to_numeric's resolved-site convention.

evaluate(eqs::AbstractMeanfieldEquations; limits=nothing, h=Int[])
evaluate(c::CorrelationFunction; limits=nothing, h=Int[])

Materialise a symbolic indexed system into a concrete, fixed-size one by inserting the index ranges and unrolling the sums. The result is ready for ModelingToolkitBase.System; its operators carry no remaining sum scope.

Keyword arguments

• limits=nothing: the concrete size of each symbolic index range, given as a Pair (N => 3), a tuple of pairs, or a Dict. Required whenever a range bound is symbolic.
• h=Int[]: restrict unrolling to these Hilbert subspaces (by acts_on index). Empty unrolls every subspace covered by limits.

scale(eqs::AbstractMeanfieldEquations; h=Int[])
scale!(eqs::AbstractMeanfieldEquations; h=Int[])

Reduce a permutation-symmetric indexed system to one representative per symmetry class, exploiting that every atom of an indexed family acts on the system identically. scale! mutates eqs in place; scale returns a new system.

Keyword arguments

• h=Int[]: the Hilbert subspaces to scale, given by their acts_on indices. Empty scales every symmetric (indexed) subspace and leaves the others untouched.

scale(c::CorrelationFunction)

Scale both the τ-equations and the underlying eqs0 so the spectrum's steady-state lookup resolves against the same scaled state names.

scale(eqs::AbstractMeanfieldEquations; h=Int[])
scale!(eqs::AbstractMeanfieldEquations; h=Int[])

Reduce a permutation-symmetric indexed system to one representative per symmetry class, exploiting that every atom of an indexed family acts on the system identically. scale! mutates eqs in place; scale returns a new system.

Keyword arguments

• h=Int[]: the Hilbert subspaces to scale, given by their acts_on indices. Empty scales every symmetric (indexed) subspace and leaves the others untouched.

NoiseMeanfieldEquations

Concrete equation set produced by meanfield when efficiencies is supplied, adding a measurement-backaction noise drift to the deterministic equations. The direction tag (Forward or Backward) drives compile-time dispatch of the noise assembly. Like MeanfieldEquations it is graph-backed: graph is the source of truth and the vector fields are a derived view.

Fields

• graph: the source of truth.
• equations: the deterministic averaged differential equations (derived view).
• noise_equations: the averaged noise (Brownian) drift per state.
• operator_equations, operator_noise_equations: the operator-level counterparts.
• states: the averages on the left-hand sides.
• operators: the operators on the left-hand sides.
• hamiltonian: the system Hamiltonian.
• jumps, jumps_dagger: the collapse operators and their adjoints.
• rates: the decay rates corresponding to jumps.
• efficiencies: the detector efficiencies per jump.
• iv: the independent (time) variable.
• order: the cumulant-expansion order, or nothing.
• direction: Forward or Backward evolution.

translate_W_to_Y(eqs::NoiseMeanfieldEquations; mix_choice=maximum)

Rewrite an SDE whose noise drift is parametrised by the underlying Wiener process dW into one parametrised by the measurement record dY instead. The substitution dW = dY - sqrt(2η)·⟨J + J†⟩·dt adds a deterministic correction to the drift; the noise drift itself is unchanged.

For each equation the RHS is augmented by cumulant_expansion(-_dY_dS_extra_term(lhs_op, J, Jdagger, rates .* efficiencies)), cumulant-expanded to eqs.order when set. Returns a fresh NoiseMeanfieldEquations of the same direction. The augmented RHS is left unsimplified; apply simplify yourself for a canonical form.

EvolutionDirection

Supertype of the singleton tags Forward and Backward used to dispatch the noise/retrodiction code path at compile time.

Forward Heisenberg evolution (positive sign on i[H, ·]).

Backward Heisenberg evolution (negative sign), used for retrodiction.

find_missing(eqs::AbstractMeanfieldEquations; filter_func=nothing)

Return the averages that appear on a right-hand side of eqs but are not yet among its tracked states: the equations still needed to close the system. A moment and its conjugate count as one: if either is already a state, neither is reported. filter_func mirrors complete!.

The report is representation-invariant: conjugate pairs are always folded onto one key, so the answer is the same whether eqs was closed with get_adjoints true or false.

See also: complete, complete!.

find_operators(h::HilbertSpace, order::Int; names=nothing) -> Vector

Generate all unique operator products up to order factors for Hilbert space h.

Starts from fundamental_operators(h) and their adjoints, forms all products with up to order factors, then filters out zero terms and adjoint-duplicates. Useful for constructing the operator basis needed for cumulant expansions or moment equations. Pass names to override auto-generated operator names.

Examples

julia> h = FockSpace(:f);

julia> length(find_operators(h, 1))
1

See also fundamental_operators, unique_ops.

complete(eqs::AbstractMeanfieldEquations; order=nothing, kw...)

Non-mutating complete!: derive equations of motion for every average that appears on a right-hand side but is not yet a state, returning a new closed system.

Keyword arguments

• order=nothing: if given, the copy is first cumulant-expanded to this order (deterministic systems only) before closing.
• kw...: forwarded to complete! (max_iter, filter_func, get_adjoints).

See also: complete!, find_missing, meanfield.

complete!(eqs::AbstractMeanfieldEquations; max_iter=100_000, filter_func=nothing,
          get_adjoints=true)

Close eqs in place by repeatedly deriving equations of motion for every average that appears on a right-hand side but is not yet a state, until the set is self-contained. Closure reads and rewrites eqs.graph, so user RHS modifications and filter substitutions already recorded there survive. Right-hand sides are left unsimplified.

Keyword arguments

• max_iter=100_000: runaway backstop; an error is raised if closure has not converged within this many iterations.
• filter_func=nothing: a predicate filter_func(avg); rejected averages are dropped (substituted by 0), e.g. to discard phase-invariant terms.
• get_adjoints=true: when false, track one representative per conjugate pair, with the partner recovered by conj at code generation. The default integrates each conjugate as its own (redundant) state; pass false for the minimal state set.

See also: complete, find_missing, meanfield.

unique_ops(ops) -> Vector

Return unique operators from ops, treating op and op' (adjoint) as the same degree of freedom. Only the first occurrence of each operator/adjoint pair is kept.

Examples

julia> h = FockSpace(:f);

julia> @qnumbers a::Destroy(h);

julia> length(unique_ops([a, a', a]))
1

See also unique_ops!, fundamental_operators.

simplify!(eqs::AbstractMeanfieldEquations; kwargs...)
simplify(eqs::AbstractMeanfieldEquations; kwargs...)

Run SymbolicUtils.simplify on every RHS in eqs (and on the noise drift RHSs of a NoiseMeanfieldEquations). simplify! mutates; simplify returns a fresh struct. The derivation pipeline leaves expressions raw so the cost of SymbolicUtils.simplify is paid only when explicitly requested.

fundamental_operators(h::HilbertSpace; names=nothing) -> Vector{Op}

Return the minimal generating set of operators for Hilbert space h.

Returns one Destroy per FockSpace; n(n+1)/2 - 1 Transition operators per NLevelSpace (upper-triangular plus diagonals, excluding the ground-state projector); three Pauli or Spin operators per PauliSpace/SpinSpace; one Position and one Momentum per PhaseSpace; and the concatenation of the above for a ProductSpace. Pass names to override the auto-generated operator names.

Examples

julia> h = FockSpace(:cavity) ⊗ NLevelSpace(:atom, 2);

julia> length(fundamental_operators(h))
3

See also find_operators, unique_ops.

to_numeric(op, basis [, d::AbstractDict{<:QSym}])
to_numeric(op, state [, d::AbstractDict{<:QSym}])

Convert a symbolic operator expression to a numeric QuantumOpticsBase operator. d substitutes individual QSyms with custom numeric operators. Throws ArgumentError if any prefactor cannot be reduced to a concrete ComplexF64.

Examples

julia> using QuantumOpticsBase

julia> h = FockSpace(:f);

julia> @qnumbers a::Destroy(h);

julia> b = FockBasis(5); ψ = fockstate(b, 2);

julia> real(QuantumOpticsBase.expect(to_numeric(a' * a, ψ), ψ)) ≈ 2
true

See also numeric_average.

numeric_average(op, state[, d::AbstractDict{<:QSym}]) -> ComplexF64
numeric_average(op, states::AbstractVector[, d::AbstractDict{<:QSym}]) -> Vector

Expectation value $\langle \psi | \hat{O} | \psi \rangle$ of a symbolic operator expression. Averaged BasicSymbolic expressions (from average) are unwrapped automatically. d substitutes symbolic operators with custom numeric operators before evaluation.

Examples

julia> using QuantumOpticsBase

julia> h = FockSpace(:f);

julia> @qnumbers a::Destroy(h);

julia> b = FockBasis(5); ψ = fockstate(b, 2);

julia> real(numeric_average(a' * a, ψ)) ≈ 2
true

See also to_numeric, average.

initial_values(eqs::AbstractMeanfieldEquations; defaults=Dict())

Map every state's u(t) variable to its initial value, defaulting to zero(ComplexF64) for averages absent from defaults.

initial_values(eqs, state)            # from a numeric quantum state
initial_values(eqs, u0::AbstractVector)  # from a value vector aligned with eqs.states

get_solution(sol, avg_or_op, eqs)

Query an ODE/SDE solution for the trajectory of an average (or a QField, averaged internally). Resolves the user's operator to a tracked state through the single structural matcher match_moment, which folds Hermitian conjugation (the side bit recovers a stored conjugate) and the system's index/symmetry treatment, so a query posed in any equivalent symbolic form hits the same state.

parameter_map(sys::MTK.System, pairs)

Keep only the entries of pairs whose key is a live unknown or parameter of the compiled system sys (MTK rejects superfluous entries).

parameter_map(eqs::AbstractMeanfieldEquations, pairs) -> Dict

Translate a user-facing Pair/Dict of parameter assignments into an MTK-compatible map. An IndexedVariable callable key (e.g. g from g(i) = IndexedVariable(:g, i)) is matched by name to the Symbolics array parameter that evaluate synthesised from the per-atom references; a scalar value is broadcast to fill that array, an AbstractArray value passes through as the explicit per-atom values. Scalar (non-array) parameters pass through.

modify_equations(eqs::AbstractMeanfieldEquations, f::Function)

Return a copy of eqs whose RHS for each equation has been rewritten by f. The function is called as f(lhs_op, rhs), receiving the operator form of the LHS (undo_average(eq.lhs)) and the symbolic RHS, and returning a new RHS.

f(lhs, rhs) = rhs + cumulant_expansion(average(commutator(1im * Hadd, lhs)), 2)
eqs_mod = modify_equations(eqs, f)

See also modify_equations!.

modify_equations!(eqs::AbstractMeanfieldEquations, f::Function)

In-place version of modify_equations. Rewrites each drift in eqs.graph with f(lhs_op, rhs), where lhs_op is the node key (the operator form of the LHS).

make_time_dependent(expr, iv) -> expr

Lift every iv-free average leaf in expr into a time-dependent Number variable name(iv) carrying the operator in AverageOperator metadata (and the VariableSource set by @variables, so it reads as a ModelingToolkit unknown). Non-average structure is rebuilt only where a child changed. The lifted node is a leaf: the walk does not descend into iv. The default name is uniqueness-only; downstream code may rename for display without changing identity.