touch up the hubbard model

docs/source/tutorial/hubbard_1d.md

@@ -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.