Coverage for tbkit/transport.py: 100%
538 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-15 06:13 +0000
1r"""
2Quantum transport through a finite Tight-Binding device attached to
3semi-infinite leads: the Landauer-Buttiker picture, in which the (linear,
4zero-temperature) conductance between two leads is the transmission
5probability at the Fermi energy,
7.. math::
9 G = \frac{e^2}{h}\,T(E_F)\, ,\qquad
10 T(E) = \mathrm{Tr}\left[\Gamma_{out}\,G^r\,\Gamma_{in}\,G^a\right]\, ,
12(the Caroli formula) with :math:`G^r = [E + i\eta - H - \sum_l\Sigma_l]^{-1}`
13the device's retarded Green's function, :math:`\Sigma_l` the self-energy of
14lead :math:`l` and :math:`\Gamma_l = i(\Sigma_l - \Sigma_l^\dagger)`.
16A lead is a semi-infinite repetition of a unit cell with Hamiltonian
17:math:`h_0`, each cell coupled to the next one (away from the device) by
18:math:`v`; its surface cell is coupled to some device sites by
19:math:`\tau`. Its surface Green's function comes from the decimation
20algorithm of Lopez Sancho, Lopez Sancho and Rubio (J. Phys. F 15, 851
21(1985)), and :math:`\Sigma = \tau\, g_s\,\tau^\dagger`.
23Example usage::
25 # a perfect chain of 10 sites between two chain leads: T = 1 in the band
26 tr = Transport(sys.ham)
27 h0, v = np.zeros((1, 1)), np.ones((1, 1))
28 tr.add_lead(h0, v, coupling=np.ones((1, 1)), sites=[0])
29 tr.add_lead(h0, v, coupling=np.ones((1, 1)), sites=[9])
30 T = tr.transmission(np.linspace(-1.9, 1.9, 50))
32*Transport.attach_lead* builds the coupling and finds the device sites
33itself, from a lead given as a 1D *KSpace* model placed in the device's
34coordinates (e.g. ``tr.attach_lead(sys, strip, 1)``).
36*Transport.smatrix* gives the scattering matrix itself. It uses the exact
37propagating and evanescent modes of the leads (*lead_modes*, with no
38broadening) and one sparse factorization of the device per energy (the
39formulation of Groth et al. 2014, as in Kwant). *Transport.wave_function*
40and *Transport.ldos* give the scattering states and the local density of
41states they carry. *Transport.transmission* uses the scattering matrix by
42default.
44With more leads, *Transport.conductance_matrix* and
45*four_terminal_resistance* solve the Landauer-Buttiker equations (Hall
46bars); *bond_currents* and *local_currents* map the current;
47*transmission_eigenvalues*, *shot_noise* and *fano_factor* give the shot
48noise. *RecursiveTransport* computes the two-terminal transmission of long
49quasi-1D devices slice by slice (recursive Green's function), without
50inverting the whole device.
51"""
52from __future__ import annotations
54from typing import NamedTuple
56import numpy as np
57from numpy.typing import ArrayLike, NDArray
58import scipy.linalg as LA
59import scipy.sparse as sp
60import scipy.sparse.linalg as spla
62import tbkit.error_handling as error_handling
65def surface_green(
66 h0: ArrayLike, v: ArrayLike, energy: float, eta: float = 1e-9, tol: float = 1e-12,
67 max_iter: int = 10000,
68) -> NDArray[np.complex128]:
69 r'''
70 Get the surface Green's function of a semi-infinite lead: cells
71 :math:`0, 1, 2, \dots`, with :math:`H_{nn} = h_0` and
72 :math:`H_{n,n+1} = v`, :math:`H_{n+1,n} = v^\dagger`, by the
73 Lopez Sancho-Rubio decimation (each iteration doubles the number of
74 cells accounted for, so it converges in a few tens of steps inside the
75 bands, as :math:`\eta\to0^+`).
77 :param h0: Square complex array. Cell Hamiltonian.
78 :param v: Square complex array, same shape. Coupling of a cell to the next one.
79 :param energy: Real number. Energy.
80 :param eta: Positive real. Default value 1e-9. Broadening.
81 :param tol: Positive real. Default value 1e-12. Convergence threshold on
82 the norm of the renormalized couplings.
83 :param max_iter: Positive integer. Default value 10000.
85 :returns:
86 * **gs** -- Complex ndarray, same shape as *h0*: :math:`g_s = [E+i\eta - h_0 - \Sigma_{bulk}]^{-1}`.
87 '''
88 h0 = np.atleast_2d(np.asarray(h0, dtype='c16'))
89 v = np.atleast_2d(np.asarray(v, dtype='c16'))
90 error_handling.lead(h0, v)
91 error_handling.real_number(energy, 'energy')
92 error_handling.positive_real(eta, 'eta')
93 error_handling.positive_real(tol, 'tol')
94 error_handling.positive_int(max_iter, 'max_iter')
95 gs, _ = _green_surface_bulk(np.array([energy + 1j * eta]), h0, v, v.conj().T, tol, max_iter,
96 'surface_green')
97 return gs[0]
100def _green_surface_bulk(
101 z: NDArray[np.complex128], h0: NDArray[np.complex128], alpha: NDArray[np.complex128],
102 beta: NDArray[np.complex128], tol: float, max_iter: int, name: str,
103) -> tuple[NDArray[np.complex128], NDArray[np.complex128]]:
104 r'''
105 Private function. The surface and bulk Green's functions of one cell
106 of the semi-infinite chain of *_decimation* (*h0*, *alpha*, *beta* of
107 shape (m, m)) at the complex energies *z*, shapes (nz, m, m).
109 The first decimation step inverts :math:`z - h_0`: at an energy equal to
110 an eigenvalue of :math:`h_0` (E = 0 for a plain chain) that is
111 :math:`\sim1/\eta`, and the precision is lost. Regrouping n = 1, 2 or 3
112 cells into one (the same chain) moves those energies: at each energy,
113 the best-conditioned grouping is used.
114 '''
115 m = len(h0)
116 # Scaling cell n by r**n is a similarity transform that leaves the
117 # Green's function of every cell unchanged: balancing the couplings of a
118 # non-reciprocal chain keeps the renormalized ones from overflowing
119 # (r = 1 when beta = alpha^dagger).
120 n_a, n_b = np.linalg.norm(alpha), np.linalg.norm(beta)
121 if n_a > 0 and n_b > 0:
122 r = np.sqrt(n_b / n_a)
123 alpha, beta = alpha * r, beta / r
124 groups, best, choice = [], None, np.zeros(len(z), int)
125 for c, n in enumerate((1, 2, 3)):
126 hn = (np.kron(np.eye(n), h0) + np.kron(np.eye(n, k=1), alpha)
127 + np.kron(np.eye(n, k=-1), beta))
128 an = np.kron(np.eye(n, k=1 - n), alpha)
129 bn = np.kron(np.eye(n, k=n - 1), beta)
130 groups.append((hn, an, bn))
131 smin = np.linalg.svd(z[:, None, None] * np.eye(n * m)[None] - hn, compute_uv=False)[:, -1]
132 if best is None:
133 best = smin
134 else:
135 better = smin > best * (1. + 1e-12)
136 choice[better], best = c, np.where(better, smin, best)
137 gs, gb = np.empty((len(z), m, m), 'c16'), np.empty((len(z), m, m), 'c16')
138 for c, (hn, an, bn) in enumerate(groups):
139 sel = choice == c
140 if np.any(sel):
141 s_, b_ = _decimation(z[sel], hn, an, bn, tol, max_iter, name)
142 gs[sel], gb[sel] = s_[:, :m, :m], b_[:, :m, :m]
143 return gs, gb
146def _decimation(
147 z: NDArray[np.complex128], h0: NDArray[np.complexfloating], alpha: NDArray[np.complexfloating],
148 beta: NDArray[np.complexfloating], tol: float, max_iter: int, name: str,
149) -> tuple[NDArray[np.complex128], NDArray[np.complex128]]:
150 r'''
151 Private function. The Lopez Sancho-Rubio decimation of a semi-infinite
152 chain of cells :math:`0, 1, 2, \dots` (:math:`H_{nn} = h_0`,
153 :math:`H_{n,n+1} = \alpha`, :math:`H_{n+1,n} = \beta`, with
154 :math:`\beta = \alpha^\dagger` unless the chain is non-reciprocal) at
155 many complex energies *z* (shape (nz,)) at once; *h0*, *alpha* and *beta*
156 have shape (m, m) or (nz, m, m). Used by *surface_green* and
157 *KSpace.surface_spectral_function*, which validate the arguments. It
158 stops when the renormalized couplings are below *tol* at every energy
159 (a RuntimeError naming *name* after *max_iter* iterations).
161 :returns:
162 * **gs** -- Complex ndarray, shape (nz, m, m). Surface Green's function (cell 0).
163 * **gb** -- Complex ndarray, shape (nz, m, m). Bulk Green's function
164 (one cell of the infinite chain).
165 '''
166 m = h0.shape[-1]
167 zi = np.asarray(z, dtype='c16')[:, None, None] * np.eye(m)[None]
168 eps_s = np.broadcast_to(h0, zi.shape).astype('c16')
169 eps = eps_s.copy()
170 alpha = np.broadcast_to(alpha, zi.shape).astype('c16')
171 beta = np.broadcast_to(beta, zi.shape).astype('c16')
172 for _ in range(max_iter):
173 g = np.linalg.inv(zi - eps)
174 ag, bg = alpha @ g, beta @ g
175 eps_s = eps_s + ag @ beta
176 eps = eps + ag @ beta + bg @ alpha
177 alpha, beta = ag @ alpha, bg @ beta
178 if np.all(np.linalg.norm(alpha, axis=(1, 2)) + np.linalg.norm(beta, axis=(1, 2)) < tol):
179 break
180 else:
181 error_handling.converged(False, name)
182 return np.linalg.inv(zi - eps_s), np.linalg.inv(zi - eps)
185def lead_from_kspace(ks, direction: int = 1) -> tuple[NDArray[np.complex128], NDArray[np.complex128]]:
186 r'''
187 Get the cell Hamiltonian :math:`h_0` and coupling :math:`v` of a lead
188 from a 1D periodic model (e.g. a ribbon, see *kspace.ribbon*), whose
189 hoppings reach at most the neighbouring cells.
191 :param ks: **KSpace** instance, 1D.
192 :param direction: +1 or -1. Default value 1. Direction, along the
193 model's primitive vector, in which the lead extends away from the
194 device: :math:`v` gathers the hoppings to the cell at
195 :math:`\mathbf{R} = \pm\mathbf{a}_1`.
197 :returns:
198 * **h0** -- Complex ndarray, shape (norb, norb).
199 * **v** -- Complex ndarray, shape (norb, norb).
200 '''
201 error_handling.dim_exact(ks.dim, 1)
202 error_handling.lead_direction(direction)
203 h0 = np.diag(ks.onsite).astype('c16') + ks._onsite_offdiag
204 v = np.zeros((ks.norb, ks.norb), 'c16')
205 for i, j, (n,), t in ks._hop_cells():
206 error_handling.nearest_cells(n)
207 if n == 0:
208 h0[i, j] += t
209 elif n == direction:
210 v[i, j] += t
211 return h0, v
214class LeadModes(NamedTuple):
215 r'''
216 The propagating modes of a lead at one energy (see *lead_modes*):
217 the incoming ones first (moving towards the device, velocity < 0), then
218 the outgoing ones, each group by block of the conservation law, then in
219 increasing momentum.
221 * **wave_functions** -- Complex ndarray, shape (m, n_modes). Mode :math:`\phi`
222 on the lead's surface cell: cell :math:`n` carries :math:`\lambda^n\phi`,
223 with :math:`\lambda = e^{ik}`. Normalized to unit current.
224 * **momenta** -- Real ndarray, shape (n_modes,). :math:`k` in :math:`(-\pi, \pi]`.
225 * **velocities** -- Real ndarray, shape (n_modes,). :math:`dE/dk`,
226 positive away from the device.
227 * **blocks** -- Integer ndarray, shape (n_modes,). The block of the
228 conservation law each mode belongs to: the index of its eigenvalue,
229 in increasing order (all 0 without a conservation law).
230 '''
231 wave_functions: NDArray[np.complex128]
232 momenta: NDArray[np.float64]
233 velocities: NDArray[np.float64]
234 blocks: NDArray[np.int64]
237def lead_modes(
238 h0: ArrayLike, v: ArrayLike, energy: float, conservation_law: ArrayLike | None = None,
239) -> LeadModes:
240 r'''
241 Get the propagating modes of a semi-infinite lead (cells
242 :math:`0, 1, 2, \dots` away from the device, :math:`H_{nn} = h_0`,
243 :math:`H_{n,n+1} = v`) at one energy. A mode :math:`\psi_n = \lambda^n\phi`
244 solves :math:`(E - h_0 - \lambda v - \lambda^{-1}v^\dagger)\phi = 0`, the
245 generalized eigenproblem of the transfer matrix
247 .. math::
249 \begin{pmatrix} 0 & 1\\ -v^\dagger & E - h_0\end{pmatrix}
250 \begin{pmatrix}\phi\\ \lambda\phi\end{pmatrix} = \lambda
251 \begin{pmatrix} 1 & 0\\ 0 & v\end{pmatrix}
252 \begin{pmatrix}\phi\\ \lambda\phi\end{pmatrix}\, ,
254 solved as an ordinary eigenproblem when :math:`v` is well conditioned,
255 and by the QZ algorithm otherwise, so that a singular :math:`v` is
256 allowed (in real arithmetic for a real lead). The
257 modes with :math:`|\lambda| = 1` propagate with the velocity
258 :math:`dE/dk = -2\,\mathrm{Im}(\lambda\phi^\dagger v\phi)` (diagonalized
259 within degenerate modes). There is no broadening :math:`\eta`, so the
260 modes are exact up to the band edges.
262 A *conservation_law* :math:`Q` (spin :math:`\sigma_z`, the electron-hole
263 :math:`\tau_z` of a normal lead, a valley) commutes with :math:`h_0` and
264 :math:`v`. The modes are then found in each eigenspace of :math:`Q`
265 separately, so that every mode carries one of its eigenvalues, even
266 when modes of different blocks are degenerate.
268 :param h0: Square complex array, Hermitian. Cell Hamiltonian.
269 :param v: Square complex array, same shape. Coupling of a cell to the
270 next one, away from the device.
271 :param energy: Real number. Energy, not at a band edge of the lead.
272 :param conservation_law: Hermitian complex array, shape of *h0*, or None.
273 Default value None.
275 :returns:
276 * **modes** -- **LeadModes**: the wave functions, momenta, velocities and blocks.
277 '''
278 h0 = np.atleast_2d(np.asarray(h0, dtype='c16'))
279 v = np.atleast_2d(np.asarray(v, dtype='c16'))
280 error_handling.lead(h0, v)
281 error_handling.real_number(energy, 'energy')
282 law = _law_blocks(h0, v, conservation_law)
283 return _lead_modes(h0, v, float(energy), law)[0]
286def _law_blocks(
287 h0: NDArray[np.complex128], v: NDArray[np.complex128], law: ArrayLike | None,
288) -> list[NDArray[np.complex128]] | None:
289 r'''
290 Private function. Validate a conservation law of a lead and return an
291 orthonormal basis (shape (m, d)) of each of its eigenspaces, in
292 increasing eigenvalue (None without a conservation law).
293 '''
294 if law is None:
295 return None
296 law = np.atleast_2d(np.asarray(law, dtype='c16'))
297 error_handling.conservation_law(law, h0, v)
298 w, u = np.linalg.eigh(law)
299 breaks = np.flatnonzero(np.diff(w) > 1e-8 * max(1., np.max(np.abs(w)))) + 1
300 return [u[:, g] for g in np.split(np.arange(len(w)), breaks)]
303def _lead_modes(
304 h0: NDArray[np.complex128], v: NDArray[np.complex128], energy: float,
305 law: list[NDArray[np.complex128]] | None = None, tol: float = 1e-6,
306) -> tuple[LeadModes, int, NDArray[np.complex128], NDArray[np.complex128]]:
307 r'''
308 Private function. The modes of a lead (see *lead_modes*), the number of
309 incoming ones, and a basis of every solution that leaves the device or
310 decays away from it: its values on cells 0 and 1, *u0* and *u1* of shape
311 (m, m). Its first columns are the outgoing modes of *LeadModes*, the
312 others span the evanescent modes (a Schur basis of the pencil, which
313 stays well conditioned when :math:`v` is singular and many
314 :math:`\lambda` vanish). Modes with :math:`||\lambda| - 1| \le` *tol*
315 propagate. With the eigenspaces *law* of a conservation law
316 (*_law_blocks*), each block :math:`P^\dagger h_0P`, :math:`P^\dagger vP`
317 is solved on its own and embedded back.
318 '''
319 if law is not None:
320 # the incoming modes of every block, then the outgoing ones, then
321 # the evanescent ones: the order of the single-block case
322 phi_parts, k_parts, vel_parts, block, out0, out1, ev0, ev1 = ([] for _ in range(8))
323 for ib, p in enumerate(law):
324 md, n, u0, u1 = _lead_modes(p.conj().T @ h0 @ p, p.conj().T @ v @ p, energy, None, tol)
325 phi_parts.append((p @ md.wave_functions[:, :n], p @ md.wave_functions[:, n:]))
326 k_parts.append((md.momenta[:n], md.momenta[n:]))
327 vel_parts.append((md.velocities[:n], md.velocities[n:]))
328 block.append(np.full(n, ib))
329 out0.append(p @ u0[:, :n])
330 out1.append(p @ u1[:, :n])
331 ev0.append(p @ u0[:, n:])
332 ev1.append(p @ u1[:, n:])
333 modes = LeadModes(np.hstack([f[0] for f in phi_parts] + [f[1] for f in phi_parts]),
334 np.concatenate([f[0] for f in k_parts] + [f[1] for f in k_parts]),
335 np.concatenate([f[0] for f in vel_parts] + [f[1] for f in vel_parts]),
336 np.concatenate(block * 2))
337 return modes, sum(len(b) for b in block), np.hstack(out0 + ev0), np.hstack(out1 + ev1)
338 error_handling.hermitian_operator(h0)
339 m = len(h0)
340 # real arithmetic for a real lead: about twice as fast
341 real = not (np.any(h0.imag) or np.any(v.imag))
342 hr, vr = (h0.real, v.real) if real else (h0, v)
343 output = 'real' if real else 'complex'
344 eye, zero = np.eye(m), np.zeros((m, m))
345 a = np.block([[zero, eye], [-vr.conj().T, energy * eye - hr]])
346 # decaying (|lambda| < 1, including lambda = 0) first; beta = 0 is lambda = infinity
347 if np.linalg.cond(v) < 1e4:
348 # invertible v: the ordinary eigenproblem of B^-1 A, an order of magnitude faster than QZ
349 a[m:] = LA.solve(vr, a[m:])
350 sort = (lambda x, y: np.hypot(x, y) < 1 - tol) if real else (lambda x: np.abs(x) < 1 - tol)
351 _, z, n_dec = LA.schur(a, output=output, sort=sort)
352 alpha, vecs = LA.eig(a)
353 beta = np.ones(2 * m)
354 else:
355 b = np.block([[eye, zero], [zero, vr]])
356 _, _, alpha, beta, _, z = LA.ordqz(a, b, sort=lambda al, be: np.abs(al) < (1 - tol) * np.abs(be),
357 output=output)
358 n_dec = int(np.sum(np.abs(alpha) < (1 - tol) * np.abs(beta)))
359 (alpha, beta), vecs = LA.eig(a, b, homogeneous_eigvals=True)
360 unit = (np.abs(np.abs(alpha) - np.abs(beta)) <= tol * np.abs(beta)) & (beta != 0)
361 lams, vecs = alpha[unit] / beta[unit], vecs[:m, unit]
362 phi_cols: list = []
363 k_list: list = []
364 vel_list: list = []
365 left = list(range(len(lams)))
366 while left:
367 # a group of degenerate lambdas: its modes span the kernel of the
368 # Hermitian E - h0 - lam v - lam^* v^dagger; diagonalize the velocity there
369 group = [i for i in left if abs(lams[i] - lams[left[0]]) < tol]
370 left = [i for i in left if i not in group]
371 lam = np.mean(lams[group])
372 lam /= abs(lam)
373 if len(group) == 1:
374 phi = vecs[:, group] / np.linalg.norm(vecs[:, group])
375 else:
376 # an orthonormal basis of the kernel, which eig does not guarantee
377 w, u = np.linalg.eigh(energy * eye - h0 - lam * v - np.conj(lam) * v.conj().T)
378 phi = u[:, np.argsort(np.abs(w))[:len(group)]]
379 vel, rot = np.linalg.eigh(1j * (lam * phi.conj().T @ v @ phi
380 - np.conj(lam) * phi.conj().T @ v.conj().T @ phi))
381 error_handling.lead_band_edge(np.min(np.abs(vel)) > tol, energy)
382 phi_cols.append(phi @ rot / np.sqrt(np.abs(vel)))
383 k_list += [np.angle(lam)] * len(group)
384 vel_list.append(vel)
385 phis = np.hstack(phi_cols) if phi_cols else np.zeros((m, 0), 'c16')
386 ks, vels = np.array(k_list), np.concatenate(vel_list) if vel_list else np.zeros(0)
387 order = np.lexsort((ks, vels > 0)) # incoming (v < 0) first, then by momentum
388 phis, ks, vels = phis[:, order], ks[order], vels[order]
389 n_in = int(np.sum(vels < 0))
390 error_handling.lead_band_edge(2 * n_in == len(vels) and n_in + n_dec == m, energy)
391 out = phis[:, n_in:]
392 u0 = np.hstack([out, z[:m, :n_dec]])
393 u1 = np.hstack([out * np.exp(1j * ks[n_in:]), z[m:, :n_dec]])
394 return LeadModes(phis, ks, vels, np.zeros(len(ks), int)), n_in, u0, u1
397class SMatrix():
398 r'''
399 The scattering matrix of a device at one energy (see
400 *Transport.smatrix*): ``data[a, b]`` is the amplitude of outgoing mode
401 *a* for an incoming mode *b* of unit current. The rows are the outgoing
402 modes of lead 0, then lead 1, ...; the columns the incoming ones, in the
403 order of *lead_info*. It is unitary for a Hermitian device.
405 A lead is an integer, or a tuple ``(lead, block)`` that keeps only the
406 modes in one block of the lead's conservation law (*Transport.add_lead*),
407 for example ``(0, 1)`` for the electrons and ``(0, 0)`` for the holes of
408 a normal lead with :math:`Q = \tau_z`.
410 :ivar data: Complex ndarray, shape (n_out, n_in).
411 :ivar lead_info: List of **LeadModes**, one per lead.
412 '''
414 def __init__(
415 self, data: NDArray[np.complex128], lead_info: list[LeadModes], n_blocks: list[int] | None = None,
416 ) -> None:
417 self.data = data
418 self.lead_info = lead_info
419 self._n_blocks = [1] * len(lead_info) if n_blocks is None else list(n_blocks)
420 counts = [int(np.sum(m.velocities < 0)) for m in lead_info]
421 self._bounds = np.concatenate([[0], np.cumsum(counts)]).astype(int)
423 def _modes(self, lead: int | tuple[int, int], outgoing: bool) -> NDArray[np.int64]:
424 r'''
425 Private method. The rows (*outgoing*) or columns of the modes of a
426 lead, or of one block of it.
427 '''
428 error_handling.lead_block(lead, self._n_blocks)
429 l = lead[0] if isinstance(lead, tuple) else lead
430 idx = np.arange(self._bounds[l], self._bounds[l + 1])
431 if isinstance(lead, tuple):
432 n = len(idx)
433 blocks = self.lead_info[l].blocks
434 idx = idx[(blocks[n:] if outgoing else blocks[:n]) == lead[1]]
435 return idx
437 def num_propagating(self, lead: int | tuple[int, int]) -> int:
438 r'''
439 Get the number of propagating modes (incoming, as many as outgoing)
440 of a lead, or of one block of it.
442 :param lead: Integer or tuple (lead, block).
444 :returns:
445 * **n** -- Integer.
446 '''
447 return len(self._modes(lead, False))
449 def submatrix(
450 self, lead_out: int | tuple[int, int], lead_in: int | tuple[int, int],
451 ) -> NDArray[np.complex128]:
452 r'''
453 Get the block of the scattering matrix from lead *lead_in* to lead
454 *lead_out*: the transmission amplitudes :math:`t`, or the reflection
455 amplitudes :math:`r` if the two are equal.
457 :param lead_out: Integer or tuple (lead, block).
458 :param lead_in: Integer or tuple (lead, block).
460 :returns:
461 * **s** -- Complex ndarray, shape (num_propagating(lead_out), num_propagating(lead_in)).
462 '''
463 return self.data[np.ix_(self._modes(lead_out, True), self._modes(lead_in, False))]
465 def transmission(self, lead_out: int | tuple[int, int], lead_in: int | tuple[int, int]) -> float:
466 r'''
467 Get the transmission :math:`\sum_{ab}|S_{ab}|^2` from lead *lead_in*
468 to lead *lead_out* (the reflection if they are equal). With blocks,
469 ``transmission((0, 0), (0, 1))`` is, for example, the Andreev
470 reflection of electrons into holes.
472 :param lead_out: Integer or tuple (lead, block).
473 :param lead_in: Integer or tuple (lead, block).
475 :returns:
476 * **T** -- Real number.
477 '''
478 return float(np.sum(np.abs(self.submatrix(lead_out, lead_in)) ** 2))
481class Transport():
482 r'''
483 A finite device, of Hamiltonian *ham* (e.g. *System.ham*), to which
484 semi-infinite leads are attached with *add_lead*; *transmission* then
485 gives the Landauer transmission between two of them.
487 The device is stored sparse. *smatrix*, *wave_function*, *ldos* and
488 *transmission* (by default) only factorize the sparse matrix. The
489 methods with a broadening *eta* use the dense Green's function, and the
490 dense matrix *ham* is built once, the first time it is needed.
492 :param ham: Square matrix (sparse or dense). Device Hamiltonian.
493 '''
495 def __init__(self, ham) -> None:
496 ham = sp.csr_matrix(ham, dtype='c16')
497 error_handling.square_matrix(ham, 'ham')
498 self._ham = ham
499 self._dense: NDArray[np.complex128] | None = None
500 self.leads: list = [] # list of (h0, v, coupling, sites)
501 self._laws: list = [] # eigenspaces of each lead's conservation law, or None
502 self._added_pos = np.zeros((0, 0)) # sites added by attach_lead(extend=True), and tags
503 self._added_tags: list[str] = []
505 @property
506 def ham(self) -> NDArray[np.complex128]:
507 r'''
508 The device Hamiltonian, as a dense complex ndarray.
509 '''
510 if self._dense is None:
511 self._dense = self._ham.toarray()
512 return self._dense
514 def add_lead(
515 self, h0: ArrayLike, v: ArrayLike, coupling: ArrayLike, sites: list[int],
516 conservation_law: ArrayLike | None = None,
517 ) -> None:
518 r'''
519 Attach a semi-infinite lead (see *surface_green*).
521 A *conservation_law* :math:`Q`, Hermitian and commuting with
522 :math:`h_0` and :math:`v`, sorts the lead's modes into the eigenspaces
523 of :math:`Q` (*lead_modes*), and *smatrix* then resolves its blocks:
524 :math:`\sigma_z` for spin-resolved transport, the electron-hole
525 :math:`\tau_z` of a normal lead for Andreev reflection. The device
526 itself need not conserve :math:`Q`.
528 :param h0: Square complex array, shape (m, m). Lead cell Hamiltonian.
529 :param v: Square complex array, shape (m, m). Coupling of a lead cell
530 to the next one, away from the device.
531 :param coupling: Complex array, shape (len(sites), m). Hoppings
532 :math:`\tau` between the device *sites* and the lead's surface cell.
533 :param sites: List of device site (row) indices the lead touches.
534 :param conservation_law: Hermitian complex array, shape (m, m), or None.
535 Default value None. Its blocks are numbered in increasing eigenvalue.
536 '''
537 h0 = np.atleast_2d(np.asarray(h0, dtype='c16'))
538 v = np.atleast_2d(np.asarray(v, dtype='c16'))
539 error_handling.lead(h0, v)
540 coupling = np.atleast_2d(np.asarray(coupling, dtype='c16'))
541 error_handling.lead_coupling(coupling, sites, len(h0), self._ham.shape[0])
542 law = _law_blocks(h0, v, conservation_law)
543 self.leads.append((h0, v, coupling, list(sites)))
544 self._laws.append(law)
546 def attach_lead(
547 self, sys, ks, direction: int = 1, conservation_law: ArrayLike | None = None,
548 extend: bool = False,
549 ) -> tuple[sp.csr_matrix, NDArray[np.float64]] | None:
550 r'''
551 Attach a semi-infinite lead given as a 1D *KSpace* model (e.g. a
552 ribbon, see *kspace.ribbon*) to the device built from *sys*, without
553 writing the coupling by hand.
555 The lead's sites sit at :math:`\boldsymbol\tau_o + n\mathbf{a}_1`
556 (its *unit_cell* positions, in the coordinates of *sys*, and its
557 primitive vector). The device sites found at these positions, with
558 the same tag, are copies of lead orbitals; the outermost cell along
559 *direction* that holds some of them is the interface. The lead
560 starts at the next cell, and couples to the interface sites through
561 its own hoppings: :math:`\tau = v` restricted to their orbitals (see
562 *lead_from_kspace*, which also gives :math:`h_0` and :math:`v`).
563 Every lead orbital with a hopping into the next cell must be present
564 in the interface cell.
566 With *extend*, a device whose edge is ragged relative to the lead
567 (a disc under a strip, a flake under a zigzag lead) is completed
568 first, as Kwant's *attach_lead* does: starting from the lead's first
569 cell, every lead site reached through the lead's hoppings and absent
570 from the device is added to it, cell by cell towards the device,
571 until the sites reached are all present. The added sites take the
572 lead's onsite energies and its hoppings to their neighbours (device
573 or added sites at lead positions). They come after the device
574 sites, so the sites of the leads already attached keep their
575 indices, and a later *attach_lead* with the same *sys* knows of them.
577 :param sys: **System** instance (one orbital per site) whose
578 Hamiltonian, in the same site order, is the device of this
579 *Transport* (followed by the sites added by earlier *extend* calls).
580 :param ks: **KSpace** instance, 1D and spinless, with hoppings
581 between neighbouring cells only.
582 :param direction: +1 or -1. Default value 1. Direction, along the
583 lead's primitive vector, in which the lead extends away from the device.
584 :param conservation_law: Hermitian complex array, shape (norb, norb),
585 or None. Default value None. See *add_lead*.
586 :param extend: Boolean. Default value False. Add the missing lead
587 sites to the device instead of refusing an incomplete interface.
589 :returns:
590 With *extend*, the enlarged device (None otherwise):
592 * **ham** -- Sparse CSR complex matrix. The device Hamiltonian,
593 the added sites last.
594 * **positions** -- Real ndarray, shape (n_sites, space_dim). The
595 positions of its sites, those of *sys* first.
597 Example usage::
599 strip = ribbon(lattices.square(), hoppings, width=5, direction=1)
600 tr = Transport(sys.ham)
601 tr.attach_lead(sys, strip, -1) # at the left end of the device
602 tr.attach_lead(sys, strip, 1) # at the right end
603 ham, pos = tr.attach_lead(disc, strip, 1, extend=True) # ragged edge
604 '''
605 from tbkit.orbital import OrbitalSystem
606 error_handling.sys(sys)
607 error_handling.not_orbital_system(sys, OrbitalSystem)
608 error_handling.empty_coor(sys.lat.coor)
609 error_handling.boolean(extend, 'extend')
610 error_handling.lead_device_size(sys.lat.sites + len(self._added_tags), self._ham.shape[0])
611 error_handling.spinless(ks.spin)
612 h0, v = lead_from_kspace(ks, direction)
613 error_handling.lead_space_dim(ks.space_dim, sys.lat.space_dim)
614 _law_blocks(h0, v, conservation_law)
615 pos = self._site_positions(sys)
616 tags = np.concatenate([sys.lat.coor['tag'], np.array(self._added_tags, dtype='U1')])
617 a = np.array(ks.lat.prim_vec[0], dtype='f8')
618 tau = ks.orbital_positions()
619 cell = np.zeros(len(pos), int)
620 orb = np.full(len(pos), -1)
621 for o in range(ks.norb):
622 rel = pos - tau[o]
623 n = np.rint(rel @ a / (a @ a)).astype(int)
624 hit = (orb < 0) & (np.linalg.norm(rel - n[:, None] * a, axis=1) < error_handling.ATOL) \
625 & (tags == ks.tags[o])
626 cell[hit], orb[hit] = n[hit], o
627 on_lead = orb >= 0
628 error_handling.lead_overlap(bool(on_lead.any()))
629 edge = direction * np.max(direction * cell[on_lead])
630 if extend:
631 added = self._extend(h0, v, direction, cell, orb, edge)
632 new = np.reshape([tau[o] + n * a for n, o in added], (-1, pos.shape[1]))
633 self._added_pos = np.concatenate([self._added_pos.reshape(-1, pos.shape[1]), new])
634 self._added_tags += [str(ks.tags[o]) for _, o in added]
635 cell = np.concatenate([cell, [n for n, _ in added]]).astype(int)
636 orb = np.concatenate([orb, [o for _, o in added]]).astype(int)
637 on_lead = orb >= 0
638 interface = np.flatnonzero(on_lead & (cell == edge))
639 bonded = np.flatnonzero(np.any(v != 0, axis=1)) # orbitals coupled to the next cell
640 error_handling.lead_interface(sorted(set(bonded) - set(orb[interface])))
641 sites = interface[np.isin(orb[interface], bonded)]
642 self.add_lead(h0, v, v[orb[sites]], [int(s) for s in sites], conservation_law)
643 if extend:
644 return self._ham.copy(), self._site_positions(sys)
645 return None
647 def _site_positions(self, sys) -> NDArray[np.float64]:
648 r'''
649 Private method. The positions of the device sites: those of *sys*,
650 then the sites added by *attach_lead* (*extend*).
651 '''
652 pos = np.stack([sys.lat.coor[f] for f in ('x', 'y', 'z')[:sys.lat.space_dim]], axis=1)
653 return np.concatenate([pos, self._added_pos.reshape(-1, pos.shape[1])])
655 def _extend(
656 self, h0: NDArray[np.complex128], v: NDArray[np.complex128], direction: int,
657 cell: NDArray[np.int64], orb: NDArray[np.int64], edge: int,
658 ) -> list[tuple[int, int]]:
659 r'''
660 Private method. Add to the device the lead sites missing between the
661 lead and the device (see *attach_lead*): a flood fill from the lead's
662 first cell (*edge* + *direction*) over the lead's hoppings, which
663 stops at the device sites (*cell*, *orb*: their lead cell and
664 orbital, -1 off the lead) and may not go past the device's innermost
665 lead cell. Returns the added sites as (cell, orbital), in the order
666 they are appended to the device.
667 '''
668 norb = len(h0)
669 present = {(int(n), int(o)): s for s, (n, o) in enumerate(zip(cell, orb)) if o >= 0}
670 inner = np.min(direction * cell[orb >= 0])
671 hop = (h0 != 0) & ~np.eye(norb, dtype=bool)
673 def neighbours(n, o):
674 return [(n, int(j)) for j in np.flatnonzero(hop[o] | hop[:, o])] \
675 + [(n + direction, int(j)) for j in np.flatnonzero(v[o])] \
676 + [(n - direction, int(i)) for i in np.flatnonzero(v[:, o])]
678 stack = [(int(edge) + direction, o) for o in range(norb)]
679 seen = set(stack)
680 added: list[tuple[int, int]] = []
681 while stack:
682 for n, o in neighbours(*stack.pop()):
683 if direction * n > direction * edge or (n, o) in seen:
684 continue
685 seen.add((n, o))
686 if (n, o) not in present:
687 error_handling.lead_extension(direction * n >= inner)
688 added.append((n, o))
689 stack.append((n, o))
690 n_dev = self._ham.shape[0]
691 index = {**present, **{site: n_dev + k for k, site in enumerate(added)}}
692 rows, cols, vals = [], [], []
693 for n, o in added:
694 s = index[(n, o)]
695 for m, p in [(n, o)] + neighbours(n, o):
696 r = index.get((m, p))
697 if r is None:
698 continue
699 t = h0[o, p] if m == n else v[o, p] if m == n + direction else np.conj(v[p, o])
700 rows += [s] if r >= n_dev else [s, r]
701 cols += [r] if r >= n_dev else [r, s]
702 vals += [t] if r >= n_dev else [t, np.conj(t)]
703 size = n_dev + len(added)
704 self._ham = sp.csr_matrix(
705 sp.block_diag([self._ham, sp.csr_matrix((len(added), len(added)))], format='csr')
706 + sp.csr_matrix((vals, (rows, cols)), shape=(size, size)), dtype='c16')
707 self._dense = None
708 return added
710 def self_energy(self, lead: int, energy: float, eta: float = 1e-9) -> NDArray[np.complex128]:
711 r'''
712 Get the self-energy :math:`\Sigma = \tau g_s \tau^\dagger` of a lead,
713 embedded in the device space.
715 :param lead: Integer. Lead index (in the order of *add_lead*).
716 :param energy: Real number.
717 :param eta: Positive real. Default value 1e-9.
719 :returns:
720 * **sigma** -- Complex ndarray, shape of *ham*.
721 '''
722 error_handling.lead_index(lead, len(self.leads))
723 h0, v, tau, sites = self.leads[lead]
724 sigma = np.zeros_like(self.ham)
725 sigma[np.ix_(sites, sites)] = tau @ surface_green(h0, v, energy, eta) @ tau.conj().T
726 return sigma
728 def get_green(self, energy: float, eta: float = 1e-9) -> NDArray[np.complex128]:
729 r'''
730 Get the retarded Green's function of the device with all its leads,
731 :math:`G^r = [E + i\eta - H - \sum_l \Sigma_l]^{-1}`.
733 :param energy: Real number.
734 :param eta: Positive real. Default value 1e-9.
736 :returns:
737 * **G** -- Complex ndarray, shape of *ham*.
738 '''
739 error_handling.real_number(energy, 'energy')
740 error_handling.positive_real(eta, 'eta')
741 sigma = sum((self.self_energy(l, energy, eta) for l in range(len(self.leads))),
742 np.zeros_like(self.ham))
743 return LA.inv((energy + 1j * eta) * np.eye(len(self.ham)) - self.ham - sigma)
745 def transmission(
746 self, energies: ArrayLike, lead_in: int = 0, lead_out: int = 1, eta: float | None = None,
747 ) -> NDArray[np.float64]:
748 r'''
749 Get the transmission from lead *lead_in* to lead *lead_out*: the
750 two-terminal conductance in units of :math:`e^2/h` (per spin, for a
751 spinless model).
753 By default (*eta* None) it is :math:`\sum_{ab}|S_{ab}|^2` over the block of
754 the scattering matrix (*smatrix*), exact up to the band edges of the
755 leads. With a broadening *eta* (and, as a fallback, for non-Hermitian
756 leads or *lead_in* = *lead_out*, with *eta* = 1e-9) it is the Caroli
757 formula :math:`T(E) = \mathrm{Tr}[\Gamma_{out} G^r \Gamma_{in} G^a]`.
758 The two agree to :math:`O(\eta)`.
760 :param energies: Real array. Energies.
761 :param lead_in: Integer. Default value 0.
762 :param lead_out: Integer. Default value 1.
763 :param eta: Positive real or None. Default value None.
765 :returns:
766 * **T** -- Real ndarray, same length as *energies*.
767 '''
768 error_handling.lead_index(lead_in, len(self.leads))
769 error_handling.lead_index(lead_out, len(self.leads))
770 energies = np.atleast_1d(np.asarray(energies, dtype='f8'))
771 hermitian = all(np.allclose(h0, h0.conj().T) for h0, _, _, _ in self.leads)
772 if eta is None and hermitian and lead_in != lead_out:
773 out = np.zeros(len(energies))
774 for n, e in enumerate(energies):
775 modes, sol = self._scattering(float(e), [lead_in], slice(self._ham.shape[0], None))
776 a = sum(len(h0) for h0, _, _, _ in self.leads[:lead_out])
777 n_out = int(np.sum(modes[lead_out].velocities > 0))
778 out[n] = np.sum(np.abs(sol[a:a + n_out]) ** 2)
779 return out
780 eta = 1e-9 if eta is None else eta
781 out = np.zeros(len(energies))
782 for n, e in enumerate(energies):
783 sig_in = self.self_energy(lead_in, e, eta)
784 sig_out = self.self_energy(lead_out, e, eta)
785 gam_in = 1j * (sig_in - sig_in.conj().T)
786 gam_out = 1j * (sig_out - sig_out.conj().T)
787 g = self.get_green(e, eta)
788 out[n] = np.trace(gam_out @ g @ gam_in @ g.conj().T).real
789 return out
791 # ------------------------------------------------------------------
792 # Scattering matrix and scattering states
793 # ------------------------------------------------------------------
795 def _scattering(self, energy: float, inject: list[int], keep: slice = slice(None)) -> tuple:
796 r'''
797 Private method. Solve the scattering problem for every incoming mode
798 of the leads *inject* (Groth et al. 2014). In lead :math:`l`, the
799 wave function on cells 0 and 1 is the incoming mode plus
800 :math:`U_0 b_l`, :math:`U_1 b_l` (*_lead_modes*: outgoing and
801 evanescent). The unknowns are the device wave function
802 :math:`\psi` and the amplitudes :math:`b_l`, and the equations are
803 the Schrodinger equation on the device and on the surface cell of
804 every lead,
806 .. math::
808 (E - H)\psi - \sum_l\tau_l U_{0,l}b_l = \sum_l\tau_l\phi^{in}_l\, ,\qquad
809 -\tau_l^\dagger\psi + [(E - h_0)U_0 - vU_1]b_l = -(E - h_0 - \lambda_{in}v)\phi^{in}_l\, ,
811 one sparse system of size n_sites + sum_l m_l, factorized once
812 (*scipy.sparse.linalg.splu*) for every right-hand side. Only the rows
813 *keep* of the solution are stored (e.g. the lead amplitudes for the
814 scattering matrix).
816 :returns:
817 * **modes** -- List of **LeadModes**, one per lead.
818 * **sol** -- Complex ndarray, shape (n_sites + sum_l m_l, n_in) before
819 *keep*: one column per incoming mode of the leads *inject*, in order.
820 '''
821 error_handling.real_number(energy, 'energy')
822 n = self._ham.shape[0]
823 found = [_lead_modes(h0, v, energy, law) for (h0, v, _, _), law in zip(self.leads, self._laws)]
824 size = n + sum(len(h0) for h0, _, _, _ in self.leads)
825 device = (energy * sp.eye(n, dtype='c16') - self._ham).tocoo()
826 coo = ([device.row], [device.col], [device.data])
827 rhs_coo: tuple = ([], [], [])
829 def block(r, c, values, out=coo):
830 rr, cc = np.meshgrid(r, c, indexing='ij')
831 out[0].append(rr.ravel())
832 out[1].append(cc.ravel())
833 out[2].append(values.ravel())
835 n_rhs, a = 0, n
836 for l, ((h0, v, tau, sites), (modes, n_in, u0, u1)) in enumerate(zip(self.leads, found)):
837 m = len(h0)
838 lead = np.arange(a, a + m)
839 e_h0 = energy * np.eye(m) - h0
840 block(sites, lead, -tau @ u0)
841 block(lead, sites, -tau.conj().T)
842 block(lead, lead, e_h0 @ u0 - v @ u1)
843 if l in inject:
844 phi = modes.wave_functions[:, :n_in]
845 cols = np.arange(n_rhs, n_rhs + n_in)
846 block(sites, cols, tau @ phi, rhs_coo)
847 block(lead, cols, -(e_h0 @ phi - v @ phi * np.exp(1j * modes.momenta[:n_in])), rhs_coo)
848 n_rhs += n_in
849 a += m
850 if n_rhs == 0:
851 return [f[0] for f in found], np.zeros((size, 0), 'c16')[keep]
852 mat, rhs = (sp.csc_matrix((np.concatenate(x), (np.concatenate(r), np.concatenate(c))), shape=(size, s))
853 for (r, c, x), s in ((coo, size), (rhs_coo, n_rhs)))
854 lu = spla.splu(mat)
855 # the right-hand sides are dense in blocks of at most 64 MB
856 chunk = max(1, 2**22 // size)
857 sol = np.hstack([lu.solve(rhs[:, c:c + chunk].toarray())[keep] for c in range(0, n_rhs, chunk)])
858 return [f[0] for f in found], sol
860 def smatrix(self, energy: float) -> SMatrix:
861 r'''
862 Get the scattering matrix at one energy: the amplitudes of the
863 outgoing modes of every lead for each incoming mode of unit
864 current, from the exact lead modes (*lead_modes*, no broadening) and
865 one sparse factorization of the device. It is unitary for a
866 Hermitian device. Every lead must be Hermitian. The blocks of the
867 leads' conservation laws (*add_lead*) are addressed as
868 ``(lead, block)`` in *SMatrix.submatrix* and *SMatrix.transmission*.
870 :param energy: Real number. Energy, not at a band edge of a lead.
872 :returns:
873 * **S** -- **SMatrix**: *data*, *lead_info*, *submatrix*, *transmission*.
875 Example usage::
877 s = tr.smatrix(0.3)
878 t = s.submatrix(1, 0) # transmission amplitudes, lead 0 to lead 1
879 T = s.transmission(1, 0) # = tr.transmission(0.3, 0, 1)
880 '''
881 modes, sol = self._scattering(energy, list(range(len(self.leads))), slice(self._ham.shape[0], None))
882 rows, a = [], 0
883 for (h0, _, _, _), mode in zip(self.leads, modes):
884 rows.append(sol[a:a + int(np.sum(mode.velocities > 0))])
885 a += len(h0)
886 data = np.vstack(rows) if rows else np.zeros((0, 0), 'c16')
887 return SMatrix(data, modes, [1 if law is None else len(law) for law in self._laws])
889 def wave_function(self, energy: float, lead: int) -> NDArray[np.complex128]:
890 r'''
891 Get the scattering states injected from a lead: the wave function on
892 the device of each incoming mode of unit current (*lead_modes*),
893 with the scattered part leaving through every lead.
895 :param energy: Real number. Energy, not at a band edge of a lead.
896 :param lead: Integer. Injecting lead.
898 :returns:
899 * **psi** -- Complex ndarray, shape (n_modes, n_sites): one row
900 per incoming mode of the lead.
901 '''
902 error_handling.lead_index(lead, len(self.leads))
903 _, sol = self._scattering(energy, [lead])
904 return sol[:self._ham.shape[0]].T
906 def ldos(self, energy: float) -> NDArray[np.float64]:
907 r'''
908 Get the local density of states carried by the scattering states,
909 :math:`\rho_i(E) = \frac{1}{2\pi}\sum_{l,a}|\psi_{la}(i)|^2` over every
910 incoming mode :math:`a` of every lead :math:`l` (unit current
911 normalization). The bound states of the device, which no lead
912 feeds, are left out.
914 :param energy: Real number. Energy, not at a band edge of a lead.
916 :returns:
917 * **rho** -- Real ndarray, shape (n_sites,).
918 '''
919 _, sol = self._scattering(energy, list(range(len(self.leads))))
920 return np.sum(np.abs(sol[:self._ham.shape[0]]) ** 2, axis=1) / (2 * np.pi)
922 # ------------------------------------------------------------------
923 # Finite temperature: conductance and thermoelectric coefficients
924 # ------------------------------------------------------------------
926 def _onsager(
927 self, mu: float, temperature: float, lead_in: int, lead_out: int, eta: float | None,
928 n_points: int,
929 ) -> tuple[float, float, float]:
930 r'''
931 Private method. The integrals :math:`L_n = \int dE\,(E - \mu)^n\,T(E)
932 (-\partial f/\partial E)`, n = 0, 1, 2.
933 '''
934 energies, weights = _fermi_window(mu, temperature, n_points)
935 t = self.transmission(energies, lead_in, lead_out, eta) * weights
936 de = energies - mu
937 return float(t.sum()), float(np.sum(t * de)), float(np.sum(t * de**2))
939 def conductance(
940 self, mu: float, temperature: float = 0., lead_in: int = 0, lead_out: int = 1,
941 eta: float | None = None, n_points: int = 101,
942 ) -> float:
943 r'''
944 Get the two-terminal linear conductance at a finite temperature, in
945 units of :math:`e^2/h`,
947 .. math::
949 G = \int dE\,\left(-\frac{\partial f}{\partial E}\right)T(E)\, ,
951 with the Fermi-Dirac :math:`f` at chemical potential :math:`\mu`. The
952 integral runs over :math:`|E - \mu| \le 36k_BT` on *n_points* equally
953 spaced energies (the trapezoidal rule, spectrally accurate for
954 a smooth :math:`T(E)`). Resolving features of :math:`T(E)` narrower than
955 :math:`k_BT` needs more points. At zero temperature it is :math:`T(\mu)`.
957 :param mu: Real number. Chemical potential.
958 :param temperature: Positive real or zero. Default value 0 (:math:`k_B = 1`).
959 :param lead_in: Integer. Default value 0.
960 :param lead_out: Integer. Default value 1.
961 :param eta: Positive real or None. Default value None (see *transmission*).
962 :param n_points: Positive integer. Default value 101.
964 :returns:
965 * **G** -- Real number.
966 '''
967 error_handling.real_number(mu, 'mu')
968 error_handling.positive_real_zero(temperature, 'temperature')
969 error_handling.positive_int(n_points, 'n_points')
970 if temperature == 0:
971 return float(self.transmission([mu], lead_in, lead_out, eta)[0])
972 return self._onsager(mu, temperature, lead_in, lead_out, eta, n_points)[0]
974 def thermoelectric(
975 self, mu: float, temperature: float, lead_in: int = 0, lead_out: int = 1,
976 eta: float | None = None, n_points: int = 101,
977 ) -> tuple[float, float, float]:
978 r'''
979 Get the two-terminal thermoelectric coefficients in linear response
980 (Sivan and Imry 1986), from :math:`L_n = \int dE\,(E - \mu)^n\,T(E)
981 (-\partial f/\partial E)`:
983 .. math::
985 G = L_0\, ,\qquad S = -\frac{L_1}{T L_0}\, ,\qquad
986 \kappa = \frac{1}{T}\left(L_2 - \frac{L_1^2}{L_0}\right)\, ,
988 the conductance (:math:`e^2/h`), the thermopower (:math:`k_B/e`, for
989 electrons of charge :math:`-e`) and the electronic thermal conductance
990 at zero current (:math:`k_B/h` times the energy unit). At low
991 temperature they obey the Mott formula
992 :math:`S = -\frac{\pi^2}{3}T\,\frac{d\ln T(E)}{dE}` and the
993 Wiedemann-Franz law :math:`\kappa = \frac{\pi^2}{3}TG`. The integrals
994 are those of *conductance*.
996 :param mu: Real number. Chemical potential.
997 :param temperature: Positive real (:math:`k_B = 1`).
998 :param lead_in: Integer. Default value 0.
999 :param lead_out: Integer. Default value 1.
1000 :param eta: Positive real or None. Default value None (see *transmission*).
1001 :param n_points: Positive integer. Default value 101.
1003 :returns:
1004 * **G** -- Real number.
1005 * **S** -- Real number (NaN when nothing is transmitted).
1006 * **kappa** -- Real number.
1007 '''
1008 error_handling.real_number(mu, 'mu')
1009 error_handling.positive_real(temperature, 'temperature')
1010 error_handling.positive_int(n_points, 'n_points')
1011 l0, l1, l2 = self._onsager(mu, temperature, lead_in, lead_out, eta, n_points)
1012 if l0 < 1e-14:
1013 return l0, np.nan, l2 / temperature
1014 return l0, -l1 / (temperature * l0), (l2 - l1**2 / l0) / temperature
1016 # ------------------------------------------------------------------
1017 # Multi-terminal conductance, currents and noise
1018 # ------------------------------------------------------------------
1020 def _broadenings(self, energy: float, eta: float) -> list[NDArray[np.complex128]]:
1021 r'''
1022 Private method. The broadening :math:`\Gamma_l = i(\Sigma_l - \Sigma_l^\dagger)`
1023 of every lead, restricted to its sites (shape (len(sites), len(sites))).
1024 '''
1025 out = []
1026 for h0, v, tau, _ in self.leads:
1027 sig = tau @ surface_green(h0, v, energy, eta) @ tau.conj().T
1028 out.append(1j * (sig - sig.conj().T))
1029 return out
1031 def transmission_matrix(self, energy: float, eta: float = 1e-9) -> NDArray[np.float64]:
1032 r'''
1033 Get the transmissions between every pair of leads at one energy,
1034 :math:`T_{pq} = \mathrm{Tr}[\Gamma_pG^r\Gamma_qG^a]` from lead
1035 :math:`q` to lead :math:`p` (the diagonal, a reflection, is set to 0).
1037 :param energy: Real number.
1038 :param eta: Positive real. Default value 1e-9.
1040 :returns:
1041 * **T** -- Real ndarray, shape (n_leads, n_leads).
1042 '''
1043 g = self.get_green(energy, eta)
1044 gam = self._broadenings(energy, eta)
1045 n = len(self.leads)
1046 out = np.zeros((n, n))
1047 for p in range(n):
1048 for q in range(n):
1049 if p != q:
1050 gpq = g[np.ix_(self.leads[p][3], self.leads[q][3])]
1051 out[p, q] = np.trace(gam[p] @ gpq @ gam[q] @ gpq.conj().T).real
1052 return out
1054 def conductance_matrix(
1055 self, energy: float, eta: float = 1e-9, temperature: float = 0., n_points: int = 101,
1056 ) -> NDArray[np.float64]:
1057 r'''
1058 Get the Landauer-Buttiker conductance matrix (Buttiker 1986), in units
1059 of :math:`e^2/h`: the linear-response currents flowing *into* the
1060 device from each lead, :math:`I_p = \sum_q G_{pq}V_q`, with
1062 .. math::
1064 G_{pq} = -T_{pq}\ (p\neq q)\, ,\qquad G_{pp} = \sum_{q\neq p}T_{qp}\, .
1066 Its rows and columns add up to zero (current conservation, and no
1067 current at equal voltages). At a finite *temperature* it is averaged
1068 over the Fermi window, :math:`\int dE\,(-\partial f/\partial E)\,G(E)`
1069 (see *conductance*).
1071 :param energy: Real number. Fermi energy (chemical potential).
1072 :param eta: Positive real. Default value 1e-9.
1073 :param temperature: Positive real or zero. Default value 0 (:math:`k_B = 1`).
1074 :param n_points: Positive integer. Default value 101. Energies in the
1075 Fermi window, at a finite temperature.
1077 :returns:
1078 * **G** -- Real ndarray, shape (n_leads, n_leads).
1079 '''
1080 error_handling.positive_real_zero(temperature, 'temperature')
1081 error_handling.positive_int(n_points, 'n_points')
1082 if temperature > 0:
1083 energies, weights = _fermi_window(energy, temperature, n_points)
1084 return sum(w * self.conductance_matrix(float(e), eta) for e, w in zip(energies, weights))
1085 t = self.transmission_matrix(energy, eta)
1086 return np.diag(t.sum(axis=0)) - t
1088 def four_terminal_resistance(
1089 self, energy: float, current: tuple[int, int], voltage: tuple[int, int], eta: float = 1e-9,
1090 temperature: float = 0., n_points: int = 101,
1091 ) -> float:
1092 r'''
1093 Get a four-terminal (or two-terminal) resistance from the
1094 Landauer-Buttiker equations: a current :math:`I` enters through lead
1095 ``current[0]`` and leaves through ``current[1]``, every other lead is
1096 an ideal voltage probe (no net current), and the result is
1097 :math:`R = (V_a - V_b)/I` with ``voltage = (a, b)``, in units of
1098 :math:`h/e^2`. On a Hall bar, a longitudinal pair gives
1099 :math:`R_{xx}` and a transverse pair :math:`R_{xy}`; in the quantum
1100 Hall regime with :math:`\nu` edge channels, :math:`|R_{xy}| = 1/\nu`
1101 and :math:`R_{xx} = 0` (Buttiker 1988).
1103 :param energy: Real number. Fermi energy.
1104 :param current: Tuple of two distinct lead indices (source, drain).
1105 :param voltage: Tuple of two lead indices (a, b).
1106 :param eta: Positive real. Default value 1e-9.
1107 :param temperature: Positive real or zero. Default value 0 (see *conductance_matrix*).
1108 :param n_points: Positive integer. Default value 101.
1110 :returns:
1111 * **R** -- Real number.
1112 '''
1113 error_handling.lead_pair(current, len(self.leads), 'current', True)
1114 error_handling.lead_pair(voltage, len(self.leads), 'voltage', False)
1115 g = self.conductance_matrix(energy, eta, temperature, n_points)
1116 source, drain = current
1117 keep = [p for p in range(len(self.leads)) if p != drain] # the drain is grounded
1118 currents = np.zeros(len(keep))
1119 currents[keep.index(source)] = 1.
1120 volts = np.zeros(len(self.leads))
1121 volts[keep] = LA.solve(g[np.ix_(keep, keep)], currents)
1122 return float(volts[voltage[0]] - volts[voltage[1]])
1124 def bond_currents(self, energy: float, lead: int, eta: float = 1e-9) -> NDArray[np.float64]:
1125 r'''
1126 Get the bond currents carried, per unit energy, by the electrons
1127 injected from one lead (the non-equilibrium current when that lead's
1128 electrochemical potential is raised),
1130 .. math::
1132 J_{i\to j} = -2\,\mathrm{Im}\left[H_{ij}\,A^{(q)}_{ji}\right]\, ,\qquad
1133 A^{(q)} = G^r\Gamma_qG^a\, ,
1135 in units of :math:`e/h` per unit energy: summed over a cross-section
1136 of a two-terminal device, it is the transmission :math:`T(E)`. The
1137 matrix is antisymmetric, and the currents are conserved at every
1138 site not touched by a lead.
1140 :param energy: Real number.
1141 :param lead: Integer. Injecting lead.
1142 :param eta: Positive real. Default value 1e-9.
1144 :returns:
1145 * **J** -- Real ndarray, shape of *ham*: J[i, j] is the current
1146 from site i to site j.
1147 '''
1148 error_handling.lead_index(lead, len(self.leads))
1149 g = self.get_green(energy, eta)
1150 sites = self.leads[lead][3]
1151 gam = self._broadenings(energy, eta)[lead]
1152 spectral = g[:, sites] @ gam @ g[:, sites].conj().T
1153 return -2 * np.imag(self.ham * spectral.T)
1155 def local_currents(
1156 self, energy: float, lead: int, positions: ArrayLike, eta: float = 1e-9,
1157 ) -> NDArray[np.float64]:
1158 r'''
1159 Get the current density vector on every site from the bond currents
1160 (*bond_currents*): half of the currents of the bonds leaving the
1161 site, along their bond vectors,
1162 :math:`\mathbf{j}_i = \frac12\sum_jJ_{i\to j}(\mathbf{r}_j - \mathbf{r}_i)`.
1164 :param energy: Real number.
1165 :param lead: Integer. Injecting lead.
1166 :param positions: Real array, shape (n_sites, d). Site positions
1167 (e.g. ``np.column_stack([lat.coor['x'], lat.coor['y']])``).
1168 :param eta: Positive real. Default value 1e-9.
1170 :returns:
1171 * **j** -- Real ndarray, shape (n_sites, d).
1172 '''
1173 positions = np.asarray(positions, dtype='f8')
1174 error_handling.site_positions(positions, len(self.ham))
1175 bonds = self.bond_currents(energy, lead, eta)
1176 return 0.5 * (bonds @ positions - bonds.sum(axis=1)[:, None] * positions)
1178 def transmission_eigenvalues(
1179 self, energy: float, lead_in: int = 0, lead_out: int = 1, eta: float = 1e-9,
1180 ) -> NDArray[np.float64]:
1181 r'''
1182 Get the transmission eigenvalues :math:`T_n`, eigenvalues of
1183 :math:`t^\dagger t` with :math:`t = \Gamma_{out}^{1/2}G^r\Gamma_{in}^{1/2}`:
1184 their sum is the transmission.
1186 :param energy: Real number.
1187 :param lead_in: Integer. Default value 0.
1188 :param lead_out: Integer. Default value 1.
1189 :param eta: Positive real. Default value 1e-9.
1191 :returns:
1192 * **T_n** -- Real ndarray, in decreasing order, of length
1193 min(len(sites_in), len(sites_out)).
1194 '''
1195 error_handling.lead_index(lead_in, len(self.leads))
1196 error_handling.lead_index(lead_out, len(self.leads))
1197 g = self.get_green(energy, eta)
1198 gam = self._broadenings(energy, eta)
1199 s_in, s_out = self.leads[lead_in][3], self.leads[lead_out][3]
1200 return _eigen_transmissions(gam[lead_out], g[np.ix_(s_out, s_in)], gam[lead_in])
1202 def shot_noise(
1203 self, energies: ArrayLike, lead_in: int = 0, lead_out: int = 1, eta: float = 1e-9,
1204 ) -> NDArray[np.float64]:
1205 r'''
1206 Get the zero-temperature shot-noise factor :math:`\sum_nT_n(1-T_n)`
1207 (Lesovik 1989; Buttiker 1990): the noise power at bias :math:`V` is
1208 :math:`S = 2\frac{e^3|V|}{h}\sum_nT_n(1-T_n)` (per spin).
1210 :param energies: Real array. Energies.
1211 :param lead_in: Integer. Default value 0.
1212 :param lead_out: Integer. Default value 1.
1213 :param eta: Positive real. Default value 1e-9.
1215 :returns:
1216 * **noise** -- Real ndarray, same length as *energies*.
1217 '''
1218 error_handling.frequencies(energies, 'energies')
1219 energies = np.atleast_1d(np.asarray(energies, dtype='f8'))
1220 return np.array([_noise(self.transmission_eigenvalues(float(e), lead_in, lead_out, eta))
1221 for e in energies])
1223 def fano_factor(
1224 self, energies: ArrayLike, lead_in: int = 0, lead_out: int = 1, eta: float = 1e-9,
1225 ) -> NDArray[np.float64]:
1226 r'''
1227 Get the Fano factor of the zero-temperature shot noise,
1229 .. math::
1231 F = \frac{\sum_nT_n(1-T_n)}{\sum_nT_n}\, ,
1233 the noise in units of the Poisson value :math:`2e|I|`: 0 for a
1234 ballistic conductor (every :math:`T_n` is 0 or 1), close to 1 for a
1235 tunnel barrier, 1/3 for a diffusive wire (Beenakker and Buttiker 1992).
1237 :param energies: Real array. Energies.
1238 :param lead_in: Integer. Default value 0.
1239 :param lead_out: Integer. Default value 1.
1240 :param eta: Positive real. Default value 1e-9.
1242 :returns:
1243 * **F** -- Real ndarray, same length as *energies* (NaN where
1244 nothing is transmitted).
1245 '''
1246 error_handling.frequencies(energies, 'energies')
1247 energies = np.atleast_1d(np.asarray(energies, dtype='f8'))
1248 return np.array([_fano(self.transmission_eigenvalues(float(e), lead_in, lead_out, eta))
1249 for e in energies])
1252def _fermi_window(
1253 mu: float, temperature: float, n_points: int,
1254) -> tuple[NDArray[np.float64], NDArray[np.float64]]:
1255 r'''
1256 Private function. Energies and weights of the trapezoidal rule for
1257 :math:`\int dE\,(-\partial f/\partial E)\,F(E)`, on *n_points* equally
1258 spaced energies with :math:`|E - \mu| \le 36T` (the weight is below
1259 1e-15 beyond). The weights add up to 1.
1260 '''
1261 x = np.linspace(-36., 36., n_points) if n_points > 1 else np.zeros(1)
1262 w = 1 / (4 * np.cosh(x / 2) ** 2)
1263 return mu + temperature * x, w / w.sum()
1266def _sqrt_psd(mat: NDArray[np.complex128]) -> NDArray[np.complex128]:
1267 '''
1268 Private function. Square root of a positive semi-definite Hermitian matrix.
1269 '''
1270 w, u = np.linalg.eigh(0.5 * (mat + mat.conj().T))
1271 return (u * np.sqrt(np.clip(w, 0., None))) @ u.conj().T
1274def _eigen_transmissions(
1275 gam_out: NDArray[np.complex128], g: NDArray[np.complex128], gam_in: NDArray[np.complex128],
1276) -> NDArray[np.float64]:
1277 r'''
1278 Private function. Transmission eigenvalues: the squared singular values
1279 of :math:`\Gamma_{out}^{1/2}G\,\Gamma_{in}^{1/2}` (*g* the block of the
1280 Green's function from the input to the output sites).
1281 '''
1282 t = _sqrt_psd(gam_out) @ g @ _sqrt_psd(gam_in)
1283 return np.linalg.svd(t, compute_uv=False) ** 2
1286def _noise(t: NDArray[np.float64]) -> float:
1287 '''
1288 Private function. Shot-noise factor sum T_n (1 - T_n).
1289 '''
1290 return float(np.sum(t * (1 - t)))
1293def _fano(t: NDArray[np.float64]) -> float:
1294 '''
1295 Private function. Fano factor from transmission eigenvalues (NaN if none transmit).
1296 '''
1297 total = t.sum()
1298 return _noise(t) / float(total) if total > 1e-12 else np.nan
1301def slices_from_positions(x: ArrayLike, tol: float = 1e-6) -> list[list[int]]:
1302 r'''
1303 Group the sites of a device into slices of equal coordinate *x* (e.g.
1304 the columns of a strip along :math:`x`), ordered by increasing *x*, for
1305 *RecursiveTransport*.
1307 :param x: Real array. One coordinate per site (e.g. ``lat.coor['x']``).
1308 :param tol: Positive real. Default value 1e-6. Coordinates closer than
1309 *tol* belong to the same slice.
1311 :returns:
1312 * **slices** -- List of lists of site indices.
1313 '''
1314 x = np.asarray(x, dtype='f8').ravel()
1315 error_handling.ndarray_empty(x, 'x')
1316 error_handling.positive_real(tol, 'tol')
1317 order = np.argsort(x, kind='stable')
1318 breaks = np.flatnonzero(np.diff(x[order]) > tol) + 1
1319 return [sorted(s.tolist()) for s in np.split(order, breaks)]
1322class RecursiveTransport():
1323 r'''
1324 Two-terminal transport through a long quasi-1D device by the recursive
1325 Green's function method (Thouless and Kirkpatrick 1981; Lee and Fisher
1326 1981; MacKinnon 1985). The device is cut into slices :math:`1\dots N`,
1327 each coupled only to its neighbours, with the left lead on slice 1 and
1328 the right lead on slice :math:`N`. The left-connected Green's functions
1330 .. math::
1332 g_1 = [z - H_{11} - \Sigma_L]^{-1}\, ,\qquad
1333 g_{i} = [z - H_{ii} - H_{i,i-1}\,g_{i-1}H_{i-1,i}]^{-1}\, ,
1335 (:math:`\Sigma_R` added in the last one, :math:`z = E + i\eta`) and the
1336 propagator :math:`G_{i1} = g_iH_{i,i-1}G_{i-1,1}` give
1337 :math:`T = \mathrm{Tr}[\Gamma_RG_{N1}\Gamma_LG_{N1}^\dagger]`. The cost
1338 grows linearly with the length, and only slice-sized matrices are
1339 inverted, never the whole device (compare *Transport*, which inverts
1340 the dense device matrix).
1342 :param ham: Square matrix, sparse (preferred for long devices) or dense.
1343 Device Hamiltonian.
1344 :param slices: List of lists of site indices: a partition of the device
1345 into slices, ordered from the left lead to the right one, with
1346 hoppings between neighbouring slices only (see *slices_from_positions*).
1347 :param left: Tuple (h0, v, coupling): the left lead (see *Transport.add_lead*),
1348 with a coupling of shape (len(slices[0]), len(h0)), its rows in the
1349 order of ``slices[0]``.
1350 :param right: Tuple (h0, v, coupling): the right lead, coupled to ``slices[-1]``.
1352 Example usage::
1354 slices = slices_from_positions(lat.coor['x'])
1355 rgf = RecursiveTransport(sys.ham, slices, (h0, v_left, tau), (h0, v_right, tau))
1356 T = rgf.transmission(energies)
1357 '''
1359 def __init__(self, ham, slices: list[list[int]], left: tuple, right: tuple) -> None:
1360 ham = sp.csr_matrix(ham, dtype='c16')
1361 error_handling.square_matrix(ham, 'ham')
1362 error_handling.slices(slices, ham.shape[0])
1363 label = np.empty(ham.shape[0], dtype=int)
1364 pos = np.empty(ham.shape[0], dtype=int)
1365 for n, s in enumerate(slices):
1366 label[s] = n
1367 pos[s] = np.arange(len(s))
1368 coo = ham.tocoo()
1369 coo.sum_duplicates()
1370 lr, lc = label[coo.row], label[coo.col]
1371 error_handling.block_tridiagonal(int(np.abs(lr - lc).max(initial=0)))
1372 self.slices = [list(s) for s in slices]
1373 sizes = [len(s) for s in self.slices]
1374 self.blocks = [np.zeros((m, m), 'c16') for m in sizes]
1375 self.couplings = [np.zeros((m, k), 'c16') for m, k in zip(sizes[:-1], sizes[1:])]
1376 # scatter every nonzero entry into its block (lc == lr) or coupling
1377 # (lc == lr + 1), grouped by slice; the lower couplings are their adjoints
1378 for target, keep in ((self.blocks, lr == lc), (self.couplings, lc == lr + 1)):
1379 order = np.argsort(lr[keep], kind='stable')
1380 rows, cols = pos[coo.row[keep]][order], pos[coo.col[keep]][order]
1381 vals, owner = coo.data[keep][order], lr[keep][order]
1382 bounds = np.searchsorted(owner, np.arange(len(target) + 1))
1383 for n, mat in enumerate(target):
1384 a, b = bounds[n], bounds[n + 1]
1385 mat[rows[a:b], cols[a:b]] = vals[a:b]
1386 self.leads = []
1387 for lead, s in ((left, self.slices[0]), (right, self.slices[-1])):
1388 error_handling.lead_tuple(lead)
1389 h0 = np.atleast_2d(np.asarray(lead[0], dtype='c16'))
1390 v = np.atleast_2d(np.asarray(lead[1], dtype='c16'))
1391 error_handling.lead(h0, v)
1392 tau = np.atleast_2d(np.asarray(lead[2], dtype='c16'))
1393 error_handling.lead_coupling(tau, list(range(len(s))), len(h0), len(s))
1394 self.leads.append((h0, v, tau))
1396 def _propagator(self, energy: float, eta: float) -> tuple:
1397 r'''
1398 Private method. :math:`G_{N1}` and the broadenings of the two leads.
1399 '''
1400 sig = [tau @ surface_green(h0, v, energy, eta) @ tau.conj().T for h0, v, tau in self.leads]
1401 gam = [1j * (s - s.conj().T) for s in sig]
1402 z = energy + 1j * eta
1403 last = len(self.blocks) - 1
1404 g = g_n1 = None
1405 for i, h in enumerate(self.blocks):
1406 a = z * np.eye(len(h)) - h
1407 if i == 0:
1408 a = a - sig[0]
1409 else:
1410 v = self.couplings[i - 1]
1411 a = a - v.conj().T @ g @ v
1412 if i == last:
1413 a = a - sig[1]
1414 g = LA.inv(a)
1415 g_n1 = g if i == 0 else g @ self.couplings[i - 1].conj().T @ g_n1
1416 return g_n1, gam
1418 def transmission(self, energies: ArrayLike, eta: float = 1e-9) -> NDArray[np.float64]:
1419 r'''
1420 Get the transmission from the left lead to the right one: the same
1421 as *Transport.transmission* for the same device and leads.
1423 :param energies: Real array. Energies.
1424 :param eta: Positive real. Default value 1e-9.
1426 :returns:
1427 * **T** -- Real ndarray, same length as *energies*.
1428 '''
1429 error_handling.frequencies(energies, 'energies')
1430 energies = np.atleast_1d(np.asarray(energies, dtype='f8'))
1431 error_handling.positive_real(eta, 'eta')
1432 out = np.zeros(len(energies))
1433 for n, e in enumerate(energies):
1434 g, (gam_l, gam_r) = self._propagator(float(e), eta)
1435 out[n] = np.trace(gam_r @ g @ gam_l @ g.conj().T).real
1436 return out
1438 def transmission_eigenvalues(self, energy: float, eta: float = 1e-9) -> NDArray[np.float64]:
1439 r'''
1440 Get the transmission eigenvalues (see *Transport.transmission_eigenvalues*).
1442 :param energy: Real number.
1443 :param eta: Positive real. Default value 1e-9.
1445 :returns:
1446 * **T_n** -- Real ndarray, in decreasing order.
1447 '''
1448 error_handling.real_number(energy, 'energy')
1449 error_handling.positive_real(eta, 'eta')
1450 g, (gam_l, gam_r) = self._propagator(energy, eta)
1451 return _eigen_transmissions(gam_r, g, gam_l)
1453 def fano_factor(self, energies: ArrayLike, eta: float = 1e-9) -> NDArray[np.float64]:
1454 r'''
1455 Get the Fano factor of the shot noise (see *Transport.fano_factor*).
1457 :param energies: Real array. Energies.
1458 :param eta: Positive real. Default value 1e-9.
1460 :returns:
1461 * **F** -- Real ndarray, same length as *energies*.
1462 '''
1463 error_handling.frequencies(energies, 'energies')
1464 energies = np.atleast_1d(np.asarray(energies, dtype='f8'))
1465 return np.array([_fano(self.transmission_eigenvalues(float(e), eta)) for e in energies])