Single Particle Fermion Operators
Multiple Particles Fermion Operators
In first quantization, we transition from classical mechanics to quantum mechanics. The classical Hamiltonian H(\{\mathbf{r}_i\}, \{\mathbf{p}_i\}) and the equations of motion for a system of particles i = 1, \ldots, N become operators. The time-dependent Schrödinger equation is:
i\hbar \frac{\mathrm{d}}{\mathrm{d}t} \Psi = \mathbf{H} \Psi
where the Hamiltonian operator \mathbf{H} is obtained by replacing momentum \mathbf{p}_i with -i\hbar \nabla_i in the classical Hamiltonian. The wavefunction \Psi(\{\mathbf{r}_i\}, t) contains all information about the quantum N-particle system as it evolves in time. Observables are expressed as expectation values:
\langle \Psi(\{\mathbf{r}_i\}, t) | \mathbf{G} | \Psi(\{\mathbf{r}_i\}, t) \rangle
Second quantization reformulates quantum mechanics to handle many interacting particles more effectively. It uses creation (\mathbf a^\dagger for bosons and \mathbf b^\dag for fermions) and annihilation (\mathbf a and \mathbf b) operators to manage particle states, obeying commutation or anticommutation rules for bosons and fermions, respectively. The formulation of these operators is available here.
The formalism extends to the occupation number representation, where states are described by the number of particles in each quantum state. Field operators \hat{\psi}^\dagger(\mathbf{r}) and \hat{\psi}(\mathbf{r}) are introduced to create and annihilate particles at position \mathbf{r}.
Starting from a of a single-fermion state the wavefunction operator for a fermion system can be expressed using the creation (\mathbf{b}^\dagger) and annihilation (\mathbf{b}) operators. For a single fermion state, the wavefunction operator \hat{\psi}(\mathbf{r}) is postulated as:
\hat{\psi}(\mathbf{r}) = \sum_j \phi_j(\mathbf{r}) \mathbf{b}_j
Here, \phi_j(\mathbf{r}) are the single-particle wavefunctions (modes), and \mathbf{b}_j are the annihilation operators for the j^{th} mode.
In this formalism, the state of the system is described by the occupation numbers of the modes, and the creation and annihilation operators modify these occupation numbers, adhering to the Pauli exclusion principle which is fundamental for fermions.
Supposing the particle is in a state k, then it can be written as:
\mathbf b^\dagger_k | \mathbf 0 \rangle
so that is creating a fermion on the k^{th} position from the vacuum state.
The wavefunction operator \hat{\psi}(\mathbf{r}) can then be used to operate on this state:
\hat\psi(\mathbf r) \mathbf b^\dagger_j | \mathbf 0 \rangle = \left( \sum_j \phi_j(\mathbf{r}) \mathbf{b}_j \right) \mathbf{b}^\dagger_k | \mathbf{0} \rangle= \sum_j \phi_j(\mathbf{r}) \mathbf{b}_j \mathbf{b}^\dagger_k | \mathbf{0} \rangle
Using the anticommutation relations, specifically \{\mathbf{b}_j, \mathbf{b}^\dagger_k\} = \delta_{jk}, since \mathbf{b}_j \mathbf{b}^\dagger_k = \delta_{jk} - \mathbf{b}^\dagger_k \mathbf{b}_j and \mathbf{b}_i | \mathbf{0} \rangle = 0, we have:
\hat\psi(\mathbf r) \mathbf b^\dagger_j | \mathbf 0 \rangle = \phi_j(\mathbf{r}) | \mathbf{0} \rangle
Thus, the wavefunction operator for the state with a single fermion in the k^{th} mode is simply the single-particle wavefunction \phi_j(\mathbf{r}) acting on the vacuum state.
Continuing from the concept of a single-particle state and extending it to a superposition state, we can write a state that is a linear superposition of the single-particle states. Suppose the state is given by:
|\Psi\rangle = \sum_k \alpha_k \mathbf{b}^\dagger_k | \mathbf{0} \rangle
where \alpha_k are the coefficients of the superposition.
Now, let’s consider the action of the wavefunction operator \hat{\psi}(\mathbf{r}) on this superposition state. We have:
\hat{\psi}(\mathbf{r}) |\Psi\rangle = \hat{\psi}(\mathbf{r}) \left( \sum_k \alpha_k \mathbf{b}^\dagger_k | \mathbf{0} \rangle \right)
Substituting the expression for \hat{\psi}(\mathbf{r}) we get:
\hat{\psi}(\mathbf{r}) \sum_k \alpha_k \mathbf{b}^\dagger_k | \mathbf{0} \rangle = \sum_j \phi_j(\mathbf{r}) \mathbf{b}_j \sum_k \alpha_k \mathbf{b}^\dagger_k | \mathbf{0} \rangle
Using the anticommutation relations, \{\mathbf{b}_j, \mathbf{b}^\dagger_k\} = \delta_{jk}, and the fact that \mathbf{b}_i | \mathbf{0} \rangle = 0:
\sum_j \phi_j(\mathbf{r}) \mathbf{b}_j \sum_k \alpha_k \mathbf{b}^\dagger_k | \mathbf{0} \rangle = \sum_j \phi_j(\mathbf{r}) \sum_k \alpha_k (\delta_{jk} - \mathbf{b}^\dagger_k \mathbf{b}_j) | \mathbf{0} \rangle
Since \mathbf{b}_j | \mathbf{0} \rangle = 0, the second term vanishes:
\hat{\psi}(\mathbf{r}) |\Psi\rangle = \sum_j \phi_j(\mathbf{r}) \sum_k \alpha_k \delta_{jk} | \mathbf{0} \rangle
This simplifies to:
\hat{\psi}(\mathbf{r}) |\Psi\rangle = \sum_j \phi_j(\mathbf{r}) \alpha_j | \mathbf{0} \rangle
Thus, the wavefunction operator acting on the superposition state |\Psi\rangle gives:
\hat{\psi}(\mathbf{r}) |\Psi\rangle = \left( \sum_j \alpha_j \phi_j(\mathbf{r}) \right) | \mathbf{0} \rangle
This shows that the result of the wavefunction operator \hat{\psi}(\mathbf{r}) acting on a superposition state is the superposition of the individual wavefunctions \phi_j(\mathbf{r}) weighted by their respective coefficients \alpha_j, all acting on the vacuum state |\mathbf{0}\rangle.
For a two-particle system, the wavefunction operator \hat{\psi}(\mathbf{r}) can be defined in terms of the creation and annihilation operators for two particles. The field operator for a two-particle state is defined as:
\hat{\psi}(\mathbf{r}_1, \mathbf{r}_2) = \frac{1}{\sqrt 2}\sum_{j,k} \phi_j(\mathbf{r}_1)\phi_k(\mathbf{r}_2) \mathbf{b}_k \mathbf{b}_j
Here, \phi_{j,k}(\mathbf{r}_1, \mathbf{r}_2) represents the two-particle wavefunctions (modes), and \mathbf{b}_j and \mathbf{b}_k are the annihilation operators for the j^{th} and k^{th} modes, respectively. The square root factor in front is to ensure that the solution is normalized.
This operator, operating on a state with two different single-particles states occupied, leads mathematically to a linear combination of product of wavefunctions that is correctly antisymmetric with respect to exchange of these two particle without the need of postulating it. Let’s consider the state:
\mathbf b_l^\dag \mathbf b_m^\dag | \mathbf 0 \rangle
And act on this vector with the wavefunction operator:
\hat{\psi}(\mathbf{r}_1, \mathbf{r}_2)\mathbf b_l^\dag \mathbf b_m^\dag | \mathbf 0 \rangle = \left(\frac{1}{\sqrt 2}\sum_{j,k} \phi_j(\mathbf{r}_1)\phi_k(\mathbf{r}_2) \mathbf{b}_k \mathbf{b}_j\right)\mathbf b_l^\dag \mathbf b_m^\dag | \mathbf 0 \rangle = \frac{1}{\sqrt 2}\sum_{j,k} \phi_j(\mathbf{r}_1)\phi_k(\mathbf{r}_2) \mathbf{b}_k \mathbf{b}_j\mathbf b_l^\dag \mathbf b_m^\dag | \mathbf 0 \rangle
It is possible to use the anticommutation relations to move the annihilator operators to the left:
\begin{aligned} \mathbf{b}_k \mathbf{b}_j\mathbf b_l^\dag \mathbf b_m^\dag | \mathbf 0 \rangle & = \mathbf{b}_k \left( \delta_{jl} - \mathbf b_l^\dag\mathbf{b}_j \right)\mathbf b_m^\dag | \mathbf 0 \rangle \\ & = \left[\delta_{jl} \mathbf{b}_k \mathbf b_m^\dag - \mathbf{b}_k \mathbf b_l^\dag\mathbf{b}_j \mathbf b_m^\dag \right]| \mathbf 0 \rangle \\ & = \left[\delta_{jl} \left(\delta_{km} - \mathbf b_m^\dag\mathbf{b}_k \right) - \left(\delta_{kl}- \mathbf b_l^\dag\mathbf{b}_k \right)\left( \delta_{jm}- \mathbf b_m^\dag\mathbf{b}_j \right) \right]| \mathbf 0 \rangle \\ & = \delta_{jl}\delta_{km} - \delta_{jl}\mathbf b_m^\dag\mathbf{b}_k - \delta_{kl}\delta_{jm} +\delta_{kl}\mathbf b_m^\dag\mathbf{b}_j + \mathbf b_l^\dag\mathbf{b}_k\delta_{jm} + \mathbf b_l^\dag\mathbf{b}_k\mathbf b_m^\dag\mathbf{b}_j | \mathbf 0 \rangle \\ & = \delta_{jl}\delta_{km} - \delta_{kl}\delta_{jm} | \mathbf 0 \rangle \end{aligned}
As the term of like \mathbf b_m^\dag\mathbf{b}_k | \mathbf 0 \rangle are all equal to zero because annihilating on an empty space.
The results it therefore:
\hat{\psi}(\mathbf{r}_1, \mathbf{r}_2)\mathbf b_l^\dag \mathbf b_m^\dag | \mathbf 0 \rangle = \frac{1}{\sqrt 2} \sum_{j,k} \phi_j(\mathbf{r}_1)\phi_k(\mathbf{r}_2) \left(\delta_{jl}\delta_{km} - \delta_{kl}\delta_{jm}\right) | \mathbf 0 \rangle
Applying the Kronecker delta, the only two term that are nonzero are j=l,k=m and j=m,k=l:
\hat{\psi}(\mathbf{r}_1, \mathbf{r}_2)\mathbf b_l^\dag \mathbf b_m^\dag | \mathbf 0 \rangle = \frac{1}{\sqrt 2} \left( \phi_l(\mathbf{r}_1)\phi_m(\mathbf{r}_2) - \phi_l(\mathbf{r}_2)\phi_m(\mathbf{r}_1) \right) | \mathbf 0 \rangle
So the antisymmetry with respect to the particle exchange is naturally emerging.
For N particles, the field operator can be extended as follows:
\hat{\psi}(\mathbf{r}_1, \mathbf{r}_2, \ldots, \mathbf{r}_N) = \frac{1}{\sqrt N} \sum_{1, 2, \ldots, N} \phi_1(\mathbf{r}_1) \phi_2(\mathbf{r}_2) \dots \phi_n(\mathbf{r}_N) \mathbf{b}_{1} \mathbf{b}_2 \dots \mathbf{b}_N
This operators will be correctly dealing with the permutation for particle exchange.
In the occupation number formalism, when dealing with multiple particles, the operator \mathbf{\hat G} corresponding to each particle must be considered. For an operator \mathbf{\hat G} acting on a single-particle state, the extension to multiple particles involves summing the contributions of \mathbf{\hat G} for each state. If \mathbf{\hat G} is an operator that acts on the j^{th} state, the multi-particle version can be expressed as:
\mathbf{\hat G}_{\mathbf r} = \sum_{a=1}^N \mathbf{G}{\mathbf r_a}
The operator \mathbf{ \hat G} is postulated to be:
\mathbf{ \hat G} = \int \hat \psi^\dag \mathbf{G}_{\mathbf r} \hat \psi \mathrm d^3 \mathbf r_1 \dots \mathrm d^3 \mathbf r_N
Substituting the wavefunction operator:
\begin{aligned} \mathbf{\hat G} & = \int \hat{\psi}^\dag(\mathbf{r}) \mathbf{\hat G}_{\mathbf{r}} \hat{\psi}(\mathbf{r}) \, \mathrm{d}^3\mathbf{r}_1 \, \dots \, \mathrm{d}^3\mathbf{r}_N \\ &= \frac{1}{N} \sum_{a=1}^N \int \left( \sum_j \bar \phi_j(\mathbf{r}) \mathbf{b}_j^\dag \right) \mathbf{\hat G}_{\mathbf r_a} \left( \sum_k \phi_k(\mathbf{r}) \mathbf{b}_k \right) \mathrm{d}^3\mathbf{r}_1\dots \mathrm{d}^3\mathbf{r}_N \\ & = \frac{1}{N} \sum_{a=1}^N \sum_{j,k} \mathbf{b}_j^\dag \mathbf{b}_k \int \bar \phi_j(\mathbf{r}) \mathbf{\hat G}_{\mathbf r_a} \phi_k(\mathbf{r}) \mathrm{d}^3\mathbf{r}_1 \dots \mathrm{d}^3\mathbf{r}_N \\ & = \frac{1}{N} \sum_{a=1}^N \sum_{j,k} \mathbf{b}_j^\dag \mathbf{b}_k G_{jk} \end{aligned}
where:
G_{jk} \equiv \int \bar \phi_j(\mathbf{r}) \mathbf{\hat G}_{\mathbf r_a} \phi_k(\mathbf{r}) \mathrm{d}^3\mathbf{r}_1 \dots \mathrm{d}^3\mathbf{r}_N
It was also used the orthogonality of the functions \phi_n to remove terms like \mathbf b_1^\dag \mathbf b_2 \int \bar \phi_1\phi_2.
It can be noted that this integral does not depend from the particle considered, so it will be always the same so summing N times will gives N G_{jk} and finally:
\mathbf{\hat G} = \frac{1}{N} \sum_{a=1}^N \sum_{j,k} \mathbf{b}_j^\dag \mathbf{b}_k G_{jk} = \frac{1}{N} N \sum_{j,k} \mathbf{b}_j^\dag \mathbf{b}_k G_{jk} = \sum_{j,k} \mathbf{b}_j^\dag \mathbf{b}_k G_{jk}
The Hamiltonian is a particular case that worth analysis. Starting from a single fermion, the Hamiltonian is:
\mathbf H_{\mathbf r} = -\frac{\hbar^2}{2m}\nabla^2 + V(\mathbf r)
Using the definition above:
\begin{aligned} \mathbf{ \hat{H}} & = \sum_{j,k} \mathbf{b}_k^\dag \mathbf{b}_j H_{ij} = \sum_{j,k} \mathbf{b}_k^\dag \mathbf{b}_j \int \bar \phi_k(\mathbf{r}) \mathbf H_{\mathbf r} \mathbf \phi_j(\mathbf{r}) \mathrm{d}^3\mathbf{r} \\ & = \sum_{j,k} \mathbf{b}_k^\dag \mathbf{b}_j\int \bar \phi_k(\mathbf{r})\left[-\frac{\hbar^2}{2m}\nabla^2 + V(\mathbf r) \right] \mathbf \phi_j(\mathbf{r}) \mathrm{d}^3\mathbf{r} \\ & = \sum_{j,k} \mathbf{b}_k^\dag \mathbf{b}_j\int \bar \phi_k(\mathbf{r})E_k \mathbf \phi_j(\mathbf{r}) \mathrm{d}^3\mathbf{r} \\ & = \sum_{j,k} \mathbf{b}_k^\dag \mathbf{b}_j E_k \int \bar \phi_k(\mathbf{r}) \mathbf \phi_j(\mathbf{r}) \mathrm{d}^3\mathbf{r} \\ & = \sum_{j,k} \mathbf{b}_k^\dag \mathbf{b}_j E_k \delta_{jk} = \sum_{j} \mathbf{b}_j^\dag \mathbf{b}_j E_j \end{aligned}
So this operator is diagonal and, as demonstrated, it is valid for one or more non-interacting particles without particular need to change it. Lets compute the expectation of the energy in a superposition of states for a single fermion:
\langle E \rangle = \langle \psi | \mathbf H | \psi \rangle
For a superposition of states, the wavefunction can be written as linear combination of basis function with some coefficient \alpha_i:
| \psi \rangle = \sum \alpha_i | \phi_i \rangle = \sum \alpha_i \mathbf b_i^\dag | \mathbf 0 \rangle
Then, using the anticommutation relations to move the annihilator operators to the right and ignoring the annihilator operators acting on the vacuum space:
\begin{aligned} \langle E \rangle & = \langle \psi | \mathbf H | \psi \rangle = \sum_{i,j,k} \langle 0 | \bar \alpha_i \mathbf b_i \mathbf{b}_j^\dag \mathbf{b}_j E_j \alpha_k \mathbf b^\dag_k | 0 \rangle \\ & = \sum_{i,j,k} \bar \alpha_i \alpha_k E_j \langle 0 | \mathbf b_i \mathbf{b}_j^\dag \mathbf{b}_j \mathbf b^\dag_k | 0 \rangle = \sum_{i,j,k} \bar \alpha_i \alpha_k E_j \langle 0 | \left[ \left(\delta_{ij} - \mathbf{b}_j^\dag\mathbf b_i\right)\left(\delta_{jk} - \mathbf b_k^\dag \mathbf{b}_j\right) \right]| 0 \rangle \\ & = \sum_{i,j,k} \bar \alpha_i \alpha_k E_j \langle 0 |\delta_{ij} \delta_{jk} | 0 \rangle = \sum_{i,j,k} \bar \alpha_i \alpha_k E_j \delta_{ij} \delta_{jk} \langle 0 | 0 \rangle \\ & = \sum_i |\alpha_i|^2 E_i \end{aligned}
which is the same result that could have been computed from the first quantization.
The power of this framework is that considering multiple particles, the Hamiltonian does not need to be changed, but we can simply apply it as is. Considering two particles in different states:
| \psi \rangle = \mathbf b_i^\dag \mathbf b_j^\dag |0\rangle
Then the energy expectation is:
\begin{aligned} \langle E \rangle & = \langle \psi | \mathbf H | \psi \rangle = \sum_{i,j,k} \langle 0 | \mathbf b_j \mathbf b_i \mathbf{b}_k^\dag \mathbf{b}_k E_k \mathbf b^\dag_i \mathbf b^\dag_j | 0 \rangle = \sum_{i,j,k} E_k \langle 0 | \mathbf b_j \mathbf b_i \mathbf{b}_k^\dag \mathbf{b}_k \mathbf b^\dag_i \mathbf b^\dag_j | 0 \rangle \\ & = \sum_{i,j,k} E_k \langle 0 | \mathbf b_j \mathbf b_i \mathbf{b}_k^\dag \mathbf{b}_k \mathbf b^\dag_i \mathbf b^\dag_j | 0 \rangle = \sum_{i,j,k} E_k \langle 0 | \left[ \mathbf b_j \left(\mathbf b_i \mathbf{b}_k^\dag \right)\left(\mathbf{b}_k \mathbf b^\dag_i \right) \mathbf b^\dag_j\right] | 0 \rangle \\ & = \sum_{i,j,k} E_k \langle 0 | \left[ \mathbf b_j \left(\delta_{ik} - \mathbf{b}_k^\dag \mathbf b_i \right)\left(\delta_{ik} - \mathbf b^\dag_i \mathbf{b}_k \right) \mathbf b^\dag_j\right] | 0 \rangle \\ & = \sum_{i,j,k} E_k \langle 0 | \left[ \mathbf b_j \left(\delta_{ik} - \mathbf{b}_k^\dag \mathbf b_i \right)\left(\delta_{ik} \mathbf b^\dag_j - \mathbf b^\dag_i\left(\delta_{jk} - \mathbf b^\dag_j \mathbf{b}_k\right) \right) \right] | 0 \rangle \\ & = \sum_{i,j,k} E_k \langle 0 | \left[ \mathbf b_j \left(\delta_{ik} - \mathbf{b}_k^\dag \mathbf b_i \right)\left(\delta_{ik} \mathbf b^\dag_j - \delta_{jk} \mathbf b^\dag_i \right) \right] | 0 \rangle \\ & = \sum_{i,j,k} E_k \langle 0 | \left[ \delta_{ik}\mathbf b_j\mathbf b^\dag_j - \delta_{ik}\delta_{jk} \mathbf b_j\mathbf b^\dag_i - \delta_{ik} \mathbf b_j\mathbf b^\dag_j \mathbf{b}_k \mathbf b^\dag_j + \delta_{jk} \mathbf b_j\mathbf b^\dag_j \mathbf{b}_k \mathbf b^\dag_i \right] | 0 \rangle \\ & = \sum_{i,j,k} E_k \langle 0 | \left[ \delta_{ik}\left(1 -\mathbf b^\dag_j \mathbf b_j\right) + \delta_{jk} \mathbf b_j\mathbf b^\dag_j \mathbf{b}_k \mathbf b^\dag_i \right] | 0 \rangle \\ & = \sum_{i,j,k} E_k \langle 0 | \left[ \delta_{ik}\left(1 -\mathbf b^\dag_j \mathbf b_j\right) + \delta_{jk} \left(1 -\mathbf b^\dag_j \mathbf b_j\right)\left(\delta_{ik}- \mathbf b^\dag_i\mathbf{b}_k \right) \right] | 0 \rangle \\ & = \sum_{i,j,k} E_k \langle 0 | \left(\delta_{ik} +\delta_{jk} \right) | 0 \rangle =\sum_{i,j,k} E_k \left(\delta_{ik} +\delta_{jk} \right) \langle 0 | 0 \rangle \\ & = E_i + E_j \end{aligned}
We used some relations like \delta_{ik}\delta_{jk}=0 as i\ne j or \mathbf b^\dag_j \mathbf b^\dag_j = 0 as it trying to create two fermions in the same state.
This formalism is useful when particles are not interacting with one another, but can be extended to multiple particles all interacting with each other. The operator \mathbf{ \hat G} is again postulated to be:
\mathbf{ \hat G} = \int \hat \psi^\dag \mathbf{G}_{\mathbf r} \hat \psi \mathrm d^3 \mathbf r_1 \dots \mathrm d^3 \mathbf r_N
where the wavefunction operator is now:
\hat{\psi}(\mathbf{r}_1, \mathbf{r}_2, \ldots, \mathbf{r}_N) = \frac{1}{\sqrt N} \sum_{1, 2, \ldots, N} \phi_1(\mathbf{r}_1) \phi_2(\mathbf{r}_2) \dots \phi_n(\mathbf{r}_N) \mathbf{b}_{1} \mathbf{b}_2 \dots \mathbf{b}_N
substituting in the operator expression:
\begin{aligned} \mathbf{ \hat G} & = \frac{1}{N} \sum_{a_1,\dots,a_n,c_1,\dots,c_n} \mathbf b_{c_n}^\dag \dots \mathbf b_{c_1}^\dag \mathbf b_{a_1}\dots \mathbf b_{a_N} \int \bar \phi_{c_1}(\mathbf r_1) \dots \bar \phi_{c_N}(\mathbf r_N) \mathbf{G}_{\mathbf r} (\mathbf r_1,\dots,\mathbf r_n) \phi_{a_1}(\mathbf r_1) \dots \phi_{a_N}(\mathbf r_N) \mathrm d^3 \mathbf r_1 \dots \mathrm d^3 \mathbf r_N \\ & = \frac{1}{N} \sum_{a_1,\dots,a_n,c_1,\dots,c_n} \mathbf b_{c_n}^\dag \dots \mathbf b_{c_1}^\dag \mathbf b_{a_1}\dots \mathbf b_{a_N} G_{a_1,\dots,a_N,c_1,\dots,c_N} \end{aligned}
where:
G_{a_1,\dots,a_N,c_1,\dots,c_N} \equiv \int \bar \phi_{c_1}(\mathbf r_1) \dots \bar \phi_{c_N}(\mathbf r_N) \mathbf{G}_{\mathbf r} (\mathbf r_1,\dots,\mathbf r_n) \phi_{a_1}(\mathbf r_1) \dots \phi_{a_N}(\mathbf r_N) \mathrm d^3 \mathbf r_1 \dots \mathrm d^3 \mathbf r_N
And this operator would not change if there are more than one set of particle non-interacting. In particular, for practical purposes, it is useful to derive fully the two-particles operator, that takes the form:
\mathbf{ \hat G} = \frac{1}{2} \sum_{a_1,a_2,c_1,c_2} \mathbf b_{c_2}^\dag \mathbf b_{c_1}^\dag \mathbf b_{a_1} \mathbf b_{a_2} G_{a_1,a_2,c_1,c_2}
and:
G_{a_1,a_2,c_1,c_2} \equiv \int \bar \phi_{c_1}(\mathbf r_1) \bar \phi_{c_2}(\mathbf r_2) \mathbf{G}_{\mathbf r}(\mathbf r_1, \mathbf r_2) \phi_{a_1}(\mathbf r_1)\phi_{a_2}(\mathbf r_2) \mathrm d^3 \mathbf r_1 \mathrm d^3 \mathbf r_2
Considering two electrons interacting through the Coulomb force. The Hamiltonian is:
\mathbf H \left(\mathbf r_1, \mathbf r_2\right) = -\frac{\hbar}{2m}\left(\nabla^2_{\mathbf r_1} +\nabla^2_{\mathbf r_2}\right) + \frac{e^2}{4\pi\varepsilon_0|\mathbf r_1 - \mathbf r_2|}
Considering two particles in different states:
| \psi \rangle = \mathbf b_i^\dag \mathbf b_j^\dag |0\rangle
Then the energy expectation is:
\begin{aligned} \langle E \rangle & = \langle \psi | \mathbf H | \psi \rangle = \frac{1}{2} \langle 0 |\sum_{a_1,a_2,c_1,c_2} H_{a_1,a_2,c_1,c_2} \mathbf b_j \mathbf b_i \mathbf b_{c_2}^\dag \mathbf b_{c_1}^\dag \mathbf b_{a_1} \mathbf b_{a_2} \mathbf b_i^\dag \mathbf b_j^\dag | 0 \rangle \\ & = \frac{1}{2} \langle 0 |\sum_{a_1,a_2,c_1,c_2} H_{a_1,a_2,c_1,c_2} \mathbf b_j \left(\mathbf b_i \mathbf b_{c_2}^\dag \right)\mathbf b_{c_1}^\dag \mathbf b_{a_1} \left(\mathbf b_{a_2} \mathbf b_i^\dag\right) \mathbf b_j^\dag | 0 \rangle\\ & = \frac{1}{2} \langle 0 |\sum_{a_1,a_2,c_1,c_2} H_{a_1,a_2,c_1,c_2}\mathbf b_j \left(\delta_{ic_2}-\mathbf b_{c_2}^\dag \mathbf b_i \right)\mathbf b_{c_1}^\dag \mathbf b_{a_1} \left(\delta_{ia_2} - \mathbf b_i^\dag \mathbf b_{a_2} \right) \mathbf b_j^\dag | 0 \rangle \\ & = \frac{1}{2} \langle 0 |\sum_{a_1,a_2,c_1,c_2} H_{a_1,a_2,c_1,c_2}\mathbf b_j \left(\delta_{ic_2}-\mathbf b_{c_2}^\dag \mathbf b_i \right)\mathbf b_{c_1}^\dag \mathbf b_{a_1} \left[\delta_{ia_2}\mathbf b_j^\dag - \mathbf b_i^\dag \left(\mathbf b_{a_2}\mathbf b_j^\dag \right)\right] | 0 \rangle \\ & = \frac{1}{2} \langle 0 |\sum_{a_1,a_2,c_1,c_2} H_{a_1,a_2,c_1,c_2} \mathbf b_j \left(\delta_{ic_2}-\mathbf b_{c_2}^\dag \mathbf b_i \right)\mathbf b_{c_1}^\dag \mathbf b_{a_1} \left[\delta_{ia_2}\mathbf b_j^\dag - \mathbf b_i^\dag \left(\delta_{ja_2} - \mathbf b_j^\dag \mathbf b_{a_2}\right)\right] | 0 \rangle \\ & = \frac{1}{2} \langle 0 |\sum_{a_1,a_2,c_1,c_2} H_{a_1,a_2,c_1,c_2} \mathbf b_j \left(\delta_{ic_2}-\mathbf b_{c_2}^\dag \mathbf b_i \right)\mathbf b_{c_1}^\dag \mathbf b_{a_1} \left(\delta_{ia_2}\mathbf b_j^\dag - \delta_{ja_2}\mathbf b_i^\dag \right) | 0 \rangle \\ & = \frac{1}{2} \langle 0 |\sum_{a_1,a_2,c_1,c_2} H_{a_1,a_2,c_1,c_2} \mathbf b_j \left(\delta_{ic_2}-\mathbf b_{c_2}^\dag \mathbf b_i \right)\mathbf b_{c_1}^\dag \left(\delta_{ia_2}\mathbf b_{a_1} \mathbf b_j^\dag - \delta_{ja_2}\mathbf b_{a_1} \mathbf b_i^\dag \right) | 0 \rangle \\ & = \frac{1}{2} \langle 0 |\sum_{a_1,a_2,c_1,c_2} H_{a_1,a_2,c_1,c_2} \mathbf b_j \left(\delta_{ic_2}-\mathbf b_{c_2}^\dag \mathbf b_i \right)\mathbf b_{c_1}^\dag \left[\delta_{ia_2}\left(\delta_{ja_1} - \mathbf b_j^\dag\mathbf b_{a_1} \right) - \delta_{ja_2}\left(\delta_{ia_1} - \mathbf b_i^\dag \mathbf b_{a_1} \right) \right] | 0 \rangle \\ & = \frac{1}{2} \langle 0 |\sum_{a_1,a_2,c_1,c_2} H_{a_1,a_2,c_1,c_2} \mathbf b_j \left(\delta_{ic_2}-\mathbf b_{c_2}^\dag \mathbf b_i \right)\mathbf b_{c_1}^\dag \left(\delta_{ia_2}\delta_{ja_1} - \delta_{ja_2}\delta_{ia_1}\right) | 0 \rangle \\ & = \frac{1}{2} \langle 0 |\sum_{a_1,a_2,c_1,c_2} H_{a_1,a_2,c_1,c_2} \mathbf b_j \left[\delta_{ic_2}\mathbf b_{c_1}^\dag-\mathbf b_{c_2}^\dag \left(\mathbf b_i\mathbf b_{c_1}^\dag\right) \right] \left(\delta_{ia_2}\delta_{ja_1} - \delta_{ja_2}\delta_{ia_1}\right) | 0 \rangle \\ & = \frac{1}{2} \langle 0 |\sum_{a_1,a_2,c_1,c_2} H_{a_1,a_2,c_1,c_2} \mathbf b_j \left[\delta_{ic_2}\mathbf b_{c_1}^\dag-\mathbf b_{c_2}^\dag \left(\delta_{ic_1} -\mathbf b_{c_1}^\dag\mathbf b_i\right) \right] \left(\delta_{ia_2}\delta_{ja_1} - \delta_{ja_2}\delta_{ia_1}\right) | 0 \rangle \\ & = \frac{1}{2} \langle 0 |\sum_{a_1,a_2,c_1,c_2} H_{a_1,a_2,c_1,c_2} \mathbf b_j \left(\delta_{ic_2}\mathbf b_{c_1}^\dag-\delta_{ic_1}\mathbf b_{c_2}^\dag\right) \left(\delta_{ia_2}\delta_{ja_1} - \delta_{ja_2}\delta_{ia_1}\right) | 0 \rangle \\ & = \frac{1}{2} \langle 0 |\sum_{a_1,a_2,c_1,c_2} H_{a_1,a_2,c_1,c_2} \left(\delta_{ic_2}\mathbf b_j \mathbf b_{c_1}^\dag-\delta_{ic_1}\mathbf b_j \mathbf b_{c_2}^\dag\right) \left(\delta_{ia_2}\delta_{ja_1} - \delta_{ja_2}\delta_{ia_1}\right) | 0 \rangle\\ & = \frac{1}{2} \langle 0 |\sum_{a_1,a_2,c_1,c_2} H_{a_1,a_2,c_1,c_2} \left[\delta_{ic_2}\left(\delta_{jc_1}- \mathbf b_{c_1}^\dag\mathbf b_j\right) - \delta_{ic_1}\left(\delta_{jc_2} -\mathbf b_{c_2}^\dag\mathbf b_j \right)\right] \left(\delta_{ia_2}\delta_{ja_1} - \delta_{ja_2}\delta_{ia_1}\right) | 0 \rangle \\ & = \frac{1}{2} \langle 0 |\sum_{a_1,a_2,c_1,c_2} H_{a_1,a_2,c_1,c_2} \left(\delta_{ic_2}\delta_{jc_1} - \delta_{ic_1}\delta_{jc_2} \right) \left(\delta_{ia_2}\delta_{ja_1} - \delta_{ja_2}\delta_{ia_1}\right) | 0 \rangle \\ & = \frac{1}{2} \langle 0 |\sum_{a_1,a_2,c_1,c_2} H_{a_1,a_2,c_1,c_2} \left(\delta_{ic_2}\delta_{jc_1}\delta_{ia_2}\delta_{ja_1} + \delta_{ic_1}\delta_{jc_2}\delta_{ja_2}\delta_{ia_1} - \delta_{ic_2}\delta_{jc_1} \delta_{ja_2}\delta_{ia_1} - \delta_{ic_1}\delta_{jc_2}\delta_{ia_2}\delta_{ja_1}\right) | 0 \rangle \end{aligned}
With this result, only a limited number of terms are non zero, the ones all the indexes are consistent (and the other three terms will be zero because at least one of the components is zero):
\langle E \rangle = \frac{1}{2} \left(H_{jiji} + H_{ijij} + H_{ijji} + H_{jiij} \right)
Explicitly:
\begin{aligned} & H_{jiji} = H_{ijij} = \int \bar \phi_j(\mathbf r_1) \bar \phi_i(\mathbf r_2) \mathbf{H} \phi_j(\mathbf r_1)\phi_i(\mathbf r_2) \mathrm d^3 \mathbf r_1 \mathrm d^3 \mathbf r_2 \\ & H_{ijji} = H_{jiij} = \int \bar \phi_j(\mathbf r_1) \bar \phi_i(\mathbf r_2) \mathbf{H} \phi_i(\mathbf r_1)\phi_j(\mathbf r_2) \mathrm d^3 \mathbf r_1 \mathrm d^3 \mathbf r_2 \end{aligned}
The first term is:
\begin{aligned} H_{jiji} = & \int \bar \phi_j(\mathbf r_1) \bar \phi_i(\mathbf r_2) \mathbf{H} \phi_j(\mathbf r_1)\phi_i(\mathbf r_2) \mathrm d^3 \mathbf r_1 \mathrm d^3 \mathbf r_2 \\ = & \int \bar \phi_j(\mathbf r_1) \bar \phi_i(\mathbf r_2) \left(-\frac{\hbar}{2m}\left(\nabla^2_{\mathbf r_1} +\nabla^2_{\mathbf r_2}\right) + \frac{e^2}{4\pi\varepsilon_0|\mathbf r_1 - \mathbf r_2|}\right) \phi_j(\mathbf r_1)\phi_i(\mathbf r_2) \mathrm d^3 \mathbf r_1 \mathrm d^3 \mathbf r_2 \\ = &- \int \bar \phi_j(\mathbf r_1) \bar \phi_i(\mathbf r_2) \frac{\hbar}{2m}\nabla^2_{\mathbf r_1} \phi_j(\mathbf r_1)\phi_i(\mathbf r_2) \mathrm d^3 \mathbf r_1 \mathrm d^3 \mathbf r_2 \\ & -\int \bar \phi_j(\mathbf r_1) \bar \phi_i(\mathbf r_2) i\frac{\hbar}{2m}\nabla^2_{\mathbf r_2} \phi_j(\mathbf r_1)\phi_i(\mathbf r_2) \mathrm d^3 \mathbf r_1 \mathrm d^3 \mathbf r_2 \\ & + \int \bar \phi_j(\mathbf r_1) \bar \phi_i(\mathbf r_2) \frac{e^2}{4\pi\varepsilon_0|\mathbf r_1 - \mathbf r_2|} \phi_j(\mathbf r_1)\phi_i(\mathbf r_2) \mathrm d^3 \mathbf r_1 \mathrm d^3 \mathbf r_2 \\ = & - \int \bar \phi_j(\mathbf r_1) \frac{\hbar}{2m}\nabla^2_{\mathbf r_1} \phi_j(\mathbf r_1) \mathrm d^3 \mathbf r_1 \left(\bar \phi_i(\mathbf r_2) \phi_i(\mathbf r_2) \mathrm d^3 \mathbf r_2 \right) \\ & -\int \bar \phi_i(\mathbf r_2) i\frac{\hbar}{2m}\nabla^2_{\mathbf r_2} \phi_i(\mathbf r_2) \mathrm d^3 \mathbf r_2 \left(\bar \phi_j(\mathbf r_1)\phi_j(\mathbf r_1) \mathrm d^3 \mathbf r_1 \right) \\ & \int \bar \phi_j(\mathbf r_1) \bar \phi_i(\mathbf r_2) \frac{e^2}{4\pi\varepsilon_0|\mathbf r_1 - \mathbf r_2|} \phi_j(\mathbf r_1)\phi_i(\mathbf r_2) \mathrm d^3 \mathbf r_1 \mathrm d^3 \mathbf r_2 \\ = & - \int \bar \phi_j(\mathbf r_1) \frac{\hbar}{2m}\nabla^2_{\mathbf r_1} \phi_j(\mathbf r_1) \mathrm d^3 \mathbf r_1 -\int \bar \phi_i(\mathbf r_2) i\frac{\hbar}{2m}\nabla^2_{\mathbf r_2} \phi_i(\mathbf r_2) \mathrm d^3 \mathbf r_2 \\ & + \int \bar \phi_j(\mathbf r_1) \bar \phi_i(\mathbf r_2) \frac{e^2}{4\pi\varepsilon_0|\mathbf r_1 - \mathbf r_2|} \phi_j(\mathbf r_1)\phi_i(\mathbf r_2) \mathrm d^3 \mathbf r_1 \mathrm d^3 \mathbf r_2 \\ = & E_{K_1} + E_{K_2} + E_{P_{12}} \end{aligned}
These terms are the standard energies, the kinetic energy of the two electron and the potential energy from the Coulomb interaction.
The remaining is an additional energy:
H_{ijji} = \int \bar \phi_j(\mathbf r_1) \bar \phi_i(\mathbf r_2) \mathbf{H} \phi_i(\mathbf r_1)\phi_j(\mathbf r_2) \mathrm d^3 \mathbf r_1 \mathrm d^3 \mathbf r_2
This is called exchange energy and arises due to the requirement that the total wavefunction of a system of fermions must be antisymmetric under particle exchange, as dictated by the Pauli exclusion principle.
This energy term represents the interaction between identical particles and is purely a quantum mechanical effect with no classical analog. It emerges from the overlap of the wavefunctions of the two electrons when their positions are exchanged. The exchange energy contributes to the total energy of the system by accounting for the correlation between the fermions due to their indistinguishability and the antisymmetric nature of their combined wavefunction.