Benalcazar-Bernevig-Hughes¶
The two-dimensional Benalcazar-Bernevig-Hughes (BBH) model is a minimal model for higher order topological insulator with quadrupolar order.
Let’s define a function that returns the BBH model with a given set of parameters.
In [1]: from pymatgen.core.structure import Lattice, Structure
In [2]: def get_BBH_model(gamma: float, lamda: float, delta: float) -> tp.TightBindingModel:
...: a = b = 1
...: c = 10
...: lattice = np.diag([a, b, c])
...: coords = [[0.5, 0.5, 0], [0, 0, 0], [0, 0.5, 0], [0.5, 0, 0]]
...: structure = Structure.from_spacegroup(sg=1, lattice=lattice, species=['C'] * 4, coords=coords)
...: model = tp.TightBindingModel(structure)
...: model.generate_couplings(0.5, 1)
...: model.invert_coupling(0)
...: model.set_coupling(0, gamma)
...: model.set_coupling(1, gamma)
...: model.set_coupling(6, gamma)
...: model.invert_coupling(7)
...: model.set_coupling(7, -gamma)
...: model.invert_coupling(2)
...: model.set_coupling(2, lamda)
...: model.set_coupling(3, lamda)
...: model.set_coupling(4, lamda)
...: model.invert_coupling(5)
...: model.set_coupling(5, -lamda)
...: model.set_onsite_scalar(0, delta)
...: model.set_onsite_scalar(1, delta)
...: model.set_onsite_scalar(2, -delta)
...: model.set_onsite_scalar(3, -delta)
...: return model
...:
Let’s print the couplings and check that everything looks like in the reference (the sign structure of the hoppings is important).
In [3]: model = get_BBH_model(0.5, 1, 0.001)
In [4]: 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.5 │ [-0. -0. -0.] │ [0. 0.5 0. ] │ 3 │ 0 │ 0 │ 0 │ 0.5 │ [0. 0. 0.] │ [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]] │
├─────────┼──────────────────┼──────────────────────┼────────────┼──────────────────┼─────────────────────┼─────────┼────────────┼─────────┼────────────┼────────────┼─────────────────────┼─────────────────────────────────────────────────────┤
│ 1 │ 1 │ x, y, z │ 0.5 │ [0. 0. 0.] │ [-0.5 0. 0. ] │ 0 │ 0 │ 2 │ 0 │ 0.5 │ [0. 0. 0.] │ [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]] │
├─────────┼──────────────────┼──────────────────────┼────────────┼──────────────────┼─────────────────────┼─────────┼────────────┼─────────┼────────────┼────────────┼─────────────────────┼─────────────────────────────────────────────────────┤
│ 2 │ 2 │ x, y, z │ 0.5 │ [-0. -1. -0.] │ [0. 0.5 0. ] │ 3 │ 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]] │
├─────────┼──────────────────┼──────────────────────┼────────────┼──────────────────┼─────────────────────┼─────────┼────────────┼─────────┼────────────┼────────────┼─────────────────────┼─────────────────────────────────────────────────────┤
│ 3 │ 3 │ x, y, z │ 0.5 │ [1. 0. 0.] │ [-0.5 0. 0. ] │ 0 │ 0 │ 2 │ 0 │ 1 │ [0. 0. 0.] │ [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]] │
├─────────┼──────────────────┼──────────────────────┼────────────┼──────────────────┼─────────────────────┼─────────┼────────────┼─────────┼────────────┼────────────┼─────────────────────┼─────────────────────────────────────────────────────┤
│ 4 │ 4 │ x, y, z │ 0.5 │ [-1. 0. 0.] │ [0.5 0. 0. ] │ 1 │ 0 │ 3 │ 0 │ 1 │ [0. 0. 0.] │ [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]] │
├─────────┼──────────────────┼──────────────────────┼────────────┼──────────────────┼─────────────────────┼─────────┼────────────┼─────────┼────────────┼────────────┼─────────────────────┼─────────────────────────────────────────────────────┤
│ 5 │ 5 │ x, y, z │ 0.5 │ [-0. 1. -0.] │ [ 0. -0.5 0. ] │ 2 │ 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]] │
├─────────┼──────────────────┼──────────────────────┼────────────┼──────────────────┼─────────────────────┼─────────┼────────────┼─────────┼────────────┼────────────┼─────────────────────┼─────────────────────────────────────────────────────┤
│ 6 │ 6 │ x, y, z │ 0.5 │ [0. 0. 0.] │ [0.5 0. 0. ] │ 1 │ 0 │ 3 │ 0 │ 0.5 │ [0. 0. 0.] │ [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]] │
├─────────┼──────────────────┼──────────────────────┼────────────┼──────────────────┼─────────────────────┼─────────┼────────────┼─────────┼────────────┼────────────┼─────────────────────┼─────────────────────────────────────────────────────┤
│ 7 │ 7 │ x, y, z │ 0.5 │ [-0. -0. -0.] │ [ 0. -0.5 0. ] │ 2 │ 0 │ 1 │ 0 │ -0.5 │ [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 │
╞═════════╪══════════════════╪══════════════════════╪════════════╪══════════════════╪═════════════════════╪═════════╪════════════╪═════════╪════════════╪════════════╪═════════════════════╪══════════╡
╘═════════╧══════════════════╧══════════════════════╧════════════╧══════════════════╧═════════════════════╧═════════╧════════════╧═════════╧════════════╧════════════╧═════════════════════╧══════════╛
Next, we create a path through the two-dimensional Brillouin zone and plot the band structure.
In [5]: path = tp.Path([[0, 0, 0], [0.5, 0, 0], [0.5, 0.5, 0], [0, 0, 0]])
In [6]: spec = tp.Spec(get_BBH_model(0.5, 1, 0.001), path)
In [7]: fig, ax = plt.subplots()
In [8]: for band in spec.energies.T:
...: ax.plot(band, c='k', ls='-')
...:
In [9]: ax.set_xticks(path.node_indices);
In [10]: ax.set_xlim(0, path.node_indices[-1]);
In [11]: ax.set_xticklabels([r'$\Gamma$', 'X', 'M', r'$\Gamma$']);
It’s an insulator. However, there is a gap-closing point that marks a topological phasetransition. We can calculate the dispersion on a finite sheet at the Gamma point and sweep through the first parameter of the model to visualize the topological phase transition. We also project the eigenstates onto the boundaries of the system to confirm that the zero energy modes are actually edge modes.
In [12]: gammas = np.linspace(-1.5, 1.5, 51)
In [13]: num_cells = 6
In [14]: scaling_factors = [num_cells, num_cells, 1]
In [15]: delta = 0.001
In [16]: supercell = tp.get_supercell(get_BBH_model(gammas[0], 1, delta), scaling_factors)
In [17]: supercell.set_open_boundaries('xy')
In [18]: spec = tp.Spec(supercell, [0, 0, 0])
In [19]: energies = spec.energies
In [20]: wavefunctions = spec.psi
In [21]: projections = tp.get_projections(spec, {'unit_cell_x': [0, num_cells - 1], 'unit_cell_y': [0, num_cells - 1]})
In [22]: for gamma in gammas[1:]:
....: supercell = tp.get_supercell(get_BBH_model(gamma, 1, delta), scaling_factors)
....: supercell.set_open_boundaries('xy')
....: spec = tp.Spec(supercell, [0, 0, 0])
....: energies = np.vstack((energies, spec.energies))
....: wavefunctions = np.vstack((wavefunctions, spec.psi))
....: projections = np.vstack((projections, tp.get_projections(spec, {'unit_cell_x': [0, num_cells - 1], 'unit_cell_y': [0, num_cells - 1]})))
....:
In [23]: fig, ax = plt.subplots()
In [24]: for band, projection in zip(energies.T, projections.T):
....: ax.plot(gammas, band, ls='-', c='midnightblue')
....: ax.scatter(gammas, band, c=projection, alpha=0.5, cmap=plt.get_cmap('Reds'), vmin=0, vmax=1)
....:
In [25]: ax.set_xlim(gammas[0], gammas[-1]);
In [26]: ax.set_xlabel(r'$\gamma / \lambda$');
In [27]: ax.set_ylabel(r'Energy');
We found surface modes! In the figure above it looks like the gap closing is slightly away from gamma / lambda = 1. However that is just because of the small size of our finite system.
If there are surface modes, we should be able to diagnose that by doing a bulk calculation and calculating the Chern number by e.g. tracking the evolution of the Wannier charge centers throughout the Brillouin zone. Let’s try that.
#line_cover = tp.get_ine
#model = tp.get_supercell(get_BBH_model(0.5, 1, delta), scaling_factors)
#supercell.set_open_boundaries('xy')
Weird. Why is there no well we have to calculate
In [28]: num_cells = 10
In [29]: scaling_factors = [num_cells, num_cells, 1]
In [30]: delta = 0.001
In [31]: supercell = tp.get_supercell(get_BBH_model(0.5, 1, delta), scaling_factors)
In [32]: supercell.set_open_boundaries('xy')
In [33]: spec = tp.Spec(supercell, [0, 0, 0])
In [34]: density = spec.get_particle_density(0)
# integrate out the sublattice degree of freedom
In [35]: density = density[:, :, 0, :, 0].sum(axis=2)
In [36]: fig, ax = plt.subplots()
In [37]: ax.imshow(density, origin='lower', cmap='bwr', extent=[1, num_cells, 1, num_cells])
Out[37]: <matplotlib.image.AxesImage at 0x7fb182a3dcd0>
In [38]: ax.set_xlabel('x')
Out[38]: Text(0.5, 0, 'x')
In [39]: ax.set_ylabel('y')
Out[39]: Text(0, 0.5, 'y')
In [40]: ax.set_title('Electron Density')
Out[40]: Text(0.5, 1.0, 'Electron Density')
