.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "api/gallery/correlations/plot_flat_band_ferromagnetism.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_correlations_plot_flat_band_ferromagnetism.py: Flat-Band Ferromagnetism: Nagaoka, Mielke and Tasaki ========================================================== Itinerant ferromagnetism is hard to get from the Hubbard model: the repulsion :math:`U` favours aligned spins, which cannot share a site, but aligning them costs kinetic energy. Three exact results say when the spins align fully anyway. * **Nagaoka (1966).** At :math:`U = \infty`, with one electron fewer than sites, the ground state is the saturated ferromagnet: the hole moves most freely through aligned spins. On a finite cluster this holds above a finite :math:`U_c`. * **Mielke (1991).** On line graphs such as the kagome lattice, the lowest band is exactly flat. Filled with electrons of one spin, it costs no kinetic energy over any other filling and no interaction energy at all (no site is doubly occupied): the saturated ferromagnet is a ground state for *every* :math:`U > 0`. * **Tasaki (1992, 1995).** Flat bands of other, decorated lattices do the same; and the ferromagnet survives when the flat band is bent into a nearly flat one, of width :math:`W`, provided :math:`U` is large enough: the threshold grows with :math:`W`. The first is checked by exact diagonalization of tiny clusters (a few lines of second quantization below), the other two with :func:`~tbkit.meanfield.hubbard_mean_field` on a :func:`~tbkit.lattices.kagome` torus. .. GENERATED FROM PYTHON SOURCE LINES 29-42 .. code-block:: Python import itertools import numpy as np import matplotlib.pyplot as plt import scipy.sparse as sparse import scipy.sparse.linalg as sla import tbkit.lattices as lattices from tbkit.meanfield import hubbard_mean_field from tbkit.plot import Plot from tbkit.system import System .. GENERATED FROM PYTHON SOURCE LINES 43-52 Nagaoka: one hole in a square cluster ---------------------------------------- With :math:`N_e = N - 1` electrons, the Hubbard Hamiltonian is diagonalized in the sector :math:`S_z = 1/2` (spin-up and spin-down Fock spaces of spinless fermions, coupled by :math:`U`). Its lowest energy equals that of the fully polarized state -- the :math:`N_e` lowest single-particle levels, filled with one spin, which feel no :math:`U` -- exactly when the ground state has the maximal spin :math:`S = N_e/2`. .. GENERATED FROM PYTHON SOURCE LINES 52-91 .. code-block:: Python def fock_hopping(h, n_el): '''Hopping of n_el spinless fermions, and the occupations of the basis.''' basis = list(itertools.combinations(range(len(h)), n_el)) index = {b: k for k, b in enumerate(basis)} occ = np.zeros((len(basis), len(h))) rows, cols, vals = [], [], [] for k, b in enumerate(basis): occ[k, list(b)] = 1 for j in b: for i in np.flatnonzero(h[:, j]): if i not in b: lo, hi = sorted((i, j)) sign = (-1) ** sum(lo < m < hi for m in b) rows.append(index[tuple(sorted(set(b) - {j} | {i}))]) cols.append(k) vals.append(sign * h[i, j]) return sparse.csr_matrix((vals, (rows, cols)), shape=(len(basis),) * 2), occ def hubbard_ground_energy(h, U, n_el): '''Lowest energy of the Hubbard model at S_z = 1/2 (or 0).''' h_up, occ_up = fock_hopping(h, (n_el + 1) // 2) h_dn, occ_dn = fock_hopping(h, n_el // 2) ham = sparse.kron(h_up, sparse.identity(h_dn.shape[0])) \ + sparse.kron(sparse.identity(h_up.shape[0]), h_dn) + U * sparse.diags((occ_up @ occ_dn.T).ravel()) return sla.eigsh(ham, k=1, which='SA')[0][0] def square_cluster(n1, n2): lat = lattices.square() lat.get_lattice(n1=n1, n2=n2) sys = System(lat) sys.set_hopping([{'n': 1, 't': -1.}]) sys.get_ham() return sys.ham.toarray().real .. GENERATED FROM PYTHON SOURCE LINES 92-94 The largest cluster, a :math:`3\times3` square of nine sites with nearest-neighbour hopping and open boundaries: .. GENERATED FROM PYTHON SOURCE LINES 94-101 .. code-block:: Python lat_vis = lattices.square() lat_vis.get_lattice(n1=3, n2=3) vis = System(lat_vis) vis.set_hopping([{'n': 1, 't': -1.}]) fig_sq = Plot(vis).lattice(plt_hop=True, ms=14, figsize=(3.5, 3.5)) .. image-sg:: /api/gallery/correlations/images/sphx_glr_plot_flat_band_ferromagnetism_001.png :alt: plot flat band ferromagnetism :srcset: /api/gallery/correlations/images/sphx_glr_plot_flat_band_ferromagnetism_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 102-103 One hole in each cluster: .. GENERATED FROM PYTHON SOURCE LINES 103-121 .. code-block:: Python nagaoka = {} for n1, n2 in ((2, 2), (3, 2), (3, 3)): h = square_cluster(n1, n2) n_el = len(h) - 1 e_ferro = np.linalg.eigvalsh(h)[:n_el].sum() gaps = [hubbard_ground_energy(h, U, n_el) - e_ferro for U in (5., 1000.)] assert gaps[0] < -1e-3 and abs(gaps[1]) < 1e-8 # low spin at U = 5, S = N_e/2 at U = 1000 if len(h) < 9: # bisection for the threshold U_c lo, hi = 5., 1000. for _ in range(30): mid = np.sqrt(lo * hi) lo, hi = (mid, hi) if hubbard_ground_energy(h, mid, n_el) < e_ferro - 1e-9 else (lo, mid) nagaoka[len(h)] = hi print('Nagaoka threshold: U_c = {:.2f}|t| (2x2), {:.1f}|t| (3x2); 3x3 saturated at U = 1000|t|.' .format(nagaoka[4], nagaoka[6])) assert abs(nagaoka[4] - 18.58) < 0.01 # the known 4-site value .. rst-class:: sphx-glr-script-out .. code-block:: none Nagaoka threshold: U_c = 18.58|t| (2x2), 51.8|t| (3x2); 3x3 saturated at U = 1000|t|. .. GENERATED FROM PYTHON SOURCE LINES 122-127 Mielke: the flat band of the kagome lattice ----------------------------------------------- On a :math:`4\times4` kagome torus (48 sites) with hopping :math:`t > 0`, the flat band at :math:`E = -2t` holds 17 states (16, plus the one where it touches the next band). Seventeen electrons of one spin fill it. .. GENERATED FROM PYTHON SOURCE LINES 127-140 .. code-block:: Python def kagome_torus(t2=0.): lat = lattices.kagome() lat.get_lattice(n1=4, n2=4) sys = System(lat, periodic=True) sys.set_hopping([{'n': 1, 't': 1.}]) if t2: sys.set_hopping([{'n': 2, 't': t2}]) sys.get_ham() return sys .. GENERATED FROM PYTHON SOURCE LINES 141-144 The kagome lattice: three sites per unit cell ('a', 'b' and 'c', one per colour) on corner-sharing triangles. The bonds that close the torus are not drawn. .. GENERATED FROM PYTHON SOURCE LINES 144-147 .. code-block:: Python fig_kag = Plot(kagome_torus()).lattice(plt_hop=True, ms=10, figsize=(5.5, 4.5)) .. image-sg:: /api/gallery/correlations/images/sphx_glr_plot_flat_band_ferromagnetism_002.png :alt: plot flat band ferromagnetism :srcset: /api/gallery/correlations/images/sphx_glr_plot_flat_band_ferromagnetism_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 148-149 Mean field from a saturated and from a paramagnetic start: .. GENERATED FROM PYTHON SOURCE LINES 149-175 .. code-block:: Python def polarized_and_paramagnetic(sys, U, n_el): '''Mean-field solutions from a saturated and from a paramagnetic start.''' n = sys.lat.sites _, vec = np.linalg.eigh(sys.ham.toarray()) up = np.sum(np.abs(vec[:, :n_el]) ** 2, axis=1) ferro = hubbard_mean_field(sys.ham, U, n_el, n_up=up, n_dn=np.zeros(n)) para = hubbard_mean_field(sys.ham, U, n_el, n_up=np.full(n, n_el / (2 * n)), n_dn=np.full(n, n_el / (2 * n))) return ferro, para sys = kagome_torus() energies = np.linalg.eigvalsh(sys.ham.toarray()) n_flat = int(np.sum(np.isclose(energies, -2.))) assert n_flat == 17 and np.all(energies[n_flat:] > -2. + 1e-6) for U in (0.01, 0.5, 2., 8.): ferro, para = polarized_and_paramagnetic(sys, U, n_flat) # saturated, with the noninteracting energy -2t per electron: a ground state assert np.isclose(ferro.total_spin, n_flat / 2) and np.isclose(ferro.energy, -2. * n_flat) assert para.energy > ferro.energy and np.isclose(para.total_spin, 0.) # random starts reach other, higher, self-consistent solutions for seed in range(2): assert hubbard_mean_field(sys.ham, U, n_flat, seed=seed).energy > ferro.energy - 1e-9 print('Kagome flat band, 17 electrons: S = 17/2 and E = -34t for every U tried, down to 0.01t.') .. rst-class:: sphx-glr-script-out .. code-block:: none Kagome flat band, 17 electrons: S = 17/2 and E = -34t for every U tried, down to 0.01t. .. GENERATED FROM PYTHON SOURCE LINES 176-182 Tasaki: a nearly flat band needs a large enough U ------------------------------------------------------ A second-neighbour hopping :math:`t_2` bends the flat band, to a width :math:`W`. The saturated ferromagnet is then a mean-field solution, and beats the paramagnet, only above a threshold :math:`U_c`, here proportional to the width. .. GENERATED FROM PYTHON SOURCE LINES 182-218 .. code-block:: Python def saturated_wins(sys, U, n_el): ferro, para = polarized_and_paramagnetic(sys, U, n_el) return np.isclose(ferro.total_spin, n_el / 2) and ferro.energy < para.energy widths, thresholds = [], [] for t2 in (0., -0.02, -0.05, -0.1, -0.15): sys = kagome_torus(t2) energies = np.linalg.eigvalsh(sys.ham.toarray()) widths.append(energies[n_flat - 1] - energies[0]) lo, hi = 0., 10. assert saturated_wins(sys, hi, n_flat) for _ in range(25): mid = (lo + hi) / 2 lo, hi = (lo, mid) if saturated_wins(sys, mid, n_flat) else (mid, hi) thresholds.append(hi) widths, thresholds = np.array(widths), np.array(thresholds) for w, u in zip(widths, thresholds): print('W = {:.3f}t: U_c = {:.3f}t'.format(w, u)) ratio = thresholds[1:] / widths[1:] print('U_c / W = {:.2f}, constant to {:.1%}.'.format(ratio.mean(), np.ptp(ratio) / ratio.mean())) assert thresholds[0] < 1e-6 and np.allclose(ratio, ratio.mean(), rtol=0.01) fig, axes = plt.subplots(1, 2, figsize=(10, 4)) axes[0].plot([4, 6], [nagaoka[4], nagaoka[6]], 'o-r') axes[0].set_xlabel('sites N (one hole: N - 1 electrons)') axes[0].set_ylabel(r'$U_c/|t|$') axes[0].set_xticks([4, 6]) axes[0].set_title('Nagaoka: saturated ferromagnet above $U_c$') axes[1].plot(widths, thresholds, 'o-b') axes[1].set_xlabel(r'width $W$ of the lowest 17 states ($t$)') axes[1].set_ylabel(r'$U_c/t$') axes[1].set_title(r'Mielke ($W = 0$) and Tasaki: kagome, $4\times4$ torus') fig.set_layout_engine('tight') .. image-sg:: /api/gallery/correlations/images/sphx_glr_plot_flat_band_ferromagnetism_003.png :alt: Nagaoka: saturated ferromagnet above $U_c$, Mielke ($W = 0$) and Tasaki: kagome, $4\times4$ torus :srcset: /api/gallery/correlations/images/sphx_glr_plot_flat_band_ferromagnetism_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none W = 0.000t: U_c = 0.000t W = 0.100t: U_c = 0.346t W = 0.251t: U_c = 0.866t W = 0.502t: U_c = 1.736t W = 0.755t: U_c = 2.610t U_c / W = 3.46, constant to 0.1%. .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 2.012 seconds) .. _sphx_glr_download_api_gallery_correlations_plot_flat_band_ferromagnetism.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_flat_band_ferromagnetism.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_flat_band_ferromagnetism.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_flat_band_ferromagnetism.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_