{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Fermions on the square lattice & perfect nesting\n",
"\n",
"This tutorial is the first in a series of tutorials on two-particle response were we will use TRIQS and the [**Two-Particle Response-Function toolbox (TPRF)**](https://triqs.github.io/tprf/latest/) to compute;\n",
"\n",
"1. the **non-interacting Green's function** of fermions on the square lattice with nearest-neighbour hopping \n",
" and study the Fermi surface [01]\n",
"\n",
"2. the **non-interacting two-particle response**, also called the susceptibility [03]\n",
"\n",
"3. the **Random-Phase Approximation (RPA)** susceptibility for weak interactions, \n",
" studying the anti-ferromagnetic divergence at ($\\pi,\\pi)$ [05]\n",
"\n",
"4. the **Two-Particle Self Conistent (TPSC)** susceptibility\n",
" and show that it satisfies the Pauli principle, while RPA does not [07]\n",
" \n",
"5. and show that TPSC obeys the **Mermin-Wagner theorem**,\n",
" since its spin susceptibility does not diverge at finite temperature [09]\n"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"execution": {
"iopub.execute_input": "2023-08-29T09:09:00.599462Z",
"iopub.status.busy": "2023-08-29T09:09:00.599177Z",
"iopub.status.idle": "2023-08-29T09:09:00.814690Z",
"shell.execute_reply": "2023-08-29T09:09:00.814365Z"
}
},
"outputs": [],
"source": [
"%matplotlib inline\n",
"%config InlineBackend.figure_format = 'svg'\n",
"from triqs.plot.mpl_interface import plt\n",
"\n",
"import numpy as np"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Square lattice with nearest-neighbour hopping\n",
"\n",
"The square lattice with nearest-neighbour hopping $t$ appeared in earier tutorials, were the dispersion relation \n",
"\n",
"\\begin{equation}\n",
" \\epsilon(\\mathbf{k})=-2t(\\cos{k_x}+\\cos{k_y}),\n",
"\\end{equation}\n",
"\n",
"was computed using TRIQS in more than one way.\n",
"\n",
"However, in TRIQS and TPRF there are a number of helper routines for lattice models that simplifies the study of general tight-binding models. Here we will therefore use these standard routines.\n",
"\n",
"Insead of constructing $\\epsilon(\\mathbf{k})$ directly in momentum space we construct a real-space tight binding lattice Hamitonian $H(\\mathbf{r})$ using [the `TBLattice` class](https://triqs.github.io/triqs/latest/documentation/python_api/triqs.lattice.tight_binding.TBLattice.html?highlight=tblattice#triqs.lattice.tight_binding.TBLattice) corresponding to the square lattice with nearest neighbour hopping $t=1$."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"execution": {
"iopub.execute_input": "2023-08-29T09:09:00.816381Z",
"iopub.status.busy": "2023-08-29T09:09:00.816294Z",
"iopub.status.idle": "2023-08-29T09:09:00.867522Z",
"shell.execute_reply": "2023-08-29T09:09:00.867316Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Warning: could not identify MPI environment!\n",
"Tight Binding Hamiltonian on Bravais Lattice with dimension 2, units \n",
"[[1,0,0]\n",
" [0,1,0]\n",
" [0,0,1]], n_orbitals 1\n",
"with hoppings [\n",
" [1,0] : \n",
"[[(-1,0)]]\n",
" [-1,0] : \n",
"[[(-1,0)]]\n",
" [0,1] : \n",
"[[(-1,0)]]\n",
" [0,-1] : \n",
"[[(-1,0)]] ]\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"Starting serial run at: 2023-08-29 11:09:00.866052\n"
]
}
],
"source": [
"from triqs.lattice.tight_binding import TBLattice\n",
"\n",
"t = 1.0 # nearest-neigbhour hopping amplitude\n",
"\n",
"H_r = TBLattice( \n",
" units=[\n",
" (1,0,0), # basis vector in the x-direction \n",
" (0,1,0), # basis vector in the y-direction\n",
" ],\n",
" hoppings={\n",
" (+1,0) : [[-t]], # hopping in the +x direction\n",
" (-1,0) : [[-t]], # hopping in the -x direction\n",
" (0,+1) : [[-t]], # hopping in the +y direction\n",
" (0,-1) : [[-t]], # hopping in the -y direction\n",
" })\n",
"\n",
"print(H_r)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The real-space Hamiltonian $H(\\mathbf{r})$ can both construct a discretized momentum mesh and evaluate the dispersion $\\epsilon(\\mathbf{k})$ on a given momentum mesh, by"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"execution": {
"iopub.execute_input": "2023-08-29T09:09:00.883989Z",
"iopub.status.busy": "2023-08-29T09:09:00.883893Z",
"iopub.status.idle": "2023-08-29T09:09:00.888623Z",
"shell.execute_reply": "2023-08-29T09:09:00.888364Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Greens Function with mesh Brillouin Zone Mesh with linear dimensions (128 128 1)\n",
" -- units = \n",
"[[0.0490874,0,0]\n",
" [0,0.0490874,0]\n",
" [0,0,6.28319]]\n",
" -- brillouin_zone: Brillouin Zone with 2 dimensions and reciprocal matrix \n",
"[[6.28319,0,0]\n",
" [0,6.28319,0]\n",
" [0,0,6.28319]] and target_shape (1, 1): \n",
"\n"
]
}
],
"source": [
"n_k = 128\n",
"kmesh = H_r.get_kmesh(n_k=n_k)\n",
"e_k = H_r.fourier(kmesh)\n",
"print(e_k)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Since the square lattice is two-dimensional it is possible to visualize the dispersion using a color plot.\n",
"\n",
"Here is an example that plots $\\epsilon(\\mathbf{k})$ in the entire Brillouin zone as well as the shape of the Fermi surface at $\\omega = 0$ (dotted line)."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"execution": {
"iopub.execute_input": "2023-08-29T09:09:00.889931Z",
"iopub.status.busy": "2023-08-29T09:09:00.889859Z",
"iopub.status.idle": "2023-08-29T09:09:01.067282Z",
"shell.execute_reply": "2023-08-29T09:09:01.067063Z"
},
"scrolled": false
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
" \n",
" \n",
" \n",
" \n",
" 2023-08-29T11:09:01.035179 \n",
" image/svg+xml \n",
" \n",
" \n",
" Matplotlib v3.7.1, https://matplotlib.org/ \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"k = np.linspace(-np.pi, np.pi, num=100)\n",
"kx, ky = np.meshgrid(k, k)\n",
"\n",
"e_k_interp = np.vectorize(lambda kx, ky : e_k((kx, ky, 0)).real)(kx, ky)\n",
"\n",
"plt.figure()\n",
"\n",
"plt.pcolormesh(kx, ky, e_k_interp, rasterized=True, cmap='RdBu')\n",
"plt.colorbar().ax.set_ylabel(r'$\\epsilon(\\mathbf{k})$'); \n",
"\n",
"plt.contour(kx, ky, e_k_interp, levels=[0], linestyles='dotted')\n",
"\n",
"plt.xlabel(r'$k_x$'); plt.ylabel(r'$k_y$');\n",
"k_ticks, k_labels = [-np.pi, 0, np.pi], [r\"$-\\pi$\", r\"0\", r\"$\\pi$\"]\n",
"plt.xticks(k_ticks, k_labels); plt.yticks(k_ticks, k_labels);\n",
"plt.axis('square');\n",
"\n",
"# -- High-symmetry path G-X-M-G\n",
"\n",
"pts = [\n",
" (0., 0., r'$\\Gamma$', 'w', 'bottom', 'right'), \n",
" (np.pi, 0., r'$X$', 'k', 'center', 'left'),\n",
" (np.pi, np.pi, r'$M$', 'r', 'bottom', 'left'),\n",
" ]\n",
"for x, y, label, color, va, ha in pts:\n",
" plt.plot(x, y, 'o', color=color, clip_on=False, zorder=110)\n",
" plt.text(x, y, label, color=color, fontsize=18, va=va, ha=ha)\n",
"\n",
"X, Y, _, _, _, _ = zip(*(pts+[pts[0]]))\n",
"plt.plot(X, Y, '-m', zorder=100, lw=4);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Momentum dependent quantities can also be visualized along high-symmetry paths in the Brillouin zone, see above for the high-symmetry points $\\Gamma$, $X$ and $M$ of the square lattice. \n",
"\n",
"Here is an example that plots the dispersion $\\epsilon(\\mathbf{k})$ along thepath $\\Gamma - X - M - \\Gamma$ in k-space using [the `triqs_tprf.lattice_utils.k_space_path` function](https://triqs.github.io/tprf/unstable/reference/python_reference.html#triqs_tprf.lattice_utils.k_space_path)."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"execution": {
"iopub.execute_input": "2023-08-29T09:09:01.068722Z",
"iopub.status.busy": "2023-08-29T09:09:01.068648Z",
"iopub.status.idle": "2023-08-29T09:09:01.126302Z",
"shell.execute_reply": "2023-08-29T09:09:01.126033Z"
},
"scrolled": false
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
" \n",
" \n",
" \n",
" \n",
" 2023-08-29T11:09:01.113436 \n",
" image/svg+xml \n",
" \n",
" \n",
" Matplotlib v3.7.1, https://matplotlib.org/ \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"G = [0.0, 0.0, 0.0]\n",
"X = [0.5, 0.0, 0.0]\n",
"M = [0.5, 0.5, 0.0]\n",
"\n",
"path = [(G, X), (X, M), (M, G)]\n",
"\n",
"from triqs_tprf.lattice_utils import k_space_path\n",
"\n",
"k_vecs, k_plot, k_ticks = k_space_path(path, num=32, bz=H_r.bz)\n",
"\n",
"e_k_interp = np.vectorize(lambda k : e_k(k).real, signature='(n)->()')\n",
"\n",
"plt.plot(k_plot, e_k_interp(k_vecs))\n",
"plt.xticks(k_ticks, labels=[r'$\\Gamma$', '$X$', '$M$', r'$\\Gamma$'])\n",
"plt.ylabel(r'$\\epsilon(\\mathbf{k})$')\n",
"plt.grid(True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In the following we will re-purpose these visualization scripts to study the one-particle and two-particle Green's functions of the square lattice model."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Non-interacting lattice Green's function\n",
"\n",
"Given the dispersion $\\epsilon(\\mathbf{k})$ the non-interacting Green's function $G_0(i\\omega_n, \\mathbf{k})$ is given by\n",
"\n",
"\\begin{equation}\n",
" G_0(i\\omega_n, \\mathbf{k}) = \\frac{1}{i\\omega_n - \\epsilon(\\mathbf{k})}\n",
" \\, .\n",
"\\end{equation}\n",
"\n",
"As shown in the Basic Tutorial it is of course possible to compute $G_0$ using a loop over frequency and momentum:\n",
"\n",
"```python\n",
"from triqs.gf import Gf, MeshImFreq, MeshProduct\n",
"\n",
"wmesh = MeshImFreq(beta=2.5, S='Fermion', n_max=128)\n",
"wkmesh = MeshProduct(wmesh, kmesh)\n",
"g0_wk = Gf(mesh=wkmesh, target_shape=e_k.target_shape)\n",
"\n",
"for w, k in wkmesh: \n",
" g0_wk[w, k] = 1/(w - e_k[k])\n",
"```\n",
"\n",
"However, TPRF has Dyson equation solvers that are OpenMP+MPI parallell and all implemented in C++, see [the TPRF documentation](https://triqs.github.io/tprf/latest/reference/python_reference.html#lattice-green-s-functions). Here we will use these fast routines!"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Exercise 1:\n",
"\n",
"Use [`triqs_tprf.lattice.lattice_dyson_g0_wk`](https://triqs.github.io/tprf/latest/reference/python_reference.html#triqs_tprf.lattice.lattice_dyson_g0_wk) to compute $G_0(i\\omega_n, \\mathbf{k})$ at inverse temperature $\\beta = 2.5$ using a [fermionic `MeshImFreq` frequency mesh](https://triqs.github.io/triqs/latest/documentation/python_api/triqs.gf.meshes.MeshImFreq.html?highlight=meshimfreq#triqs.gf.meshes.MeshImFreq) with 128 Matsubara frequencies and name the resulting Green's function `g0_wk`.\n",
"\n",
"Check the properties of `g0_wk` by printing it, i.e.\n",
"```python\n",
"print(g0_wk)\n",
"```"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"execution": {
"iopub.execute_input": "2023-08-29T09:09:01.127961Z",
"iopub.status.busy": "2023-08-29T09:09:01.127871Z",
"iopub.status.idle": "2023-08-29T09:09:01.320874Z",
"shell.execute_reply": "2023-08-29T09:09:01.320627Z"
},
"scrolled": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Greens Function with mesh Imaginary Freq Mesh with beta = 2.5, statistic = Fermion, n_iw = 32, positive_only = false, Brillouin Zone Mesh with linear dimensions (128 128 1)\n",
" -- units = \n",
"[[0.0490874,0,0]\n",
" [0,0.0490874,0]\n",
" [0,0,6.28319]]\n",
" -- brillouin_zone: Brillouin Zone with 2 dimensions and reciprocal matrix \n",
"[[6.28319,0,0]\n",
" [0,6.28319,0]\n",
" [0,0,6.28319]] and target_shape (1, 1): \n",
"\n"
]
}
],
"source": [
"# Write your code here\n",
"\n",
"from triqs.gf import MeshImFreq\n",
"wmesh = MeshImFreq(beta=2.5, S='Fermion', n_iw=32)\n",
"\n",
"from triqs_tprf.lattice import lattice_dyson_g0_wk\n",
"g0_wk = lattice_dyson_g0_wk(mu=0., e_k=e_k, mesh=wmesh)\n",
"\n",
"print(g0_wk)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Questions**\n",
"\n",
"- How many meshes does the Green's function have?\n",
"- What is the k-space discretization?\n",
"- How is the reciprocal basis vectors of the Brillouin zone related to the lattice (`units`) vectors of the tight binding lattice `H_r` above?"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Exercise 2:\n",
"\n",
"Save the Green's function `g0_wk` into a hdf5 file named `g0_wk.h5` using [`h5.HDFArchive`](https://triqs.github.io/triqs/latest/documentation/manual/triqs/hdf5/ref.html), so that we can read `g0_wk` from file in the following tutorials. Do the same with the dispersion `e_k` and save it to a file with the name `e_k.h5`."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"execution": {
"iopub.execute_input": "2023-08-29T09:09:01.322296Z",
"iopub.status.busy": "2023-08-29T09:09:01.322205Z",
"iopub.status.idle": "2023-08-29T09:09:01.459973Z",
"shell.execute_reply": "2023-08-29T09:09:01.459631Z"
}
},
"outputs": [],
"source": [
"# Write your code here\n",
"\n",
"from h5 import HDFArchive\n",
"\n",
"with HDFArchive(\"g0_wk.h5\", \"w\") as R:\n",
" R['g0_wk'] = g0_wk\n",
" \n",
"with HDFArchive(\"e_k.h5\", \"w\") as R:\n",
" R['e_k'] = e_k"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Fermi surface nesting\n",
"\n",
"We will now study the Fermi surface of Fermions on the square lattice with nearest neighbour hopping, which has a special property called *perfect nesting*. A Fermi surface is said to be *nested* if parts of the Fermi surface map to each other by a single momentum vector $\\mathbf{Q}$, called the *nesting vector*.\n",
"\n",
"For a non-interacting system the Fermi surface is the surface in k-space defined by\n",
"\n",
"$$ \\epsilon(\\mathbf{k}) - \\mu = 0 \\, ,$$\n",
"\n",
"where $\\mu$ is the chemical potential. In terms of the spectral function $A(\\omega, \\mathbf{k})$ this corresponds to large values of $A$ at $\\omega=0$\n",
"\n",
"$$ A(\\omega = 0, \\mathbf{k}) = \\frac{1}{\\pi} \\text{Im} \n",
"\\left[ \\frac{1}{ 0 - \\epsilon(\\mathbf{k}) + \\mu - i\\delta } \\right] \\gg 1 \\, ,$$\n",
"\n",
"which also generalizes to interacting systems."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Exercise 3:\n",
"\n",
"Make a color plot of the zero-frequency spectral function $A(k, \\omega=0)$ over the Brillouin zone, using the approximation\n",
"\n",
"$$ A(k, \\omega=0) \\approx -\\frac{1}{\\pi} \\text{Im}[ G_0(\\mathbf{k}, i\\omega_0) ] \\, ,$$\n",
"\n",
"where we neglect the fact that the first fermionic Matsubara frequency $i\\omega_0$ is not exactly $0$.\n",
"\n",
"The right hand side can be evaluated using the Triqs Green's function `g0_wk` and the interpolation feature:\n",
"\n",
"```python\n",
"n = 0 # Matsubara frequency index\n",
"kx, ky, kz = 0., 0., 0.\n",
"k_vec = (kx, ky, kz)\n",
"g0_wk(n, k_vec)\n",
"```"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"execution": {
"iopub.execute_input": "2023-08-29T09:09:01.461536Z",
"iopub.status.busy": "2023-08-29T09:09:01.461453Z",
"iopub.status.idle": "2023-08-29T09:09:01.715649Z",
"shell.execute_reply": "2023-08-29T09:09:01.715396Z"
}
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
" \n",
" \n",
" \n",
" \n",
" 2023-08-29T11:09:01.673693 \n",
" image/svg+xml \n",
" \n",
" \n",
" Matplotlib v3.7.1, https://matplotlib.org/ \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Write your code here\n",
"\n",
"kgrid1d = np.linspace(-np.pi, np.pi, n_k + 1, endpoint=True)\n",
"kx, ky = np.meshgrid(kgrid1d, kgrid1d)\n",
"\n",
"A_k = np.vectorize(lambda kx, ky: -g0_wk( 0, (kx, ky, 0) ).imag / np.pi)\n",
"A_inv_k = np.vectorize(lambda kx, ky: (1 / g0_wk( 0, (kx, ky, 0) )).real)\n",
"\n",
"plt.contour(kx, ky, A_inv_k(kx, ky), levels=[0], colors='white')\n",
"plt.pcolormesh(kx, ky, A_k(kx, ky), rasterized=True)\n",
"\n",
"plt.colorbar().ax.set_ylabel(r\"$G_0(i\\omega_0, \\mathbf{k})$\")\n",
"plt.xticks([-np.pi, 0, np.pi],[r\"$-\\pi$\", r\"0\", r\"$\\pi$\"]) \n",
"plt.yticks([-np.pi, 0, np.pi],[r\"$-\\pi$\", r\"0\", r\"$\\pi$\"])\n",
"plt.axis('square'); plt.xlabel(r\"$k_x$\"); plt.ylabel(r\"$k_y$\");"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Hint**: Re-purpose the code for the color plot of $\\epsilon(\\mathbf{k})$ above.\n",
"\n",
"**Questions**\n",
"\n",
" * How can we see from the plot that the Fermi surface is nested?\n",
" * What is the nesting vector?\n",
" * Actually the Fermi surface is **perfectly** nested. What do you think is the difference between *nesting* and *perfect nesting*?"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Exercise 4:\n",
"\n",
"For non-interacting systems the Fermi surface can also be observed in the Fermionic density distribution in k-space $\\rho(\\mathbf{k})$, which can be computed from the Matsubara Green's function $G_0(i\\omega_n, \\mathbf{k})$.\n",
"\n",
"In TRIQS this is available through the `.density()` method of Matsubara frequency Green's functions. For the lattice Green's function it can be evaluated for a given k-vector using\n",
"\n",
"```python\n",
"kx, ky, kz = 0., 0., 0.\n",
"k_vec = (kx, ky, kz)\n",
"rho_k = g0_wk(all, k_vec).density().real\n",
"```\n",
"\n",
"Plot the momentum distribution $\\rho(\\mathbf{k})$ along the Brillouin zone high-symmetry path $\\Gamma - X - M - \\Gamma$."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"execution": {
"iopub.execute_input": "2023-08-29T09:09:01.717211Z",
"iopub.status.busy": "2023-08-29T09:09:01.717131Z",
"iopub.status.idle": "2023-08-29T09:09:01.772482Z",
"shell.execute_reply": "2023-08-29T09:09:01.772191Z"
}
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
" \n",
" \n",
" \n",
" \n",
" 2023-08-29T11:09:01.760870 \n",
" image/svg+xml \n",
" \n",
" \n",
" Matplotlib v3.7.1, https://matplotlib.org/ \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"rho_k_interp = np.vectorize(lambda k: g0_wk(all, k).density(), signature='(n)->()')\n",
"\n",
"plt.plot(k_plot, rho_k_interp(k_vecs).real)\n",
"plt.xticks(k_ticks, labels=[r'$\\Gamma$', '$X$', '$M$', r'$\\Gamma$'])\n",
"plt.ylabel(r'$\\rho(\\mathbf{k})$')\n",
"plt.grid(True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Hint**: Re-purpose the code above plotting $\\epsilon(\\mathbf{k})$ along the same path.\n",
"\n",
"**Questions**\n",
"\n",
"- What is the value of $\\rho(\\mathbf{k})$ at the Fermi surface?\n",
"- What is the sign of $\\epsilon(\\mathbf{k})$ in the regions of k-space where $\\rho(\\mathbf{k}) \\approx 1$ and $\\approx 0$, respectively?\n",
"- How is this related to the Fermi distribution function $$f(\\epsilon) = \\frac{1}{1 + e^{\\beta \\epsilon}}$$ plotted below?"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"execution": {
"iopub.execute_input": "2023-08-29T09:09:01.773973Z",
"iopub.status.busy": "2023-08-29T09:09:01.773898Z",
"iopub.status.idle": "2023-08-29T09:09:01.825794Z",
"shell.execute_reply": "2023-08-29T09:09:01.825547Z"
}
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
" \n",
" \n",
" \n",
" \n",
" 2023-08-29T11:09:01.813545 \n",
" image/svg+xml \n",
" \n",
" \n",
" Matplotlib v3.7.1, https://matplotlib.org/ \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"beta = 2.5\n",
"e = np.linspace(-4., 4.)\n",
"f = lambda e : 1/(1 + np.exp(beta * e))\n",
"plt.plot(e, f(e))\n",
"plt.xlabel(r'$\\epsilon$'); plt.ylabel(r'$f(\\epsilon)$'); plt.grid(True);"
]
}
],
"metadata": {
"@webio": {
"lastCommId": null,
"lastKernelId": null
},
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.11.5"
}
},
"nbformat": 4,
"nbformat_minor": 4
}