Implementation

This page walks through the main concepts in SecondQuantizedAlgebra.jl: defining Hilbert spaces, creating operators, building algebraic expressions, and converting to numerics.

Hilbert spaces

Every system starts with a Hilbert space. The package provides five concrete spaces:

hf = FockSpace(:cavity)

NLevelSpace requires the number of levels. Levels can be integers or symbolic names, and a ground state can be specified (defaults to the first level):

ha = NLevelSpace(:atom, 2)                  # levels 1, 2; ground state = 1
ha_sym = NLevelSpace(:atom, (:g, :e))       # symbolic levels; ground state = :g
ha_gs = NLevelSpace(:atom, 4, 2)            # 4 levels; ground state = 2

The ground-state projector $|g\rangle\langle g|$ is eliminated via completeness $\sum_j |j\rangle\langle j| = 1$. Same-site composition that produces one (e.g. $\sigma^{12} \cdot \sigma^{21}$) triggers this eagerly inside *. User-constructed ground-state projectors (e.g. Transition(h, :σ, 1, 1) typed directly, never multiplied) stay atomic until they enter a * or until you call expand_completeness explicitly.

Composite systems are built with ⊗ (\otimes<tab>) or tensor:

h = hf ⊗ ha
h3 = hf ⊗ ha ⊗ NLevelSpace(:spin, 3)

Operators

Operators are the symbolic building blocks. Each operator type lives on a specific Hilbert space:

All operators are subtypes of QSym. Here are some examples:

hf = FockSpace(:cavity)
a = Destroy(hf, :a)
a'   # Create — the adjoint of Destroy

ha = NLevelSpace(:atom, (:g, :e))
σge = Transition(ha, :σ, :g, :e)   # |g⟩⟨e|

hp = PauliSpace(:qubit)
σx = Pauli(hp, :σ, 1)   # axis: 1=x, 2=y, 3=z

In composite systems, specify which subspace the operator acts on by index:

h = FockSpace(:c1) ⊗ FockSpace(:c2) ⊗ NLevelSpace(:atom, 2)
a = Destroy(h, :a, 1)             # acts on first FockSpace
b = Destroy(h, :b, 2)             # acts on second FockSpace
σ12 = Transition(h, :σ, 1, 2, 3)  # acts on the NLevelSpace

For convenience, the @qnumbers macro creates operators in one line:

h = FockSpace(:cavity) ⊗ NLevelSpace(:atom, 2, 1)
@qnumbers a::Destroy(h, 1)
σ(i, j) = Transition(h, :σ, i, j, 2)

Symbolic parameters

Symbolic parameters (c-numbers) are created using @variables from Symbolics.jl, which is re-exported by SecondQuantizedAlgebra. Variables are real by default. The symtype annotation controls how conjugation is handled:

h = FockSpace(:cavity)
@qnumbers a::Destroy(h)
@variables ω η::Number      # ω is real, η is complex

H = ω * a' * a + η * a + conj(η) * a'

Real variables satisfy conj(ω) == ω, which simplifies adjoint expressions:

conj(ω)

\[ \begin{equation} \omega \end{equation} \]

A ::Number parameter is complex-valued but stays a single symbol, so its conjugate is the symbolic conj(η):

conj(η)

\[ \begin{equation} \mathrm{conj}\left( \eta \right) \end{equation} \]

Use ::Number for a complex parameter you want kept atomic (e.g. a coupling amplitude); use ::Complex when you want the parameter decomposed into independent real and imaginary unknowns. ::Number keeps coefficient arithmetic on a single symbol, while ::Complex carries both parts through every product.

Algebraic expressions and commutation relations

All operator expressions are stored as QAdd, a dictionary mapping operator sequences to prefactors. Commutation rules are applied eagerly at construction time: every * immediately normal-orders the result, applies algebraic identities (Transition composition, Pauli products), and expands ground-state completeness. The dict key always reflects the canonical form.

a * a'   # immediately gives a†a + 1

\[1 + a^{\dagger}a\]

normal_order can be called explicitly as an idempotent finalizer (useful in tests and for hand-constructed expressions that bypass *):

normal_order(a * a')   # a†a + 1

\[1 + a^{\dagger}a\]

simplify runs normal_order first, then walks the resulting terms applying Symbolics.simplify to each coefficient and dropping any summation indices no surviving term depends on. Since eager * already produces canonical form, calling simplify on the output of * is typically idempotent at the operator level, but it can still reduce symbolic coefficients (e.g. collapsing sin²(ω) + cos²(ω) to 1).

Substitution

Use substitute for both scalar parameters and operator leaves. Mixed rule dictionaries are split by key: QSym keys replace operators, while all other keys are applied to coefficients through SymbolicUtils.

Operator substitutions are simultaneous and single-pass. The replacement expression is not searched again for operator keys, so transformations such as a => g*a + h*b are well-defined even though the right-hand side contains a. Missing adjoint rules are generated by default; pass replace_adjoint=false to match only the keys supplied explicitly.

h = FockSpace(:cavity)
@qnumbers a::Destroy(h)
@variables Δ g

substitute(Δ * a, Dict(a => g * a + 1, Δ => 2.0))

\[2 + 2 ~ g a\]

Averaging

In many-body theory, equations of motion are typically written for expectation values rather than operators. The average function converts an operator expression into a symbolic scalar representing its expectation value $\langle \hat{O} \rangle$:

average(a)

\[ \begin{equation} \langle a \rangle \end{equation} \]

Averaging is linear — it distributes over sums and pulls out c-number prefactors:

average(ω * a' * a + a + a')

\[ \begin{equation} \langle a^{\dagger} \rangle + \langle a \rangle + \langle a^{\dagger}a \rangle ~ \omega \end{equation} \]

The result is a BasicSymbolic (from SymbolicUtils.jl) that participates in standard symbolic arithmetic. To recover the underlying operator expression, use undo_average:

avg = average(a' * a)
undo_average(avg)

\[a^{\dagger}a\]

Averaged expressions also preserve summation metadata when working with indexed operators, so that symbolic sums carry through the averaging process (see Symbolic Sums and Indices).

Numeric conversion

Operators can be converted to numeric representations using QuantumOpticsBase.jl.

Direct conversion

Use to_numeric to convert an operator expression to a numeric operator:

using SecondQuantizedAlgebra, QuantumOpticsBase
h = FockSpace(:cavity)
@qnumbers a::Destroy(h)

b = FockBasis(7)
a_num = to_numeric(a, b)
@assert a_num == destroy(b)

Scalar parameters and custom numeric operators can be supplied directly at the numeric boundary:

@variables Δ::Real
H = Δ * a' * a
H_num = to_numeric(H, b; parameter = Dict(Δ => 2.0))
@assert H_num == 2.0 * create(b) * destroy(b)

For time-dependent parameters, pass numbers or functions in time_parameter. The result is a callable t -> op(t).

Numeric averages

Use numeric_average to compute expectation values for a given state:

α = 0.3 + 0.1im
ψ = coherentstate(b, α)
numeric_average(a' * a, ψ)

0.09999999999818197 - 4.265150033476549e-25im

For symbolic scalar expressions such as average(a), call numeric_average directly. The expect alias is intentionally kept for operator expressions (QField) only.

Symbolic levels for NLevelSpace

Symbolic level names are resolved to integer basis indices at Transition construction time, using the order of the levels tuple passed to NLevelSpace. The first level maps to basis index 1, the second to 2, and so on:

using SecondQuantizedAlgebra, QuantumOpticsBase
h = NLevelSpace(:atom, (:g, :e))
bn = NLevelBasis(2)
s = Transition(h, :s, :g, :e)
@assert to_numeric(s, bn) == transition(bn, 1, 2)

Composite systems

For product spaces, to_numeric returns LazyTensor operators, which efficiently handle large tensor products:

using SecondQuantizedAlgebra, QuantumOpticsBase
hc = FockSpace(:cavity)
ha = NLevelSpace(:atom, (:a, :b, :c))
h = hc ⊗ ha

@qnumbers a::Destroy(h, 1)
s(i, j) = Transition(h, :s, i, j, 2)

bf = FockBasis(10)
bn = NLevelBasis(3)
bc = bf ⊗ bn

a_num = to_numeric(a, bc)
s_num = to_numeric(s(:a, :c), bc)

For very large systems, use LazyKet to avoid materializing the full state vector:

ψ = LazyKet(bc, (coherentstate(bf, 0.3), (nlevelstate(bn, 1) + nlevelstate(bn, 3)) / sqrt(2)))
avg = average(a' * s(:a, :c))
numeric_average(avg, ψ)

0.15000000000000002 + 0.0im