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)\]

where the 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)\]

where \(␣\) is used in number-nonconserving operators to point out the empty “slots”.

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\)