Second quantization#

Fock space#

Consider a system of electrons in \(M\) spin orbitals. An occupation number vector lists the spin-orbital occupation

\[ | \mathbf{k} \rangle = | k_1 \ldots k_p\ldots k_M \rangle \]

where the occupation numbers, \(k_p\), are restricted to 0 or 1.

The Fock space, \(F(M)\), is defined as the vector space of all possible occupation vectors.

There is a unique mapping between occupation number vectors and Slater determinants, but the Fock space includes states with different number of electrons.

The Fock space is a linear vector space with the following properties:

  • occupation number vectors are orthonormal

    \[ \langle \mathbf{k} | \mathbf{m} \rangle = \prod_{p=1}^{M} \delta_{k_p m_p} \]
  • the identity operator can be written as

    \[ \hat{I} = \sum_\mathbf{k} | \mathbf{k} \rangle \langle \mathbf{k} | \]
  • it can be partitioned into sub-spaces according to

    \[ F(M) = F(0,M) \oplus F(1,M) \oplus \cdots \oplus F(M,M) \]

    where \(F(N, M)\) is the set of occupation number vectors with \(N\) occupied spin-orbitals.

  • subspace \(F(0, M)\) corresponds to a completely unoccupied Slater determinant known as the vacuum state, \(| \mathrm{vac} \rangle\). For the vacuum state, we have

    \[\begin{equation*} | \mathrm{vac} \rangle = | 0_1 0_2 \ldots 0_M \rangle ; \qquad \langle \mathrm{vac} | \mathrm{vac} \rangle = 1 \end{equation*}\]

Creation and annihilation operators#

Second quantization introduces a convenient way to manipulate states by means of creation and annihilation operators. As their name imply, these operators create and destroy electrons in specific spin-orbitals as follows

  • creation operator \(\hat{a}_p^\dagger\):

\[ \hat{a}_p^\dagger | \mathbf{k} \rangle = \delta_{k_p 0} \prod_{i = 1}^{p-1} (-1)^{k_p} | k_1 k_2...1_p ...k_M \rangle \]
  • annihilation operator \(\hat{a}_p\):

\[ \hat{a}_p | \mathbf{k} \rangle = \delta_{k_p 1} \prod_{i = 1}^{p-1} (-1)^{k_p} | k_1 k_2...0_p ...k_M \rangle \]

where the factor \(\delta_{k_p 0}\) prevents the addition of an electron to a state that is already populated and \(\delta_{k_p 1}\) prevents the destruction of a non-existing electron.

The phase factor \(\prod_{i = 1}^{p-1} (-1)^{k_p}\) ensures the anti-symmetry of the associated wave function and it is a direct consequence of the anti-commutator relationships:

\[\begin{eqnarray*} \hat{a}_p^\dagger \hat{a}_q + \hat{a}_q \hat{a}_p^\dagger &=& \delta_{pq} \\ \hat{a}_p^\dagger \hat{a}_q^\dagger + \hat{a}_q^\dagger \hat{a}_p^\dagger &=& 0 \\ \hat{a}_p \hat{a}_q + \hat{a}_q \hat{a}_p &=& 0 \end{eqnarray*}\]

Ground state#

The ground state of an \(N\)-electron system can be formed from the vacuum state

\[ | 0 \rangle = \prod_{i = 1}^{N} \hat{a}^\dagger_{i} |\mathrm{vac}\rangle \]

Slater determinants#

There is a one-to-one correspondence between occupation number vectors and Slater determinants. For \(F(N,M)\), let us collect the \(N\) positions in a given occupation number vector that have nonzero occupations in an ordered index set \(\{i_1, i_2, \ldots,i_N\}\). We then have

\[ | \mathbf{k} \rangle = | k_1 \ldots k_p\ldots k_M \rangle \leftrightarrow | \psi_{i_1}, \psi_{i_2}, \ldots, \psi_{i_N}\rangle \]

and, out of convenience, we may at times denote the occupation number vector as equal to the associated Slater determinant.

We can write an electron excitation from occupied orbital \(i_p\) to unoccupied orbital \(a\) with direct reference to Slater determinants

\[ \hat{a}^\dagger_{a}\hat{a}_{i_p} | \psi_{i_1}, \ldots, \psi_{i_p}, \ldots, \psi_{i_N}\rangle = | \psi_{i_1}, \ldots, \psi_{a}, \ldots, \psi_{i_N}\rangle \]

where the phase factors associated with the annihilation and creation operators in general cancel at the same time as the energy ordering of orbitals in the resulting Slater determinant is lost.

The actions of individual creation and annihilation operators can also be defined for Slater determinants

\[\begin{align*} \hat{a}^\dagger_{a} | \psi_{1}, \psi_{2}, \ldots, \psi_{N}\rangle & = | \psi_{a}, \psi_{1}, \psi_{2}, \ldots, \psi_{N}\rangle \\ \hat{a}_{i} | \psi_{i}, \psi_{1}, \psi_{2}, \ldots, \psi_{N}\rangle & = | \psi_{1}, \psi_{2}, \ldots, \psi_{N}\rangle \\ \end{align*}\]

The explicit reference to the associated phase factors is also here avoided by a presumed ordering of orbitals in the Slater determinant.

Illustration for Fock space F(4,10)#

Let us consider the closed-shell ground state for a system with four electrons in five molecular orbitals, or ten spin orbitals,

\[ |0\rangle = |1 1 1 1 0 0 0 0 0 0\rangle = |\psi_1, \psi_\bar{1}, \psi_2, \psi_\bar{2}\rangle \]
import numpy as np

M = 10  # number of spin orbitals in our Fock space

We implement a class of occupation number vectors.

class OccupationNumberVector:
    def __init__(self, occ_orbs=[]):
        self.phase = 1
        self.occ_num_vec = np.zeros(M, dtype=int)  # vacuum state

        for i in occ_orbs:
            self.occ_num_vec[i] = 1

    def __mul__(self, other):
        if np.array_equal(self.occ_num_vec, other.occ_num_vec):
            inner_prod = self.phase * other.phase
        else:
            inner_prod = 0

        return inner_prod

    def __str__(self):
        return (
            "phase: "
            + str(self.phase)
            + "\noccupation number vector: "
            + str(self.occ_num_vec)
        )

    def creation(self, p):
        if len(self.occ_num_vec) == 0:
            return

        if self.occ_num_vec[p] == 1:
            self.occ_num_vec = []
        else:
            self.phase *= np.prod((-1) ** self.occ_num_vec[:p])
            self.occ_num_vec[p] = 1

    def annihilation(self, p):
        if len(self.occ_num_vec) == 0:
            return

        if self.occ_num_vec[p] == 0:
            self.occ_num_vec = []
        else:
            self.phase *= np.prod((-1) ** self.occ_num_vec[:p])
            self.occ_num_vec[p] = 0

The ground state can formed by applying a series of creation operators on the vacuum state

\[ | 0 \rangle = \prod_{i}^\mathrm{occ} \hat{a}^\dagger_i |\mathrm{vac}\rangle \]
ket = OccupationNumberVector()

for i in [3, 2, 1, 0]:
    ket.creation(i)

print(ket)
phase: 1
occupation number vector: [1 1 1 1 0 0 0 0 0 0]

One-electron operators#

A one-electron operator takes the form

\[ \hat{\Omega} = \sum_{p,q} \omega_{pq} \, \hat{a}_p^\dagger \hat{a}_q \]

where

\[ \omega_{pq} = \langle \psi_p | \hat{\omega} | \psi_q \rangle \]

One-particle density#

The one-particle density matrix becomes

\[ D_{pq} = \langle 0| \hat{a}_p^\dagger \hat{a}_q| 0\rangle \]
bra = OccupationNumberVector([0, 1, 2, 3])

D = np.zeros((M, M), dtype=int)

for p in range(M):
    for q in range(M):
        ket = OccupationNumberVector([0, 1, 2, 3])

        ket.annihilation(q)
        ket.creation(p)

        D[p, q] = bra * ket

print("One-particle denisty matrix:\n", D)
One-particle denisty matrix:
 [[1 0 0 0 0 0 0 0 0 0]
 [0 1 0 0 0 0 0 0 0 0]
 [0 0 1 0 0 0 0 0 0 0]
 [0 0 0 1 0 0 0 0 0 0]
 [0 0 0 0 0 0 0 0 0 0]
 [0 0 0 0 0 0 0 0 0 0]
 [0 0 0 0 0 0 0 0 0 0]
 [0 0 0 0 0 0 0 0 0 0]
 [0 0 0 0 0 0 0 0 0 0]
 [0 0 0 0 0 0 0 0 0 0]]

Two-electron operators#

A two-electron operator takes the form

\[ \hat{\Omega} = \sum_{p,q,r,s} \omega_{pqrs} \, \hat{a}_p^\dagger \hat{a}_q^\dagger \hat{a}_s \hat{a}_r \]

where, in the physicist’s notation,

\[ \omega_{pqrs} = \langle \psi_p \psi_q | \hat{\omega} | \psi_r \psi_s \rangle \]

Hamiltonian#

The Hamiltonian becomes

\[ \hat H = \sum_{p,q} h_{pq} \hat{a}_p^\dagger \hat{a}_q + \frac{1}{2} \sum_{p,q,r,s} g_{pqrs} \hat{a}_p^\dagger \hat{a}_r^\dagger \hat{a}_s \hat{a}_q \]

where, in the chemist’s notation,

\[ g_{pqrs} = (\phi_p \phi_q | \hat{\omega} | \phi_r \phi_s ) \delta_{\sigma_p\sigma_q} \delta_{\sigma_r\sigma_r} \]

Here \(\phi_p\) is the spatial part (MO) of spin orbital \(\psi_p\) and the Kronecker deltas are the result of the spin integration.

Energy#

The energy of a single-determinant state equals

\[ \langle 0| \hat H | 0\rangle = \sum_{p,q} h_{pq} \langle 0| \hat{a}_p^\dagger \hat{a}_q| 0\rangle + \frac{1}{2} \sum_{p,q,r,s} g_{pqrs} \langle 0| \hat{a}_p^\dagger \hat{a}_r^\dagger \hat{a}_s \hat{a}_q | 0\rangle \]

We thus need to evaluate the one- and two-particle density matrices.

One-electron term#

For the one-electron term, we notice that orbital \(q\) has to be occupied to obtain a nonzero result as the annihilation operator acts to the right on \(| 0\rangle\). Using the anti-commutator relation, we get

\[ \langle 0| \hat{a}_p^\dagger \hat{a}_i | 0 \rangle = \delta_{pi} \langle 0| 0 \rangle - \langle 0| \hat{a}_i \hat{a}_p^\dagger | 0 \rangle = \delta_{pi} \]

This results in

\[ \sum_{p,q} h_{pq} \langle 0| \hat{a}_p^\dagger \hat{a}_q| 0\rangle = \sum_i^\mathrm{occ} h_{ii} \]

Two-electron term#

bra = OccupationNumberVector([0, 1, 2, 3])

for p, q, r, s in [(0, 0, 1, 1), (0, 1, 1, 0), (0, 1, 0, 1)]:
    ket = OccupationNumberVector([0, 1, 2, 3])

    ket.annihilation(q)
    ket.annihilation(s)
    ket.creation(r)
    ket.creation(p)

    print("Two-particle density matrix element:", bra * ket)
Two-particle density matrix element: 1
Two-particle density matrix element: -1
Two-particle density matrix element: 0

For the two-electron term, we notice that orbitals \(q\) and \(s\) both have to be occupied to obtain a nonzero result. We focus on \(s\) and will use the same trick to move \(\hat{a}_i\) to the left

\[\begin{align*} \hat{a}_p^\dagger \hat{a}_r^\dagger \hat{a}_i \hat{a}_q & = \delta_{ri} \hat{a}_p^\dagger \hat{a}_q - \hat{a}_p^\dagger \hat{a}_i \hat{a}_r^\dagger \hat{a}_q \\ &= \delta_{ri} \hat{a}_p^\dagger \hat{a}_q - \delta_{pi} \hat{a}_r^\dagger \hat{a}_q + \hat{a}_i \hat{a}_p^\dagger \hat{a}_r^\dagger \hat{a}_q \end{align*}\]

The last term will vanish in the expectation value and we are left with only terms involving only pairs of creation or annihilation operators. We get

\[\begin{align*} \sum_{p,q,r,s} g_{pqrs} \langle 0| \hat{a}_p^\dagger \hat{a}_r^\dagger \hat{a}_s \hat{a}_q | 0\rangle &= \sum_{i,p,q} g_{pqii} \langle 0|\hat{a}_p^\dagger \hat{a}_q | 0\rangle - \sum_{i,q,r} g_{iqri} \langle 0| \hat{a}_r^\dagger \hat{a}_q | 0\rangle \\ & = \sum_{i,j} \big( g_{jjii} - g_{ijji} \big) \end{align*}\]

and the total energy thus becomes equal to

\[ \langle 0| \hat H | 0\rangle = \sum_i h_{ii} + \frac{1}{2} \sum_{i,j} \big( g_{jjii} - g_{ijji} \big) \]