Coverage for tbkit/exceptional.py: 100%

1r""" 

2Band touchings of two-dimensional (non-Hermitian) Bloch bands: 

3**diabolical points** (DPs), Hermitian degeneracies such as graphene's 

4Dirac points, where two eigenvalues meet but their eigenvectors stay 

5distinct; and **exceptional points** (EPs), where the eigenvalues *and* 

6the eigenvectors coalesce, so that :math:`H(\mathbf{k})` stops being 

7diagonalizable (a Jordan block). 

8 

9The tools below tell them apart and compute their invariants. Each takes a 

10*model*: a 2D **KSpace** (the plane is the Brillouin zone, in the k 

11coordinates of *get_ham*), or any callable ``model(p)`` returning a square 

12matrix for a point ``p`` of a two-parameter plane. 

13 

14* **Eigenvalue vorticity** of a pair of bands around a closed loop 

15 :math:`\Gamma` (Shen, Zhen and Fu, Phys. Rev. Lett. 120, 146402 (2018)), 

16 

17 .. math:: 

18 

19 \nu_{mn}(\Gamma) = -\frac{1}{2\pi}\oint_\Gamma 

20 \nabla_{\mathbf{k}}\arg\left[E_m(\mathbf{k}) - E_n(\mathbf{k})\right]\cdot d\mathbf{k}\, . 

21 

22 Band labels are not single-valued around an EP (the two eigenvalues swap 

23 after one loop), so *vorticity* follows the eigenvalues by continuation 

24 along the loop (*track_eigenvalues*), never by sorting. A second-order EP 

25 has :math:`\nu = \pm1/2`, a DP :math:`\nu = 0`. 

26* **Discriminant winding** (Yang, Schnyder, Hu and Chiu, Phys. Rev. Lett. 

27 126, 086401 (2021)): the discriminant 

28 :math:`\Delta(\mathbf{k}) = \prod_{m<n}(E_m - E_n)^2` is a polynomial in the 

29 matrix elements, computed from the characteristic polynomial without any 

30 eigenvalue, and its winding 

31 :math:`W(\Gamma) = \frac{1}{2\pi}\oint_\Gamma \nabla_{\mathbf{k}}\arg\Delta\cdot d\mathbf{k}` 

32 counts the EPs inside :math:`\Gamma` with their charges. With the two sign 

33 conventions above, :math:`W = -2\sum_{m<n}\nu_{mn}`: :math:`W = -2\nu = \mp1` 

34 for a second-order EP of a two-band model. 

35* **EP finder** (*find_exceptional_points*): zeros of :math:`\Delta` with 

36 :math:`W \neq 0` on a plaquette mesh of the Brillouin zone, refined by 

37 Newton's method, with their charge and their order (the size of the Jordan 

38 block, from the Petermann factors, which diverge at an EP). Since 

39 :math:`\Delta` is periodic, the charges sum to zero over the zone: EPs come in 

40 pairs of opposite charge (the doubling theorem). When a symmetry makes 

41 :math:`\Delta` real, EPs form lines instead (exceptional rings, Zhen et al., 

42 Nature 525, 354 (2015)), which the finder also returns; *fermi_arcs* draws 

43 the bulk Fermi arcs :math:`\mathrm{Re}(E_m - E_n) = 0` that join a pair of 

44 EPs (Kozii and Fu, arXiv:1708.05841 (2017); Zhou et al., Science 359, 1009 

45 (2018)). 

46* **Encircling** (*encircle*): parallel transport of an eigenvector around 

47 a loop. Around a DP one loop returns the state with a Berry phase 

48 :math:`\pi`; around a second-order EP one loop swaps the two eigenstates, 

49 two loops return the state with a sign flip, and four loops restore it 

50 (Heiss, Phys. Rev. E 61, 929 (2000); Dembowski et al., Phys. Rev. Lett. 

51 86, 787 (2001)). 

52""" 

53from __future__ import annotations 

54 

55from typing import Callable 

56 

57import numpy as np 

58from numpy.typing import ArrayLike, NDArray 

59import scipy.linalg as LA 

60from scipy.optimize import linear_sum_assignment 

61 

62import tbkit.error_handling as error_handling 

63from tbkit.kspace import KSpace 

64 

65 

66PI = np.pi 

67_HALVINGS = 30 # a continuation step may be halved this many times 

68_ARG_STEP = PI / 4 # largest phase change of an eigenvalue difference per step 

69_MESH_SHIFT = (0.3183098861837907, 0.2718281828459045) # keeps mesh edges off symmetry lines 

70 

71 

72################################# 

73# MODELS AND LOOPS 

74################################# 

75 

76 

77def _model(model) -> Callable[[NDArray[np.float64]], NDArray[np.complex128]]: 

78 ''' 

79 Private function. The matrix-valued function p -> H(p) of a validated model. 

80 ''' 

81 error_handling.ham_model(model, KSpace) 

82 if isinstance(model, KSpace): 

83 error_handling.no_overlap(model._overlap_hop) 

84 return model.get_ham 

85 

86 def ham(p): 

87 h = np.asarray(model(np.asarray(p, dtype='f8'))) 

88 error_handling.ham_matrix(h) 

89 return h.astype('c16') 

90 return ham 

91 

92 

93def circle(center: ArrayLike, radius: float) -> Callable[[float], NDArray[np.float64]]: 

94 r''' 

95 Get a circular loop, run counterclockwise: 

96 :math:`\mathbf{k}(s) = \mathbf{k}_0 + r(\cos 2\pi s, \sin 2\pi s)`, 

97 :math:`s\in[0, 1]`. 

98 

99 :param center: Tuple/list/ndarray of 2 real numbers. :math:`\mathbf{k}_0`. 

100 :param radius: Positive real. :math:`r`. 

101 

102 :returns: 

103 * **loop** -- Callable of s, returning a real ndarray of shape (2,). 

104 ''' 

105 error_handling.k_vector(center, 'center', 2) 

106 error_handling.positive_real(radius, 'radius') 

107 k0 = np.asarray(center, dtype='f8') 

108 

109 def loop(s): 

110 return k0 + radius * np.array([np.cos(2 * PI * s), np.sin(2 * PI * s)]) 

111 return loop 

112 

113 

114def _loop(loop) -> Callable[[float], NDArray[np.float64]]: 

115 ''' 

116 Private function. The callable s -> k(s), s in [0, 1], of a validated loop 

117 (a callable, or the vertices of a polygon, each side covering an equal 

118 range of s). 

119 ''' 

120 error_handling.loop(loop) 

121 if callable(loop): 

122 def traced(s): 

123 p = np.asarray(loop(float(s)), dtype='f8') 

124 error_handling.loop_point(p) 

125 return p 

126 return traced 

127 vert = np.asarray(loop, dtype='f8') 

128 closed = np.vstack([vert, vert[:1]]) 

129 m = len(vert) 

130 

131 def path(s): 

132 x = (s % 1.) * m 

133 i = min(int(x), m - 1) 

134 return closed[i] + (x - i) * (closed[i + 1] - closed[i]) 

135 return path 

136 

137 

138def _closed(ham, path) -> int: 

139 ''' 

140 Private function. Check that the loop is closed (H at its end equals H at 

141 its start, which allows loops across the zone), and return norb. 

142 ''' 

143 h0 = ham(path(0.)) 

144 error_handling.closed_loop(h0, ham(path(1.))) 

145 return h0.shape[0] 

146 

147 

148################################# 

149# EIGENVALUE CONTINUATION 

150################################# 

151 

152 

153def _eig(h: NDArray[np.complex128], vectors: bool): 

154 ''' 

155 Private function. Eigenvalues, right and left eigenvectors (as kets: 

156 columns of vl with vl^H H = E vl^H) of the general eigenproblem. 

157 ''' 

158 if vectors: 

159 w, vl, vr = LA.eig(h, left=True, right=True) 

160 return w, vr, vl 

161 return LA.eigvals(h), None, None 

162 

163 

164def _match(pred, w_old, w_new, vr_old, vr_new): 

165 ''' 

166 Private function. Permutation of *w_new* continuing *w_old* (predicted 

167 at *pred*), or None if the step is ambiguous: a matched eigenvalue must 

168 be four times closer to its prediction than to any other eigenvalue, 

169 no eigenvalue difference may turn by more than pi/4, and (with vectors) 

170 each eigenvector must stay close to its predecessor. Eigenvalues equal to 

171 rounding precision form a cluster and do not make the step ambiguous. 

172 ''' 

173 _, perm = linear_sum_assignment(np.abs(pred[:, None] - w_new[None, :])) 

174 w = w_new[perm] 

175 deg = 1e-9 * max(1., float(np.max(np.abs(w)))) 

176 dist = np.abs(w[:, None] - w[None, :]) 

177 np.fill_diagonal(dist, np.inf) 

178 single = dist.min(axis=1) >= deg # no partner equal to rounding precision 

179 dist[dist < deg] = np.inf 

180 if np.any(np.abs(w - pred) >= 0.25 * dist.min(axis=1)): 

181 return None 

182 d_old = w_old[:, None] - w_old[None, :] 

183 d_new = w[:, None] - w[None, :] 

184 keep = (np.abs(d_old) > deg) & (np.abs(d_new) > deg) 

185 if np.any(np.abs(np.angle(d_new[keep] / d_old[keep])) > _ARG_STEP): 

186 return None 

187 if vr_old is not None: 

188 # within a degenerate cluster the eigenvectors are arbitrary 

189 overlap = np.abs(np.sum(vr_old.conj() * vr_new[:, perm], axis=0)) 

190 if np.any(overlap[single] < 0.9): 

191 return None 

192 return perm 

193 

194 

195def _track(ham, path, s_end: float, nk: int, vectors: bool = False): 

196 ''' 

197 Private function. Follow the eigenvalues of H(path(s)) by continuation 

198 from s = 0 to *s_end* (loops are evaluated at s modulo 1): steps of at 

199 most 1/nk, halved while a step is ambiguous (see *_match*), with a 

200 linear prediction from the previous step. The labels are those of the 

201 eigenvalues at s = 0 sorted by real part. 

202 

203 Returns s (n,), en (n, norb), and, with *vectors*, the right and left 

204 eigenvectors (n, norb, norb), column m following label m. 

205 ''' 

206 ds_max = 1. / nk 

207 ds_min = ds_max / 2 ** _HALVINGS 

208 w, vr, vl = _eig(ham(path(0.)), vectors) 

209 order = np.argsort(w.real, kind='stable') 

210 w = w[order] 

211 if vectors: 

212 vr, vl = vr[:, order], vl[:, order] 

213 s_list, w_list, vr_list, vl_list = [0.], [w], [vr], [vl] 

214 slope = np.zeros_like(w) 

215 s, ds = 0., ds_max 

216 while s < s_end: 

217 last = ds >= s_end - s 

218 s_new = s_end if last else s + ds 

219 step = s_new - s 

220 k = path(s_new % 1. if s_new != s_end or s_end % 1. else 1.) 

221 w_new, vr_new, vl_new = _eig(ham(k), vectors) 

222 perm = _match(w + slope * step, w, w_new, vr, vr_new) 

223 if perm is None: 

224 if step <= ds_min: 

225 error_handling.tracking(k) 

226 ds = step / 2 

227 continue 

228 slope = (w_new[perm] - w) / step 

229 w, s = w_new[perm], s_new 

230 if vectors: 

231 vr, vl = vr_new[:, perm], vl_new[:, perm] 

232 s_list.append(s) 

233 w_list.append(w) 

234 vr_list.append(vr) 

235 vl_list.append(vl) 

236 ds = min(2 * step, ds_max) 

237 s_arr, en = np.array(s_list), np.array(w_list) 

238 if vectors: 

239 return s_arr, en, np.array(vr_list), np.array(vl_list) 

240 return s_arr, en, None, None 

241 

242 

243def track_eigenvalues( 

244 model, loop, n_loops: int = 1, nk: int = 100, 

245) -> tuple[NDArray[np.float64], NDArray[np.float64], NDArray[np.complex128]]: 

246 r''' 

247 Follow the (complex) eigenvalues of :math:`H(\mathbf{k})` continuously 

248 around a closed loop, *n_loops* times. Each step (at most 1/nk of the 

249 loop) matches the new eigenvalues to a linear prediction from the 

250 previous ones, and is halved while that matching is ambiguous; nothing 

251 depends on how the eigensolver orders its output. Label :math:`m` is the 

252 :math:`m`-th eigenvalue at the start of the loop, sorted by real part. 

253 After one loop around a second-order exceptional point the labels come 

254 back swapped: ``en[-1]`` is then a permutation of ``en[0]``. 

255 

256 :param model: 2D **KSpace**, or callable of a point of a plane returning 

257 a square matrix. 

258 :param loop: Closed loop: a callable of :math:`s\in[0, 1]` returning a 

259 point (e.g. *circle*), or an array of polygon vertices, shape (m, 2). 

260 For a **KSpace**, points are k in the coordinates of *get_ham*, and a 

261 loop may close up to a reciprocal lattice vector. 

262 :param n_loops: Positive integer. Default value 1. Number of turns. 

263 :param nk: Positive integer. Default value 100. Steps per turn before 

264 refinement. 

265 

266 :returns: 

267 * **s** -- Real ndarray, shape (n,). Loop parameter, from 0 to *n_loops*. 

268 * **k** -- Real ndarray, shape (n, 2). Points of the loop. 

269 * **en** -- Complex ndarray, shape (n, norb). Eigenvalues, en[:, m] 

270 continuous in s. 

271 ''' 

272 ham = _model(model) 

273 path = _loop(loop) 

274 error_handling.positive_int(n_loops, 'n_loops') 

275 error_handling.positive_int(nk, 'nk') 

276 _closed(ham, path) 

277 s, en, _, _ = _track(ham, path, float(n_loops), nk) 

278 return s, np.array([path(x % 1.) for x in s]), en 

279 

280 

281def vorticity(model, loop, bands: tuple[int, int] = (0, 1), nk: int = 100) -> float: 

282 r''' 

283 Get the eigenvalue vorticity of a pair of bands around a closed loop 

284 :math:`\Gamma` (Shen, Zhen and Fu, Phys. Rev. Lett. 120, 146402 (2018)), 

285 

286 .. math:: 

287 

288 \nu_{mn}(\Gamma) = -\frac{1}{2\pi}\oint_\Gamma 

289 \nabla_{\mathbf{k}}\arg\left[E_m(\mathbf{k}) - E_n(\mathbf{k})\right] 

290 \cdot d\mathbf{k}\, , 

291 

292 with the eigenvalues followed by continuation (see *track_eigenvalues*). 

293 :math:`\nu_{mn} = \pm1/2` around a second-order exceptional point, where 

294 :math:`E_m - E_n \propto \sqrt{\mathbf{k}-\mathbf{k}_{EP}}` changes sign 

295 after one loop; :math:`\pm1/3` for each pair of a third-order one; 

296 0 around a diabolical point, or a loop enclosing no degeneracy. 

297 :math:`\nu_{mn} = \nu_{nm}`. 

298 

299 :param model: See *track_eigenvalues*. 

300 :param loop: See *track_eigenvalues*. Counterclockwise for the sign above. 

301 :param bands: Pair of distinct integers. Default value (0, 1). Band labels 

302 :math:`(m, n)` at the start of the loop, sorted by real part. 

303 :param nk: Positive integer. Default value 100. See *track_eigenvalues*. 

304 

305 :returns: 

306 * **nu** -- Real number, a multiple of 1/2 around second-order EPs. 

307 ''' 

308 ham = _model(model) 

309 path = _loop(loop) 

310 error_handling.positive_int(nk, 'nk') 

311 norb = _closed(ham, path) 

312 error_handling.band_pair(bands, norb) 

313 s, en, _, _ = _track(ham, path, 1., nk) 

314 diff = en[:, bands[0]] - en[:, bands[1]] 

315 if np.any(diff == 0.): 

316 error_handling.tracking(path(s[np.argmin(np.abs(diff))] % 1.)) 

317 return float(-np.sum(np.angle(diff[1:] / diff[:-1])) / (2 * PI)) 

318 

319 

320################################# 

321# DISCRIMINANT 

322################################# 

323 

324 

325def _charpoly(h: NDArray[np.complex128]) -> NDArray[np.complex128]: 

326 r''' 

327 Private function. Coefficients, highest degree first, of the monic 

328 characteristic polynomial :math:`p(z) = \det(z - H)`, from its values 

329 at N+1 points of a circle of radius :math:`\|H\|_F` (a discrete Fourier 

330 transform; no eigenvalues). 

331 ''' 

332 n = h.shape[0] 

333 r = LA.norm(h) or 1. 

334 z = r * np.exp(2j * PI * np.arange(n + 1) / (n + 1)) 

335 vals = np.array([LA.det(zz * np.eye(n) - h) for zz in z]) 

336 coef = np.fft.fft(vals) / (n + 1) / r ** np.arange(n + 1) 

337 coef[n] = 1. 

338 return coef[::-1] 

339 

340 

341def _discriminant(h: NDArray[np.complex128]) -> complex: 

342 r''' 

343 Private function. :math:`\Delta = (-1)^{N(N-1)/2}\,\mathrm{Res}(p, p')` 

344 of the monic characteristic polynomial, from the determinant of the 

345 Sylvester matrix, computed for :math:`H/\|H\|_F` and rescaled. 

346 ''' 

347 n = h.shape[0] 

348 if n < 2: 

349 return 1. + 0j 

350 r = LA.norm(h) 

351 if r == 0.: 

352 return 0j 

353 c = _charpoly(h / r) 

354 dc = c[:-1] * np.arange(n, 0, -1) 

355 syl = np.zeros((2 * n - 1, 2 * n - 1), 'c16') 

356 for i in range(n - 1): 

357 syl[i, i:i + n + 1] = c 

358 for i in range(n): 

359 syl[n - 1 + i, i:i + n] = dc 

360 return complex((-1) ** (n * (n - 1) // 2) * LA.det(syl) * r ** (n * (n - 1))) 

361 

362 

363def discriminant(model, k: ArrayLike) -> complex: 

364 r''' 

365 Get the discriminant of :math:`H(\mathbf{k})`, 

366 

367 .. math:: 

368 

369 \Delta(\mathbf{k}) = \prod_{m<n}\left[E_m(\mathbf{k}) - E_n(\mathbf{k})\right]^2 

370 = (-1)^{N(N-1)/2}\,\mathrm{Res}(p, p')\, , 

371 

372 computed as the resultant of the characteristic polynomial 

373 :math:`p(E) = \det[E - H(\mathbf{k})]` and its derivative (the determinant 

374 of their Sylvester matrix), without the eigenvalues. It is a polynomial 

375 in the matrix elements, hence smooth and single-valued, and it vanishes 

376 exactly where two eigenvalues meet. For a two-band model 

377 :math:`H = d_0 + \mathbf{d}\cdot\boldsymbol\sigma`, 

378 :math:`\Delta = 4\,\mathbf{d}\cdot\mathbf{d}`. 

379 

380 :param model: See *track_eigenvalues*. 

381 :param k: Tuple/list/ndarray of 2 real numbers. Point of the plane. 

382 

383 :returns: 

384 * **delta** -- Complex number, in units of energy^(N(N-1)) (1 if N = 1). 

385 ''' 

386 ham = _model(model) 

387 error_handling.k_vector(k, 'k', 2) 

388 return _discriminant(ham(np.asarray(k, dtype='f8'))) 

389 

390 

391def _arg_change(disc, path, s_a, s_b, d_a, d_b, depth=0) -> float: 

392 ''' 

393 Private function. Change of arg(Delta) from path(s_a) to path(s_b), 

394 bisecting until each piece turns by at most pi/4. 

395 ''' 

396 if d_a == 0. or d_b == 0. or depth > 40: 

397 error_handling.winding_path(path(s_b if d_b == 0. else s_a)) 

398 step = float(np.angle(d_b / d_a)) 

399 if abs(step) <= _ARG_STEP: 

400 return step 

401 s_m = (s_a + s_b) / 2 

402 d_m = disc(path(s_m)) 

403 return (_arg_change(disc, path, s_a, s_m, d_a, d_m, depth + 1) 

404 + _arg_change(disc, path, s_m, s_b, d_m, d_b, depth + 1)) 

405 

406 

407def _winding(ham, path, nk: int) -> float: 

408 ''' 

409 Private function. Winding of the discriminant along path(s), s in [0, 1]. 

410 ''' 

411 def disc(p): 

412 return _discriminant(ham(p)) 

413 s = np.linspace(0., 1., nk + 1) 

414 d = [disc(path(x)) for x in s] 

415 total = sum(_arg_change(disc, path, s[i], s[i + 1], d[i], d[i + 1]) for i in range(nk)) 

416 return total / (2 * PI) 

417 

418 

419def discriminant_winding(model, loop, nk: int = 100) -> float: 

420 r''' 

421 Get the winding number of the discriminant (see *discriminant*) around a 

422 closed loop :math:`\Gamma` (Yang, Schnyder, Hu and Chiu, Phys. Rev. Lett. 

423 126, 086401 (2021)), 

424 

425 .. math:: 

426 

427 W(\Gamma) = \frac{1}{2\pi}\oint_\Gamma \nabla_{\mathbf{k}}\arg\Delta(\mathbf{k}) 

428 \cdot d\mathbf{k} = -2\sum_{m<n}\nu_{mn}(\Gamma)\, , 

429 

430 the total charge of the exceptional points inside :math:`\Gamma` 

431 (:math:`\pm1` for a second-order EP, where :math:`\Delta` vanishes 

432 linearly). It needs no eigenvalue tracking. A diabolical point has 

433 :math:`W = 0`: :math:`\Delta \geq 0` for a Hermitian model. The relation to 

434 the vorticities *vorticity* holds with the signs as written, the minus 

435 sign coming from Shen, Zhen and Fu's definition of :math:`\nu`. 

436 

437 :param model: See *track_eigenvalues*. 

438 :param loop: See *track_eigenvalues*. Counterclockwise for the sign above. 

439 :param nk: Positive integer. Default value 100. Points on the loop; a 

440 step along which :math:`\arg\Delta` turns by more than :math:`\pi/4` 

441 is bisected. 

442 

443 :returns: 

444 * **W** -- Real number, close to an integer. 

445 ''' 

446 ham = _model(model) 

447 path = _loop(loop) 

448 error_handling.positive_int(nk, 'nk') 

449 _closed(ham, path) 

450 return float(_winding(ham, path, nk)) 

451 

452 

453################################# 

454# EXCEPTIONAL POINTS 

455################################# 

456 

457 

458def _petermann(h: NDArray[np.complex128]): 

459 ''' 

460 Private function. Eigenvalues (sorted by real part) and Petermann factors. 

461 ''' 

462 w, vl, vr = LA.eig(h, left=True, right=True) 

463 order = np.argsort(w.real, kind='stable') 

464 w, vl, vr = w[order], vl[:, order], vr[:, order] 

465 norms = np.sum(np.abs(vl) ** 2, axis=0) * np.sum(np.abs(vr) ** 2, axis=0) 

466 pair = np.abs(np.sum(vl.conj() * vr, axis=0)) ** 2 

467 return w, norms / np.maximum(pair, np.finfo(float).tiny) 

468 

469 

470def petermann_factors(model, k: ArrayLike) -> tuple[NDArray[np.complex128], NDArray[np.float64]]: 

471 r''' 

472 Get the eigenvalues of :math:`H(\mathbf{k})` and their Petermann factors 

473 

474 .. math:: 

475 

476 K_n = \frac{\langle L_n|L_n\rangle\langle R_n|R_n\rangle}{|\langle L_n|R_n\rangle|^2} 

477 \geq 1\, , 

478 

479 with :math:`|R_n\rangle`, :math:`|L_n\rangle` the right and left 

480 eigenvectors: 1 for a normal (e.g. Hermitian) matrix, diverging at an 

481 exceptional point, where :math:`\langle L_n|R_n\rangle \to 0` 

482 (as :math:`|\mathbf{k}-\mathbf{k}_{EP}|^{-1}` near a second-order one). 

483 

484 :param model: See *track_eigenvalues*. 

485 :param k: Tuple/list/ndarray of 2 real numbers. 

486 

487 :returns: 

488 * **en** -- Complex ndarray, shape (norb,). Eigenvalues, sorted by real part. 

489 * **K** -- Real ndarray, shape (norb,). Petermann factors. 

490 ''' 

491 ham = _model(model) 

492 error_handling.k_vector(k, 'k', 2) 

493 return _petermann(ham(np.asarray(k, dtype='f8'))) 

494 

495 

496def _order(h: NDArray[np.complex128], k_min: float = 1e4) -> tuple[int, float, complex]: 

497 ''' 

498 Private function. Order of the exceptional point at H: the number of 

499 eigenvalues whose Petermann factor exceeds sqrt(max K) (the coalescing 

500 ones, all of order max K), 1 if max K < k_min (no EP); with max K and 

501 the mean coalescing eigenvalue. 

502 ''' 

503 w, kp = _petermann(h) 

504 k_max = float(np.max(kp)) 

505 if k_max < k_min: 

506 return 1, k_max, complex(w[np.argmax(kp)]) 

507 big = kp > np.sqrt(k_max) 

508 return int(np.sum(big)), k_max, complex(np.mean(w[big])) 

509 

510 

511class ExceptionalPoints(): 

512 r''' 

513 Exceptional points of a 2D model, see *find_exceptional_points*. 

514 

515 :ivar k: Real ndarray, shape (n, 2). Isolated exceptional points. 

516 :ivar charge: Integer ndarray, shape (n,). Discriminant winding around each 

517 one (see *discriminant_winding*): :math:`\pm1` for a second-order EP. 

518 :ivar order: Integer ndarray, shape (n,). Number of coalescing eigenvalues 

519 (the size of the Jordan block): 2, 3, ... 

520 :ivar petermann: Real ndarray, shape (n,). Largest Petermann factor at 

521 each point (divergent at an EP: about :math:`10^{8}` or more once refined). 

522 :ivar energy: Complex ndarray, shape (n,). The coalescing eigenvalue. 

523 :ivar lines: Real ndarray, shape (m, 2). Points on exceptional lines 

524 (rings), when a symmetry makes the discriminant real; empty otherwise. 

525 ''' 

526 

527 def __init__(self, k, charge, order, petermann, energy, lines) -> None: 

528 self.k = k 

529 self.charge = charge 

530 self.order = order 

531 self.petermann = petermann 

532 self.energy = energy 

533 self.lines = lines 

534 

535 @property 

536 def total_charge(self) -> int: 

537 ''' 

538 Sum of the charges: zero over a whole Brillouin zone (doubling theorem). 

539 ''' 

540 return int(np.sum(self.charge)) 

541 

542 

543def _mesh(ks: KSpace, nk: int): 

544 ''' 

545 Private function. Shifted nk x nk mesh of the zone: node k-points, 

546 shape (nk + 1, nk + 1, 2) (the last row/column repeat the first ones 

547 shifted by a reciprocal vector), and the plaquette size. 

548 ''' 

549 fr = (np.arange(nk + 1)[:, None] + np.array(_MESH_SHIFT)[None, :]) / nk 

550 fracs = np.stack(np.meshgrid(fr[:, 0], fr[:, 1], indexing='ij'), axis=-1) 

551 size = min(np.linalg.norm(b) for b in ks.rec_vec_k) / nk 

552 return fracs @ ks.rec_vec_k, size 

553 

554 

555def _edges(nk: int): 

556 ''' 

557 Private function. Mesh edges ((i, j), (i', j')) along b1 then along b2. 

558 ''' 

559 return ([((i, j), (i + 1, j)) for i in range(nk) for j in range(nk)] 

560 + [((i, j), (i, j + 1)) for i in range(nk) for j in range(nk)]) 

561 

562 

563def _bisect(func, k_a, k_b, f_a, n_iter=60): 

564 ''' 

565 Private function. Root of the real function *func* on the segment 

566 [k_a, k_b], across which it changes sign (f_a = func(k_a)). 

567 ''' 

568 lo, hi = 0., 1. 

569 for _ in range(n_iter): 

570 mid = (lo + hi) / 2 

571 f_m = func(k_a + mid * (k_b - k_a)) 

572 if np.sign(f_m) == np.sign(f_a): 

573 lo, f_a = mid, f_m 

574 else: 

575 hi = mid 

576 return k_a + (lo + hi) / 2 * (k_b - k_a) 

577 

578 

579def _newton(disc, k0, mult, h, tol, max_iter): 

580 r''' 

581 Private function. Newton's method on (Re Delta, Im Delta) = 0, with a 

582 central-difference Jacobian; the step is multiplied by the multiplicity 

583 *mult* of the zero (Delta ~ L(k)^mult, L complex linear, as at a 

584 higher-order exceptional point). 

585 ''' 

586 k = np.array(k0, dtype='f8') 

587 for _ in range(max_iter): 

588 f = disc(k) 

589 jac = np.zeros((2, 2)) 

590 for mu in range(2): 

591 e = np.zeros(2) 

592 e[mu] = h 

593 d = (disc(k + e) - disc(k - e)) / (2 * h) 

594 jac[:, mu] = d.real, d.imag 

595 step = mult * np.linalg.lstsq(jac, np.array([f.real, f.imag]), rcond=None)[0] 

596 k = k - step 

597 if np.linalg.norm(step) < tol: 

598 break 

599 return k 

600 

601 

602def find_exceptional_points( 

603 ks: KSpace, nk: int = 40, tol: float = 1e-12, max_iter: int = 50, 

604) -> ExceptionalPoints: 

605 r''' 

606 Locate the exceptional points of a 2D model over its Brillouin zone. 

607 

608 1. The discriminant :math:`\Delta(\mathbf{k})` (see *discriminant*) is 

609 computed on an :math:`nk\times nk` mesh of the zone (shifted off the 

610 high-symmetry lines), and its winding around every plaquette from the 

611 phase changes along the mesh edges (bisected where :math:`\arg\Delta` 

612 turns by more than :math:`\pi/4`). A plaquette with a nonzero winding 

613 :math:`w` holds an EP. 

614 2. Each one is refined by Newton's method on 

615 :math:`(\mathrm{Re}\,\Delta, \mathrm{Im}\,\Delta) = 0` from the plaquette 

616 centre (the step multiplied by :math:`|w|`, the multiplicity of the zero). 

617 3. Its charge is the winding of :math:`\Delta` around a small circle 

618 (a quarter of a plaquette), and its order the number of eigenvalues 

619 whose Petermann factor diverges there (see *petermann_factors*). 

620 

621 If a symmetry makes :math:`\Delta` real everywhere (up to a constant 

622 phase), EPs are not isolated: they form lines, exceptional rings, where 

623 :math:`\Delta` changes sign. The lines are then sampled where they cross 

624 the mesh edges (bisection to machine precision), and no isolated point is 

625 returned. A Hermitian model has :math:`\Delta \geq 0`, and no EP. 

626 

627 Since :math:`\Delta` is periodic, the charges add up to zero 

628 (*total_charge*): EPs come in pairs of opposite charge. The mesh must be 

629 fine enough that no plaquette holds two EPs of opposite charge. 

630 

631 :param ks: **KSpace** instance with two primitive vectors. 

632 :param nk: Positive integer. Default value 40. Mesh points along each 

633 reciprocal vector. 

634 :param tol: Positive real. Default value 1e-12. Newton's method stops 

635 when a step is shorter than *tol* (units of k). 

636 :param max_iter: Positive integer. Default value 50. Newton iterations. 

637 

638 :returns: 

639 * **eps** -- **ExceptionalPoints** instance (positions *k*, *charge*, 

640 *order*, *petermann*, *energy*, and the exceptional *lines*). 

641 ''' 

642 error_handling.kspace(ks, KSpace) 

643 error_handling.dim_exact(ks.dim, 2) 

644 ham = _model(ks) 

645 error_handling.positive_int(nk, 'nk') 

646 error_handling.positive_real(tol, 'tol') 

647 error_handling.positive_int(max_iter, 'max_iter') 

648 

649 def disc(k): 

650 return _discriminant(ham(k)) 

651 

652 kmesh, size = _mesh(ks, nk) 

653 d = np.array([[disc(kmesh[i, j]) for j in range(nk)] for i in range(nk)]) 

654 d = np.pad(d, ((0, 1), (0, 1)), mode='wrap') 

655 scale = float(np.max(np.abs(d))) 

656 error_handling.discriminant_nonzero(scale) 

657 phase = np.exp(-0.5j * np.angle(np.sum(d ** 2))) 

658 empty = np.zeros((0, 2)) 

659 if np.max(np.abs((phase * d).imag)) <= 1e-8 * scale: 

660 # a real discriminant: exceptional lines where it changes sign 

661 lines = [_bisect(lambda k: (phase * disc(k)).real, kmesh[a], kmesh[b], (phase * d[a]).real) 

662 for a, b in _edges(nk) if (phase * d[a]).real * (phase * d[b]).real < 0.] 

663 return ExceptionalPoints(empty, np.zeros(0, int), np.zeros(0, int), np.zeros(0), 

664 np.zeros(0, 'c16'), np.array(lines).reshape(-1, 2)) 

665 # phase change along the edge leaving node (i, j) along b1 (axis 0) or b2 

666 # (axis 1); Delta is periodic, so the last row/column repeat the first ones 

667 edge = np.zeros((2, nk, nk)) 

668 for a, b in _edges(nk): 

669 edge[int(b[1] != a[1])][a] = _arg_change( 

670 disc, lambda s, a=a, b=b: kmesh[a] + s * (kmesh[b] - kmesh[a]), 0., 1., d[a], d[b]) 

671 b1, b2 = ks.rec_vec_k 

672 orient = np.sign(b1[0] * b2[1] - b1[1] * b2[0]) 

673 points: list = [] 

674 for i in range(nk): 

675 for j in range(nk): 

676 w = (edge[0, i, j] + edge[1, (i + 1) % nk, j] 

677 - edge[0, i, (j + 1) % nk] - edge[1, i, j]) / (2 * PI) 

678 w = int(np.rint(orient * w)) 

679 if w == 0: 

680 continue 

681 center = (kmesh[i, j] + kmesh[i + 1, j + 1]) / 2 

682 k_ep = _newton(disc, center, abs(w), 1e-5 * size, tol, max_iter) 

683 if all(np.linalg.norm(k_ep - p) > 1e-6 * size for p in points): 

684 points.append(k_ep) 

685 charge, order, peter, energy = [], [], [], [] 

686 for p in points: 

687 charge.append(int(np.rint(_winding(ham, circle(p, 0.25 * size), 16)))) 

688 o, kp, e = _order(ham(p)) 

689 order.append(o) 

690 peter.append(kp) 

691 energy.append(e) 

692 return ExceptionalPoints(np.array(points).reshape(-1, 2), np.array(charge, int), 

693 np.array(order, int), np.array(peter), 

694 np.array(energy, 'c16'), empty) 

695 

696 

697def fermi_arcs(ks: KSpace, nk: int = 60) -> NDArray[np.float64]: 

698 r''' 

699 Get the bulk Fermi arcs of a 2D non-Hermitian model: the curves where 

700 two eigenvalues have equal real parts but different imaginary parts, 

701 :math:`\mathrm{Re}(E_m - E_n) = 0`, :math:`\mathrm{Im}(E_m - E_n) \neq 0`. 

702 Such an arc joins two exceptional points of opposite charge (for a 

703 two-band model it is where :math:`\Delta` is real and negative, see 

704 *discriminant*); it is what a Hermitian Dirac point turns into when a 

705 non-Hermitian term splits it into an EP pair (Kozii and Fu 2017; Zhou et 

706 al., Science 359, 1009 (2018)). 

707 

708 Along every edge of an :math:`nk\times nk` mesh of the zone the 

709 eigenvalues are followed by continuation (see *track_eigenvalues*), and 

710 where :math:`\mathrm{Re}(E_m - E_n)` changes sign, the crossing is found by 

711 bisection. Regions where the real parts coincide over an area (inside an 

712 exceptional ring) show no sign change and are not reported, and neither 

713 are edges that cross an exceptional point or line, along which the 

714 eigenvalues cannot be followed. 

715 

716 :param ks: **KSpace** instance with two primitive vectors. 

717 :param nk: Positive integer. Default value 60. Mesh points along each 

718 reciprocal vector. 

719 

720 :returns: 

721 * **points** -- Real ndarray, shape (m, 2). Points on the arcs. 

722 ''' 

723 error_handling.kspace(ks, KSpace) 

724 error_handling.dim_exact(ks.dim, 2) 

725 ham = _model(ks) 

726 error_handling.positive_int(nk, 'nk') 

727 kmesh, _ = _mesh(ks, nk) 

728 points = [] 

729 for a, b in _edges(nk): 

730 k_a, k_b = kmesh[a], kmesh[b] 

731 try: 

732 s, en, _, _ = _track(ham, lambda x: k_a + x * (k_b - k_a), 1., 1) 

733 except ValueError: 

734 continue # the edge crosses an exceptional point (or line) 

735 noise = 1e-10 * max(1., float(np.max(np.abs(en)))) 

736 norb = en.shape[1] 

737 for m in range(norb): 

738 for n in range(m + 1, norb): 

739 g = (en[:, m] - en[:, n]).real 

740 for i in np.nonzero((g[:-1] * g[1:] < 0.) & (np.abs(g[:-1]) > noise) 

741 & (np.abs(g[1:]) > noise))[0]: