Haldane Model¶
The Haldane model, a time-reversal symmetry breaking model of Graphene, is the fruitfly of (band) topological systems. Let us set up a simple tight-binding model and investigate its topological properties.
We use a pymatgen Structure as as input. We start by constructing a hexagonal unit cell. Pymatgen structures are always three-dimensional, so we just put some vacuum along the c-direction of the unit cell. We also set the space group symmetry to P6/mmm (#191).
In [1]: from pymatgen.core.structure import Lattice, Structure
In [2]: space_group = 191
In [3]: lattice = Lattice.hexagonal(1.42, 10)
In [4]: structure = Structure.from_spacegroup(sg=space_group, lattice=lattice, species=['C'], coords=[[1 / 3, 2 / 3, 0]])
In [5]: print(structure)
Full Formula (C2)
Reduced Formula: C
abc : 1.420000 1.420000 10.000000
angles: 90.000000 90.000000 120.000000
pbc : True True True
Sites (2)
# SP a b c
--- ---- -------- -------- ---
0 C 0.333333 0.666667 0
1 C 0.666667 0.333333 0
Next, we instantiate a tight-binding model using that structure.
In [6]: import topwave as tp
In [7]: model = tp.TightBindingModel(structure)
Graphene has hybrid orbitals sp2 on its sites with three major inplane lobes at 120 degree and the pz orbital sticking out of plane. The inplane orbitals form strong sigma bonds, whereas the pz orbitals form pi bonds. These bonds are responsible for graphene’s low energy properties, so we only need to consider a single orbital per site. CITE?
We can confirm that a single orbital per site is topwave’s default setting by printing the site properties.
In [8]: model.show_site_properties()
╒═════════╤═══════════╤═════════╤════════════╤════════════════════════════════════╤═══════════════════════════════════════════════════╤══════════╤═════════════════╤═════════════════╤═════════════════════════════════════════════════════╤═══════════════════╤════════════════════╤═════════╕
│ index │ species │ label │ orbitals │ coordinates (latt.) │ coordinates (cart.) │ magmom │ onsite scalar │ onsite vector │ onsite matrix │ unit cell index │ supercell vector │ layer │
╞═════════╪═══════════╪═════════╪════════════╪════════════════════════════════════╪═══════════════════════════════════════════════════╪══════════╪═════════════════╪═════════════════╪═════════════════════════════════════════════════════╪═══════════════════╪════════════════════╪═════════╡
│ 0 │ C1 │ │ 1 │ [0.33333333 0.66666667 0. ] │ [-7.09999093e-11 8.19837382e-01 8.69499227e-17] │ │ 0 │ [0. 0. 0.] │ [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]] │ │ │ │
├─────────┼───────────┼─────────┼────────────┼────────────────────────────────────┼───────────────────────────────────────────────────┼──────────┼─────────────────┼─────────────────┼─────────────────────────────────────────────────────┼───────────────────┼────────────────────┼─────────┤
│ 1 │ C1 │ │ 1 │ [0.66666667 0.33333333 0. ] │ [7.10000000e-01 4.09918691e-01 8.69499227e-17] │ │ 0 │ [0. 0. 0.] │ [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]] │ │ │ │
╘═════════╧═══════════╧═════════╧════════════╧════════════════════════════════════╧═══════════════════════════════════════════════════╧══════════╧═════════════════╧═════════════════╧═════════════════════════════════════════════════════╧═══════════════════╧════════════════════╧═════════╛
Spin Polarized: False
Zeeman: [0. 0. 0.]
Supercell Size: None
Great, the orbitals column shows one orbital for each site. Next, we generate the possible couplings up to next-nearest neighbors and output them.
In [9]: model.generate_couplings(max_distance=1.5, space_group=space_group)
In [10]: model.show_couplings()
╒═════════╤══════════════════╤══════════════════════╤════════════╤══════════════════╤═══════════════════════════════════════╤═════════╤════════════╤═════════╤════════════╤════════════╤═════════════════════╤═════════════════════════════════════════════════════╕
│ index │ symmetry index │ symmetry operation │ distance │ lattice_vector │ sublattice_vector │ site1 │ orbital1 │ site2 │ orbital2 │ strength │ spin-orbit vector │ matrix │
╞═════════╪══════════════════╪══════════════════════╪════════════╪══════════════════╪═══════════════════════════════════════╪═════════╪════════════╪═════════╪════════════╪════════════╪═════════════════════╪═════════════════════════════════════════════════════╡
│ 0 │ 0 │ x, y, z │ 0.819837 │ [-1. 0. 0.] │ [ 0.33333333 -0.33333333 0. ] │ 0 │ 0 │ 1 │ 0 │ 0j │ [0. 0. 0.] │ [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]] │
├─────────┼──────────────────┼──────────────────────┼────────────┼──────────────────┼───────────────────────────────────────┼─────────┼────────────┼─────────┼────────────┼────────────┼─────────────────────┼─────────────────────────────────────────────────────┤
│ 1 │ 0 │ x-y, x, z │ 0.819837 │ [-0. -1. -0.] │ [-0.33333333 0.33333333 0. ] │ 1 │ 0 │ 0 │ 0 │ 0j │ [0. 0. 0.] │ [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]] │
├─────────┼──────────────────┼──────────────────────┼────────────┼──────────────────┼───────────────────────────────────────┼─────────┼────────────┼─────────┼────────────┼────────────┼─────────────────────┼─────────────────────────────────────────────────────┤
│ 2 │ 0 │ x-y, -y, z │ 0.819837 │ [0. 0. 0.] │ [ 0.33333333 -0.33333333 0. ] │ 0 │ 0 │ 1 │ 0 │ 0j │ [0. 0. 0.] │ [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]] │
├─────────┼──────────────────┼──────────────────────┼────────────┼──────────────────┼───────────────────────────────────────┼─────────┼────────────┼─────────┼────────────┼────────────┼─────────────────────┼─────────────────────────────────────────────────────┤
│ 3 │ 1 │ x, y, z │ 1.42 │ [0. 1. 0.] │ [0. 0. 0.] │ 0 │ 0 │ 0 │ 0 │ 0j │ [0. 0. 0.] │ [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]] │
├─────────┼──────────────────┼──────────────────────┼────────────┼──────────────────┼───────────────────────────────────────┼─────────┼────────────┼─────────┼────────────┼────────────┼─────────────────────┼─────────────────────────────────────────────────────┤
│ 4 │ 1 │ -y, -x, z │ 1.42 │ [-1. 0. 0.] │ [0. 0. 0.] │ 0 │ 0 │ 0 │ 0 │ 0j │ [0. 0. 0.] │ [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]] │
├─────────┼──────────────────┼──────────────────────┼────────────┼──────────────────┼───────────────────────────────────────┼─────────┼────────────┼─────────┼────────────┼────────────┼─────────────────────┼─────────────────────────────────────────────────────┤
│ 5 │ 1 │ -x+y, y, z │ 1.42 │ [ 1. 1. -0.] │ [0. 0. 0.] │ 0 │ 0 │ 0 │ 0 │ 0j │ [0. 0. 0.] │ [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]] │
├─────────┼──────────────────┼──────────────────────┼────────────┼──────────────────┼───────────────────────────────────────┼─────────┼────────────┼─────────┼────────────┼────────────┼─────────────────────┼─────────────────────────────────────────────────────┤
│ 6 │ 1 │ y, x, z │ 1.42 │ [ 1. -0. -0.] │ [0. 0. 0.] │ 1 │ 0 │ 1 │ 0 │ 0j │ [0. 0. 0.] │ [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]] │
├─────────┼──────────────────┼──────────────────────┼────────────┼──────────────────┼───────────────────────────────────────┼─────────┼────────────┼─────────┼────────────┼────────────┼─────────────────────┼─────────────────────────────────────────────────────┤
│ 7 │ 1 │ -x, -x+y, z │ 1.42 │ [-0. 1. -0.] │ [0. 0. 0.] │ 1 │ 0 │ 1 │ 0 │ 0j │ [0. 0. 0.] │ [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]] │
├─────────┼──────────────────┼──────────────────────┼────────────┼──────────────────┼───────────────────────────────────────┼─────────┼────────────┼─────────┼────────────┼────────────┼─────────────────────┼─────────────────────────────────────────────────────┤
│ 8 │ 1 │ x-y, -y, z │ 1.42 │ [-1. -1. 0.] │ [0. 0. 0.] │ 1 │ 0 │ 1 │ 0 │ 0j │ [0. 0. 0.] │ [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]] │
╘═════════╧══════════════════╧══════════════════════╧════════════╧══════════════════╧═══════════════════════════════════════╧═════════╧════════════╧═════════╧════════════╧════════════╧═════════════════════╧═════════════════════════════════════════════════════╛
Onsite (Inter-Orbital) Couplings:
╒═════════╤══════════════════╤══════════════════════╤════════════╤══════════════════╤═════════════════════╤═════════╤════════════╤═════════╤════════════╤════════════╤═════════════════════╤══════════╕
│ index │ symmetry index │ symmetry operation │ distance │ lattice_vector │ sublattice_vector │ site1 │ orbital1 │ site2 │ orbital2 │ strength │ spin-orbit vector │ matrix │
╞═════════╪══════════════════╪══════════════════════╪════════════╪══════════════════╪═════════════════════╪═════════╪════════════╪═════════╪════════════╪════════════╪═════════════════════╪══════════╡
╘═════════╧══════════════════╧══════════════════════╧════════════╧══════════════════╧═════════════════════╧═════════╧════════════╧═════════╧════════════╧════════════╧═════════════════════╧══════════╛
There are three nearest neighbors bonds, about 0.82 Angstrom apart. We now couple these sites by some hopping amplitude that accounts for the mobile electrons along the pi bonds. We can select the bonds from the list in different ways. Because we input the space group symmetry, topwave automatically classifies the bonds based on symmetry. Graphene has rotational symmetry and all the nearest neighbors are symmetrically equivalent, which is why they share the same symmetry index. Let’s assign the hopping using that symmetry index.
In [11]: t = -1
In [12]: model.set_coupling(attribute_value=0, strength=t, attribute='symmetry_id')
In [13]: model.show_couplings()
╒═════════╤══════════════════╤══════════════════════╤════════════╤══════════════════╤═══════════════════════════════════════╤═════════╤════════════╤═════════╤════════════╤════════════╤═════════════════════╤═════════════════════════════════════════════════════╕
│ index │ symmetry index │ symmetry operation │ distance │ lattice_vector │ sublattice_vector │ site1 │ orbital1 │ site2 │ orbital2 │ strength │ spin-orbit vector │ matrix │
╞═════════╪══════════════════╪══════════════════════╪════════════╪══════════════════╪═══════════════════════════════════════╪═════════╪════════════╪═════════╪════════════╪════════════╪═════════════════════╪═════════════════════════════════════════════════════╡
│ 0 │ 0 │ x, y, z │ 0.819837 │ [-1. 0. 0.] │ [ 0.33333333 -0.33333333 0. ] │ 0 │ 0 │ 1 │ 0 │ -1 │ [0. 0. 0.] │ [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]] │
├─────────┼──────────────────┼──────────────────────┼────────────┼──────────────────┼───────────────────────────────────────┼─────────┼────────────┼─────────┼────────────┼────────────┼─────────────────────┼─────────────────────────────────────────────────────┤
│ 1 │ 0 │ x-y, x, z │ 0.819837 │ [-0. -1. -0.] │ [-0.33333333 0.33333333 0. ] │ 1 │ 0 │ 0 │ 0 │ -1 │ [0. 0. 0.] │ [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]] │
├─────────┼──────────────────┼──────────────────────┼────────────┼──────────────────┼───────────────────────────────────────┼─────────┼────────────┼─────────┼────────────┼────────────┼─────────────────────┼─────────────────────────────────────────────────────┤
│ 2 │ 0 │ x-y, -y, z │ 0.819837 │ [0. 0. 0.] │ [ 0.33333333 -0.33333333 0. ] │ 0 │ 0 │ 1 │ 0 │ -1 │ [0. 0. 0.] │ [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]] │
├─────────┼──────────────────┼──────────────────────┼────────────┼──────────────────┼───────────────────────────────────────┼─────────┼────────────┼─────────┼────────────┼────────────┼─────────────────────┼─────────────────────────────────────────────────────┤
│ 3 │ 1 │ x, y, z │ 1.42 │ [0. 1. 0.] │ [0. 0. 0.] │ 0 │ 0 │ 0 │ 0 │ 0j │ [0. 0. 0.] │ [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]] │
├─────────┼──────────────────┼──────────────────────┼────────────┼──────────────────┼───────────────────────────────────────┼─────────┼────────────┼─────────┼────────────┼────────────┼─────────────────────┼─────────────────────────────────────────────────────┤
│ 4 │ 1 │ -y, -x, z │ 1.42 │ [-1. 0. 0.] │ [0. 0. 0.] │ 0 │ 0 │ 0 │ 0 │ 0j │ [0. 0. 0.] │ [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]] │
├─────────┼──────────────────┼──────────────────────┼────────────┼──────────────────┼───────────────────────────────────────┼─────────┼────────────┼─────────┼────────────┼────────────┼─────────────────────┼─────────────────────────────────────────────────────┤
│ 5 │ 1 │ -x+y, y, z │ 1.42 │ [ 1. 1. -0.] │ [0. 0. 0.] │ 0 │ 0 │ 0 │ 0 │ 0j │ [0. 0. 0.] │ [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]] │
├─────────┼──────────────────┼──────────────────────┼────────────┼──────────────────┼───────────────────────────────────────┼─────────┼────────────┼─────────┼────────────┼────────────┼─────────────────────┼─────────────────────────────────────────────────────┤
│ 6 │ 1 │ y, x, z │ 1.42 │ [ 1. -0. -0.] │ [0. 0. 0.] │ 1 │ 0 │ 1 │ 0 │ 0j │ [0. 0. 0.] │ [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]] │
├─────────┼──────────────────┼──────────────────────┼────────────┼──────────────────┼───────────────────────────────────────┼─────────┼────────────┼─────────┼────────────┼────────────┼─────────────────────┼─────────────────────────────────────────────────────┤
│ 7 │ 1 │ -x, -x+y, z │ 1.42 │ [-0. 1. -0.] │ [0. 0. 0.] │ 1 │ 0 │ 1 │ 0 │ 0j │ [0. 0. 0.] │ [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]] │
├─────────┼──────────────────┼──────────────────────┼────────────┼──────────────────┼───────────────────────────────────────┼─────────┼────────────┼─────────┼────────────┼────────────┼─────────────────────┼─────────────────────────────────────────────────────┤
│ 8 │ 1 │ x-y, -y, z │ 1.42 │ [-1. -1. 0.] │ [0. 0. 0.] │ 1 │ 0 │ 1 │ 0 │ 0j │ [0. 0. 0.] │ [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]] │
╘═════════╧══════════════════╧══════════════════════╧════════════╧══════════════════╧═══════════════════════════════════════╧═════════╧════════════╧═════════╧════════════╧════════════╧═════════════════════╧═════════════════════════════════════════════════════╛
Onsite (Inter-Orbital) Couplings:
╒═════════╤══════════════════╤══════════════════════╤════════════╤══════════════════╤═════════════════════╤═════════╤════════════╤═════════╤════════════╤════════════╤═════════════════════╤══════════╕
│ index │ symmetry index │ symmetry operation │ distance │ lattice_vector │ sublattice_vector │ site1 │ orbital1 │ site2 │ orbital2 │ strength │ spin-orbit vector │ matrix │
╞═════════╪══════════════════╪══════════════════════╪════════════╪══════════════════╪═════════════════════╪═════════╪════════════╪═════════╪════════════╪════════════╪═════════════════════╪══════════╡
╘═════════╧══════════════════╧══════════════════════╧════════════╧══════════════════╧═════════════════════╧═════════╧════════════╧═════════╧════════════╧════════════╧═════════════════════╧══════════╛
We can see in the table that first three bonds now have strength -1. Next, we want to compute the
band structure of our model. As an input we just need a list of k-points. We
could create the list manually, or use the functionality of topwave.set_of_kpoints.Path.
In [14]: path = tp.Path([[0, 0, 0],
....: [1 / 2, 0, 0],
....: [1 / 3, 1 / 3, 0],
....: [0, 0, 0]])
....:
In [15]: spec = tp.Spec(model, path)
# Let's plot it
In [16]: import matplotlib.pyplot as plt
In [17]: fig, ax = plt.subplots()
In [18]: for band in spec.energies.T:
....: ax.plot(band, ls='-', lw=2.5, c='deeppink')
....:
In [19]: ax.set_xticks(path.node_indices);
In [20]: ax.set_ylabel('Energy');
In [21]: ax.set_xticklabels([r'$\Gamma$', 'M', 'K', r'$\Gamma$']);
We can see graphene’s famous Dirac point at the K-point. It’s time reversal partner, the K’-point, sits at opposite momentum. Let us construct the Wilson loop around the K-point using class:topwave.set_of_kpoints.Circle and calculate the Berry phase for the valence band. We also add a small staggered onsite potential to the sublattices. This so-called mass gap is trivial, and makes the Berry curvature at the K- and K’-point well defined.
In [22]: mass_term = 0.005
In [23]: model.set_onsite_scalar(0, mass_term)
In [24]: model.set_onsite_scalar(1, -mass_term)
In [25]: circle = tp.Circle(radius=0.02, center=[1 / 3, 1 / 3, 0], normal=[0, 0, 1])
In [26]: spec = tp.Spec(model, circle)
In [27]: print(spec.get_berry_phase(band_indices=[0]))
2.9995524533587985
As expected, the Berry phase is quantized to (-)pi. Because the system is time-reversal invariant, the Berry phase around the K’-point must have the opposite sign. There is no net flux of Berry curvature. Let us have a look at the distribution of Berry flux over the entire two-dimensional Brillouin zone. We cover the plane with small plaquettes and compute the flux through each plaquette.
# Create a cover of plaquettes of the xy-plane (with the z-direction as the norm).
In [28]: num_x, num_y = 13, 13
In [29]: plaquettes = tp.get_plaquette_cover('z', num_x, num_y)
# Calculate the spectrum for each plaquette.
In [30]: spectra = [tp.Spec(model, plaquette) for plaquette in plaquettes]
# Compute the Berry phase of each spectrum and plot it.
In [31]: fluxes = np.array([spectrum.get_berry_phase(band_indices=[0]) for spectrum in spectra])
In [32]: fig, ax = plt.subplots()
In [33]: im = ax.imshow(fluxes.reshape((num_x, num_y)), origin='lower', cmap='PiYG', vmin=-np.pi, vmax=np.pi, extent=[-0.5, 0.5, -0.5, 0.5])
In [34]: cbar = fig.colorbar(im, ax=ax)
In [35]: ax.set_xlabel(r'$k_x$')
Out[35]: Text(0.5, 0, '$k_x$')
In [36]: ax.set_ylabel(r'$k_y$');
# Compute the Chern number by integrating out the Berry curvature flux.
In [37]: chern_number = fluxes.sum() / (2 * np.pi)
In [38]: print('Chern number: %.4f' % chern_number)
Chern number: -0.0000
There are two peaks of Berry flux with opposite sign at K and K’. The Berry flux exactly compensates (it is an odd function under time-reversal, so the integral over the whole Brillouin zone vanishes). To get an overall Berry curvature, let us break time-reversal symmetry. We add spin-orbit coupling terms that point out-of the plane to all the bonds. We do not use the symmetry functionality of topwave in this case, because we want to explicitly break time-reversal symmetry by introducing fluxes à la Haldane (which regular spin-orbit coupling does not do).
The orientation of the couplings is generated automatically and arbitrary. We could just stare at the table with the couplings above to see how to choose the sign of the spin-orbit terms. We want to make sure the next-nearest neighbor couplings form closed flux loops. Let us plot how the lattice vectors connect the sublattices to see how to choose the signs of the spin-orbit terms.
In [39]: fig, ax = plt.subplots()
In [40]: origin_A = model.structure[0].frac_coords[:2]
In [41]: origin_B = model.structure[1].frac_coords[:2]
In [42]: for index_A, index_B in zip(range(3,6), range(6, 9)):
....: coupling_A, coupling_B = model.couplings[index_A], model.couplings[index_B]
....: arrow_A, arrow_B = coupling_A.lattice_vector[:2], coupling_B.lattice_vector[:2]
....: ax.arrow(*origin_A, *arrow_A, head_width=0.1, color='deeppink', length_includes_head=True)
....: ax.arrow(*origin_B, *arrow_B, head_width=0.1, color='red', length_includes_head=True)
....: origin_A += arrow_A
....: origin_B += arrow_B
....:
In [43]: ax.set_aspect('equal')
The pink triangle does not form a closed loop (no finite flux/winding. We can either flip the orientation of the coupling, or just account for the orientation of the bond when putting the spin-orbit interaction. Let us do the former and plot the triangles again.
In [44]: model.invert_coupling(3)
In [45]: fig, ax = plt.subplots()
In [46]: origin_A = model.structure[0].frac_coords[:2]
In [47]: origin_B = model.structure[1].frac_coords[:2]
In [48]: for index_A, index_B in zip(range(3,6), range(6, 9)):
....: coupling_A, coupling_B = model.couplings[index_A], model.couplings[index_B]
....: arrow_A, arrow_B = coupling_A.lattice_vector[:2], coupling_B.lattice_vector[:2]
....: ax.arrow(*origin_A, *arrow_A, head_width=0.1, color='deeppink', length_includes_head=True)
....: ax.arrow(*origin_B, *arrow_B, head_width=0.1, color='red', length_includes_head=True)
....: origin_A += arrow_A
....: origin_B += arrow_B
....:
In [49]: ax.set_aspect('equal')
Perfect! We the loops are closed and have the same sense of rotation. So let us assign some spin-orbit terms that point out of plane to all these bonds and plot the spectrum, and the flux again. We could choose the next-nearest neighbors based on their symmetry classification again, but let us select them based on distance this time (which we just read off the tables from above).
# this sets the strength of the complex hopping term
In [50]: lamda = 0.03
In [51]: model.set_spin_orbit(attribute_value=1.42,
....: vector=[0, 0, 1],
....: strength=lamda,
....: attribute='distance')
....:
In [52]: spec = tp.Spec(model, path)
In [53]: fig, (ax1, ax2) = plt.subplots(1, 2)
In [54]: for band in spec.energies.T:
....: ax1.plot(band, ls='-', lw=2.5, c='deeppink')
....:
In [55]: ax1.set_xticks(path.node_indices);
In [56]: ax1.set_xticklabels([r'$\Gamma$', 'M', 'K', r'$\Gamma$']);
In [57]: ax1.set_ylabel('Energy');
In [58]: fluxes = np.array([tp.Spec(model, plaquette).get_berry_phase([0]) for plaquette in plaquettes])
In [59]: im = ax2.imshow(fluxes.reshape((num_x, num_y)), origin='lower', cmap='PiYG', vmin=-np.pi, vmax=np.pi, extent=[-0.5, 0.5, -0.5, 0.5])
In [60]: cbar = fig.colorbar(im, ax=ax2)
In [61]: ax2.set_xlabel(r'$k_x$');
In [62]: ax2.set_ylabel(r'$k_y$');
In [63]: fig.set_size_inches(10, 3)
In [64]: plt.tight_layout()
Let’s focus on the plot on the left. A small gap has openend. Seems like we did everything right! Why does the Berry flux look like a QR-code though? With the spin orbit coupling, we introduced a spin-dependent term, so the size size of our Hilbert space was doubled. We are looking at four pairwise degenerate bands, not at two anymore. We can confirm by checking the number of eigenvalues at any k-point (or by calling the check_if_spinful-method).
In [65]: model.check_if_spinful()
Out[65]: True
In [66]: spec.energies[0].shape
Out[66]: (4,)
The (abelian) Berry curvature is not well-defined for degenerate bands. In the Haldane model, the spinless or rather spin-polarized case is considered. What we are looking at now is actually two copies of the Haldane model, the so-called Kane-Mele model. We can spin-polarize our system by applying a strong external magnetic field and only computing the Berry curvature for the lowest energy band, or we can use the set_spin_polarized method of the model to only consider the spin-up block of the Hamiltoninan.
In [67]: model.set_spin_polarized()
In [68]: spec = tp.Spec(model, path)
In [69]: fig, ax = plt.subplots()
In [70]: fluxes = np.array([tp.Spec(model, plaquette).get_berry_phase([0]) for plaquette in plaquettes])
In [71]: im = ax.imshow(fluxes.reshape((num_x, num_y)), origin='lower', cmap='PiYG', vmin=-np.pi, vmax=np.pi, extent=[-0.5, 0.5, -0.5, 0.5])
In [72]: cbar = fig.colorbar(im, ax=ax)
In [73]: ax.set_xlabel(r'$k_x$');
In [74]: ax.set_ylabel(r'$k_y$');
In [75]: chern_number = fluxes.sum() / (2 * np.pi)
In [76]: print('Chern number: %.4f' % chern_number)
Chern number: 1.0000
Instead of computing the Berry phase for a set of plaquettes that cover the two-dimensional Brillouin zone, we can also evaluate the Berry Curvature directly at each k-point. Because the Berry Curvature is often an ill-behaved function - it has very sharp peaks close to avoided crossings - it is often preferred to use the Wilson loop method as above. Another option is to track the evolution of the Wannier charge centers through the Brillouin zone, which is equivalent to computing the Berry phase for lines that we use to cover the Brillouin zone.
In [77]: nk = 13
In [78]: plane = tp.Plane(normal=[0, 0, 1], num_x=nk, num_y=nk)
In [79]: spec = tp.Spec(model, plane)
In [80]: berry_curvature_all_bands = tp.get_berry_curvature(spec, component='z')
In [81]: berry_curvature = berry_curvature_all_bands[:, 0].reshape(plane.shape)
In [82]: vmax = np.abs(berry_curvature).max()
In [83]: fig, (ax1, ax2) = plt.subplots(1, 2)
In [84]: im = ax1.imshow(berry_curvature, origin='lower', cmap='PiYG', vmin=-vmax, vmax=vmax, extent=[-0.5, 0.5, -0.5, 0.5])
In [85]: ax1.set_xlabel(r'$k_x$');
In [86]: ax1.set_ylabel(r'$k_y$');
In [87]: line_cover = tp.get_line_cover('z', 'x', 50, 100)
In [88]: spectra = [tp.Spec(model, line) for line in line_cover]
In [89]: charge_centers = np.array([spectrum.get_berry_phase(band_indices=[0]) for spectrum in spectra])
In [90]: ax2.scatter(line_cover[:, 0, 0], charge_centers, c='deeppink')
Out[90]: <matplotlib.collections.PathCollection at 0x7fb17f7000d0>
In [91]: ax2.set_aspect(1 / (2 * np.pi))
In [92]: ax2.set_xlabel(r'$k_x$')
Out[92]: Text(0.5, 0, '$k_x$')
In [93]: ax2.set_ylabel('Berry phase')
Out[93]: Text(0, 0.5, 'Berry phase')
In [94]: plt.tight_layout()
In [95]: chern_number = berry_curvature.sum() / (2 * np.pi) / nk**2
In [96]: print('Chern number: %.4f' % chern_number)
Chern number: -1.0918
Note that even though we used the same grid size as we had numbers of plaquettes, the Chern number for integrating out the Berry curvature directly is not as close to its correct value as the Wilson loop method. We can also nicely see that charge is ‘pumped’ throughout one adiabatic cycle, which means taking kx from -pi to pi.
We have confirmed, using three different methods, that the model is in a Chern Hall insulating phase. This means that even though our model is an insulator, we expect a quantized Hall current on the boundaries of the system. To see this, we perform a slab calculation. We create a supercell in one spatial direction and open the boundaries along that direction. We then use the get_projections function to project the wave functions onto the first and last unit cell of the slab that we created and plot the color of the slab bands as a function of that projector.
In [97]: num_cells = 10
In [98]: supercell_model = tp.get_supercell(model, [1, num_cells, 1])
In [99]: supercell_model.set_open_boundaries('y')
In [100]: path = tp.Path([[0, 0, 0], [1, 0, 0]], 101)
In [101]: spec = tp.Spec(supercell_model, path)
In [102]: projections = tp.get_projections(spec, {'unit_cell_y' : [0, num_cells - 1]})
In [103]: fig, ax = plt.subplots()
In [104]: for band, projection in zip(spec.energies.T, projections.T):
.....: ax.scatter(path.kpoints[:, 0], band, alpha=0.4, c=projection, cmap=plt.get_cmap('RdPu'), vmin=0, vmax=1)
.....:
In [105]: ax.set_xlabel(r'$k_x$')
Out[105]: Text(0.5, 0, '$k_x$')
In [106]: ax.set_ylabel(r'Energy')
Out[106]: Text(0, 0.5, 'Energy')
In [107]: plt.tight_layout()
A gapless mode has appeared! Even though our bulk model was an insulator, the boundary has become metallic. We can also see that these boundary modes are located at the ends of the slab that we created. Let’s compute the electron density on a finite sheet with open boundary conditions in x, y and both directions.
In [108]: import matplotlib.transforms as mtransforms
In [109]: num_cells = 10
In [110]: supercell_model_x = tp.get_supercell(model, [num_cells, num_cells, 1])
In [111]: supercell_model_x.set_open_boundaries('x')
In [112]: supercell_model_y = tp.get_supercell(model, [num_cells, num_cells, 1])
In [113]: supercell_model_y.set_open_boundaries('y')
In [114]: supercell_model_xy = tp.get_supercell(model, [num_cells, num_cells, 1])
In [115]: supercell_model_xy.set_open_boundaries('xy')
In [116]: spec_x = tp.Spec(supercell_model_x, [[0, 0, 0]])
In [117]: spec_y = tp.Spec(supercell_model_y, [[0, 0, 0]])
In [118]: spec_xy = tp.Spec(supercell_model_xy, [[0, 0, 0]])
In [119]: density_x = spec_x.get_particle_density(0)[:, :, 0, :, 0].sum(axis=2)
In [120]: density_y = spec_y.get_particle_density(0)[:, :, 0, :, 0].sum(axis=2)
In [121]: density_xy = spec_xy.get_particle_density(0)[:, :, 0, :, 0].sum(axis=2)
In [122]: fig, (ax1, ax2, ax3) = plt.subplots(1, 3)
In [123]: fig.set_size_inches(6, 2)
In [124]: for ax, density in zip([ax1, ax2, ax3], [density_x, density_y, density_xy]):
.....: ax.imshow(density, origin='lower', cmap='PiYG', extent=[-0.5, 0.5, -0.5, 0.5],
.....: transform=mtransforms.Affine2D().skew(np.pi / 6, 0) + ax.transData)
.....: ax.set_xlim(-0.8, 0.8)
.....: ax.set_ylim(-0.8, 0.8)
.....: ax.spines[['left', 'right', 'bottom', 'top']].set_visible(False)
.....: ax.set_xticks([])
.....: ax.set_yticks([])
.....:
In [125]: plt.tight_layout()
