Coverage for tbkit/higher_order.py: 100%
160 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-10-30 11:04 +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 Parameters
72 ----------
73 gamma : Real number, default 0.5
74 Intra-cell hopping.
75 lam : Real number, default 1
76 Inter-cell hopping.
77 delta : Real number, default 0
78 Onsite :math:`\pm\delta`.
80 Returns
81 -------
82 bbh : **KSpace** instance, 4 orbitals, lattice constant 1
83 '''
84 error_handling.real_number(gamma, 'gamma')
85 error_handling.real_number(lam, 'lam')
86 error_handling.real_number(delta, 'delta')
87 s0, sx = np.eye(2), np.array([[0., 1.], [1., 0.]])
88 sy, sz = np.array([[0., -1j], [1j, 0.]]), np.diag([1., -1.])
89 g1, g2, g3 = (-np.kron(sy, s) for s in (sx, sy, sz))
90 g4 = np.kron(sx, s0)
91 onsite = gamma * (g4 + g2)
92 bx = lam / 2 * (g4 - 1j * g3) # lam cos(kx) G4 + lam sin(kx) G3
93 by = lam / 2 * (g2 - 1j * g1)
94 r = 0.25
95 # orbital 0 top right, 1 bottom left, 2 top left, 3 bottom right
96 cell = [{'tag': t, 'r0': p} for t, p in zip('abcd', [(r, r), (-r, -r), (-r, r), (r, -r)])]
97 bbh = KSpace(Lattice(unit_cell=cell, prim_vec=[(1., 0.), (0., 1.)]))
98 hops = []
99 for mat, R in ((onsite, (0, 0)), (bx, (1, 0)), (by, (0, 1))):
100 for i in range(4):
101 for j in range(4):
102 if abs(mat[i, j]) > 1e-14 and (R != (0, 0) or i < j):
103 t = complex(mat[i, j])
104 hops.append({'i': i, 'j': j, 'R': R, 't': t.real if t.imag == 0 else t})
105 bbh.set_hopping(hops)
106 bbh.set_onsite({'a': float(delta), 'b': float(delta), 'c': -float(delta), 'd': -float(delta)})
107 return bbh
110def _check_model(ks: KSpace, bands) -> list[int]:
111 '''
112 Private function. Validate a 2D Hermitian KSpace and its band list.
113 '''
114 error_handling.kspace(ks, KSpace)
115 error_handling.dim_exact(ks.dim, 2)
116 error_handling.hermitian_kspace(ks.is_hermitian(), ks._overlap_hop)
117 if isinstance(bands, int):
118 bands = [bands]
119 error_handling.band_indices(bands, ks.norb)
120 return bands
123def _occupied_mesh(ks: KSpace, bands: list[int], nk: tuple[int, int]) -> NDArray[np.complex128]:
124 '''
125 Private function. Eigenvectors of *bands* on the uniform mesh, shape
126 (n1, n2, norb, len(bands)), from one vectorized Bloch sum.
127 '''
128 n1, n2 = nk
129 fracs = np.stack(np.meshgrid(np.arange(n1) / n1, np.arange(n2) / n2, indexing='ij'), axis=-1)
130 k_cart = (fracs.reshape(-1, 2) @ ks.rec_vec_k) @ ks.k_basis.T
131 _, vec = np.linalg.eigh(ks._bloch_ham(k_cart))
132 return vec[:, :, bands].reshape(n1, n2, ks.norb, len(bands))
135def _link_phase(ks: KSpace, direction: int, n: int, positions: bool) -> NDArray[np.complex128]:
136 '''
137 Private function. The phases exp(-i dk.tau) carried by each link of a
138 loop of n steps along b_direction (ones without positions).
139 '''
140 if not positions:
141 return np.ones(ks.norb, 'c16')
142 dk = ks.k_basis @ (ks.rec_vec_k[direction] / n)
143 return np.exp(-1j * ks.orbital_positions() @ dk)
146def _unitary(mat: NDArray[np.complex128]) -> NDArray[np.complex128]:
147 '''
148 Private function. Unitary part u vh of each matrix (last two axes).
149 '''
150 u, _, vh = np.linalg.svd(mat)
151 return u @ vh
154def _wilson_all(vec: NDArray[np.complex128], phase: NDArray[np.complex128]) -> NDArray[np.complex128]:
155 r'''
156 Private function. Wilson loops along the first axis of *vec* (shape
157 (n, m, norb, nb)) from every base point: W[i, j] = F_i F_{i+1} ... F_{i-1},
158 F_i = <u_i|e^{-i dk tau}|u_{i+1}> (unitarized). Shape (n, m, nb, nb).
159 '''
160 n = vec.shape[0]
161 nxt = np.roll(vec, -1, axis=0)
162 links = _unitary(vec.conj().transpose(0, 1, 3, 2) @ (phase[:, None] * nxt))
163 out = np.zeros(links.shape, 'c16')
164 for i in range(n):
165 w = np.broadcast_to(np.eye(links.shape[-1], dtype='c16'), links.shape[1:]).copy()
166 for m in range(n):
167 w = w @ links[(i + m) % n]
168 out[i] = w
169 return out
172def wannier_bands(
173 ks: KSpace, bands: int | list[int], nk: int | tuple[int, int] = 40, direction: int = 0,
174 positions: bool = True,
175) -> tuple[NDArray[np.float64], NDArray[np.complex128]]:
176 r'''
177 Get the Wannier bands :math:`\nu_{d}^j(\mathbf{k})` of a group of bands
178 of a 2D model, and the Wannier-sector states, from the Wilson loop
179 :math:`W_{d,\mathbf{k}}` along :math:`\mathbf{b}_d` starting from every
180 point :math:`\mathbf{k}` of the mesh. The Wannier bands do not depend on
181 the base point along :math:`d` (only the eigenvectors do); they are the
182 hybrid Wannier centres of *KSpace.wannier_centers*, in units of
183 :math:`\mathbf{a}_d`, here in :math:`(-1/2, 1/2]` and sorted.
185 Parameters
186 ----------
187 ks : **KSpace** instance, 2D, Hermitian, without overlap
188 bands : Band index or list of band indices (e.g. the occupied bands)
189 nk : Positive integer, or tuple of two, default 40
190 Mesh
191 along :math:`\mathbf{b}_1, \mathbf{b}_2`.
192 direction : 0 or 1, default 0
193 Direction :math:`d` of the Wilson loop.
194 positions : Boolean, default True
195 Include the orbital
196 positions in the links (see *KSpace.berry_phase*).
198 Returns
199 -------
200 nu : Real ndarray, shape (n1, n2, len(bands))
201 Wannier bands at
202 every base point, sorted, in :math:`(-1/2, 1/2]`.
203 w : Complex ndarray, shape (n1, n2, norb, len(bands))
204 The
205 Wannier-sector states :math:`|w^j_{d,\mathbf{k}}\rangle`, column *j*
206 for ``nu[..., j]``.
207 '''
208 bands = _check_model(ks, bands)
209 error_handling.nk(nk, 2)
210 error_handling.direction(direction, 2)
211 error_handling.boolean(positions, 'positions')
212 if isinstance(nk, int):
213 nk = (nk, nk)
214 vec = _occupied_mesh(ks, bands, nk)
215 phase = _link_phase(ks, direction, nk[direction], positions)
216 if direction == 1:
217 vec = vec.transpose(1, 0, 2, 3)
218 wil = _wilson_all(vec, phase)
219 val, eig = np.linalg.eig(wil)
220 nu = -np.angle(val) / (2 * PI)
221 nu = np.where(nu <= -0.5, nu + 1., nu)
222 order = np.argsort(nu, axis=-1, kind='stable')
223 nu = np.take_along_axis(nu, order, axis=-1)
224 eig = np.take_along_axis(eig, order[..., None, :], axis=-1)
225 # W is unitary: orthonormalize the eigenvectors (degenerate ones included)
226 eig = np.linalg.qr(eig)[0]
227 w = vec @ eig
228 if direction == 1:
229 nu, w = nu.transpose(1, 0, 2), w.transpose(1, 0, 2, 3)
230 return nu, w
233def wannier_sector_polarization(
234 ks: KSpace, bands: int | list[int], sector: int | list[int] | None = None,
235 nk: int | tuple[int, int] = 40, direction: int = 0, positions: bool = True,
236) -> float:
237 r'''
238 Get the polarization of a Wannier sector, :math:`p_{d'}^{\nu_d}`, by the
239 nested Wilson loop of Benalcazar, Bernevig and Hughes (see the module
240 docstring): the Wannier-sector states of *sector* (from *wannier_bands*
241 along *direction* :math:`d`) are carried around the Brillouin zone along
242 the other direction :math:`d'`,
244 .. math::
246 p_{d'}^{\nu_d} = -\frac{1}{N_d}\sum_{k_d}\frac{1}{2\pi}\arg\det
247 \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
248 \pmod 1\, ,
250 in units of :math:`\mathbf{a}_{d'}`. The sector must be separated from
251 the other Wannier bands by a Wannier gap at every :math:`\mathbf{k}`
252 (a ValueError is raised otherwise). For the BBH model it is 1/2 in the
253 quadrupole phase and 0 in the trivial one.
255 Parameters
256 ----------
257 ks : **KSpace** instance, 2D, Hermitian, without overlap
258 bands : Band index or list of band indices (the occupied bands)
259 sector : Integer or list of integers
260 Default value None (the
261 lower half, :math:`\nu_d^-`). Indices of the Wannier bands (sorted
262 in :math:`(-1/2, 1/2]`, see *wannier_bands*) forming the sector.
263 nk : Positive integer, or tuple of two, default 40
264 Mesh.
265 direction : 0 or 1, default 0
266 Direction :math:`d` of the
267 first Wilson loop; the polarization is along the other one.
268 positions : Boolean, default True
269 See *wannier_bands*.
271 Returns
272 -------
273 p : Real number in :math:`[0, 1)`
274 '''
275 bands = _check_model(ks, bands)
276 if sector is None:
277 sector = list(range(max(1, len(bands) // 2)))
278 if isinstance(sector, int):
279 sector = [sector]
280 error_handling.band_indices(sector, len(bands))
281 nu, w = wannier_bands(ks, bands, nk, direction, positions)
282 error_handling.wannier_gap(nu, sector)
283 other = 1 - direction
284 w = w[..., sector]
285 n_other = w.shape[other]
286 phase = _link_phase(ks, other, n_other, positions)
287 if other == 0:
288 w = w.transpose(1, 0, 2, 3)
289 # w: (n_d, n_other, norb, ns); links along the second axis
290 nxt = np.roll(w, -1, axis=1)
291 links = _unitary(w.conj().transpose(0, 1, 3, 2) @ (phase[:, None] * nxt))
292 loop = np.broadcast_to(np.eye(len(sector), dtype='c16'), (w.shape[0], len(sector), len(sector))).copy()
293 for m in range(n_other):
294 loop = loop @ links[:, m]
295 # the phases are continuous in k_d: unwrap them before averaging
296 p = -np.unwrap(np.angle(np.linalg.det(loop))) / (2 * PI)
297 return _mod1(float(np.mean(p)))
300def _mod1(x: float) -> float:
301 '''
302 Private function. x modulo 1, in [0, 1), with 1 - 1e-10 < x < 1 taken as 0.
303 '''
304 x = x % 1.
305 return 0. if x > 1. - 1e-10 else x
308def quadrupole_moment(
309 ks: KSpace, bands: int | list[int], nk: int | tuple[int, int] = 40, positions: bool = True,
310) -> float:
311 r'''
312 Get the bulk quadrupole moment of a 2D insulator with reflection
313 symmetries, from its Wannier-sector polarizations (Benalcazar, Bernevig
314 and Hughes 2017):
316 .. math::
318 q_{xy} = 2\,p_y^{\nu_x^-}\,p_x^{\nu_y^-} \pmod 1\, ,
320 in units of :math:`e` per unit cell area. With the reflections
321 :math:`M_x, M_y` (and :math:`C_4`, or the chiral symmetry, of the BBH
322 model) both polarizations are quantized to 0 or 1/2, and so is
323 :math:`q_{xy}`. The occupied bands must have a Wannier gap in both
324 directions (a ValueError is raised otherwise, see
325 *wannier_sector_polarization*).
327 Parameters
328 ----------
329 ks : **KSpace** instance, 2D, Hermitian, without overlap
330 bands : List of band indices (the occupied bands, at least two)
331 nk : Positive integer, or tuple of two, default 40
332 Mesh.
333 positions : Boolean, default True
334 See *wannier_bands*.
336 Returns
337 -------
338 q_xy : Real number in :math:`[0, 1)`
339 '''
340 p_y = wannier_sector_polarization(ks, bands, None, nk, 0, positions)
341 p_x = wannier_sector_polarization(ks, bands, None, nk, 1, positions)
342 # each p is taken in (-1/2, 1/2] first, so that 1 - eps counts as 0
343 p_y, p_x = (p - 1. if p > 0.5 + 1e-9 else p for p in (p_y, p_x))
344 return _mod1(2 * p_y * p_x)
347def flake_positions(ks: KSpace, n_cells: int | tuple[int, ...]) -> NDArray[np.float64]:
348 r'''
349 Get the positions of the orbitals of the finite sample of
350 *KSpace.finite_ham*, in its order: orbital *o* of cell
351 :math:`(n_1, n_2, \dots)` at :math:`\sum_i n_i\mathbf{a}_i + \boldsymbol\tau_o`.
353 Parameters
354 ----------
355 ks : **KSpace** instance
356 n_cells : Positive integer, or tuple of *dim* positive integers
358 Returns
359 -------
360 positions : Real ndarray, shape (N*norb, space_dim)
361 '''
362 error_handling.kspace(ks, KSpace)
363 error_handling.nk(n_cells, ks.dim)
364 if isinstance(n_cells, int):
365 n_cells = (n_cells,) * ks.dim
366 # cell (n1, n2, n3) is row ((n3 N2 + n2) N1 + n1): n1 runs fastest
367 grids = np.meshgrid(*[np.arange(n) for n in n_cells[::-1]], indexing='ij')
368 cells = np.stack([g.ravel() for g in grids[::-1]], axis=1)
369 origin = cells @ np.array(ks.lat.prim_vec, dtype='f8')
370 return (origin[:, None, :] + ks.orbital_positions()[None]).reshape(-1, ks.space_dim)
373class CornerCharges():
374 r'''
375 Charges of the four corners of a finite 2D flake, see *corner_charges*.
377 Attributes
378 ----------
379 energies : Real ndarray
380 Spectrum of the flake.
381 density : Real ndarray
382 Electron density :math:`\rho_i` on every orbital.
383 positions : Real ndarray, shape (N*norb, 2)
384 Orbital positions.
385 background : Real number
386 Uniform background :math:`\bar\rho`
387 (electrons per orbital) that neutralizes the flake.
388 e_fermi : Real number
389 Fermi level.
390 charges : Real ndarray, shape (2, 2)
391 Excess electron number
392 :math:`\sum_{i\in Q}(\rho_i - \bar\rho)` in each quadrant *Q*:
393 ``charges[a, b]`` for the half :math:`a` (0: low, 1: high) along
394 :math:`\mathbf{a}_1` and :math:`b` along :math:`\mathbf{a}_2`.
395 '''
397 def __init__(self, energies, density, positions, background, e_fermi, charges) -> None:
398 self.energies = energies
399 self.density = density
400 self.positions = positions
401 self.background = background
402 self.e_fermi = e_fermi
403 self.charges = charges
406def corner_charges(
407 ks: KSpace, n_cells: int | tuple[int, int], n_electrons: float | None = None,
408 temperature: float = 0.,
409) -> CornerCharges:
410 r'''
411 Get the corner charges of a finite flake of :math:`N_1\times N_2` cells
412 of a 2D model (*KSpace.finite_ham*, open boundaries): the electron
413 density :math:`\rho_i = \sum_n f_n|\psi_n(i)|^2` of the lowest states,
414 minus a uniform neutralizing background :math:`\bar\rho = N_e/(N\,\mathrm{norb})`,
415 summed over each quadrant of the flake (the halves being split at the
416 middle of the sample along :math:`\mathbf{a}_1` and
417 :math:`\mathbf{a}_2`). A quadrant charge is the corner charge of
418 Benalcazar, Bernevig and Hughes,
419 :math:`Q^{\mathrm{corner}} = p^{\mathrm{edge}}_x + p^{\mathrm{edge}}_y - q_{xy}`:
420 :math:`\pm 1/2` for the BBH quadrupole insulator (with
421 :math:`\delta \ne 0` to fill two of its four corner states), 0 for a
422 trivial insulator. The sign is that of the *electron* number (charge
423 in units of :math:`-e`).
425 At :math:`T = 0`, the electrons left for a degenerate level at the
426 Fermi energy are shared equally among its states (e.g. the four exact
427 zero modes of the BBH model at :math:`\delta = 0`, which then gives
428 zero corner charges): split the degeneracy, or use a temperature.
430 Parameters
431 ----------
432 ks : **KSpace** instance, 2D, Hermitian, without overlap
433 n_cells : Positive integer, or tuple of two
434 Cells along
435 :math:`\mathbf{a}_1, \mathbf{a}_2`.
436 n_electrons : Positive real
437 Default value None (half filling).
438 Number of electrons in the flake. An integer at T = 0.
439 temperature : Positive real or zero, default 0
441 Returns
442 -------
443 result : :class:`CornerCharges`
444 '''
445 _check_model(ks, 0)
446 error_handling.nk(n_cells, 2)
447 if isinstance(n_cells, int):
448 n_cells = (n_cells, n_cells)
449 n_states = int(np.prod(n_cells)) * ks.norb
450 if n_electrons is None:
451 n_electrons = n_states // 2
452 error_handling.electrons(n_electrons, n_states)
453 error_handling.positive_real_zero(temperature, 'temperature')
454 if temperature == 0:
455 error_handling.integer_electrons(n_electrons)
456 en, vec = LA.eigh(ks.finite_ham(n_cells))
457 f, _, mu = _occupy(en, en[:0], n_electrons, temperature)
458 density = (np.abs(vec) ** 2) @ f
459 background = n_electrons / n_states
460 positions = flake_positions(ks, n_cells)
461 # fractional coordinates along a1, a2, split at the middle of the sample
462 a = np.array(ks.lat.prim_vec, dtype='f8')
463 frac = np.linalg.solve(a.T, positions.T).T
464 mid = (np.array(n_cells) - 1) / 2 + np.mean(np.linalg.solve(a.T, ks.orbital_positions().T).T, axis=0)
465 high = frac > mid[None, :]
466 excess = density - background
467 charges = np.array([[excess[(high[:, 0] == bool(s1)) & (high[:, 1] == bool(s2))].sum()
468 for s2 in (0, 1)] for s1 in (0, 1)])
469 return CornerCharges(en, density, positions, background, float(mu), charges)