Coverage for tbkit/bdg.py: 100%
130 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-15 06:13 +0000
1r"""
2Superconductivity in the Bogoliubov-de Gennes (BdG) formalism. With a
3pairing term :math:`\sum_{ij}\Delta_{ij}c^\dagger_ic^\dagger_j + h.c.`,
4the Hamiltonian is :math:`\frac12\Psi^\dagger H_{BdG}\Psi` in the Nambu
5basis :math:`\Psi = (c_1, \dots, c_N, c_1^\dagger, \dots, c_N^\dagger)^T`,
7.. math::
9 H_{BdG} = \begin{pmatrix} h - \mu & \Delta \\ \Delta^\dagger & -(h-\mu)^T \end{pmatrix}\, ,
11with :math:`\Delta^T = -\Delta` (Fermi statistics). Its spectrum is
12symmetric, :math:`\pm E_n`, by the particle-hole symmetry
13:math:`\tau_x H_{BdG}^*\tau_x = -H_{BdG}`; the quasiparticle energies are
14the :math:`E_n \ge 0`, and a zero mode is a Majorana mode. Real-space
15models come from *bdg_ham*, periodic ones from *bdg_kspace* (a
16:class:`tbkit.kspace.KSpace` with the holes as extra orbitals, so its bands,
17Berry phases, *finite_ham*, ... all apply).
18"""
19from __future__ import annotations
21import numpy as np
22from numpy.typing import NDArray
23import scipy.sparse as sparse
25import tbkit.error_handling as error_handling
26from tbkit.lattice import Lattice
29def bdg_ham(ham, pairing, mu: float = 0.) -> sparse.csr_matrix:
30 r'''
31 Get the real-space BdG Hamiltonian (see the module docstring).
33 :param ham: Square Hermitian matrix, shape (N, N) (e.g. *System.ham*).
34 :param pairing: Square antisymmetric matrix, shape (N, N): :math:`\Delta`
35 (see *pairing_bonds*, *pairing_s_wave*).
36 :param mu: Real number. Default value 0. Chemical potential.
38 :returns:
39 * **ham_bdg** -- Sparse CSR matrix, shape (2N, 2N).
40 '''
41 h = sparse.csr_matrix(ham, dtype='c16')
42 delta = sparse.csr_matrix(pairing, dtype='c16')
43 error_handling.square_matrix(h, 'ham')
44 error_handling.hermitian(h)
45 error_handling.pairing(delta, h.shape[0])
46 error_handling.real_number(mu, 'mu')
47 h = h - mu * sparse.identity(h.shape[0], dtype='c16', format='csr')
48 return sparse.csr_matrix(sparse.bmat([[h, delta], [delta.conj().T, -h.T]]))
51def pairing_bonds(sys, delta: complex) -> sparse.csr_matrix:
52 r'''
53 Get a p-wave (odd) pairing on the bonds of *sys.hop*:
54 :math:`\Delta_{ij} = \Delta`, :math:`\Delta_{ji} = -\Delta` for every
55 stored hopping :math:`i\to j` (oriented as *System.set_hopping* stores
56 them: angle in :math:`[0, 180)`). On a chain it is Kitaev's
57 :math:`\Delta\,c^\dagger_jc^\dagger_{j+1}`.
59 :param sys: **System** instance, with hoppings.
60 :param delta: Complex number. :math:`\Delta`.
62 :returns:
63 * **pairing** -- Sparse CSR matrix, shape (sites, sites).
64 '''
65 error_handling.empty_hop(sys.hop)
66 error_handling.number(delta, 'delta')
67 i, j = sys.hop['i'].astype(int), sys.hop['j'].astype(int)
68 n = sys.lat.sites
69 vals = np.full(len(i), delta, dtype='c16')
70 return sparse.csr_matrix((np.concatenate([vals, -vals]), (np.concatenate([i, j]), np.concatenate([j, i]))),
71 shape=(n, n))
74def pairing_s_wave(n_sites: int, delta: complex) -> sparse.csr_matrix:
75 r'''
76 Get an onsite singlet (s-wave) pairing
77 :math:`\Delta\sum_i c^\dagger_{i\uparrow}c^\dagger_{i\downarrow} + h.c.`
78 for a spinful model with rows ordered site by site, spin up then down
79 (as *KSpace* with ``spin=True`` and *OrbitalSystem* with one orbital).
81 :param n_sites: Positive integer. Number of sites.
82 :param delta: Complex number. :math:`\Delta`.
84 :returns:
85 * **pairing** -- Sparse CSR matrix, shape (2 n_sites, 2 n_sites).
86 '''
87 error_handling.positive_int(n_sites, 'n_sites')
88 error_handling.number(delta, 'delta')
89 block = np.array([[0., delta], [-delta, 0.]], dtype='c16')
90 return sparse.csr_matrix(sparse.kron(sparse.identity(n_sites), block))
93def bdg_kspace(ks, pairing: list[dict], mu: float = 0.):
94 r'''
95 Get the BdG Bloch Hamiltonian of a periodic model,
97 .. math::
99 H_{BdG}(\mathbf{k}) = \begin{pmatrix} h(\mathbf{k})-\mu & \Delta(\mathbf{k}) \\
100 \Delta(\mathbf{k})^\dagger & -[h(-\mathbf{k})-\mu]^T\end{pmatrix}\, ,\qquad
101 \Delta(\mathbf{k}) = \sum_{\mathbf{R}}\Delta_{ij}(\mathbf{R})\,e^{i\mathbf{k}\cdot\mathbf{R}}\, ,
103 as a new **KSpace** whose orbitals are those of *ks* (particles,
104 ``0 ... norb-1``) followed by their holes (``norb ... 2 norb-1``, at the
105 same positions). Each pairing :math:`\Delta\,c^\dagger_{i,\mathbf{0}}c^\dagger_{j,\mathbf{R}} + h.c.`
106 enters antisymmetrized, :math:`\Delta_{ij}(\mathbf{R}) = -\Delta_{ji}(-\mathbf{R})`.
108 :param ks: **KSpace** instance (Hermitian, orthogonal basis).
109 :param pairing: List of dictionaries with keys ('i', 'j', 'R', 'delta'):
110 orbital indices of *ks* (``0 ... norb-1``; with ``spin=True``,
111 ``2*site + spin``), lattice vector, amplitude.
112 :param mu: Real number. Default value 0. Chemical potential.
114 :returns:
115 * **ks_bdg** -- **KSpace** instance, 2 norb orbitals.
117 Example usage::
119 # Kitaev chain: E(k) = +-sqrt((2t cos k - mu)^2 + 4 Delta^2 sin^2 k)
120 chain = KSpace(lattices.chain())
121 chain.set_hopping([{'i': 0, 'j': 0, 'R': (1,), 't': t}])
122 kit = bdg_kspace(chain, [{'i': 0, 'j': 0, 'R': (1,), 'delta': delta}], mu=mu)
123 '''
124 from tbkit.kspace import KSpace
125 error_handling.hermitian_kspace(ks.is_hermitian(), ks._overlap_hop)
126 error_handling.pairing_kspace(pairing, ks.norb, ks.dim)
127 error_handling.real_number(mu, 'mu')
128 norb = ks.norb
129 tau = ks.orbital_positions()
130 tags = np.repeat(ks.tags, 2) if ks.spin else ks.tags
131 cell = [{'tag': str(t), 'r0': tuple(float(c) for c in r)} for t, r in zip(tags, tau)]
132 new = KSpace(Lattice(unit_cell=cell + [dict(c) for c in cell], prim_vec=ks.lat.prim_vec))
133 h0 = np.diag(ks.onsite) + ks._onsite_offdiag - mu * np.eye(norb)
134 new.onsite = np.concatenate([np.diag(h0), -np.diag(h0)]).astype('c16')
135 off = h0 - np.diag(np.diag(h0))
136 new._onsite_offdiag = np.block([[off, np.zeros((norb, norb))], [np.zeros((norb, norb)), -off.T]]).astype('c16')
137 hops = [(i, j, R, t) for i, j, R, t in ks._hop]
138 hops += [(j + norb, i + norb, -R, -t) for i, j, R, t in ks._hop]
139 for dic in pairing:
140 R = np.zeros(ks.space_dim)
141 for n, a in zip(dic['R'], ks.lat.prim_vec):
142 R += n * np.array(a)
143 i, j, d = dic['i'], dic['j'], complex(dic['delta'])
144 hops += [(i, j + norb, R, d), (j, i + norb, -R, -d),
145 (j + norb, i, -R, np.conj(d)), (i + norb, j, R, -np.conj(d))]
146 new._hop = hops
147 return new
150def particle_hole(n: int) -> NDArray[np.float64]:
151 r'''
152 Get the unitary part :math:`\tau_x` of the BdG particle-hole symmetry
153 :math:`C = \tau_x\mathcal{K}` (:math:`C^2 = +1`), for *n* particle
154 orbitals: pass it to *KSpace.symmetry_error* (``antiunitary=True``,
155 ``anti=True``) or *KSpace.tenfold_class* (``particle_hole=``).
157 :param n: Positive integer. Number of particle orbitals.
159 :returns:
160 * **tau_x** -- Real ndarray, shape (2n, 2n).
161 '''
162 error_handling.positive_int(n, 'n')
163 return np.kron(np.array([[0., 1.], [1., 0.]]), np.eye(n))
166def pfaffian(a) -> complex:
167 r'''
168 Get the Pfaffian of an antisymmetric matrix,
169 :math:`\mathrm{Pf}(A)^2 = \det A`, by the Parlett-Reid reduction to
170 tridiagonal form with pivoting (Wimmer, ACM Trans. Math. Softw. 38, 30
171 (2012)). For a 2x2 block :math:`\begin{pmatrix}0 & a\\ -a & 0\end{pmatrix}`
172 it is :math:`a`; it is 0 for an odd dimension.
174 :param a: Square antisymmetric matrix, shape (n, n).
176 :returns:
177 * **pf** -- Complex number (a float for a real matrix).
178 '''
179 real = not np.iscomplexobj(a)
180 a = np.array(a, dtype='c16')
181 error_handling.square_matrix(a, 'a')
182 error_handling.antisymmetric(a)
183 n = len(a)
184 if n % 2:
185 return 0.
186 pf = 1. + 0j
187 for k in range(0, n - 1, 2):
188 # bring the largest entry of column k below the diagonal to row k + 1
189 p = k + 1 + int(np.argmax(np.abs(a[k + 1:, k])))
190 if p != k + 1:
191 a[[k + 1, p], :] = a[[p, k + 1], :]
192 a[:, [k + 1, p]] = a[:, [p, k + 1]]
193 pf = -pf
194 if a[k + 1, k] == 0.:
195 return 0.
196 pf *= a[k, k + 1]
197 if k + 2 < n:
198 tau = a[k, k + 2:] / a[k, k + 1]
199 a[k + 2:, k + 2:] += np.outer(tau, a[k + 2:, k + 1]) - np.outer(a[k + 2:, k + 1], tau)
200 return float(pf.real) if real else complex(pf)
203def majorana_number(ks) -> int:
204 r'''
205 Get Kitaev's Majorana number of a BdG model (Kitaev, Physics-Uspekhi
206 44, 131 (2001)): the product, over the time-reversal-invariant momenta
207 :math:`\Gamma_i`, of the signs of the Pfaffians of the Hamiltonian in
208 the Majorana basis,
210 .. math::
212 \mathcal{M} = \prod_i \mathrm{sign}\,\mathrm{Pf}\,A(\Gamma_i)\, ,\qquad
213 A = -i\,U H_{BdG} U^\dagger\, ,\qquad
214 U = \frac{1}{\sqrt2}\begin{pmatrix}1 & 1\\ -i & i\end{pmatrix}\, ,
216 :math:`U` taking the Nambu spinor :math:`(c, c^\dagger)` to the Majorana
217 operators :math:`(c + c^\dagger, -i(c - c^\dagger))`. Particle-hole
218 symmetry makes :math:`A(\Gamma_i)` real and antisymmetric. In 1D,
219 :math:`\mathcal{M} = \mathrm{sign}\,\mathrm{Pf}A(0)\,\mathrm{Pf}A(\pi) = -1`
220 in the topological phase (a Majorana zero mode at each end of an open
221 chain), +1 in the trivial one: the :math:`\mathbb{Z}_2` invariant of
222 class D. In 2D the product over the four momenta is
223 :math:`(-1)^C`, the parity of the Chern number (Ghosh et al., Phys. Rev.
224 B 82, 184525 (2010)). The invariant changes only where the gap closes
225 at one of the :math:`\Gamma_i`.
227 :param ks: **KSpace** instance of a BdG model, 2n orbitals: n particles
228 then their n holes (see *bdg_kspace*).
230 :returns:
231 * **m** -- Integer, -1 (topological) or +1 (trivial).
232 '''
233 from tbkit.kspace import KSpace
234 error_handling.kspace(ks, KSpace)
235 error_handling.bdg_orbitals(ks.norb)
236 n = ks.norb // 2
237 u = np.kron(np.array([[1., 1.], [-1j, 1j]]) / np.sqrt(2), np.eye(n))
238 sign = 1
239 for k in ks._trims():
240 a = -1j * u @ ks.get_ham(k) @ u.conj().T
241 error_handling.majorana_form(a)
242 pf = pfaffian(a.real)
243 error_handling.gapped_trim(pf)
244 sign *= int(np.sign(pf))
245 return sign
248class GapResult():
249 r'''
250 Self-consistent s-wave pairing of *s_wave_gap*.
252 :ivar delta: Complex ndarray. Gap :math:`\Delta_i` on every site.
253 :ivar energies: Real ndarray, shape (2N,). Spectrum of the reduced BdG
254 Hamiltonian (symmetric, :math:`\pm E_n`).
255 :ivar states: Complex ndarray, shape (2N, 2N). Its eigenvectors
256 :math:`(u_n, v_n)` (columns).
257 :ivar iterations: Integer. Number of iterations.
258 '''
260 def __init__(self, delta, energies, states, iterations) -> None:
261 self.delta = delta
262 self.energies = energies
263 self.states = states
264 self.iterations = iterations
266 @property
267 def mean_gap(self) -> float:
268 r'''
269 Average :math:`|\Delta_i|` over the sites.
270 '''
271 return float(np.mean(np.abs(self.delta)))
274def s_wave_gap(
275 ham, V: float | NDArray, mu: float = 0., temperature: float = 0., delta0: complex | NDArray = 0.1,
276 mixing: float = 1., tol: float = 1e-8, max_iter: int = 10000,
277) -> GapResult:
278 r'''
279 Solve the real-space gap equation of an s-wave superconductor
280 self-consistently (the attractive Hubbard model
281 :math:`-\sum_iV_in_{i\uparrow}n_{i\downarrow}` in the BCS mean-field
282 approximation, de Gennes 1966). For a spin-independent *ham* the singlet
283 BdG Hamiltonian reduces to the Nambu basis
284 :math:`(c_{i\uparrow}, c^\dagger_{i\downarrow})`,
286 .. math::
288 H = \begin{pmatrix} h-\mu & \Delta\\ \Delta^* & -(h-\mu)^*\end{pmatrix}\, ,
289 \qquad \Delta = \mathrm{diag}(\Delta_i)\, ,
291 whose eigenvectors :math:`(u_n, v_n)`, :math:`E_n > 0`, give the gap
292 equation
294 .. math::
296 \Delta_i = V_i\,\langle c_{i\downarrow}c_{i\uparrow}\rangle
297 = V_i\sum_{E_n>0}u_n(i)\,v_n^*(i)\tanh\frac{E_n}{2T}\, ,
299 iterated (with linear mixing) until self-consistent. On a clean lattice
300 it is the BCS gap equation :math:`1 = \frac{V}{N}\sum_{\mathbf{k}}
301 \tanh(E_{\mathbf{k}}/2T)/2E_{\mathbf{k}}`, with
302 :math:`\Delta(T_c) = 0` at :math:`T_c \approx \Delta(0)/1.764` at weak
303 coupling. The chemical potential is fixed (the Hartree shift
304 :math:`-V_in_i/2` is left out, i.e. absorbed in :math:`\mu`). Near
305 :math:`T_c` the iteration slows down (the linearized map has an
306 eigenvalue close to 1): raise *max_iter* there.
308 :param ham: Square Hermitian matrix, (N, N), spin independent (e.g. *System.ham*).
309 :param V: Positive real, or array of N non-negative reals. Attraction.
310 :param mu: Real number. Default value 0. Chemical potential.
311 :param temperature: Positive real or zero. Default value 0.
312 :param delta0: Number, or array of N numbers. Default value 0.1.
313 Initial gap.
314 :param mixing: Real in (0, 1]. Default value 1 (plain iteration).
315 :param tol: Positive real. Default value 1e-8. Convergence threshold
316 on the largest change of :math:`\Delta_i`.
317 :param max_iter: Positive integer. Default value 10000.
319 :returns:
320 * **result** -- :class:`GapResult`.
321 '''
322 h = ham.toarray() if hasattr(ham, 'toarray') else np.asarray(ham)
323 error_handling.square_matrix(h, 'ham')
324 error_handling.hermitian_dense(h)
325 n = len(h)
326 error_handling.interaction(V, n)
327 error_handling.real_number(mu, 'mu')
328 error_handling.positive_real_zero(temperature, 'temperature')
329 error_handling.pairing_amplitudes(delta0, n)
330 error_handling.mixing(mixing)
331 error_handling.positive_real(tol, 'tol')
332 error_handling.positive_int(max_iter, 'max_iter')
333 V = np.broadcast_to(np.asarray(V, dtype='f8'), (n,))
334 delta = np.broadcast_to(np.asarray(delta0, dtype='c16'), (n,)).copy()
335 xi = h.astype('c16') - mu * np.eye(n)
336 big = np.zeros((2 * n, 2 * n), 'c16')
337 big[:n, :n], big[n:, n:] = xi, -xi.conj()
338 diag = np.arange(n)
339 for it in range(1, max_iter + 1):
340 big[diag, diag + n] = delta
341 big[diag + n, diag] = delta.conj()
342 en, vec = np.linalg.eigh(big)
343 pos = en > 0
344 weight = np.tanh(en[pos] / (2 * temperature)) if temperature > 0 else np.ones(pos.sum())
345 new = V * np.einsum('in,in,n->i', vec[:n, pos], vec[n:, pos].conj(), weight)
346 change = np.max(np.abs(new - delta))
347 delta = (1 - mixing) * delta + mixing * new
348 if change < tol:
349 break
350 else:
351 error_handling.converged(False, 's_wave_gap')
352 return GapResult(new, en, vec, it)