Coverage for tbkit/floquet.py: 100%

1r""" 

2Floquet theory of periodically driven Tight-Binding models, 

3:math:`H(t + T) = H(t)`, :math:`\omega = 2\pi/T` (:math:`\hbar = 1`). 

4 

5The evolution over one period, :math:`U(T) = \mathcal{T}e^{-i\int_0^TH(t)dt} 

6= e^{-iH_FT}`, defines the Floquet Hamiltonian :math:`H_F`, whose 

7eigenvalues -- the quasienergies :math:`\epsilon`, defined modulo 

8:math:`\omega` -- replace the energies of a static model: the driven 

9system, observed once per period, behaves as if it were static with 

10Hamiltonian :math:`H_F`. Driving can therefore *engineer* band structures 

11and their topology: circularly polarized light turns graphene into a 

12Chern insulator (Oka and Aoki, Phys. Rev. B 79, 081406 (2009); Kitagawa, 

13Oka, Brataas, Fu and Demler, Phys. Rev. B 84, 235108 (2011)). 

14 

15Two routes give the quasienergies, checked against each other in the tests: 

16 

17* the time-ordered product of short-time propagators (*evolution_operator*, 

18 *quasienergies*, *effective_hamiltonian*); 

19* the Sambe (extended) space of the harmonics 

20 :math:`H(t) = \sum_m H_m e^{-im\omega t}`, a static eigenvalue problem 

21 (*harmonics*, *sambe_hamiltonian*). 

22 

23:class:`FloquetKSpace` makes a driven **KSpace** model: its *get_ham* is 

24the Floquet Hamiltonian :math:`H_F(\mathbf{k})`, so its bands, Berry 

25curvature and Chern numbers are the driven ones. 

26 

27Anomalous Floquet phases. The quasienergies are defined modulo 

28:math:`\omega`, so :math:`H_F` depends on where the branch cut of the 

29logarithm is put (the *epsilon* arguments). The Chern numbers of the 

30Floquet bands no longer count the edge states: every band can have 

31:math:`C = 0` while chiral edge states cross every gap (Kitagawa, Berg, 

32Rudner and Demler, Phys. Rev. B 82, 235114 (2010)). The winding number of 

33Rudner, Lindner, Berg and Levin (Phys. Rev. X 3, 031005 (2013)), which 

34uses the evolution at *all* times of the period, does 

35(*DrivenKSpace.winding_number*): 

36 

37* :class:`DrivenKSpace` makes a driven Bloch model out of any 

38 time-dependent :math:`H(\mathbf{k}, t)`; 

39* *step_drive* builds a piecewise-constant drive from a list of models 

40 (**KSpace** models, ribbons made by *kspace.ribbon*, real-space 

41 **System** instances, or matrices), evolved exactly as a product of 

42 step exponentials (:class:`StepDrive` in real space); 

43* *DrivenKSpace.edge_state_count* counts the chiral edge states of a 

44 driven ribbon, independently of the bulk winding number. 

45""" 

46from __future__ import annotations 

47 

48import copy 

49from typing import TYPE_CHECKING, Callable 

50 

51import numpy as np 

52from numpy.typing import ArrayLike, NDArray 

53import scipy.linalg as LA 

54 

55import tbkit.error_handling as error_handling 

56 

57if TYPE_CHECKING: 

58 from tbkit.kspace import KSpace 

59 from tbkit.lattice import Lattice 

60 

61 

62def evolution_operator(ham_t: Callable, period: float, n_steps: int = 100) -> NDArray[np.complex128]: 

63 r''' 

64 Get the one-period evolution operator 

65 :math:`U(T) = \prod_m e^{-iH(t_m)\,\delta t}` (later times on the left, 

66 midpoints :math:`t_m = (m + 1/2)\delta t`, error :math:`O(\delta t^2)`). 

67 

68 :param ham_t: Callable. ``ham_t(t)`` returns the Hamiltonian at time t (a square matrix). 

69 :param period: Positive real. :math:`T`. 

70 :param n_steps: Positive integer. Default value 100. Time steps per period. 

71 

72 :returns: 

73 * **U** -- Complex ndarray, the shape of the Hamiltonian. 

74 ''' 

75 error_handling.is_callable(ham_t, 'ham_t') 

76 error_handling.positive_real(period, 'period') 

77 error_handling.positive_int(n_steps, 'n_steps') 

78 dt = period / n_steps 

79 u = LA.expm(-1j * dt * np.asarray(ham_t(0.5 * dt), dtype='c16')) 

80 for m in range(1, n_steps): 

81 h = np.asarray(ham_t((m + 0.5) * dt), dtype='c16') 

82 u = LA.expm(-1j * dt * h) @ u 

83 return u 

84 

85 

86def _fold(eps: NDArray[np.float64], period: float, epsilon: float) -> NDArray[np.float64]: 

87 ''' 

88 Private function. Quasienergies folded into [epsilon - 2 pi/T, epsilon). 

89 ''' 

90 omega = 2 * np.pi / period 

91 rest = np.mod(eps - epsilon, omega) 

92 rest = np.where(rest >= omega, rest - omega, rest) # np.mod may round up to omega 

93 return epsilon - omega + rest 

94 

95 

96def _floquet( 

97 u: NDArray[np.complex128], period: float, epsilon: float | None = None, 

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

99 ''' 

100 Private function. Quasienergies (sorted, in [-pi/T, pi/T), or in 

101 [epsilon - 2 pi/T, epsilon)) and Floquet states of a unitary 

102 U = exp(-i H_F T), from its complex Schur form. 

103 ''' 

104 tri, z = LA.schur(u, output='complex') 

105 eps = -np.angle(np.diag(tri)) / period 

106 if epsilon is not None: 

107 eps = _fold(eps, period, epsilon) 

108 order = np.argsort(eps, kind='stable') 

109 return eps[order], z[:, order] 

110 

111 

112def _heff(u: NDArray[np.complex128], period: float, epsilon: float | None = None) -> NDArray[np.complex128]: 

113 ''' 

114 Private function. The Floquet Hamiltonian (i/T) ln U, branch cut at epsilon. 

115 ''' 

116 eps, z = _floquet(u, period, epsilon) 

117 h_f = z @ np.diag(eps) @ z.conj().T 

118 return (h_f + h_f.conj().T) / 2 

119 

120 

121def quasienergies( 

122 ham_t: Callable, period: float, n_steps: int = 100, epsilon: float | None = None, 

123) -> NDArray[np.float64]: 

124 r''' 

125 Get the quasienergies :math:`\epsilon`, eigenvalues of the Floquet 

126 Hamiltonian, in the first Floquet zone :math:`[-\omega/2, \omega/2)`, 

127 or in :math:`[\epsilon_{cut} - \omega, \epsilon_{cut})`. 

128 

129 :param ham_t: Callable. See *evolution_operator*. 

130 :param period: Positive real. 

131 :param n_steps: Positive integer. Default value 100. 

132 :param epsilon: Real number. Default value None (the zone 

133 :math:`[-\omega/2, \omega/2)`). Branch cut, see *effective_hamiltonian*. 

134 

135 :returns: 

136 * **eps** -- Real ndarray, sorted. 

137 ''' 

138 error_handling.is_callable(ham_t, 'ham_t') 

139 error_handling.positive_real(period, 'period') 

140 error_handling.positive_int(n_steps, 'n_steps') 

141 error_handling.branch_cut(epsilon) 

142 return _floquet(evolution_operator(ham_t, period, n_steps), period, epsilon)[0] 

143 

144 

145def effective_hamiltonian( 

146 ham_t: Callable, period: float, n_steps: int = 100, epsilon: float | None = None, 

147) -> NDArray[np.complex128]: 

148 r''' 

149 Get the Floquet Hamiltonian :math:`H_F = \frac{i}{T}\ln U(T)`, with its 

150 quasienergies in the first Floquet zone (Hermitian for a Hermitian drive). 

151 

152 The logarithm is defined up to multiples of :math:`\omega` on each 

153 eigenvalue: *epsilon* puts its branch cut at the quasienergy 

154 :math:`\epsilon_{cut}` (the phase :math:`e^{-i\epsilon_{cut}T}` of 

155 :math:`U`), so that the quasienergies lie in 

156 :math:`[\epsilon_{cut} - \omega, \epsilon_{cut})`. The default, the zone 

157 :math:`[-\omega/2, \omega/2)`, is the cut :math:`\epsilon_{cut} = \omega/2`. 

158 Rudner et al. (2013) take :math:`(\epsilon_{cut}, \epsilon_{cut} + \omega]` 

159 instead: :math:`H_F` shifted by :math:`\omega` (the same projectors, 

160 so the same Chern numbers and winding numbers). 

161 

162 :param ham_t: Callable. See *evolution_operator*. 

163 :param period: Positive real. 

164 :param n_steps: Positive integer. Default value 100. 

165 :param epsilon: Real number. Default value None (:math:`\epsilon_{cut} = \omega/2`, 

166 computed exactly as before this option existed). Branch cut, 

167 best placed in a quasienergy gap. 

168 

169 :returns: 

170 * **H_F** -- Complex ndarray. 

171 ''' 

172 error_handling.is_callable(ham_t, 'ham_t') 

173 error_handling.positive_real(period, 'period') 

174 error_handling.positive_int(n_steps, 'n_steps') 

175 error_handling.branch_cut(epsilon) 

176 return _heff(evolution_operator(ham_t, period, n_steps), period, epsilon) 

177 

178 

179def harmonics(ham_t: Callable, period: float, n_harmonics: int, n_samples: int | None = None) -> dict: 

180 r''' 

181 Get the Fourier components :math:`H_m = \frac1T\int_0^T H(t)\,e^{im\omega t}\,dt` 

182 of the drive, :math:`H(t) = \sum_m H_m e^{-im\omega t}`, for 

183 :math:`|m| \le` *n_harmonics* (exact for a drive with no higher harmonics). 

184 

185 :param ham_t: Callable. See *evolution_operator*. 

186 :param period: Positive real. 

187 :param n_harmonics: Positive integer or zero. 

188 :param n_samples: Positive integer. Default value None (4 n_harmonics + 16). 

189 Time samples per period. 

190 

191 :returns: 

192 * **H_m** -- Dictionary {m: complex ndarray}. 

193 ''' 

194 error_handling.is_callable(ham_t, 'ham_t') 

195 error_handling.positive_real(period, 'period') 

196 error_handling.positive_int_zero(n_harmonics, 'n_harmonics') 

197 if n_samples is None: 

198 n_samples = 4 * n_harmonics + 16 

199 error_handling.positive_int(n_samples, 'n_samples') 

200 omega = 2 * np.pi / period 

201 ts = period * np.arange(n_samples) / n_samples 

202 hams = np.array([np.asarray(ham_t(t), dtype='c16') for t in ts]) 

203 return {m: np.mean(hams * np.exp(1j * m * omega * ts)[:, None, None], axis=0) 

204 for m in range(-n_harmonics, n_harmonics + 1)} 

205 

206 

207def sambe_hamiltonian(harm: dict, omega: float, n_max: int) -> NDArray[np.complex128]: 

208 r''' 

209 Get the Sambe-space Floquet Hamiltonian, truncated to the photon 

210 sectors :math:`n = -n_{max}, \dots, n_{max}`: 

211 

212 .. math:: 

213 

214 (\mathcal{H})_{nn'} = H_{n-n'} + n\omega\,\delta_{nn'}\, . 

215 

216 Its eigenvalues are the quasienergies, repeated in every Floquet zone 

217 (shifted by multiples of :math:`\omega`); those of the central zone 

218 converge fastest as *n_max* grows. 

219 

220 :param harm: Dictionary {m: square complex array}, see *harmonics*. 

221 :param omega: Positive real. Driving frequency. 

222 :param n_max: Positive integer or zero. 

223 

224 :returns: 

225 * **H** -- Complex ndarray, shape ((2 n_max + 1) d, (2 n_max + 1) d). 

226 ''' 

227 error_handling.harmonics(harm) 

228 error_handling.positive_real(omega, 'omega') 

229 error_handling.positive_int_zero(n_max, 'n_max') 

230 d = len(harm[0]) 

231 size = 2 * n_max + 1 

232 out = np.zeros((size * d, size * d), 'c16') 

233 for a, n in enumerate(range(-n_max, n_max + 1)): 

234 for b, n2 in enumerate(range(-n_max, n_max + 1)): 

235 block = np.asarray(harm.get(n - n2, np.zeros((d, d))), dtype='c16') 

236 if n == n2: 

237 block = block + n * omega * np.eye(d) 

238 out[a*d:(a+1)*d, b*d:(b+1)*d] = block 

239 return out 

240 

241 

242# ---------------------------------------------------------------------- 

243# Piecewise-constant (step) drives 

244# ---------------------------------------------------------------------- 

245 

246 

247def _step_evolution(steps: list, t: float | None = None) -> NDArray[np.complex128]: 

248 r''' 

249 Private function. The exact evolution of a piecewise-constant drive, 

250 *steps* = [(H_1, dt_1), (H_2, dt_2), ...] in time order: 

251 U(t) = exp(-i H_n (t - t_{n-1})) ... exp(-i H_1 dt_1), U(T) if t is None. 

252 ''' 

253 u = np.eye(len(steps[0][0]), dtype='c16') 

254 elapsed = 0. 

255 for h, dt in steps: 

256 tau = dt if t is None else min(dt, t - elapsed) 

257 if tau <= 0: 

258 break 

259 u = LA.expm(-1j * tau * h) @ u 

260 elapsed += dt 

261 return u 

262 

263 

264class StepDrive: 

265 r''' 

266 A finite (real-space) system driven by a piecewise-constant 

267 Hamiltonian: :math:`H(t) = H_n` for 

268 :math:`t_{n-1} \le t < t_n = \delta t_1 + \dots + \delta t_n`, repeated 

269 with period :math:`T = \sum_n \delta t_n`. The evolution is exact, a 

270 product of step exponentials, 

271 

272 .. math:: 

273 

274 U(t) = e^{-iH_n(t - t_{n-1})}\,e^{-iH_{n-1}\delta t_{n-1}}\cdots e^{-iH_1\delta t_1}\, , 

275 

276 with no time grid. Built by *step_drive* from **System** instances or 

277 matrices (see *DrivenKSpace* for Bloch models and ribbons). 

278 

279 :param hams: List of square complex ndarrays, one Hermitian matrix per step (all of one shape). 

280 :param durations: List of positive reals, the step durations :math:`\delta t_n`. 

281 :param epsilon: Real number. Default value None. Default branch cut of 

282 *quasienergies* and *effective_hamiltonian*. 

283 ''' 

284 

285 def __init__(self, hams: list, durations: ArrayLike, epsilon: float | None = None) -> None: 

286 error_handling.step_matrices(hams) 

287 error_handling.durations(durations, len(hams)) 

288 error_handling.branch_cut(epsilon) 

289 self.hams = [np.asarray(h, dtype='c16') for h in hams] 

290 self.durations = np.array(durations, dtype='f8') 

291 self.period = float(self.durations.sum()) 

292 self.epsilon = epsilon 

293 self.n_sites = len(self.hams[0]) 

294 

295 def get_ham(self, t: float) -> NDArray[np.complex128]: 

296 r''' 

297 Get the Hamiltonian :math:`H(t)` (periodic in :math:`T`). 

298 

299 :param t: Real number. Time. 

300 

301 :returns: 

302 * **H** -- Complex ndarray. 

303 ''' 

304 error_handling.real_number(t, 't') 

305 n = int(np.searchsorted(np.cumsum(self.durations), t % self.period, side='right')) 

306 return self.hams[min(n, len(self.hams) - 1)] 

307 

308 def evolution_operator(self, t: float | None = None) -> NDArray[np.complex128]: 

309 r''' 

310 Get the exact evolution operator :math:`U(t)` from time 0. 

311 

312 :param t: Real number in :math:`[0, T]`. Default value None (one period, :math:`U(T)`). 

313 

314 :returns: 

315 * **U** -- Complex ndarray. 

316 ''' 

317 error_handling.time_in_period(t, self.period) 

318 return _step_evolution(list(zip(self.hams, self.durations)), t) 

319 

320 def quasienergies(self, epsilon: float | None = None) -> NDArray[np.float64]: 

321 r''' 

322 Get the quasienergies, sorted, in :math:`[-\omega/2, \omega/2)` or 

323 :math:`[\epsilon_{cut} - \omega, \epsilon_{cut})` (see *effective_hamiltonian*). 

324 

325 :param epsilon: Real number. Default value None (the *epsilon* of the drive). 

326 

327 :returns: 

328 * **eps** -- Real ndarray. 

329 ''' 

330 error_handling.branch_cut(epsilon) 

331 epsilon = self.epsilon if epsilon is None else epsilon 

332 return _floquet(self.evolution_operator(), self.period, epsilon)[0] 

333 

334 def effective_hamiltonian(self, epsilon: float | None = None) -> NDArray[np.complex128]: 

335 r''' 

336 Get the Floquet Hamiltonian :math:`H_F = \frac{i}{T}\ln U(T)`, branch 

337 cut at *epsilon* (see *tbkit.floquet.effective_hamiltonian*). 

338 

339 :param epsilon: Real number. Default value None (the *epsilon* of the drive). 

340 

341 :returns: 

342 * **H_F** -- Complex ndarray. 

343 ''' 

344 error_handling.branch_cut(epsilon) 

345 epsilon = self.epsilon if epsilon is None else epsilon 

346 return _heff(self.evolution_operator(), self.period, epsilon) 

347 

348 def evolve(self, psi: NDArray, n_periods: int) -> NDArray[np.complex128]: 

349 r''' 

350 Evolve a state stroboscopically, :math:`\psi(nT) = U(T)^n\psi(0)`. 

351 

352 :param psi: Complex ndarray, the initial state (one amplitude per site). 

353 :param n_periods: Positive integer or zero. 

354 

355 :returns: 

356 * **psi_n** -- Complex ndarray, shape (n_periods + 1, n_sites): 

357 the state after 0, 1, ..., n_periods periods. 

358 ''' 

359 error_handling.ndarray(psi, 'psi', self.n_sites) 

360 error_handling.positive_int_zero(n_periods, 'n_periods') 

361 u = self.evolution_operator() 

362 out = np.zeros((n_periods + 1, self.n_sites), 'c16') 

363 out[0] = psi 

364 for n in range(n_periods): 

365 out[n + 1] = u @ out[n] 

366 return out 

367 

368 

369def step_drive(models: list, durations: ArrayLike, epsilon: float | None = None): 

370 r''' 

371 Build a piecewise-constant drive: the Hamiltonian of ``models[n]`` acts 

372 for ``durations[n]``, in order, and the cycle repeats with period 

373 :math:`T = \sum_n \delta t_n`. The evolution is exact (a product of step 

374 exponentials, no time grid, see *StepDrive*). 

375 

376 :param models: List of models of one size, all of one kind: 

377 

378 * **KSpace** instances (including ribbons made by *kspace.ribbon*) 

379 -- returns a *DrivenKSpace*, whose *get_ham* is the Floquet 

380 Hamiltonian :math:`H_F(\mathbf{k})`; 

381 * **System** instances (their *get_ham* is called) or square 

382 matrices (dense or scipy.sparse) -- returns a *StepDrive*. 

383 

384 Each must be Hermitian; KSpace models need an orthogonal basis. 

385 :param durations: List of positive reals. :math:`\delta t_n`, one per model. 

386 :param epsilon: Real number. Default value None. Branch cut of the 

387 Floquet Hamiltonian (see *effective_hamiltonian*). 

388 

389 :returns: 

390 * **drive** -- *DrivenKSpace* or *StepDrive*. 

391 

392 Example usage:: 

393 

394 # Rudner et al. (2013): hopping along one bond direction per step 

395 drive = step_drive([ks1, ks2, ks3, ks4, ks5], [T/5] * 5) 

396 w = drive.winding_number(epsilon=np.pi / T) 

397 ''' 

398 from tbkit.kspace import KSpace 

399 from tbkit.system import System 

400 kind = error_handling.step_models(models, KSpace, System) 

401 error_handling.durations(durations, len(models)) 

402 error_handling.branch_cut(epsilon) 

403 if kind == 'kspace': 

404 return DrivenKSpace._from_steps(models, durations, epsilon) 

405 hams = [] 

406 for m in models: 

407 if isinstance(m, System): 

408 m.get_ham() 

409 m = m.ham 

410 hams.append(m.toarray() if hasattr(m, 'toarray') else np.asarray(m)) 

411 return StepDrive(hams, durations, epsilon) 

412 

413 

414# ---------------------------------------------------------------------- 

415# Driven Bloch models 

416# ---------------------------------------------------------------------- 

417 

418 

419def _bloch_grid(model, ks: NDArray[np.float64], a: ArrayLike | None = None) -> NDArray[np.complex128]: 

420 r''' 

421 Private function. The Bloch Hamiltonians of a static **KSpace** at many 

422 k-points at once, shape (nk, norb, norb): *get_ham* (or, with a vector 

423 potential *a*, *get_ham_peierls*) vectorized over k. A model whose class 

424 redefines these methods is evaluated point by point. 

425 ''' 

426 from tbkit.kspace import KSpace 

427 if type(model).get_ham is not KSpace.get_ham or type(model).get_ham_peierls is not KSpace.get_ham_peierls: 

428 if a is None: 

429 return np.array([model.get_ham(k) for k in ks]) 

430 return np.array([model.get_ham_peierls(k, a) for k in ks]) 

431 return model._bloch_ham(ks @ model.k_basis.T, None if a is None else np.asarray(a, dtype='f8')) 

432 

433 

434def _spectral_derivative(u: NDArray[np.complex128], axis: int) -> NDArray[np.complex128]: 

435 r''' 

436 Private function. Derivative of a function sampled on a uniform grid 

437 of one period [0, 1) along *axis* (exact for a trigonometric 

438 polynomial, exponentially accurate for an analytic periodic function). 

439 ''' 

440 n = u.shape[axis] 

441 m = np.fft.fftfreq(n, 1. / n) 

442 if n % 2 == 0: 

443 m[n // 2] = 0. # the Nyquist mode has no odd derivative 

444 shape = [1] * u.ndim 

445 shape[axis] = n 

446 return np.fft.ifft(np.fft.fft(u, axis=axis) * (2j * np.pi * m).reshape(shape), axis=axis) 

447 

448 

449def _winding_density(u: NDArray[np.complex128], a_t: NDArray[np.complexfloating]) -> float: 

450 r''' 

451 Private function. Brillouin-zone average of 

452 Tr(A_t [A_1, A_2]), A_i = U^-1 d U / d f_i along the fractional 

453 coordinates f_1, f_2 (u has shape (n1, n2, N, N)). 

454 ''' 

455 u_inv = np.conj(np.swapaxes(u, -1, -2)) 

456 a1 = u_inv @ _spectral_derivative(u, 0) 

457 a2 = u_inv @ _spectral_derivative(u, 1) 

458 return float(np.mean(np.trace(a_t @ (a1 @ a2 - a2 @ a1), axis1=-2, axis2=-1)).real) 

459 

460 

461class DrivenKSpace: 

462 r''' 

463 A periodic model driven by any time-periodic Bloch Hamiltonian 

464 :math:`H(\mathbf{k}, t + T) = H(\mathbf{k}, t)`. 

465 

466 It is a copy of the **KSpace** model whose *get_ham* returns the Floquet 

467 Hamiltonian :math:`H_F(\mathbf{k}) = \frac{i}{T}\ln U(\mathbf{k}, T)` 

468 (branch cut at *epsilon*, see *effective_hamiltonian*), so *get_bands*, 

469 *k_path*, *berry_curvature*, *chern_number*, ... describe the driven 

470 system, and *winding_number* gives its anomalous topology. The 

471 evolution :math:`U(\mathbf{k}, T)` is the product of *n_steps* 

472 midpoint propagators (see *evolution_operator*); drives built by 

473 *step_drive* are evolved exactly. Real-space methods that need static 

474 hoppings (*finite_ham*, *get_ham_beta*, *gbz*) are not available. 

475 

476 :param ks: **KSpace** instance. The model driven: its lattice, orbitals 

477 and k-space conventions (periodic gauge) are those of the driven model. 

478 :param ham_kt: Callable. ``ham_kt(k, t)`` returns the Hermitian Bloch 

479 Hamiltonian :math:`H(\mathbf{k}, t)`, shape (norb, norb), in the 

480 periodic gauge of *ks* (:math:`H(\mathbf{k} + \mathbf{G}, t) = H(\mathbf{k}, t)`). 

481 :param period: Positive real. :math:`T`. 

482 :param n_steps: Positive integer. Default value 100. Time steps per period. 

483 :param epsilon: Real number. Default value None (quasienergies in 

484 :math:`[-\omega/2, \omega/2)`). Branch cut of *get_ham*. 

485 

486 Example usage:: 

487 

488 # a sublattice mass switched on and off 

489 drive = DrivenKSpace(ks, lambda k, t: ks.get_ham(k) + m * np.cos(w * t) * sz, 2 * np.pi / w) 

490 ''' 

491 

492 # Set by _copy (and _from_steps or __new__) on an instance whose class 

493 # also derives from type(ks): the rest is the copied KSpace's. 

494 static: KSpace 

495 period: float 

496 epsilon: float | None 

497 steps: list | None 

498 durations: NDArray[np.float64] 

499 n_steps: int 

500 ham_kt: Callable 

501 lat: Lattice 

502 dim: int 

503 norb: int 

504 spin: bool 

505 rec_vec_k: NDArray[np.float64] 

506 _plane_orientation: Callable 

507 

508 def __new__(cls, ks, ham_kt: Callable, period: float, n_steps: int = 100, epsilon: float | None = None): 

509 from tbkit.kspace import KSpace 

510 error_handling.kspace(ks, KSpace) 

511 error_handling.is_callable(ham_kt, 'ham_kt') 

512 error_handling.positive_real(period, 'period') 

513 error_handling.positive_int(n_steps, 'n_steps') 

514 error_handling.branch_cut(epsilon) 

515 error_handling.drive_matrix(ham_kt(np.zeros(ks.dim), 0.), ks.norb) 

516 obj = cls._copy(ks, period, epsilon) 

517 obj.ham_kt = ham_kt 

518 obj.n_steps = n_steps 

519 return obj 

520 

521 def __init__(self, ks, ham_kt, period, n_steps=100, epsilon=None) -> None: 

522 pass 

523 

524 @classmethod 

525 def _copy(cls, ks, period: float, epsilon: float | None): 

526 ''' 

527 Private method. A copy of *ks* whose class also derives from *cls*. 

528 ''' 

529 driven_cls = type(cls.__name__, (cls, type(ks)), {}) 

530 obj: DrivenKSpace = object.__new__(driven_cls) 

531 obj.__dict__.update(copy.deepcopy(ks.__dict__)) 

532 obj.static = ks 

533 obj.period = period 

534 obj.epsilon = epsilon 

535 obj.steps = None 

536 return obj 

537 

538 @classmethod 

539 def _from_steps(cls, models: list, durations: ArrayLike, epsilon: float | None): 

540 ''' 

541 Private method. The drive of *step_drive*: models[n] for durations[n]. 

542 ''' 

543 durations = np.array(durations, dtype='f8') 

544 obj = cls._copy(models[0], float(durations.sum()), epsilon) 

545 obj.steps = list(models) 

546 obj.durations = durations 

547 obj.n_steps = len(models) 

548 obj.ham_kt = obj._step_ham 

549 return obj 

550 

551 def _step_ham(self, k: ArrayLike, t: float) -> NDArray[np.complex128]: 

552 ''' 

553 Private method. H(k, t) of a step drive. 

554 ''' 

555 steps = self.steps or [] # set by _from_steps, the only caller 

556 n = int(np.searchsorted(np.cumsum(self.durations), t % self.period, side='right')) 

557 return steps[min(n, len(steps) - 1)].get_ham(k) 

558 

559 def _step_list(self, k: NDArray[np.float64]) -> list: 

560 ''' 

561 Private method. The drive at k as steps [(H_n(k), dt_n)]: exact for 

562 a step drive, the midpoint propagators otherwise. 

563 ''' 

564 if self.steps is not None: 

565 return [(m.get_ham(k), d) for m, d in zip(self.steps, self.durations)] 

566 dt = self.period / self.n_steps 

567 return [(np.asarray(self.ham_kt(k, (m + 0.5) * dt), dtype='c16'), dt) for m in range(self.n_steps)] 

568 

569 def _grid_steps(self, ks: NDArray[np.float64]) -> tuple[NDArray[np.complex128], list]: 

570 ''' 

571 Private method. The steps of *_step_list* at many k-points at once: 

572 the Hamiltonians, shape (nk, n_steps, norb, norb), and the durations. 

573 ''' 

574 if self.steps is not None: 

575 return np.stack([_bloch_grid(m, ks) for m in self.steps], axis=1), list(self.durations) 

576 dt = self.period / self.n_steps 

577 hams = np.array([[s[0] for s in self._step_list(k)] for k in ks]) 

578 return hams, [dt] * self.n_steps 

579 

580 def get_evolution(self, k: ArrayLike, t: float | None = None) -> NDArray[np.complex128]: 

581 r''' 

582 Get the evolution operator :math:`U(\mathbf{k}, t) = \mathcal{T}e^{-i\int_0^tH(\mathbf{k},t')dt'}`. 

583 

584 :param k: k point (see *KSpace.get_ham*). 

585 :param t: Real number in :math:`[0, T]`. Default value None (one 

586 period: the Floquet operator :math:`U(\mathbf{k}, T)`). 

587 

588 :returns: 

589 * **U** -- Complex ndarray, shape (norb, norb). 

590 ''' 

591 error_handling.k_vector(k, 'k', self.dim) 

592 error_handling.time_in_period(t, self.period) 

593 k = np.asarray(k, dtype='f8') 

594 if t is None and self.steps is None: 

595 return evolution_operator(lambda tt: self.ham_kt(k, tt), self.period, self.n_steps) 

596 return _step_evolution(self._step_list(k), t) 

597 

598 def get_ham(self, k: ArrayLike) -> NDArray[np.complex128]: 

599 r''' 

600 Get the Floquet Hamiltonian :math:`H_F(\mathbf{k})`, quasienergies in 

601 :math:`[\epsilon_{cut} - \omega, \epsilon_{cut})` (default 

602 :math:`[-\omega/2, \omega/2)`). 

603 

604 :param k: k point (see *KSpace.get_ham*). 

605 

606 :returns: 

607 * **H_F** -- Complex ndarray, shape (norb, norb). 

608 ''' 

609 error_handling.k_vector(k, 'k', self.dim) 

610 return _heff(self.get_evolution(k), self.period, self.epsilon) 

611 

612 def is_hermitian(self) -> bool: 

613 ''' 

614 The Floquet Hamiltonian of a Hermitian drive is Hermitian. 

615 ''' 

616 return True 

617 

618 def _hop_cells(self): 

619 error_handling.not_floquet() 

620 

621 def _check_static(self): 

622 error_handling.not_floquet() 

623 

624 def winding_number(self, epsilon: float, nk: int | tuple[int, int] = 48, n_t: int = 40) -> float: 

625 r''' 

626 Get the winding number of Rudner, Lindner, Berg and Levin (Phys. Rev. 

627 X 3, 031005 (2013)), the number of chiral edge states in the 

628 quasienergy gap at *epsilon*: 

629 

630 .. math:: 

631 

632 W[U_\epsilon] = \frac{1}{8\pi^2}\int_0^T dt\int_{BZ} dk_x dk_y\, 

633 \mathrm{Tr}\left(U_\epsilon^{-1}\partial_tU_\epsilon 

634 \left[U_\epsilon^{-1}\partial_{k_x}U_\epsilon, U_\epsilon^{-1}\partial_{k_y}U_\epsilon\right]\right), 

635 

636 of the loop :math:`U_\epsilon(\mathbf{k}, t) = U(\mathbf{k}, 2t)` for 

637 :math:`t \le T/2`, then :math:`V_\epsilon(\mathbf{k}, 2T - 2t)`, with 

638 the return map :math:`V_\epsilon(\mathbf{k}, t) = e^{-iH_F^\epsilon(\mathbf{k})t}` 

639 and :math:`H_F^\epsilon` the Floquet Hamiltonian with its branch cut 

640 at *epsilon* (see *effective_hamiltonian*). As the integrand is a 

641 3-form, this is the winding of :math:`U` over :math:`[0, T]` minus 

642 that of :math:`V_\epsilon` over :math:`[0, T]`. 

643 

644 Signs. :math:`W` counts the chiral edge states in the gap: 

645 :math:`W = n_{lower} = -n_{upper}` of *edge_state_count*, so 

646 :math:`W > 0` means states moving along :math:`+\mathbf{a}_1` on 

647 the edge of a ribbon on the side of :math:`-\hat{z}\times\mathbf{a}_1` 

648 (the bottom edge of a ribbon along x). Across a group of Floquet bands, 

649 :math:`W[U_{\epsilon'}] - W[U_\epsilon] = -\sum_n C_n` for 

650 :math:`\epsilon < \epsilon'`, the sum over the bands between the two 

651 gaps, with the Chern numbers :math:`C_n` of *chern_number* (the 

652 paper's :math:`\mathcal{C}_n = -\frac{1}{2\pi}\int\nabla\times\mathcal{A}_n` 

653 is minus ours, so there :math:`W' - W = +\sum\mathcal{C}`). A 

654 static model (whose bands fit in the zone below the cut) has 

655 :math:`W = 0` at the zone edge, and :math:`W = -\sum C` of the 

656 bands below any other gap. The anomalous phases are those with 

657 :math:`W \ne 0` in every gap while all the :math:`C_n = 0`. 

658 

659 Discretization: :math:`U(\mathbf{k}, t)` on a uniform 

660 :math:`n_1\times n_2` Brillouin-zone grid, differentiated spectrally 

661 in :math:`\mathbf{k}` (exponentially accurate: both :math:`U` and, 

662 while the gap at *epsilon* is open, :math:`V_\epsilon` are analytic 

663 and periodic in :math:`\mathbf{k}`), and :math:`\partial_tU = -iHU` 

664 exactly, integrated by Gauss-Legendre quadrature on each time step 

665 (*n_t* nodes per period in total, at least 2 per step) and over 

666 :math:`[0, T]` for the return map. Step drives are exact; a 

667 continuous drive is replaced by its *n_steps* midpoint steps, whose 

668 loop has the same :math:`W` as long as the gap stays open. 

669 

670 Tolerance. The result is returned raw. Its error falls off 

671 exponentially with *nk*, at a rate set by the k-width of the 

672 region where :math:`V_\epsilon` varies fast, gap / band velocity 

673 near the smallest gap at *epsilon* (the Berry curvature is 

674 concentrated there). For the five-step model 

675 (:math:`\delta_{AB} = 0.5\pi/T`) the error is below 

676 :math:`10^{-12}` at perfect transfer; with a gap of 

677 :math:`0.067\omega` it is :math:`2\times10^{-5}` at the default 

678 *nk* = 48; with a gap of :math:`0.021\omega` it is 

679 :math:`10^{-2}` at *nk* = 48 and :math:`4\times10^{-4}` at 

680 *nk* = 96. Circularly driven graphene (gap :math:`0.01\omega` at K) 

681 needs *nk* = 96 for :math:`3\times10^{-3}`. Rounding to the nearest 

682 integer is trustworthy when :math:`|W - \mathrm{round}(W)| < 0.01` and the 

683 value is unchanged within that tolerance when *nk* is doubled; a 

684 value far from an integer means the gap is too small for the grid 

685 (or closed). 

686 

687 :param epsilon: Real number. Quasienergy of the gap (branch cut). 

688 :param nk: Positive integer (at least 4), or tuple of 2. Default 

689 value 48. k-points along each reciprocal vector. 

690 :param n_t: Positive integer. Default value 40. Gauss-Legendre 

691 nodes per period (at least 2 per step; 2 n_t, and at least 24, 

692 for the return map). 

693 

694 :returns: 

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

696 ''' 

697 error_handling.dim_exact(self.dim, 2) 

698 error_handling.real_number(epsilon, 'epsilon') 

699 error_handling.nk(nk, 2) 

700 if isinstance(nk, int): 

701 nk = (nk, nk) 

702 error_handling.nk_min(nk, 4) 

703 error_handling.positive_int(n_t, 'n_t') 

704 n1, n2 = nk 

705 fracs = np.stack(np.meshgrid(np.arange(n1) / n1, np.arange(n2) / n2, indexing='ij'), -1) 

706 ks = fracs @ self.rec_vec_k 

707 hams, durations = self._grid_steps(ks.reshape(-1, 2)) 

708 hams = hams.reshape(n1, n2, len(durations), self.norb, self.norb) 

709 x, w = np.polynomial.legendre.leggauss(max(2, -(-n_t // len(durations)))) 

710 # the driven evolution U(k, t), 0 <= t <= T 

711 u = np.broadcast_to(np.eye(self.norb, dtype='c16'), (n1, n2, self.norb, self.norb)) 

712 integral = 0. 

713 for s, dt in enumerate(durations): 

714 h = hams[:, :, s] 

715 e, q = np.linalg.eigh(h) 

716 qh = np.conj(np.swapaxes(q, -1, -2)) 

717 for xj, wj in zip(x, w): 

718 ut = (q * np.exp(-1j * e * dt * (1 + xj) / 2)[..., None, :]) @ qh @ u 

719 a_t = -1j * np.conj(np.swapaxes(ut, -1, -2)) @ h @ ut 

720 integral += wj * dt / 2 * _winding_density(ut, a_t) 

721 u = (q * np.exp(-1j * e * dt)[..., None, :]) @ qh @ u 

722 # the return map V(k, t) = exp(-i H_F t), 0 <= t <= T, run backwards 

723 eps = np.zeros((n1, n2, self.norb)) 

724 z = np.zeros((n1, n2, self.norb, self.norb), 'c16') 

725 for i1 in range(n1): 

726 for i2 in range(n2): 

727 eps[i1, i2], z[i1, i2] = _floquet(u[i1, i2], self.period, epsilon) 

728 zh = np.conj(np.swapaxes(z, -1, -2)) 

729 h_f = (z * eps[..., None, :]) @ zh 

730 x, w = np.polynomial.legendre.leggauss(max(24, 2 * n_t)) 

731 for xj, wj in zip(x, w): 

732 tau = self.period * (1 + xj) / 2 

733 v = (z * np.exp(-1j * eps * tau)[..., None, :]) @ zh 

734 integral -= wj * self.period / 2 * _winding_density(v, -1j * h_f) 

735 return self._plane_orientation((0, 1)) * integral / (8 * np.pi**2) 

736 

737 def edge_state_count( 

738 self, epsilon: float, nk: int = 200, window: float | None = None, 

739 ) -> tuple[float, float]: 

740 r''' 

741 Count the chiral edge states of a driven ribbon (a 1D model, e.g. 

742 made by *step_drive* from *kspace.ribbon* models) crossing the 

743 quasienergy gap at *epsilon*, on each of its two edges. 

744 

745 Following the quasienergies :math:`\epsilon_n(k)` over the 

746 Brillouin zone (:math:`k` along :math:`\mathbf{a}_1`, *nk* points), a 

747 state crossing *epsilon* upwards (:math:`d\epsilon/dk > 0`) counts 

748 +1, downwards -1, weighted by its probability on each half of the 

749 ribbon (the halves below and above its middle line, along 

750 :math:`\hat{z}\times\mathbf{a}_1`). States are followed from one k to 

751 the next by their largest overlap. Only states within *window* of 

752 *epsilon* are followed; bulk states do not cross a gap, so only 

753 edge states contribute. The counts are integers, :math:`\pm1` per 

754 chiral mode, up to the tail of the edge states on the other half 

755 (exponentially small in the width) and the two edges sum to zero. 

756 With the conventions of *winding_number*, 

757 ``n_upper = -n_lower = W`` (checked in the tests and examples). 

758 

759 :param epsilon: Real number. Quasienergy of the gap. 

760 :param nk: Positive integer (at least 4). Default value 200. k-points, 

761 fine enough that a state moves by less than *window* between two. 

762 :param window: Positive real. Default value None (:math:`\omega/8`). 

763 

764 :returns: 

765 * **n_lower** -- Real number. Signed count on the lower edge. 

766 * **n_upper** -- Real number. Signed count on the upper edge. 

767 ''' 

768 error_handling.dim_exact(self.dim, 1) 

769 error_handling.real_number(epsilon, 'epsilon') 

770 error_handling.positive_int(nk, 'nk') 

771 error_handling.nk_min((nk,), 4) 

772 omega = 2 * np.pi / self.period 

773 if window is None: 

774 window = omega / 8 

775 error_handling.positive_real(window, 'window') 

776 a1 = np.asarray(self.lat.prim_vec[0], dtype='f8')[:2] 

777 normal = np.array([-a1[1], a1[0]]) / np.linalg.norm(a1) 

778 across = np.array([dic['r0'][:2] for dic in self.lat.unit_cell], dtype='f8') @ normal 

779 if self.spin: 

780 across = np.repeat(across, 2) 

781 middle = (across.min() + across.max()) / 2 

782 lower = np.where(across < middle, 1., np.where(across > middle, 0., 0.5)) 

783 ks = ((np.arange(nk) + 0.5) / nk)[:, None] * self.rec_vec_k[0][None, :]