{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Finite temperature antiferromagnetism in two dimensional systems\n",
"\n",
"## Mermin-Wagner theorem\n",
"\n",
"The Mermin-Wagner theorem [[Mermin and Wagner, Phys. Rev. Lett. 17, 1133 (1966)]](https://doi.org/10.1103/PhysRevLett.17.1133) states that a continuous-symmetry cannot be broken at finite temperature in two dimensions. So, finite temperature antiferromagnetism is impossible in two dimensions, contrary to the prediction of RPA. \n",
"\n",
"Here is a heuristic proof. Assume that the system has long range order and that the spins are collinear in the $z$-direction. Then the free energy density contains a term proportional to $(\\nabla S^z)^2/2$, where $S^z$ is the spin density operator in the $z$-direction. In Fourier space, this becomes $-q^2|S^z(\\mathbf{q})|^2/2$. In the long wavelength limit, where fluctuations are slow, we can use the classical equipartition theorem which gives \n",
"$$\n",
"|S^z(\\mathbf{q})|^2=\\frac{k_BT}{q^2} \\, .\n",
"$$\n",
"However, the denominator gives rise to a so called *infrared divergence* which cause the local moment to diverge \n",
"$$\n",
"\\left<(S^z)^2\\right>\\sim \\int d^2q |S^z(\\mathbf{q})|^2\\sim \\int d^2q\\frac{k_BT}{q^2}=\\infty \\, .\n",
"$$\n",
"We come to an absurdity, which proves that our initial hypothesis of the existence of long range order is wrong. Hence, there is no long-range order. \n",
"\n",
"## TPSC and the Mermin-Wagner theorem\n",
"\n",
"To see that TPSC satisfies the Mermin-Wagner theorem, we first note that the spin susceptibility has the following spectral representation\n",
"\n",
"$$\\chi_{sp}(\\mathbf{q},i\\omega_n)=\\int\\frac{d\\omega}{\\pi}\\frac{\\chi_{sp}''(\\mathbf{q},\\omega)}{\\omega-i\\omega_n}=\\int\\frac{d\\omega}{\\pi}\\frac{\\chi_{sp}''(\\mathbf{q},\\omega)\\omega}{\\omega^2+(\\omega_n)^2}.$$\n",
"\n",
"The last equality follows from the fact $\\chi_{sp}''(\\mathbf{q},\\omega)$ is odd in frequency. This last result shows that the finite Matsubara frequencies should be regular and that the largest contribution is at the zeroth-Matsubara frequency. This allows us to give a rough idea of why the theorem is satisfied by focusing on the zero Matsubara frequency contribution. \n",
"\n",
"Let us write the self-consistency condition for $U_{sp}$ as follows\n",
"\n",
"\\begin{equation}\n",
"\\frac{T}{N}\\sum_{\\mathbf{q}} \\frac{\\chi_0(\\mathbf{q},0)}{1-\\frac{U_{sp}}{2}\\chi_0(\\mathbf{q},0)}=n-2\\left< n_\\uparrow n_\\downarrow\\right>-C(T)\n",
"\\end{equation}\n",
"\n",
"where $C(T)$ contains the non-singular contribution of the finite Matsubara frequencies. \n",
"\n",
"Calling the right-hand side $C'(T)$, expanding the denominator around the maximum at $\\mathbf{Q}=(\\pi,\\pi)$ and shifting the origin of the wave vector integration to $\\mathbf{Q}=(\\pi,\\pi)$, the self-consistency condition can be written as \n",
"\n",
"\\begin{equation}\n",
"\\frac{T}{N}\\sum_{\\mathbf{q}} \\frac{A}{\\xi^{-2}+q^2}=C'(T)\n",
"\\end{equation}\n",
"where $A$ is a constant and $\\xi$ the correlation length contains the value of $U_{sp}$. Since the right-hand side is finite, $\\xi$ adjusts itself not to become infinite, otherwise the left-hand side diverges. The divergence of the susceptibility can occur only at $T=0$ where we cannot treat the non-zero Matsubara frequencies separately. \n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Code from previous notebooks\n",
"\n",
"To study the temperature dependence we will reuse the code from the previous TPSC notebook. Please look through the functions and make sure that they are familiar."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"execution": {
"iopub.execute_input": "2023-08-29T09:09:45.339540Z",
"iopub.status.busy": "2023-08-29T09:09:45.339071Z",
"iopub.status.idle": "2023-08-29T09:09:45.613081Z",
"shell.execute_reply": "2023-08-29T09:09:45.612838Z"
}
},
"outputs": [],
"source": [
"%matplotlib inline\n",
"%config InlineBackend.figure_format = 'svg'\n",
"from triqs.plot.mpl_interface import plt, oplot\n",
"\n",
"import numpy as np\n",
"from triqs.gf import *"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## $\\chi_0$ calculator for arbitrary $\\beta$\n",
"\n",
"Since we have to recompute the bare bubble $\\chi_0$ for every temperature we also need an implementation `get_chi0` of the bubble calculation, see below."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"execution": {
"iopub.execute_input": "2023-08-29T09:09:45.615131Z",
"iopub.status.busy": "2023-08-29T09:09:45.614607Z",
"iopub.status.idle": "2023-08-29T09:09:45.621718Z",
"shell.execute_reply": "2023-08-29T09:09:45.621523Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Warning: could not identify MPI environment!\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"Starting serial run at: 2023-08-29 11:09:45.617709\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",
"kmesh = H_r.get_kmesh(n_k=16)\n",
"e_k = H_r.fourier(kmesh)\n",
"\n",
"from triqs.gf import MeshImFreq\n",
"from triqs_tprf.lattice import lattice_dyson_g0_wk\n",
"from triqs_tprf.lattice_utils import imtime_bubble_chi0_wk\n",
"\n",
"def get_chi0(beta, n_iw=128):\n",
" wmesh = MeshImFreq(beta=beta, S='Fermion', n_iw=64)\n",
" g0_wk = lattice_dyson_g0_wk(mu=0., e_k=e_k, mesh=wmesh)\n",
" chi0_wk = 2 * imtime_bubble_chi0_wk(g0_wk, nw=32, verbose=False)\n",
" return chi0_wk"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here are the functions used for the TPSC self consistency."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"execution": {
"iopub.execute_input": "2023-08-29T09:09:45.638123Z",
"iopub.status.busy": "2023-08-29T09:09:45.638026Z",
"iopub.status.idle": "2023-08-29T09:09:45.726470Z",
"shell.execute_reply": "2023-08-29T09:09:45.726133Z"
}
},
"outputs": [],
"source": [
"from scipy.optimize import brentq\n",
"\n",
"from triqs_tprf.lattice import solve_rpa_PH\n",
"\n",
"def solve_rpa(chi0_wk, U):\n",
" I = np.ones([1, 1, 1, 1], dtype=complex)\n",
" chi_wk = solve_rpa_PH(chi0_wk, U/2 * I)\n",
" return chi_wk\n",
"\n",
"def trace_chi(chi_wk):\n",
" \"\"\"Compute the sum, \\sum_k \\sum_\\nu \\chi(k,\\nu)\"\"\" \n",
" wmesh, kmesh = chi_wk.mesh.components\n",
" chi_w = Gf(mesh=wmesh, target_shape=[])\n",
" chi_w.data[:] = np.sum(np.squeeze(chi_wk.data), axis=1) / len(kmesh)\n",
" return -chi_w.density().real\n",
"\n",
"def Usp_root(Usp, chi0_wk, n, U):\n",
" tr_chi_sp = trace_chi(solve_rpa(chi0_wk, U=Usp))\n",
" diff = tr_chi_sp + 0.5 * Usp/U * n**2 - n\n",
" return diff\n",
"\n",
"def Uch_root(Uch, chi0_wk, n, U, docc):\n",
" tr_chi = trace_chi(solve_rpa(chi0_wk, U=-Uch))\n",
" diff = tr_chi - 2 * docc - n + n**2\n",
" return diff\n",
"\n",
"def solve_tpsc(chi0_wk, U, n):\n",
" Uc = 2/np.squeeze(chi0_wk(0, [np.pi, np.pi, 0])).real\n",
" Usp = brentq(Usp_root, 0, Uc, args=(chi0_wk, n, U), xtol=1e-2)\n",
" docc = 0.25 * Usp / U * n**2\n",
" Uch = brentq(Uch_root, 0, 100, args=(chi0_wk, n, U, docc), xtol=1e-2)\n",
" return Usp, Uch, docc, Uc"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Exercise\n",
"\n",
"Solve the TPCS equations for a range of temperatures at $n=1$ (half-filling) and $U=4$ and examine the validity of the Mermin-Wagner theorem in the TPSC approximation. In particular study the spin structure factor $S(\\mathbf{q})$ defined as \n",
"\n",
"$$\n",
"S(\\mathbf{q})\\equiv T\\sum_n \\chi_{sp}(\\mathbf{q},i\\omega_n)\n",
"$$\n",
"\n",
"and reproduce the results of following figure from __[[Vilk and Tremblay, J. Phys. I France 7 (1997) 1309-1368]](https://jp1.journaldephysique.org/articles/jp1/abs/1997/11/jp1v7p1309/jp1v7p1309.html)__ __[[arXiv version]](https://arxiv.org/abs/cond-mat/9702188v3)__:\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## RPA spin structure factor $S_{RPA}$ as a function of temperature $T$\n",
"\n",
"For comparison we compute the RPA spin structure factor $S_{RPA}$ for a range of temperatures. Note that $T_c^{(RPA)} \\approx 0.75$."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"execution": {
"iopub.execute_input": "2023-08-29T09:09:45.728211Z",
"iopub.status.busy": "2023-08-29T09:09:45.728125Z",
"iopub.status.idle": "2023-08-29T09:09:46.076876Z",
"shell.execute_reply": "2023-08-29T09:09:46.076646Z"
},
"scrolled": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"| T | beta | S_rpa |\n",
"-----------------------------------------\n",
"| 4.0000E+00 | 2.5000E-01 | 1.7196E-01 |\n",
"| 3.0000E+00 | 3.3333E-01 | 2.5925E-01 |\n",
"| 2.8000E+00 | 3.5714E-01 | 2.8783E-01 |\n",
"| 2.6000E+00 | 3.8462E-01 | 3.2302E-01 |\n",
"| 2.4000E+00 | 4.1667E-01 | 3.6725E-01 |\n",
"| 2.2000E+00 | 4.5455E-01 | 4.2425E-01 |\n",
"| 2.0000E+00 | 5.0000E-01 | 5.0006E-01 |\n",
"| 1.8000E+00 | 5.5556E-01 | 6.0503E-01 |\n",
"| 1.6000E+00 | 6.2500E-01 | 7.5858E-01 |\n",
"| 1.4000E+00 | 7.1429E-01 | 1.0020E+00 |\n",
"| 1.2000E+00 | 8.3333E-01 | 1.4436E+00 |\n",
"| 1.0000E+00 | 1.0000E+00 | 2.5028E+00 |\n",
"| 8.0000E-01 | 1.2500E+00 | 9.9329E+00 |\n"
]
}
],
"source": [
"U = 4.0\n",
"n = 1.0\n",
"\n",
"T_rpa_vec = np.concatenate((np.arange(4., 3., -1.), np.arange(3, 0.75, -0.2)))\n",
"S_rpa_vec = np.zeros_like(T_rpa_vec)\n",
"\n",
"print(''.join('| %-11s' % s for s in ['T', 'beta', 'S_rpa']), '|')\n",
"print('-'*41)\n",
"\n",
"for idx, T in enumerate(T_rpa_vec):\n",
"\n",
" beta = 1. / T\n",
" chi0_wk = get_chi0(beta)\n",
" chi_wk = solve_rpa(chi0_wk, U)\n",
" \n",
" S_rpa = -chi_wk(all, (np.pi, np.pi, 0))[0,0,0,0].density().real * beta\n",
" S_rpa_vec[idx] = S_rpa\n",
" \n",
" print(''.join('| %4.4E ' % x for x in [T, beta, S_rpa_vec[idx]]), '|') "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## TPSC spin structure factor $S_{TPSC}$ as a function of temperature\n",
"\n",
"Using the ansatz $U_{sp}\\left \\left=U\\left$, the spin susceptibility obeys\n",
"\n",
"\\begin{equation}\n",
"\\frac{T}{N}\\sum_{\\mathbf{q},i\\omega_n} \\frac{\\chi_0(\\mathbf{q},i\\omega_n)}{1-\\frac{U\\left}{2\\left \\left}\\chi_0(\\mathbf{q},i\\omega_n)}=n-2\\left< n_\\uparrow n_\\downarrow\\right>\n",
"\\end{equation}\n",
"\n",
"When the susceptibility increases, $\\left$ on the right-hand side decreases, but then the denominator of the spin susceptibility will lead to a decrease in susceptibility.\n",
"\n",
"More rigorously, we can see that dimension is important here. Let us repeat the argument at the beginning of the notebook. The right-hand side of the equation cannot diverge. Also, on the left-hand side, note that the most divergent contribution is the zero Matsubara frequency, as one can see from the spectral representation and $\\chi''(\\mathbf{q},\\omega)=-\\chi''(\\mathbf{q},-\\omega)$\n",
"\n",
"\\begin{equation}\n",
"\\chi(\\mathbf{q},i\\omega_n)=\\int \\frac{d\\omega}{\\pi}\\frac{\\chi''(\\mathbf{q},\\omega)}{\\omega-i\\omega_n}=\\int \\frac{d\\omega}{\\pi}\\frac{\\omega\\chi''(\\mathbf{q},\\omega)}{\\omega^2+\\omega_n^2}.\n",
"\\end{equation}\n",
"\n",
"Using these results, the non-singular finite Matsubara frequency terms can be put on the right-hand side of the sum rule and all that is left is \n",
"\n",
"\\begin{equation}\n",
"T\\int d^2q \\frac{a}{\\xi^{2}+q^2}\\sim C'(T)\n",
"\\end{equation}\n",
"\n",
"where we have expanded the susceptibility around $(\\pi,\\pi)$, gone from sum to integral and shifted the origin of integration so that now $\\mathbf{q}$ is the deviation from $(\\pi,\\pi)$. On dimensional grounds, the left-hand side is logarithmic in two dimensions so that the correlation length scales like $\\exp(C'(T)/T)$.\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Exercise: Spin structure factor\n",
"\n",
"Compute the TPSC spin structure factor for a range of temperatures $T \\in [0.25, 4]$ and plot $S_{TPSC}$ and $S_{RPA}$ and determine whether the Mermin-Wagner theorem holds."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"execution": {
"iopub.execute_input": "2023-08-29T09:09:46.078314Z",
"iopub.status.busy": "2023-08-29T09:09:46.078233Z",
"iopub.status.idle": "2023-08-29T09:09:47.375999Z",
"shell.execute_reply": "2023-08-29T09:09:47.375724Z"
},
"scrolled": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"| T | beta | Usp | Uch | docc | S_tpsc |\n",
"--------------------------------------------------------------------------------\n",
"| 4.0000E+00 | 2.5000E-01 | 3.1151E+00 | 4.9735E+00 | 1.9469E-01 | 1.6108E-01 |\n",
"| 3.0000E+00 | 3.3333E-01 | 2.9065E+00 | 5.3054E+00 | 1.8166E-01 | 2.3189E-01 |\n",
"| 2.5000E+00 | 4.0000E-01 | 2.7718E+00 | 5.5801E+00 | 1.7324E-01 | 2.9464E-01 |\n",
"| 2.0000E+00 | 5.0000E-01 | 2.6148E+00 | 6.0225E+00 | 1.6343E-01 | 3.9851E-01 |\n",
"| 1.5000E+00 | 6.6667E-01 | 2.4421E+00 | 6.8411E+00 | 1.5263E-01 | 5.9599E-01 |\n",
"| 1.2000E+00 | 8.3333E-01 | 2.3406E+00 | 7.7133E+00 | 1.4629E-01 | 8.2155E-01 |\n",
"| 1.0000E+00 | 1.0000E+00 | 2.2802E+00 | 8.5589E+00 | 1.4251E-01 | 1.0735E+00 |\n",
"| 8.0000E-01 | 1.2500E+00 | 2.2301E+00 | 9.6547E+00 | 1.3938E-01 | 1.5004E+00 |\n",
"| 6.0000E-01 | 1.6667E+00 | 2.1949E+00 | 1.0900E+01 | 1.3718E-01 | 2.3575E+00 |\n",
"| 4.0000E-01 | 2.5000E+00 | 2.1657E+00 | 1.2228E+01 | 1.3535E-01 | 4.8974E+00 |\n",
"| 3.5000E-01 | 2.8571E+00 | 2.1527E+00 | 1.2658E+01 | 1.3454E-01 | 6.6014E+00 |\n",
"| 3.0000E-01 | 3.3333E+00 | 2.1327E+00 | 1.3218E+01 | 1.3330E-01 | 1.0230E+01 |\n",
"| 2.5000E-01 | 4.0000E+00 | 2.0869E+00 | 1.4279E+01 | 1.3043E-01 | 2.2999E+01 |\n"
]
}
],
"source": [
"# Write your code here\n",
"\n",
"T_tpsc_vec = np.array([ 4., 3., \n",
" 2.5, 2.0, 1.5, 1.2, 1.0, \n",
" 0.8, 0.6, 0.4, 0.35, 0.3, 0.25, \n",
" ])\n",
"\n",
"S_tpsc_vec = np.zeros_like(T_tpsc_vec)\n",
"U_sp_vec = np.zeros_like(T_tpsc_vec)\n",
"U_ch_vec = np.zeros_like(T_tpsc_vec)\n",
"docc_vec = np.zeros_like(T_tpsc_vec)\n",
"\n",
"print(''.join('| %-11s' % s for s in ['T', 'beta', 'Usp', 'Uch', 'docc', 'S_tpsc']), '|')\n",
"print('-'*80)\n",
"\n",
"for idx, T in enumerate(T_tpsc_vec):\n",
"\n",
" beta = 1. / T \n",
" chi0_wk = get_chi0(beta) \n",
" \n",
" Usp, Uch, docc, UcRPA = solve_tpsc(chi0_wk, U, n)\n",
" \n",
" chi_sp_wk = solve_rpa(chi0_wk, Usp)\n",
" S_tpsc = -chi_sp_wk(all, (np.pi, np.pi, 0))[0,0,0,0].density().real * beta\n",
"\n",
" S_tpsc_vec[idx], U_sp_vec[idx], U_ch_vec[idx], docc_vec[idx] = S_tpsc, Usp, Uch, docc \n",
" print(''.join('| %4.4E ' % x for x in [T, beta, Usp, Uch, docc, S_tpsc]), '|') "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Excercise: Critical temperature and double occupancy\n",
"\n",
"To see the divergencies it is useful to also study the inverse spin structure factor $S^{-1}$. Plot $S^{-1}$ and see at what temperatures the curves intercept $S^{-1}=0$ to determine the critical temperatures $T_c$ of RPA and TPSC. Also plot the double occupancy and explain its behaviour as a function of temperature."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"execution": {
"iopub.execute_input": "2023-08-29T09:09:47.377434Z",
"iopub.status.busy": "2023-08-29T09:09:47.377361Z",
"iopub.status.idle": "2023-08-29T09:09:47.678024Z",
"shell.execute_reply": "2023-08-29T09:09:47.677790Z"
}
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
"