Free fermions

In this tutorial we use the quantum toolkit, specifically the operator framework and the free fermion features, to study the dynamics and eigenstates of a free fermion system. For a comprehensive introduction to free fermion systems and their representation, see the concept Non-interacting fermions.

note

In the context of non-interacting fermionic models, we make use of the word free as a shorthand for the single-particle sector, in analogy to a single fermionic particle free to roam in a potential landscape.

Below we will see how the free fermion features can be used to efficiently simulate and study the time evolution of a single fermion bouncing in a closed boxed containing a centered potential well. Furthermore, we will extract the eigenstates of the system and study their structure by analyzing the probability distribution of the particle. More precisely, we will:

• Study the dynamics of a free fermion initially confined to one edge of the box.
• Study the probability distribution of the ground-state and show that it is localized near the impurity.

Constructing the Hamiltonian

We want to describe the dynamics of an electron that is free to move in a box with closed boundary conditions in the presence of an impurity. This can be done by first constructing a square lattice of size $L \times L$. In Aleph (see the Lattice tutorial for more details on how to use the lattice functionalities) we can construct a square lattice as follows

var L = 50; // Specify the system size
var bcs = [boundary_condition.open, boundary_condition.open]
var square_lattice = lattice("square", [L,L], bcs); // Create the lattice

Next, we add nearest neighbor hopping terms at every site along both the $x$ and $y$ axis. The hopping terms capture the tendency of the electronic wave-function to delocalize in free space. The hopping Hamiltonian is given by

$$H_{\text{hopping}} = -\gamma\sum_{x=0,y=0}^{L-2}(c_{x,y}^{\dagger}c_{x+1,y} + c_{x,y}^{\dagger}c_{x,y+1}) + \text{h.c}$$

where $\gamma$ sets the energy scale of the hopping terms, chosen here to be $\gamma = 1$. In Aleph, we can construct the hopping hamiltonian as follows

var H_hopping = operator_sum(as_real) // Create an empty initial operator sum

// Create displacement vectors along the x and y direction
var delta_x = lattice_coordinate([1,0]);
var delta_y = lattice_coordinate([0,1]);
var gamma = 1.0;

// Go through every coordinate in the lattice
for (coord : lattice_range(square_lattice)) { 

    // Adding the hopping terms
    var x_coords = [coord, coord + delta_x]; // The coordinates [(x,y), (x+1,y)]
    var y_coords = [coord, coord + delta_y]; // The coordinates [(x,y), (x,y+1)]

    if (square_lattice.contains(coord + delta_x)) // Check if (x + 1, y) is in the lattice
    {
        // Convert the x hopping coordinates to the corresponding matrix indices
        var x_indices = [square_lattice.to_index(x_coords[0]), square_lattice.to_index(x_coords[1])]

        H_hopping += -1.0*gamma*Hop(x_indices);
    }

    if (square_lattice.contains(coord + delta_y)) // Check if (x, y + 1) is in the lattice
    {
        // Convert the y hopping coordinates to the corresponding matrix indices
        var y_indices = [square_lattice.to_index(y_coords[0]), square_lattice.to_index(y_coords[1])]

        H_hopping += -1.0*gamma*Hop(y_indices);
    }
}

Additionally, we add a an impurity at the center of the lattice in the form of a Gaussian potential well

$$H_{\text{impurity}} = U_{\text{center}}\sum_{x=0,y=0}^{L-1} e^{-\frac{1}{2\sigma^2}(( x - \mu_x)^2 + (y - \mu_y)^2)}(n_{x,y})$$

where $\sigma = L/4$ sets the the width of the potential well, $U_{\text{center}} = -6\gamma = -6$ is its maximum depth and $(\mu_x, \mu_y)$ is the position around which it is centered, chosen here to be $(\mu_x, \mu_y) = (L/2,L/2)$.

note

The parameter $\sigma$ is best chosen when used to approximate the surface area that contains the potential well. For a lattice of size $L \times L$ for which the lattice spacing is set to $1$, then the potential well approximately occupies a fraction $\sigma^2/L^2$ of the total surface area. Hence the choice $\sigma = L/4$ implies the potential well occupies about $1/16$ of the total system.

In Aleph, we can construct the impurity Hamiltonian as follows:

var H_impurity = operator_sum(as_real) // Create an empty initial operator sum
var U_center = -6.0
var sigma = L/4.0
var mu_x = (L-1)/2.0  // (L-1) to account for zero indexing
var mu_y = (L-1)/2.0

// Go through every coordinate in the lattice
for (coord : lattice_range(square_lattice)) { 

    var x_site = coord.primitive_vector()[0];
    var y_site = coord.primitive_vector()[1];
    var local_amplitude = U_center*exp(-(1.0/(2.0*sigma**2))*((x_site - mu_x)**2 + (y_site - mu_y)**2));
        
    // Get the corresponding matrix index 
    var n_index = square_lattice.to_index(coord)
    H_impurity += Number(n_index)*local_amplitude;
}

// Construct the total Hamiltonian
var H_tot = H_hopping + H_impurity

where the last line constructs the Hamiltonian of the system

$$H_{\text{tot}} = H_{\text{hopping}} + H_{\text{impurity}}$$

In its current form, the obtained Hamiltonian has not been restricted to the single-particle sector, and so by default will expect to interact with state vectors of size $2^L$. Given our goal we would like to cast the obtained Hamiltonian to its free representation. This can be achieved either by making use of the matrix factory with specified options or by making use of the operator_free factory. The former directly returns the matrix representation of the Hamiltonian in the single-particle sector whilst the second returns an operator representing the free matrix representation of the operator (see the concept Non-interacting fermions for more details). Both these representations can be constructed as follows in Aleph

//First specify the options needed to build a free operator
var free_option = as_fermion | as_singleparticle

//Extract the bare free matrix representation
var free_matrix = matrix(H_tot, free_option)

//Extract the operator representation
var free_operator = operator_free(H_tot, free_option)

note

The choice of constructing either a matrix or an operator depends on the problem at hand. The matrix factory returns a matrix which gives the user a lot of control over the data at the cost of having to ensure that custom made interactions with other objects behave as intended. Making use of `operator_free produces an operator which comes with specialized methods, safety checks and optimizations at the cost of loosing precise control over how it interacts with other objects.

Studying the dynamics

Now that we have constructed the Hamiltonian, our next task is to study the time evolution of a state initially localized to one edge of the system. Here we may do so by concentrating the weight of the electronic probability distribution uniformly on the top edge of the box.

$$\vert \psi(t = 0) \rangle = \frac{1}{\sqrt{L}}\sum_{x = 0}^{L-1}c_{x,0}^{\dagger}\ket{0}$$

In Aleph, we may achieve this as follows

//Create a single particle state
var free_state = state_vector(L**2, as_fermion | as_singleparticle); // The fermionic particle can occupy L**2 possible sites

var initial_coord = lattice_coordinate([0,0])
for (var x_pos = 0; x_pos < L; ++x_pos) //Loop over all the lattice coordinates
{ 
    var coord = lattice_coordinate([x_pos,0])
    var vector_index = square_lattice.to_index(coord);
    free_state[vector_index] = 1.0;
}
//Normalize the state
free_state.normalize()

Next, we want to time-evolve the electronic wave function in time. From the solution to the time-independent Schrödinger equation, we have that

$$\ket{\psi(t = T)} = e^{-iH_{\text{tot}}T}\ket{\psi(0)} = (e^{-iH_{\text{tot}}(T/N_{T})})^{N_{T}}\ket{\psi(0)} = U_{dT}^{N_{T}}\ket{\psi(0)}$$

where $U_{dT} = e^{-iH_{\text{tot}}dT}$, $dT = T/N_{T}$ and $T$ is the total time for which we want to evolve the system. We have also introduced $N_{T}$ which is the total number of time steps we would like to take between $t = 0$ and $t = T$ to study how the wave function evolves in time. As can be observed from the equation above, this can be achieved by applying the unitary $U_{dT}$ repetitively on the initial state $\ket{\psi(0)}$. To do so, we first construct the unitary $U_{dT}$ and then apply it recursively on the state. Below, we show two ways in which this can be achieved in Aleph. The first method makes use of a direct calculation of the matrix exponential. The second method makes use of operator_function which leverages efficient Krylov methods to accelerate the calculation of matrix on vector multiplications. Crucially it avoids the explicit calculation of the matrix exponential. (see the Krylov evolution and operator matrix tutorial for more details on how to use Krylov methods)

Method 1

//Setup
var T = 20.0
var N_T = 200
var dT = T/N_T

// Construct the unitary matrix
var U_dT = operator_free(exp(-1i*dT*free_matrix), as_fermion | as_singleparticle)
var krylov_method = false
for (var time = 0; time < N_T; ++time){
    U_dT >> free_state
}

Method 2

//Setup
var T = 20.0
var N_T = 200
var dT = T/N_T

// Krylov setup
var func = fun[dT](x) { return exp(-1i * dT * x) } // Specify the function that generate the unitary matrix
var Nkrylov = 4

// Construct an operator function
var U_dT = operator_function(free_matrix, func, ["krylov_dimension" : Nkrylov])
var free_coefficients = free_state.coefficients()
for (var time = 0; time < N_T; ++time){
    U_dT >> free_coefficients
}

note

Performing time evolution by making use of the first method will in general be more precise, but it can be orders of magnitude slower than the second method. The appropriate choice is determined by the error tolerance for the problem at hand which is left to the discretion of the user.

The following animation displays the obtained probability distribution of finding the the electron on any site of the lattice as a function of time. We see that the attractive impurity initially focuses the incoming probability distribution near the center. Then, the particle scatters of the potential well. Some of the probability distribution is transmitted whilst some of it is reflected. Interestingly the reflected part of the probability distribution is concentrated around the $x = L/2$ axis, showing that an impurity can act as a sort of quantum lense.

Time evolution

Eigenmode analysis

As a final task in this tutorial, we study the ground state of the single-particle Hamiltonian $H_{\text{tot}}$. The ground state of $H_{\text{tot}}$ is the eigenstate of $H_{\text{tot}}$ with the smallest possible eigenvalue. It can be found by making use of the iterative_eigensolver as follows in Aleph (see the ground-state tutorial for more details on how to use Krylov methods)

var num_eigenpairs = 5
var eigen_pairs = iterative_eigensolver(free_matrix, num_eigenpairs)
var ground_state = eigen_pairs[1][0:L**2,0:1]
var excited_state = eigen_pairs[1][0:L**2,4:5]

The following figures display the obtained probability distribution of finding the the electron on any site of the lattice for the ground state and the fifth excited stationary state. We see that the ground-state is localized around the the center as expected. The fifth excited stationary state exhibits four symmetrically opposite regions of high probability. These are reminiscent of the well known electronic orbitals of the Hydrogen atom. We also observe regions for which the probability distribution vanishes. This fits the well known fact that high energy solutions of the Schrödinger equation have an increasingly larger number of nodes as the energy increases.

Ground state	Fifth excited state

What's next?

• Can this be used to study the energy levels of the hydrogen atom? What would need to be changed?
• Which time evolution method is faster? Which is more precise? What is the minimum Krylov subspace dimension needed before the state start to diverge form one another?
• What happens when we add more impurities? Change the sign of $U_{\text{center}}$?

Documentation Contributors

Pierre-Gabriel Rozon