Many-body systems

An ensemble of identical atoms or an array of coupled modes is written once per site and summed. Consider $N$ two-level atoms in a cavity, the Tavis–Cummings model

\[H = \Delta a^\dagger a + \sum_{i=1}^{N} g_i \left( a^\dagger \sigma_i^{12} + a \sigma_i^{21} \right).\]

The site label $i$ and the number of atoms $N$ both stay symbolic, so equations are derived once for any $N$ and any distribution of couplings. When the atoms are indistinguishable, collective operators describe the same ensemble without any site label.

Indices
Indexed operators and parameters
Sums
Explicit sites
Collective systems

Indices

A site label is an Index. It belongs to one subsystem of the product space and runs over a range that may be a number or a symbol. For the Tavis–Cummings model the index runs over the atoms, the second factor.

using SecondQuantizedAlgebra@variables N Δh = FockSpace(:cavity) ⊗ NLevelSpace(:atoms, 2)i = Index(h, :i, N, 2)j = Index(h, :j, N, 2)

SecondQuantizedAlgebra.Index — Type

Index(h::HilbertSpace, name::Symbol, range, space::HilbertSpace)
Index(h::HilbertSpace, name::Symbol, range, subsystem::Int)
i(k::Integer)

Symbolic label for one of many identical sites.

An index runs from 1 to range, an integer or a symbolic parameter such as N, over the factor of h given either by the space itself or by its position. Attached to operators with IndexedOperator and to parameters with IndexedVariable, it is summed over with Σ. Two indices are equal only if they have the same name, range and subsystem.

Calling an index with an integer, i(3), gives the index of a concrete site. It is named i_3, keeps the range and subsystem of i, and addresses that site after a sum has been written out site by site, for instance in numerical conversion or in cumulant equations.

Examples

julia> @variables N;julia> h = FockSpace(:cavity) ⊗ NLevelSpace(:atom, 2);julia> i = Index(h, :i, N, NLevelSpace(:atom, 2));julia> i == Index(h, :i, N, 2)truejulia> σ = Transition(h, :σ, 2, 1);julia> IndexedOperator(σ, i), IndexedOperator(σ, i(3))(σ_i₂₁, σ_i_3₂₁)

source

Replaces `ClusterSpace`

Earlier versions inherited ClusterSpace from QuantumCumulants.jl, which wrapped a fixed number of identical copies of a space. It has been removed, together with cluster_expand and insert_index, in favor of an Index with a symbolic range. The number of copies then stays symbolic, and the truncation order of a cumulant expansion is set in the solver rather than in the Hilbert space.

Indexed operators and parameters

An operator is placed on a site with IndexedOperator. Wrapping it in a function of the levels and the index keeps the notation close to $\sigma_i^{kl}$.

@qnumbers a::Destroy(h, 1)σ(k, l, m) = IndexedOperator(Transition(h, :σ, k, l), m)σ(2, 1, i)

\[{\sigma}_{i}^{{21}}\]

Operators on different sites commute, while operators on the same site obey the local operator algebra. To retain individual atom labels, use indexed Transition operators as above and Σ to form their collective sums.

A CollectiveTransition instead represents an entire ensemble; see Collective systems. An index on such an operator could label separate ensembles, each with its own collective algebra, but IndexedOperator does not currently support this. For a fixed number of ensembles, give each its own CollectiveNLevelSpace in a product space.

Site-dependent parameters such as the couplings $g_i$ are made with IndexedVariable, and pair parameters with DoubleIndexedVariable.

g = IndexedVariable(:g, i)

\[ \begin{equation} g\left( i \right) \end{equation} \]

SecondQuantizedAlgebra.IndexedOperator — Function

IndexedOperator(op, i::Index)

Copy of the operator op placed on site i.

Operators on different sites act independently; operators on the same site follow the algebra of op. Any operator except a CollectiveTransition can be indexed.

Examples

julia> h = NLevelSpace(:atom, 2); i = Index(h, :i, 10, h);julia> σ = Transition(h, :σ, 2, 1);julia> IndexedOperator(σ, i) * IndexedOperator(σ', i)σ_i₂₂

See also Index, Σ, operator_index.

source

SecondQuantizedAlgebra.IndexedVariable — Function

IndexedVariable(name::Symbol, i::Index)

Site-dependent parameter name(i), such as a local detuning or coupling.

The parameter is real. Under a sum it follows the index, so it is evaluated site by site once the sum is unrolled.

Examples

julia> h = NLevelSpace(:atom, 2); i = Index(h, :i, 10, h);julia> g = IndexedVariable(:g, i)g(i)julia> Σ(g * IndexedOperator(Transition(h, :σ, 2, 1), i), i)Σ(i=1:10) g(i) * σ_i₂₁

source

SecondQuantizedAlgebra.DoubleIndexedVariable — Function

DoubleIndexedVariable(name::Symbol, i::Index, j::Index; identical = true)

Pair parameter name(i, j), such as a hopping amplitude or a dipole-dipole coupling.

With identical = false the parameter vanishes whenever both indices refer to the same site, so the diagonal terms of a double sum drop out.

Examples

julia> h = NLevelSpace(:atom, 2); i = Index(h, :i, 10, h); j = Index(h, :j, 10, h);julia> σ(k, l, m) = IndexedOperator(Transition(h, :σ, k, l), m);julia> J = DoubleIndexedVariable(:J, i, j; identical = false);julia> Σ(J * σ(2, 1, i) * σ(1, 2, j), i, j)Σ(i=1:10)Σ(j=1:10)(i≠j) J(i, j) * σ_i₂₁ * σ_j₁₂

Sums

With these pieces the Hamiltonian is a single Σ.

H = Δ * a' * a + Σ(g * (a' * σ(1, 2, i) + a * σ(2, 1, i)), i)

\[\Delta a^{\dagger}a + \underset{i}{\overset{N}{\sum}}\left( g\left( i \right) a{\sigma}_{i}^{{21}} + g\left( i \right) a^{\dagger}{\sigma}_{i}^{{12}} \right)\]

A sum takes part in the algebra like any other expression. The Heisenberg equation of the cavity field collects the contribution of every atom.

1im * commutator(H, a)

\[ - \Delta ~ \mathit{i} a + \underset{i}{\overset{N}{\sum}} - g\left( i \right) ~ \mathit{i} {\sigma}_{i}^{{12}}\]

SecondQuantizedAlgebra.Σ — Function

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

Symbolic sum of expr over all sites of the index i.

Several indices give a nested sum. The sum stays symbolic and follows the algebra like any other expression. Where its index meets another index on the same subsystem, the terms in which the two sites coincide are split off, so the remaining sum carries the constraint i ≠ j. The optional non_equal adds such constraints by hand, excluding those sites from the sum. A summand that does not depend on i is multiplied by the number of sites.

∑ (typed \sum<tab>) is the same function.

Examples

julia> @variables N;julia> h = NLevelSpace(:atom, 2); i = Index(h, :i, N, h); j = Index(h, :j, N, h);julia> σ(k, l, m) = IndexedOperator(Transition(h, :σ, k, l), m);julia> Σ(σ(2, 1, i) * σ(1, 2, j), i, j)Σ(j=1:N) σ_j₂₂ + Σ(i=1:N)Σ(j=1:N)(i≠j) σ_i₂₁ * σ_j₁₂julia> Σ(σ(2, 2, j), i, [j])(-1 + N) * σ_j₂₂

source

Coinciding sites

An index under a sum also runs over the site of every other index on the same subsystem. On that site the two operators act on one atom and compose, so the term is split off and the remaining sum excludes it.

Σ(σ(2, 2, i), i) * σ(2, 1, j)

\[{\sigma}_{j}^{{21}} + \underset{i{\neq}j}{\overset{N}{\sum}}{\sigma}_{i}^{{22}}{\sigma}_{j}^{{21}}\]

A ground-state projector produced on the shared site stays as it is, as everywhere else.

Σ(σ(1, 2, i), i) * σ(2, 1, j)

\[{\sigma}_{j}^{{11}} + \underset{i{\neq}j}{\overset{N}{\sum}}{\sigma}_{i}^{{12}}{\sigma}_{j}^{{21}}\]

Free indices

Two indices that are not summed over may or may not label the same site. Their operators are therefore left in the order they were written, and no rule is applied between them.

σ(2, 2, j) * σ(1, 2, i)

\[{\sigma}_{j}^{{22}}{\sigma}_{i}^{{12}}\]

When the indices are known to label different sites, assume_distinct_index says so, and the product is brought into canonical form.

assume_distinct_index(σ(2, 2, j) * σ(1, 2, i), [(i, j)])

\[{\sigma}_{i}^{{12}}{\sigma}_{j}^{{22}}\]

SecondQuantizedAlgebra.assume_distinct_index — Function

assume_distinct_index(expr, pairs::Vector{Tuple{Index, Index}})

Treat each pair of free indices in pairs as different sites.

Two free indices on the same subsystem may or may not refer to the same site, so their operators are left as they are. Outside a Σ, which settles this by itself, this function supplies the missing information: the operators are then sorted into canonical order, and ground-state projectors in the result are eliminated with expand_completeness.

Examples

julia> h = NLevelSpace(:atom, 2); i = Index(h, :i, 10, h); j = Index(h, :j, 10, h);julia> σ(k, l, m) = IndexedOperator(Transition(h, :σ, k, l), m);julia> σ(2, 1, j) * σ(1, 2, i)σ_j₂₁ * σ_i₁₂julia> assume_distinct_index(σ(2, 1, j) * σ(1, 2, i), [(i, j)])σ_i₁₂ * σ_j₂₁

Renaming indices

An equation derived for one site holds for every other site under a new index name, and a correlation between two sites is obtained from a single-site expression by renaming. change_index renames an index throughout and brings the result back into canonical form. The coherence of atom $i$ becomes that of atom $j$.

change_index(a' * σ(1, 2, i), i, j)

\[a^{\dagger}{\sigma}_{j}^{{12}}\]

SecondQuantizedAlgebra.change_index — Function

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

Replace the index from by to in operators and coefficients of expr.

The result is brought back into canonical form, so operators that now share a site compose, and a DoubleIndexedVariable declared with identical = false vanishes when its two indices become equal. The second form renames several indices at once, which makes swaps such as Dict(i => j, j => i) possible.

Examples

julia> h = NLevelSpace(:atom, 2); i = Index(h, :i, 10, h); j = Index(h, :j, 10, h);julia> σ(k, l, m) = IndexedOperator(Transition(h, :σ, k, l), m);julia> change_index(σ(2, 1, i) * σ(1, 2, j), i, j)σ_j₂₂julia> change_index(σ(2, 1, i) * σ(1, 2, j), Dict(i => j, j => i))σ_j₂₁ * σ_i₁₂

Inspecting indices

The indices of an expression and their properties are read off directly.

get_indices(a' * σ(1, 2, i) * σ(2, 1, j))

2-element Vector{Index}:
 i
 j

index_name(i), index_range(i)

(:i, N)

SecondQuantizedAlgebra.get_indices — Function

get_indices(expr) -> Vector{Index}

Indices occurring in expr, both summed and free.

Examples

julia> h = NLevelSpace(:atom, 2); i = Index(h, :i, 10, h); j = Index(h, :j, 10, h);julia> σ(k, l, m) = IndexedOperator(Transition(h, :σ, k, l), m);julia> get_indices(Σ(σ(2, 1, i) * σ(1, 2, j), i))2-element Vector{Index}: i j

Explicit sites

An index treats all of its sites alike. A lattice with nearest-neighbour couplings, open boundaries or a definite number of sites is written instead with one subsystem per site, as in the driven Bose–Hubbard chain of $M$ coupled cavities

\[H = \sum_{k=1}^{M} \left[ -\Delta\, a_k^\dagger a_k + \frac{U}{2}\, a_k^\dagger a_k^\dagger a_k a_k + \Omega\,(a_k + a_k^\dagger) \right] - J \sum_{k=1}^{M-1} \left( a_k^\dagger a_{k+1} + a_{k+1}^\dagger a_k \right).\]

Each site gets its own FockSpace, and both sums are a sum over a vector of terms.

using SecondQuantizedAlgebra@variables Δ U Ω JM = 6hc = reduce(⊗, [FockSpace(Symbol(:c, k)) for k in 1:M])a = [Destroy(hc, Symbol(:a, k), k) for k in 1:M]onsite = sum([-Δ * a[k]' * a[k] + (U / 2) * a[k]' * a[k]' * a[k] * a[k] + Ω * (a[k] + a[k]') for k in 1:M])hopping = sum([-J * (a[k]' * a[k+1] + a[k+1]' * a[k]) for k in 1:(M-1)])H = onsite + hoppinglength(H)

Summing a vector builds the result in place. A chain of + gives the same expression, but each + copies the sum accumulated so far, so $n$ terms cost $O(n^2)$ instead of $O(n)$.

terms = [Ω * (a[k] + a[k]') for k in 1:M]sum(terms) == foldl(+, terms)

true

For a few sites the difference is negligible. For long chains the copies outweigh the work of forming the terms themselves:

Sites $M$	Terms	Forming the terms	`foldl(+, terms)`	`sum(terms)`
24	142	0.4 ms	0.18 ms	0.008 ms
400	2398	6.3 ms	63 ms	0.22 ms

Use brackets, not a bare generator or a `+=` loop

The fast path is taken by sum, reduce(+, …) and MutableArithmetics.operate(sum, …) over an array. A bare generator and an accumulation loop both fall back to repeated +:

H = sum([Ω * (a[k] + a[k]') for k in 1:M])            # fastH = sum(Ω * (a[k] + a[k]') for k in 1:M)              # slow: generatorH = 0; for k in 1:M; H += Ω * (a[k] + a[k]'); end     # slow: copies on every +=

Wrapping the generator in brackets, or collecting the terms into a vector first, is enough. The extra vector of length $M$ is negligible next to the quadratic saving.

MutableArithmetics.@rewrite gives the same result, but its macro expansion adds overhead and it does not take the fast path for generators. It suits code that already composes a larger MutableArithmetics expression.

import MutableArithmetics as MAfast = sum([Ω * (a[k] + a[k]') for k in 1:M])fast == MA.operate(sum, [Ω * (a[k] + a[k]') for k in 1:M]) == MA.@rewrite(sum(Ω * (a[k] + a[k]') for k in 1:M))

The binary +, - and scalar * copy because the caller still holds both operands, so neither may be modified in place. Only a reduction over a whole collection can build its result in one buffer. The same holds wherever many terms are added up:

equations of motion and cumulant systems, where each derived $\frac{d}{dt}\langle O \rangle$ is a sum of many contributions;
ensembles written with explicit sites, for instance collective operators $S_\pm = \sum_i \sigma_i^\pm$;
Lindbladians with one dissipator per channel;
truncated Magnus, Dyson or BCH expansions.

Products are not affected; their cost lies in the reordering itself. A symbolic Σ stays a single sum however many sites it covers, so the question does not arise there.

Collective systems

Collective operators $S^{ij} = \sum_k \sigma_k^{ij}$ act on the ensemble as a whole and form a closed algebra without individual site labels. This algebra holds on the full atomic Hilbert space as well as within its permutation-symmetric sector. The numerical example below chooses a basis for that sector; this restriction belongs to the numerical representation, not to the symbolic commutation relations. The two descriptions suit different problems:

	Indexed sites	Collective operators
Operators	$\sum_k \sigma_k^{ij}$ with `Index` and `Σ`	`CollectiveTransition` $S^{ij}$
Number of atoms	symbolic	fixed by the numerical basis
Site-dependent parameters	yes	no
Site-resolved moments and correlations	yes	no
Numerical basis	full product space for a concrete $N$	permutation-symmetric sector in the example below

They encode different physical roles and live on different Hilbert-space types, so they do not mix on one subsystem.

SecondQuantizedAlgebra.CollectiveNLevelSpace — Type

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

Hilbert space of an ensemble of identical n-level systems, home of the CollectiveTransition operators.

The number of systems is not part of the space; it is fixed only by the basis chosen for numerics. Levels are numbered or named as for NLevelSpace. There is no ground state to eliminate, since the collective projectors add up to the number of systems rather than to one.

Examples

julia> h = CollectiveNLevelSpace(:atoms, (:g, :e))ℋ(atoms)julia> S = CollectiveTransition(h, :S, :e, :g);julia> commutator(S, S')-S₁₁ + S₂₂

source

SecondQuantizedAlgebra.CollectiveTransition — Function

CollectiveTransition(h::CollectiveNLevelSpace, name::Symbol, i, j)
CollectiveTransition(h::ProductSpace, name::Symbol, i, j, [subsystem::Int])

Collective transition from level j to level i, summed over all systems of the ensemble.

The operator is $S^{ij} = \sum_k |i\rangle_k\langle j|$. Collective transitions obey $[S^{ij}, S^{kl}] = \delta_{jk}S^{il} - \delta_{li}S^{kj}$ and are normal ordered like ladder operators. Unlike a single Transition, a product of two of them does not compose, since it includes pairs of different systems. On a ProductSpace, subsystem selects the CollectiveNLevelSpace factor and may be omitted when there is only one.

Examples

julia> h = CollectiveNLevelSpace(:atoms, (:g, :e));julia> @qnumbers S::CollectiveTransition(h, :e, :g);julia> S' * SS₁₁ - S₂₂ + S₂₁ * S₁₂

source

The number of systems enters only through the numerical basis. A bosonic many-body basis over an NLevelBasis with $N$ particles is the permutation-symmetric sector, of dimension $\binom{N+n-1}{n-1}$ for $n$ levels, so $N + 1$ for two.

using SecondQuantizedAlgebra, QuantumOpticsBaseh = CollectiveNLevelSpace(:atoms, (:g, :e))S(k, l) = CollectiveTransition(h, :S, k, l)N = 10b1 = NLevelBasis(2)b = ManyBodyBasis(b1, bosonstates(b1, N))length(b)

In this basis the collective projectors add up to the number of atoms.

to_numeric(S(:g, :g) + S(:e, :e), b) ≈ N * identityoperator(b)

true

For two levels this is the collective spin of length $N/2$. With $e$ as the upper level, $S^{eg} = S_+$, $S^{ge} = S_-$ and $S_z = (S^{ee} - S^{gg})/2$. In the ordering of the SpinBasis of QuantumOpticsBase, where level 1 is spin up, the same identification reads $S^{12} = S_+$.

get_sum_indices(expr) reads only the bound summation indices. It returns an empty vector for an unsummed expression.

SecondQuantizedAlgebra.get_sum_indices — Method

get_sum_indices(q::QAdd) -> Vector{Index}
get_sum_indices(q::QSym) -> Vector{Index}

The summation indices bound by q, in the order of its sum scope. A QSym and an unsummed QAdd bind none. The result is a copy, so mutating it does not change q.

julia> h = FockSpace(:site);julia> @qnumbers a::Destroy(h);julia> i = Index(h, :i, 3, h);julia> SecondQuantizedAlgebra.get_sum_indices(Σ(IndexedOperator(a, i), i)) == [i]true

source