Coverage for tbkit/transport.py: 100%
578 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-10-30 11:04 +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`.
23Examples
24--------
25Example usage::
27 # a perfect chain of 10 sites between two chain leads: T = 1 in the band
28 tr = Transport(sys.ham)
29 h0, v = np.zeros((1, 1)), np.ones((1, 1))
30 tr.add_lead(h0, v, coupling=np.ones((1, 1)), sites=[0])
31 tr.add_lead(h0, v, coupling=np.ones((1, 1)), sites=[9])
32 T = tr.transmission(np.linspace(-1.9, 1.9, 50))
34*Transport.attach_lead* builds the coupling and finds the device sites
35itself, from a lead given as a 1D *KSpace* model placed in the device's
36coordinates (e.g. ``tr.attach_lead(sys, strip, 1)``).
38*Transport.smatrix* gives the scattering matrix itself. It uses the exact
39propagating and evanescent modes of the leads (*lead_modes*, with no
40broadening) and one sparse factorization of the device per energy (the
41formulation of Groth et al. 2014, as in Kwant). *Transport.wave_function*
42and *Transport.ldos* give the scattering states and the local density of
43states they carry. *Transport.transmission* uses the scattering matrix by
44default.
46With more leads, *Transport.conductance_matrix* and
47*four_terminal_resistance* solve the Landauer-Buttiker equations (Hall
48bars); *bond_currents* and *local_currents* map the current;
49*transmission_eigenvalues*, *shot_noise* and *fano_factor* give the shot
50noise. By default they too come from the scattering matrix and the
51scattering states, on the same sparse factorization. *RecursiveTransport* computes the two-terminal transmission of long
52quasi-1D devices slice by slice (recursive Green's function), without
53inverting the whole device.
54"""
55from __future__ import annotations
57from concurrent.futures import ThreadPoolExecutor
58from typing import Callable, NamedTuple, TypeVar
60import numpy as np
61from numpy.typing import ArrayLike, NDArray
62import scipy.linalg as LA
63import scipy.sparse as sp
64import scipy.sparse.linalg as spla
66import tbkit.error_handling as error_handling
68_T = TypeVar('_T')
70def surface_green(
71 h0: ArrayLike, v: ArrayLike, energy: float, eta: float = 1e-9, tol: float = 1e-12,
72 max_iter: int = 10000,
73) -> NDArray[np.complex128]:
74 r'''
75 Get the surface Green's function of a semi-infinite lead: cells
76 :math:`0, 1, 2, \dots`, with :math:`H_{nn} = h_0` and
77 :math:`H_{n,n+1} = v`, :math:`H_{n+1,n} = v^\dagger`, by the
78 Lopez Sancho-Rubio decimation (each iteration doubles the number of
79 cells accounted for, so it converges in a few tens of steps inside the
80 bands, as :math:`\eta\to0^+`).
82 Parameters
83 ----------
84 h0 : Square complex array
85 Cell Hamiltonian.
86 v : Square complex array, same shape
87 Coupling of a cell to the next one.
88 energy : Real number
89 Energy.
90 eta : Positive real, default 1e-9
91 Broadening.
92 tol : Positive real, default 1e-12
93 Convergence threshold on
94 the norm of the renormalized couplings.
95 max_iter : Positive integer, default 10000
97 Returns
98 -------
99 gs
100 Complex ndarray, same shape as *h0*: :math:`g_s = [E+i\eta - h_0 - \Sigma_{bulk}]^{-1}`.
101 '''
102 h0 = np.atleast_2d(np.asarray(h0, dtype='c16'))
103 v = np.atleast_2d(np.asarray(v, dtype='c16'))
104 error_handling.lead(h0, v)
105 error_handling.real_number(energy, 'energy')
106 error_handling.positive_real(eta, 'eta')
107 error_handling.positive_real(tol, 'tol')
108 error_handling.positive_int(max_iter, 'max_iter')
109 gs, _ = _green_surface_bulk(np.array([energy + 1j * eta]), h0, v, v.conj().T, tol, max_iter,
110 'surface_green')
111 return gs[0]
114def _green_surface_bulk(
115 z: NDArray[np.complex128], h0: NDArray[np.complex128], alpha: NDArray[np.complex128],
116 beta: NDArray[np.complex128], tol: float, max_iter: int, name: str,
117) -> tuple[NDArray[np.complex128], NDArray[np.complex128]]:
118 r'''
119 Private function. The surface and bulk Green's functions of one cell
120 of the semi-infinite chain of *_decimation* (*h0*, *alpha*, *beta* of
121 shape (m, m)) at the complex energies *z*, shapes (nz, m, m).
123 The first decimation step inverts :math:`z - h_0`: at an energy equal to
124 an eigenvalue of :math:`h_0` (E = 0 for a plain chain) that is
125 :math:`\sim1/\eta`, and the precision is lost. Regrouping n = 1, 2 or 3
126 cells into one (the same chain) moves those energies: at each energy,
127 the best-conditioned grouping is used.
128 '''
129 m = len(h0)
130 # Scaling cell n by r**n is a similarity transform that leaves the
131 # Green's function of every cell unchanged: balancing the couplings of a
132 # non-reciprocal chain keeps the renormalized ones from overflowing
133 # (r = 1 when beta = alpha^dagger).
134 n_a, n_b = np.linalg.norm(alpha), np.linalg.norm(beta)
135 if n_a > 0 and n_b > 0:
136 r = np.sqrt(n_b / n_a)
137 alpha, beta = alpha * r, beta / r
138 groups, best, choice = [], None, np.zeros(len(z), int)
139 for c, n in enumerate((1, 2, 3)):
140 hn = (np.kron(np.eye(n), h0) + np.kron(np.eye(n, k=1), alpha)
141 + np.kron(np.eye(n, k=-1), beta))
142 an = np.kron(np.eye(n, k=1 - n), alpha)
143 bn = np.kron(np.eye(n, k=n - 1), beta)
144 groups.append((hn, an, bn))
145 smin = np.linalg.svd(z[:, None, None] * np.eye(n * m)[None] - hn, compute_uv=False)[:, -1]
146 if best is None:
147 best = smin
148 else:
149 better = smin > best * (1. + 1e-12)
150 choice[better], best = c, np.where(better, smin, best)
151 gs, gb = np.empty((len(z), m, m), 'c16'), np.empty((len(z), m, m), 'c16')
152 for c, (hn, an, bn) in enumerate(groups):
153 sel = choice == c
154 if np.any(sel):
155 s_, b_ = _decimation(z[sel], hn, an, bn, tol, max_iter, name)
156 gs[sel], gb[sel] = s_[:, :m, :m], b_[:, :m, :m]
157 return gs, gb
160def _decimation(
161 z: NDArray[np.complex128], h0: NDArray[np.complexfloating], alpha: NDArray[np.complexfloating],
162 beta: NDArray[np.complexfloating], tol: float, max_iter: int, name: str,
163) -> tuple[NDArray[np.complex128], NDArray[np.complex128]]:
164 r'''
165 Private function. The Lopez Sancho-Rubio decimation of a semi-infinite
166 chain of cells :math:`0, 1, 2, \dots` (:math:`H_{nn} = h_0`,
167 :math:`H_{n,n+1} = \alpha`, :math:`H_{n+1,n} = \beta`, with
168 :math:`\beta = \alpha^\dagger` unless the chain is non-reciprocal) at
169 many complex energies *z* (shape (nz,)) at once; *h0*, *alpha* and *beta*
170 have shape (m, m) or (nz, m, m). Used by *surface_green* and
171 *KSpace.surface_spectral_function*, which validate the arguments. It
172 stops when the renormalized couplings are below *tol* at every energy
173 (a RuntimeError naming *name* after *max_iter* iterations).
175 Returns
176 -------
177 gs : Complex ndarray, shape (nz, m, m)
178 Surface Green's function (cell 0).
179 gb : Complex ndarray, shape (nz, m, m)
180 Bulk Green's function
181 (one cell of the infinite chain).
182 '''
183 m = h0.shape[-1]
184 zi = np.asarray(z, dtype='c16')[:, None, None] * np.eye(m)[None]
185 eps_s = np.broadcast_to(h0, zi.shape).astype('c16')
186 eps = eps_s.copy()
187 alpha = np.broadcast_to(alpha, zi.shape).astype('c16')
188 beta = np.broadcast_to(beta, zi.shape).astype('c16')
189 for _ in range(max_iter):
190 g = np.linalg.inv(zi - eps)
191 ag, bg = alpha @ g, beta @ g
192 eps_s = eps_s + ag @ beta
193 eps = eps + ag @ beta + bg @ alpha
194 alpha, beta = ag @ alpha, bg @ beta
195 if np.all(np.linalg.norm(alpha, axis=(1, 2)) + np.linalg.norm(beta, axis=(1, 2)) < tol):
196 break
197 else:
198 error_handling.converged(False, name)
199 return np.linalg.inv(zi - eps_s), np.linalg.inv(zi - eps)
202def lead_from_kspace(ks, direction: int = 1) -> tuple[NDArray[np.complex128], NDArray[np.complex128]]:
203 r'''
204 Get the cell Hamiltonian :math:`h_0` and coupling :math:`v` of a lead
205 from a 1D periodic model (e.g. a ribbon, see *kspace.ribbon*), whose
206 hoppings reach at most the neighbouring cells.
208 Parameters
209 ----------
210 ks : **KSpace** instance, 1D
211 direction : +1 or -1, default 1
212 Direction, along the
213 model's primitive vector, in which the lead extends away from the
214 device: :math:`v` gathers the hoppings to the cell at
215 :math:`\mathbf{R} = \pm\mathbf{a}_1`.
217 Returns
218 -------
219 h0 : Complex ndarray, shape (norb, norb)
220 v : Complex ndarray, shape (norb, norb)
221 '''
222 error_handling.dim_exact(ks.dim, 1)
223 error_handling.lead_direction(direction)
224 h0 = np.diag(ks.onsite).astype('c16') + ks._onsite_offdiag
225 v = np.zeros((ks.norb, ks.norb), 'c16')
226 for i, j, (n,), t in ks._hop_cells():
227 error_handling.nearest_cells(n)
228 if n == 0:
229 h0[i, j] += t
230 elif n == direction:
231 v[i, j] += t
232 return h0, v
235class LeadModes(NamedTuple):
236 r'''
237 The propagating modes of a lead at one energy (see *lead_modes*):
238 the incoming ones first (moving towards the device, velocity < 0), then
239 the outgoing ones, each group by block of the conservation law, then in
240 increasing momentum.
242 * **wave_functions** -- Complex ndarray, shape (m, n_modes). Mode :math:`\phi`
243 on the lead's surface cell: cell :math:`n` carries :math:`\lambda^n\phi`,
244 with :math:`\lambda = e^{ik}`. Normalized to unit current.
245 * **momenta** -- Real ndarray, shape (n_modes,). :math:`k` in :math:`(-\pi, \pi]`.
246 * **velocities** -- Real ndarray, shape (n_modes,). :math:`dE/dk`,
247 positive away from the device.
248 * **blocks** -- Integer ndarray, shape (n_modes,). The block of the
249 conservation law each mode belongs to: the index of its eigenvalue,
250 in increasing order (all 0 without a conservation law).
251 '''
252 wave_functions: NDArray[np.complex128]
253 momenta: NDArray[np.float64]
254 velocities: NDArray[np.float64]
255 blocks: NDArray[np.int64]
258def lead_modes(
259 h0: ArrayLike, v: ArrayLike, energy: float, conservation_law: ArrayLike | None = None,
260) -> LeadModes:
261 r'''
262 Get the propagating modes of a semi-infinite lead (cells
263 :math:`0, 1, 2, \dots` away from the device, :math:`H_{nn} = h_0`,
264 :math:`H_{n,n+1} = v`) at one energy. A mode :math:`\psi_n = \lambda^n\phi`
265 solves :math:`(E - h_0 - \lambda v - \lambda^{-1}v^\dagger)\phi = 0`, the
266 generalized eigenproblem of the transfer matrix
268 .. math::
270 \begin{pmatrix} 0 & 1\\ -v^\dagger & E - h_0\end{pmatrix}
271 \begin{pmatrix}\phi\\ \lambda\phi\end{pmatrix} = \lambda
272 \begin{pmatrix} 1 & 0\\ 0 & v\end{pmatrix}
273 \begin{pmatrix}\phi\\ \lambda\phi\end{pmatrix}\, ,
275 solved as an ordinary eigenproblem when :math:`v` is well conditioned,
276 and by the QZ algorithm otherwise, so that a singular :math:`v` is
277 allowed (in real arithmetic for a real lead). The
278 modes with :math:`|\lambda| = 1` propagate with the velocity
279 :math:`dE/dk = -2\,\mathrm{Im}(\lambda\phi^\dagger v\phi)` (diagonalized
280 within degenerate modes). There is no broadening :math:`\eta`, so the
281 modes are exact up to the band edges.
283 A *conservation_law* :math:`Q` (spin :math:`\sigma_z`, the electron-hole
284 :math:`\tau_z` of a normal lead, a valley) commutes with :math:`h_0` and
285 :math:`v`. The modes are then found in each eigenspace of :math:`Q`
286 separately, so that every mode carries one of its eigenvalues, even
287 when modes of different blocks are degenerate.
289 Parameters
290 ----------
291 h0 : Square complex array, Hermitian
292 Cell Hamiltonian.
293 v : Square complex array, same shape
294 Coupling of a cell to the
295 next one, away from the device.
296 energy : Real number
297 Energy, not at a band edge of the lead.
298 conservation_law : Hermitian complex array, shape of *h0*, or None, default None
300 Returns
301 -------
302 modes : **LeadModes**: the wave functions, momenta, velocities and blocks
303 '''
304 h0 = np.atleast_2d(np.asarray(h0, dtype='c16'))
305 v = np.atleast_2d(np.asarray(v, dtype='c16'))
306 error_handling.lead(h0, v)
307 error_handling.real_number(energy, 'energy')
308 law = _law_blocks(h0, v, conservation_law)
309 return _lead_modes(h0, v, float(energy), law)[0]
312def _law_blocks(
313 h0: NDArray[np.complex128], v: NDArray[np.complex128], law: ArrayLike | None,
314) -> list[NDArray[np.complex128]] | None:
315 r'''
316 Private function. Validate a conservation law of a lead and return an
317 orthonormal basis (shape (m, d)) of each of its eigenspaces, in
318 increasing eigenvalue (None without a conservation law).
319 '''
320 if law is None:
321 return None
322 law = np.atleast_2d(np.asarray(law, dtype='c16'))
323 error_handling.conservation_law(law, h0, v)
324 w, u = np.linalg.eigh(law)
325 breaks = np.flatnonzero(np.diff(w) > 1e-8 * max(1., np.max(np.abs(w)))) + 1
326 return [u[:, g] for g in np.split(np.arange(len(w)), breaks)]
329def _lead_modes(
330 h0: NDArray[np.complex128], v: NDArray[np.complex128], energy: float,
331 law: list[NDArray[np.complex128]] | None = None, tol: float = 1e-6,
332) -> tuple[LeadModes, int, NDArray[np.complex128], NDArray[np.complex128]]:
333 r'''
334 Private function. The modes of a lead (see *lead_modes*), the number of
335 incoming ones, and a basis of every solution that leaves the device or
336 decays away from it: its values on cells 0 and 1, *u0* and *u1* of shape
337 (m, m). Its first columns are the outgoing modes of *LeadModes*, the
338 others span the evanescent modes (a Schur basis of the pencil, which
339 stays well conditioned when :math:`v` is singular and many
340 :math:`\lambda` vanish). Modes with :math:`||\lambda| - 1| \le` *tol*
341 propagate. With the eigenspaces *law* of a conservation law
342 (*_law_blocks*), each block :math:`P^\dagger h_0P`, :math:`P^\dagger vP`
343 is solved on its own and embedded back.
344 '''
345 if law is not None:
346 # the incoming modes of every block, then the outgoing ones, then
347 # the evanescent ones: the order of the single-block case
348 phi_parts, k_parts, vel_parts, block, out0, out1, ev0, ev1 = ([] for _ in range(8))
349 for ib, p in enumerate(law):
350 md, n, u0, u1 = _lead_modes(p.conj().T @ h0 @ p, p.conj().T @ v @ p, energy, None, tol)
351 phi_parts.append((p @ md.wave_functions[:, :n], p @ md.wave_functions[:, n:]))
352 k_parts.append((md.momenta[:n], md.momenta[n:]))
353 vel_parts.append((md.velocities[:n], md.velocities[n:]))
354 block.append(np.full(n, ib))
355 out0.append(p @ u0[:, :n])
356 out1.append(p @ u1[:, :n])
357 ev0.append(p @ u0[:, n:])
358 ev1.append(p @ u1[:, n:])
359 modes = LeadModes(np.hstack([f[0] for f in phi_parts] + [f[1] for f in phi_parts]),
360 np.concatenate([f[0] for f in k_parts] + [f[1] for f in k_parts]),
361 np.concatenate([f[0] for f in vel_parts] + [f[1] for f in vel_parts]),
362 np.concatenate(block * 2))
363 return modes, sum(len(b) for b in block), np.hstack(out0 + ev0), np.hstack(out1 + ev1)
364 error_handling.hermitian_operator(h0)
365 m = len(h0)
366 # real arithmetic for a real lead: about twice as fast
367 real = not (np.any(h0.imag) or np.any(v.imag))
368 hr, vr = (h0.real, v.real) if real else (h0, v)
369 output = 'real' if real else 'complex'
370 eye, zero = np.eye(m), np.zeros((m, m))
371 a = np.block([[zero, eye], [-vr.conj().T, energy * eye - hr]])
372 # decaying (|lambda| < 1, including lambda = 0) first; beta = 0 is lambda = infinity
373 if np.linalg.cond(v) < 1e4:
374 # invertible v: the ordinary eigenproblem of B^-1 A, an order of magnitude faster than QZ
375 a[m:] = LA.solve(vr, a[m:])
376 sort = (lambda x, y: np.hypot(x, y) < 1 - tol) if real else (lambda x: np.abs(x) < 1 - tol)
377 _, z, n_dec = LA.schur(a, output=output, sort=sort)
378 alpha, vecs = LA.eig(a)
379 beta = np.ones(2 * m)
380 else:
381 b = np.block([[eye, zero], [zero, vr]])
382 _, _, alpha, beta, _, z = LA.ordqz(a, b, sort=lambda al, be: np.abs(al) < (1 - tol) * np.abs(be),
383 output=output)
384 n_dec = int(np.sum(np.abs(alpha) < (1 - tol) * np.abs(beta)))
385 (alpha, beta), vecs = LA.eig(a, b, homogeneous_eigvals=True)
386 unit = (np.abs(np.abs(alpha) - np.abs(beta)) <= tol * np.abs(beta)) & (beta != 0)
387 lams, vecs = alpha[unit] / beta[unit], vecs[:m, unit]
388 phi_cols: list = []
389 k_list: list = []
390 vel_list: list = []
391 left = list(range(len(lams)))
392 while left:
393 # a group of degenerate lambdas: its modes span the kernel of the
394 # Hermitian E - h0 - lam v - lam^* v^dagger; diagonalize the velocity there
395 group = [i for i in left if abs(lams[i] - lams[left[0]]) < tol]
396 left = [i for i in left if i not in group]
397 lam = np.mean(lams[group])
398 lam /= abs(lam)
399 if len(group) == 1:
400 phi = vecs[:, group] / np.linalg.norm(vecs[:, group])
401 else:
402 # an orthonormal basis of the kernel, which eig does not guarantee
403 w, u = np.linalg.eigh(energy * eye - h0 - lam * v - np.conj(lam) * v.conj().T)
404 phi = u[:, np.argsort(np.abs(w))[:len(group)]]
405 vel, rot = np.linalg.eigh(1j * (lam * phi.conj().T @ v @ phi
406 - np.conj(lam) * phi.conj().T @ v.conj().T @ phi))
407 error_handling.lead_band_edge(np.min(np.abs(vel)) > tol, energy)
408 phi_cols.append(phi @ rot / np.sqrt(np.abs(vel)))
409 k_list += [np.angle(lam)] * len(group)
410 vel_list.append(vel)
411 phis = np.hstack(phi_cols) if phi_cols else np.zeros((m, 0), 'c16')
412 ks, vels = np.array(k_list), np.concatenate(vel_list) if vel_list else np.zeros(0)
413 order = np.lexsort((ks, vels > 0)) # incoming (v < 0) first, then by momentum
414 phis, ks, vels = phis[:, order], ks[order], vels[order]
415 n_in = int(np.sum(vels < 0))
416 error_handling.lead_band_edge(2 * n_in == len(vels) and n_in + n_dec == m, energy)
417 out = phis[:, n_in:]
418 u0 = np.hstack([out, z[:m, :n_dec]])
419 u1 = np.hstack([out * np.exp(1j * ks[n_in:]), z[m:, :n_dec]])
420 return LeadModes(phis, ks, vels, np.zeros(len(ks), int)), n_in, u0, u1
423class SMatrix():
424 r'''
425 The scattering matrix of a device at one energy (see
426 *Transport.smatrix*): ``data[a, b]`` is the amplitude of outgoing mode
427 *a* for an incoming mode *b* of unit current. The rows are the outgoing
428 modes of lead 0, then lead 1, ...; the columns the incoming ones, in the
429 order of *lead_info*. It is unitary for a Hermitian device.
431 A lead is an integer, or a tuple ``(lead, block)`` that keeps only the
432 modes in one block of the lead's conservation law (*Transport.add_lead*),
433 for example ``(0, 1)`` for the electrons and ``(0, 0)`` for the holes of
434 a normal lead with :math:`Q = \tau_z`.
436 Attributes
437 ----------
438 data : Complex ndarray, shape (n_out, n_in)
439 lead_info : List of **LeadModes**, one per lead
440 '''
442 def __init__(
443 self, data: NDArray[np.complex128], lead_info: list[LeadModes], n_blocks: list[int] | None = None,
444 ) -> None:
445 self.data = data
446 self.lead_info = lead_info
447 self._n_blocks = [1] * len(lead_info) if n_blocks is None else list(n_blocks)
448 counts = [int(np.sum(m.velocities < 0)) for m in lead_info]
449 self._bounds = np.concatenate([[0], np.cumsum(counts)]).astype(int)
451 def _modes(self, lead: int | tuple[int, int], outgoing: bool) -> NDArray[np.int64]:
452 r'''
453 Private method. The rows (*outgoing*) or columns of the modes of a
454 lead, or of one block of it.
455 '''
456 error_handling.lead_block(lead, self._n_blocks)
457 l = lead[0] if isinstance(lead, tuple) else lead
458 idx = np.arange(self._bounds[l], self._bounds[l + 1])
459 if isinstance(lead, tuple):
460 n = len(idx)
461 blocks = self.lead_info[l].blocks
462 idx = idx[(blocks[n:] if outgoing else blocks[:n]) == lead[1]]
463 return idx
465 def num_propagating(self, lead: int | tuple[int, int]) -> int:
466 r'''
467 Get the number of propagating modes (incoming, as many as outgoing)
468 of a lead, or of one block of it.
470 Parameters
471 ----------
472 lead : Integer or tuple (lead, block)
474 Returns
475 -------
476 n : Integer
477 '''
478 return len(self._modes(lead, False))
480 def submatrix(
481 self, lead_out: int | tuple[int, int], lead_in: int | tuple[int, int],
482 ) -> NDArray[np.complex128]:
483 r'''
484 Get the block of the scattering matrix from lead *lead_in* to lead
485 *lead_out*: the transmission amplitudes :math:`t`, or the reflection
486 amplitudes :math:`r` if the two are equal.
488 Parameters
489 ----------
490 lead_out : Integer or tuple (lead, block)
491 lead_in : Integer or tuple (lead, block)
493 Returns
494 -------
495 s
496 Complex ndarray, shape (num_propagating(lead_out), num_propagating(lead_in)).
497 '''
498 return self.data[np.ix_(self._modes(lead_out, True), self._modes(lead_in, False))]
500 def transmission(self, lead_out: int | tuple[int, int], lead_in: int | tuple[int, int]) -> float:
501 r'''
502 Get the transmission :math:`\sum_{ab}|S_{ab}|^2` from lead *lead_in*
503 to lead *lead_out* (the reflection if they are equal). With blocks,
504 ``transmission((0, 0), (0, 1))`` is, for example, the Andreev
505 reflection of electrons into holes.
507 Parameters
508 ----------
509 lead_out : Integer or tuple (lead, block)
510 lead_in : Integer or tuple (lead, block)
512 Returns
513 -------
514 T : Real number
515 '''
516 return float(np.sum(np.abs(self.submatrix(lead_out, lead_in)) ** 2))
519class _EnergyWorkers():
520 r'''
521 Private mixin. The threads over which *Transport* and
522 *RecursiveTransport* spread the energies of their many-energy methods.
523 '''
524 workers = 1 # threads over the energies, see set_workers
526 def set_workers(self, workers: int) -> None:
527 r'''
528 Set the number of threads over which the methods taking many energies
529 spread them: *transmission* and *fano_factor*, and for *Transport*
530 also *shot_noise*, and *conductance*, *thermoelectric*,
531 *conductance_matrix* and *four_terminal_resistance* at a finite
532 temperature (the energies of the Fermi window). Each
533 energy is one sparse factorization or dense inversion, mostly in
534 LAPACK and SuperLU, which release the GIL, so the energies run in
535 parallel. The results do not depend on *workers*.
537 Parameters
538 ----------
539 workers : Positive integer
540 Number of threads (default 1, no threads).
542 Examples
543 --------
544 Example usage::
546 tr.set_workers(4)
547 T = tr.transmission(np.linspace(-2, 2, 400))
548 '''
549 error_handling.positive_int(workers, 'workers')
550 self.workers = workers
552 def _map(self, func: Callable[[float], _T], energies: NDArray[np.float64]) -> list[_T]:
553 r'''
554 Private method. *func* at every energy, in order, spread over
555 *workers* threads.
556 '''
557 if self.workers > 1 and len(energies) > 1:
558 with ThreadPoolExecutor(min(self.workers, len(energies))) as pool:
559 return list(pool.map(func, (float(e) for e in energies)))
560 return [func(float(e)) for e in energies]
563class Transport(_EnergyWorkers):
564 r'''
565 A finite device, of Hamiltonian *ham* (e.g. *System.ham*), to which
566 semi-infinite leads are attached with *add_lead*; *transmission* then
567 gives the Landauer transmission between two of them.
569 The device is stored sparse. *smatrix*, *wave_function*, *ldos*, and by
570 default (*eta* None) the transmissions, conductance matrix, currents and
571 noise, only factorize the sparse matrix. With a broadening *eta* (and,
572 as a fallback, for non-Hermitian leads, with *eta* = 1e-9) these use the
573 dense Green's function, and the dense matrix *ham* is built once, the
574 first time it is needed.
576 Parameters
577 ----------
578 ham : Square matrix (sparse or dense)
579 Device Hamiltonian.
580 '''
582 def __init__(self, ham) -> None:
583 ham = sp.csr_matrix(ham, dtype='c16')
584 error_handling.square_matrix(ham, 'ham')
585 self._ham = ham
586 self._dense: NDArray[np.complex128] | None = None
587 self.leads: list = [] # list of (h0, v, coupling, sites)
588 self._laws: list = [] # eigenspaces of each lead's conservation law, or None
589 self._added_pos = np.zeros((0, 0)) # sites added by attach_lead(extend=True), and tags
590 self._added_tags: list[str] = []
592 @property
593 def ham(self) -> NDArray[np.complex128]:
594 r'''
595 The device Hamiltonian, as a dense complex ndarray.
596 '''
597 if self._dense is None:
598 self._dense = self._ham.toarray()
599 return self._dense
601 def add_lead(
602 self, h0: ArrayLike, v: ArrayLike, coupling: ArrayLike, sites: list[int],
603 conservation_law: ArrayLike | None = None,
604 ) -> None:
605 r'''
606 Attach a semi-infinite lead (see *surface_green*).
608 A *conservation_law* :math:`Q`, Hermitian and commuting with
609 :math:`h_0` and :math:`v`, sorts the lead's modes into the eigenspaces
610 of :math:`Q` (*lead_modes*), and *smatrix* then resolves its blocks:
611 :math:`\sigma_z` for spin-resolved transport, the electron-hole
612 :math:`\tau_z` of a normal lead for Andreev reflection. The device
613 itself need not conserve :math:`Q`.
615 Parameters
616 ----------
617 h0 : Square complex array, shape (m, m)
618 Lead cell Hamiltonian.
619 v : Square complex array, shape (m, m)
620 Coupling of a lead cell
621 to the next one, away from the device.
622 coupling : Complex array, shape (len(sites), m)
623 Hoppings
624 :math:`\tau` between the device *sites* and the lead's surface cell.
625 sites : List of device site (row) indices the lead touches
626 conservation_law : Hermitian complex array, shape (m, m), or None, default None
627 Its blocks are numbered in increasing eigenvalue.
628 '''
629 h0 = np.atleast_2d(np.asarray(h0, dtype='c16'))
630 v = np.atleast_2d(np.asarray(v, dtype='c16'))
631 error_handling.lead(h0, v)
632 coupling = np.atleast_2d(np.asarray(coupling, dtype='c16'))
633 error_handling.lead_coupling(coupling, sites, len(h0), self._ham.shape[0])
634 law = _law_blocks(h0, v, conservation_law)
635 self.leads.append((h0, v, coupling, list(sites)))
636 self._laws.append(law)
638 def attach_lead(
639 self, sys, ks, direction: int = 1, conservation_law: ArrayLike | None = None,
640 extend: bool = False,
641 ) -> tuple[sp.csr_matrix, NDArray[np.float64]] | None:
642 r'''
643 Attach a semi-infinite lead given as a 1D *KSpace* model (e.g. a
644 ribbon, see *kspace.ribbon*) to the device built from *sys*, without
645 writing the coupling by hand. Returns, with *extend*, the enlarged
646 device (None otherwise).
648 The lead's sites sit at :math:`\boldsymbol\tau_o + n\mathbf{a}_1`
649 (its *unit_cell* positions, in the coordinates of *sys*, and its
650 primitive vector). The device sites found at these positions, with
651 the same tag, are copies of lead orbitals; the outermost cell along
652 *direction* that holds some of them is the interface. The lead
653 starts at the next cell, and couples to the interface sites through
654 its own hoppings: :math:`\tau = v` restricted to their orbitals (see
655 *lead_from_kspace*, which also gives :math:`h_0` and :math:`v`).
656 Every lead orbital with a hopping into the next cell must be present
657 in the interface cell.
659 With *extend*, a device whose edge is ragged relative to the lead
660 (a disc under a strip, a flake under a zigzag lead) is completed
661 first, as Kwant's *attach_lead* does: starting from the lead's first
662 cell, every lead site reached through the lead's hoppings and absent
663 from the device is added to it, cell by cell towards the device,
664 until the sites reached are all present. The added sites take the
665 lead's onsite energies and its hoppings to their neighbours (device
666 or added sites at lead positions). They come after the device
667 sites, so the sites of the leads already attached keep their
668 indices, and a later *attach_lead* with the same *sys* knows of them.
670 Parameters
671 ----------
672 sys
673 **System** instance (one orbital per site) whose
674 Hamiltonian, in the same site order, is the device of this
675 *Transport* (followed by the sites added by earlier *extend* calls).
676 ks
677 **KSpace** instance, 1D and spinless, with hoppings
678 between neighbouring cells only.
679 direction : +1 or -1, default 1
680 Direction, along the
681 lead's primitive vector, in which the lead extends away from the device.
682 conservation_law : Hermitian complex array, shape (norb, norb), or None, default None
683 See *add_lead*.
684 extend : Boolean, default False
685 Add the missing lead
686 sites to the device instead of refusing an incomplete interface.
688 Returns
689 -------
690 ham : Sparse CSR complex matrix
691 The device Hamiltonian,
692 the added sites last.
693 positions : Real ndarray, shape (n_sites, space_dim)
694 The
695 positions of its sites, those of *sys* first.
697 Examples
698 --------
699 Example usage::
701 strip = ribbon(lattices.square(), hoppings, width=5, direction=1)
702 tr = Transport(sys.ham)
703 tr.attach_lead(sys, strip, -1) # at the left end of the device
704 tr.attach_lead(sys, strip, 1) # at the right end
705 ham, pos = tr.attach_lead(disc, strip, 1, extend=True) # ragged edge
706 '''
707 from tbkit.orbital import OrbitalSystem
708 error_handling.sys(sys)
709 error_handling.not_orbital_system(sys, OrbitalSystem)
710 error_handling.empty_coor(sys.lat.coor)
711 error_handling.boolean(extend, 'extend')
712 error_handling.lead_device_size(sys.lat.sites + len(self._added_tags), self._ham.shape[0])
713 error_handling.spinless(ks.spin)
714 h0, v = lead_from_kspace(ks, direction)
715 error_handling.lead_space_dim(ks.space_dim, sys.lat.space_dim)
716 _law_blocks(h0, v, conservation_law)
717 pos = self._site_positions(sys)
718 tags = np.concatenate([sys.lat.coor['tag'], np.array(self._added_tags, dtype='U1')])
719 a = np.array(ks.lat.prim_vec[0], dtype='f8')
720 tau = ks.orbital_positions()
721 cell = np.zeros(len(pos), int)
722 orb = np.full(len(pos), -1)
723 for o in range(ks.norb):
724 rel = pos - tau[o]
725 n = np.rint(rel @ a / (a @ a)).astype(int)
726 hit = (orb < 0) & (np.linalg.norm(rel - n[:, None] * a, axis=1) < error_handling.ATOL) \
727 & (tags == ks.tags[o])
728 cell[hit], orb[hit] = n[hit], o
729 on_lead = orb >= 0
730 error_handling.lead_overlap(bool(on_lead.any()))
731 edge = direction * np.max(direction * cell[on_lead])
732 if extend:
733 added = self._extend(h0, v, direction, cell, orb, edge)
734 new = np.reshape([tau[o] + n * a for n, o in added], (-1, pos.shape[1]))
735 self._added_pos = np.concatenate([self._added_pos.reshape(-1, pos.shape[1]), new])
736 self._added_tags += [str(ks.tags[o]) for _, o in added]
737 cell = np.concatenate([cell, [n for n, _ in added]]).astype(int)
738 orb = np.concatenate([orb, [o for _, o in added]]).astype(int)
739 on_lead = orb >= 0
740 interface = np.flatnonzero(on_lead & (cell == edge))
741 bonded = np.flatnonzero(np.any(v != 0, axis=1)) # orbitals coupled to the next cell
742 error_handling.lead_interface(sorted(set(bonded) - set(orb[interface])))
743 sites = interface[np.isin(orb[interface], bonded)]
744 self.add_lead(h0, v, v[orb[sites]], [int(s) for s in sites], conservation_law)
745 if extend:
746 return self._ham.copy(), self._site_positions(sys)
747 return None
749 def _site_positions(self, sys) -> NDArray[np.float64]:
750 r'''
751 Private method. The positions of the device sites: those of *sys*,
752 then the sites added by *attach_lead* (*extend*).
753 '''
754 pos = np.stack([sys.lat.coor[f] for f in ('x', 'y', 'z')[:sys.lat.space_dim]], axis=1)
755 return np.concatenate([pos, self._added_pos.reshape(-1, pos.shape[1])])
757 def _extend(
758 self, h0: NDArray[np.complex128], v: NDArray[np.complex128], direction: int,
759 cell: NDArray[np.int64], orb: NDArray[np.int64], edge: int,
760 ) -> list[tuple[int, int]]:
761 r'''
762 Private method. Add to the device the lead sites missing between the
763 lead and the device (see *attach_lead*): a flood fill from the lead's
764 first cell (*edge* + *direction*) over the lead's hoppings, which
765 stops at the device sites (*cell*, *orb*: their lead cell and
766 orbital, -1 off the lead) and may not go past the device's innermost
767 lead cell. Returns the added sites as (cell, orbital), in the order
768 they are appended to the device.
769 '''
770 norb = len(h0)
771 present = {(int(n), int(o)): s for s, (n, o) in enumerate(zip(cell, orb)) if o >= 0}
772 inner = np.min(direction * cell[orb >= 0])
773 hop = (h0 != 0) & ~np.eye(norb, dtype=bool)
775 def neighbours(n, o):
776 return [(n, int(j)) for j in np.flatnonzero(hop[o] | hop[:, o])] \
777 + [(n + direction, int(j)) for j in np.flatnonzero(v[o])] \
778 + [(n - direction, int(i)) for i in np.flatnonzero(v[:, o])]
780 stack = [(int(edge) + direction, o) for o in range(norb)]
781 seen = set(stack)
782 added: list[tuple[int, int]] = []
783 while stack:
784 for n, o in neighbours(*stack.pop()):
785 if direction * n > direction * edge or (n, o) in seen:
786 continue
787 seen.add((n, o))
788 if (n, o) not in present:
789 error_handling.lead_extension(direction * n >= inner)
790 added.append((n, o))
791 stack.append((n, o))
792 n_dev = self._ham.shape[0]
793 index = {**present, **{site: n_dev + k for k, site in enumerate(added)}}
794 rows, cols, vals = [], [], []
795 for n, o in added:
796 s = index[(n, o)]
797 for m, p in [(n, o)] + neighbours(n, o):
798 r = index.get((m, p))
799 if r is None:
800 continue
801 t = h0[o, p] if m == n else v[o, p] if m == n + direction else np.conj(v[p, o])
802 rows += [s] if r >= n_dev else [s, r]
803 cols += [r] if r >= n_dev else [r, s]
804 vals += [t] if r >= n_dev else [t, np.conj(t)]
805 size = n_dev + len(added)
806 self._ham = sp.csr_matrix(
807 sp.block_diag([self._ham, sp.csr_matrix((len(added), len(added)))], format='csr')
808 + sp.csr_matrix((vals, (rows, cols)), shape=(size, size)), dtype='c16')
809 self._dense = None
810 return added
812 def self_energy(self, lead: int, energy: float, eta: float = 1e-9) -> NDArray[np.complex128]:
813 r'''
814 Get the self-energy :math:`\Sigma = \tau g_s \tau^\dagger` of a lead,
815 embedded in the device space.
817 Parameters
818 ----------
819 lead : Integer
820 Lead index (in the order of *add_lead*).
821 energy : Real number
822 eta : Positive real, default 1e-9
824 Returns
825 -------
826 sigma : Complex ndarray, shape of *ham*
827 '''
828 error_handling.lead_index(lead, len(self.leads))
829 h0, v, tau, sites = self.leads[lead]
830 sigma = np.zeros_like(self.ham)
831 sigma[np.ix_(sites, sites)] = tau @ surface_green(h0, v, energy, eta) @ tau.conj().T
832 return sigma
834 def get_green(self, energy: float, eta: float = 1e-9) -> NDArray[np.complex128]:
835 r'''
836 Get the retarded Green's function of the device with all its leads,
837 :math:`G^r = [E + i\eta - H - \sum_l \Sigma_l]^{-1}`.
839 Parameters
840 ----------
841 energy : Real number
842 eta : Positive real, default 1e-9
844 Returns
845 -------
846 G : Complex ndarray, shape of *ham*
847 '''
848 error_handling.real_number(energy, 'energy')
849 error_handling.positive_real(eta, 'eta')
850 sigma = sum((self.self_energy(l, energy, eta) for l in range(len(self.leads))),
851 np.zeros_like(self.ham))
852 return LA.inv((energy + 1j * eta) * np.eye(len(self.ham)) - self.ham - sigma)
854 def transmission(
855 self, energies: ArrayLike, lead_in: int = 0, lead_out: int = 1, eta: float | None = None,
856 ) -> NDArray[np.float64]:
857 r'''
858 Get the transmission from lead *lead_in* to lead *lead_out*: the
859 two-terminal conductance in units of :math:`e^2/h` (per spin, for a
860 spinless model).
862 By default (*eta* None) it is :math:`\sum_{ab}|S_{ab}|^2` over the block of
863 the scattering matrix (*smatrix*), exact up to the band edges of the
864 leads. With a broadening *eta* (and, as a fallback, for non-Hermitian
865 leads or *lead_in* = *lead_out*, with *eta* = 1e-9) it is the Caroli
866 formula :math:`T(E) = \mathrm{Tr}[\Gamma_{out} G^r \Gamma_{in} G^a]`.
867 The two agree to :math:`O(\eta)`.
869 Parameters
870 ----------
871 energies : Real array
872 Energies.
873 lead_in : Integer, default 0
874 lead_out : Integer, default 1
875 eta : Positive real or None, default None
877 Returns
878 -------
879 T : Real ndarray, same length as *energies*
880 '''
881 error_handling.lead_index(lead_in, len(self.leads))
882 error_handling.lead_index(lead_out, len(self.leads))
883 energies = np.atleast_1d(np.asarray(energies, dtype='f8'))
884 if self._exact(eta) and lead_in != lead_out:
885 return np.array(self._map(lambda e: np.sum(np.abs(self._amplitudes(e, lead_in, lead_out)) ** 2),
886 energies), dtype='f8')
887 eta = 1e-9 if eta is None else eta
888 self.ham # built once, before the threads
890 def caroli(e):
891 sig_in = self.self_energy(lead_in, e, eta)
892 sig_out = self.self_energy(lead_out, e, eta)
893 gam_in = 1j * (sig_in - sig_in.conj().T)
894 gam_out = 1j * (sig_out - sig_out.conj().T)
895 g = self.get_green(e, eta)
896 return np.trace(gam_out @ g @ gam_in @ g.conj().T).real
897 return np.array(self._map(caroli, energies), dtype='f8')
899 # ------------------------------------------------------------------
900 # Scattering matrix and scattering states
901 # ------------------------------------------------------------------
903 def _exact(self, eta: float | None) -> bool:
904 r'''
905 Private method. Whether a method with a broadening *eta* uses the
906 scattering matrix and states (*eta* None and Hermitian leads) rather
907 than the dense Green's function.
908 '''
909 return eta is None and all(np.allclose(h0, h0.conj().T) for h0, _, _, _ in self.leads)
911 def _amplitudes(self, energy: float, lead_in: int, lead_out: int) -> NDArray[np.complex128]:
912 r'''
913 Private method. The block of the scattering matrix from lead
914 *lead_in* to lead *lead_out*, from the modes injected by *lead_in*
915 only: shape (outgoing modes of *lead_out*, incoming modes of *lead_in*).
916 '''
917 modes, sol = self._scattering(energy, [lead_in], slice(self._ham.shape[0], None))
918 a = sum(len(h0) for h0, _, _, _ in self.leads[:lead_out])
919 return sol[a:a + int(np.sum(modes[lead_out].velocities > 0))]
921 def _scattering(self, energy: float, inject: list[int], keep: slice = slice(None)) -> tuple:
922 r'''
923 Private method. Solve the scattering problem for every incoming mode
924 of the leads *inject* (Groth et al. 2014). In lead :math:`l`, the
925 wave function on cells 0 and 1 is the incoming mode plus
926 :math:`U_0 b_l`, :math:`U_1 b_l` (*_lead_modes*: outgoing and
927 evanescent). The unknowns are the device wave function
928 :math:`\psi` and the amplitudes :math:`b_l`, and the equations are
929 the Schrodinger equation on the device and on the surface cell of
930 every lead,
932 .. math::
934 (E - H)\psi - \sum_l\tau_l U_{0,l}b_l = \sum_l\tau_l\phi^{in}_l\, ,\qquad
935 -\tau_l^\dagger\psi + [(E - h_0)U_0 - vU_1]b_l = -(E - h_0 - \lambda_{in}v)\phi^{in}_l\, ,
937 one sparse system of size n_sites + sum_l m_l, factorized once
938 (*scipy.sparse.linalg.splu*) for every right-hand side. Only the rows
939 *keep* of the solution are stored (e.g. the lead amplitudes for the
940 scattering matrix).
942 Returns
943 -------
944 modes : List of **LeadModes**, one per lead
945 sol
946 Complex ndarray, shape (n_sites + sum_l m_l, n_in) before
947 *keep*: one column per incoming mode of the leads *inject*, in order.
948 '''
949 error_handling.real_number(energy, 'energy')
950 n = self._ham.shape[0]
951 found = [_lead_modes(h0, v, energy, law) for (h0, v, _, _), law in zip(self.leads, self._laws)]
952 size = n + sum(len(h0) for h0, _, _, _ in self.leads)
953 device = (energy * sp.eye(n, dtype='c16') - self._ham).tocoo()
954 coo = ([device.row], [device.col], [device.data])
955 rhs_coo: tuple = ([], [], [])
957 def block(r, c, values, out=coo):
958 rr, cc = np.meshgrid(r, c, indexing='ij')
959 out[0].append(rr.ravel())
960 out[1].append(cc.ravel())
961 out[2].append(values.ravel())
963 n_rhs, a = 0, n
964 for l, ((h0, v, tau, sites), (modes, n_in, u0, u1)) in enumerate(zip(self.leads, found)):
965 m = len(h0)
966 lead = np.arange(a, a + m)
967 e_h0 = energy * np.eye(m) - h0
968 block(sites, lead, -tau @ u0)
969 block(lead, sites, -tau.conj().T)
970 block(lead, lead, e_h0 @ u0 - v @ u1)
971 if l in inject:
972 phi = modes.wave_functions[:, :n_in]
973 cols = np.arange(n_rhs, n_rhs + n_in)
974 block(sites, cols, tau @ phi, rhs_coo)
975 block(lead, cols, -(e_h0 @ phi - v @ phi * np.exp(1j * modes.momenta[:n_in])), rhs_coo)
976 n_rhs += n_in
977 a += m
978 if n_rhs == 0:
979 return [f[0] for f in found], np.zeros((size, 0), 'c16')[keep]
980 mat, rhs = (sp.csc_matrix((np.concatenate(x), (np.concatenate(r), np.concatenate(c))), shape=(size, s))
981 for (r, c, x), s in ((coo, size), (rhs_coo, n_rhs)))
982 lu = spla.splu(mat)
983 # the right-hand sides are dense in blocks of at most 64 MB
984 chunk = max(1, 2**22 // size)
985 sol = np.hstack([lu.solve(rhs[:, c:c + chunk].toarray())[keep] for c in range(0, n_rhs, chunk)])
986 return [f[0] for f in found], sol
988 def smatrix(self, energy: float) -> SMatrix:
989 r'''
990 Get the scattering matrix at one energy: the amplitudes of the
991 outgoing modes of every lead for each incoming mode of unit
992 current, from the exact lead modes (*lead_modes*, no broadening) and
993 one sparse factorization of the device. It is unitary for a
994 Hermitian device. Every lead must be Hermitian. The blocks of the
995 leads' conservation laws (*add_lead*) are addressed as
996 ``(lead, block)`` in *SMatrix.submatrix* and *SMatrix.transmission*.
998 Parameters
999 ----------
1000 energy : Real number
1001 Energy, not at a band edge of a lead.
1003 Returns
1004 -------
1005 S : **SMatrix**: *data*, *lead_info*, *submatrix*, *transmission*
1007 Examples
1008 --------
1009 Example usage::
1011 s = tr.smatrix(0.3)
1012 t = s.submatrix(1, 0) # transmission amplitudes, lead 0 to lead 1
1013 T = s.transmission(1, 0) # = tr.transmission(0.3, 0, 1)
1014 '''
1015 modes, sol = self._scattering(energy, list(range(len(self.leads))), slice(self._ham.shape[0], None))
1016 rows, a = [], 0
1017 for (h0, _, _, _), mode in zip(self.leads, modes):
1018 rows.append(sol[a:a + int(np.sum(mode.velocities > 0))])
1019 a += len(h0)
1020 data = np.vstack(rows) if rows else np.zeros((0, 0), 'c16')
1021 return SMatrix(data, modes, [1 if law is None else len(law) for law in self._laws])
1023 def wave_function(self, energy: float, lead: int) -> NDArray[np.complex128]:
1024 r'''
1025 Get the scattering states injected from a lead: the wave function on
1026 the device of each incoming mode of unit current (*lead_modes*),
1027 with the scattered part leaving through every lead.
1029 Parameters
1030 ----------
1031 energy : Real number
1032 Energy, not at a band edge of a lead.
1033 lead : Integer
1034 Injecting lead.
1036 Returns
1037 -------
1038 psi
1039 Complex ndarray, shape (n_modes, n_sites): one row
1040 per incoming mode of the lead.
1041 '''
1042 error_handling.lead_index(lead, len(self.leads))
1043 _, sol = self._scattering(energy, [lead])
1044 return sol[:self._ham.shape[0]].T
1046 def ldos(self, energy: float) -> NDArray[np.float64]:
1047 r'''
1048 Get the local density of states carried by the scattering states,
1049 :math:`\rho_i(E) = \frac{1}{2\pi}\sum_{l,a}|\psi_{la}(i)|^2` over every
1050 incoming mode :math:`a` of every lead :math:`l` (unit current
1051 normalization). The bound states of the device, which no lead
1052 feeds, are left out.
1054 Parameters
1055 ----------
1056 energy : Real number
1057 Energy, not at a band edge of a lead.
1059 Returns
1060 -------
1061 rho : Real ndarray, shape (n_sites,)
1062 '''
1063 _, sol = self._scattering(energy, list(range(len(self.leads))))
1064 return np.sum(np.abs(sol[:self._ham.shape[0]]) ** 2, axis=1) / (2 * np.pi)
1066 # ------------------------------------------------------------------
1067 # Finite temperature: conductance and thermoelectric coefficients
1068 # ------------------------------------------------------------------
1070 def _onsager(
1071 self, mu: float, temperature: float, lead_in: int, lead_out: int, eta: float | None,
1072 n_points: int,
1073 ) -> tuple[float, float, float]:
1074 r'''
1075 Private method. The integrals :math:`L_n = \int dE\,(E - \mu)^n\,T(E)
1076 (-\partial f/\partial E)`, n = 0, 1, 2.
1077 '''
1078 energies, weights = _fermi_window(mu, temperature, n_points)
1079 t = self.transmission(energies, lead_in, lead_out, eta) * weights
1080 de = energies - mu
1081 return float(t.sum()), float(np.sum(t * de)), float(np.sum(t * de**2))
1083 def conductance(
1084 self, mu: float, temperature: float = 0., lead_in: int = 0, lead_out: int = 1,
1085 eta: float | None = None, n_points: int = 101,
1086 ) -> float:
1087 r'''
1088 Get the two-terminal linear conductance at a finite temperature, in
1089 units of :math:`e^2/h`,
1091 .. math::
1093 G = \int dE\,\left(-\frac{\partial f}{\partial E}\right)T(E)\, ,
1095 with the Fermi-Dirac :math:`f` at chemical potential :math:`\mu`. The
1096 integral runs over :math:`|E - \mu| \le 36k_BT` on *n_points* equally
1097 spaced energies (the trapezoidal rule, spectrally accurate for
1098 a smooth :math:`T(E)`). Resolving features of :math:`T(E)` narrower than
1099 :math:`k_BT` needs more points. At zero temperature it is :math:`T(\mu)`.
1101 Parameters
1102 ----------
1103 mu : Real number
1104 Chemical potential.
1105 temperature : Positive real or zero
1106 Default value 0 (:math:`k_B = 1`).
1107 lead_in : Integer, default 0
1108 lead_out : Integer, default 1
1109 eta : Positive real or None
1110 Default value None (see *transmission*).
1111 n_points : Positive integer, default 101
1113 Returns
1114 -------
1115 G : Real number
1116 '''
1117 error_handling.real_number(mu, 'mu')
1118 error_handling.positive_real_zero(temperature, 'temperature')
1119 error_handling.positive_int(n_points, 'n_points')
1120 if temperature == 0:
1121 return float(self.transmission([mu], lead_in, lead_out, eta)[0])
1122 return self._onsager(mu, temperature, lead_in, lead_out, eta, n_points)[0]
1124 def thermoelectric(
1125 self, mu: float, temperature: float, lead_in: int = 0, lead_out: int = 1,
1126 eta: float | None = None, n_points: int = 101,
1127 ) -> tuple[float, float, float]:
1128 r'''
1129 Get the two-terminal thermoelectric coefficients in linear response
1130 (Sivan and Imry 1986), from :math:`L_n = \int dE\,(E - \mu)^n\,T(E)
1131 (-\partial f/\partial E)`:
1133 .. math::
1135 G = L_0\, ,\qquad S = -\frac{L_1}{T L_0}\, ,\qquad
1136 \kappa = \frac{1}{T}\left(L_2 - \frac{L_1^2}{L_0}\right)\, ,
1138 the conductance (:math:`e^2/h`), the thermopower (:math:`k_B/e`, for
1139 electrons of charge :math:`-e`) and the electronic thermal conductance
1140 at zero current (:math:`k_B/h` times the energy unit). At low
1141 temperature they obey the Mott formula
1142 :math:`S = -\frac{\pi^2}{3}T\,\frac{d\ln T(E)}{dE}` and the
1143 Wiedemann-Franz law :math:`\kappa = \frac{\pi^2}{3}TG`. The integrals
1144 are those of *conductance*.
1146 Parameters
1147 ----------
1148 mu : Real number
1149 Chemical potential.
1150 temperature : Positive real (:math:`k_B = 1`)
1151 lead_in : Integer, default 0
1152 lead_out : Integer, default 1
1153 eta : Positive real or None
1154 Default value None (see *transmission*).
1155 n_points : Positive integer, default 101
1157 Returns
1158 -------
1159 G : Real number
1160 S : Real number (NaN when nothing is transmitted)
1161 kappa : Real number
1162 '''
1163 error_handling.real_number(mu, 'mu')
1164 error_handling.positive_real(temperature, 'temperature')
1165 error_handling.positive_int(n_points, 'n_points')
1166 l0, l1, l2 = self._onsager(mu, temperature, lead_in, lead_out, eta, n_points)
1167 if l0 < 1e-14:
1168 return l0, np.nan, l2 / temperature
1169 return l0, -l1 / (temperature * l0), (l2 - l1**2 / l0) / temperature
1171 # ------------------------------------------------------------------
1172 # Multi-terminal conductance, currents and noise
1173 # ------------------------------------------------------------------
1175 def _broadenings(self, energy: float, eta: float) -> list[NDArray[np.complex128]]:
1176 r'''
1177 Private method. The broadening :math:`\Gamma_l = i(\Sigma_l - \Sigma_l^\dagger)`
1178 of every lead, restricted to its sites (shape (len(sites), len(sites))).
1179 '''
1180 out = []
1181 for h0, v, tau, _ in self.leads:
1182 sig = tau @ surface_green(h0, v, energy, eta) @ tau.conj().T
1183 out.append(1j * (sig - sig.conj().T))
1184 return out
1186 def transmission_matrix(self, energy: float, eta: float | None = None) -> NDArray[np.float64]:
1187 r'''
1188 Get the transmissions :math:`T_{pq}` between every pair of leads at
1189 one energy, from lead :math:`q` to lead :math:`p` (the diagonal, a
1190 reflection, is set to 0). By default (*eta* None) they come from one
1191 scattering matrix (*smatrix*), :math:`T_{pq} = \sum_{ab}|S_{ab}|^2`
1192 over its block from :math:`q` to :math:`p`; with a broadening *eta*,
1193 from the Caroli formula :math:`\mathrm{Tr}[\Gamma_pG^r\Gamma_qG^a]`
1194 (see *transmission*).
1196 Parameters
1197 ----------
1198 energy : Real number
1199 eta : Positive real or None, default None
1201 Returns
1202 -------
1203 T : Real ndarray, shape (n_leads, n_leads)
1204 '''
1205 n = len(self.leads)
1206 if self._exact(eta):
1207 s = self.smatrix(energy)
1208 return np.array([[s.transmission(p, q) if p != q else 0. for q in range(n)] for p in range(n)])
1209 eta = 1e-9 if eta is None else eta
1210 g = self.get_green(energy, eta)
1211 gam = self._broadenings(energy, eta)
1212 out = np.zeros((n, n))
1213 for p in range(n):
1214 for q in range(n):
1215 if p != q:
1216 gpq = g[np.ix_(self.leads[p][3], self.leads[q][3])]
1217 out[p, q] = np.trace(gam[p] @ gpq @ gam[q] @ gpq.conj().T).real
1218 return out
1220 def conductance_matrix(
1221 self, energy: float, eta: float | None = None, temperature: float = 0., n_points: int = 101,
1222 ) -> NDArray[np.float64]:
1223 r'''
1224 Get the Landauer-Buttiker conductance matrix (Buttiker 1986), in units
1225 of :math:`e^2/h`: the linear-response currents flowing *into* the
1226 device from each lead, :math:`I_p = \sum_q G_{pq}V_q`, with
1228 .. math::
1230 G_{pq} = -T_{pq}\ (p\neq q)\, ,\qquad G_{pp} = \sum_{q\neq p}T_{qp}\, .
1232 Its rows and columns add up to zero (current conservation, and no
1233 current at equal voltages). At a finite *temperature* it is averaged
1234 over the Fermi window, :math:`\int dE\,(-\partial f/\partial E)\,G(E)`
1235 (see *conductance*). The transmissions come from *transmission_matrix*.
1237 Parameters
1238 ----------
1239 energy : Real number
1240 Fermi energy (chemical potential).
1241 eta : Positive real or None, default None
1242 See *transmission_matrix*.
1243 temperature : Positive real or zero
1244 Default value 0 (:math:`k_B = 1`).
1245 n_points : Positive integer, default 101
1246 Energies in the
1247 Fermi window, at a finite temperature.
1249 Returns
1250 -------
1251 G : Real ndarray, shape (n_leads, n_leads)
1252 '''
1253 error_handling.positive_real_zero(temperature, 'temperature')
1254 error_handling.positive_int(n_points, 'n_points')
1255 if temperature > 0:
1256 energies, weights = _fermi_window(energy, temperature, n_points)
1257 if not self._exact(eta):
1258 self.ham # built once, before the threads
1259 return sum(w * g for g, w in zip(self._map(lambda e: self.conductance_matrix(e, eta), energies),
1260 weights))
1261 t = self.transmission_matrix(energy, eta)
1262 return np.diag(t.sum(axis=0)) - t
1264 def four_terminal_resistance(
1265 self, energy: float, current: tuple[int, int], voltage: tuple[int, int], eta: float | None = None,
1266 temperature: float = 0., n_points: int = 101,
1267 ) -> float:
1268 r'''
1269 Get a four-terminal (or two-terminal) resistance from the
1270 Landauer-Buttiker equations: a current :math:`I` enters through lead
1271 ``current[0]`` and leaves through ``current[1]``, every other lead is
1272 an ideal voltage probe (no net current), and the result is
1273 :math:`R = (V_a - V_b)/I` with ``voltage = (a, b)``, in units of
1274 :math:`h/e^2`. On a Hall bar, a longitudinal pair gives
1275 :math:`R_{xx}` and a transverse pair :math:`R_{xy}`; in the quantum
1276 Hall regime with :math:`\nu` edge channels, :math:`|R_{xy}| = 1/\nu`
1277 and :math:`R_{xx} = 0` (Buttiker 1988).
1279 Parameters
1280 ----------
1281 energy : Real number
1282 Fermi energy.
1283 current : Tuple of two distinct lead indices (source, drain)
1284 voltage : Tuple of two lead indices (a, b)
1285 eta : Positive real or None, default None
1286 See *transmission_matrix*.
1287 temperature : Positive real or zero
1288 Default value 0 (see *conductance_matrix*).
1289 n_points : Positive integer, default 101
1291 Returns
1292 -------
1293 R : Real number
1294 '''
1295 error_handling.lead_pair(current, len(self.leads), 'current', True)
1296 error_handling.lead_pair(voltage, len(self.leads), 'voltage', False)
1297 g = self.conductance_matrix(energy, eta, temperature, n_points)
1298 source, drain = current
1299 keep = [p for p in range(len(self.leads)) if p != drain] # the drain is grounded
1300 currents = np.zeros(len(keep))
1301 currents[keep.index(source)] = 1.
1302 volts = np.zeros(len(self.leads))
1303 volts[keep] = LA.solve(g[np.ix_(keep, keep)], currents)
1304 return float(volts[voltage[0]] - volts[voltage[1]])
1306 def _bond_list(
1307 self, energy: float, lead: int, eta: float | None,
1308 ) -> tuple[NDArray[np.int64], NDArray[np.int64], NDArray[np.float64]]:
1309 r'''
1310 Private method. The bond currents of *bond_currents* on the bonds
1311 :math:`i \neq j` of the device (the nonzero :math:`H_{ij}`): rows,
1312 columns and currents. By default :math:`A^{(q)}_{ji}` comes from the
1313 scattering states, otherwise from the dense Green's function.
1314 '''
1315 error_handling.lead_index(lead, len(self.leads))
1316 h = self._ham.tocoo()
1317 bond = h.row != h.col
1318 rows, cols, hij = h.row[bond], h.col[bond], h.data[bond]
1319 if self._exact(eta):
1320 psi = self.wave_function(energy, lead)
1321 a_ji = np.sum(psi[:, cols] * psi[:, rows].conj(), axis=0)
1322 else:
1323 eta = 1e-9 if eta is None else eta
1324 g = self.get_green(energy, eta)
1325 sites = self.leads[lead][3]
1326 gam = self._broadenings(energy, eta)[lead]
1327 a_ji = np.sum((g[cols][:, sites] @ gam) * g[rows][:, sites].conj(), axis=1)
1328 return rows, cols, -2 * np.imag(hij * a_ji)
1330 def bond_currents(self, energy: float, lead: int, eta: float | None = None) -> NDArray[np.float64]:
1331 r'''
1332 Get the bond currents carried, per unit energy, by the electrons
1333 injected from one lead (the non-equilibrium current when that lead's
1334 electrochemical potential is raised),
1336 .. math::
1338 J_{i\to j} = -2\,\mathrm{Im}\left[H_{ij}\,A^{(q)}_{ji}\right]\, ,\qquad
1339 A^{(q)} = G^r\Gamma_qG^a = \sum_a\psi_a\psi_a^\dagger\, ,
1341 in units of :math:`e/h` per unit energy: summed over a cross-section
1342 of a two-terminal device, it is the transmission :math:`T(E)`. The
1343 matrix is antisymmetric, and the currents are conserved at every
1344 site not touched by a lead. By default (*eta* None) :math:`A^{(q)}`
1345 comes from the scattering states :math:`\psi_a` of the lead
1346 (*wave_function*, unit current), on the bonds of the sparse device
1347 only; with a broadening *eta*, from the dense Green's function.
1349 Parameters
1350 ----------
1351 energy : Real number
1352 lead : Integer
1353 Injecting lead.
1354 eta : Positive real or None, default None
1356 Returns
1357 -------
1358 J
1359 Real ndarray, shape of *ham*: J[i, j] is the current
1360 from site i to site j.
1361 '''
1362 rows, cols, j = self._bond_list(energy, lead, eta)
1363 n = self._ham.shape[0]
1364 return sp.csr_matrix((j, (rows, cols)), shape=(n, n)).toarray()
1366 def local_currents(
1367 self, energy: float, lead: int, positions: ArrayLike, eta: float | None = None,
1368 ) -> NDArray[np.float64]:
1369 r'''
1370 Get the current density vector on every site from the bond currents
1371 (*bond_currents*): half of the currents of the bonds leaving the
1372 site, along their bond vectors,
1373 :math:`\mathbf{j}_i = \frac12\sum_jJ_{i\to j}(\mathbf{r}_j - \mathbf{r}_i)`.
1374 By default (*eta* None) it never builds a dense matrix of the device.
1376 Parameters
1377 ----------
1378 energy : Real number
1379 lead : Integer
1380 Injecting lead.
1381 positions : Real array, shape (n_sites, d)
1382 Site positions
1383 (e.g. ``np.column_stack([lat.coor['x'], lat.coor['y']])``).
1384 eta : Positive real or None, default None
1386 Returns
1387 -------
1388 j : Real ndarray, shape (n_sites, d)
1389 '''
1390 positions = np.asarray(positions, dtype='f8')
1391 error_handling.site_positions(positions, self._ham.shape[0])
1392 rows, cols, bonds = self._bond_list(energy, lead, eta)
1393 out = np.zeros_like(positions)
1394 np.add.at(out, rows, 0.5 * bonds[:, None] * (positions[cols] - positions[rows]))
1395 return out
1397 def transmission_eigenvalues(
1398 self, energy: float, lead_in: int = 0, lead_out: int = 1, eta: float | None = None,
1399 ) -> NDArray[np.float64]:
1400 r'''
1401 Get the transmission eigenvalues :math:`T_n`, eigenvalues of
1402 :math:`t^\dagger t`: their sum is the transmission. By default (*eta*
1403 None) :math:`t` is the block of the scattering matrix from
1404 *lead_in* to *lead_out* (*smatrix*), and the list is padded with the
1405 zeros of the closed channels; with a broadening *eta* it is
1406 :math:`t = \Gamma_{out}^{1/2}G^r\Gamma_{in}^{1/2}` (and, as a
1407 fallback, with *eta* = 1e-9 for *lead_in* = *lead_out*).
1409 Parameters
1410 ----------
1411 energy : Real number
1412 lead_in : Integer, default 0
1413 lead_out : Integer, default 1
1414 eta : Positive real or None, default None
1416 Returns
1417 -------
1418 T_n
1419 Real ndarray, in decreasing order, of length
1420 min(len(sites_in), len(sites_out)).
1421 '''
1422 error_handling.lead_index(lead_in, len(self.leads))
1423 error_handling.lead_index(lead_out, len(self.leads))
1424 s_in, s_out = self.leads[lead_in][3], self.leads[lead_out][3]
1425 if self._exact(eta) and lead_in != lead_out:
1426 t = self._amplitudes(energy, lead_in, lead_out)
1427 vals = np.linalg.svd(t, compute_uv=False) ** 2 if t.size else np.zeros(0)
1428 # at most min(len(s_in), len(s_out)) channels transmit: the others are zeros
1429 out = np.zeros(min(len(s_in), len(s_out)))
1430 out[:min(len(out), len(vals))] = vals[:len(out)]
1431 return out
1432 eta = 1e-9 if eta is None else eta
1433 g = self.get_green(energy, eta)
1434 gam = self._broadenings(energy, eta)
1435 return _eigen_transmissions(gam[lead_out], g[np.ix_(s_out, s_in)], gam[lead_in])
1437 def shot_noise(
1438 self, energies: ArrayLike, lead_in: int = 0, lead_out: int = 1, eta: float | None = None,
1439 ) -> NDArray[np.float64]:
1440 r'''
1441 Get the zero-temperature shot-noise factor :math:`\sum_nT_n(1-T_n)`
1442 (Lesovik 1989; Buttiker 1990): the noise power at bias :math:`V` is
1443 :math:`S = 2\frac{e^3|V|}{h}\sum_nT_n(1-T_n)` (per spin).
1445 Parameters
1446 ----------
1447 energies : Real array
1448 Energies.
1449 lead_in : Integer, default 0
1450 lead_out : Integer, default 1
1451 eta : Positive real or None, default None
1452 See *transmission_eigenvalues*.
1454 Returns
1455 -------
1456 noise : Real ndarray, same length as *energies*
1457 '''
1458 error_handling.frequencies(energies, 'energies')
1459 energies = np.atleast_1d(np.asarray(energies, dtype='f8'))
1460 if not (self._exact(eta) and lead_in != lead_out):
1461 self.ham # built once, before the threads
1462 return np.array(self._map(lambda e: _noise(self.transmission_eigenvalues(e, lead_in, lead_out, eta)),
1463 energies), dtype='f8')
1465 def fano_factor(
1466 self, energies: ArrayLike, lead_in: int = 0, lead_out: int = 1, eta: float | None = None,
1467 ) -> NDArray[np.float64]:
1468 r'''
1469 Get the Fano factor of the zero-temperature shot noise,
1471 .. math::
1473 F = \frac{\sum_nT_n(1-T_n)}{\sum_nT_n}\, ,
1475 the noise in units of the Poisson value :math:`2e|I|`: 0 for a
1476 ballistic conductor (every :math:`T_n` is 0 or 1), close to 1 for a
1477 tunnel barrier, 1/3 for a diffusive wire (Beenakker and Buttiker 1992).
1479 Parameters
1480 ----------
1481 energies : Real array
1482 Energies.
1483 lead_in : Integer, default 0
1484 lead_out : Integer, default 1
1485 eta : Positive real or None, default None
1486 See *transmission_eigenvalues*.
1488 Returns
1489 -------
1490 F
1491 Real ndarray, same length as *energies* (NaN where
1492 nothing is transmitted).
1493 '''
1494 error_handling.frequencies(energies, 'energies')
1495 energies = np.atleast_1d(np.asarray(energies, dtype='f8'))
1496 if not (self._exact(eta) and lead_in != lead_out):
1497 self.ham # built once, before the threads
1498 return np.array(self._map(lambda e: _fano(self.transmission_eigenvalues(e, lead_in, lead_out, eta)),
1499 energies), dtype='f8')
1502def _fermi_window(
1503 mu: float, temperature: float, n_points: int,
1504) -> tuple[NDArray[np.float64], NDArray[np.float64]]:
1505 r'''
1506 Private function. Energies and weights of the trapezoidal rule for
1507 :math:`\int dE\,(-\partial f/\partial E)\,F(E)`, on *n_points* equally
1508 spaced energies with :math:`|E - \mu| \le 36T` (the weight is below
1509 1e-15 beyond). The weights add up to 1.
1510 '''
1511 x = np.linspace(-36., 36., n_points) if n_points > 1 else np.zeros(1)
1512 w = 1 / (4 * np.cosh(x / 2) ** 2)
1513 return mu + temperature * x, w / w.sum()
1516def _sqrt_psd(mat: NDArray[np.complex128]) -> NDArray[np.complex128]:
1517 '''
1518 Private function. Square root of a positive semi-definite Hermitian matrix.
1519 '''
1520 w, u = np.linalg.eigh(0.5 * (mat + mat.conj().T))
1521 return (u * np.sqrt(np.clip(w, 0., None))) @ u.conj().T
1524def _eigen_transmissions(
1525 gam_out: NDArray[np.complex128], g: NDArray[np.complex128], gam_in: NDArray[np.complex128],
1526) -> NDArray[np.float64]:
1527 r'''
1528 Private function. Transmission eigenvalues: the squared singular values
1529 of :math:`\Gamma_{out}^{1/2}G\,\Gamma_{in}^{1/2}` (*g* the block of the
1530 Green's function from the input to the output sites).
1531 '''
1532 t = _sqrt_psd(gam_out) @ g @ _sqrt_psd(gam_in)
1533 return np.linalg.svd(t, compute_uv=False) ** 2
1536def _noise(t: NDArray[np.float64]) -> float:
1537 '''
1538 Private function. Shot-noise factor sum T_n (1 - T_n).
1539 '''
1540 return float(np.sum(t * (1 - t)))
1543def _fano(t: NDArray[np.float64]) -> float:
1544 '''
1545 Private function. Fano factor from transmission eigenvalues (NaN if none transmit).
1546 '''
1547 total = t.sum()
1548 return _noise(t) / float(total) if total > 1e-12 else np.nan
1551def slices_from_positions(x: ArrayLike, tol: float = 1e-6) -> list[list[int]]:
1552 r'''
1553 Group the sites of a device into slices of equal coordinate *x* (e.g.
1554 the columns of a strip along :math:`x`), ordered by increasing *x*, for
1555 *RecursiveTransport*.
1557 Parameters
1558 ----------
1559 x : Real array
1560 One coordinate per site (e.g. ``lat.coor['x']``).
1561 tol : Positive real, default 1e-6
1562 Coordinates closer than
1563 *tol* belong to the same slice.
1565 Returns
1566 -------
1567 slices : List of lists of site indices
1568 '''
1569 x = np.asarray(x, dtype='f8').ravel()
1570 error_handling.ndarray_empty(x, 'x')
1571 error_handling.positive_real(tol, 'tol')
1572 order = np.argsort(x, kind='stable')
1573 breaks = np.flatnonzero(np.diff(x[order]) > tol) + 1
1574 return [sorted(s.tolist()) for s in np.split(order, breaks)]
1577class RecursiveTransport(_EnergyWorkers):
1578 r'''
1579 Two-terminal transport through a long quasi-1D device by the recursive
1580 Green's function method (Thouless and Kirkpatrick 1981; Lee and Fisher
1581 1981; MacKinnon 1985). The device is cut into slices :math:`1\dots N`,
1582 each coupled only to its neighbours, with the left lead on slice 1 and
1583 the right lead on slice :math:`N`. The left-connected Green's functions
1585 .. math::
1587 g_1 = [z - H_{11} - \Sigma_L]^{-1}\, ,\qquad
1588 g_{i} = [z - H_{ii} - H_{i,i-1}\,g_{i-1}H_{i-1,i}]^{-1}\, ,
1590 (:math:`\Sigma_R` added in the last one, :math:`z = E + i\eta`) and the
1591 propagator :math:`G_{i1} = g_iH_{i,i-1}G_{i-1,1}` give
1592 :math:`T = \mathrm{Tr}[\Gamma_RG_{N1}\Gamma_LG_{N1}^\dagger]`. The cost
1593 grows linearly with the length, and only slice-sized matrices are
1594 inverted, never the whole device (compare *Transport*, which inverts
1595 the dense device matrix).
1597 Parameters
1598 ----------
1599 ham : Square matrix, sparse (preferred for long devices) or dense
1600 Device Hamiltonian.
1601 slices
1602 List of lists of site indices: a partition of the device
1603 into slices, ordered from the left lead to the right one, with
1604 hoppings between neighbouring slices only (see *slices_from_positions*).
1605 left
1606 Tuple (h0, v, coupling): the left lead (see *Transport.add_lead*),
1607 with a coupling of shape (len(slices[0]), len(h0)), its rows in the
1608 order of ``slices[0]``.
1609 right : Tuple (h0, v, coupling): the right lead, coupled to ``slices[-1]``
1611 Examples
1612 --------
1613 Example usage::
1615 slices = slices_from_positions(lat.coor['x'])
1616 rgf = RecursiveTransport(sys.ham, slices, (h0, v_left, tau), (h0, v_right, tau))
1617 T = rgf.transmission(energies)
1618 '''
1620 def __init__(self, ham, slices: list[list[int]], left: tuple, right: tuple) -> None:
1621 ham = sp.csr_matrix(ham, dtype='c16')
1622 error_handling.square_matrix(ham, 'ham')
1623 error_handling.slices(slices, ham.shape[0])
1624 label = np.empty(ham.shape[0], dtype=int)
1625 pos = np.empty(ham.shape[0], dtype=int)
1626 for n, s in enumerate(slices):
1627 label[s] = n
1628 pos[s] = np.arange(len(s))
1629 coo = ham.tocoo()
1630 coo.sum_duplicates()
1631 lr, lc = label[coo.row], label[coo.col]
1632 error_handling.block_tridiagonal(int(np.abs(lr - lc).max(initial=0)))
1633 self.slices = [list(s) for s in slices]
1634 sizes = [len(s) for s in self.slices]
1635 self.blocks = [np.zeros((m, m), 'c16') for m in sizes]
1636 self.couplings = [np.zeros((m, k), 'c16') for m, k in zip(sizes[:-1], sizes[1:])]
1637 # scatter every nonzero entry into its block (lc == lr) or coupling
1638 # (lc == lr + 1), grouped by slice; the lower couplings are their adjoints
1639 for target, keep in ((self.blocks, lr == lc), (self.couplings, lc == lr + 1)):
1640 order = np.argsort(lr[keep], kind='stable')
1641 rows, cols = pos[coo.row[keep]][order], pos[coo.col[keep]][order]
1642 vals, owner = coo.data[keep][order], lr[keep][order]
1643 bounds = np.searchsorted(owner, np.arange(len(target) + 1))
1644 for n, mat in enumerate(target):
1645 a, b = bounds[n], bounds[n + 1]
1646 mat[rows[a:b], cols[a:b]] = vals[a:b]
1647 self.leads = []
1648 for lead, s in ((left, self.slices[0]), (right, self.slices[-1])):
1649 error_handling.lead_tuple(lead)
1650 h0 = np.atleast_2d(np.asarray(lead[0], dtype='c16'))
1651 v = np.atleast_2d(np.asarray(lead[1], dtype='c16'))
1652 error_handling.lead(h0, v)
1653 tau = np.atleast_2d(np.asarray(lead[2], dtype='c16'))
1654 error_handling.lead_coupling(tau, list(range(len(s))), len(h0), len(s))
1655 self.leads.append((h0, v, tau))
1657 def _propagator(self, energy: float, eta: float) -> tuple:
1658 r'''
1659 Private method. :math:`G_{N1}` and the broadenings of the two leads.
1660 '''
1661 sig = [tau @ surface_green(h0, v, energy, eta) @ tau.conj().T for h0, v, tau in self.leads]
1662 gam = [1j * (s - s.conj().T) for s in sig]
1663 z = energy + 1j * eta
1664 last = len(self.blocks) - 1
1665 g = g_n1 = None
1666 for i, h in enumerate(self.blocks):
1667 a = z * np.eye(len(h)) - h
1668 if i == 0:
1669 a = a - sig[0]
1670 else:
1671 v = self.couplings[i - 1]
1672 a = a - v.conj().T @ g @ v
1673 if i == last:
1674 a = a - sig[1]
1675 g = LA.inv(a)
1676 g_n1 = g if i == 0 else g @ self.couplings[i - 1].conj().T @ g_n1
1677 return g_n1, gam
1679 def transmission(self, energies: ArrayLike, eta: float = 1e-9) -> NDArray[np.float64]:
1680 r'''
1681 Get the transmission from the left lead to the right one: the same
1682 as *Transport.transmission* for the same device and leads.
1684 Parameters
1685 ----------
1686 energies : Real array
1687 Energies.
1688 eta : Positive real, default 1e-9
1690 Returns
1691 -------
1692 T : Real ndarray, same length as *energies*
1693 '''
1694 error_handling.frequencies(energies, 'energies')
1695 energies = np.atleast_1d(np.asarray(energies, dtype='f8'))
1696 error_handling.positive_real(eta, 'eta')
1698 def caroli(e):
1699 g, (gam_l, gam_r) = self._propagator(e, eta)
1700 return np.trace(gam_r @ g @ gam_l @ g.conj().T).real
1701 return np.array(self._map(caroli, energies), dtype='f8')
1703 def transmission_eigenvalues(self, energy: float, eta: float = 1e-9) -> NDArray[np.float64]:
1704 r'''
1705 Get the transmission eigenvalues (see *Transport.transmission_eigenvalues*).
1707 Parameters
1708 ----------
1709 energy : Real number
1710 eta : Positive real, default 1e-9
1712 Returns
1713 -------
1714 T_n : Real ndarray, in decreasing order
1715 '''
1716 error_handling.real_number(energy, 'energy')
1717 error_handling.positive_real(eta, 'eta')
1718 g, (gam_l, gam_r) = self._propagator(energy, eta)
1719 return _eigen_transmissions(gam_r, g, gam_l)
1721 def fano_factor(self, energies: ArrayLike, eta: float = 1e-9) -> NDArray[np.float64]:
1722 r'''
1723 Get the Fano factor of the shot noise (see *Transport.fano_factor*).
1725 Parameters
1726 ----------
1727 energies : Real array
1728 Energies.
1729 eta : Positive real, default 1e-9
1731 Returns
1732 -------
1733 F : Real ndarray, same length as *energies*
1734 '''
1735 error_handling.frequencies(energies, 'energies')
1736 energies = np.atleast_1d(np.asarray(energies, dtype='f8'))
1737 return np.array(self._map(lambda e: _fano(self.transmission_eigenvalues(e, eta)), energies), dtype='f8')