.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "api/gallery/correlations/plot_triangular_graphene_nanoislands.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_api_gallery_correlations_plot_triangular_graphene_nanoislands.py>`
        to download the full example code.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_api_gallery_correlations_plot_triangular_graphene_nanoislands.py:


Magnetism of Triangular Graphene Nanoislands
==================================================

J. Fernández-Rossier and J. J. Palacios (2007) asked which graphene
flakes are magnets, and found the answer in their shape. A triangle with
zigzag edges has all three edges on the same sublattice, :math:`A`, and
more :math:`A` than :math:`B` sites. A bipartite Hamiltonian then has
(at least) :math:`N_Z = |A| - |B|` states at exactly zero energy, all on
sublattice :math:`A` and concentrated on the edges. At half filling they
hold :math:`N_Z` electrons: like the open shell of an atom, they are
filled with parallel spins (Hund's rule), and the flake carries a spin
:math:`S = N_Z/2` -- for *any* repulsion :math:`U`, however small,
while a hexagonal flake needs :math:`U` above a threshold. This is
Lieb's theorem made concrete: the moment grows with the size of the
triangle. Triangular nanographenes (the "triangulenes") were
synthesized and seen to be magnetic in 2017-2021.

:meth:`~tbkit.graphene.GrapheneLattice.triangle_zigzag` cuts the flakes;
:func:`~tbkit.meanfield.hubbard_mean_field` solves the Hubbard model.

.. GENERATED FROM PYTHON SOURCE LINES 22-50

.. code-block:: Python

    import numpy as np
    import matplotlib.pyplot as plt

    from tbkit.graphene import GrapheneLattice
    from tbkit.meanfield import hubbard_mean_field
    from tbkit.plot import Plot
    from tbkit.system import System


    def flake(shape, n):
        lat = GrapheneLattice()
        getattr(lat, shape)(n)
        sys = System(lat)
        sys.set_hopping([{'n': 1, 't': -1.}])
        sys.get_ham()
        return lat, sys


    def ground_state(lat, sys, U, starts=2):
        '''Lowest-energy mean-field solution: a start polarized on sublattice
        A, and random ones.'''
        a = lat.coor['tag'] == 'a'
        up = np.where(a, 0.7, 0.3)
        results = [hubbard_mean_field(sys.ham, U, lat.sites, n_up=up, n_dn=1 - up, max_iter=100000)]
        results += [hubbard_mean_field(sys.ham, U, lat.sites, seed=seed, max_iter=100000) for seed in range(starts)]
        return min(results, key=lambda r: r.energy)









.. GENERATED FROM PYTHON SOURCE LINES 51-56

The two flakes
--------------
A zigzag triangle (left): all three edges end on sublattice A (tag 'a'),
so A has more sites than B. A zigzag hexagon (right) alternates A and B
edges and keeps the two sublattices balanced.

.. GENERATED FROM PYTHON SOURCE LINES 56-61

.. code-block:: Python


    fig_lat, axes_lat = plt.subplots(1, 2, figsize=(10, 4.5))
    Plot(flake('triangle_zigzag', 5)[1]).lattice(plt_hop=True, ms=8, ax=axes_lat[0])
    Plot(flake('hexagon_zigzag', 3)[1]).lattice(plt_hop=True, ms=8, ax=axes_lat[1])




.. image-sg:: /api/gallery/correlations/images/sphx_glr_plot_triangular_graphene_nanoislands_001.png
   :alt: plot triangular graphene nanoislands
   :srcset: /api/gallery/correlations/images/sphx_glr_plot_triangular_graphene_nanoislands_001.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    <Figure size 1000x450 with 2 Axes>



.. GENERATED FROM PYTHON SOURCE LINES 62-64

Zero modes: one per extra A site
---------------------------------------

.. GENERATED FROM PYTHON SOURCE LINES 64-77

.. code-block:: Python


    print('  n  sites  |A|-|B|  zero modes  S (U = 0.1|t|)  S (U = 2|t|)')
    for n in (3, 5, 7, 9):
        lat, sys = flake('triangle_zigzag', n)
        a = lat.coor['tag'] == 'a'
        imbalance = int(a.sum() - (~a).sum())
        energies, states = np.linalg.eigh(sys.ham.toarray())
        zero = np.abs(energies) < 1e-8
        assert zero.sum() == imbalance and np.allclose(states[~a][:, zero], 0.)  # all on A
        spins = [abs(ground_state(lat, sys, U).total_spin) for U in (0.1, 2.)]
        assert np.allclose(spins, imbalance / 2, atol=1e-6)
        print('{:3d} {:6d} {:8d} {:11d} {:15.2f} {:13.2f}'.format(n, lat.sites, imbalance, zero.sum(), *spins))





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

      n  sites  |A|-|B|  zero modes  S (U = 0.1|t|)  S (U = 2|t|)
      3      9        3           3            1.50          1.50
      5     25        5           5            2.50          2.50
      7     49        7           7            3.50          3.50
      9     81        9           9            4.50          4.50




.. GENERATED FROM PYTHON SOURCE LINES 78-84

A moment for any U: the zero modes, half filled
---------------------------------------------------
At small :math:`U` the moment is the density of the zero modes, each
filled by one spin-up electron: :math:`m_i \to \frac12\sum_z|\psi_z(i)|^2`.
A zigzag hexagon of similar size has no zero mode and stays
nonmagnetic at the same :math:`U`.

.. GENERATED FROM PYTHON SOURCE LINES 84-112

.. code-block:: Python


    lat, sys = flake('triangle_zigzag', 9)
    energies, states = np.linalg.eigh(sys.ham.toarray())
    rho = np.sum(np.abs(states[:, np.abs(energies) < 1e-8]) ** 2, axis=1)
    m_small = ground_state(lat, sys, 0.05).magnetization
    deviation = np.max(np.abs(np.abs(m_small) - rho / 2))
    print('Triangle, U = 0.05|t|: |m_i| = zero-mode density / 2 to {:.4f}.'.format(deviation))
    assert deviation < 0.005
    m_large = ground_state(lat, sys, 2.).magnetization

    hexagon, sys_hex = flake('hexagon_zigzag', 4)
    assert not np.any(np.abs(np.linalg.eigvalsh(sys_hex.ham.toarray())) < 1e-8)
    m_hex = ground_state(hexagon, sys_hex, 0.1).magnetization
    print('Hexagon ({} sites), U = 0.1|t|: largest |m_i| = {:.1e}.'.format(hexagon.sites, np.abs(m_hex).max()))
    assert np.abs(m_hex).max() < 1e-6

    fig, axes = plt.subplots(1, 3, figsize=(14, 4.5))
    panels = ((rho / 2, 'zero-mode density / 2'),
                  (m_small * np.sign(m_small.sum()), r'$m_i$, $U = 0.05|t|$'),
                  (m_large * np.sign(m_large.sum()), r'$m_i$, $U = 2|t|$'))
    for ax, (values, title) in zip(axes, panels):
        vmax = np.abs(values).max()
        sc = ax.scatter(lat.coor['x'], lat.coor['y'], c=values, cmap='RdBu_r', vmin=-vmax, vmax=vmax, s=40)
        fig.colorbar(sc, ax=ax, shrink=0.8)
        ax.set_aspect('equal')
        ax.axis('off')
        ax.set_title(title)
    fig.suptitle('Triangular nanoisland, {} sites: S = {:.1f}'.format(lat.sites, abs(m_large.sum())))



.. image-sg:: /api/gallery/correlations/images/sphx_glr_plot_triangular_graphene_nanoislands_002.png
   :alt: Triangular nanoisland, 81 sites: S = 4.5, zero-mode density / 2, $m_i$, $U = 0.05|t|$, $m_i$, $U = 2|t|$
   :srcset: /api/gallery/correlations/images/sphx_glr_plot_triangular_graphene_nanoislands_002.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    Triangle, U = 0.05|t|: |m_i| = zero-mode density / 2 to 0.0022.
    Hexagon (96 sites), U = 0.1|t|: largest |m_i| = 5.7e-10.

    Text(0.5, 0.98, 'Triangular nanoisland, 81 sites: S = 4.5')




.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 7.515 seconds)


.. _sphx_glr_download_api_gallery_correlations_plot_triangular_graphene_nanoislands.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: plot_triangular_graphene_nanoislands.ipynb <plot_triangular_graphene_nanoislands.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: plot_triangular_graphene_nanoislands.py <plot_triangular_graphene_nanoislands.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: plot_triangular_graphene_nanoislands.zip <plot_triangular_graphene_nanoislands.zip>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_