Coverage for tbkit/higher_order.py: 100%
160 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-15 06:13 +0000
1r"""
2Higher-order topology: nested Wilson loops, Wannier-sector polarizations,
3the quadrupole moment and corner charges (Benalcazar, Bernevig and Hughes,
4Science 357, 61 (2017); Phys. Rev. B 96, 245115 (2017)).
6A two-dimensional insulator can have no dipole moment (no Chern number, no
7bulk polarization) and still be topological: its hybrid Wannier functions
8split into two gapped *Wannier sectors*, each carrying a quantized
9polarization of its own, and the bulk quadrupole moment
11.. math::
13 q_{xy} = 2\,p_y^{\nu_x^-}\,p_x^{\nu_y^-} \pmod 1
15(1/2 or 0 with reflection symmetries) shows up as a fractional charge
16:math:`\pm e/2` bound to each corner of a finite flake, with zero-energy
17corner states at the corners in the chiral-symmetric limit.
19The recipe of the nested Wilson loop, for occupied bands :math:`|u^n_{\mathbf{k}}\rangle`:
211. The Wilson loop along :math:`x` from every base point :math:`\mathbf{k}`,
22 :math:`W_{x,\mathbf{k}}|\nu^j_{x,\mathbf{k}}\rangle = e^{-2\pi i\nu_x^j(k_y)}|\nu^j_{x,\mathbf{k}}\rangle`;
23 its phases :math:`\nu_x^j(k_y)` are the Wannier bands (see
24 :meth:`tbkit.kspace.KSpace.wannier_flow`).
252. The Wannier-sector states :math:`|w^j_{x,\mathbf{k}}\rangle = \sum_n|u^n_{\mathbf{k}}\rangle[\nu^j_{x,\mathbf{k}}]^n`
26 of the sector :math:`\nu_x^-` (the Wannier bands below the Wannier gap).
273. The Wilson loop of those states along :math:`y` (the *nested* Wilson
28 loop), whose phase, averaged over :math:`k_x`, is the Wannier-sector
29 polarization :math:`p_y^{\nu_x^-}`.
31The Bloch matrices of the whole mesh come at once from
32:meth:`KSpace._bloch_ham <tbkit.kspace.KSpace>`; the links carry the orbital
33positions (``positions=True``), as in *KSpace.berry_phase*.
34"""
35from __future__ import annotations
37import numpy as np
38from numpy.typing import NDArray
39import scipy.linalg as LA
41import tbkit.error_handling as error_handling
42from tbkit.kspace import KSpace
43from tbkit.meanfield import _occupy
44from tbkit.lattice import Lattice
47PI = np.pi
50def bbh_model(gamma: float = 0.5, lam: float = 1., delta: float = 0.) -> KSpace:
51 r'''
52 Get the Benalcazar-Bernevig-Hughes (BBH) quadrupole insulator: four
53 orbitals per square cell, intra-cell hoppings :math:`\gamma`, inter-cell
54 hoppings :math:`\lambda`, and a :math:`\pi` flux through every plaquette,
56 .. math::
58 H(\mathbf{k}) = [\gamma+\lambda\cos k_x]\Gamma_4 + \lambda\sin k_x\,\Gamma_3
59 + [\gamma+\lambda\cos k_y]\Gamma_2 + \lambda\sin k_y\,\Gamma_1 + \delta\,\Gamma_0\, ,
61 with :math:`\Gamma_0 = \tau_3\sigma_0`, :math:`\Gamma_k = -\tau_2\sigma_k`
62 (:math:`k = 1, 2, 3`), :math:`\Gamma_4 = \tau_1\sigma_0` (orbital
63 :math:`2\tau + \sigma`). The hoppings are read off :math:`H(\mathbf{k})`,
64 and the four orbitals sit at the corners of a square of side 1/2 centred
65 on the cell origin, each next to the orbitals it is bonded to by
66 :math:`\gamma`. For :math:`|\gamma| < |\lambda|` (and :math:`\delta = 0`) the model is a
67 quadrupole insulator, :math:`q_{xy} = 1/2`; for :math:`|\gamma| > |\lambda|` it is trivial.
68 :math:`\delta` breaks the chiral symmetry and splits the corner states
69 to :math:`\pm\delta`.
71 :param gamma: Real number. Default value 0.5. Intra-cell hopping.
72 :param lam: Real number. Default value 1. Inter-cell hopping.
73 :param delta: Real number. Default value 0. Onsite :math:`\pm\delta`.
75 :returns:
76 * **bbh** -- **KSpace** instance, 4 orbitals, lattice constant 1.
77 '''
78 error_handling.real_number(gamma, 'gamma')
79 error_handling.real_number(lam, 'lam')
80 error_handling.real_number(delta, 'delta')
81 s0, sx = np.eye(2), np.array([[0., 1.], [1., 0.]])
82 sy, sz = np.array([[0., -1j], [1j, 0.]]), np.diag([1., -1.])
83 g1, g2, g3 = (-np.kron(sy, s) for s in (sx, sy, sz))
84 g4 = np.kron(sx, s0)
85 onsite = gamma * (g4 + g2)
86 bx = lam / 2 * (g4 - 1j * g3) # lam cos(kx) G4 + lam sin(kx) G3
87 by = lam / 2 * (g2 - 1j * g1)
88 r = 0.25
89 # orbital 0 top right, 1 bottom left, 2 top left, 3 bottom right
90 cell = [{'tag': t, 'r0': p} for t, p in zip('abcd', [(r, r), (-r, -r), (-r, r), (r, -r)])]
91 bbh = KSpace(Lattice(unit_cell=cell, prim_vec=[(1., 0.), (0., 1.)]))
92 hops = []
93 for mat, R in ((onsite, (0, 0)), (bx, (1, 0)), (by, (0, 1))):
94 for i in range(4):
95 for j in range(4):
96 if abs(mat[i, j]) > 1e-14 and (R != (0, 0) or i < j):
97 t = complex(mat[i, j])
98 hops.append({'i': i, 'j': j, 'R': R, 't': t.real if t.imag == 0 else t})
99 bbh.set_hopping(hops)
100 bbh.set_onsite({'a': float(delta), 'b': float(delta), 'c': -float(delta), 'd': -float(delta)})
101 return bbh
104def _check_model(ks: KSpace, bands) -> list[int]:
105 '''
106 Private function. Validate a 2D Hermitian KSpace and its band list.
107 '''
108 error_handling.kspace(ks, KSpace)
109 error_handling.dim_exact(ks.dim, 2)
110 error_handling.hermitian_kspace(ks.is_hermitian(), ks._overlap_hop)
111 if isinstance(bands, int):
112 bands = [bands]
113 error_handling.band_indices(bands, ks.norb)
114 return bands
117def _occupied_mesh(ks: KSpace, bands: list[int], nk: tuple[int, int]) -> NDArray[np.complex128]:
118 '''
119 Private function. Eigenvectors of *bands* on the uniform mesh, shape
120 (n1, n2, norb, len(bands)), from one vectorized Bloch sum.
121 '''
122 n1, n2 = nk
123 fracs = np.stack(np.meshgrid(np.arange(n1) / n1, np.arange(n2) / n2, indexing='ij'), axis=-1)
124 k_cart = (fracs.reshape(-1, 2) @ ks.rec_vec_k) @ ks.k_basis.T
125 _, vec = np.linalg.eigh(ks._bloch_ham(k_cart))
126 return vec[:, :, bands].reshape(n1, n2, ks.norb, len(bands))
129def _link_phase(ks: KSpace, direction: int, n: int, positions: bool) -> NDArray[np.complex128]:
130 '''
131 Private function. The phases exp(-i dk.tau) carried by each link of a
132 loop of n steps along b_direction (ones without positions).
133 '''
134 if not positions:
135 return np.ones(ks.norb, 'c16')
136 dk = ks.k_basis @ (ks.rec_vec_k[direction] / n)
137 return np.exp(-1j * ks.orbital_positions() @ dk)
140def _unitary(mat: NDArray[np.complex128]) -> NDArray[np.complex128]:
141 '''
142 Private function. Unitary part u vh of each matrix (last two axes).
143 '''
144 u, _, vh = np.linalg.svd(mat)
145 return u @ vh
148def _wilson_all(vec: NDArray[np.complex128], phase: NDArray[np.complex128]) -> NDArray[np.complex128]:
149 r'''
150 Private function. Wilson loops along the first axis of *vec* (shape
151 (n, m, norb, nb)) from every base point: W[i, j] = F_i F_{i+1} ... F_{i-1},
152 F_i = <u_i|e^{-i dk tau}|u_{i+1}> (unitarized). Shape (n, m, nb, nb).
153 '''
154 n = vec.shape[0]
155 nxt = np.roll(vec, -1, axis=0)
156 links = _unitary(vec.conj().transpose(0, 1, 3, 2) @ (phase[:, None] * nxt))
157 out = np.zeros(links.shape, 'c16')
158 for i in range(n):
159 w = np.broadcast_to(np.eye(links.shape[-1], dtype='c16'), links.shape[1:]).copy()
160 for m in range(n):
161 w = w @ links[(i + m) % n]
162 out[i] = w
163 return out
166def wannier_bands(
167 ks: KSpace, bands: int | list[int], nk: int | tuple[int, int] = 40, direction: int = 0,
168 positions: bool = True,
169) -> tuple[NDArray[np.float64], NDArray[np.complex128]]:
170 r'''
171 Get the Wannier bands :math:`\nu_{d}^j(\mathbf{k})` of a group of bands
172 of a 2D model, and the Wannier-sector states, from the Wilson loop
173 :math:`W_{d,\mathbf{k}}` along :math:`\mathbf{b}_d` starting from every
174 point :math:`\mathbf{k}` of the mesh. The Wannier bands do not depend on
175 the base point along :math:`d` (only the eigenvectors do); they are the
176 hybrid Wannier centres of *KSpace.wannier_centers*, in units of
177 :math:`\mathbf{a}_d`, here in :math:`(-1/2, 1/2]` and sorted.
179 :param ks: **KSpace** instance, 2D, Hermitian, without overlap.
180 :param bands: Band index or list of band indices (e.g. the occupied bands).
181 :param nk: Positive integer, or tuple of two. Default value 40. Mesh
182 along :math:`\mathbf{b}_1, \mathbf{b}_2`.
183 :param direction: 0 or 1. Default value 0. Direction :math:`d` of the Wilson loop.
184 :param positions: Boolean. Default value True. Include the orbital
185 positions in the links (see *KSpace.berry_phase*).
187 :returns:
188 * **nu** -- Real ndarray, shape (n1, n2, len(bands)). Wannier bands at
189 every base point, sorted, in :math:`(-1/2, 1/2]`.
190 * **w** -- Complex ndarray, shape (n1, n2, norb, len(bands)). The
191 Wannier-sector states :math:`|w^j_{d,\mathbf{k}}\rangle`, column *j*
192 for ``nu[..., j]``.
193 '''
194 bands = _check_model(ks, bands)
195 error_handling.nk(nk, 2)
196 error_handling.direction(direction, 2)
197 error_handling.boolean(positions, 'positions')
198 if isinstance(nk, int):
199 nk = (nk, nk)
200 vec = _occupied_mesh(ks, bands, nk)
201 phase = _link_phase(ks, direction, nk[direction], positions)
202 if direction == 1:
203 vec = vec.transpose(1, 0, 2, 3)
204 wil = _wilson_all(vec, phase)
205 val, eig = np.linalg.eig(wil)
206 nu = -np.angle(val) / (2 * PI)
207 nu = np.where(nu <= -0.5, nu + 1., nu)
208 order = np.argsort(nu, axis=-1, kind='stable')
209 nu = np.take_along_axis(nu, order, axis=-1)
210 eig = np.take_along_axis(eig, order[..., None, :], axis=-1)
211 # W is unitary: orthonormalize the eigenvectors (degenerate ones included)
212 eig = np.linalg.qr(eig)[0]
213 w = vec @ eig
214 if direction == 1:
215 nu, w = nu.transpose(1, 0, 2), w.transpose(1, 0, 2, 3)
216 return nu, w
219def wannier_sector_polarization(
220 ks: KSpace, bands: int | list[int], sector: int | list[int] | None = None,
221 nk: int | tuple[int, int] = 40, direction: int = 0, positions: bool = True,
222) -> float:
223 r'''
224 Get the polarization of a Wannier sector, :math:`p_{d'}^{\nu_d}`, by the
225 nested Wilson loop of Benalcazar, Bernevig and Hughes (see the module
226 docstring): the Wannier-sector states of *sector* (from *wannier_bands*
227 along *direction* :math:`d`) are carried around the Brillouin zone along
228 the other direction :math:`d'`,
230 .. math::
232 p_{d'}^{\nu_d} = -\frac{1}{N_d}\sum_{k_d}\frac{1}{2\pi}\arg\det
233 \prod_{k_{d'}}\langle w_{d,\mathbf{k}}|e^{-i\delta\mathbf{k}\cdot\boldsymbol\tau}|w_{d,\mathbf{k}+\delta\mathbf{k}_{d'}}\rangle
234 \pmod 1\, ,
236 in units of :math:`\mathbf{a}_{d'}`. The sector must be separated from
237 the other Wannier bands by a Wannier gap at every :math:`\mathbf{k}`
238 (a ValueError is raised otherwise). For the BBH model it is 1/2 in the
239 quadrupole phase and 0 in the trivial one.
241 :param ks: **KSpace** instance, 2D, Hermitian, without overlap.
242 :param bands: Band index or list of band indices (the occupied bands).
243 :param sector: Integer or list of integers. Default value None (the
244 lower half, :math:`\nu_d^-`). Indices of the Wannier bands (sorted
245 in :math:`(-1/2, 1/2]`, see *wannier_bands*) forming the sector.
246 :param nk: Positive integer, or tuple of two. Default value 40. Mesh.
247 :param direction: 0 or 1. Default value 0. Direction :math:`d` of the
248 first Wilson loop; the polarization is along the other one.
249 :param positions: Boolean. Default value True. See *wannier_bands*.
251 :returns:
252 * **p** -- Real number in :math:`[0, 1)`.
253 '''
254 bands = _check_model(ks, bands)
255 if sector is None:
256 sector = list(range(max(1, len(bands) // 2)))
257 if isinstance(sector, int):
258 sector = [sector]
259 error_handling.band_indices(sector, len(bands))
260 nu, w = wannier_bands(ks, bands, nk, direction, positions)
261 error_handling.wannier_gap(nu, sector)
262 other = 1 - direction
263 w = w[..., sector]
264 n_other = w.shape[other]
265 phase = _link_phase(ks, other, n_other, positions)
266 if other == 0:
267 w = w.transpose(1, 0, 2, 3)
268 # w: (n_d, n_other, norb, ns); links along the second axis
269 nxt = np.roll(w, -1, axis=1)
270 links = _unitary(w.conj().transpose(0, 1, 3, 2) @ (phase[:, None] * nxt))
271 loop = np.broadcast_to(np.eye(len(sector), dtype='c16'), (w.shape[0], len(sector), len(sector))).copy()
272 for m in range(n_other):
273 loop = loop @ links[:, m]
274 # the phases are continuous in k_d: unwrap them before averaging
275 p = -np.unwrap(np.angle(np.linalg.det(loop))) / (2 * PI)
276 return _mod1(float(np.mean(p)))
279def _mod1(x: float) -> float:
280 '''
281 Private function. x modulo 1, in [0, 1), with 1 - 1e-10 < x < 1 taken as 0.
282 '''
283 x = x % 1.
284 return 0. if x > 1. - 1e-10 else x
287def quadrupole_moment(
288 ks: KSpace, bands: int | list[int], nk: int | tuple[int, int] = 40, positions: bool = True,
289) -> float:
290 r'''
291 Get the bulk quadrupole moment of a 2D insulator with reflection
292 symmetries, from its Wannier-sector polarizations (Benalcazar, Bernevig
293 and Hughes 2017):
295 .. math::
297 q_{xy} = 2\,p_y^{\nu_x^-}\,p_x^{\nu_y^-} \pmod 1\, ,
299 in units of :math:`e` per unit cell area. With the reflections
300 :math:`M_x, M_y` (and :math:`C_4`, or the chiral symmetry, of the BBH
301 model) both polarizations are quantized to 0 or 1/2, and so is
302 :math:`q_{xy}`. The occupied bands must have a Wannier gap in both
303 directions (a ValueError is raised otherwise, see
304 *wannier_sector_polarization*).
306 :param ks: **KSpace** instance, 2D, Hermitian, without overlap.
307 :param bands: List of band indices (the occupied bands, at least two).
308 :param nk: Positive integer, or tuple of two. Default value 40. Mesh.
309 :param positions: Boolean. Default value True. See *wannier_bands*.
311 :returns:
312 * **q_xy** -- Real number in :math:`[0, 1)`.
313 '''
314 p_y = wannier_sector_polarization(ks, bands, None, nk, 0, positions)
315 p_x = wannier_sector_polarization(ks, bands, None, nk, 1, positions)
316 # each p is taken in (-1/2, 1/2] first, so that 1 - eps counts as 0
317 p_y, p_x = (p - 1. if p > 0.5 + 1e-9 else p for p in (p_y, p_x))
318 return _mod1(2 * p_y * p_x)
321def flake_positions(ks: KSpace, n_cells: int | tuple[int, ...]) -> NDArray[np.float64]:
322 r'''
323 Get the positions of the orbitals of the finite sample of
324 *KSpace.finite_ham*, in its order: orbital *o* of cell
325 :math:`(n_1, n_2, \dots)` at :math:`\sum_i n_i\mathbf{a}_i + \boldsymbol\tau_o`.
327 :param ks: **KSpace** instance.
328 :param n_cells: Positive integer, or tuple of *dim* positive integers.
330 :returns:
331 * **positions** -- Real ndarray, shape (N*norb, space_dim).
332 '''
333 error_handling.kspace(ks, KSpace)
334 error_handling.nk(n_cells, ks.dim)
335 if isinstance(n_cells, int):
336 n_cells = (n_cells,) * ks.dim
337 # cell (n1, n2, n3) is row ((n3 N2 + n2) N1 + n1): n1 runs fastest
338 grids = np.meshgrid(*[np.arange(n) for n in n_cells[::-1]], indexing='ij')
339 cells = np.stack([g.ravel() for g in grids[::-1]], axis=1)
340 origin = cells @ np.array(ks.lat.prim_vec, dtype='f8')
341 return (origin[:, None, :] + ks.orbital_positions()[None]).reshape(-1, ks.space_dim)
344class CornerCharges():
345 r'''
346 Charges of the four corners of a finite 2D flake, see *corner_charges*.
348 :ivar energies: Real ndarray. Spectrum of the flake.
349 :ivar density: Real ndarray. Electron density :math:`\rho_i` on every orbital.
350 :ivar positions: Real ndarray, shape (N*norb, 2). Orbital positions.
351 :ivar background: Real number. Uniform background :math:`\bar\rho`
352 (electrons per orbital) that neutralizes the flake.
353 :ivar e_fermi: Real number. Fermi level.
354 :ivar charges: Real ndarray, shape (2, 2). Excess electron number
355 :math:`\sum_{i\in Q}(\rho_i - \bar\rho)` in each quadrant *Q*:
356 ``charges[a, b]`` for the half :math:`a` (0: low, 1: high) along
357 :math:`\mathbf{a}_1` and :math:`b` along :math:`\mathbf{a}_2`.
358 '''
360 def __init__(self, energies, density, positions, background, e_fermi, charges) -> None:
361 self.energies = energies
362 self.density = density
363 self.positions = positions
364 self.background = background
365 self.e_fermi = e_fermi
366 self.charges = charges
369def corner_charges(
370 ks: KSpace, n_cells: int | tuple[int, int], n_electrons: float | None = None,
371 temperature: float = 0.,
372) -> CornerCharges:
373 r'''
374 Get the corner charges of a finite flake of :math:`N_1\times N_2` cells
375 of a 2D model (*KSpace.finite_ham*, open boundaries): the electron
376 density :math:`\rho_i = \sum_n f_n|\psi_n(i)|^2` of the lowest states,
377 minus a uniform neutralizing background :math:`\bar\rho = N_e/(N\,\mathrm{norb})`,
378 summed over each quadrant of the flake (the halves being split at the
379 middle of the sample along :math:`\mathbf{a}_1` and
380 :math:`\mathbf{a}_2`). A quadrant charge is the corner charge of
381 Benalcazar, Bernevig and Hughes,
382 :math:`Q^{\mathrm{corner}} = p^{\mathrm{edge}}_x + p^{\mathrm{edge}}_y - q_{xy}`:
383 :math:`\pm 1/2` for the BBH quadrupole insulator (with
384 :math:`\delta \ne 0` to fill two of its four corner states), 0 for a
385 trivial insulator. The sign is that of the *electron* number (charge
386 in units of :math:`-e`).
388 At :math:`T = 0`, the electrons left for a degenerate level at the
389 Fermi energy are shared equally among its states (e.g. the four exact
390 zero modes of the BBH model at :math:`\delta = 0`, which then gives
391 zero corner charges): split the degeneracy, or use a temperature.
393 :param ks: **KSpace** instance, 2D, Hermitian, without overlap.
394 :param n_cells: Positive integer, or tuple of two. Cells along
395 :math:`\mathbf{a}_1, \mathbf{a}_2`.
396 :param n_electrons: Positive real. Default value None (half filling).
397 Number of electrons in the flake. An integer at T = 0.
398 :param temperature: Positive real or zero. Default value 0.
400 :returns:
401 * **result** -- :class:`CornerCharges`.
402 '''
403 _check_model(ks, 0)
404 error_handling.nk(n_cells, 2)
405 if isinstance(n_cells, int):
406 n_cells = (n_cells, n_cells)
407 n_states = int(np.prod(n_cells)) * ks.norb
408 if n_electrons is None:
409 n_electrons = n_states // 2
410 error_handling.electrons(n_electrons, n_states)
411 error_handling.positive_real_zero(temperature, 'temperature')
412 if temperature == 0:
413 error_handling.integer_electrons(n_electrons)
414 en, vec = LA.eigh(ks.finite_ham(n_cells))
415 f, _, mu = _occupy(en, en[:0], n_electrons, temperature)
416 density = (np.abs(vec) ** 2) @ f
417 background = n_electrons / n_states
418 positions = flake_positions(ks, n_cells)
419 # fractional coordinates along a1, a2, split at the middle of the sample
420 a = np.array(ks.lat.prim_vec, dtype='f8')
421 frac = np.linalg.solve(a.T, positions.T).T
422 mid = (np.array(n_cells) - 1) / 2 + np.mean(np.linalg.solve(a.T, ks.orbital_positions().T).T, axis=0)
423 high = frac > mid[None, :]
424 excess = density - background
425 charges = np.array([[excess[(high[:, 0] == bool(s1)) & (high[:, 1] == bool(s2))].sum()
426 for s2 in (0, 1)] for s1 in (0, 1)])
427 return CornerCharges(en, density, positions, background, float(mu), charges)