GitHub

Averages and numerics

Measured quantities are expectation values, and the equations of motion of an open or closed system are usually written for them. The Heisenberg equation turns into $\frac{d}{dt}\langle X \rangle = i \langle [H, X] \rangle$, a relation between averages of operator products. Consider a driven Jaynes–Cummings model,

\[H = \Delta a^\dagger a + g \left(a^\dagger \sigma^{ge} + \sigma^{eg} a\right) + \eta \left(a + a^\dagger\right).\]

Averages

average maps an operator expression to its expectation value. Averaging is linear: it distributes over sums and pulls out the coefficients, so the right-hand side of the Heisenberg equation for $a$ becomes a sum of averages.

using SecondQuantizedAlgebrah = FockSpace(:cavity) ⊗ NLevelSpace(:atom, (:g, :e))a = Destroy(h, :a, 1)σ(i, j) = Transition(h, :σ, i, j, 2)@variables Δ g η tH = Δ * a' * a + g * (a' * σ(:g, :e) + σ(:e, :g) * a) + η * (a + a')da = average(1im * commutator(H, a))

\[ \begin{equation} - \mathtt{im} ~ \eta - \langle {\sigma}^{{12}} \rangle ~ g ~ \mathtt{im} - \langle a \rangle ~ \mathtt{im} ~ \Delta \end{equation} \]

An average is an ordinary symbolic scalar, so it takes part in Symbolics arithmetic. undo_average recovers the operator expression.

undo_average(da)

\[ - \eta ~ \mathit{i} - \Delta ~ \mathit{i} a - g ~ \mathit{i} {\sigma}^{{12}}\]

SecondQuantizedAlgebra.average — Function
average(expr)

Expectation value of an operator expression, as a symbolic scalar.

The average is linear: it distributes over sums and takes coefficients out, so the result is a sum of averages of operator products, printed as ⟨…⟩. It takes part in ordinary symbolic arithmetic with Symbolics.jl. The average of a Hermitian operator is real. A sum over an index becomes a sum of averages over the same index.

Examples

julia> a = Destroy(FockSpace(:cavity), :a);julia> @variables ω;julia> average(ω * a' * a + a)⟨a⟩ + ⟨a' * a⟩*ω

See also undo_average, numeric_average, make_time_dependent.

source

Conjugation

The equation for $\langle a^\dagger \rangle$ is the complex conjugate of the one for $\langle a \rangle$. Taken with conj, the conjugate of an average stays formal. inner_adjoint moves it inside, where it becomes the average of the adjoint operator, $\langle X \rangle^* = \langle X^\dagger \rangle$, so the conjugate equation is again written in averages.

inner_adjoint(da)

\[ \begin{equation} \mathtt{im} ~ \eta + \langle {\sigma}^{{21}} \rangle ~ g ~ \mathtt{im} + \langle a^{\dagger} \rangle ~ \mathtt{im} ~ \Delta \end{equation} \]

SecondQuantizedAlgebra.inner_adjoint — Function
inner_adjoint(x)

Complex conjugate of an expression of averages, written again as averages.

The conjugate of $\langle X \rangle$ is written as $\langle X^\dagger \rangle$ rather than as conj(⟨X⟩), so that equations for averages stay expressed in averages of operators. Coefficients are conjugated as by qadjoint.

Examples

julia> a = Destroy(FockSpace(:cavity), :a);julia> @variables g::Number;julia> inner_adjoint(g * average(a))⟨a'⟩*conj(g)

See also qadjoint, average.

source

Time dependence

To integrate the equations, every average becomes an unknown function of time.

make_time_dependent(da, t)

\[ \begin{equation} - \mathtt{im} ~ \eta - \langle a \rangle\left( t \right) ~ \mathtt{im} ~ \Delta - \langle {\sigma}^{{12}} \rangle\left( t \right) ~ g ~ \mathtt{im} \end{equation} \]

SecondQuantizedAlgebra.make_time_dependent — Function
make_time_dependent(expr, t)

Replace every average in expr by an unknown function of time t.

Each average $\langle X \rangle$ becomes $\langle X \rangle(t)$, which can be differentiated with respect to t and serves as an unknown of a differential equation, for instance in ModelingToolkit.jl. The unknown is real when $X$ is Hermitian and complex otherwise. undo_average recovers the operators.

Examples

julia> a = Destroy(FockSpace(:cavity), :a);julia> @variables t;julia> make_time_dependent(average(a' * a) + average(a), t)⟨a⟩(t) + ⟨a' * a⟩(t)

See also average, undo_average.

source

Operator bases

The equation for $\langle a \rangle$ involves $\langle \sigma^{ge} \rangle$, whose own equation brings in products such as $\langle a \sigma^{ee} \rangle$. Closing this hierarchy starts from a set of operators to track. fundamental_operators gives the generators of the algebra on a space, and find_operators all products up to a given order, each counted once together with its adjoint.

fundamental_operators(h)

\[ \left[ \begin{array}{c} a \\ {\sigma}^{{12}} \\ {\sigma}^{{22}} \\ \end{array} \right] \]

find_operators(h, 2)

\[ \left[ \begin{array}{c} a \\ {\sigma}^{{12}} \\ {\sigma}^{{22}} \\ aa \\ a{\sigma}^{{12}} \\ a{\sigma}^{{22}} \\ a^{\dagger}a \\ a{\sigma}^{{21}} \\ \end{array} \right] \]

unique_up_to_adjoint applies the same reduction to a list of your own.

unique_up_to_adjoint([a, a', a' * a, σ(:g, :e), σ(:e, :g)])

\[ \left[ \begin{array}{c} a \\ a^{\dagger}a \\ {\sigma}^{{12}} \\ \end{array} \right] \]

SecondQuantizedAlgebra.fundamental_operators — Function
fundamental_operators(h::HilbertSpace; names = Symbol[])

Smallest set of operators from which every operator on h is built by products and adjoints.

A bosonic mode contributes its annihilation operator; a finite-level system its transitions $\sigma^{ij}$ with $i < j$ and its projectors except the ground state; a CollectiveNLevelSpace all of its collective transitions; Pauli and spin systems their three components; and a pair of quadratures position and momentum. On a product space the sets of the factors are joined. names sets the operator names, one per factor; otherwise they are generated.

Examples

julia> h = FockSpace(:cavity) ⊗ NLevelSpace(:atom, 3);julia> fundamental_operators(h; names = [:b, :s])6-element Vector{Op}: b s₁₂ s₁₃ s₂₂ s₂₃ s₃₃

See also find_operators, unique_up_to_adjoint.

source
SecondQuantizedAlgebra.find_operators — Function
find_operators(h::HilbertSpace, order::Int; names = Symbol[])

All distinct products of up to order operators on h, counting adjoints only once.

The products are built from fundamental_operators and their adjoints, brought into canonical form, and filtered with unique_up_to_adjoint. They are the averages a set of moment or cumulant equations of that order is written in. names is passed on to fundamental_operators.

Examples

julia> find_operators(FockSpace(:cavity), 2)3-element Vector{Union{Op, SecondQuantizedAlgebra.QAdd}}: a a * a a' * a

See also fundamental_operators, unique_up_to_adjoint.

source
SecondQuantizedAlgebra.unique_up_to_adjoint — Function
unique_up_to_adjoint(ops)

Operators of ops with duplicates removed, counting an operator and its adjoint as the same.

The first of each pair is kept. Equations for averages need only one of $\langle X \rangle$ and $\langle X^\dagger \rangle$, since the other is its complex conjugate.

Examples

julia> a = Destroy(FockSpace(:cavity), :a);julia> unique_up_to_adjoint([a, a', a' * a])2-element Vector{SecondQuantizedAlgebra.QField}: a a' * a

See also find_operators, fundamental_operators.

source

Many emitters

With $N$ atoms in the cavity, each with its own coupling $g_i$, the interaction becomes a Σ over the atoms. The average of a sum is a sum of averages over the same index, with the same excluded sites, so the equations keep $N$ symbolic.

using SecondQuantizedAlgebrah = FockSpace(:cavity) ⊗ NLevelSpace(:atom, 2)a = Destroy(h, :a, 1)@variables Ni = Index(h, :i, N, 2)σ(k, l) = IndexedOperator(Transition(h, :σ, k, l, 2), i)g = IndexedVariable(:g, i)H = Σ(g * (a' * σ(1, 2) + a * σ(2, 1)), i)average(H)

\[ \begin{equation} \underset{i}{\overset{N}{\sum}}\left( \langle a{\sigma}_{i}^{{21}} \rangle ~ g\left( i \right) + \langle a^{\dagger}{\sigma}_{i}^{{12}} \rangle ~ g\left( i \right) \right) \end{equation} \]

Numerical conversion

to_numeric turns an operator expression into a matrix of QuantumOpticsBase.jl or QuantumToolbox.jl, ready for their solvers. The numerical side is an extension: it becomes available once one of the two packages is loaded, and the backend is chosen with QuantumOpticsBackend or QuantumToolboxBackend. The portable form takes the symbolic Hilbert space and one truncation per factor, a photon cutoff for a mode. Every parameter needs a value.

using QuantumOpticsBasebackend = QuantumOpticsBackend()parameter = Dict(Δ => 0.5, g => 0.1, η => 0.2)Hn = to_numeric(H, h, (10, 2); backend, parameter)size(Hn)
(22, 22)
SecondQuantizedAlgebra.to_numeric — Function
to_numeric(op, h::HilbertSpace, dims; backend, parameter, time_parameter, operators, adjoint_ops, op_type)
to_numeric(op, basis; kwargs...)
to_numeric(op, state; kwargs...)
to_numeric(op, basis_or_state, operators)
to_numeric(ops::AbstractVector, args...; kwargs...)

Matrix of an operator expression on a truncated basis, for a numerical backend.

The first form builds the basis from the symbolic space h and works with every backend. dims holds one entry per factor: the photon cutoff for a mode, so 5 gives six states, and the spin length for a spin. The backend is QuantumOpticsBackend() or QuantumToolboxBackend(), and may be left out when only one of the two packages is loaded. The other forms take a basis of the backend, or a state whose basis is used; QuantumToolbox also accepts its own dimension tuple, which counts states rather than photons. The third argument operators replaces operators by given matrices, as the keyword does with adjoint_ops = false. A vector of expressions gives a vector of matrices.

Keywords:

  • parameter: values of the symbolic parameters. A coefficient that does not reduce to a number throws an ArgumentError.
  • time_parameter: values that depend on time, as functions t -> value or numbers; the time variable itself is given as t => identity. A non-empty time_parameter makes the result time dependent, a TimeDependentSum for QuantumOptics and a QobjEvo for QuantumToolbox. With QuantumToolbox, a function (p, t) -> value of the solver parameters p makes the QobjEvo differentiable in p, for optimal control, and may be mixed with the other entries.
  • operators: matrices on the full basis to use for given operators in place of the standard ones.
  • adjoint_ops = true: also replace the adjoint of each operator in operators by the adjoint matrix.
  • op_type: the default gives a sparse matrix. identity keeps the backend's lazy representation, which saves memory on large product spaces, and another function, such as dense or QuantumToolbox.to_dense, is applied to the result. A time-dependent result is always the backend's time-dependent operator, so only the default and identity are accepted.

Formal operator functions throw an ArgumentError; expand them with taylor first.

Examples

julia> using QuantumOpticsBasejulia> h = FockSpace(:cavity) ⊗ NLevelSpace(:atom, 2);julia> @qnumbers a::Destroy(h, 1) σ::Transition(h, 1, 2, 2);julia> @variables g;julia> Hn = to_numeric(g * (a' * σ + σ' * a), h, (5, 2); backend = QuantumOpticsBackend(), parameter = Dict(g => 0.1));julia> size(Hn)(12, 12)

See also numeric_average, QuantumOpticsBackend, QuantumToolboxBackend.

source
SecondQuantizedAlgebra.QuantumOpticsBackend — Type
QuantumOpticsBackend()

Numerical backend built on QuantumOpticsBase.jl.

Available after using QuantumOpticsBase. Operators become Operators on QuantumOptics bases, and time-dependent ones TimeDependentSums, ready for the solvers of QuantumOptics.jl.

Examples

julia> using QuantumOpticsBasejulia> h = FockSpace(:cavity); a = Destroy(h, :a);julia> to_numeric(a' * a, h, 4; backend = QuantumOpticsBackend()) ≈ number(FockBasis(4))true

See also QuantumToolboxBackend, to_numeric.

source

Bases and levels

A basis of the backend can be passed in place of the space and its truncation, and then also fixes the backend. Symbolic level names map to basis indices in the order of the levels of the NLevelSpace, so $g$ is the first basis state and $e$ the second.

bf = FockBasis(10)bn = NLevelBasis(2)b = bf ⊗ bnto_numeric(σ(:g, :e), b) == embed(b, 2, transition(bn, 1, 2))
true

Time-dependent parameters

Parameters that change in time go into time_parameter, as functions of time. The result is then the backend's time-dependent operator, a TimeDependentSum for QuantumOptics.

Ht = to_numeric(H, h, (10, 2); backend, parameter = Dict(Δ => 0.5, g => 0.1), time_parameter = Dict(η => t -> 0.2 * sin(t)))Ht isa TimeDependentSum
true

An explicit time dependence of the coefficients, such as a phase from a rotating frame, needs the time variable itself, given as t => identity. A constant may stand in time_parameter next to it.

Hd = η * (expim(Δ * t) * a + expim(-Δ * t) * a')to_numeric(Hd, h, (10, 2); backend, parameter = Dict(Δ => 0.5), time_parameter = Dict(t => identity, η => 0.2)) isa TimeDependentSum
true

Custom operators and lazy output

operators replaces the matrix of an operator by one of your choice, given on the full basis, and the adjoint is replaced along with it. Replacing $a$ by $a + \alpha$ displaces the mode numerically.

α = 0.5a_displaced = embed(b, 1, destroy(bf)) + α * identityoperator(b)to_numeric(a' * a, b; operators = Dict(a => a_displaced)) ≈ to_numeric((a' + α) * (a + α), b)
true

The default result is a sparse matrix. On a large product space, op_type = identity keeps the backend's lazy sum of tensor products and never forms the full matrix.

to_numeric(H, b; parameter, op_type = identity) isa LazySum
true

Indexed systems

A sum over an index is written out over explicit copies of the subsystem, so the basis fixes the number of atoms. It has one factor per copy, and sites maps the position of each subsystem in the symbolic space to the positions of its copies. Parameters that depend on the site are keyed by site, obtained by substituting a number for the index.

using QuantumOpticsBaseb = FockBasis(3) ⊗ NLevelBasis(2) ⊗ NLevelBasis(2)couplings = Dict(substitute(g, Dict(index_sym(i) => k)) => 0.1k for k in 1:2)Hn = to_numeric(H, b, Dict(1 => [1], 2 => [2, 3]); parameter = couplings)size(Hn)
(16, 16)
SecondQuantizedAlgebra.to_numeric — Method
to_numeric(expr, basis, sites::Dict{Int, Vector{Int}}; parameter, operators)
to_numeric(expr, basis, sites::Dict{Int, Vector{Int}}[, operators[, parameter]])

Matrix of an expression with sums over indices, on a basis with one factor per site.

Each sum is written out over the sites. sites maps the position of each subsystem in the symbolic space to the positions of its copies in basis; a subsystem without an index has a single copy. Pairs of indices that must differ skip their coinciding sites. parameter gives the values of all parameters. A site-dependent parameter $g_i$ is keyed by site, with the value of $g_1$ under substitute(g, Dict(index_sym(i) => 1)). operators maps operators to matrices as for the other forms. This form is available for QuantumOpticsBase only.

Examples

julia> using QuantumOpticsBasejulia> h = FockSpace(:cavity) ⊗ NLevelSpace(:atom, 2);julia> a = Destroy(h, :a, 1); i = Index(h, :i, 2, 2);julia> σ(k, l) = IndexedOperator(Transition(h, :σ, k, l, 2), i);julia> g = IndexedVariable(:g, i); H = Σ(g * (a' * σ(1, 2) + a * σ(2, 1)), i);julia> b = FockBasis(3) ⊗ NLevelBasis(2) ⊗ NLevelBasis(2);julia> gs = Dict(substitute(g, Dict(index_sym(i) => k)) => 0.1k for k in 1:2);julia> size(to_numeric(H, b, Dict(1 => [1], 2 => [2, 3]); parameter = gs))(16, 16)

See also Σ, index_sym.

source

Expectation values

numeric_average evaluates an operator expression, or an expression of averages, in a state. For a coherent state with amplitude $\alpha$ and the atom in its ground state, $\langle \sigma^{ge} \rangle$ vanishes and the right-hand side of the equation for $\langle a \rangle$ is $-i(\Delta\alpha + \eta) = -0.45i$.

ψ = coherentstate(bf, α) ⊗ nlevelstate(bn, 1)numeric_average(substitute(da, parameter), ψ)
0.0 - 0.4499999999999477im

A vector of states, such as the output of a time evolution, gives a vector of values. expect is the same for operator expressions, under the name used by QuantumOptics.jl and QuantumToolbox.jl.

ψs = [coherentstate(bf, β) ⊗ nlevelstate(bn, 1) for β in (0.1, 0.5)]numeric_average(a' * a, ψs)
2-element Vector{ComplexF64}:
 0.009999999999999998 + 0.0im
  0.24999999999994768 + 0.0im
expect(a' * a, ψ)
0.24999999999994768 + 0.0im

On a large product space, a LazyKet holds the state of each factor separately and avoids the full state vector.

ψl = LazyKet(b, (coherentstate(bf, 0.5), nlevelstate(bn, 2)))numeric_average(average(a' * a + σ(:e, :e)), ψl)
1.2499999999999432 + 0.0im
SecondQuantizedAlgebra.numeric_average — Function
numeric_average(x, state[, operators])
numeric_average(x, states::AbstractVector[, operators])

Expectation value of an operator expression or of an expression of averages in a numerical state.

x is either an operator expression or a result of average, with every parameter already replaced by a number. The state is a ket or a density matrix of the backend, and its basis fixes the truncation; a vector of states gives a vector of values. operators maps operators to the matrices to use in place of the standard ones.

Examples

julia> using QuantumOpticsBasejulia> a = Destroy(FockSpace(:cavity), :a);julia> ψ = coherentstate(FockBasis(20), 0.5);julia> numeric_average(a' * a, ψ) ≈ 0.25truejulia> numeric_average(average(a) + average(a'), ψ) ≈ 1true

See also expect, to_numeric, average.

source
QuantumInterface.expect — Method
expect(op, state[, operators])

Expectation value of an operator expression, the same as numeric_average.

The name matches expect of QuantumOptics.jl and QuantumToolbox.jl. Unlike numeric_average, it does not accept expressions of averages.

Examples

julia> using QuantumOpticsBasejulia> a = Destroy(FockSpace(:cavity), :a);julia> expect(a' * a, fockstate(FockBasis(5), 2)) ≈ 2true

See also numeric_average.

source