Wick’s Theorem

The getting-started page covers using sequant::WickTheorem — building operator products and reading its normal-ordered/overlap output. This page covers its configuration surface and the contraction algorithm behind it, for anyone extending or debugging that algorithm.

Configuration

sequant::WickTheorem is built from a sequant::NormalOperatorSequence or a general sequant::Expr; either way an index that appears once in the input is external and a repeated one is a dummy, unless the Context names the external indices (see Canonicalization options). It is then tuned via a handful of fluent setters before calling compute():

• full_contractions(bool) (default true): full contractions only, versus all contractions including partial ones.

• use_topology(bool) (default true): when the input is an Expr, treats the Ops within a bra or ket of a NormalOperator that can be swapped as a symmetry of the input (see below), and NormalOperator objects attached to the same tensor label, as topologically equivalent, so that contractions related by this equivalence are not separately enumerated.

• set_nop_connections() / set_nop_avoided_connections(): force, or forbid, contraction between specific pairs of normal-operator ordinals. Under a Vacuum::MultiProduct vacuum a required connection is enforced after cumulant expansion rather than during the contraction recursion, since a cumulant connects operators that no pair does; an avoided pair is still rejected the moment a contraction between it is attempted (see below).

• set_nop_partitions() / set_op_partitions() / make_default_op_partitions(): declare explicit equivalence groups of normal operators, or of individual Ops, so that contractions related by permuting within a group are counted once with a combinatorial degeneracy factor rather than enumerated redundantly — the general form of what use_topology() infers automatically.

compute(count_only, skip_input_canonicalization) then applies the theorem and returns a Constant, Product, or Sum. It is not reentrant, but is optionally threaded internally (a Sum input has its summands processed concurrently, one WickTheorem instance per summand, merged into a shared accumulator under a mutex — see below); it throws if the input’s vacuum does not match the current sequant::Context’s.

The contraction algorithm

compute() (SeQuant/core/wick.impl.hpp) expands a general expression input into normal-operator-sequence form and fully canonicalizes it, unless skip_input_canonicalization is set; for a Sum it then recurses per summand in parallel before merging results. Canonicalization leaves an expression it already canonicalized alone (see Tensor Network Canonicalizer), so an input the caller has simplified costs nothing to canonicalize again; a Product the caller has not simplified is canonicalized topologically, which is slower than the rapid canonicalization it used to get. After canonicalizing the input, compute() scopes a context that registers null canonicalizers for normal operators; a mark recorded outside that context is not valid within it, and is valid again once compute() returns. The actual enumeration happens in compute_nontensor_wick/recursive_nontensor_wick: it walks pairs of Op s across the flattened operator sequence, left to right. Under full_contractions_ (the default) it only ever extends a contraction starting from the leftmost still-free operator — the standard recursive formulation of full-contraction enumeration; with it disabled, all pairs are considered.

Because canonicalizing the produced normal operators would undo their normalization, a top-level compute() installs, after canonicalizing its input, a scoped sequant::Context in which the normal-operator labels map to NullTensorCanonicalizer for every statistics; a canonicalizer registered by the user for those labels is shadowed for the duration, not removed. The parallel workers of the per-summand WickTheorem instances see that scope.

A candidate pair (left, right) contracts (can_contract) iff left is a quasiparticle annihilator, right is a quasiparticle creator, and their quasiparticle spaces intersect (is_qpannihilator/is_qpcreator/IndexBasisRegistry::intersection, from SeQuant/core/op.hpp); under Vacuum::MultiProduct their actions must also differ (see below). The contraction value (contract) depends on whether those quasiparticle spaces are pure hole/particle subspaces or not:

• both pure: the contraction is a plain overlap \(s(L,R)\);

• neither pure: \(\delta(L,l)\, s(l,r)\, \delta(R,r)\);

• only the left space pure: \(s(L,r)\, \delta(R,r)\);

• only the right space pure: \(\delta(L,l)\, s(l,R)\),

where \(l = L \cap H\) and \(r = R \cap P\) (\(H\)/\(P\) the hole/particle subspaces), materialized via temporary indices from Index::make_tmp_index when a space needs projecting. The middle factor comes from WickTheorem::contraction_value, which returns \(s\) except under Vacuum::MultiProduct (see below). For fermionic statistics, each contraction additionally picks up a sign of \(-1\) for every operator that sat strictly between left and right at the time of contraction — the usual anticommutation rule.

Topological pruning and degeneracy

When use_topology_ is set, is_topologically_unique() restricts contraction attempts within a partition group (from use_topology()’s automatic grouping, or an explicit set_op_partitions()/set_nop_partitions()) to the group’s first free member, skipping contractions that would only reproduce one already tried by symmetry. The automatic grouping puts two Ops in one partition only if swapping them alone — their index vertices together with the tensor slots these occupy, every other vertex of the tensor network graph fixed — is an automorphism; an automorphism that also moves other indices, tensors or operators does not make the Ops independently interchangeable. The corresponding combinatorial weight — equivalent to a multinomial coefficient over how many operators from each partition have already been paired off — is recovered by op_permutational_degeneracy(), so the pruned enumeration and the exhaustive one agree on the final coefficient. They also agree for an input that vanishes by symmetry, i.e. has an automorphism of phase -1 (see Tensor Network Canonicalizer), because such an input is returned as zero before any contraction is attempted: by the up-front canonicalization of a Sum input, or by the topology analysis for a Product.

For the spin-free case over a fermionic vacuum, an additional factor of \(2^{n}\) is folded in per completed contraction, where \(n\) is the number of creation/annihilation “partner” cycles formed by that contraction — the generalized Wick’s theorem for spin-free operators (attributed, in the surrounding comment, to Kutzelnigg).

The extended Wick theorem (Vacuum::MultiProduct)

Relative to a general (multiconfigurational) reference the quasiparticle picture is lost: an operator on a partially occupied (active) orbital is both a quasiparticle creator and annihilator. SeQuant follows Kutzelnigg and Mukherjee’s generalized normal order: strings are normal-ordered so that their reference expectation value vanishes, and the theorem acquires, besides the pair contractions (now valued \(\gamma\) for cre·ann and \(\eta = \delta - \gamma\) for ann·cre), multi-leg contractions of \(k \ge 2\) creators and \(k\) annihilators valued by the \(k\)-body density cumulant \(\kappa_k\).

The reference must commute with the number operator (a state of definite particle number, or an ensemble of such states): only then does every string with unequal numbers of creators and annihilators have a vanishing reference expectation value and cumulant, so that the blocks above are the only multi-leg contractions (a reference that superposes particle numbers, such as a Bogoliubov vacuum, would need anomalous contractions such as \(\langle a_p a_q \rangle\) as well). Kutzelnigg and Mukherjee define generalized normal order for number-conserving strings; SeQuant extends it to every string by the same rule, under which a string is the sum, over all sets of disjoint internal contractions (pairs and cumulant blocks), of their values times the normal-ordered remainder. Hence \(\{a_p\} = a_p\), the reference expectation value of every normal-ordered string vanishes, and the theorem holds for products of such strings; the GNO products match the core-vacuum path test in tests/unit/test_mbpt.cpp checks it, built from that rule, against the standard theorem under the core vacuum.

Rather than a second engine, this is layered on the standard one, which changes in three places:

• the classifiers in SeQuant/core/op.hpp gain a MultiProduct branch in which the “hole” space is the registry’s reference-occupied space \(R\) and the “particle” space its vacuum-unoccupied space \(U\); the active space is \(R \cap U\), and an active Op is classified as both. can_contract therefore also requires opposite actions, so that only cre·ann (over \(R\)) and ann·cre (over \(U\)) contract;

• WickTheorem::contraction_value returns γ (left a creator) or η (left an annihilator), from SeQuant/core/density.hpp, as the middle factor of contract. A γ over \(R\) thus stands for \(\delta\) on the core plus \(\gamma\) on the active space, and likewise an η over \(U\) for \(\delta\) on the virtual space plus \(\eta\) on the active space. An index with protoindices must not reach the active space (nonorthogonal active orbitals are not supported);

• the engine skips the required-connectivity filter (WickTheorem::pairwise_connectivity()), since a cumulant can connect operators that no pair does; it keeps rejecting contractions between avoided pairs early, because no cumulant can undo a direct contraction, and cumulant_expand filters the cumulant-mediated connections of both kinds.

Under this vacuum WickTheorem::compute hands its input and options to detail::extended_wick (SeQuant/core/wick_extended.hpp), which drives the rest; count_only and bosonic statistics are not supported there. For each input term it applies the term’s own \(\delta\)s and overlaps over summed indices, since identifying two indices is not a contraction (one may leave an index shared by two operators, and one left carries the indices of at most one operator), canonicalizes, records which input NormalOperator every operator index comes from (its provenance; an index shared by two operators is renamed apart, the \(\delta\) binding the names multiplying the result), and runs the standard theorem (the private WickTheorem::compute_contractions) with partial contractions on the bare operator sequence, multiplying the c-number factors back in afterwards. In that run every operator index is external, so no term is canonicalized and no surviving index is renamed, which keeps the provenance valid. Topological equivalence needs dummy indices, though, so with use_topology the partitions are computed on the canonicalized term before its shared indices are renamed (renaming leaves every operator at its ordinal) and passed to that run with set_nop_partitions()/set_op_partitions(); the statistics of the run accumulate in the dispatching WickTheorem’s stats(). The terms of a Sum input are handled one after another, each with its own topology analysis, where the standard path contracts them in parallel. Each resulting term then

• has its mixed-space γ/η and surviving operators split into pure pieces: a γ into a core \(\delta\) plus an active γ, an η into a virtual \(\delta\) plus an active η (a one-body γ or η of the input is split the same way, whatever space its indices range over; the \(\delta\) binds the bra to a core or virtual index in its basis, so between two bases, e.g. indices with different protoindices or in two basis instances, it is followed by an overlap; as for an operator, an index with protoindices must not reach the active part), a survivor onto its active part and, with partial contractions, its core and virtual parts. A part that is not a registered space is split over its base spaces. The pieces are reduced with the operator indices kept fixed, so each projected index stays \(\delta\)-bound to an input index and inherits its provenance;

• is handed to detail::cumulant_expand, which groups the surviving active operators into disjoint blocks of \(k\) creators and \(k\) annihilators, \(2 \le k \le\) WickTheorem::max_cumulant_rank, with legs from at least two input operators. Each block is valued \(\kappa_k\) (detail::block_value), the term takes the parity of the permutation that pulls each block’s legs, in order, to the front of the survivor string, and any remaining operators (partial contractions only) stay as a MultiProduct-vacuum NormalOperator. Terms that miss a required nop_connections pair or realize a nop_avoided_connections pair are dropped; a factor the theorem produces (a \(\gamma\), \(\eta\), \(\kappa\), \(\delta\) or overlap) connects all the input operators its indices came from, any other tensor (e.g. a coefficient) none.

Finally every one-body η is optionally rewritten as \(\delta - \gamma\) (WickTheorem::eta_as_delta_minus_gamma; a multi-body η of the input is kept, and the \(\delta\) is an overlap unless that is a Kronecker delta: between indices with different protoindices, in different basis instances or in a non-orthonormal basis, see sequant::is_kronecker_equivalent()), the \(\delta\)s over summed indices are applied, and the result is simplified. In it γ{ann;cre} and η{ann;cre} are one-body, Hermitian and column-symmetric, κ{ann…;cre…} is of rank \(\ge 2\), antisymmetric, Hermitian and column-symmetric, and all their indices are active. These symmetries take part in the tensor hash, so a density spelled differently would not merge with the engine’s; hence γ, η, κ and the spin-free Γ are reserved labels (reserved::density_labels()): a Tensor carrying one must have the symmetries density::symmetries() gives for its label and rank, or its constructor (and set_label) throws, and the parser supplies them; the permutational symmetry of a one-body density is void, so any spelled-out one is accepted and normalized away. The one alternative is a perm-nonsymmetric multi-body γ, η or κ, a spin component such as spin tracing produces. The factories in SeQuant/core/density.hpp build densities with them.

The no-double-counting rule is that a block is only ever built from operators that survive the standard theorem: each extended term descends from exactly one partial-contraction term (the one carrying its pair contractions), so no \(1/k!\) weights are needed. Spin-free evaluation is not supported for this vacuum (WickTheorem::compute throws); the hooks for it are WickTheorem::contraction_value and detail::block_value/detail::term_weight in SeQuant/core/wick_extended.hpp.

mbpt::ref_av has no branch of its own for this vacuum, only forcing full contractions and so accepting connectivity lists, under which a density or cumulant carrying indices of two operators connects them. There the mbpt operators (ã) are normal-ordered relative to the reference, so it evaluates a different quantity than the SingleProduct path, which normal-orders them relative to the core; the two agree for products of elementary operators, which the tests verify for results with cumulants up to \(\kappa_4\). sequant::mbpt::decompositions::cumulants_to_densities() converts cumulants of every rank to densities. The tests are in tests/unit/test_wick_extended.cpp, tests/unit/test_wick.cpp and the MRSO-MultiProduct section of tests/unit/test_mbpt.cpp.

Reducing the result

Raw contraction output can contain chains of Kronecker deltas and overlaps introduced by the space-projection cases above. WickTheorem::reduce() turns these into an index-replacement map and substitutes it through the expression, collapsing delta chains and eliminating deltas wherever the internal/external status of the indices they bind allows it. An overlap is applied only if one of its indices is covariant, a dummy that appears exactly twice and neither has nor is a protoindex, since such an index can be rotated into the basis of the other; between two noncovariant indices it stands. The exception is an overlap that is itself a Kronecker delta (sequant::is_kronecker_equivalent(): both indices in one orthonormal basis, i.e. with the same protoindices, not in two different basis instances and neither in a basis of general metric): it identifies its indices, through the indices they are protoindices of as well, so it is applied unless both are unpaired dummies (not appearing exactly twice); unlike a Kronecker delta between them, which is applied, it then stands. These rules apply only if the input of Wick’s theorem has noncovariant indices; otherwise every overlap is applied. The result is then, like any other sequant::Product/sequant::Sum, put into canonical form by the tensor-network canonicalizer so that like terms collect correctly — Wick’s-theorem correctness therefore also rests on canonicalization being correct. The reduction treats a basis instance (sequant::IndexBasis) like a subspace: an index without one runs over the space’s own basis, which includes every instance of the space, so a delta or overlap between a generic and a specific index replaces the generic one by the specific one, and a specific index is never replaced by a generic one. Two different instances are not related by an identity: an overlap between them (or between indices that the rules collected so far reduce to different instances) stands and yields no rule, and a Kronecker delta between them throws sequant::Exception, since basis functions of different bases cannot be compared for equality; either is zero if the spaces of its indices (as narrowed by the rules collected so far) are disjoint, whatever the bases. When a basis instance is present the deltas are applied first, then the overlaps, and an overlap between two generic indices last, once every index that will be specific is, so the result does not depend on the order of the factors.

Debugging and tests

Logger::instance().wick_harness, .wick_contract, and .wick_reduce (SeQuant/core/logger.hpp) trace, respectively, the top-level expand/canonicalize/recurse flow, individual contraction attempts, and the delta/overlap reduction step. The algorithm’s test coverage lives in tests/unit/test_wick.cpp, that of the basis-instance rules in tests/unit/test_basis_wick.cpp.