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 Parameters 

100 ---------- 

101 center : Tuple/list/ndarray of 2 real numbers 

102 :math:`\mathbf{k}_0`. 

103 radius : Positive real 

104 :math:`r`. 

105 

106 Returns 

107 ------- 

108 loop : Callable of s, returning a real ndarray of shape (2,) 

109 ''' 

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

111 error_handling.positive_real(radius, 'radius') 

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

113 

114 def loop(s): 

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

116 return loop 

117 

118 

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

120 ''' 

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

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

123 range of s). 

124 ''' 

125 error_handling.loop(loop) 

126 if callable(loop): 

127 def traced(s): 

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

129 error_handling.loop_point(p) 

130 return p 

131 return traced 

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

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

134 m = len(vert) 

135 

136 def path(s): 

137 x = (s % 1.) * m 

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

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

140 return path 

141 

142 

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

144 ''' 

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

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

147 ''' 

148 h0 = ham(path(0.)) 

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

150 return h0.shape[0] 

151 

152 

153################################# 

154# EIGENVALUE CONTINUATION 

155################################# 

156 

157 

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

159 ''' 

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

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

162 ''' 

163 if vectors: 

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

165 return w, vr, vl 

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

167 

168 

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

170 ''' 

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

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

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

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

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

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

177 ''' 

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

179 w = w_new[perm] 

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

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

182 np.fill_diagonal(dist, np.inf) 

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

184 dist[dist < deg] = np.inf 

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

186 return None 

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

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

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

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

191 return None 

192 if vr_old is not None: 

193 # within a degenerate cluster the eigenvectors are arbitrary 

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

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

196 return None 

197 return perm 

198 

199 

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

201 ''' 

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

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

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

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

206 eigenvalues at s = 0 sorted by real part. 

207 

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

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

210 ''' 

211 ds_max = 1. / nk 

212 ds_min = ds_max / 2 ** _HALVINGS 

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

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

215 w = w[order] 

216 if vectors: 

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

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

219 slope = np.zeros_like(w) 

220 s, ds = 0., ds_max 

221 while s < s_end: 

222 last = ds >= s_end - s 

223 s_new = s_end if last else s + ds 

224 step = s_new - s 

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

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

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

228 if perm is None: 

229 if step <= ds_min: 

230 error_handling.tracking(k) 

231 ds = step / 2 

232 continue 

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

234 w, s = w_new[perm], s_new 

235 if vectors: 

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

237 s_list.append(s) 

238 w_list.append(w) 

239 vr_list.append(vr) 

240 vl_list.append(vl) 

241 ds = min(2 * step, ds_max) 

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

243 if vectors: 

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

245 return s_arr, en, None, None 

246 

247 

248def track_eigenvalues( 

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

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

251 r''' 

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

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

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

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

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

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

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

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

260 

261 Parameters 

262 ---------- 

263 model 

264 2D **KSpace**, or callable of a point of a plane returning 

265 a square matrix. 

266 loop 

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

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

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

270 loop may close up to a reciprocal lattice vector. 

271 n_loops : Positive integer, default 1 

272 Number of turns. 

273 nk : Positive integer, default 100 

274 Steps per turn before 

275 refinement. 

276 

277 Returns 

278 ------- 

279 s : Real ndarray, shape (n,) 

280 Loop parameter, from 0 to *n_loops*. 

281 k : Real ndarray, shape (n, 2) 

282 Points of the loop. 

283 en : Complex ndarray, shape (n, norb) 

284 Eigenvalues, en[:, m] 

285 continuous in s. 

286 ''' 

287 ham = _model(model) 

288 path = _loop(loop) 

289 error_handling.positive_int(n_loops, 'n_loops') 

290 error_handling.positive_int(nk, 'nk') 

291 _closed(ham, path) 

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

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

294 

295 

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

297 r''' 

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

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

300 

301 .. math:: 

302 

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

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

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

306 

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

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

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

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

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

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

313 

314 Parameters 

315 ---------- 

316 model 

317 See *track_eigenvalues*. 

318 loop 

319 See *track_eigenvalues*. Counterclockwise for the sign above. 

320 bands : Pair of distinct integers 

321 Default value (0, 1). Band labels 

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

323 nk : Positive integer, default 100 

324 See *track_eigenvalues*. 

325 

326 Returns 

327 ------- 

328 nu : Real number, a multiple of 1/2 around second-order EPs 

329 ''' 

330 ham = _model(model) 

331 path = _loop(loop) 

332 error_handling.positive_int(nk, 'nk') 

333 norb = _closed(ham, path) 

334 error_handling.band_pair(bands, norb) 

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

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

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

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

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

340 

341 

342################################# 

343# DISCRIMINANT 

344################################# 

345 

346 

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

348 r''' 

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

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

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

352 transform; no eigenvalues). 

353 ''' 

354 n = h.shape[0] 

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

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

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

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

359 coef[n] = 1. 

360 return coef[::-1] 

361 

362 

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

364 r''' 

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

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

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

368 ''' 

369 n = h.shape[0] 

370 if n < 2: 

371 return 1. + 0j 

372 r = LA.norm(h) 

373 if r == 0.: 

374 return 0j 

375 c = _charpoly(h / r) 

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

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

378 for i in range(n - 1): 

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

380 for i in range(n): 

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

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

383 

384 

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

386 r''' 

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

388 

389 .. math:: 

390 

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

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

393 

394 computed as the resultant of the characteristic polynomial 

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

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

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

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

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

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

401 

402 Parameters 

403 ---------- 

404 model 

405 See *track_eigenvalues*. 

406 k : Tuple/list/ndarray of 2 real numbers 

407 Point of the plane. 

408 

409 Returns 

410 ------- 

411 delta : Complex number, in units of energy^(N(N-1)) (1 if N = 1) 

412 ''' 

413 ham = _model(model) 

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

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

416 

417 

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

419 ''' 

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

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

422 ''' 

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

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

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

426 if abs(step) <= _ARG_STEP: 

427 return step 

428 s_m = (s_a + s_b) / 2 

429 d_m = disc(path(s_m)) 

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

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

432 

433 

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

435 ''' 

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

437 ''' 

438 def disc(p): 

439 return _discriminant(ham(p)) 

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

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

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

443 return total / (2 * PI) 

444 

445 

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

447 r''' 

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

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

450 126, 086401 (2021)), 

451 

452 .. math:: 

453 

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

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

456 

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

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

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

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

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

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

463 

464 Parameters 

465 ---------- 

466 model 

467 See *track_eigenvalues*. 

468 loop 

469 See *track_eigenvalues*. Counterclockwise for the sign above. 

470 nk : Positive integer, default 100 

471 Points on the loop; a 

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

473 is bisected. 

474 

475 Returns 

476 ------- 

477 W : Real number, close to an integer 

478 ''' 

479 ham = _model(model) 

480 path = _loop(loop) 

481 error_handling.positive_int(nk, 'nk') 

482 _closed(ham, path) 

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

484 

485 

486################################# 

487# EXCEPTIONAL POINTS 

488################################# 

489 

490 

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

492 ''' 

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

494 ''' 

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

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

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

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

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

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

501 

502 

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

504 r''' 

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

506 

507 .. math:: 

508 

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

510 \geq 1\, , 

511 

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

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

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

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

516 

517 Parameters 

518 ---------- 

519 model 

520 See *track_eigenvalues*. 

521 k : Tuple/list/ndarray of 2 real numbers 

522 

523 Returns 

524 ------- 

525 en : Complex ndarray, shape (norb,) 

526 Eigenvalues, sorted by real part. 

527 K : Real ndarray, shape (norb,) 

528 Petermann factors. 

529 ''' 

530 ham = _model(model) 

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

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

533 

534 

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

536 ''' 

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

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

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

540 the mean coalescing eigenvalue. 

541 ''' 

542 w, kp = _petermann(h) 

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

544 if k_max < k_min: 

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

546 big = kp > np.sqrt(k_max) 

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

548 

549 

550class ExceptionalPoints(): 

551 r''' 

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

553 

554 Attributes 

555 ---------- 

556 k : Real ndarray, shape (n, 2) 

557 Isolated exceptional points. 

558 charge : Integer ndarray, shape (n,) 

559 Discriminant winding around each 

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

561 order : Integer ndarray, shape (n,) 

562 Number of coalescing eigenvalues 

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

564 petermann : Real ndarray, shape (n,) 

565 Largest Petermann factor at 

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

567 energy : Complex ndarray, shape (n,) 

568 The coalescing eigenvalue. 

569 lines : Real ndarray, shape (m, 2) 

570 Points on exceptional lines 

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

572 ''' 

573 

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

575 self.k = k 

576 self.charge = charge 

577 self.order = order 

578 self.petermann = petermann 

579 self.energy = energy 

580 self.lines = lines 

581 

582 @property 

583 def total_charge(self) -> int: 

584 ''' 

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

586 ''' 

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

588 

589 

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

591 ''' 

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

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

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

595 ''' 

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

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

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

599 return fracs @ ks.rec_vec_k, size 

600 

601 

602def _edges(nk: int): 

603 ''' 

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

605 ''' 

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

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

608 

609 

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

611 ''' 

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

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

614 ''' 

615 lo, hi = 0., 1. 

616 for _ in range(n_iter): 

617 mid = (lo + hi) / 2 

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

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

620 lo, f_a = mid, f_m 

621 else: 

622 hi = mid 

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

624 

625 

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

627 r''' 

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

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

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

631 higher-order exceptional point). 

632 ''' 

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

634 for _ in range(max_iter): 

635 f = disc(k) 

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

637 for mu in range(2): 

638 e = np.zeros(2) 

639 e[mu] = h 

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

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

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

643 k = k - step 

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

645 break 

646 return k 

647 

648 

649def find_exceptional_points( 

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

651) -> ExceptionalPoints: 

652 r''' 

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

654 

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

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

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

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

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

660 :math:`w` holds an EP. 

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

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

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

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

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

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

667 

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

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

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

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

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

673 

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

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

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

677 

678 Parameters 

679 ---------- 

680 ks : **KSpace** instance with two primitive vectors 

681 nk : Positive integer, default 40 

682 Mesh points along each 

683 reciprocal vector. 

684 tol : Positive real, default 1e-12 

685 Newton's method stops 

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

687 max_iter : Positive integer, default 50 

688 Newton iterations. 

689 

690 Returns 

691 ------- 

692 eps 

693 **ExceptionalPoints** instance (positions *k*, *charge*, 

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

695 ''' 

696 error_handling.kspace(ks, KSpace) 

697 error_handling.dim_exact(ks.dim, 2) 

698 ham = _model(ks) 

699 error_handling.positive_int(nk, 'nk') 

700 error_handling.positive_real(tol, 'tol') 

701 error_handling.positive_int(max_iter, 'max_iter') 

702 

703 def disc(k): 

704 return _discriminant(ham(k)) 

705 

706 kmesh, size = _mesh(ks, nk) 

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

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

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

710 error_handling.discriminant_nonzero(scale) 

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

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

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

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

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

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

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

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

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

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

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

722 for a, b in _edges(nk): 

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

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

725 b1, b2 = ks.rec_vec_k 

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

727 points: list = [] 

728 for i in range(nk): 

729 for j in range(nk): 

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

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

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

733 if w == 0: 

734 continue 

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

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

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

738 points.append(k_ep) 

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

740 for p in points: 

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

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

743 order.append(o) 

744 peter.append(kp) 

745 energy.append(e) 

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

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

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

749 

750 

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

752 r''' 

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

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

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

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

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

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

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

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

761 

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

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

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

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

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

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

768 eigenvalues cannot be followed. 

769 

770 Parameters 

771 ---------- 

772 ks : **KSpace** instance with two primitive vectors 

773 nk : Positive integer, default 60 

774 Mesh points along each 

775 reciprocal vector.