Open quantum systems

A cavity that leaks photons or an atom that emits into free space does not evolve unitarily. For a memoryless environment the density matrix obeys a master equation $\dot\rho = \mathcal{L}\rho$, whose generator has the Gorini–Kossakowski–Sudarshan–Lindblad form

\[\mathcal{L}\rho = -i[H, \rho] + \sum \gamma\, \mathcal{D}[J, K]\rho, \qquad \mathcal{D}[J, K]\rho = J \rho K^\dagger - \tfrac12 \left\{K^\dagger J, \rho\right\},\]

with $\mathcal{D}[J] = \mathcal{D}[J, J]$ for an independent decay channel. Consider a two-level atom in a driven cavity,

\[H = \omega_c a^\dagger a + \omega_a \sigma^{ee} + g\left(a^\dagger \sigma^{ge} + \sigma^{eg} a\right) + \eta\left(e^{-i\omega_d t} a^\dagger + e^{i\omega_d t} a\right),\]

with cavity loss $\kappa$ and spontaneous emission $\gamma$.

Lindblad generators
Rotating frame
Equations of motion
Numerics
Displacement frames
Superoperators
GKSL form of a map
Correlated and collective decay
Formal operator functions
Inspecting and building maps
Secular projection

Lindblad generators

hamiltonian_action builds the coherent part $-i[H, \cdot]$ and dissipator a decay channel, with the rate as a scalar prefactor. The generator keeps its structure and prints as the Hamiltonian and its channels.

using SecondQuantizedAlgebrah = FockSpace(:cavity) ⊗ NLevelSpace(:atom, (:g, :e))a = Destroy(h, :a, 1)σ(i, j) = Transition(h, :σ, i, j, 2)@variables t::Real ωc::Real ωa::Real ωd::Real g::Real η::Real κ::Real γ::RealH = ωc * a' * a + ωa * σ(:e, :e) + g * (a' * σ(:g, :e) + σ(:e, :g) * a) +    η * (expim(-ωd * t) * a' + expim(ωd * t) * a)L = hamiltonian_action(H) + κ * dissipator(a) + γ * dissipator(σ(:g, :e))

\[-i\left[e^{i t ~ \mathtt{{\omega}d}} ~ \eta a + e^{-i t ~ \mathtt{{\omega}d}} ~ \eta a^{\dagger} + \mathtt{{\omega}a} {\sigma}^{{22}} + g a{\sigma}^{{21}} + \mathtt{{\omega}c} a^{\dagger}a + g a^{\dagger}{\sigma}^{{12}},\cdot\right] + \kappa\,\mathcal{D}\left[a\right] + \gamma\,\mathcal{D}\left[{\sigma}^{{12}}\right]\]

SecondQuantizedAlgebra.Lindbladian — Type

Lindbladian{O}

A Lindblad generator that keeps its Hamiltonian and decay channels.

The generator is L = -i[H, ·] + Σ γ 𝒟[J, K]. O is QAdd for a polynomial model and QExpr when the Hamiltonian or a jump contains a formal operator function. Build a generator with hamiltonian_action, dissipator and liouvillian, and combine generators with +, scalar * and Σ. L(ρ) applies the generator, and transpose(L)(A) is the Heisenberg drift of the observable A.

A Lindbladian is a form, not a claim of physicality: H may be non-Hermitian and the rates complex or negative. It always preserves the trace. It preserves Hermiticity when H and the rate matrix are Hermitian, and it generates a completely positive evolution when the rate matrix is also positive semidefinite.

Equality is structural. Equal spellings of one channel compare equal, but a cross term is not expanded, so two generators that act alike are compared with SuperOperator(L₁) == SuperOperator(L₂).

Examples

julia> h = FockSpace(:c); a = Destroy(h, :a);julia> L = hamiltonian_action(a' * a) + 2 * dissipator(a)-i[a' * a, ·] + 2 * 𝒟[a]julia> transpose(L)(a)(-1 + -1im) * a

source

SecondQuantizedAlgebra.hamiltonian_action — Function

hamiltonian_action(H)

The coherent part of a Lindblad generator with Hamiltonian H.

It acts as ρ ↦ -i[H, ρ]. The identity part of H is dropped, since it commutes with every state. A summed H keeps its summation scope.

Examples

julia> h = FockSpace(:c); a = Destroy(h, :a);julia> hamiltonian_action(a' * a + 1)-i[a' * a, ·]

Rotating frame

The drive oscillates at $\omega_d$. transform moves the generator into the frame rotating with $\omega_d (a^\dagger a + \sigma^{ee})$: the Hamiltonian is transformed and every jump is conjugated. The drive becomes static, the frequencies become the detunings $\Delta_c = \omega_c - \omega_d$ and $\Delta_a = \omega_a - \omega_d$, and the jumps are unchanged, since a loss channel only picks up a phase.

Lrot = transform(L, RotatingFrame(ωd * (a' * a + σ(:e, :e)), t))

\[-i\left[\eta a + \eta a^{\dagger} + \left(\mathtt{{\omega}a} - \mathtt{{\omega}d}\right) {\sigma}^{{22}} + g a{\sigma}^{{21}} + \left(\mathtt{{\omega}c} - \mathtt{{\omega}d}\right) a^{\dagger}a + g a^{\dagger}{\sigma}^{{12}},\cdot\right] + \kappa\,\mathcal{D}\left[a\right] + \gamma\,\mathcal{D}\left[{\sigma}^{{12}}\right]\]

Equations of motion

The transpose of a generator is its Heisenberg-picture form, $\operatorname{Tr}(A\,\mathcal{L}\rho) = \operatorname{Tr}(\mathcal{L}^T\!A\,\rho)$, so $\frac{d}{dt}\langle A \rangle = \langle \mathcal{L}^T A \rangle$. The field is damped at $\kappa/2$, driven by $\eta$ and coupled to the atomic coherence.

average(transpose(Lrot)(a))

\[ \begin{equation} - \frac{1}{2} ~ \langle a \rangle ~ \kappa - \mathtt{im} ~ \eta - \langle {\sigma}^{{12}} \rangle ~ g ~ \mathtt{im} + \langle a \rangle ~ \mathtt{im} ~ \left( - \mathtt{{\omega}c} + \mathtt{{\omega}d} \right) \end{equation} \]

The coherence decays at $\gamma/2$ and is driven by the field times the inversion.

average(transpose(Lrot)(σ(:g, :e)))

\[ \begin{equation} - \frac{1}{2} ~ \langle {\sigma}^{{12}} \rangle ~ \gamma + \langle a{\sigma}^{{22}} \rangle ~ g ~ \mathtt{im} - \langle a{\sigma}^{{11}} \rangle ~ g ~ \mathtt{im} + \langle {\sigma}^{{12}} \rangle ~ \mathtt{im} ~ \left( - \mathtt{{\omega}a} + \mathtt{{\omega}d} \right) \end{equation} \]

For a weak drive the atom stays near its ground state, $\langle a\sigma^{gg}\rangle \approx \langle a \rangle$ and $\langle a\sigma^{ee}\rangle \approx 0$, and the two equations close. Their steady state is the vacuum-Rabi response of the cavity field,

\[\langle a \rangle = \frac{-i\eta}{\kappa/2 + i\Delta_c + g^2/(\gamma/2 + i\Delta_a)}.\]

Trace preservation reads $\mathcal{L}^T 1 = 0$.

transpose(Lrot)(1)

\[0\]

Base.transpose — Method

transpose(S::SuperOperator)

The Heisenberg-picture map of a map on operators.

It satisfies Tr(A * S(ρ)) == Tr(transpose(S)(A) * ρ), and turns each term c * L ρ R into c * R A L. For a generator L, transpose(L)(A) is the drift of the observable A, and average(transpose(L)(A)) its equation of motion.

Examples

julia> h = FockSpace(:c); a = Destroy(h, :a);julia> transpose(2 * dissipator(a))(a' * a)-2 * a' * a

Numerics

to_numeric converts a generator to the backend's superoperator, with the same parameter, time_parameter and indexed sites arguments as for operators. The steady state of the rotating-frame Liouvillian at weak drive reproduces the vacuum-Rabi amplitude.

import QuantumToolbox as QTp = Dict(ωc => 1.0, ωa => 1.1, ωd => 0.95, g => 0.2, η => 0.002, κ => 0.1, γ => 0.05)backend = QuantumToolboxBackend()Lnum = to_numeric(Lrot, h, (10, 2); backend, parameter = p)anum = to_numeric(a, h, (10, 2); backend)Δc, Δa = 0.05, 0.15QT.expect(anum, QT.steadystate(Lnum)), -0.002im / (0.05 + Δc * im + 0.2^2 / (0.025 + Δa * im))

(0.0079710666045887 - 0.003550538017207381im, 0.007969151670951156 - 0.003547557840616967im)

Solvers that take a Hamiltonian and collapse operators get them from to_numeric_master, here $\sqrt\kappa\, a$ and $\sqrt\gamma\, \sigma^{ge}$.

using SecondQuantizedAlgebra: to_numeric_masterHnum, c_ops = to_numeric_master(Lrot, h, (10, 2); backend, parameter = p)QT.expect(anum, QT.steadystate(Hnum, c_ops))

0.0079710666045887 - 0.003550538017207381im

SecondQuantizedAlgebra.to_numeric_master — Function

to_numeric_master(L::Lindbladian, h, dims; backend, parameter, time_parameter, sites, scalar_subs, tol)

The numerical Hamiltonian and collapse operators of a Lindbladian, as inputs for a master-equation solver.

Returns (hamiltonian = H, collapse_operators = c_ops), which feed QuantumOptics timeevolution.master and QuantumToolbox mesolve without separate rates. backend, parameter, time_parameter and sites work as for to_numeric, and scalar_subs resolves coefficients of unrolled indices, such as g(k). A single channel γ 𝒟[J] gives √γ J. The rate matrix of every block of correlated channels is diagonalized, c_k = √λ_k Σ_m U[m, k] M_m, where M_m are the operator monomials of the jumps.

The generator must be physical: the rate matrices must be Hermitian and positive semidefinite and H Hermitian, or the conversion throws. Eigenvalues with |λ| ≤ tol * max|λ| are set to zero, with tol = √eps by default. A time-dependent H is checked at a few sample times.

With a time_parameter, which needs QuantumToolbox, H becomes a QobjEvo, and a single channel with a rate γ(t) gives the collapse operator √γ(t) J; the rate must stay real and nonnegative during the solve. A correlated block with time-dependent rates throws, since the eigenvectors of Γ(t) move with t; to_numeric(L, h, dims; time_parameter) gives the full generator instead.

Examples

julia> using SecondQuantizedAlgebra: to_numeric_masterjulia> import QuantumToolboxjulia> h = FockSpace(:c); a = Destroy(h, :a);julia> L = hamiltonian_action(a' * a) + 4 * dissipator(a);julia> H, c_ops = to_numeric_master(L, h, 4; backend = QuantumToolboxBackend());julia> only(c_ops) ≈ 2 * to_numeric(a, h, 4; backend = QuantumToolboxBackend())true

Displacement frames

A damped linear mode is displaced to its steady state. DisplacementFrame takes the amplitude from the fixed point of the full drift, so loss and thermal heating both enter it, and the drive is gone in the displaced frame.

using SecondQuantizedAlgebrab = Destroy(FockSpace(:c), :b)Lth = hamiltonian_action(2 * b' * b + 3 * (b + b')) + (1 // 2) * dissipator(b) + (1 // 10) * dissipator(b')U = DisplacementFrame(b, Lth)

UnitaryTransform with 2 rules:
  b ↦ (-150//101 + -15//101im) + b
  b' ↦ (-150//101 + 15//101im) + b'

transpose(transform(Lth, U))(b)

\[\left(\frac{-1}{5}-2\mathit{i}\right) b\]

Superoperators

A Lindbladian is one kind of linear map on operators. The general polynomial map $\rho \mapsto \sum c\, A \rho B$ is a SuperOperator, built from the sandwich sprepost and its one-sided forms spre and spost. SuperOperator(L) lowers a Lindbladian to its sandwiches.

using SecondQuantizedAlgebra: SuperOperator, spre, spost, sprepostD = SuperOperator(dissipator(a))

\[-\frac{1}{2} \cdot\,a^{\dagger} a + a\,\cdot\,a^{\dagger} - \frac{1}{2} a^{\dagger} a\,\cdot\]

D == sprepost(a, a') - (1 // 2) * spre(a' * a) - (1 // 2) * spost(a' * a)

true

These names are public but not exported, because QuantumOpticsBase and QuantumToolbox export the same names for their numerical superoperators. With a backend loaded, its spre on a symbolic operator forwards to the symbolic one.

Maps compose with ∘ or *, and the right factor acts first. Composing, transposing or applying a Lindbladian gives a SuperOperator.

(spre(a) ∘ spost(a'))(σ(:e, :e))

\[{\sigma}^{{22}} + a^{\dagger}a{\sigma}^{{22}}\]

Every Lindbladian preserves the trace. The identity map does not, so it is not a Lindbladian: L^0, like a Lindbladian plus a scalar, is a general map. It prints as the bare placeholder of $\rho \mapsto \rho$.

L^0

\[\cdot\]

transpose and adjoint are different maps. The adjoint belongs to the Hilbert–Schmidt product $\operatorname{Tr}(X^\dagger S\rho)$ and takes the adjoint of both sides, so the two agree only for maps that preserve Hermiticity.

transpose(spre(a))(σ(:e, :e)), adjoint(spre(a))(σ(:e, :e))

(a * σ₂₂, a' * σ₂₂)

SecondQuantizedAlgebra.SuperOperator — Type

SuperOperator
SuperOperator(L::Lindbladian)
SuperOperator(pairs, indices)

A polynomial linear map on operators.

The map acts as ρ ↦ Σ c * A * ρ * B. Build maps with sprepost, spre and spost, or lower a Lindbladian with SuperOperator(L). SuperOperator(pairs, indices) rebuilds a map from SuperTerm => coefficient pairs and a summation scope. Maps compose with ∘ or *, the right factor acting first, and S(ρ) applies a map. transpose(S) satisfies Tr(A * S(ρ)) == Tr(transpose(S)(A) * ρ) and is the Heisenberg-picture map, while adjoint(S) is the Hilbert–Schmidt adjoint. Iterating a map gives its SuperTerm => coefficient pairs.

A map is not an operator, so it does not enter operator products or averages. The name is public but not exported, since the numeric backends export a SuperOperator of their own.

Examples

julia> using SecondQuantizedAlgebra: SuperOperatorjulia> h = FockSpace(:c); a = Destroy(h, :a);julia> SuperOperator(dissipator(a))-1//2 * (·) ⋅ a' * a + a ⋅ (·) ⋅ a' - 1//2 * a' * a ⋅ (·)

source

SecondQuantizedAlgebra.sprepost — Function

sprepost(A, B)

The map that multiplies an operator by A from the left and by B from the right.

It acts as ρ ↦ A * ρ * B. A and B are operators or scalars. A side with a formal operator function gives a SuperExpr. A side that carries its own Σ throws; sum the map instead, Σ(sprepost(Aᵢ, Bᵢ), i).

Examples

julia> using SecondQuantizedAlgebra: sprepostjulia> h = FockSpace(:c); a = Destroy(h, :a);julia> S = sprepost(a, a')a ⋅ (·) ⋅ a'julia> S(a' * a)1 + 3 * a' * a + a' * a' * a * a

See also spre, spost, SuperOperator.

source

SecondQuantizedAlgebra.spre — Function

spre(A)

The map that multiplies an operator by A from the left, ρ ↦ A * ρ.

Examples

julia> using SecondQuantizedAlgebra: sprejulia> h = FockSpace(:c); a = Destroy(h, :a);julia> spre(a)(a')1 + a' * a

GKSL form of a map

A map written in sandwiches, or one that comes out of a composition or a frame change, is read back as a Lindbladian. The extraction is exact in the operator algebra and needs no Hilbert-space cutoff.

Lindbladian(κ * (sprepost(a, a') - (1 // 2) * spre(a' * a) - (1 // 2) * spost(a' * a)))

\[\kappa\,\mathcal{D}\left[a\right]\]

Lindbladian(SuperOperator(Lrot)) == Lrot

true

The Kossakowski matrix depends on the basis of jumps. In the basis $(a, a + a^\dagger a)$, independent loss and dephasing become a correlated matrix.

kossakowski(κ * dissipator(a) + γ * dissipator(a' * a), JumpBasis(a, a + a' * a))

\[ \left[ \begin{array}{cc} \gamma + \kappa & - \gamma \\ - \gamma & \gamma \\ \end{array} \right] \]

Basis operators are taken modulo the identity and with the completeness relation expanded.

JumpBasis(σ(:g, :g), a + 1)

JumpBasis with 2 operators:
  1: -σ₂₂
  2: a

SecondQuantizedAlgebra.Lindbladian — Method

Lindbladian(S::SuperOperator[, basis::JumpBasis]; rtol = 0)

A map in GKSL form, read back as a structured generator. The result has the Hamiltonian of S and one channel d_{μν} 𝒟[F_μ, F_ν] for every nonzero entry of the Kossakowski matrix in basis, so SuperOperator(Lindbladian(S, basis)) equals S up to completeness expansion (and up to the terms a positive rtol drops). Without a basis the support_basis of S is used. The Hamiltonian is defined modulo the identity, which hamiltonian_action drops.

Throws and handles floating-point coefficients like kossakowski. The rates are not checked for positivity.

Examples

julia> @qnumbers a::Destroy(FockSpace(:c));julia> using SecondQuantizedAlgebra: SuperOperatorjulia> Lindbladian(SuperOperator(dissipator(a + 3)))-i[3//2im * a + -3//2im * a', ·] + 𝒟[a]

source

SecondQuantizedAlgebra.JumpBasis — Type

JumpBasis(operators...)
JumpBasis(operators::AbstractVector)

An ordered basis of jump operators $(F_1, …, F_q)$ in which a map is written in GKSL form,

L(ρ) = -i[H, ρ] + Σ_{μν} d_{μν} (F_μ ρ F_ν' - {F_ν' F_μ, ρ}/2),

with d the Kossakowski matrix returned by kossakowski. Each operator is canonicalized with expand_completeness and taken modulo the identity, since a scalar shift of a jump changes only the Hamiltonian. The operators must be linearly independent modulo the identity, and an ArgumentError is thrown otherwise, including for an operator proportional to the identity; they need not be orthogonal, normalized or complete. Their order fixes the order of d, so bases with the same span in a different order are different. An empty basis is the basis of a coherent map.

A JumpBasis iterates and indexes like a vector of its canonical operators.

Examples

julia> @qnumbers a::Destroy(FockSpace(:c));julia> JumpBasis(a + 2, a' * a)JumpBasis with 2 operators:  1: a  2: a' * a

source

SecondQuantizedAlgebra.kossakowski — Function

kossakowski(S::SuperOperator, basis::JumpBasis; rtol = 0) -> Matrix{Coeff}
kossakowski(L::Lindbladian, basis::JumpBasis; rtol = 0) -> Matrix{Coeff}

The Hermitian Kossakowski matrix d of a map in the ordered basis:

S(ρ) = -i[H, ρ] + Σ_{μν} d_{μν} (F_μ ρ F_ν' - {F_ν' F_μ, ρ}/2).

The extraction works in the operator algebra and introduces no Hilbert-space cutoff. It throws an ArgumentError when the map has a jump direction outside basis, and a GKSLFormError when the map is not of GKSL form. A Lindbladian is read through its SuperOperator, so its own channels need not lie in basis.

The extraction is exact. Floating-point coefficients carry the roundoff of building the map, so a float map is usually of GKSL form only up to that roundoff and is rejected; write rates as exact rationals such as 1//10 or as symbolic parameters instead. A positive rtol opts into a heuristic: a leftover float coefficient counts as roundoff when, for every parameter monomial, it is at most rtol times the size of the coefficients it was computed from. No single rtol suits every basis, since roundoff grows as the basis operators approach linear dependence, and genuine terms below rtol are dropped. Exact coefficients are always compared exactly.

Examples

julia> @qnumbers a::Destroy(FockSpace(:c));julia> using Symbolics: @variablesjulia> @variables κ::Real γ::Real;julia> kossakowski(κ * dissipator(a) + γ * dissipator(a' * a), JumpBasis(a, a + a' * a))2×2 Matrix{SecondQuantizedAlgebra.Coeff}: γ + κ  -γ -γ     γ

source

SecondQuantizedAlgebra.support_basis — Function

support_basis(S::SuperOperator) -> JumpBasis
support_basis(L::Lindbladian) -> JumpBasis

The JumpBasis of the operator monomials on the two-sided sandwiches of a map, after expand_completeness, in term_order_key order. Every jump direction of the map lies in its span, so Lindbladian(S) reads S in this basis.

Examples

julia> using SecondQuantizedAlgebra: support_basisjulia> @qnumbers b::Destroy(FockSpace(:c));julia> support_basis(dissipator(b) + dissipator(b' * b))JumpBasis with 2 operators:  1: b  2: b' * b

source

SecondQuantizedAlgebra.independent_basis — Function

independent_basis(candidates::AbstractVector) -> JumpBasis

The JumpBasis of those candidates that are linearly independent modulo the identity of the candidates before them, kept in order. A candidate proportional to the identity is skipped.

Examples

julia> using SecondQuantizedAlgebra: independent_basisjulia> @qnumbers a::Destroy(FockSpace(:c));julia> independent_basis([a, 2a + 1, one(a), a' * a])JumpBasis with 2 operators:  1: a  2: a' * a

Correlated and collective decay

A cross channel $\mathcal{D}[J, K]$ with $J \neq K$ describes two decay paths into a shared bath, such as the cavity and the atom radiating into the same outgoing mode. liouvillian builds a generator from a list of jumps and a rate matrix in the QuantumOptics convention, where $\Gamma_{ij}$ multiplies $J_i \rho J_j^\dagger$. The off-diagonal entries become cross channels.

using SecondQuantizedAlgebra: liouvillian@variables Γcliouvillian(0, [a, σ(:g, :e)]; rates = [κ Γc; Γc γ])

\[\kappa\,\mathcal{D}\left[a\right] + \mathtt{{\Gamma}c}\,\mathcal{D}\left[a,{\sigma}^{{12}}\right] + \mathtt{{\Gamma}c}\,\mathcal{D}\left[{\sigma}^{{12}},a\right] + \gamma\,\mathcal{D}\left[{\sigma}^{{12}}\right]\]

SecondQuantizedAlgebra.liouvillian — Function

liouvillian(H, J; rates = nothing)

The Lindblad generator of the Hamiltonian H and the jump operators J, in the QuantumOptics convention.

The generator is -i[H, ·] + Σᵢⱼ Γ[i, j] 𝒟[J[i], J[j]]. rates is nothing for unit rates, a vector of rates, or a matrix Γ whose entry Γ[i, j] multiplies J[i] ρ J[j]'. An empty J gives the coherent part only. The name is public but not exported, since QuantumOpticsBase exports a numerical liouvillian.

Examples

julia> using SecondQuantizedAlgebra: liouvillianjulia> h = FockSpace(:c); a = Destroy(h, :a);julia> liouvillian(a' * a, [a, a']; rates = [1, 2])-i[a' * a, ·] + 𝒟[a] + 2 * 𝒟[a']

source

For an ensemble of emitters, a generator is summed over sites with Σ. Independent emission at rate $\gamma$ on every atom drains the total excitation at rate $\gamma$.

using SecondQuantizedAlgebraha = NLevelSpace(:atom, (:g, :e))@variables N γi = Index(ha, :i, N, ha)j = Index(ha, :j, N, ha)k = Index(ha, :k, N, ha)s(x, y, n) = IndexedOperator(Transition(ha, :s, x, y), n)Lind = Σ(γ * dissipator(s(:g, :e, i)), i)transpose(Lind)(Σ(s(:e, :e, j), j))

\[\underset{i}{\overset{N}{\sum}} - \gamma {s}_{i}^{{22}}\]

Emitters closer than a wavelength radiate collectively, with a rate matrix $\Gamma_{ij}$ over pairs of sites. dissipator(J, Γ, i => j) builds $\sum_{i,j} \Gamma_{ij} \mathcal{D}[J_i, J_j]$ and splits off the diagonal $i = j$, so the generator shows its local and correlated parts.

Lcol = dissipator(s(:g, :e, i), DoubleIndexedVariable(:Γ, i, j), i => j)

\[\underset{i}{\overset{N}{\sum}} \Gamma\left( i, i \right)\,\mathcal{D}\left[{s}_{i}^{{12}}\right] + \underset{i}{\overset{N}{\sum}} \underset{j,i{\neq}j}{\overset{N}{\sum}} \Gamma\left( i, j \right)\,\mathcal{D}\left[{s}_{i}^{{12}},{s}_{j}^{{12}}\right]\]

The population of one emitter $k$ then couples to the coherences between $k$ and every other emitter.

average(transpose(Lcol)(s(:e, :e, k)))

\[ \begin{equation} \underset{i{\neq}k}{\overset{N}{\sum}} - \frac{1}{2} ~ \langle {s}_{i}^{{12}}{s}_{k}^{{21}} \rangle ~ \Gamma\left( i, k \right) + \underset{j{\neq}k}{\overset{N}{\sum}} - \frac{1}{2} ~ \langle {s}_{j}^{{21}}{s}_{k}^{{12}} \rangle ~ \Gamma\left( k, j \right) - \langle {s}_{k}^{{22}} \rangle ~ \Gamma\left( k, k \right) \end{equation} \]

A map shares one summation scope between its two sides, as in $\sum_i J_i \rho J_i^\dagger$. Constructions that would make the scope ambiguous throw instead of guessing, and the algebra never invents an index to resolve them:

composing two maps that bind the same index; rename one with change_index;
applying a map to an operator in which its bound index is free, since the sum would capture it, which is why the drift above is taken on s(:e, :e, k);
a jump that carries its own Σ, since $J \rho J^\dagger$ with $J = \sum_i J_i$ needs two independent copies of the index;
powers of a map with bound indices, for the same reason.

julia> Lind ∘ LindERROR: ArgumentError: summation index i is bound by both operands. Bound indices of composed maps, or of a map and the operator it acts on, must be distinct; rename one side with `change_index`.julia> dissipator(Σ(s(:g, :e, i), i))ERROR: ArgumentError: a dissipator cannot take an operator with its own bound sum: `J ρ J'` with `J = Σᵢ Jᵢ` needs two independent copies of the summation index. Use `dissipator(Jᵢ, Γ, i => j)` for collective decay or `Σ(dissipator(Jᵢ), i)` for independent decay, and `Σ(sprepost(Aᵢ, Bᵢ), i)` for a summed sandwich.

Formal operator functions

A Hamiltonian or a jump with a formal function such as cos, sin or expim gives a Lindbladian{QExpr}, as for a transmon with its cosine potential and dephasing of the phase. taylor lowers it to a polynomial generator.

using SecondQuantizedAlgebraφ = Position(PhaseSpace(:transmon), :φ)@variables γφLt = hamiltonian_action(-cos(φ)) + γφ * dissipator(φ)

\[-i\left[-\cos\left( \hat{\varphi} \right),\cdot\right] + \mathtt{\gamma\varphi}\,\mathcal{D}\left[\hat{\varphi}\right]\]

taylor(Lt, 0:2)

\[-i\left[\frac{1}{2} \hat{\varphi}\hat{\varphi},\cdot\right] + \mathtt{\gamma\varphi}\,\mathcal{D}\left[\hat{\varphi}\right]\]

A general map with a formal side is a SuperExpr, and taylor lowers it in the same way.

using SecondQuantizedAlgebra: spreposttaylor(sprepost(cos(φ), φ), 0:2)

\[\cdot\,\hat{\varphi} - \frac{1}{2} \hat{\varphi} \hat{\varphi}\,\cdot\,\hat{\varphi}\]

SecondQuantizedAlgebra.SuperExpr — Type

SuperExpr
SuperExpr(L::Lindbladian)

A linear map on operators whose sides may contain formal operator functions.

A SuperExpr is a polynomial SuperOperator plus sandwiches c * f ρ g with a QExpr side. It arises from sprepost with a formal side, from operations on a Lindbladian{QExpr}, and from SuperExpr(L). It supports the arithmetic, composition, transpose and application of a SuperOperator, and keeps its type even when no formal sandwich remains. taylor(X, 0:n) lowers every formal side and returns a SuperOperator, and SuperOperator(X) converts an X without formal sandwiches.

Examples

julia> using SecondQuantizedAlgebra: sprepostjulia> h = FockSpace(:c); a = Destroy(h, :a);julia> X = sprepost(cos(a' * a), a)cos(a' * a) ⋅ (·) ⋅ ajulia> X isa SecondQuantizedAlgebra.SuperExprtrue

See also SuperOperator, taylor, poly.

source

Inspecting and building maps

The parts of a generator are read back with accessors. hamiltonian and channels split a Lindbladian into its Hamiltonian and its Channel => rate pairs, and jump, partner and scope read a channel.

using SecondQuantizedAlgebra: hamiltonian, channels, jump, partner, scopech, rate = first(channels(Lrot))jump(ch), partner(ch), rate

(a, a, κ)

Iterating a SuperOperator gives its SuperTerm => coefficient pairs, whose sides are read with left_operator and right_operator. A SuperExpr splits into its polynomial part and its formal sandwiches with poly and formal_terms.

using SecondQuantizedAlgebra: left_operator, right_operatorterm, c = first(D)left_operator(term), right_operator(term), c

(a' * a, 1, (-1//2))

sum over an array of Lindbladians or SuperOperators accumulates in place, like sum over operators, so a generator with many channels costs linear time. A loop uses a builder directly.

using SecondQuantizedAlgebra: LindbladianBuilder, QAdd, accumulate!, buildhm = reduce(⊗, [FockSpace(Symbol(:mode, n)) for n in 1:3])modes = [Destroy(hm, Symbol(:a, n), n) for n in 1:3]builder = LindbladianBuilder{QAdd}()for m in modes    accumulate!(builder, κ * dissipator(m))endbuild(builder) == sum([κ * dissipator(m) for m in modes])

true

SecondQuantizedAlgebra.hamiltonian — Function

hamiltonian(L::Lindbladian)
channels(L::Lindbladian)

The Hamiltonian of a Lindbladian, and a vector of its Channel => rate pairs in channel_order_key order.

Examples

julia> using SecondQuantizedAlgebra: hamiltonian, channels, jumpjulia> h = FockSpace(:c); a = Destroy(h, :a);julia> L = hamiltonian_action(a' * a) + 2 * dissipator(a);julia> hamiltonian(L)a' * ajulia> ch, rate = first(channels(L)); jump(ch), rate(a, 2)

Secular projection

When phases remain after the transformation, secular keeps the time-independent part of the Hamiltonian, the rates and the jumps, and returns a Lindbladian. For a quadrature jump this is the rotating-wave dissipator:

using SecondQuantizedAlgebra@qnumbers b::Destroy(FockSpace(:secular_mode))@variables ω::Real κ::Real ωd::Real t::Realimport SecondQuantizedAlgebra: secularLx = hamiltonian_action(ω * b' * b) + κ * dissipator(b + b')secular(transform(Lx, Rotation(b, ωd * t, t)), ωd, t)

\[-i\left[\left(\omega - \mathtt{{\omega}d}\right) b^{\dagger}b,\cdot\right] + \kappa\,\mathcal{D}\left[b\right] + \kappa\,\mathcal{D}\left[b^{\dagger}\right]\]