{ "cells": [ { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Nanostructures in Equilibrium with Python\n", "\n", "© Prof. Branislav K. Nikolić & Dr. Marko D. Petrović, University of Delaware\n", "\n", "[PHYS824: Nanophysics & Nanotechnology](https://wiki.physics.udel.edu/phys824) " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## What is covered in this notebook:\n", "- Charge and spin density in 1D nanowire with or without impurities\n", "- DOS of disordered 1D nanowire via eigenvalue binning and Anderson localization of its eigenstates\n", "- DOS of 1D closed or open nanowire via Green functions\n", "- Equilibrium density matrix of open 1D nanowire via Green functions and Ozaki contour complex integration\n", "- Hofstadter butterfly energy spectrum of electrons on a square lattice in magnetic field\n", "- DOS of graphene flake via Green functions" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import scipy.linalg as la\n", "import matplotlib.pyplot as plt\n", "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1. Charge and spin density in 1D nanostructures" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### 1.1 Charge density\n", "Compute electron charge density for a 1D nanowire via equilibrium grand canonical density matrix. The nanowire is modeled using [discretized continuous Hamiltonian](https://wiki.physics.udel.edu/phys824/Discretization_of_1D_continuous_Hamiltonian) with possible impurity potential on its diagonal:\n", "\n", "$\\hat{H} = \\sum_{i=1}^n \\big[(U_i+2t)|i\\rangle \\langle i| - t|i\\rangle \\langle i+1| -t|i+1 \\rangle \\langle i| \\big]$\n", "\n", "Identical tight-binding Hamiltonian (without $2t$ on the diagonal) appears also in modeling of [1D chain of atoms](https://wiki.physics.udel.edu/wiki_phys824/images/8/87/PHYS824_lecture3_atom_to_1D_nanowires.pdf) with one orbital included per atom.\n", "\n", "__NOTE: In order for this example to run correctely, all the cells have to be executed one after another.__\n", "\n", "Define physical constants" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "HBAR = 1.055e-34 # Planck constant [J*s]\n", "M_EL = 9.110e-31 # Electron mass [kg]\n", "Q_EL = 1.602e-19 # Electron charge [C]\n", "mu = 0.25 # Chemical potential \n", "kbT = 0.025; # Temperature in eV; useful to remember is that room temperature T=300 K is 0.025eV" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "and lattice properties" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "a = 2e-10 # Lattice spacing in meters \n", "n = 100 # Number of discrete lattice points \n", "x = np.arange(0, n*a, a) # Lattice site positions\n", "t = (HBAR**2)/(2*M_EL*(a**2))/Q_EL # Hopping energy " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Load Hamiltonian matrix representation\n", "\n", "$\\hat{H} = \\begin{bmatrix} 2t & -t & 0 & \\cdots & 0 \\\\\n", " -t & 2t & -t & \\cdots & 0 \\\\\n", " 0 & -t & 2t & \\cdots & 0 \\\\\n", " \\cdots & \\cdots & \\cdots & 2t & -t \\\\\n", " 0 & \\cdots & \\cdots & -t & 2t\n", " \\end{bmatrix}$\n", "\n", "into NumPy 2D array via shortcut offered by `numpy.diag`" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "ham = 2*t*np.diag(np.ones(n)) # Create diagonal elements \n", "ham += -t*np.diag(np.ones(n-1), -1) # Add lower diagonal hoppings\n", "ham += -t*np.diag(np.ones(n-1), 1) # Add upper diagonal hoppings" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Check shape" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ham.shape" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Introduce possible onsite potential and/or periodic boundary conditions" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "# Uncomment to add periodic boundary conditions\n", "# ham[0, -1] = -t\n", "# ham[-1, 0] = -t\n", "pot = np.zeros((n, n))\n", "# Uncomment to add impurity in the middle\n", "# pot_impurity = np.zeros(n)\n", "# pot_impurity[int(n/2)] += 0.1\n", "# pot = np.diag(pot_impurity)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Check that Hamiltonian is Hermitian matrix" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "la.norm(ham-np.conjugate(ham.T))" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Compute eigenvalues and eigenvectors" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "evals, evecs = np.linalg.eigh(ham + pot)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Find equilibrium density matrix\n", "\n", "$\\hat{\\rho}_\\mathrm{eq} = f(\\hat{H} - \\mu \\hat{I})$\n", "\n", "as the Fermi function of the Hamiltonian" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "rho = 1. / (1 + np.exp((evals-mu)/kbT))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from functools import reduce\n", "rho1 = reduce(np.matmul, [evecs, np.diag(rho), np.conj(evecs.T)]) \n", "density = np.diag(rho1) / a " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plot the results" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [], "source": [ "import matplotlib as mpl # These two lines are for setting DPI \n", "mpl.rcParams['figure.dpi'] = 100\n", "\n", "plt.plot(x, density)\n", "plt.xlabel(r'$x$ (m)')\n", "plt.ylabel(r'Electron density (/$m^3$)')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "For running this as standalone Python code, we can put all cells above into a function" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "from functools import reduce\n", "\n", "HBAR = 1.055e-34 # Planck constant [J*s]\n", "M_EL = 9.110e-31 # Electron mass [kg]\n", "Q_EL = 1.602e-19 # Electron charge [C]\n", "\n", "def adj(ar):\n", " return np.conj(ar.T)\n", "\n", "def density(n=100, mu=0.025, u_impurity=0.0, \n", " kbT=0.025, periodic=False):\n", " a = 2e-10 # Lattice spacing in meters\n", " x = np.arange(0, n*a, a) # Lattice site positions\n", " # t = (HBAR**2)/(2*M_EL*(a**2))/Q_EL # Hopping energy \n", " t = 1\n", " # Define Hamiltonian\n", " ham = 2*t*np.diag(np.ones(n)) # Create diagonal elements \n", " ham += -t*np.diag(np.ones(n-1), -1) # Add lower diagonal\n", " ham += -t*np.diag(np.ones(n-1), 1) # Add upper diagonal\n", " \n", " if periodic:\n", " ham[0, -1] = -t\n", " ham[-1, 0] = -t\n", " \n", " pot_impurity = np.zeros(n)\n", " pot_impurity[int(n/2)] += u_impurity\n", " pot = np.diag(pot_impurity)\n", " \n", " evals, evecs = np.linalg.eigh(ham + pot)\n", " rho = 1. / (1 + np.exp((evals-mu)/kbT))\n", " \n", " rho1 = reduce(np.matmul, [evecs, np.diag(rho), np.conj(evecs.T)]) \n", " density = np.diag(rho1) / a \n", " return x, evals, density" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Define plotting function " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "def plot_density(x, density):\n", " %matplotlib inline\n", " from matplotlib import pyplot as plt \n", " import matplotlib as mpl # These two lines are for setting DPI \n", " mpl.rcParams['figure.dpi'] = 100\n", "\n", " plt.plot(x, density)\n", " plt.xlabel(r'$x$ (m)')\n", " plt.ylabel(r'Electron density (/$m^3$)')\n", " plt.show()" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Test the function by comparing results with those obtained by executing individual cells above" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x, evals, dens = density(n=100, mu=0.25, u_impurity=0.0,\n", " kbT=0.0025, periodic=False)\n", "plot_density(x, dens)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### 1.2 Spin density \n", "\n", "For the same 1D nanowire, let us compute equilibrium electron spin density in the presence of __classical magnetic impurity__ within the wire. First, our tight-binding Hamiltonian has to be extended into spinfull Hilbert space, $\\mathcal{H}_\\mathrm{orbital} \\otimes \\mathcal{H}_\\mathrm{spin}$, of a single electron\n", "\n", "$\\hat{H}= \\big[-t\\sum_{i=1}^n \\big[|i\\rangle\\langle i+1| + |i+1\\rangle\\langle i|)\\big] \\otimes \\hat{I}_\\mathrm{spin}$\n", "\n", "Here the unit operator $\\hat{I}_\\mathrm{spin}$ in spin space will become just `np.eye(2)` in Python. Second, we add magnetic impurity at one site, as described by its classical magnetization vector $\\mathbf{m}=(\\sin\\theta\\cos\\phi,\\sin \\theta \\sin\\phi,\\cos \\theta)$, which interacts with quantum spin of an electron via the so-called $sd$ exchange interaction \n", "\n", "$\\hat{H} = \\left[-t\\sum_{i=1}^n (|i\\rangle\\langle i+1| + |i+1\\rangle\\langle i|)\\right] \\otimes \\hat{I}_\\mathrm{spin} - |p\\rangle\\langle p| \\otimes J_{sd} \\hat{\\boldsymbol{\\sigma}} \\cdot \\mathbf{m}$\n", "\n", "Treating magnetic impurity spin as classical vector means we are assuming its spin is large. Note that the eigenvalue of spin operatore squared in quantum mechanics $\\hat{\\mathbf{S}}^2$ is $S(S+1)=S^2(1+1/S)$ which deviates from \"classical spin\" $S^2$ by that term $1/S$ which becomes smaller and smaller correction as $S \\rightarrow \\infty$. " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "code_folding": [], "slideshow": { "slide_type": "slide" } }, "outputs": [], "source": [ "# Define Pauli matrices for spin-1/2 of electron:\n", "SIGMA_X = np.array([[0, 1],\n", " [1, 0]])\n", "SIGMA_Y = np.array([[0, -1j],\n", " [1j, 0]])\n", "SIGMA_Z = np.array([[1, 0],\n", " [0, -1]])\n", "UNITY = np.array([[1, 0],\n", " [0, 1]])\n", "\n", "\n", "def density_matrix(n=100, mu=0.025, kbT=0.025, periodic=False, \n", " moment=(1, 0, 0), jsd=0.1,\n", " mag_impurity_index=49):\n", " \n", " a = 2e-10 # Lattice spacing in meters\n", " x = np.arange(0, n*a, a) # Lattice site positions\n", " \n", " # Define Hamiltonian\n", " ham = 2*t*np.diag(np.ones(n, dtype=complex)) # Create diagonal elements \n", " ham += -t*np.diag(np.ones(n-1), -1) # Add lower diagonal\n", " ham += -t*np.diag(np.ones(n-1), 1) # Add upper diagonal\n", " if periodic: # Add periodic boundary conditions\n", " ham[0, -1] = -t\n", " ham[-1, 0] = -t\n", " ham = np.kron(ham, UNITY) # Make tight-binding Hamiltonian spinfull\n", " \n", " # Add magnetic impurity potential with impurity magnetization pointing along (mx,my,mz) vector\n", " mm = np.array(moment) / np.linalg.norm(moment) # Impurity magnetization vector should be a unit vector\n", " mx, my, mz = mm\n", " pot = -jsd * (mx*SIGMA_X + my*SIGMA_Y + mz*SIGMA_Z)\n", " ind = 2*mag_impurity_index\n", " ham[ind:ind+2, ind:ind+2] += pot\n", " \n", " # Find equilibrium density matrix\n", " evals, evecs = np.linalg.eigh(ham)\n", " rho_zero = 1. / (1 + np.exp((evals-mu)/kbT))\n", " rho = reduce(np.matmul, [evecs, np.diag(rho_zero), np.conj(evecs.T)]) \n", " rho /= a\n", " \n", " return x, rho" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Function to compute spin and charge density" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "def get_density(rho):\n", " n_sites = int(rho.shape[0]/2)\n", " Tr = np.trace # If a function name is too long\n", " mul = np.matmul # you can rename it\n", " \n", " charge = []\n", " sx = []\n", " sy = []\n", " sz = []\n", " for i in range(n_sites):\n", " rho_site = rho[2*i:2*i+2, 2*i:2*i+2]\n", " charge.append(Tr(mul(rho_site, UNITY)))\n", " sx.append(Tr(mul(rho_site, SIGMA_X)))\n", " sy.append(Tr(mul(rho_site, SIGMA_Y)))\n", " sz.append(Tr(mul(rho_site, SIGMA_Z)))\n", " \n", " # Convert lists into numpy arrays\n", " charge = np.array(charge)\n", " sx = np.array(sx)\n", " sy = np.array(sy)\n", " sz = np.array(sz)\n", " return charge.real, sx.real, sy.real, sz.real" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Function to plot multiplot figures with matplotlib" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "import matplotlib as mpl # These two lines are for setting DPI \n", "mpl.rcParams['figure.dpi'] = 400\n", "#%matplotlib inline\n", "\n", "def plot_4panel(x, n, sx, sy, sz):\n", " fig, axs = plt.subplots(2, 2)\n", " axs[0, 0].plot(x, n, color='black', lw=0.7)\n", " axs[0, 0].set_title(r'$n$ vs. $x$')\n", " axs[0, 0].set_xlabel(r'$x$ (m)')\n", " axs[0, 0].set_ylabel(r'Electron density ($m^{-1}$)')\n", "\n", " axs[0, 1].plot(x, sx, color='C0', lw=0.7)\n", " axs[0, 1].set_title(r'${\\rm S}_{\\rm x}$ vs. $x$')\n", " axs[0, 1].set_xlabel(r'$x$ (m)')\n", " axs[0, 1].set_ylabel(r'Spin density ($m^{-1}$)')\n", "\n", " axs[1, 0].plot(x, sy, color='C2', lw=0.7)\n", " axs[1, 0].set_title(r'${\\rm S}_{\\rm y}$ vs. $x$')\n", " axs[1, 0].set_xlabel(r'$x$ (m)')\n", " axs[1, 0].set_ylabel(r'Spin density ($m^{-1}$)')\n", "\n", " axs[1, 1].plot(x, sz, color='C3', lw=0.7)\n", " axs[1, 1].set_title(r'${\\rm S}_{\\rm z}$ vs. $x$')\n", " axs[1, 1].set_xlabel(r'$x$ (m)')\n", " axs[1, 1].set_ylabel(r'Spin density ($m^{-1}$)')\n", "\n", " fig.tight_layout(h_pad=1, w_pad=0.2) # Use this to control the spacing\n", "\n", " plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Test function for different orientations of impurity magnetization $\\mathbf{m}$" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [], "source": [ "x, rho = density_matrix(n=100, mu=0.25, kbT=0.025, periodic=False, \n", " moment=(1, 1, 1), jsd=-2.0, mag_impurity_index=49)\n", "n, sx, sy, sz = get_density(rho)\n", "plot_4panel(x, n, sx, sy, sz)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## 2. DOS of 1D disordered nanowire via eigenvalue binning\n", "\n", "We use random on-site potential $U_i \\in [-W/2,W/2]$ on every site of tight-binding Hamiltonian:\n", "\n", "$\\hat{H} = \\sum_{i=1}^n \\big[U_i|i\\rangle \\langle i| -t |i\\rangle \\langle i+1| -t |i+1 \\rangle \\langle i| \\big]$\n", "\n", "\n", "where increasing $W$ yields an example of [Anderson localization](https://physicstoday.scitation.org/doi/10.1063/1.3206091) of wavefunctions in 1D lattice which is manifestation of quantum intereference effects. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from numpy.random import rand\n", "\n", "def disorder_1d(n=500, t=1, periodic=True, w=0.1):\n", " \n", "\n", " # Hamiltonian matrix for clean wire\n", " ham = np.diag(np.zeros(n)) # Create diagonal elements \n", " ham += -t*np.diag(np.ones(n-1), -1) # Add lower diagonal\n", " ham += -t*np.diag(np.ones(n-1), 1) # Add upper diagonal\n", " \n", " if periodic:\n", " ham[0, -1] = -t\n", " ham[-1, 0] = -t\n", "\n", " # Add Anderson disorder = uniform random variable as\n", " # on-site potential on every atom\n", " pot = w*t # Anderson disorder strength\n", "\n", " # Impurity potential is uniform random variable in the interval [-W/2,W/2]\n", " u = -0.5*pot + pot*rand(n)\n", " ham = ham + np.diag(u)\n", " \n", " evals, evecs = np.linalg.eigh(ham)\n", "\n", " return evals, evecs\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "evals, evects = disorder_1d(w=5.0)\n", "\n", "fig, ax = plt.subplots(1, 3)\n", "ax[0].plot(evals)\n", "ax[0].set_title('Energy Eigenvalues')\n", "ax[0].set_xlabel('Eigenvalue No.')\n", "ax[0].set_ylabel('Energy (eV)')\n", "\n", "ax[1].hist(evals, bins=80)\n", "ax[1].set_ylabel('DOS')\n", "ax[1].set_xlabel('Energy (eV)')\n", "\n", "ax[2].plot(evects[:, 220], color='C3', lw=0.3) # Eigenfunction for state 220\n", "ax[2].set_title('$\\Psi(x)$ for Eigenvalue No. 220')\n", "ax[2].set_xlabel('$x$')\n", "ax[2].set_ylabel(r'$\\Psi(x)$')\n", "\n", "fig.tight_layout(h_pad=1, w_pad=0.2) # Use this to control the spacing\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## 3. DOS of 1D nanowire via Green functions\n", "\n", "### 3.1 Closed 1D nanowire" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The following function computes DOS in a finite or infinite 1D nanowire via Green functions. The DOS \n", "$${\\rm DOS}(E) = -\\frac{1}{\\pi} {\\rm Tr}\\left[{\\rm Im}\\mathbf{G}^{r}(E)\\right]$$\n", "is computed using retarded Green function\n", "$$ \\mathbf{G}^{r}(E) = \\left[(E + i\\eta)\\mathbf{I} - \\mathbf{H}\\right]^{-1}$$\n", "whose imaginary part is\n", "$$ \\mathrm{Im}\\mathbf{G}^r = (\\mathbf{G}^{r}-\\mathbf{G}^{a})/2i$$" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "code_folding": [] }, "outputs": [], "source": [ "from numpy.random import rand\n", "import numpy as np\n", "\n", "\n", "def dos_closed(n=100, t=1.0, periodic=True, eta=1.e-3,\n", " emin=-5, emax=5, ne=100):\n", " \"\"\"Compute density of states in the nanowire\n", " \n", " Parameters\n", " ----------\n", " n : int, optional\n", " Number of sites in the nanowire\n", " t : float, optional\n", " Hopping constant between neighboring sites\n", " periodic : boolean, optional\n", " If True, set hopping between the first and the last site.\n", " eta : float, optional\n", " Small number used to obtain retarded Green's function.\n", " emin : float, optional\n", " Lower energy limit for DOS computation\n", " emax : float, optional\n", " Upper energy limit for DOS computation\n", " ne : int, optional\n", " Number of energy points used in DOS computation\n", " \n", " Returns\n", " -------\n", " ens : numpy 1D array of floats\n", " Energies at which the function computes DOS.\n", " dos : numpy 1D array of floats\n", " DOS for the given energies.\n", " \"\"\" \n", "\n", " # Hamiltonian matrix for a clean wire\n", " ham = np.zeros((n, n), dtype=complex)\n", " ham += -t*np.diag(np.ones(n-1), -1) # Add lower diagonal\n", " ham += -t*np.diag(np.ones(n-1), 1) # Add upper diagonal\n", " \n", " if periodic:\n", " ham[0, -1] = -t\n", " ham[-1, 0] = -t\n", " \n", " energies = np.linspace(emin, emax, int(ne))\n", " \n", " dos = []\n", " for energy in energies:\n", " green = np.linalg.inv((energy + 1j*eta)*np.eye(n) - ham)\n", " dos_e = (-1/np.pi)*np.trace((green - np.conj(green.T))/2j)\n", " dos.append(dos_e)\n", " \n", " return energies, np.array(dos)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "\n", "en, dos = dos_closed(eta=1.e-3, n=100, periodic=False, ne=500)\n", "plt.plot(en, dos)\n", "plt.xlabel('Energy (t)')\n", "plt.ylabel('DOS')\n", "plt.title('DOS in finite 1D nanowire')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 3.2 Open 1D nanowire" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First we define a function that computes self-energy of a semi-infinite 1D lead from an analytical expression\n", "$$\\boldsymbol{\\Sigma}^r(E) = \n", " \\left\\{ \n", " \\begin{array}{lr}\n", " \\frac{t_{\\textrm{C}}^2}{2t_{\\textrm{L}}^2}\n", " \\left(E - i\\sqrt{4t_{\\textrm{L}}^2 - E^2}\\right), & \n", " \\textrm{for }|E|\\lt 2|t_{\\textrm{L}}| \\\\\n", " \\frac{t_{\\textrm{C}}^2}{2t_{\\textrm{L}}^2}\n", " \\left(E - \\textrm{sgn}(E)\n", " \\sqrt{E^2 - 4t_{\\textrm{L}}^2}\\right), & \n", " \\textrm{for }|E|\\ge 2|t_{\\textrm{L}}| \n", " \\end{array}\n", " \\right. $$\n", "Here, $t_{\\textrm{L}}$ is the hopping between neighboring sites inside of the semi-infinite lead and $t_{\\textrm{C}}$ is the hopping between the first site of the lead and the central region as finit length 1D nanowire. The retarded Green function of open wire is then given by\n", "$$\\mathbf{G}^{r}(E) = \\left[E\\mathbf{I} - \\mathbf{H} - \\boldsymbol{\\Sigma}^r_L(E) -\\boldsymbol{\\Sigma}^r_R(E) \\right]^{-1}$$\n", "if both left (L) and right (R) leads are attached. Note that $\\boldsymbol{\\Sigma}^r_L(E)=\\boldsymbol{\\Sigma}^r_R(E)=\\boldsymbol{\\Sigma}^r(E)$ from the expression above, assuming that $t_{\\textrm{L}}$ and $t_{\\textrm{C}}$ are the same for both L and R leads." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def sigma(energy, t_couple=1.0, t_lead=1.0):\n", " \"\"\"Compute selfenergy of a semi-infinite lead\n", " \n", " Parameters\n", " ----------\n", " energy : float\n", " Energy at which to compute the selefenergy.\n", " t_couple : float, optional\n", " Hopping between the 1D wire and the\n", " semi-infinite lead (the default is 1.0)\n", " t_lead : float, optional\n", " Hopping between sites withing the 1D \n", " semi-infinite lead (the default is 1.0)\n", " \n", " Returns\n", " -------\n", " se : float\n", " The selfenergy for the semi-infinite lead\n", " \"\"\"\n", " \n", " rad = 4*t_lead**2 - energy**2\n", " coef = t_couple**2 / (2*t_lead**2)\n", " \n", " if (rad > 0):\n", " se = coef * (energy - 1j*np.sqrt(rad)) \n", " else:\n", " se = coef * (energy - np.sign(energy)*np.sqrt(-rad))\n", " \n", " return se" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Test sigma function by printing its value for energy inside and outside of the band" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(sigma(2.0))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from numpy.random import rand\n", "import numpy as np\n", "\n", "\n", "def dos_open(n=100, t=1.0, t_lead=1.0, t_couple=1.0,\n", " emin=-5, emax=5, ne=100, attach='both'):\n", " \"\"\"Compute density of states in an open 1D nano-wire\n", " \n", " Parameters\n", " ----------\n", " n : int, optional\n", " Number of sites in the nanowire.\n", " t : float, optional\n", " Hopping constant between neighboring sites.\n", " t_lead : float, optional\n", " Hopping in the leads. Both leads have identical hopping.\n", " t_couple : float, optional\n", " Hopping in between the leads and the wire.\n", " attach : string, optional\n", " Can be set to \n", " \"left\" - attach only left lead, \n", " \"right\" - attach only right lead\n", " \"both\" - attach both leads.\n", " emin : float, optional\n", " Lower energy limit for DOS computation.\n", " emax : float, optional\n", " Upper energy limit for DOS computation.\n", " ne : int, optional\n", " Number of energy points used in DOS computation.\n", " \n", " Returns\n", " -------\n", " ens : numpy 1D array of floats\n", " Energies at which the function computes DOS.\n", " dos : numpy 1D array of floats\n", " DOS for the given energies.\n", " ldos : numpy 1D array of floats\n", " LDOS at the central site.\n", " \"\"\" \n", " # Hamiltonian matrix for a clean wire\n", " ham = np.zeros((n, n), dtype=complex)\n", " n_half = int(n/2)\n", " ham += -t*np.diag(np.ones(n-1), -1) # Add lower diagonal\n", " ham += -t*np.diag(np.ones(n-1), 1) # Add upper diagonal\n", " \n", " energies = np.linspace(emin, emax, int(ne))\n", " \n", " dos = []\n", " ldos = []\n", " \n", " for energy in energies:\n", " \n", " selfenergy = sigma(energy, t_couple, t_lead)\n", " \n", " # Pay attention, we are making a copy of ham matrix \n", " ham_eff = np.array(ham) \n", " \n", " if attach == 'left':\n", " ham_eff[0, 0] += selfenergy\n", " elif attach == 'right':\n", " ham_eff[-1, -1] += selfenergy\n", " elif attach == 'both':\n", " ham_eff[0, 0] += selfenergy\n", " ham_eff[-1, -1] += selfenergy\n", " else:\n", " msg = 'Unknown option %s' % attach\n", " raise ValueError(msg)\n", " \n", " green = np.linalg.inv(energy*np.eye(n) - ham_eff)\n", " im_g = -(green - np.conj(green.T)) / (2j*np.pi)\n", " ldos.append(im_g[n_half, n_half])\n", " dos_e = (-1/np.pi)*np.trace((green - np.conj(green.T))/2j)\n", " dos.append(dos_e)\n", " \n", " return energies, np.array(dos), np.array(ldos)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def plot_2panel(en, dos, ldos):\n", " fig, ax = plt.subplots(2, 1)\n", " ax[0].plot(en, ldos)\n", " ax[0].set_xlabel('Energy (t)')\n", " ax[0].set_ylabel('LDOS')\n", " ax[0].set_title('LDOS in open 1D nanowire')\n", " \n", " ax[1].plot(en, dos)\n", " ax[1].set_xlabel('Energy (t)')\n", " ax[1].set_ylabel('DOS')\n", " ax[1].set_title('DOS in open 1D nanowire')\n", "\n", " fig.tight_layout(h_pad=1, w_pad=0.2) \n", " plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "\n", "en, dos, ldos = dos_open(n=100, ne=500, attach='right')\n", "plot_2panel(en, dos, ldos)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 4. Equilibrium density matrix of open 1D nanowire via Green functions and Ozaki contour complex integration" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The equilibrium density matrix can be expressed as\n", "\\begin{equation}\n", " \\label{egRho}\n", " \\boldsymbol{\\rho}_{\\rm eq}(E_F) = -{\\rm Im}\\left[\n", " \\frac{1}{\\pi} \\int_{-\\infty}^{\\infty}\n", " \\mathbf{G}^r(E)f(E, E_F)\\,dE \\right],\n", "\\end{equation}\n", "where $G^r(E)$ is the same retarded Green function of an open active region attached to two semi-infinite leads. \n", "To compute the given integral, one can use the so-called [Ozaki contour](https://doi.org/10.1103/PhysRevB.75.035123). The main idea behind this method is to approximate the Fermi-Dirac function\n", "\\begin{equation}\n", " f(x) = \\frac{1}{1 + \\exp(x)},\n", "\\end{equation}\n", "where the variable $x$ is expressed in terms of the\n", "temperature and the Fermi energy $x = {(E - E_F)} / k_{B}T$.\n", "After aproximation, one has to compute the poles and residues\n", "of the approximate function. The function is approximated by\n", "continued fraction representation\n", "\\begin{equation}\n", " \\frac{1}{1+\\exp(x)} = \\frac{1}{2} - \\frac{x}{4}\n", " \\left(\\frac{1}{1+\\frac{\\frac{x^2}{2}}{3+\n", " \\frac{\\frac{x^2}{2}}{5 + \\ldots}}} \\right).\n", "\\end{equation}\n", "\n", "In the Ozaki method, the equilibrium density matrix is expressed through poles\n", "$\\alpha_{p}$ and residues $R_{p}$ of the approximate Fermi-Dirac function is\n", "\\begin{equation}\n", " \\boldsymbol{\\rho}_{\\rm eq}(E_F) = -\\frac{1}{\\pi}{\\rm Im}\\left[\n", " 2\\pi i k_{B}T \\sum_{p=1}^{N}\n", " \\mathbf{G}(\\alpha_{p})R_{p} \\right] + iR\\mathbf{G}(iR),\n", "\\end{equation}\n", "For continued fraction approximation function, the poles and residues can be computed numericaly. \n", "The first $N$ poles of function are calculated by\n", "solving the eigenvalue problem of an $N\\times N$ matrix $B$ \n", "\\begin{equation}\n", " B_{n,n+1} = B_{n+1,n} =\n", " \\frac{1}{\\sqrt{(2n+1)(2n-1)}},\\,\\,n\\geq1.\n", "\\end{equation}\n", "The positive eigenvalues of $B$\n", "\\begin{equation}\n", " B |b_{p}\\rangle = b_{p}\n", " |b_{p}\\rangle, \\quad b_{p} > 0\n", "\\end{equation}\n", "determine the positions of poles in the upper-half\n", "of the complex plane\n", "\\begin{equation}\n", " z_{p} = \\frac{1}{b_p}i,\n", "\\end{equation}\n", "while the first elements of the corresponding eigenvectors\n", "$|b_{p}\\rangle$ determine the residues\n", "\\begin{equation}\n", " R_{p} = \\frac{|\\langle\n", " 1|b_{p}\\rangle|^{2}}{4b_{p}^{2}}.\n", "\\end{equation}\n", "When integrating over the energy axis, we must perform a\n", "substitution\n", "\\begin{equation}\n", " z_p = \\frac{\\alpha_{p} - E_{F}}{k_{B}T},\n", "\\end{equation}\n", "from which we obtain\n", "\\begin{eqnarray}\n", " \\alpha_{p} & = & E_F + z_{p}k_{B}T \\nonumber \\\\\n", " & = & E_{F} + \\frac{k_{B}T}{b_{p}}i.\n", "\\end{eqnarray}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First, we want to define a function that computes Ozaki poles and residues" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "from scipy.sparse import diags\n", "from scipy.sparse.linalg import eigsh\n", "\n", "def ozaki_poles(n_poles):\n", " '''Compute the first n poles of the Ozaki continued\n", " fraction representation of the Fermi-Dirac function\n", " [Phys. Rev. B 75, 035123 (2007)]. The implementation\n", " presented here is adopted from [Phys. Rev. B 82, 125114 (2010)]\n", "\n", " Parameters\n", " ----------\n", " n_poles : int\n", " Number of poles to compute\n", "\n", " Returns\n", " -------\n", " zeros : numpy array of doubles\n", " An array of n_poles elements with the singular points of the\n", " Ozaki continued fraction representation\n", " res : numpy array of doubles \n", " An array of residues at positions of zeros.\n", " '''\n", "\n", " n = n_poles*2 + 1\n", " x = np.arange(1, n, 1)\n", " x = 1.0 / (2*np.sqrt(4*x**2 - 1))\n", " diag_elements = [x, x]\n", "\n", " b = diags(diag_elements, [-1, 1])\n", " evals, evects = eigsh(b, k=n-1)\n", "\n", " # Function eigsh returns the eigenvectors in columnwise\n", " # order, so evects[0] is an array of the first elements\n", " # and not the first eigenvector.\n", " res = np.abs(evects[0])**2 / (4*evals**2)\n", "\n", " res = res[np.where(evals > 0)]\n", " evals = evals[np.where(evals > 0)]\n", "\n", " zeros = 1.0 / evals\n", "\n", " # The strange indexing in this part is used just to\n", " # reverse the ordering in these two arrays\n", " return zeros[::-1], res[::-1]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we define Green function for 1D chain" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def green_1d(sites=100, energy=0.0, t=1., t_lead=1.,\n", " t_ls=1., onsite=0.0, u=0.25, imp_index=50):\n", "\n", " h = (-t*(np.eye(sites, k=1) + \n", " np.eye(sites, k=-1))).astype(complex)\n", " \n", " pot_diag = np.zeros(sites)\n", " pot_diag[imp_index] = u\n", " h += np.diag(pot_diag)\n", " rad = 4*t_lead**2 - energy**2\n", " self_energy = 0*1j\n", " if abs(energy.real) < 2*t_lead:\n", " self_energy = energy - 1j*np.sqrt(rad)\n", " else:\n", " if energy.real > 0:\n", " sgn = +1.0\n", " else:\n", " sgn = -1.0\n", " self_energy = energy - sgn*np.sqrt(-rad)\n", " self_energy *= (t_ls**2 / (2*t_lead**2))\n", "\n", " h[0, 0] += self_energy\n", " h[-1, -1] += self_energy\n", "\n", " h += (onsite * np.eye(sites))\n", " return (np.linalg.inv(energy*np.eye(sites) - h))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def ozaki(green_func, npoles=20, rd=1.e+20, energy=-3., kbt=0.25):\n", " \n", " beta = 1. / (kbt)\n", " \n", " poles, res = ozaki_poles(npoles)\n", "\n", " sum_ozaki = 0.0\n", " for z, r in zip(poles, res): # zip function is used to iterate over pairs\n", " pole = energy + 1j*z/beta\n", " sum_ozaki += (r*green_func(energy=pole))\n", " sum_ozaki = 2j*sum_ozaki/beta\n", " rho = (sum_ozaki - np.conjugate(sum_ozaki.T)) / 2j\n", " moment = 0.5j * rd * green_func(energy=1j*rd)\n", " \n", " return rho + moment" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "rho_ozaki = ozaki(green_1d, npoles=500, energy=-1.75, kbt=0.025)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "a = 2.e-10\n", "x, evals, dens = density(n=100, mu=0.25, u_impurity=0.25,\n", " kbT=0.025, periodic=True)\n", "# Compare Ozaki density with density in periodic closed system\n", "plt.plot(np.diag(rho_ozaki)/a, color='C1', lw=3.5);\n", "plt.plot(dens, color='C0', ls=(1, (2, 2)), lw=3.5);\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 5. Hofstadter butterfly energy spectrum of electrons on square lattice in magnetic field" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import math\n", "from numpy import arange\n", "t = 1 # Nearest-neighbor hopping between atoms\n", "hamhop = np.zeros((3,3)) # Hopping submatrix between different columns of atoms \n", "hamhop = -t*np.diag(np.ones(3)) # Add hoppings on diagonal" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x_data_plot = []\n", "y_data_plot = []\n", "for phi in arange(-math.pi,math.pi,0.1):\n", "\n", " ham11 = np.zeros((3,3),dtype=complex) # Submatrix of the first row of atoms\n", " ham11 += -t*np.exp(-2*1j*phi)*np.diag(np.ones(n-1), -1) # Add lower diagonal of complex hoppings\n", " ham11 += -t*np.exp(2*1j*phi)*np.diag(np.ones(n-1), 1) # Add upper diagonal of complex hoppings\n", " \n", " ham22 = np.zeros((3,3),dtype=complex) # Submatrix of the second row of atoms\n", " ham22 += -t*np.exp(-4*1j*phi)*np.diag(np.ones(n-1), -1) # Add lower diagonal of complex hoppings\n", " ham22 += -t*np.exp(4*1j*phi)*np.diag(np.ones(n-1), 1) # Add upper diagonal of complex hoppings\n", " \n", " ham33 = np.zeros((3,3),dtype=complex) # submatrix of the third row of atoms\n", " ham33 += -t*np.exp(-6*1j*phi)*np.diag(np.ones(n-1), -1) # Add lower diagonal of complex hoppings\n", " ham33 += -t*np.exp(6*1j*phi)*np.diag(np.ones(n-1), 1) # Add upper diagonal of complex hoppings\n", " \n", " # load submatrices ham11,ham22,ham33,hamhop into full Hamiltonian matrix\n", " ham = np.zeros((9,9),dtype=complex)\n", " ham[0:3,0:3] += ham11\n", " ham[3:6,3:6] += ham22\n", " ham[6:9,6:9] += ham33\n", " ham[0:3,3:6] += hamhop\n", " ham[3:6,6:9] += hamhop\n", " ham[3:6,0:3] += hamhop\n", " ham[6:9,3:6] += hamhop\n", " \n", " eigvals, eigvecs = la.eig(ham)\n", " \n", " for j in range(0,9):\n", " x_data_plot.append(phi)\n", " y_data_plot.append(eigvals[j].real)\n", "\n", "# Add features to our figure\n", "plt.scatter(x_data_plot,y_data_plot,s=10) # s defines the size of points in scatter plot\n", "plt.title('Eigenenergies of 3x3 tight-binding Hamiltonian in magnetic field')\n", "plt.xlabel('Magnetic Flux Ba$^2$ (h/e)')\n", "plt.ylabel('Eigenenergy (eV)')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 6. DOS of graphene flake via Green functions" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "def dos_graphene(nx=7, ny=4, gamma=1.0, v=0.0,\n", " periodic=False, eta=1.e-3, emin=-3.5, emax=3.5, ne=100):\n", " \"\"\"Compute density of states in graphene ribbon\n", " \n", " Parameters\n", " ----------\n", " nx : integer, optional\n", " Number of carbon atoms along the x-axis.\n", " ny : integer, optional\n", " Number of carbon atoms along the y-axis.\n", " gamma : float, optional\n", " Hopping between nearest neighbors. The basic unit of energy.\n", " periodic : boolean, optional\n", " Introduce periodic boundry conditions along the y-axis.\n", " v : float, optional\n", " Onsite potential energy in units of hopping.\n", " eta : float, optional\n", " Small number used to obtain retarded Green's function.\n", " emin : float, optional\n", " Lower energy limit for DOS computation in the units of hopping.\n", " emax : float, optional\n", " Upper energy limit for DOS computation in the units of hopping.\n", " ne : integer\n", " Number of points in energy for DOS computation.\n", " \"\"\"\n", "\n", " # Define hamiltonians of individual rows of atoms\n", " # h_11 or h_22\n", " # C C\n", " # | \n", " # C C\n", " # | \n", " # C C ^ y\n", " # | | \n", " # C C |\n", " # ------> x\n", " \n", " h_size = nx * ny # Size of the full Hamiltonian matrix\n", " ham = np.zeros((h_size, h_size), dtype=complex)\n", " t = -gamma\n", " v = v*gamma\n", " \n", " h_11 = (v*np.diag(np.ones(ny)) +\n", " t*np.diag(np.ones(ny-1), 1) +\n", " t*np.diag(np.ones(ny-1), -1))\n", " h_22 = np.array(h_11) # Make a copy. Caution! h_22 = h_11 will create an alias\n", " \n", " # Remove odd hoppings from h_11\n", " for i in range(1, ny-1, 2):\n", " h_11[i, i+1] = 0.\n", " h_11[i+1, i] = 0.\n", " \n", " # Remove even hoppings from h_22\n", " for i in range(0, ny-1, 2):\n", " h_22[i, i+1] = 0.\n", " h_22[i+1, i] = 0.\n", " \n", " # Matrices for coupling between two rows of atoms\n", " h_ij = t*np.diag(np.ones(ny)) \n", " \n", " if periodic:\n", " # Wrap system along the y-axis\n", " h_11[0, -1] = t\n", " h_11[-1, 0] = t\n", " h_22[0, -1] = t\n", " h_22[-1, 0] = t\n", " \n", " # Wrap system along the x-axis\n", " ham[:ny, -ny:] = h_ij\n", " ham[-ny:, :ny] = h_ij\n", " \n", " # Put h_11 and h_22 as diagonal blocks of the main hamiltonian\n", " for i in range(0, nx, 2):\n", " ham[i*ny:(i+1)*ny, i*ny:(i+1)*ny] = h_11\n", " \n", " for i in range(1, nx, 2):\n", " ham[i*ny:(i+1)*ny, i*ny:(i+1)*ny] = h_22\n", " \n", " # Put the coupling matrices as the main off-diagonal blocks\n", " hx = (np.diag(np.ones(nx-1), 1) \n", " + np.diag(np.ones(nx-1), -1))\n", " ham = ham + np.kron(hx, h_ij)\n", " \n", " energies = np.linspace(emin, emax, ne)\n", " \n", " dos = []\n", " for energy in energies:\n", " green = np.linalg.inv((energy + 1j*eta)*np.eye(h_size) - ham)\n", " dos_e = (-1/np.pi)*np.trace((green - np.conj(green.T))/2j)\n", " dos.append(dos_e)\n", " \n", " return energies, np.array(dos)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ens, dos = dos_graphene(nx=7, ny=4, eta=0.1, v=0, periodic=True, ne=500)\n", "plt.plot(ens, dos)\n", "plt.xlabel('Energy (t)')\n", "plt.ylabel('DOS')\n", "plt.title('DOS of graphene flake')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Send corrections to bnikolic@udel.edu" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "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.8.3" } }, "nbformat": 4, "nbformat_minor": 4 }