API reference

FermionicHilbertSpaces

Symbolic fermions, bosons and spins, Hilbert spaces with constraints, and the tensor products, partial traces and matrix representations that connect them.

hc

Adding this is equivalent to adding the hermitian conjugate.

AdditiveConstraint(allowed_values; spaces=missing, maps=missing)
AdditiveConstraint(allowed_values, maps)

Constraint enforcing that the sum of per-subspace contributions lies in allowed_values. The contribution of each subspace in spaces is the value of its map applied to the subspace's state. A single function in maps is applied to all spaces; a collection of maps supplies one function per subspace.

If spaces is missing, tensor_product uses the spaces it combines, and constrain_space treats the whole space as a single subspace.

Axes(H1, H2, ...)

The Hilbert spaces of an array's axes, in axis order: Axes(Hout, Hin) for an operator from Hin to Hout, Axes(H) for a vector. A bare space H means H on every axis. Public, not exported.

Examples

using FermionicHilbertSpacesusing FermionicHilbertSpaces: Axes@fermions fH2 = hilbert_space(f, 1:2, NumberConservation(2))H1 = hilbert_space(f, 1:2, NumberConservation(1))representation(f[1], Axes(H1, H2)) # a dim(H1) × dim(H2) matrix

BosonField

A family of bosonic modes, created by @bosons. b[i] is the annihilation operator of mode i.

BosonicState(n)

The basis state with n bosons in a single mode, as returned by basisstates of a bosonic space.

BranchConstraint(f)

Constraint for pruning branches while generating basis states. f(partial_state, depth, spaces) -> Bool is evaluated during backtracking by generate_states, after the state at position depth has been assigned. Entries after depth in partial_state have not yet been assigned.

Returning true continues exploring the branch, while returning false prunes it. The callback should therefore accept partial states and reject a branch only when no valid completion can satisfy the constraint.

BranchConstraint is intended for use with tensor_product or more directly with generate_states.

FilterConstraint(reducer; maps=missing, spaces=missing)
FilterConstraint(reducer, maps, spaces)

Constraint that keeps complete basis states for which reducer returns true.

With only reducer, the reducer is called directly as reducer(state). When maps and spaces are provided, the state is split into the selected spaces, each map is applied to its corresponding subspace, and the reducer is called on the resulting values. A single function is applied to all spaces; a collection of maps supplies one function per subspace.

FockNumber

A type representing a Fock state as the bitstring of an integer. Bit k - 1 (counting from the least significant bit) is the occupation of the k-th mode of the space, in the order of the labels given to hilbert_space. In an unconstrained space the basis state at index i is FockNumber(i - 1).

GenericHilbertSpace(label, basisstates[, state_index])
GenericHilbertSpace{B}(label, basisstates[, state_index])

Atomic Hilbert space whose basis states are given explicitly by basisstates, indexed by integers. label is returned by symbolic_group and atomic_id and identifies the space in symbolic expressions. state_index gives the integer index of a basis state. It is either an AbstractDict from basis states to indices, or a function that takes a basis state and returns its index. By default it is a dictionary built from basisstates.

B is the basis-state type and defaults to eltype(basisstates). When B is given and basisstates is an AbstractVector or Tuple, its elements are converted to B.

Kets(H)

Callable that turns a label into a symbolic ket of H, as in Kets(H)("10"). The label is parsed by interpret_state. Kets(H)[i] is the ket of basisstate(i, H).

julia> using FermionicHilbertSpaces; @fermions f;julia> k = Kets(hilbert_space(f, 1:2));julia> k[2] == k("10")true

NoSymmetry()

Constraint that leaves a Hilbert space unchanged.

NumberConservation(total=missing; spaces=missing, weights=missing)
NumberConservation(total, spaces; weights=missing)
NumberConservation(total, spaces, weights)
NumberConservation(space)

Constraint enforcing conservation of a (possibly weighted) particle number, optionally on selected spaces. total can be a single value or collection of allowed values; missing allows any value. spaces restricts the count to the given subspaces (default: all of them). weights gives one weight per subspace in spaces (default: all 1), and the conserved quantity is the weighted sum of their particle numbers.

ParityConservation(parities=[-1, 1]; spaces=missing)
ParityConservation(parities, spaces)
ParityConservation(space)

Constraint enforcing allowed fermion parities, optionally on selected spaces. parities is a parity or a collection of parities, each 1 (even) or -1 (odd). Sectors are ordered by increasing quantum number, so the odd sector comes first.

ProductState(states::Tuple)

A basis state of a tensor product space, holding one basis state per factor, e.g. ProductState((SpinState(1//2), BosonicState(0))).

SectorConstraint(reducer; maps=missing, spaces=missing)
SectorConstraint(reducer, maps, spaces)

Constraint that groups complete basis states into sectors according to a custom rule. It has the same reducer, maps, and spaces interface as FilterConstraint, but reducer must return the sector label for a state. Returning missing discards the state; states with the same non-missing label are placed in the same sector.

When maps and spaces are provided, the state is split into the selected spaces, the maps are applied to them, and the reducer receives the resulting values. A single function is applied to all spaces; a collection of maps supplies one function per subspace.

SectorHilbertSpace

Hilbert space whose basis is grouped into sectors labeled by quantum numbers. Use quantumnumbers, sector, and indices to access individual sectors.

SpinField

A family of spins, created by @spins. s[i] is the spin on site i, e.g. s[i][:z].

SpinState(m)

The basis state of a single spin with s[:z] eigenvalue m, e.g. SpinState(-1//2).

Ket(label, H)

The symbolic ket of H with the given label, such as the bitstring "10". See Kets for the multi-state form.

julia> @fermions f;julia> Ket("10", hilbert_space(f, 1:2))|10⟩

basisstate(i, H)

Return the i-th basis state of the Hilbert space H, in the order of basisstates(H). Throws a BoundsError if i is not in 1:dim(H). A string or character that spells a state of H is accepted in place of an index, and an ArgumentError is thrown if it is not a state of H.

basisstates(H)

Return the basis states of the Hilbert space H, in the order used by matrix representations. The result has length(basisstates(H)) == dim(H) and can be iterated and indexed, with basisstates(H)[i] == basisstate(i, H). Use collect to get a Vector.

constrain_space(space, constraint[, states])
constrain_space(space, states)

Build a constrained Hilbert space from space.

With a constraint, keep the basis states of space that satisfy it. If states is given, only those states are candidates. Constraints that group states into sectors, such as NumberConservation, return a space with sectors, with states ordered by sector. constrain_space(space, NoSymmetry()) returns space unchanged.

With only states, return the subspace spanned by exactly those basis states, in the given order; this is equivalent to constrain_space(space, NoSymmetry(), states).

dim(H)

Return the Hilbert-space dimension of H, i.e. the number of basis states.

dissipator(L)

The symbolic Lindblad dissipator ρ ↦ LρL† - ½{L†L, ρ} for the jump operator L, acting on the vectorized density matrix of a space returned by open_system. Public, not exported.

Examples

using FermionicHilbertSpacesusing FermionicHilbertSpaces: open_system, dissipator@fermions cHlr, Hl, Hr, left, right = open_system(c, 1:1)representation(dissipator(c[1]), Hlr)

embed(m, Hsub => H; complement=complementary_subsystem(H, Hsub), kwargs...)

Compute the embedding of a matrix m in the basis Hsub into the basis H. Fermionic phase factors are included if the two spaces are fermionic Hilbert spaces. Hsub => H reads input => output; see Pair arguments.

The result is the compression P (m ⊗ I) P onto H, so that tr(embed(A, Hsub => H) * ρ) == tr(A * partial_trace(ρ, H => Hsub)) for ρ on H. Unless projection=true, throws an ArgumentError if m ⊗ I maps a state of H outside H with nonzero amplitude. Columns of m for states of Hsub that no state of H contains are dropped, and states of H whose part on Hsub is not in Hsub get zero.

If m is a vector, Hsub and H must have the same modes, as for a sector Hsub of H. The result is m placed at the indices indices(Hsub, H) and zero elsewhere.

julia> using FermionicHilbertSpacesjulia> @fermions f;julia> H = hilbert_space(f, 1:2, NumberConservation());julia> embed([1, 2], sector(1, H) => H)4-element Vector{Int64}: 0 1 2 0

embed(Hsub => H; kwargs...)

Create a callable embedding operation from Hsub into H.

The operation can be applied directly as op(m) and converted to a sparse matrix map with sparse(op). Hsub => H reads input => output; see Pair arguments. It is the compression onto H and, unlike embed, does not check its input.

expectation(op, state, H)

Return the expectation value of the operator op in state on the space H. A vector ψ gives ψ' * M * ψ and a density matrix ρ gives tr(ρ * M), where M = representation(op, H). The state is assumed to be normalized.

On hilbert_space(f, labels, FreeFermions.SingleParticle()), state is the correlation matrix G[i, j] = ⟨f_i' f_j⟩, and on hilbert_space(f, labels, FreeFermions.BdG()) it is the Nambu correlation matrix G[i, j] = ⟨Ψ_i' Ψ_j⟩, with Ψ = (f; f') in the order of the BdG space. On the Majorana space hilbert_space(γ, labels, FreeFermions.SingleParticle()) it is G[i, j] = ⟨γ_i γ_j⟩ / 2. The state is then the Gaussian state with that correlation matrix, and op of any degree is evaluated with Wick's theorem. Terms that do not conserve the particle number (single-particle space) or the parity (BdG and Majorana spaces) throw an error, unless projection=true, which drops them.

Examples

@fermions fH = hilbert_space(f, 1:2)ψ = representation("10", H)expectation(f[1]' * f[1] + 2, ψ, H) == 3

factors(H)

Split H one level down: a product space into its groups, a grouped space (such as a fermionic space) into its atoms, and an atomic space into itself. Use this when building partitions, permutations and constraints. Internal code that needs the state layout should use groups instead.

generalized_kron(ms, Hs, H::AbstractHilbertSpace=tensor_product(Hs))

Compute the fermionic tensor product of the matrices or vectors in ms on the spaces Hs, as an array on H, which defaults to the tensor_product of Hs.

For vectors, the result is the product state, the same as tensor_product(ms, Hs => H). For matrices, write each local basis operator as a product of single-mode factors (one of c c†, c†, c, c†c per mode) in the mode order of its space; generalized_kron multiplies the factors of all spaces together in the mode order of H, and extends this linearly. The result does not depend on the order of ms and Hs, and is in general not a product of embeddings: for A = c₁†c₃ on modes [1, 3] and B = c₂†c₄ on modes [2, 4] it gives c₁†c₂†c₃c₄ = -(c₁†c₃)(c₂†c₄). It equals tensor_product(ms, Hs => H) when each space in Hs covers consecutive modes of H and Hs is listed in that order. Otherwise generalized_kron of density matrices need not be a density matrix. Use tensor_product for product states, local observables, and anything that should equal a product of embedded operators.

generalized_kron(ms, Hs => H) is also accepted, and generalized_kron(Hs => H) returns a function of ms. Hs => H reads input => output; see Pair arguments.

The result is the compression of the tensor product onto H. Unless projection=true, throws an ArgumentError if the tensor product of matrices maps a state of H outside H with nonzero amplitude, or if the tensor product of vectors has a nonzero amplitude outside H.

hilbert_space(b::BosonSym, dimension[, constraint])
hilbert_space(b::BosonField, labels, dimension[, constraint])

The truncated space of the bosonic mode b, with basis states |0⟩, …, |dimension-1⟩, so the maximal occupation is dimension - 1. With labels, the tensor product of one such space for each mode b[l], l in labels. An optional constraint restricts the basis states.

Operator terms that would raise the occupation above the cutoff are dropped without error, while leaving a constrained space throws unless projection=true. Products are normal ordered before truncation, so b * b' is represented as b'b + 1.

hilbert_space(f, labels, [constraint])
hilbert_space(b, labels, dimension, [constraint])
hilbert_space(s, labels, [J], [constraint])

Construct the Hilbert space of the fermions f, bosons b (truncated at dimension) or spins s (spin J) with the given labels, optionally restricted by a constraint, e.g. NumberConservation(2).

julia> @fermions f; @bosons b; @spins s 1//2;julia> dim(hilbert_space(f, 1:3, NumberConservation(1)))3julia> dim(hilbert_space(b, 1:2, 3))9julia> dim(hilbert_space(s, 1:2))4

indices(Hsub, H)
indices(qn, H)

Return the basis-state indices in H belonging to a given sector, specified either by a sector Hilbert space Hsub or by a quantum number qn, given as for sector. Throws an ArgumentError if qn is not one of quantumnumbers(H).

julia> using FermionicHilbertSpaces; @fermions f;julia> indices(1, hilbert_space(f, 1:2, NumberConservation()))2:3

matrix_representation(op, space::AbstractHilbertSpace, [repr]; kwargs...)
matrix_representation(op, Axes(Hout, Hin), [repr]; kwargs...)

Return the matrix representation of symbolic operator op in Hilbert space space. This is the operator-specific form of representation.

With Axes(Hout, Hin), return the dim(Hout) × dim(Hin) matrix of op from Hin to Hout. The two spaces must have the same groups, for example two sectors of the same modes. Two different spaces support only eager formats without chunking.

The optional repr selects the format, as in representation: :sparse (the default), :dense, :lazy, or a format object (EagerSparseRepr(), EagerDenseRepr(), LazyRepr()). A numeric element type also works as repr: ComplexF64 gives a sparse matrix with that element type, and Matrix{Float64} gives a dense matrix with that element type.

Keyword arguments include projection (default false) which if true will ignore the error thrown when the operator maps outside the space (outside Hout with Axes). Projection drops matrix elements to states outside the space, so an operator that changes a conserved quantity gives zero on a single sector; use Axes(Hout, Hin) with the two sectors instead. Use chunking to control NCAdd construction strategy: NoChunking(), TermChunking(scheduler), or StateChunking(scheduler).

Examples

@fermions fH = hilbert_space(f, 1:2)op = f[1]' * f[2] + hcM = FermionicHilbertSpaces.matrix_representation(op, H)size(M) == (dim(H), dim(H))

numberoperator(H)

Return the number operator Σ_j n_j for the Hilbert space H.

open_system(sym, args...; kwargs...)
open_system(H; kwargs...)

Set up the doubled space for density matrices of a system with Hilbert space H, or of hilbert_space(sym, args...; kwargs...). For H, kwargs are passed to tensor_product. Returns (Hlr, Hl, Hr, left, right):

• Hl and Hr are copies of the system space tagged :left and :right. Hr is transposed, so that representation(right(A * B), Hr) == representation(right(B), Hr) * representation(right(A), Hr).
• Hlr = tensor_product(Hl, Hr) is the space of vectorized density matrices.
• left(op) and right(op) map a symbolic operator of the system to the copy acting on Hl or Hr, so left(A) * right(B) represents ρ -> A ρ B and -1im * (left(H) - right(H)) is the Liouvillian of a Hamiltonian H.

Use reshape(ρ, (Hl, Hr) => Hlr) and reshape(v, Hlr => (Hl, Hr)) to convert between density matrices and vectors.

julia> using FermionicHilbertSpaces: dissipator; @fermions f;julia> Hlr, Hl, Hr, left, right = open_system(f, 1:1);julia> n = f[1]' * f[1];julia> L = representation(-1im * (left(n) - right(n)) + dissipator(f[1]), Hlr);julia> size(L)(4, 4)

parityoperator(H)

Return the fermionic parity operator (-1)^N for the Hilbert space H: 1 on even and -1 on odd basis states.

partial_trace(m, H => Hsub; complement=complementary_subsystem(H, Hsub))
partial_trace(H => Hsub; kwargs...)  # -> PartialTraceMap

Compute the partial trace of m from H to Hsub. Fermionic phase factors are included if both H and Hsub are Fermionic, unless specified otherwise in kwargs. H => Hsub reads input => output; see Pair arguments. A vector m is treated as a pure state ψ, and the result is the reduced density matrix of ψψ†.

The trace drops the states of H whose part on Hsub is not in Hsub, or whose part on the complement is not in complement. Unless projection=true, throws an ArgumentError if m has a nonzero element in the row or column of such a state. The map partial_trace(H => Hsub) does not check its input.

julia> using FermionicHilbertSpaces; @fermions f;julia> H = hilbert_space(f, 1:2); Hsub = hilbert_space(f, 1:1);julia> partial_trace(representation(f[1]' * f[1], H), H => Hsub) == [0 0; 0 2]true

partial_trace(H => Hsub; kwargs...)

Create a callable partial-trace operation from H to Hsub.

The operation can be applied directly as op(m) and can be converted to a sparse matrix map with sparse(op). It skips the states of H that partial_trace drops and, unlike partial_trace, does not check its input.

partial_transpose(m, H, Hsub; complement=complementary_subsystem(H, Hsub), phase_factors=true, projection=false)
partial_transpose(m, H => Hout, Hsub; kwargs...)
partial_transpose(H, Hsub; kwargs...)  # -> PartialTransposeMap
partial_transpose(H => Hout, Hsub; kwargs...)  # -> PartialTransposeMap

Partial transpose of m (an operator on H, a pure state, or a vectorized density matrix) on the subsystem Hsub, which may be any subset of the modes of H (contiguous or not). The result is an operator on Hout, which defaults to H. H => Hout reads input => output; see Pair arguments.

Without m, a callable PartialTransposeMap is returned. It takes the same keyword arguments.

For fermions this is the fermionic partial transpose (partial time reversal) R_A(Ẽ^{ν,ν'}) = i^{k_A} Ẽ^{(ν'_A,ν_B),(ν_A,ν'_B)} in the fermionic operator basis, which is independent of the mode ordering. With phase_factors=false, or for spaces without phase factors, it is the ordinary partial transpose.

The partial transpose does not preserve number or parity constraints. For constrained H, pass a larger output space on the same modes as H => Hout. For fermions the result is not Hermitian: (ρ^{R_A})' == p_A ρ^{R_A} p_A for parity-even ρ, with p_A the parity of Hsub.

particle_number(state)

Return the number of particles in the basis state state: the number of occupied modes for fermions, the occupation for a boson, and the total over the factors of a product state. It is not defined for spins. This is the quantity NumberConservation conserves.

quantumnumbers(H)

Return the quantum-number labels that define the sectors of H, sorted in increasing order when they can be compared. For parity this is [-1, 1]: odd before even. For spaces without sector structure, this returns (nothing,).

representation(op_or_state, space::AbstractHilbertSpace, [repr]; kwargs...)
representation(op, Axes(Hout, Hin), [repr]; kwargs...)
representation(state_expression; kwargs...)

Return a concrete representation of op_or_state in Hilbert space space. This is the main entry point, which calls matrix_representation for operators and vector_representation for states.

• For symbolic operators and numbers, returns a sparse matrix by default. The optional repr selects the format: :dense gives a dense matrix, and :lazy (or LazyRepr()) gives a matrix-free SciMLOperators.FunctionOperator that wraps a LazyOperator. A number c gives c times the identity.
• For AbstractBasisState and strings such as "10", returns a one-hot sparse column vector.
• For SymbolicState kets and bras, and sums of them, returns the corresponding column or row vector.
• For ketbras, and operators that contain ketbras, returns a matrix.

Expressions that combine kets or bras with operators or with states of other spaces have no representation and throw an ArgumentError. Without the space argument, only expressions of symbolic states are accepted, and the space is taken from the states. The single-particle space hilbert_space(f, labels, FreeFermions.SingleParticle()) and the Bogoliubov-de Gennes space hilbert_space(f, labels, FreeFermions.BdG()) represent quadratic fermionic operators only, and a number is represented by the zero matrix.

With Axes(Hout, Hin), op is represented as a map from Hin to Hout.

julia> @fermions f;julia> representation(f[1]' * f[1], hilbert_space(f, 1:1)) == [0 0; 0 1]true

sector(qn, H)

Return the sector of H corresponding to quantum number qn, for example a particle number for NumberConservation or 1 (even) / -1 (odd) for ParityConservation. For qn === nothing, this returns H for non-sector spaces. The quantum number of a product of constraints is a vector, which can also be given as a tuple: sector((1, -1), H). Throws an ArgumentError if qn is not one of quantumnumbers(H).

julia> using FermionicHilbertSpaces; @fermions f;julia> dim(sector(-1, hilbert_space(f, 1:3, ParityConservation())))4

sectors(H)

Return all sectors of H in the same order as quantumnumbers(H).

state_index(state, H)

Return the index of the basis state state in the Hilbert space H, the inverse of basisstate. A string or character that spells a state of H is accepted in place of the state. Returns 0 for a state of the same kind as the states of H (the same type up to type parameters) that is not in H, and throws a MethodError for a value of another kind.

Spaces implement the lookup of a basis state as a method of _state_index(state, H).

subregion(Hs, H::AbstractHilbertSpace)

Return the subsystem of H spanned by the factors Hs, given as a space, a symbolic operator, or a list/tuple of spaces or symbolic operators (e.g. [f[1], b[2]]).

If Hs contains only symbolic operators, their order does not matter: the modes are ordered as in H. If any element of Hs is a space, the modes must keep their order within each group of H (e.g. the Jordan-Wigner order of fermions), otherwise an ArgumentError is thrown.

Examples

julia> using FermionicHilbertSpaces; @fermions f;julia> H = hilbert_space(f, 1:3);julia> subregion([f[3], f[1]], H) == hilbert_space(f, [1, 3])true

symmetric_sector(H, Hs, sector=:symmetric, T=default; orth_method=:qr, check=true, normalize=true, sparse=true)

Return a basis for the sector of H selected by sector, under permutations of the factors Hs.

sector is either a tuple (perms, weights), with permutations of 1:length(Hs) and one weight per permutation, or :symmetric or :antisymmetric, which stand for all permutations with weight one or the permutation sign. T is the element type of the result. It defaults to Float64 for :symmetric and :antisymmetric, and for a tuple (perms, weights) to the floating point version of the element type of weights, for example ComplexF64 for complex weights.

For :symmetric and :antisymmetric with Hs a partition of H and orth_method not nothing, the basis is built directly from the orbits of the basis states under permutations, with one column per orbit that contributes to the sector. The basis then has element type float(T), and the keywords normalize and sparse have no effect. In the remaining cases, the projector is built with permutation_projector(H, Hs, perms, weights, T; normalize, sparse), and :symmetric and :antisymmetric require Combinatorics.jl to be loaded. orth_method selects what is returned:

• :qr: a matrix with orthonormal columns spanning the sector, from a QR decomposition of the projector.
• :eig: a matrix with orthonormal columns spanning the sector, from the eigenvectors of the projector with nonzero eigenvalue.
• nothing: the projector itself.

:qr and :eig find the numerical rank with a tolerance relative to the scale of the projector, so they assume that the projector is proportional to an orthogonal projector. This holds when perms form a group and weights are a positive multiple of a character of it or of its complex conjugate.

With check=true and a tuple (perms, weights), symmetric_sector throws an ArgumentError if perms does not contain the identity permutation or its weight is not real and positive, and, for orth_method not nothing, if the projector is not Hermitian. check=false skips these checks.

Hs may be either a full partition of H or a proper subsystem list; in the latter case the sector is constructed on tensor_product(Hs...) and embedded back into H.

tensor_product(spaces...; constraint=NoSymmetry())
tensor_product(spaces; constraint=NoSymmetry())

Construct the composite Hilbert space from the spaces in spaces, with optional keyword argument constraint.

Spaces of the same group, such as fermions of one @fermions basis, merge into one factor where the group first appears. The first factor varies fastest.

julia> using FermionicHilbertSpaces; @fermions f;julia> H = tensor_product(hilbert_space(f, 1:1), hilbert_space(f, 2:2));julia> state_index("10", H)2

tensor_product(ms, Hs, H::AbstractHilbertSpace; kwargs...)
tensor_product(ms, Hs => H; kwargs...)
tensor_product(Hs => H; kwargs...)

Compute the ordered product of the fermionic embeddings of the matrices ms in the spaces Hs into the space H. If H is unconstrained, tensor_product((A, B), (HA, HB) => H) == embed(A, HA => H) * embed(B, HB => H). If H is constrained, the result is the ordered product on the unconstrained space restricted to H, as representation of the corresponding symbolic product. Terms that take a state of H outside H throw an ArgumentError unless projection=true, which drops them. A matrix on a constrained space in Hs acts as zero on states outside that space. For vectors, the result is the product state: if each ms[k] = ψₖ|0⟩ with ψₖ an operator on the modes of Hs[k], the result is ψ₁ψ₂⋯|0⟩. Swapping two odd states flips the sign. kwargs can be passed a bool phase_factors. The last form returns a function of ms. Hs => H reads input => output; see Pair arguments.

vector_representation(state, space, [repr]; kwargs...) -> vector or transposed vector

Return the vector of state in space. state is a basis state, a symbolic ket or bra, or a product or sum of symbolic kets or of symbolic bras. Kets give a column vector, and bras give the transpose of one. This is the state-specific form of representation. The optional repr selects the format: :sparse (the default), :dense, or a format object (EagerSparseRepr(), EagerDenseRepr()).

add!!(a, b)

Return a + b, reusing the storage of a when possible, so use the returned value. foldl(add!!, terms; init=0) sums many operators in time linear in their number, where sum takes quadratic time.

@boson b

Create a symbolic bosonic annihilation operator b. Use b' for the corresponding creation operator.

@bosons b

Create a bosonic basis b. Index into it to get bosonic mode operators: b[i] returns the annihilation operator for mode i, b[i]' for the creation operator.

@fermions a b ...

Create one or more fermion species with the given names. Indexing into fermions species gives a concrete fermion. Fermions in one @fermions block anticommute with each other, and commute with fermions in other @fermions blocks. a[l] is the annihilation operator of the mode labelled l, a[l]' its creation operator, and a[i, σ] is the mode labelled (i, σ).

Examples:

• @fermions a b creates two species of fermions that anticommute:
  • a[1]' * a[1] + a[1] * a[1]' == 1
  • a[1]' * b[1] + b[1] * a[1]' == 0
• @fermions a; @fermions b creates two species of fermions that commute with each other:
  • a[1]' * a[1] + a[1] * a[1]' == 1
  • a[1] * b[1] - b[1] * a[1] == 0

See also @majoranas.

@majoranas a b ...

Create one or more Majorana species with the given names. Returns the species, or a tuple of them if several names are given. Indexing into Majorana species gives a concrete Majorana. Majoranas in one @majoranas block anticommute with each other, and commute with Majoranas in other @majoranas blocks.

Majoranas are Hermitian and normalised to square to one: γᵢγⱼ + γⱼγᵢ = 2δᵢⱼ.

Examples:

• @majoranas a b creates two species of Majoranas that anticommute:
  • a[1] * a[1] == 1
  • a[1] * b[1] + b[1] * a[1] == 0
• @majoranas a; @majoranas b creates two species of Majoranas that commute with each other:
  • a[1] * a[1] == 1
  • a[1] * b[1] - b[1] * a[1] == 0

See also @fermions.

@spin s

@spin s J

Create a symbolic spin basis s for constructing spin operators such as s[:x], s[:y], s[:z], s[:+], and s[:-]. The spin J can be given here or later to hilbert_space.

The operators are spin operators, not Pauli matrices: s[:z] has eigenvalues -J, …, J, so for spin 1/2 s[:x] is σˣ/2, and s[:±] = s[:x] ± i s[:y]. The basis states of hilbert_space(s, J) are ordered by increasing s[:z] eigenvalue, starting at -J.

For J ≥ 1, == and iszero on symbolic expressions compare their stored form, which is not unique, so equal operators can compare unequal: iszero(s[:z]*s[:+]^2 - s[:+]^2) is false for J = 1. Matrix representations are exact.

@spins s 

@spins s J

Create a spin basis s. Index into it to get site-specific spin bases: s[1][:z] gives the z-operator on site 1. Operators and basis states follow the conventions of @spin.