.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "api/gallery/disorder/plot_lifshitz_tails.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_api_gallery_disorder_plot_lifshitz_tails.py: Lifshitz Tails in the Density of States of a Disordered Lattice ===================================================================== Disorder moves states beyond the band edge of the clean crystal, but very few. Lifshitz (1964) explained why: a state deep in the tail needs a large region of size :math:`\ell` where the disorder happens to conspire (every site low, or every site present), which costs a probability exponentially small in the volume :math:`\ell^d`, while the kinetic energy of a state confined to that region only falls as :math:`1/\ell^2`. Balancing the two gives the integrated density of states near the lowest possible energy :math:`E_0`, .. math:: N(E) \sim \exp\!\left[-c\,(E - E_0)^{-d/2}\right]\, , an essential singularity, smaller than any power of :math:`E - E_0`. In 1D, with sites missing at random, the law is exact: a chain where each site is present with probability :math:`p` splits into segments, a segment of :math:`\ell` sites occurs with weight :math:`p^\ell`, and its lowest level is :math:`E_\ell = -2t\cos[\pi/(\ell+1)] \approx -2t + t\pi^2/\ell^2`, so that :math:`\ln N(E) \approx -\pi\ln(1/p)\,\sqrt{t/(E + 2t)}`. In 2D the tail is just as real but much harder to reach: the asymptotic law sets in only where :math:`N(E)` is far below one state per sample. What a finite lattice, with any method, can show is its hallmark, :math:`N(E)` falling faster than every power of :math:`E - E_0`, and where the sample's own tail stops. .. GENERATED FROM PYTHON SOURCE LINES 33-44 .. code-block:: Python import numpy as np import matplotlib.pyplot as plt from scipy.integrate import cumulative_trapezoid from scipy.linalg import eigvalsh_tridiagonal import tbkit.kpm as kpm import tbkit.lattices as lattices from tbkit.plot import Plot from tbkit.system import System .. GENERATED FROM PYTHON SOURCE LINES 45-52 The diluted chain: the exact Lifshitz law --------------------------------------------- :meth:`~tbkit.lattice.Lattice.remove_sites` removes 40% of the sites of a 200 000-site chain at random. The Hamiltonian is tridiagonal, with a zero bond wherever a site is missing, so its full spectrum is cheap. The integrated density of states just above :math:`E_\ell` falls by :math:`p` per added site, :math:`\ln N(E_\ell) \simeq \ell\ln p`. .. GENERATED FROM PYTHON SOURCE LINES 52-91 .. code-block:: Python p, n = 0.6, 200000 lat = lattices.chain() lat.get_lattice(n1=n) rng = np.random.default_rng(0) lat.remove_sites(np.flatnonzero(rng.random(n) > p).tolist()) sys = System(lat) sys.set_hopping([{'n': 1, 't': -1.}]) sys.get_ham() en = eigvalsh_tridiagonal(sys.ham.diagonal().real, sys.ham.diagonal(1).real) assert en.min() > -2. # no state below the clean band edge ... ells = np.arange(4, 15) e_ell = -2. * np.cos(np.pi / (ells + 1)) n_int = np.array([np.mean(en < e + 1e-9) for e in e_ell]) # N(E) just above E_ell slope = np.polyfit(ells, np.log(n_int), 1)[0] print('ln N(E_l) = {:.3f} l + const (expected ln p = {:.3f})'.format(slope, np.log(p))) assert abs(slope / np.log(p) - 1.) < 0.05 # ... but an exponentially thin tail just above it def square(n, p=1., w=0.): '''An n x n square torus, with a fraction 1 - p of sites removed or onsite disorder w.''' lat = lattices.square() lat.get_lattice(n1=n, n2=n) if p < 1.: rng = np.random.default_rng(0) lat.remove_sites(np.flatnonzero(rng.random(lat.sites) > p).tolist()) sys = System(lat, periodic=True) sys.set_hopping([{'n': 1, 't': -1.}]) sys.set_onsite({'a': 0.}) if w: state = np.random.get_state() np.random.seed(0) sys.set_onsite_dis(alpha=w) np.random.set_state(state) sys.get_ham() return sys .. rst-class:: sphx-glr-script-out .. code-block:: none ln N(E_l) = -0.521 l + const (expected ln p = -0.511) .. GENERATED FROM PYTHON SOURCE LINES 92-95 The diluted lattices, drawn: 30 sites of the chain with 40% removed (top), and a :math:`16\times16` square with 5% removed (bottom). A missing site cuts every bond it carried. .. GENERATED FROM PYTHON SOURCE LINES 95-112 .. code-block:: Python lat_vis = lattices.chain() lat_vis.get_lattice(n1=30) lat_vis.remove_sites(np.flatnonzero(np.random.default_rng(0).random(30) > p).tolist()) vis = System(lat_vis) vis.set_hopping([{'n': 1, 't': -1.}]) fig_lat, axes_lat = plt.subplots(2, 1, figsize=(7, 7.5), height_ratios=(1, 6)) Plot(vis).lattice(plt_hop=True, ms=8, ax=axes_lat[0]) Plot(square(16, p=0.95)).lattice(plt_hop=True, ms=8, ax=axes_lat[1]) def kpm_dos(sys, e, bounds): '''Density of states per site, and integrated density of states in states.''' rho = kpm.dos(sys.ham, n_moments=1024, n_random=8, e_grid=e, seed=1, bounds=bounds)[1] return rho / sys.lat.sites, cumulative_trapezoid(rho, e, initial=0.) .. image-sg:: /api/gallery/disorder/images/sphx_glr_plot_lifshitz_tails_001.png :alt: plot lifshitz tails :srcset: /api/gallery/disorder/images/sphx_glr_plot_lifshitz_tails_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 113-129 The diluted square lattice: faster than any power ----------------------------------------------------- With 5% of the sites of a :math:`400 \times 400` square lattice removed, the lowest possible energy is still the clean edge :math:`E_0 = -4t`, reached only inside ever larger vacancy-free disks. The kernel polynomial method (:func:`~tbkit.kpm.dos`) gives :math:`N(E)`. Near a clean 2D band edge the density of states is a constant :math:`1/4\pi t`, so :math:`N \propto E - E_0` and the local exponent :math:`d\ln N/d\ln(E - E_0)` is 1. In the diluted lattice it keeps growing towards :math:`E_0`. The asymptotic constant itself, :math:`(E - E_0)\ln N \to -\pi j_0^2\ln(1/p)` (:math:`j_0` the first zero of the Bessel function :math:`J_0`, from the ground state of a disk), takes over only for :math:`E - E_0 \lesssim 0.01t`, where :math:`N \sim e^{-90}`: no sample holds such a state, which is why the quantitative test is the 1D one. .. GENERATED FROM PYTHON SOURCE LINES 129-144 .. code-block:: Python xs = np.array([0.05, 0.1, 0.2, 0.4, 0.8]) # E - E_0 e_edge = np.linspace(-4., -3., 401) clean, diluted = square(200), square(400, p=0.95) n_clean = kpm_dos(clean, e_edge, (-4.01, 4.01))[1] / clean.lat.sites n_dil = kpm_dos(diluted, e_edge, (-4.01, 4.01))[1] / diluted.lat.sites exp_clean = np.diff(np.log(np.interp(xs - 4., e_edge, n_clean))) / np.diff(np.log(xs)) exp_dil = np.diff(np.log(np.interp(xs - 4., e_edge, n_dil))) / np.diff(np.log(xs)) print('local exponents d ln N / d ln(E - E_0), from E - E_0 = 0.8t down to 0.05t:') print(' clean {}'.format(np.round(exp_clean[::-1], 2))) print(' diluted {}'.format(np.round(exp_dil[::-1], 2))) assert np.allclose(exp_clean, 1., atol=0.1) # N = (E - E_0)/4 pi assert np.all(np.diff(exp_dil) < 0) and exp_dil[0] > 5. # steeper and steeper towards E_0 print('clean edge: N ~ (E - E_0)^1; diluted tail: local exponent rises past 5: OK') .. rst-class:: sphx-glr-script-out .. code-block:: none local exponents d ln N / d ln(E - E_0), from E - E_0 = 0.8t down to 0.05t: clean [1.04 1.03 1.03 1.02] diluted [1.27 1.51 2.38 7.23] clean edge: N ~ (E - E_0)^1; diluted tail: local exponent rises past 5: OK .. GENERATED FROM PYTHON SOURCE LINES 145-156 Uniform disorder: where a finite sample's tail stops -------------------------------------------------------- With onsite energies uniform in :math:`[-W, W]` (:meth:`~tbkit.system.System.set_onsite_dis`), :math:`E_0 = -4t - W`. The density of states of a :math:`200 \times 200` sample spills past the clean edge :math:`-4t`, then collapses at the sample's lowest eigenvalue (:meth:`~tbkit.system.System.get_eig_sparse`), more than :math:`t` above :math:`E_0`. The kernel polynomial method finds no state below it, as it should: the cutoff is the sample's, not the method's. Further down the integrated density of states is less than one state per 40 000 sites, so the rest of the tail exists only in larger samples. .. GENERATED FROM PYTHON SOURCE LINES 156-171 .. code-block:: Python W = 2. dis = square(200, w=W) dis.get_eig_sparse(n_eig=1, sigma=-4. - W) e_min = dis.en[0].real e_dis = np.linspace(-6., -3., 241) bounds = (-4. - W - 0.01, 4. + W + 0.01) rho_dis, count_dis = kpm_dos(dis, e_dis, bounds) rho_clean = kpm_dos(clean, e_dis, bounds)[0] assert e_min < -4.5 and e_min - (-4. - W) > 1. # past the clean edge, far above E_0 assert np.interp(e_min - 0.05, e_dis, count_dis) < 1. # no state below the lowest eigenvalue ... assert np.interp(e_min + 0.2, e_dis, count_dis) > 10. # ... and the sample's tail just above it print('W = {:g}: lowest eigenvalue {:.3f}, {:.2f}t above E_0; KPM count below it < 1 state: OK' .format(W, e_min, e_min + 4. + W)) .. rst-class:: sphx-glr-script-out .. code-block:: none W = 2: lowest eigenvalue -4.852, 1.15t above E_0; KPM count below it < 1 state: OK .. GENERATED FROM PYTHON SOURCE LINES 172-174 The three tails ------------------- .. GENERATED FROM PYTHON SOURCE LINES 174-206 .. code-block:: Python fig, axes = plt.subplots(1, 3, figsize=(15, 4.2)) e_fine = np.linspace(-1.97, -1.6, 400) axes[0].semilogy(1 / np.sqrt(e_fine + 2.), [np.mean(en < e) for e in e_fine], 'C0', label='diluted chain') axes[0].semilogy(1 / np.sqrt(e_ell + 2.), np.exp(np.polyval(np.polyfit(ells, np.log(n_int), 1), ells)), 'k--', label=r'fit through the $E_\ell$') axes[0].set_xlim(1.6, 6.) axes[0].set_xlabel(r'$(E + 2t)^{-1/2}$') axes[0].set_ylabel('integrated DOS $N(E)$') axes[0].set_title('1D, $p = {}$: $\\ln N$ linear in $(E-E_0)^{{-1/2}}$'.format(p)) axes[0].legend() x = e_edge[1:] + 4. axes[1].loglog(x, n_clean[1:], 'C7', label='clean') axes[1].loglog(x, n_dil[1:], 'C0', label='5% vacancies') axes[1].loglog(x, x / (4 * np.pi), 'k:', label=r'$(E-E_0)/4\pi$') axes[1].set_xlim(0.03, 1.) axes[1].set_ylim(1e-6, 0.2) axes[1].set_xlabel('$(E - E_0)/t$') axes[1].set_ylabel('integrated DOS $N(E)$') axes[1].set_title('2D diluted: faster than any power') axes[1].legend() axes[2].semilogy(e_dis, np.maximum(rho_clean, 1e-9), 'C7', label='clean') axes[2].semilogy(e_dis, np.maximum(rho_dis, 1e-9), 'C3', label='$W = {:g}t$'.format(W)) axes[2].axvline(e_min, c='C3', ls='--', label='lowest eigenvalue') for x0 in (-4., -4. - W): axes[2].axvline(x0, c='k', ls=':') axes[2].set_ylim(1e-9, 1) axes[2].set_xlabel('$E/t$') axes[2].set_ylabel(r'$\rho(E)$ per site') axes[2].set_title('2D uniform disorder: the sample ends before $E_0$') axes[2].legend() fig.tight_layout() .. image-sg:: /api/gallery/disorder/images/sphx_glr_plot_lifshitz_tails_002.png :alt: 1D, $p = 0.6$: $\ln N$ linear in $(E-E_0)^{-1/2}$, 2D diluted: faster than any power, 2D uniform disorder: the sample ends before $E_0$ :srcset: /api/gallery/disorder/images/sphx_glr_plot_lifshitz_tails_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 11.068 seconds) .. _sphx_glr_download_api_gallery_disorder_plot_lifshitz_tails.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_lifshitz_tails.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_lifshitz_tails.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_lifshitz_tails.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_