Skip to content
Snippets Groups Projects

Examples

Merged Examples
examples into main
Merged Kostas Vilkelis requested to merge examples into main
Viewing commit 57aab479
Prev Next
Show latest version
1 file
+ 89
− 44
Compare changes
  • Side-by-side
  • Inline
docs/source/tutorial/hubbard_1d.md
+ 89
− 44
@@ -10,39 +10,60 @@ kernelspec:
language: python
name: python3
---
# 1D Hubbard model
# 1D Hubbard on spin chain
## Background physics
To simulate how the package can be used to simulate infinite systems, we show how to use it with a tight-binding model in 1 dimension.
We exemplify this by computing the ground state of an infinite spinful chain with onsite interactions.
First, the basic imports are done.
To show the basic functionality of the package, we consider a simple interacting electronic system: a 1D chain of sites that allow nearest-neighbor tunneling with strength $t$ and on-site repulsion $U$ between two electrons if they are on the same site.
Such a model is known as the 1D [Hubbard model](https://en.wikipedia.org/wiki/Hubbard_model) and is useful for understanding the onset of insulating phases in interacting metals.
```{code-cell} ipython3
import numpy as np
import matplotlib.pyplot as plt
import pymf
```
To begin, we first consider the second quantized form of the non-interacting Hamiltonian.
Because we expect the interacting ground state to be antiferromagnetic, we build a two-atom cell and name the two sublattices $A$ and $B$.
These sublattices are identical to each other in the non-interacting case $U=0$.
The non-interacting Hamiltonian reads:
$$
\hat{H_0} = - t \sum_\sigma \sum_i \left(c_{i, B, \sigma}^{\dagger}c_{i, A, \sigma} + c_{i, A, \sigma}^{\dagger}c_{i+1, B, \sigma} + \textrm{h.c}\right).
$$
After this, we start by constructing the non-interacting Hamiltonian. As we expect the ground state to be an antiferromagnet, we build a two-atom cell. We name the two sublattices, $A$ and $B$. The Hamiltonian is then:
where $\textrm{h.c}$ is the hermitian conjugate, $\sigma$ denotes spin ($\uparrow$ or $\downarrow$) and $c_{i, A, \sigma}^{\dagger}$ creates an electron with spin $\sigma$ in unit cell $i$ of sublattice $A$.
Next up, is the interacting part of the Hamiltonian:
$$
H_0 = \sum_i c_{i, B}^{\dagger}c_{i, A} + c_{i, A}^{\dagger}c_{i+1, B} + h.c.
\hat{V} = U \sum_i \left(n_{i, A, \uparrow} n_{i, A, \downarrow} + n_{i, B, \uparrow} n_{i, B, \downarrow}\right).
$$
We write down the spinful part by simply taking $H_0(k) \otimes \mathbb{1}$.
where $n_{i, A, \sigma} = c_{i, A, \sigma}^{\dagger}c_{i, A, \sigma}$ is the number operator for sublattice $A$ and spin $\sigma$.
The total Hamiltonian is then $\hat{H} = \hat{H_0} + \hat{V}$.
With the model defined, we can now proceed to input the Hamiltonian into the package and solve it using the mean-field approximation.
## Problem definition
To translat ethis Hamiltonian into a tight-binding model, we specify the hopping vectors together with the hopping amplitudes. We ensure that the Hamiltonian is hermitian by letting the hopping amplitudes from $A$ to $B$ be the complex conjugate of the hopping amplitudes from $B$ to $A$.
### Non-interacting Hamiltonian
First, lets get the basic imports out of the way.
```{code-cell} ipython3
import numpy as np
import matplotlib.pyplot as plt
import pymf
```
Now lets translate the non-interacting Hamiltonian $\hat{H_0}$ defined above into a basic input format for the package: a **tight-binding dictionary**.
The tight-binding dictionary is a python dictionary where the keys are tuples of integers representing the hopping vectors and the values are the hopping matrices.
For example, a key `(0,)` represents the onsite term and a key `(1,)` represents the hopping a single unit cell to the right.
In the case of our 1D Hubbard model, non-interacting Hamiltonian is:
```{code-cell} ipython3
hopp = np.kron(np.array([[0, 1], [0, 0]]), np.eye(2))
h_0 = {(0,): hopp + hopp.T.conj(), (1,): hopp, (-1,): hopp.T.conj()}
```
Here `hopp` is the hopping matrix which we define as a kronecker product between sublattice and spin degrees of freedom: `np.array([[0, 1], [0, 0]])` corresponds to the hopping between sublattices and `np.eye(2)` leaves the spin degrees of freedom unchanged.
In the corresponding tight-binding dictionary `h_0`, the key `(0,)` contains hopping within the unit cell and the keys `(1,)` and `(-1,)` correspond to the hopping between the unit cells to the right and left respectively.
We verify this tight-binding model by plotting the band structure and observing the two bands due to the Brillouin zone folding. In order to do this we transform the tight-binding model into a Hamiltonian on a k-point grid using the `tb_to_kgrid` and then diagonalize it.
We verify the validity of `h_0`, we evaluate it in the reciprocal space using the {autolink}`~pymf.tb.transforms.tb_to_kgrid`, diagonalize it and plot the band structure:
```{code-cell} ipython3
# Set number of k-points
nk = 100
nk = 50 # number of k-points
ks = np.linspace(0, 2*np.pi, nk, endpoint=False)
hamiltonians_0 = pymf.tb_to_kgrid(h_0, nk)
@@ -53,29 +74,41 @@ plt.xlim(0, 2 * np.pi)
plt.ylabel("$E - E_F$")
plt.xlabel("$k / a$")
plt.show()
```
After confirming that the non-interacting part is correct, we can set up the interacting Hamiltonian. We define the interaction similarly as a tight-binding dictionary. To keep the physics simple, we let the interaction be onsite only, which gives us the following interaction matrix:
which seems metallic as expected.
$$
H_{int} =
\left(\begin{array}{cccc}
U & U & 0 & 0\\
U & U & 0 & 0\\
0 & 0 & U & U\\
0 & 0 & U & U
\end{array}\right)~.
$$
### Interaction Hamiltonian
This is simply constructed by writing:
We now proceed to define the interaction Hamiltonian $\hat{V}$.
To achieve this, we utilize the same tight-binding dictionary format as before.
Because the interaction Hamiltonian is on-site, it must be defined only for the key `(0,)` and only for electrons on the same sublattice with opposite spins.
Based on the kronecker product structure we defined earlier, the interaction Hamiltonian is:
```{code-cell} ipython3
U=0.5
h_int = {(0,): U * np.kron(np.eye(2), np.ones((2,2))),}
U=2
s_x = np.array([[0, 1], [1, 0]])
h_int = {(0,): U * np.kron(np.eye(2), s_x),}
```
Here `s_x` is the Pauli matrix acting on the spin degrees of freedom, which ensures that the interaction is only between electrons with opposite spins whereas the `np.eye(2)` ensures that the interaction is only between electrons on the same sublattice.
In order to find a mean-field solution, we combine the non interacting Hamiltonian with the interaction Hamiltonian and the relevant filling into a `Model` object. We then generate a starting guess for the mean-field solution and solve the model using the `solver` function. It is important to note that the guess will influence the possible solutions which the `solver` can find in the mean-field procedure. The `generate_guess` function generates a random Hermitian tight-binding dictionary, with the keys provided as hopping vectors and the values of the size as specified. We specifically choose the keys, meaning the hopping vectors, for the `guess` to be the same as the hopping vectors for the interaction Hamiltonian. This is because we do not expect the mean-field solution to contain terms more than the hoppings from the interacting part.
### Putting it all together
To combine the non-interacting and interaction Hamiltonians, we use the {autolink}`~pymf.model.Model` class.
In addition to the Hamiltonians, we also need to specify the filling of the system --- the number of electrons per unit cell.
```{code-cell} ipython3
filling = 2
full_model = pymf.Model(h_0, h_int, filling)
```
The object `full_model` now contains all the information needed to solve the mean-field problem.
## Solving the mean-field problem
To find a mean-field solution, we first require a starting guess.
In cases where the non-interacting Hamiltonian is highly degenerate, there exists several possible mean-field solutions, many of which are local and not global minima of the energy landscape.
Here the problem is simple enough that we can generate a random guess for the mean-field solution through the {autolink}`~pymf.tb.utils.generate_guess` function.
Finally, to solve the model, we use the {autolink}`~pymf.solvers.solver` function which by default employes a root-finding algorithm to find a self-consistent mean-field solution.
```{code-cell} ipython3
filling = 2
@@ -84,7 +117,26 @@ guess = pymf.generate_guess(frozenset(h_int), ndof=4)
mf_sol = pymf.solver(full_model, guess, nk=nk)
```
The `solver` function returns only the mean-field correction to the non-interacting Hamiltonian. To get the full Hamiltonian, we add the mean-field correction to the non-interacting Hamiltonian. To take a look at whether the result is correct, we first do the mean-field computation for a wider range of $U$ values and then plot the gap as a function of $U$.
The {autolink}`~pymf.solvers.solver` function returns only the mean-field correction to the non-interacting Hamiltonian in the same tight-binding dictionary format.
To get the full Hamiltonian, we add the mean-field correction to the non-interacting Hamiltonian and plot the band structure just as before:
```{code-cell} ipython3
h_mf = pymf.add_tb(h_0, mf_sol)
hamiltonians = pymf.tb_to_kgrid(h_mf, nk)
vals, vecs = np.linalg.eigh(hamiltonians)
plt.plot(ks, vals, c="k")
plt.xticks([0, np.pi, 2 * np.pi], ["$0$", "$\pi$", "$2\pi$"])
plt.xlim(0, 2 * np.pi)
plt.ylabel("$E - E_F$")
plt.xlabel("$k / a$")
plt.show()
```
the band structure now shows a gap at the Fermi level, indicating that the system is in an insulating phase!
We can go further and compute the gap for a wider range of $U$ values:
```{code-cell} ipython3
def compute_sol(U, h_0, nk, filling=2):
@@ -120,21 +172,14 @@ def compute_phase_diagram(
vals.append(_vals)
return np.asarray(gaps, dtype=float), np.asarray(vals)
```
Now we run this function over a larger range of $U$ values and plot the gap as a function of $U$.
Us = np.linspace(0, 4, 40, endpoint=True)
gap, vals = compute_phase_diagram(Us=Us, nk=20, nk_dense=100)
```{code-cell} ipython3
Us = np.linspace(0, 8, 40, endpoint=True)
nk, nk_dense = 40, 100
gap, vals = compute_phase_diagram(Us=Us, nk=nk, nk_dense=nk_dense)
```
```{code-cell} ipython3
plt.plot(Us, gap, c="k")
plt.xlabel("$U$")
plt.ylabel("Gap")
plt.xlabel("$U / t$")
plt.ylabel("$\Delta{E}/t$")
plt.show()
```
As expected, there is a critical $U$ around $U=1$ where the gap opens. Furthermore, at large $U$, the gap scales linearly with $U$.
We see that at around $U=1$ the gap opens up and the system transitions from a metal to an insulator.