Non-interacting fermions

Overview

In many-body physics, we often look for good approximations of complicated interacting models so that we can obtain analytical tools for calculating relevant physical properties. In the context of fermionic many-body systems, non-interacting fermionic models are often used as a basis to construct such approximations.

A prelude on fermions

Fermionic particles

Fermions are a class of quantum particles defined by their obedience to the Pauli exclusion principle.

note

The Pauli exclusion principle states that identical fermionic particles can never occupy the same physical quantum state. This contrasts with bosons, which can accumulate in the same state without restriction.

An example of a fermion is the electron which has total spin $1/2$.

note

Any fermionic particle must have a half-integer total spin. The quantum spin of a particle is an intrinsic property, just as electric charge is an intrinsic property.

In first quantization, to describe a system of fermions one first defines a relevant basis of functions $\phi_s(\vec{r})$ where $s$ labels the different states the fermion can occupy. The wave function corresponding to having a fermion in each state $\phi_{s_1}, \phi_{s_2} ..., \phi_{s_L}$ can be written as

$$\psi(\vec{r}_1,...,\vec{r}_L) = F(\phi_{s_1},...,\phi_{s_L})$$

where $L$ is the number of fermions, $\vec{r}_i$ is the position coordinate of the $i^{th}$ particle and $F$ is a function whose structure is discussed below.

Identical fundamental particles are indistinguishable from one another. As a consequence the coordinates $\vec{r}_i$ cannot be tied to any specific fermionic particle. The function $F$ must be chosen such that exchanging a pair of position labels $\vec{r}_j$ keeps the wave function invariant up to an unobservable global phase.

Taking $P_{i,j}$ to be the operator which exchanges the position labels $\vec{r}_i$ and $\vec{r}_j$, it then holds that for fermions

$$P_{i,j} \psi(...,\vec{r}_i, ...,\vec{r}_j, ...) = - \psi(...,\vec{r}_j, ...,\vec{r}_i, ...)$$

where we see that the acquired phase is $-1$.

note

The reason for the negative phase is tied to the Pauli exclusion principle. Indeed, if fermion $i$ and $j$ were in the same state, then exchanging the fermions keeps the wave function invariant which then implies $\psi = -\psi$ and thus $\psi = 0$, so no fermions can occupy the same state as expected.

The only possible form of $F(\phi_{s_1},...,\phi_{s_L})$ which satisfies the above constraints is given by

$$F(\phi_{s_1},...,\phi_{s_L}) = \frac{1}{\sqrt{L!}} \sum_{P}\text{sgn}(P)\phi_{s_{P(1)}}(\vec{r}_1)...\phi_{s_{P(L)}}(\vec{r}_L)$$

where $P$ labels the possible permutations of the state indices and $\text{sgn}(P)$ is the sign of the permutation ($1$ if it can be constructed from an even number of pairwise permutations, $-1$ otherwise). More complicated fermionic states are obtained by taking linear superpositions of the function $F$ evaluated with different states $\phi_{s_j}$ and potentially different number of particles if the state describes a situation in which the number of fermions is not fixed.

The Fock space

Albeit possible to work exclusively with wave functions, a more convenient way to describe systems with multiple fermions is instead to use a framework that captures the possible states of the fermions together with the occupation of each available state. This can be achieved by defining the so-called Fock space. In this representation, the available states are labeled $s_1, s_2, ..., s_N$ where $N$ is the total number of states. A many-body fermionic configuration can then be represented as $\vert n_1n_2...n_N \rangle$ where $n_k \in \{1,0\}$ (correspondingly occupied or unoccupied) is the occupation number of the corresponding state $s_k$. The states $\vert n_1n_2...n_N \rangle$ are in one to one correspondence with the states $F(\phi_{s_1},...,\phi_{s_L})$ discussed above. As an example, suppose $N = 4$ and $(n_1 = 0, n_2 = 1, n_3 = 0, n_4 = 1)$, then we have that

$$F(\phi_{s_2},\phi_{s_4}) \equiv \ket{0,1,0,1}$$

Fermionic operators

Now that we have a way of describing the available states the fermions can be in, the natural next step is to introduce operators that can change the occupancy of the available states. Since the only relevant information to describe a system of fermions is the occupancy of the available states, we introduce corresponding creation and annihilation operators $f^{\dagger}_{k}$ and $f_{k}$ which respectively add or remove a fermion from the state $s_k$. Sums of products of these operators alone is sufficient to describe the dynamics of any system of fermions since any physical process can be viewed as an update rule for the population of the available fermionic states. Any Fock state can be represented in terms of these operators as

$$\ket{n_1, n_2, ...} = ((f_1^{\dagger})^{n_1}(f_2^{\dagger})^{n_2}...)\ket{0}$$

where $\ket{0}$ is the vacuum state (the state for which $n_j = 0$ $\forall j$). The action of the creation and annihilation operators on any Fock states is given by

$$f_j^{\dagger}\vert ...n_j... \rangle = (-1)^{\sum_{r < j}n_r}(1 - n_j)\vert ...(1 + n_j)... \rangle \quad f_j\vert ...n_j... \rangle = (-1)^{\sum_{r < j}n_r}n_j\vert ...(n_j-1)... \rangle.$$

The creation and annihilation operators correspond to the physical process of adding a fermion in the specified state $s_j$ to the many-body wave function in such a way that the resulting state continues to respect the Pauli exclusion principle and does not break the indistinguishably property. This in turn enforces the following anti-commutation relations

$$f_if_j + f_jf_i = 0$$ $$f_i^{\dagger}f_j^{\dagger} + f_j^{\dagger}f_i^{\dagger} = 0$$ $$f_i^{\dagger}f_j + f_jf_i^{\dagger} = \delta_{i \neq j}$$

Quadratic fermionic models

A system of fermions is said to be quadratic whenever the Hamiltonian can be written down as a sum of pairs of creation and annihilation operators.

$$H_{\text{quadratic}} = \sum_{j,k} A_{j,k} f_j^{\dagger}f_k + \frac{1}{2}(B_{j,k}f_j^{\dagger}f_k^{\dagger} + \text{h.c}).$$

note

An interacting Hamiltonian would include terms such as $f_j^\dagger f_k^\dagger f_{j + \delta} f_{k + \delta} + \text{h.c.}$ which here represents an interaction between two fermionic particles.

Quadratic Hamiltonians can always be diagonalized with a Bogoliubov–de Gennes transformation. Aleph provides several operators that can be used to represent either the totality or part of $H_{\text{quadratic}}$. Before discussing the available Aleph functionalities it is useful to consider the more general quadratic operator

$$O_{\text{quadratic}} = \sum_{j,k} M^{+,-}_{j,k} f_j^{\dagger}f_k + M^{-,+}_{j,k} f_jf_k^{\dagger} + M^{+,+}_{j,k}f_j^{\dagger}f_k^{\dagger} + M^{-,-}_{j,k}f_jf_k + \alpha I$$

where $\alpha$ is some constant.

note

We do not enforce $O_{\text{quadratic}}$ to be Hermitian. As such, the matrices $M^{+,+},M^{+,-},M^{-,+},M^{-,-}$ can be chosen arbitrarily.

Any operator of this form can be compactly represented by making use of the so-called Nambu spinor. Given

$$\Psi^- = \begin{pmatrix} f_1 \\ \vdots \\ f_L \end{pmatrix}, \quad \Psi^+ = \begin{pmatrix} f_1^\dagger \\ \vdots \\ f_L^\dagger \end{pmatrix}$$

the Nambu spinor can be written as

$$\Psi = \begin{pmatrix} \Psi^- \\ \Psi^+ \end{pmatrix}.$$

The operator $O_{\text{quadratic}}$ can be written compactly as

$$O_{\text{quadratic}} = \Psi^\dagger M \Psi + \alpha I$$

where we have defined

$$M = \begin{pmatrix} M^{+,-}, M^{+,+} \\ M^{-,-}, M^{-,+} \end{pmatrix}.$$

The matrix $M$ is not unique due to the fermionic anti-commutation relationships $\{f_i,f_j\} = 0$, $\{f_i^{\dagger},f_j^{\dagger}\} = 0$, $\{f_i^{\dagger},f_j\} = \delta_{i,j}$

The representation of $O_{\text{quadratic}}$ can be fixed by instead considering the matrix

$$Q = \begin{pmatrix} Q^{+,-}, Q^{+,+} \\ Q^{-,-}, Q^{-,+} \end{pmatrix}$$

where

$$\begin{aligned} Q^{+,+} &= (M^{+,+} - (M^{+,+})^T)/2 \\ Q^{-,-} &= (M^{-,-} - (M^{-,-})^T)/2 \\ Q^{+,-} &= (M^{+,-} - (M^{-,+})^T)/2 \\ Q^{-,+} &= (M^{-,+} - (M^{+,-})^T)/2 \\ E_{\text{shift}} &= \text{Tr}(M)/2 + \alpha \end{aligned}$$

which allows us to write the following compact representation of any quadratic fermionic operator

$$O_{\text{quadratic}} = \Psi^\dagger Q \Psi + E_{\text{shift}} I.$$

We note that in this form it holds true that

$$Q^{-,+} = -\left(Q^{+,-}\right)^T,$$ $$Q^{+,+} = -\left(Q^{+,+}\right)^T,$$ $$Q^{-,-} = -\left(Q^{-,-}\right)^T.$$

When the operator $O_{\text{quadratic}}$ is Hermitian the matrix $\frac{1}{2}Q$ corresponds to the well known Bogoliubov-de Gennes representation $H_{\text{BDG}}$ of the Hamiltonian. Given the constraints imposed by the structure of the matrix $Q$ we see that out of the $2L \times 2L$ parameters of the matrix $Q$ only $2 \times L \times L - L$ of them are non-redundant. Indeed $Q^{+,-}$ can be used to deduce $Q^{-,+}$ and the upper/lower triangular part of $Q^{+,+}/Q^{-,-}$ can be used to find the lower/upper triangular part of $Q^{+,+}/Q^{-,-}$. Furthermore, the diagonal part of $Q^{+,+}$ and $Q^{-,-}$ are fixed to zero. These constraints motivate the introduction of two operators. The first one TransferSum is tied to the $Q^{+,-}$ block. Explicitly we have that

$\text{TransferSum}(A) = \Psi^+ A \Psi^- = \sum_{j,k}A_{j,k}f_j^{\dagger}f_k$

for some arbitrary matrix $A$. This operator is compatible with the single-particle functionalities of Aleph (see the next section for more details), as it conserves the total particle number. The second relevant operator is the PairingSum operator. Explicitly we have that

$\text{PairingSum}(P) = \Psi^+ P_{<} \Psi^+ + \Psi^- P_{>} \Psi^- = \sum_{j < k}P_{j,k}f_j^\dagger f_k^\dagger + \sum_{j > k}P_{j,k}f_j f_k$ where we have used the notation $P_{<}$ and $P_{>}$ to denote the upper and lower triangular part of the matrix $P$ respectively. A generic $O_{\text{quadratic}}$ can thus be written as

$O_{\text{quadratic}} =\text{TransferSum}(A) + \text{PairingSum}(P) + \beta I$

which implicitly defines the relationships

$$\begin{aligned} &Q^{+,-} = \frac{1}{2}A \\ &Q^{-,+} = -\frac{1}{2}A^T \\ &Q^{+,+} = \frac{1}{2}(P_{<} - P_{<}^T) \\ &Q^{-,-} = \frac{1}{2}(P_{>} - P_{>}^T) \\ &E_{\text{shift}} = \beta + \mathrm{Tr}(A)/2 \end{aligned}$$

For convenience the operators TransferSum and PairingSum can be constructed either by providing their explicit parameter $A$ and $P$ or simply by passing a sum of compatible fermionic operators.

Signature	Description
TransferSum(matrix)	Pass the A matrix explicitly.
TransferSum(operator [, options])	Pass a valid sum of quadratic fermionic terms and optionally a scalar option.
PairingSum(matrix)	Pass the P matrix explicitly.
PairingSum(operator [, options])	Pass a valid sum of quadratic fermionic terms and optionally a scalar option.

note

Both the TransferSum and the PairingSum operator do not admit the identity as part of their valid set of operators. To construct an operator $O_{\text{quadratic}}$ with an explicitly added energy shift simply add a multiple of the Identity() to TransferSum(A) + PairingSum(P).

For convenience the operators TransferSum and PairingSum provide the methods extract_bdg_matrix and extract_bdg_shift which return the corresponding matrix $\frac{1}{2}Q$ and the energy shift $E_{\text{shift}}$ of the operator. Finally, aleph also provides the operator QuadraticSum which is a convenient way of representing $O_{\text{quadratic}}$ as a whole. First we note that it is always possible to write

$$O_{\text{quadratic}} = \Psi^{\dagger} S \Psi = \sum_{j,k}S_{j,k}^{+,-}f_j^{\dagger}f_k + S_{j,k}^{-,+}f_jf_k^{\dagger} + S_{j,k}^{+,+}f_j^{\dagger}f_k^{\dagger} + S_{j,k}^{-,-}f_jf_k$$

by incorporating the identity term in the matrix elements of $S$ by making use of the identity $I = c^\dagger c + c c^\dagger$. In this form one obtains the relationship

$$\begin{aligned} Q^{+,+} &= (S^{+,+} - (S^{+,+})^T)/2 \\ Q^{-,-} &= (S^{-,-} - (S^{-,-})^T)/2 \\ Q^{+,-} &= (S^{+,-} - (S^{-,+})^T)/2 \\ Q^{-,+} &= (S^{-,+} - (S^{+,-})^T)/2 \\ E_{\text{shift}} &= \text{Tr}(S)/2 \end{aligned}$$

Just like TransferSum and PairingSum, the QuadraticSum operator provides the methods extract_bdg_matrix and extract_bdg_shift. A QuadraticSum operator can be constructed either by providing the explicit matrix $Q$ or by passing a sum of valid quadratic operators.

Signature	Description
QuadraticSum(matrix)	Pass the Q matrix explicitly.
QuadraticSum(operator [, options])	Pass a valid sum of quadratic fermionic terms and optionally a scalar option.

As an example a typical workflow to solve problems involving a quadratic Hamiltonian may look like

// Define the Hamiltonian parameters
var t = -1.0 // hopping
var mu = 0.5 // chemical potential
var delta = 0.4 // pairing terms
var alpha = 0.7 // explicit energy shift

// Group terms of interest
var transfer_terms = t*(Hop(0,1) +  Hop(1,2)) + mu*(Number(0) + Number(1)) // transfer terms
var pairing_terms = delta*(CreateCreate(0,1) + CreateCreate(1,2) + DestroyDestroy(1,0) + DestroyDestroy(2,1)) // pairing terms
var identity_term = alpha*Identity()
var quadratic_terms = transfer_terms + pairing_terms + identity_term // quadratic terms

// Build the relevant operators
var transfer_sum = TransferSum(transfer_terms)
var pairing_sum = PairingSum(pairing_terms)
var quadratic_sum = QuadraticSum(quadratic_terms)

// Obtain quantities of interest
var bdg_matrix_hopping = transfer_sum.get_bdg_matrix()
var bdg_shift_hopping = transfer_sum.get_bdg_shift()

var bdg_matrix_pairing = pairing_sum.get_bdg_matrix()
var bdg_shift_pairing = pairing_sum.get_bdg_shift() // Always zero

var bdg_matrix_quadratic = quadratic_sum.get_bdg_matrix()
var bdg_shift_quadratic = quadratic_sum.get_bdg_shift()

In the code snippet above it holds true that bdg_matrix_hopping + bdg_matrix_pairing == bdg_matrix_quadratic and that bdg_shift_quadratic == bdg_shift_pairing + bdg_shift_hopping + alpha. The Bogoliubov-de Gennes matrix and the corresponding energy shift can then be used to find the quasiparticle energies $\epsilon_k$ and quasiparticle eigenmodes $\vec{v}_k$.

The single-particle sector

Often, a relevant first step in studying a many-body fermionic quantum system is to solve the single-particle problem by using some heuristics to replace interacting terms by an effective potential in the original Hamiltonian. A generic single-particle fermionic Hamiltonian can be expressed as

$$H_{\text{free}} = \left(\sum_{j,k} A_{j,k} f_j^{\dagger}f_k\right)$$

note

In Aleph and everywhere in the documentation, we refer to operators that only involve hopping terms as free, in analogy to a particle which is free to roam around in a potential landscape.

The matrix $A_{j,k}$ is a Hermitian matrix when the above representation represents a Hamiltonian. By relaxing the hermitian constraint, the above equation can be used to represent any operator whose action on states is non-trivial only in the single-particle sector.

note

The single-particle sector is characterized by the property $\sum_{i}n_j = 1$. Terms such as $f_j^{\dagger}f_k^{\dagger}$ violate particle conservation and so do not appear in the single-particle Hamiltonian. Furthermore, any fermionic operator involving a product of terms acting on more than one fermion vanishes in the single-particle sector.

From the form of $H_{\text{free}}$, it is clear that single-particle Hamiltonians are a subset of the more general quadratic fermionic models. The additional structure provided by the single-particle constraint can be leveraged to represent both the states and the Hamiltonian efficiently. First, we note that any fermionic state containing only one fermion can be written down as

$$\psi_{\text{free}} = \sum_{j = 0}^{L-1} \alpha_j f_j^{\dagger} \ket{0}$$

where $\sum_j |\alpha_j|^2 = 1$. As such, one may represent a single-particle state with a vector of coefficients $\alpha_0, ..., \alpha_{L-1}$. In Aleph, such a representation can be created as follows

// 12 possible modes for the fermionic particle
var L = 12;

// specify the state vector options to produce a single-particle state vector
var options = as_fermion | as_singleparticle | as_real

//Set a randomized normalized state
var sp_state = state_vector(L, options)
sp_state.set_random()
sp_state.normalize()

//Extract the probability to find the state in a given mode
var mode = 3
var prob_occupied = abs(sp_state[mode])**2 // Yields the probability to find the state in the third fermionic mode

The resulting state vector stores $L$ coefficients whose amplitudes squared correspond to the probability to find the fermionic particle in the specified mode.

Single-particle state vectors are most useful when used with operators that represent a set of single-particle operations. As discussed above, any single-particle operator can be described with a matrix $A_{j,k}$ of size $L$ by $L$. In Aleph, a single-particle operator can be created either by directly specifying the matrix $A_{j,k}$ or by constructing the single-particle operator from a sum of operators which individually act non-trivially on the single-particle sector.

// Similar set-up as above
var L = 4;
var options = as_fermion | as_singleparticle | as_real
var sp_state = state_vector(L, options)
sp_state.set_random()
sp_state.normalize()

// Create a single-particle operator from a matrix
{
var sp_representation = operator_free(matrix([[1,1,0,0],[1,1,0,0],[0,0,1,0],[0,0,0,1]]), options) // Represents Number(0) + Number(1) + Number(3) + Hop(0,1) 
// Basic operations on state vectors
sp_representation >> sp_state
sp_representation * sp_state
}

// Create a single-particle operator from an operator sum
{
var op_sum = Number(0) + Number(1) + Hop(1,2) + Number(3) + Hop(0,1)
var sp_representation = operator_free(op_sum, options) 
// Basic operations on state vectors
sp_representation >> sp_state
sp_representation * sp_state
}

Further Reading

• For a more in depth example of what can be done with single-particle operators, stay tuned for the upcoming tutorial on how to use the functionalities of operator_free.

Documentation Contributors

Pierre-Gabriel Rozon