Wick’s Theorem¶

To get started let’s use SeQuant to apply Wick’s theorem to a simple product of elementary (creation and annihilation) fermionic operators:

\[a_{p_3} a_{p_4} a^\dagger_{p_1} a^\dagger_{p_2}.\]

This is achieved by the following SeQuant program:

Index p1(L"p_1"), p2(L"p_2"), p3(L"p_3"), p4(L"p_4");

auto cp1 = fcrex(p1), cp2 = fcrex(p2);
auto ap3 = fannx(p3), ap4 = fannx(p4);

std::wcout << to_latex(ap3 * ap4 * cp1 * cp2) << " = "
           << to_latex(FWickTheorem{ap3 * ap4 * cp1 * cp2}
                           .full_contractions(false)
                           .compute())
           << std::endl;

Hint

All core user-facing SeQuant code lives in C++ namespace sequant, and the shown code assumes this namespace has been imported via using namespace sequant.

Running this program should produce a LaTeX expression for this formula:

\[{{a_ {{p_3}}}{a_ {{p_4}}}{a^{{p_1}}}{a^{{p_2}}}} = \bigl( - {{a^{{p_1}{p_2}}_ {{p_3}{p_4}}}} - {{s^{{p_1}}_ {{p_3}}}{s^{{p_2}}_ {{p_4}}}} + {{s^{{p_1}}_ {{p_3}}}{a^{{p_2}}_ {{p_4}}}} - {{s^{{p_1}}_ {{p_4}}}{a^{{p_2}}_ {{p_3}}}} + {{s^{{p_2}}_ {{p_3}}}{s^{{p_1}}_ {{p_4}}}} - {{s^{{p_2}}_ {{p_3}}}{a^{{p_1}}_ {{p_4}}}} + {{s^{{p_2}}_ {{p_4}}}{a^{{p_1}}_ {{p_3}}}}\bigr)\]

Wick’s theorem reduces this to a sum of terms in which the elementary operators have been reordered so that every creation operator sits to the left of every annihilation operator — the normal order, or normal-ordered form — plus terms in which some operator pairs have been contracted into overlaps instead. Normal-ordered operators are useful because their vacuum expectation value vanishes trivially, which is why SeQuant expresses Wick’s-theorem output in terms of them. The following tensor notation is used to denote elementary and composite normal-ordered (or, shortly, normal) operators:

\[a^p \equiv a_p^\dagger\]

\[a^{p_1 p_2 \dots p_c}_ {q_1 q_2 \dots q_a} \equiv a_ {p_1}^\dagger a_ {p_2}^\dagger \dots a_ {p_c}^\dagger a_ {q_a} \dots a_ {q_2} a_ {q_1}\]

\(s^p_q \equiv \langle q | p \rangle\) denotes inner products (“overlaps”) of 1-particle states. Wick’s theorem can of course be applied directly to products of normal composite operators, e.g,

auto nop1 = ex<FNOperator>(cre(p1, p2), ann(p3, p4));
// OR
// auto nop1 = ex<FNOperator>(cre({p1, p2}), ann({p3, p4}));
// OR
// auto nop1 = ex<FNOperator>(cre(std::vector{p1, p2}),
//                            ann(std::array{L"p3", L"p4"}));
// OR
// auto nop1 = ex<FNOperator>(cre(std::set{"p1", "p2"}),
//                            ann(std::vector{L"p3", L"p4"}));
auto nop2 = ex<FNOperator>(cre({L"p_5"}), ann({L"p_6", L"p_7"}));

std::wcout
    << to_latex(nop1 * nop2) << " = "
    << to_latex(FWickTheorem{nop1 * nop2}.full_contractions(false).compute())
    << std::endl;

produces

\[{{a^{{p_1}{p_2}}_ {{p_3}{p_4}}}{a^{␣\,{p_5}}_ {{p_6}{p_7}}}} = \bigl({a^{␣\,{p_1}{p_2}{p_5}}_ {{p_3}{p_4}{p_6}{p_7}}} - {{s^{{p_5}}_ {{p_4}}}{a^{␣\,{p_1}{p_2}}_ {{p_3}{p_6}{p_7}}}} + {{s^{{p_5}}_ {{p_3}}}{a^{␣\,{p_1}{p_2}}_ {{p_4}{p_6}{p_7}}}}\bigr)\]

A number-nonconserving operator has an unequal number of upper (creation) and lower (annihilation) indices; \(␣\) marks the resulting empty “slots” that would otherwise misalign the tensor’s rows.

Same algebra can be performed for bosons:

auto nop3 = ex<BNOperator>(cre({p1, p2}), ann({p3, p4}));
auto nop4 = ex<BNOperator>(cre({L"p_5", L"p_6"}), ann({L"p_7"}));

std::wcout
    << to_latex(nop3 * nop4) << " = "
    << to_latex(BWickTheorem{nop3 * nop4}.full_contractions(false).compute())
    << std::endl;

\[{{b^{{p_1}{p_2}}_ {{p_3}{p_4}}}{b^{{p_5}{p_6}}_ {␣\,{p_7}}}} = \bigl( {b^{{p_5}{p_1}{p_2}{p_6}}_ {␣\,{p_3}{p_4}{p_7}}} + {{s^{{p_6}}_ {{p_3}}}{b^{{p_5}{p_2}{p_1}}_{␣\,{p_4}{p_7}}}} + {{s^{{p_5}}_{{p_3}}}{b^{{p_1}{p_2}{p_6}}_{␣\,{p_4}{p_7}}}} + {{s^{{p_6}}_{{p_4}}}{b^{{p_5}{p_1}{p_2}}_{␣\,{p_3}{p_7}}}} + {{s^{{p_5}}_{{p_4}}}{b^{{p_2}{p_1}{p_6}}_{␣\,{p_3}{p_7}}}} + {{s^{{p_5}}_{{p_3}}}{s^{{p_6}}_{{p_4}}}{b^{{p_1}{p_2}}_{␣\,{p_7}}}} + {{s^{{p_6}}_{{p_3}}}{s^{{p_5}}_{{p_4}}}{b^{{p_2}{p_1}}_{␣\,{p_7}}}} \bigr)\]

where \(b\) denotes normal bosonic operators constructed analogously with the normal fermionic operators \(a\).

Extended Wick’s theorem¶

With Vacuum::MultiProduct the reference is a general (e.g. multiconfigurational) state that commutes with the number operator, and operators, which need not conserve particle number, are normal-ordered relative to it (generalized normal order). The same WickTheorem then evaluates products of such operators into one-body densities \(\gamma\), one-body hole densities \(\eta\) and density cumulants \(\kappa_k\) over the active (partially occupied) orbitals; max_cumulant_rank bounds the rank of the cumulants and eta_as_delta_minus_gamma spells \(\eta\) as \(\delta - \gamma\):

// a multireference vocabulary of index spaces: core (o, i), active (u),
// virtual (a, g); the reference is a general state in the active space
set_default_context(
    Context({.index_basis_registry_shared_ptr = mbpt::make_mr_spaces(),
             .vacuum = Vacuum::MultiProduct}));

// a product of two generalized-normal-ordered one-body operators
auto expr = ex<FNOperator>(cre({L"u_1"}), ann({L"u_2"})) *
            ex<FNOperator>(cre({L"u_3"}), ann({L"u_4"}));

// its reference expectation value: a γη pair plus a 2-body cumulant κ
auto vev = FWickTheorem{expr}.compute();
std::wcout << to_latex(vev) << std::endl;

// its generalized-normal-ordered form, cumulants truncated at rank 2
auto gno = FWickTheorem{expr}
               .full_contractions(false)
               .max_cumulant_rank(2)
               .compute();
std::wcout << to_latex(gno) << std::endl;

This prints the reference expectation value

\[\bar{\kappa}^{u_1 u_3}_{u_2 u_4} + \gamma^{u_1}_{u_4} \eta^{u_3}_{u_2}\]

followed by the generalized-normal-ordered form of the product, in which \(\tilde{a}\) denotes operators normal-ordered relative to the reference:

\[\tilde{a}^{u_1 u_3}_{u_2 u_4} - \gamma^{u_1}_{u_4} \tilde{a}^{u_3}_{u_2} + \eta^{u_3}_{u_2} \tilde{a}^{u_1}_{u_4} + \bar{\kappa}^{u_1 u_3}_{u_2 u_4} + \gamma^{u_1}_{u_4} \eta^{u_3}_{u_2}\]

The labels γ, η and κ are reserved: build such tensors with the factories in SeQuant/core/density.hpp (or parse them, e.g. γ{u_1;u_2}), which give them their fixed symmetries. Note that parsing reads the context in effect: a parsed ã takes its vacuum, and a parsed density takes the bra-ket symmetry its field implies, so parse an expression under the context it is meant for.