physicskit.condensed #

Returns:

ndarray, shape (n_sites, n_sites) – Real, symmetric Hamiltonian matrix.

Examples

>>> H = anderson_chain_hamiltonian(n_sites=4, disorder_strength=0.0)
>>> H
array([[ 0., -1.,  0.,  0.],
       [-1.,  0., -1.,  0.],
       [ 0., -1.,  0., -1.],
       [ 0.,  0., -1.,  0.]])

physicskit.condensed.apply_peierls_phase(positions, hoppings, flux_quanta_per_plaquette, area_per_plaquette)[source]#

Apply a Peierls substitution phase to real-space hoppings for a uniform magnetic field.

In the Landau gauge \(\mathbf{A} = (-By, 0)\), a hopping from site at r_i to site at r_j acquires the phase \(\exp\!\big(i\frac{2\pi\Phi}{\Phi_0}\,\bar{y}\,(x_j-x_i)/a^2\big)\) where \(\Phi/\Phi_0\) is the flux per plaquette in units of the flux quantum, evaluated at the bond midpoint \(\bar y\).

Parameters:

positions (ndarray) – Real-space Cartesian coordinates of every site in a finite lattice.
hoppings (list of tuple(int, int, complex)) – (i, j, amplitude) real-space bonds (indices into positions).
flux_quanta_per_plaquette (float) – Magnetic flux per unit-cell plaquette, in units of the flux quantum \(\Phi_0 = h/e\).
area_per_plaquette (float) – Real-space area of one plaquette, used to convert flux density to the vector potential prefactor.

Returns:

list of tuple(int, int, complex) – The same bonds with Peierls phases multiplied into the amplitude.

Examples

>>> import numpy as np
>>> positions = np.array([[0.0, 0.0], [1.0, 0.0]])
>>> bonds = [(0, 1, 1.0)]
>>> out = apply_peierls_phase(positions, bonds, flux_quanta_per_plaquette=0.25, area_per_plaquette=1.0)
>>> bool(abs(out[0][2] - 1.0) < 1e-12)
True

physicskit.condensed.bdg_bcs_hamiltonian(k, mu=0.0, t=1.0, delta=0.5)[source]#

Mean-field Bogoliubov-de Gennes Hamiltonian for a 1D s-wave BCS superconductor.

Parameters:

k (float) – Reduced crystal momentum, periodic on \([0, 2\pi)\).
mu (float) – Chemical potential.
t (float) – Nearest-neighbor hopping amplitude (sets the normal-state band \(\xi(k) = -2t\cos k - \mu\)).
delta (complex) – s-wave (momentum-independent) pairing amplitude.

Return type:

Returns:

ndarray, shape (2, 2) – BdG Hamiltonian \(H(k) = \xi(k)\tau_z + \mathrm{Re}(\Delta)\tau_x - \mathrm{Im}(\Delta)\tau_y\) in the Nambu basis \((c_k, c_{-k}^\dagger)\).

See also

bdg_spectrum: Quasiparticle energies over a grid of k, exposing the gap.

Examples

>>> import numpy as np
>>> H = bdg_bcs_hamiltonian(k=np.pi / 2, mu=0.0, t=1.0, delta=0.5)
>>> np.round(np.linalg.eigvalsh(H), 8)
array([-0.5,  0.5])

physicskit.condensed.bdg_spectrum(mu=0.0, t=1.0, delta=0.5, n_k=200)[source]#

Quasiparticle (positive-energy) BdG spectrum over the 1D Brillouin zone.

Parameters:

mu (float) – Chemical potential.
t (float) – Nearest-neighbor hopping amplitude.
delta (complex) – s-wave pairing amplitude.
n_k (int) – Number of momentum points.

Return type:

Returns:

k_grid (ndarray, shape (n_k,)) – Momentum grid over \([0, 2\pi)\).
energies (ndarray, shape (n_k,)) – Positive quasiparticle branch \(E(k) = \sqrt{\xi(k)^2 + |\Delta|^2}\).

Examples

The minimum quasiparticle energy equals the pairing gap \(|\Delta|\):

>>> k, E = bdg_spectrum(mu=0.0, t=1.0, delta=0.5, n_k=400)
>>> round(float(E.min()), 3)
0.5

physicskit.condensed.bethe_ansatz_xxx_ground_energy(N, J=1.0)[source]#

Ground-state energy of the periodic spin-1/2 Heisenberg antiferromagnet from the Bethe ansatz.

With \(M = N/2\) magnons of rapidities \(\lambda_j\), the logarithmic Bethe equations for \(H = J\sum_i \mathbf{S}_i\cdot \mathbf{S}_{i+1}\) read

\[2N \arctan(2\lambda_j) = 2\pi I_j + 2\sum_{k \ne j} \arctan(\lambda_j - \lambda_k),\]

with the ground state given by the consecutive quantum numbers \(I_j = -\tfrac{M-1}{2}, \dots, \tfrac{M-1}{2}\), and

\[E_0 = \frac{JN}{4} - \frac{J}{2}\sum_{j=1}^{M} \frac{1}{\lambda_j^2 + 1/4}\]

(Bethe 1931; Karbach and Muller, Computers in Physics 11, 36 (1997), Eqs. 20-24). As \(N \to \infty\), \(E_0/N \to J(1/4 - \ln 2)\) (Hulthen 1938).

Parameters:

N (int) – Even number of sites, \(N \ge 4\).
J (float) – Exchange coupling.

Return type:

Returns:

float

Examples

>>> round(bethe_ansatz_xxx_ground_energy(4), 12)
-2.0

physicskit.condensed.bhz_hamiltonian(kx, ky, A=1.0, B=1.0, M=1.0, D=0.0)[source]#

Bloch Hamiltonian of the Bernevig-Hughes-Zhang (BHZ) model for HgTe quantum wells.

A minimal 4-band \(\mathbb{Z}_2\) topological insulator model, block-diagonal in a time-reversed pair of 2x2 Dirac-like blocks. Basis order is (E up, H up, E down, H down).

Parameters:

kx (float) – Reduced crystal momenta, each periodic on \([0, 2\pi)\).
ky (float) – Reduced crystal momenta, each periodic on \([0, 2\pi)\).
A (float) – Model parameters controlling the Dirac velocity and quadratic band curvature.
B (float) – Model parameters controlling the Dirac velocity and quadratic band curvature.
M (float) – Band inversion mass. With \(d_z = M - 2B(2-\cos k_x-\cos k_y)\) the model is topological (band-inverted, spin Chern number \(\pm1\)) for \(0 < M/B < 8\) and trivial otherwise (the gap closes at \(\Gamma\), \(X\)/\(Y\), and \(M\) for \(M/B = 0, 4, 8\)).
D (float) – Particle-hole asymmetry parameter, entering as \(\epsilon(k) = -2D(2-\cos k_x-\cos k_y)\), the lattice regularization of the continuum \(-Dk^2\) (Qi, Hughes & Zhang 2008), matching \(B\)’s.

Return type:

Returns:

ndarray, shape (4, 4)

Examples

>>> import numpy as np
>>> H = bhz_hamiltonian(0.1, -0.2, M=1.0, B=1.0)
>>> np.allclose(H, H.conj().T)
True
>>> np.allclose(H[:2, :2], H[2:, 2:].conj())
True

physicskit.condensed.bhz_ribbon_hamiltonian(kx, n_cells, A=1.0, B=1.0, M=1.0, D=0.0)[source]#

Real-space BHZ ribbon: periodic along x, open (finite) along y.

Truncating the y direction exposes the pair of helical, spin-locked edge states that make the BHZ model a quantum spin Hall insulator – the effect Konig et al. (2007) measured directly in HgTe/CdTe quantum wells, by two-terminal conductance quantized at \(2e^2/h\).

Parameters:

kx (float) – Reduced crystal momentum along the periodic (x) direction.
n_cells (int) – Number of unit cells stacked along the open (y) direction.
A (float) – Model parameters controlling the Dirac velocity and quadratic band curvature (see bhz_hamiltonian()).
B (float) – Model parameters controlling the Dirac velocity and quadratic band curvature (see bhz_hamiltonian()).
M (float) – Band inversion mass. Topological (edge-state-carrying) for \(0 < M/B < 8\).
D (float) – Particle-hole asymmetry parameter.

Return type:

Returns:

ndarray, shape (4 * n_cells, 4 * n_cells) – Hermitian, open-boundary ribbon Hamiltonian, block-diagonal in the (E up, H up) / (E down, H down) time-reversed sectors.

Examples

In the topological regime (\(0 < M/B < 8\)), the ribbon has a pair of near-zero-energy states crossing at \(k_x = 0\) – the helical edge modes – absent in the trivial regime (e.g. \(M/B < 0\)):

>>> import numpy as np
>>> spectrum = np.linalg.eigvalsh(bhz_ribbon_hamiltonian(kx=0.0, n_cells=40, M=1.0, B=1.0))
>>> bool(np.any(np.abs(spectrum) < 1e-6))
True
>>> spectrum_trivial = np.linalg.eigvalsh(bhz_ribbon_hamiltonian(kx=0.0, n_cells=40, M=-1.0, B=1.0))
>>> bool(np.any(np.abs(spectrum_trivial) < 1e-6))
False

physicskit.condensed.boltzmann_transport(energies, velocities, cell_volume, mu, T, tau=1.0, q=-1.0, spin_degeneracy=2.0)[source]#

Boltzmann transport coefficients in the constant relaxation-time approximation.

For a band \(\varepsilon_{\mathbf k}\) with velocities \(\mathbf v_{\mathbf k} = \nabla_{\mathbf k}\varepsilon\), define the moments

\[\mathcal K_n = \tau \int \frac{g_s\,d^dk}{(2\pi)^d} \left(-\frac{\partial f}{\partial\varepsilon}\right) \mathbf v_{\mathbf k}\mathbf v_{\mathbf k}^T (\varepsilon_{\mathbf k} - \mu)^n.\]

Then (Ashcroft and Mermin, Eqs. 13.37-13.47, for general charge \(q\))

\[\sigma = q^2\mathcal K_0, \qquad S = \frac{1}{qT}\,\mathcal K_0^{-1}\mathcal K_1, \qquad \kappa = \frac{1}{T}\left(\mathcal K_2 - \mathcal K_1\mathcal K_0^{-1}\mathcal K_1\right).\]

For a degenerate metal the Sommerfeld expansion gives the Wiedemann-Franz law \(\kappa/(\sigma T) = \pi^2/(3q^2)\) and Mott’s formula \(S = \frac{\pi^2 T}{3q}\,\frac{d\ln\sigma(\varepsilon)}{d\varepsilon}\big|_\mu\); for a parabolic band \(\sigma = nq^2\tau/m\), Drude’s result with \(n\) the actual carrier density.

Parameters:

energies (ndarray) – Band energies on a uniform \(k\)-grid, shape (M,) (flatten multi-band or multi-dimensional grids).
velocities (ndarray) – Group velocities at the same points, shape (M, d).
cell_volume (float) – \(k\)-space volume per grid point, \(\Delta^d k\).
mu (float) – Chemical potential.
T (float) – Temperature, \(T > 0\).
tau (float) – Relaxation time.
q (float) – Carrier charge.
spin_degeneracy (float) – \(g_s\).

Return type:

BoltzmannTransport

Returns:

BoltzmannTransport

Examples

A 1D parabolic band at low temperature reproduces Drude’s \(\sigma = nq^2\tau/m\):

>>> import numpy as np
>>> k = np.linspace(-4, 4, 40001)
>>> res = boltzmann_transport(k**2 / 2, k[:, None], k[1] - k[0], mu=1.0, T=0.01)
>>> bool(np.isclose(res.sigma[0, 0], res.density, rtol=1e-3))
True

physicskit.condensed.build_ribbon(hamiltonian, open_direction, n_cells)[source]#

Build a ribbon/slab: periodic in-plane, open (finite) along open_direction.

Truncates the periodic boundary condition along one primitive-lattice direction, producing a quasi-1D (2D lattice) or quasi-2D (3D lattice) strip Hamiltonian as a function of the remaining reduced momenta. This exposes edge/surface states localized at the two open boundaries.

Parameters:

hamiltonian (Hamiltonian) – A Hamiltonian built on a 2D (or higher) Lattice.
open_direction (int) – Index of the primitive-lattice direction to truncate.
n_cells (int) – Number of unit cells stacked along open_direction.

Returns:

callable – A function H_ribbon(k_parallel) returning the (n_cells * n_orbitals, n_cells * n_orbitals) ribbon Hamiltonian, where k_parallel is a reduced-momentum vector with one fewer component than the bulk lattice (the remaining periodic directions).

Examples

SSH chain cut open into a finite 20-site wire; the topological phase (v < w) hosts a mid-gap zero mode:

>>> import numpy as np
>>> from physicskit.condensed.models import ssh_lattice_hamiltonian
>>> H = ssh_lattice_hamiltonian(v=0.5, w=1.0)
>>> H_wire = build_ribbon(H, open_direction=0, n_cells=30)
>>> spectrum = np.linalg.eigvalsh(H_wire(np.array([])))
>>> bool(np.any(np.abs(spectrum) < 1e-6))
True

physicskit.condensed.calabrese_cardy_entropy(ell, N, c, const=0.0, periodic=True)[source]#

Calabrese-Cardy entanglement entropy of an interval in a critical (CFT) chain.

\[S(\ell) = \frac{c}{3\eta}\,\ln\!\left[\frac{\eta N}{\pi} \sin\frac{\pi\ell}{N}\right] + c',\]

with \(\eta = 1\) for a periodic chain and \(\eta = 2\) for an open chain cut at distance \(\ell\) from an edge (Calabrese and Cardy, J. Stat. Mech. P06002 (2004), Eqs. 3 and 14). A gapped chain instead obeys an area law: \(S(\ell)\) saturates.

Parameters:

ell (ndarray) – Subsystem sizes.
N (int) – Chain length.
c (float) – Central charge (\(1/2\) for the critical Ising chain, \(1\) for the Heisenberg chain).
const (float) – Non-universal constant \(c'\).
periodic (bool) – Boundary conditions.

Return type:

Returns:

numpy.ndarray

Examples

>>> import numpy as np
>>> S = calabrese_cardy_entropy(np.array([4, 8]), 16, c=3.0)
>>> bool(np.isclose(S[1], np.log(16 / np.pi)))
True

physicskit.condensed.compute_berry_curvature(hamiltonian_func, grid_size=50, band_index=0)[source]#

Compute the discretized Berry curvature of one band over the Brillouin zone.

Uses the Fukui-Hatsugai-Suzuki (FHS) link-variable formula, which is gauge invariant plaquette by plaquette and needs no smooth choice of eigenvector phase.

Parameters:

hamiltonian_func (callable) – A function H(k1, k2) returning the (N, N) complex Bloch Hamiltonian at reduced crystal momentum (k1, k2), each periodic on \([0, 2\pi)\) (see physicskit.condensed.tight_binding).
grid_size (int) – Number of plaquettes along each of \(k_1, k_2\).
band_index (int) – Band index (0 = lowest energy) to compute the curvature for.

Return type:

Returns:

ndarray, shape (grid_size, grid_size) – Berry flux through each plaquette, in \((-\pi, \pi]\). Summing and dividing by \(2\pi\) gives the band’s Chern number.

See also

compute_chern_number: Integrates this curvature over the full zone for every band.

Examples

>>> import numpy as np
>>> from physicskit.condensed.models import haldane_model
>>> H_func = lambda k1, k2: haldane_model(k1, k2, t=1.0, t2=0.2, phi=np.pi / 2, M=0.0)
>>> F = compute_berry_curvature(H_func, grid_size=30, band_index=0)
>>> F.shape
(30, 30)
>>> round(float(F.sum() / (2 * np.pi)))
1

physicskit.condensed.compute_chern_number(hamiltonian_func, grid_size=50)[source]#

Compute the Chern number of every band using the Fukui-Hatsugai-Suzuki method.

Parameters:

hamiltonian_func (callable) – A function H(k1, k2) returning an (N, N) complex Bloch Hamiltonian, using the reduced-momentum convention of physicskit.condensed.tight_binding.
grid_size (int) – Discretization resolution for \(k_1\) and \(k_2\) across \([0, 2\pi)\).

Return type:

list

Returns:

list of int – Chern integer for each energy band, ordered from lowest to highest energy.

Examples

>>> import numpy as np
>>> from physicskit.condensed.models import haldane_model
>>> H_func = lambda k1, k2: haldane_model(k1, k2, t=1.0, t2=0.2, phi=np.pi / 2, M=0.0)
>>> chern_numbers = compute_chern_number(H_func, grid_size=30)
>>> print(chern_numbers)
[1, -1]

physicskit.condensed.cyclotron_frequency(B, m=1.0, e=1.0)[source]#

Cyclotron frequency \(\omega_c = eB/m\).

Parameters:

B (float) – Magnetic field strength.
m (float) – Particle mass.
e (float) – Particle charge magnitude.

Return type:

Returns:

float

Examples

>>> cyclotron_frequency(B=2.0, m=0.5)
4.0

physicskit.condensed.debye_heat_capacity(T, theta_D)[source]#

Debye heat capacity per atom of a 3D solid, in units of \(k_B\).

\[\frac{C}{Nk_B} = 9\left(\frac{T}{\Theta_D}\right)^3 \int_0^{\Theta_D/T}\frac{x^4 e^x}{(e^x - 1)^2}\,dx\]

(Debye 1912; Kittel, Eq. 5.35), interpolating between the Dulong-Petit value 3 for \(T \gg \Theta_D\) and the Debye \(T^3\) law \(\frac{12\pi^4}{5}(T/\Theta_D)^3\) for \(T \ll \Theta_D\).

Parameters:

T (ndarray) – Temperature (same units as \(\Theta_D\)).
theta_D (float) – Debye temperature.

Return type:

Returns:

numpy.ndarray

Examples

>>> import numpy as np
>>> T = 0.01
>>> bool(np.isclose(debye_heat_capacity(T, 1.0), 12 * np.pi**4 / 5 * T**3, rtol=1e-6))
True

physicskit.condensed.debye_temperature(sound_speed, number_density)[source]#

Debye temperature \(\Theta_D = c\,(6\pi^2 n)^{1/3}\) (units \(\hbar = k_B = 1\)).

Cutting the linear dispersion \(\omega = ck\) off at the Debye wavenumber \(k_D = (6\pi^2 n)^{1/3}\), which holds exactly \(N\) modes per polarization (Kittel, Eq. 5.30), with \(c\) the (suitably averaged) sound speed.

Parameters:

sound_speed (float) – Sound speed \(c\).
number_density (float) – Atoms per unit volume \(n\).

Return type:

Returns:

float

Examples

>>> import numpy as np
>>> bool(np.isclose(debye_temperature(1.0, 1 / (6 * np.pi**2)), 1.0))
True

physicskit.condensed.diatomic_chain_dispersion(k, K=1.0, m1=1.0, m2=2.0, a=1.0)[source]#

Acoustic and optical branches of a diatomic harmonic chain.

Alternating masses \(m_1, m_2\) joined by identical springs \(K\), with unit-cell length \(a\):

\[\omega_\pm^2 = K\left(\frac{1}{m_1} + \frac{1}{m_2}\right) \pm K\sqrt{\left(\frac{1}{m_1} + \frac{1}{m_2}\right)^2 - \frac{4\sin^2(ka/2)}{m_1 m_2}}\]

(Kittel, Introduction to Solid State Physics, 8th ed., Eq. 4.20). The acoustic branch has sound speed \(c = a\sqrt{K/(2(m_1+m_2))}\); the optical branch starts at \(\sqrt{2K(1/m_1+1/m_2)}\), and a gap opens at the zone boundary between \(\sqrt{2K/m_{\max}}\) and \(\sqrt{2K/m_{\min}}\).

Parameters:

k (ndarray) – Wavenumber.
K (float) – Spring constant.
m1 (float) – The two atomic masses.
m2 (float) – The two atomic masses.
a (float) – Unit-cell length.

Return type:

tuple[ndarray, ndarray]

Returns:

acoustic, optical (numpy.ndarray)

Examples

>>> import numpy as np
>>> ac, op = diatomic_chain_dispersion(np.pi, K=1.0, m1=1.0, m2=2.0)
>>> bool(np.isclose(ac, 1.0) and np.isclose(op, np.sqrt(2.0)))
True

physicskit.condensed.drude_ac_conductivity(omega, n, tau, q=-1.0, m=1.0)[source]#

Drude AC conductivity \(\sigma(\omega) = \sigma_0/(1 - i\omega\tau)\).

For fields \(\propto e^{-i\omega t}\) (Ashcroft and Mermin, Eq. 1.29). The real part is a Lorentzian of width \(1/\tau\), whose integral obeys the f-sum rule \(\int_0^\infty \mathrm{Re}\,\sigma\,d\omega = \pi nq^2/(2m)\).

Parameters:

omega (ndarray) – Angular frequency.
n (float) – As in drude_conductivity().
tau (float) – As in drude_conductivity().
q (float) – As in drude_conductivity().
m (float) – As in drude_conductivity().

Return type:

Returns:

numpy.ndarray – Complex conductivity.

Examples

>>> complex(drude_ac_conductivity(1.0, n=1.0, tau=1.0))
(0.5+0.5j)

physicskit.condensed.drude_conductivity(n, tau, q=-1.0, m=1.0)[source]#

Drude DC conductivity \(\sigma_0 = nq^2\tau/m\).

Carriers of charge \(q\) and mass \(m\), accelerated by the field and randomized by collisions every \(\tau\) on average (Drude 1900; Ashcroft and Mermin, Eq. 1.6).

Parameters:

n (float) – Carrier density.
tau (float) – Relaxation (mean free) time.
q (float) – Carrier charge.
m (float) – Carrier mass.

Return type:

Returns:

float

Examples

>>> drude_conductivity(n=2.0, tau=3.0)
6.0

physicskit.condensed.drude_conductivity_tensor(B, n, tau, q=-1.0, m=1.0)[source]#

In-plane Drude conductivity tensor in a perpendicular magnetic field \(B\hat z\).

Solving \(m\dot{\mathbf v} = q(\mathbf E + \mathbf v\times\mathbf B) - m\mathbf v/\tau\) in steady state gives

\[\begin{split}\sigma = \frac{\sigma_0}{1 + (\omega_c\tau)^2} \begin{pmatrix} 1 & \omega_c\tau \\ -\omega_c\tau & 1 \end{pmatrix}, \qquad \omega_c = \frac{qB}{m},\end{split}\]

whose inverse, the resistivity, has \(\rho_{xx} = 1/\sigma_0\) independent of \(B\) (no magnetoresistance) and \(\rho_{yx} = B/(nq)\) (Ashcroft and Mermin, Eqs. 1.40-1.44).

Parameters:

B (float) – Magnetic field along \(z\).
n (float) – As in drude_conductivity().
tau (float) – As in drude_conductivity().
q (float) – As in drude_conductivity().
m (float) – As in drude_conductivity().

Return type:

Returns:

numpy.ndarray – Shape (2, 2).

Examples

>>> import numpy as np
>>> rho = np.linalg.inv(drude_conductivity_tensor(B=2.0, n=0.5, tau=1.0))
>>> bool(np.isclose(rho[1, 0], 2.0 / (0.5 * -1.0)))
True

physicskit.condensed.embed_state(psi, basis, N)[source]#

Embed a sector-basis state vector into the full \(2^N\)-dimensional space.

Parameters:

psi (ndarray) – Amplitudes on basis.
basis (ndarray) – Sector basis states.
N (int) – Number of sites.

Return type:

Returns:

numpy.ndarray – Shape (2**N,).

Examples

>>> import numpy as np
>>> embed_state(np.array([1.0, 0.0]), np.array([1, 2]), 2)
array([0., 1., 0., 0.])

physicskit.condensed.energy_gap(H)[source]#

Gap \(E_1 - E_0\) between the two lowest eigenvalues of H.

Parameters:: H (spmatrix | ndarray) – Hamiltonian (typically restricted to one symmetry sector).
Return type:: float
Returns:: float

Examples

The critical TFIM’s gap in the even sector is the two-quasiparticle energy \(\varepsilon_{\pi/N} + \varepsilon_{-\pi/N}\):

>>> import numpy as np
>>> H, _ = tfim_hamiltonian(8, h=1.0, parity=+1)
>>> bool(np.isclose(energy_gap(H), 8 * np.sin(np.pi / 16)))
True

physicskit.condensed.entanglement_entropy(psi, N, n_A, basis=None)[source]#

Von Neumann entanglement entropy of the first n_A sites.

\[S_A = -\operatorname{tr}\rho_A \ln \rho_A = -\sum_\alpha s_\alpha^2 \ln s_\alpha^2,\]

with \(s_\alpha\) the Schmidt coefficients – the singular values of \(\psi\) reshaped to a \(2^{n_A} \times 2^{N - n_A}\) matrix (Nielsen and Chuang, Sec. 2.5).

Parameters:

psi (ndarray) – Normalized state, on the full space or on basis.
N (int) – Number of sites.
n_A (int) – Number of leading sites in subsystem \(A\).
basis (ndarray | None) – Sector basis of psi, if it is not a full-space vector.

Return type:

Returns:

float – \(S_A\) in nats.

Examples

A singlet carries \(\ln 2\) of entanglement:

>>> import numpy as np
>>> singlet = np.array([0, 1, -1, 0]) / np.sqrt(2)
>>> bool(np.isclose(entanglement_entropy(singlet, 2, 1), np.log(2)))
True

physicskit.condensed.entanglement_profile(psi, N, basis=None)[source]#

Entanglement entropy \(S(\ell)\) for every cut \(\ell = 1, \dots, N-1\).

Parameters:

psi (ndarray) – Normalized state.
N (int) – Number of sites.
basis (ndarray | None) – Sector basis of psi.

Return type:

Returns:

numpy.ndarray – Shape (N - 1,).

Examples

A product state is unentangled across every cut:

>>> import numpy as np
>>> psi = np.zeros(8); psi[0] = 1.0
>>> entanglement_profile(psi, 3)
array([0., 0.])

physicskit.condensed.fermi_window(energy, mu, T)[source]#

The Fermi window \(-\partial f/\partial\varepsilon = \dfrac{1}{4T\cosh^2[(\varepsilon - \mu)/2T]}\).

Parameters:

energy (ndarray) – Energies.
mu (float) – Chemical potential.
T (float) – Temperature, \(T > 0\).

Return type:

Returns:

numpy.ndarray – Integrates to 1 over energy.

Examples

>>> float(fermi_window(0.0, mu=0.0, T=0.25))
1.0

physicskit.condensed.filling_factor(density, B, hbar=1.0, e=1.0)[source]#

Landau-level filling factor \(\nu = n_e / n_B = n_e h/(eB)\).

The number of filled Landau levels (generally non-integer) for a 2D electron density density. Integer \(\nu\) is the condition for an incompressible quantum Hall plateau in physicskit.condensed.topology.

Parameters:

density (float) – 2D electron number density (particles per unit area).
B (float) – Magnetic field strength.
hbar (float) – Reduced Planck constant.
e (float) – Particle charge magnitude.

Return type:

Returns:

float

Examples

>>> round(filling_factor(density=2.0, B=1.0), 6)
12.566371

physicskit.condensed.fit_central_charge(profile, N, periodic=True, trim=1)[source]#

Least-squares fit of calabrese_cardy_entropy() to an entanglement profile.

Parameters:

profile (ndarray) – \(S(\ell)\) for \(\ell = 1, \dots, N-1\) (from entanglement_profile()).
N (int) – Chain length.
periodic (bool) – Boundary conditions.
trim (int) – Number of cuts dropped at each end, where lattice corrections are largest.

Return type:

tuple[float, float]

Returns:

c (float) – Fitted central charge.
const (float) – Fitted non-universal constant.

Examples

>>> import numpy as np
>>> ell = np.arange(1, 20)
>>> c, _ = fit_central_charge(calabrese_cardy_entropy(ell, 20, c=0.5, const=0.3), 20)
>>> round(c, 10)
0.5

physicskit.condensed.ginzburg_landau_parameter(coherence_length, penetration_depth)[source]#

Ginzburg-Landau parameter \(\kappa = \lambda/\xi\).

The single dimensionless number that decides a superconductor’s response to a magnetic field: \(\kappa < 1/\sqrt2\) is Type I (the normal-superconducting interface has positive surface energy, and the field is excluded entirely below \(H_c\)); \(\kappa > 1/\sqrt2\) is Type II (negative surface energy favors flux penetrating as an Abrikosov vortex lattice between \(H_{c1}\) and \(H_{c2}\)).

Parameters:

coherence_length (float) – \(\xi\), from gl_coherence_length().
penetration_depth (float) – \(\lambda\), from gl_penetration_depth().

Return type:

Returns:

float

Examples

>>> round(ginzburg_landau_parameter(coherence_length=1.0, penetration_depth=1.0), 4)
1.0

physicskit.condensed.gl_coherence_length(a, hbar=1.0, m=1.0)[source]#

Ginzburg-Landau coherence length \(\xi = \hbar/\sqrt{2m|a|}\).

The length scale over which the order parameter heals back to its bulk value after being suppressed at a boundary or a vortex core.

Parameters:

a (float) – Quadratic coefficient (only its magnitude matters).
hbar (float) – Reduced Planck constant.
m (float) – Effective mass of the condensate.

Return type:

Returns:

float

Examples

>>> gl_coherence_length(a=-0.5)
1.0

physicskit.condensed.gl_equilibrium_order_parameter(a, b)[source]#

Equilibrium (uniform, field-free) order parameter magnitude \(|\psi_0|=\sqrt{-a/b}\).

Parameters:

a (float) – Quadratic coefficient.
b (float) – Quartic coefficient, b > 0.

Return type:

Returns:

float – \(\sqrt{-a/b}\) if a < 0 (ordered phase), else 0.0 (disordered phase, only the normal state minimizes \(f\)).

Examples

>>> gl_equilibrium_order_parameter(a=-2.0, b=2.0)
1.0
>>> gl_equilibrium_order_parameter(a=1.0, b=2.0)
0.0

physicskit.condensed.gl_free_energy_density(psi, a, b, grad_psi=0.0, hbar=1.0, m=1.0)[source]#

Ginzburg-Landau free energy density \(f = a|\psi|^2 + \tfrac{b}{2}|\psi|^4 + \tfrac{\hbar^2}{2m}|\nabla\psi|^2\).

Parameters:

psi (complex or array_like) – Order parameter value(s).
a (float) – Quadratic coefficient. Changes sign at the transition (a < 0 in the ordered phase, a > 0 in the disordered phase).
b (float) – Quartic coefficient, b > 0 for stability.
grad_psi (complex or array_like, default=0.0) – Gradient \(\nabla\psi\), same shape as psi.
hbar (float) – Reduced Planck constant.
m (float) – Effective mass of the condensate.

Returns:

float or ndarray

Examples

>>> gl_free_energy_density(psi=0.0, a=-1.0, b=1.0)
0.0
>>> psi0 = gl_equilibrium_order_parameter(a=-1.0, b=1.0)
>>> round(gl_free_energy_density(psi0, a=-1.0, b=1.0), 6)
-0.5

physicskit.condensed.gl_order_parameter_profile(x, xi)[source]#

Order parameter healing profile \(\psi(x)/\psi_0 = \tanh(x/(\sqrt2\,\xi))\).

The exact solution of the dimensionless Ginzburg-Landau equation \(\xi^2\psi'' = \psi^3-\psi\) for a condensate pinned to zero at a boundary (x = 0, e.g. a normal-superconducting interface) and recovering its bulk value far from it – the direct, textbook illustration of the coherence length as a healing length.

Parameters:

x (array_like) – Distance from the boundary.
xi (float) – Coherence length, from gl_coherence_length().

Return type:

Returns:

ndarray – \(\psi(x)/\psi_0\), in \([0, 1)\) for \(x \geq 0\).

Examples

>>> import numpy as np
>>> round(float(gl_order_parameter_profile(0.0, xi=1.0)), 8)
0.0
>>> bool(gl_order_parameter_profile(np.array([10.0]), xi=1.0)[0] > 0.999)
True

physicskit.condensed.gl_penetration_depth(psi0, e=1.0, m=1.0)[source]#

London penetration depth \(\lambda = \sqrt{m/(e^2|\psi_0|^2)}\).

The length scale over which an external magnetic field is screened from a superconductor’s interior by the supercurrent it induces in the condensate of density \(|\psi_0|^2\).

Parameters:

psi0 (float) – Equilibrium order parameter magnitude (see gl_equilibrium_order_parameter()).
e (float) – Cooper-pair charge magnitude.
m (float) – Effective mass of the condensate.

Return type:

Returns:

float

Examples

>>> gl_penetration_depth(psi0=1.0)
1.0

physicskit.condensed.graphene_hamiltonian(k1, k2, t=1.0)[source]#

Bloch Hamiltonian of nearest-neighbor graphene (honeycomb lattice).

Parameters:

k1 (float) – Reduced crystal momenta conjugate to the two honeycomb primitive vectors, each periodic on \([0, 2\pi)\).
k2 (float) – Reduced crystal momenta conjugate to the two honeycomb primitive vectors, each periodic on \([0, 2\pi)\).
t (float) – Nearest-neighbor hopping amplitude.

Return type:

Returns:

ndarray, shape (2, 2) – Bloch Hamiltonian on the A/B sublattice basis.

Notes

The Dirac points sit at \((k_1, k_2) = (2\pi/3, 4\pi/3)\) and its time-reversed partner \((4\pi/3, 2\pi/3)\), where the gap closes and the dispersion is linear (massless Dirac cone), as verified in the Examples.

Examples

>>> import numpy as np
>>> K = np.array([2 * np.pi / 3, 4 * np.pi / 3])
>>> np.round(np.linalg.eigvalsh(graphene_hamiltonian(*K)), 8)
array([-0.,  0.])
>>> eigs = np.linalg.eigvalsh(graphene_hamiltonian(*(K + [1e-4, 0])))
>>> round(float(eigs[1] / 1e-4), 4)
1.0

physicskit.condensed.graphene_lattice_hamiltonian(t=1.0)[source]#

Nearest-neighbor honeycomb (graphene) tight-binding Hamiltonian.

Parameters:: t (float) – Nearest-neighbor hopping amplitude.
Return type:: Hamiltonian
Returns:: Hamiltonian

See also

graphene_hamiltonian: Closed-form 2x2 Bloch Hamiltonian equivalent to this builder.

physicskit.condensed.haldane_lattice_hamiltonian(t=1.0, t2=0.2, phi=1.5707963267948966, M=0.0)[source]#

Haldane model (quantum anomalous Hall) as a Hamiltonian.

Parameters:

t (float) – Nearest-neighbor hopping amplitude.
t2 (float) – Next-nearest-neighbor hopping magnitude.
phi (float) – Next-nearest-neighbor hopping phase, breaking time-reversal symmetry.
M (float) – Sublattice (Semenoff) mass, breaking inversion symmetry.

Return type:

Returns:

Hamiltonian

See also

haldane_model: Closed-form 2x2 Bloch Hamiltonian equivalent to this builder.

physicskit.condensed.haldane_model(kx, ky, t=1.0, t2=0.2, phi=1.5707963267948966, M=0.0)[source]#

Bloch Hamiltonian of the Haldane model on the honeycomb lattice.

The Haldane model realizes the quantum anomalous Hall effect: a Chern insulator with zero net magnetic flux per unit cell, arising from a complex next-nearest-neighbor hopping \(t_2 e^{i\phi}\) that breaks time-reversal symmetry while preserving the lattice translational symmetry.

Parameters:

kx (float) – Reduced crystal momenta (see physicskit.condensed.tight_binding), each periodic on \([0, 2\pi)\).
ky (float) – Reduced crystal momenta (see physicskit.condensed.tight_binding), each periodic on \([0, 2\pi)\).
t (float) – Nearest-neighbor hopping amplitude.
t2 (float) – Next-nearest-neighbor hopping magnitude.
phi (float) – Next-nearest-neighbor hopping phase.
M (float) – Sublattice mass. The model is a Chern insulator (\(C = \mathrm{sgn}(\sin\phi)\)) for \(|M| < 3\sqrt{3}\,t_2|\sin\phi|\), and a trivial insulator otherwise.

Return type:

Returns:

ndarray, shape (2, 2)

See also

physicskit.condensed.topology.compute_chern_number: Computes the Chern number of this model.

Examples

>>> import numpy as np
>>> from physicskit.condensed.topology import compute_chern_number
>>> H_func = lambda k1, k2: haldane_model(k1, k2, t=1.0, t2=0.2, phi=np.pi / 2, M=0.0)
>>> compute_chern_number(H_func, grid_size=30)
[1, -1]
>>> H_trivial = lambda k1, k2: haldane_model(k1, k2, t=1.0, t2=0.2, phi=np.pi / 2, M=2.0)
>>> compute_chern_number(H_trivial, grid_size=30)
[0, 0]

physicskit.condensed.hall_coefficient(n, q=-1.0)[source]#

Drude Hall coefficient \(R_H = E_y/(j_x B) = 1/(nq)\).

Negative for electrons; its measured sign revealed hole conduction in metals such as Be and Al that the free-electron model could not explain (Ashcroft and Mermin, Table 1.4).

Parameters:

n (float) – Carrier density.
q (float) – Carrier charge.

Return type:

Returns:

float

Examples

>>> hall_coefficient(0.25)
-4.0

physicskit.condensed.harper_hofstadter_hamiltonian(k1, k2, p, q, t=1.0)[source]#

Bloch Hamiltonian of the square-lattice Harper-Hofstadter model at flux p/q.

A charged particle on a square lattice threaded by a uniform magnetic flux \(p/q\) (in units of the flux quantum) per plaquette, in Landau gauge \(\mathbf{A} = (0, Bx)\). Translational symmetry along x survives only in steps of q lattice constants, so the q inequivalent sublattice sites m = 0, ..., q-1 within one magnetic unit cell become an internal band index, coupled by

\[H_{mm}(k_1, k_2) = -2t\cos\!\left(k_2 + 2\pi \frac{p}{q} m\right), \qquad H_{m, m+1} = -t,\]

plus the boundary hopping \(H_{q-1, 0} = -t\,e^{ik_1}\) that closes the magnetic unit cell, with \(k_1\) the reduced momentum conjugate to translation by one magnetic cell (period \(q\) sites) and \(k_2\) the ordinary reduced momentum along the unbroken direction. This is the lattice (Bloch) route to the same physics as physicskit.condensed.landau_levels, and, fed one band at a time into compute_chern_number(), the standard numerical verification of the TKNN integer quantum Hall formula: summing the Chern numbers of the lowest \(r\) (non-touching) bands gives the exactly quantized Hall conductance \(\sigma_{xy} = C\,e^2/h\) at the filling between band \(r\) and \(r+1\).

Parameters:

k1 (float) – Reduced crystal momenta, each periodic on \([0, 2\pi)\). k1 is conjugate to the enlarged (q-site) magnetic unit cell.
k2 (float) – Reduced crystal momenta, each periodic on \([0, 2\pi)\). k1 is conjugate to the enlarged (q-site) magnetic unit cell.
p (int) – Flux per plaquette \(p/q\) (in lowest terms), \(0 < p < q\). Odd q avoids the exact band touchings that occur for some even q, where individual-band Chern numbers become ill-defined.
q (int) – Flux per plaquette \(p/q\) (in lowest terms), \(0 < p < q\). Odd q avoids the exact band touchings that occur for some even q, where individual-band Chern numbers become ill-defined.
t (float) – Nearest-neighbor hopping amplitude.

Return type:

Returns:

ndarray, shape (q, q) – Hermitian Bloch Hamiltonian.

See also

physicskit.condensed.topology.compute_chern_number: Chern number of each resulting band.

Examples

At flux 1/3 the three bands carry Chern numbers -1, 2, -1 (summing to zero, as any complete set of bands of a lattice Hamiltonian must), so the Hall conductance is quantized to \(-1\) and \(+1\) (in units of \(e^2/h\)) at the two gaps:

>>> import numpy as np
>>> from physicskit.condensed.topology import compute_chern_number
>>> H_func = lambda k1, k2: harper_hofstadter_hamiltonian(k1, k2, p=1, q=3)
>>> chern_numbers = compute_chern_number(H_func, grid_size=30)
>>> chern_numbers
[-1, 2, -1]
>>> sum(chern_numbers)
0

physicskit.condensed.hubbard_1d_exact_diagonalization(n_sites, n_up, n_dn, t=1.0, U=4.0, pbc=False, return_eigenvectors=False)[source]#

Exact diagonalization of the 1D Fermi-Hubbard model on a small cluster.

\[H = -t \sum_{\langle ij\rangle, \sigma} c_{i\sigma}^\dagger c_{j\sigma} + U \sum_i n_{i\uparrow} n_{i\downarrow}\]

Builds the Fock-space Hamiltonian in the fixed \((n_\uparrow, n_\downarrow)\) sector via full (dense) diagonalization; suitable for small clusters (\(n_{sites} \lesssim 8\)).

Parameters:

n_sites (int) – Number of lattice sites.
n_up (int) – Number of spin-up and spin-down electrons.
n_dn (int) – Number of spin-up and spin-down electrons.
t (float) – Nearest-neighbor hopping amplitude.
U (float) – Onsite Coulomb repulsion.
pbc (bool) – Periodic (True) or open (False) boundary conditions.
return_eigenvectors (bool) – If True, also return the ground-state vector and the occupation- number bases it is expanded in, e.g. for hubbard_spin_correlations().

Return type:

dict

Returns:

dict – {"ground_state_energy": float, "eigenvalues": ndarray, "dimension": int}, plus {"ground_state_vector": ndarray, "states_up": list, "states_dn": list} when return_eigenvectors=True. ground_state_vector[a * len(states_dn) + b] is the amplitude of the product basis state (states_up[a], states_dn[b]), each an integer bitmask with bit i set if site i is occupied.

Notes

In the large-\(U\) limit at half filling, double occupancy is suppressed and the ground state approaches the Mott-insulating, singly-occupied regime characteristic of spin-charge separation in 1D.

Examples

The 2-site Hubbard dimer at half filling has the exact ground-state energy \(E_0 = \tfrac{1}{2}\big(U - \sqrt{U^2 + 16t^2}\big)\):

>>> import numpy as np
>>> result = hubbard_1d_exact_diagonalization(n_sites=2, n_up=1, n_dn=1, t=1.0, U=4.0)
>>> E_exact = 0.5 * (4.0 - np.sqrt(4.0**2 + 16 * 1.0**2))
>>> bool(round(result["ground_state_energy"], 8) == round(E_exact, 8))
True

physicskit.condensed.hubbard_spin_correlations(n_sites, ground_state_vector, states_up, states_dn)[source]#

Equal-time z-spin correlations \(\langle S_i^z S_j^z\rangle\) from a Hubbard ground state.

\(S_i^z = \tfrac{1}{2}(n_{i\uparrow} - n_{i\downarrow})\) is diagonal in the occupation-number basis hubbard_1d_exact_diagonalization() builds its Hamiltonian in, so \(\langle S_i^z S_j^z\rangle\) reduces to a weighted sum over basis-state probabilities, \(\sum_\alpha |\psi_\alpha|^2\, s_i(\alpha)\, s_j(\alpha)\), with no off-diagonal matrix elements to track.

Parameters:

n_sites (int) – Number of lattice sites.
ground_state_vector (ndarray) – Ground-state amplitudes, as returned by hubbard_1d_exact_diagonalization() with return_eigenvectors=True.
states_up (list) – Occupation-number bases for each spin species, as returned alongside ground_state_vector.
states_dn (list) – Occupation-number bases for each spin species, as returned alongside ground_state_vector.

Return type:

Returns:

ndarray, shape (n_sites,) – correlations[r] is \(\langle S_i^z S_j^z\rangle\) averaged over every site pair with \(|i-j|=r\); correlations[0] is the onsite \(\langle (S_i^z)^2\rangle\), suppressed toward 0 by double/empty occupancy and toward its \(1/4\) ceiling by strong onsite repulsion.

Notes

A Mott-insulating antiferromagnet’s short-range Neel order shows up as correlations alternating in sign with r – negative (antialigned) at odd separations, positive (aligned) at even ones – the same short-range magnetic correlations believed to survive doping into the cuprate superconductors’ metallic phase.

Examples

At half filling and strong coupling, nearest-neighbor spins are antialigned:

>>> import numpy as np
>>> result = hubbard_1d_exact_diagonalization(n_sites=6, n_up=3, n_dn=3, t=1.0, U=8.0, pbc=True, return_eigenvectors=True)
>>> corr = hubbard_spin_correlations(6, result["ground_state_vector"], result["states_up"], result["states_dn"])
>>> bool(corr[1] < 0)
True

physicskit.condensed.inverse_participation_ratio(psi)[source]#

Inverse participation ratio \(\text{IPR} = \sum_i|\psi_i|^4 / \left(\sum_i|\psi_i|^2\right)^2\).

A dimensionless measure of how many sites an eigenstate is spread over: \(\text{IPR}\sim 1/N\) for a state extended over all \(N\) sites, and \(\text{IPR} = O(1)\), independent of \(N\), for a state localized on a handful of sites – the diagnostic Anderson used to distinguish the two regimes.

Parameters:: psi (array_like) – Eigenvector components (need not be pre-normalized).
Return type:: float
Returns:: float

Examples

A state spread equally over N sites has IPR = 1/N:

>>> import numpy as np
>>> psi = np.ones(10) / np.sqrt(10)
>>> round(inverse_participation_ratio(psi), 6)
0.1
>>> psi_localized = np.zeros(10); psi_localized[0] = 1.0
>>> inverse_participation_ratio(psi_localized)
1.0

physicskit.condensed.kane_mele_hamiltonian(kx, ky, t=1.0, lambda_so=0.06, lambda_v=0.0, lambda_r=0.0)[source]#

Bloch Hamiltonian of the Kane-Mele quantum spin Hall model.

Two time-reversed Haldane copies (intrinsic spin-orbit coupling with opposite Chern-number-generating phase for each spin), realizing a \(\mathbb{Z}_2\) topological insulator with helical edge states. Basis order is (A up, B up, A down, B down).

Parameters:

kx (float) – Reduced crystal momenta, each periodic on \([0, 2\pi)\).
ky (float) – Reduced crystal momenta, each periodic on \([0, 2\pi)\).
t (float) – Nearest-neighbor hopping amplitude.
lambda_so (float) – Intrinsic spin-orbit coupling strength.
lambda_v (float) – Sublattice (staggered) potential.
lambda_r (float) – Rashba spin-orbit coupling strength. A nonzero value mixes the spin blocks (breaking \(s_z\) conservation) but preserves overall time-reversal symmetry.

Return type:

Returns:

ndarray, shape (4, 4)

Notes

When lambda_r == 0, \(s_z\) is conserved and the Hamiltonian is block-diagonal in spin; each block is a Haldane model with \(t_2 = \lambda_{so}\), \(\phi = \pm\pi/2\). In that case the \(\mathbb{Z}_2\) invariant equals the spin-up Chern number modulo 2 (see physicskit.condensed.topology.z2_invariant()).

Examples

>>> import numpy as np
>>> H = kane_mele_hamiltonian(0.3, 0.7, lambda_so=0.06)
>>> np.allclose(H, H.conj().T)
True
>>> up = H[:2, :2]
>>> down_at_minus_k = kane_mele_hamiltonian(-0.3, -0.7, lambda_so=0.06)[2:, 2:]
>>> np.allclose(down_at_minus_k, up.conj())  # time reversal: H_dn(-k) = H_up(k)*
True

physicskit.condensed.kitaev_chain_bdg_real_space(n_sites, mu=0.0, t=1.0, delta=1.0)[source]#

Real-space open-boundary BdG Hamiltonian of the Kitaev chain.

Parameters:

n_sites (int) – Number of lattice sites in the finite, open chain.
mu (float) – Chemical potential.
t (float) – Nearest-neighbor hopping amplitude.
delta (float) – p-wave pairing amplitude.

Return type:

Returns:

ndarray, shape (2*n_sites, 2*n_sites) – BdG Hamiltonian in the Nambu basis \((c_1, \dots, c_N, c_1^\dagger, \dots, c_N^\dagger)\).

Examples

In the topological phase (\(|\mu| < 2t\)), diagonalizing this matrix yields a pair of near-zero-energy eigenvalues (Majorana end modes), exponentially localized at the two ends of the chain:

>>> import numpy as np
>>> H = kitaev_chain_bdg_real_space(n_sites=60, mu=0.0, t=1.0, delta=1.0)
>>> eigs = np.linalg.eigvalsh(H)
>>> bool(np.abs(eigs[np.argmin(np.abs(eigs))]) < 1e-6)
True

physicskit.condensed.kitaev_chain_hamiltonian(k, mu=0.0, t=1.0, delta=1.0)[source]#

Bogoliubov-de Gennes Bloch Hamiltonian of the Kitaev p-wave superconducting chain.

Parameters:

k (float) – Reduced crystal momentum, periodic on \([0, 2\pi)\).
mu (float) – Chemical potential.
t (float) – Nearest-neighbor hopping amplitude.
delta (float) – p-wave pairing amplitude.

Return type:

Returns:

ndarray, shape (2, 2) – BdG Hamiltonian \(H(k) = \xi(k)\tau_z + \Delta(k)\tau_y\) in the Nambu basis, with \(\xi(k) = -2t\cos k - \mu\) and \(\Delta(k) = 2\Delta\sin k\).

Notes

The chain is a topological superconductor, hosting unpaired Majorana zero modes at the ends of an open chain, for \(|\mu| < 2t\); it is topologically trivial for \(|\mu| > 2t\).

Examples

>>> import numpy as np
>>> H = kitaev_chain_hamiltonian(k=np.pi / 2, mu=0.0, t=1.0, delta=1.0)
>>> np.round(np.linalg.eigvalsh(H), 8)
array([-2.,  2.])

physicskit.condensed.landau_degeneracy(area, B, hbar=1.0, e=1.0)[source]#

Number of degenerate single-particle states per Landau level.

Each Landau level holds \(n_B \cdot \text{area}\) states, where \(n_B = 1/(2\pi\ell_B^2) = eB/(2\pi\hbar) = B/\Phi_0\) is one state per flux quantum \(\Phi_0 = 2\pi\hbar/e\) threading the sample – the origin of the macroscopic degeneracy behind Landau diamagnetism.

Parameters:

area (float) – Real-space sample area.
B (float) – Magnetic field strength.
hbar (float) – Reduced Planck constant.
e (float) – Particle charge magnitude.

Return type:

Returns:

float

Examples

>>> landau_degeneracy(area=2 * np.pi, B=1.0)
1.0

physicskit.condensed.landau_density_of_states(energies, B, m=1.0, e=1.0, hbar=1.0, n_max=20, broadening=0.05)[source]#

Disorder-broadened density of states per unit area.

Replaces each infinitely sharp level \(n_B\,\delta(E - E_n)\) with a Gaussian of width broadening, the standard phenomenological model for the disorder- or finite-lifetime-broadened Landau levels seen in a real 2D electron gas (and needed to have a finite conductivity between the exactly quantized Hall plateaus).

Parameters:

energies (array_like) – Energies at which to evaluate the density of states.
B (float) – Magnetic field strength.
m (float) – Particle mass.
e (float) – Particle charge magnitude.
hbar (float) – Reduced Planck constant.
n_max (int) – Highest Landau index included in the sum.
broadening (float) – Gaussian standard deviation of each broadened level.

Return type:

Returns:

ndarray – Density of states per unit area, same shape as energies.

Examples

>>> dos = landau_density_of_states([0.5], B=1.0, n_max=0, broadening=0.1)
>>> round(float(dos[0]), 4)
0.6349

physicskit.condensed.landau_level_energies(n_max, B, m=1.0, e=1.0, hbar=1.0)[source]#

Landau energy levels \(E_n = \hbar\omega_c(n+\tfrac12)\) for \(n=0,\dots,n_{max}\).

Parameters:

n_max (int) – Highest Landau index to include.
B (float) – Magnetic field strength.
m (float) – Particle mass.
e (float) – Particle charge magnitude.
hbar (float) – Reduced Planck constant.

Return type:

Returns:

ndarray, shape (n_max + 1,) – Equally spaced energies \(E_0 < E_1 < \dots < E_{n_{max}}\), spacing \(\hbar\omega_c\).

See also

landau_degeneracy: Number of degenerate single-particle states per level.

Examples

>>> landau_level_energies(n_max=2, B=1.0)
array([0.5, 1.5, 2.5])

physicskit.condensed.lattice_heat_capacity(frequencies, T)[source]#

Average heat capacity per phonon mode, \(C/(N_\mathrm{modes} k_B)\).

Each harmonic mode of frequency \(\omega\) contributes Einstein’s

\[c(\omega, T) = \left(\frac{\omega}{T}\right)^2 \frac{e^{\omega/T}}{(e^{\omega/T} - 1)^2}\]

(Einstein, Ann. Phys. 327, 180 (1907)), tending to 1 (equipartition) at high \(T\). Summing over a Brillouin-zone sample of a phonon dispersion gives the lattice heat capacity without Debye’s continuum approximation.

Parameters:

frequencies (ndarray) – Mode frequencies (any shape; zero-frequency modes contribute 1).
T (ndarray) – Temperatures.

Return type:

Returns:

numpy.ndarray – Shape of T.

Examples

>>> import numpy as np
>>> round(float(lattice_heat_capacity(np.array([1.0]), 1000.0)), 6)
1.0

physicskit.condensed.laughlin_metropolis_sweep(z, m, step=1.0)[source]#

Perform one Metropolis sweep sampling \(|\Psi_m|^2\), in place.

Proposes a small random displacement of one particle at a time and accepts or rejects it with the exact Metropolis ratio of \(|\Psi_m|^2\) – Laughlin’s “plasma analogy” Boltzmann weight, \(\prod_{i<j}|z_i-z_j|^{2m}\,e^{-\sum_i|z_i|^2/2}\) – with no explicit normalization needed since only the ratio between configurations ever enters.

Parameters:

z (ndarray) – Particle positions \(z_j=x_j+iy_j\). Modified in place.
m (int) – Laughlin exponent (filling \(\nu=1/m\)); odd for fermions.
step (float) – Maximum displacement of a trial move along each of \(x, y\), drawn uniformly from [-step, step].

Return type:

Returns:

ndarray of shape (N,) – The same array passed in.

Notes

Trial moves draw from Numba’s internal RNG, which neither np.random.seed nor a np.random.default_rng generator affects when called from ordinary Python. For a reproducible chain, seed it first with physicskit.statphys.core.monte_carlo.seed_numba_random().

Examples

A two-particle system’s Metropolis chain never proposes an exactly coincident configuration acceptable (\(|\Psi_m|^2 = 0\) there), so repeated sweeps keep the two particles strictly apart:

>>> import numpy as np
>>> from physicskit.statphys.core.monte_carlo import seed_numba_random
>>> seed_numba_random(0)
>>> z = np.array([0.1 + 0.0j, -0.1 + 0.0j])
>>> for _ in range(200):
...     _ = laughlin_metropolis_sweep(z, m=3, step=0.5)
>>> bool(np.abs(z[0] - z[1]) > 1e-3)
True

physicskit.condensed.laughlin_pair_correlation(z_samples, m, r_max, n_bins=60)[source]#

Pair correlation function \(g(r)\) from sampled Laughlin configurations.

Bins every pairwise separation across a stack of equilibrium samples (drawn with laughlin_metropolis_sweep()) and normalizes by the Laughlin state’s exact bulk density \(\rho=1/(2\pi m \ell_B^2)\), so that \(g(r) \to 1\) far from any correlation the wavefunction itself imposes.

Parameters:

z_samples (ndarray) – Stack of particle-position snapshots, ideally decorrelated by several sweeps between samples.
m (int) – Laughlin exponent used to generate the samples.
r_max (float) – Largest separation to bin, in units of \(\ell_B\).
n_bins (int) – Number of radial bins.

Return type:

Returns:

r (ndarray of shape (n_bins,)) – Bin-center separations.
g (ndarray of shape (n_bins,)) – Pair correlation function; g[0] \approx 0 is the correlation hole, g \to 1 at large separation.

See also

laughlin_metropolis_sweep: Generates the equilibrium samples this function bins.

Examples

>>> import numpy as np
>>> from physicskit.statphys.core.monte_carlo import seed_numba_random
>>> seed_numba_random(0)
>>> rng = np.random.default_rng(0)
>>> N, m = 12, 3
>>> R0 = np.sqrt(2 * m * N)
>>> r0 = R0 * np.sqrt(rng.random(N))
>>> th0 = rng.uniform(0, 2 * np.pi, N)
>>> z = (r0 * np.cos(th0) + 1j * r0 * np.sin(th0)).astype(np.complex128)
>>> for _ in range(300):
...     _ = laughlin_metropolis_sweep(z, m=m, step=1.0)
>>> samples = np.array([laughlin_metropolis_sweep(z, m=m, step=1.0).copy() for _ in range(200)])
>>> r, g = laughlin_pair_correlation(samples, m=m, r_max=6.0, n_bins=30)
>>> bool(g[0] < 0.1)
True

physicskit.condensed.laughlin_radial_density(z_samples, r_max, n_bins=60)[source]#

Radially averaged particle density profile from sampled Laughlin configurations.

Parameters:

z_samples (ndarray) – Stack of particle-position snapshots (see laughlin_pair_correlation()).
r_max (float) – Largest radius to bin, in units of \(\ell_B\).
n_bins (int) – Number of radial bins.

Return type:

Returns:

r (ndarray of shape (n_bins,)) – Bin-center radii.
density (ndarray of shape (n_bins,)) – Mean particle density in each annulus, in units of \(\ell_B^{-2}\).

Notes

A large-\(N\) Laughlin droplet is an incompressible liquid: this profile is approximately flat at the bulk density \(\rho=1/(2\pi m \ell_B^2)\) out to the droplet radius, falling to zero over a few magnetic lengths at the edge – the real-space counterpart of the sharp gap that makes the plateau itself possible.

Examples

>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> N, m = 12, 3
>>> R0 = np.sqrt(2 * m * N)
>>> r0 = R0 * np.sqrt(rng.random(N))
>>> th0 = rng.uniform(0, 2 * np.pi, N)
>>> z = (r0 * np.cos(th0) + 1j * r0 * np.sin(th0)).astype(np.complex128)
>>> for _ in range(300):
...     _ = laughlin_metropolis_sweep(z, m=m, step=1.0)
>>> samples = np.array([laughlin_metropolis_sweep(z, m=m, step=1.0).copy() for _ in range(200)])
>>> r, density = laughlin_radial_density(samples, r_max=2 * R0, n_bins=20)
>>> bool(density[-1] < density[2])
True

physicskit.condensed.localization_length(psi, positions=None)[source]#

Estimate an eigenstate’s localization length from the decay of its envelope.

Fits \(\log|\psi_i|^2\) linearly against position on either side of the state’s peak and returns \(\xi = -2/\text{slope}\), averaged over both sides – the length scale over which a localized state’s probability density decays as \(|\psi(x)|^2 \sim e^{-2|x-x_0|/\xi}\). Returns inf for a state too extended (flat/non-monotonic envelope) for the fit to detect exponential decay.

Parameters:

psi (array_like) – Eigenvector components, indexed by site.
positions (array_like, optional) – Site coordinates. Defaults to integer indices 0, 1, ..., N-1.

Return type:

Returns:

float

Examples

An exactly exponential envelope recovers its localization length:

>>> import numpy as np
>>> x = np.arange(200)
>>> xi_true = 5.0
>>> psi = np.exp(-np.abs(x - 100) / xi_true)
>>> round(localization_length(psi), 4)
5.0

physicskit.condensed.lowest_eigenstates(H, k=2)[source]#

The \(k\) lowest eigenpairs of a Hermitian (sparse) Hamiltonian.

Uses dense numpy.linalg.eigh() for small matrices and the Lanczos method (scipy.sparse.linalg.eigsh(), which="SA") otherwise.

Parameters:

H (spmatrix | ndarray) – Hamiltonian, shape (dim, dim).
k (int) – Number of eigenpairs.

Return type:

tuple[ndarray, ndarray]

Returns:

energies (numpy.ndarray) – Ascending, shape (k,).
vectors (numpy.ndarray) – Eigenvectors as columns, shape (dim, k).

Examples

>>> H, _ = xxz_hamiltonian(4, n_up=2)
>>> E, _ = lowest_eigenstates(H, k=1)
>>> round(float(E[0]), 12)
-2.0

physicskit.condensed.magnetic_length(B, hbar=1.0, e=1.0)[source]#

Magnetic length \(\ell_B = \sqrt{\hbar/(eB)}\).

The natural length scale of the lowest Landau level’s cyclotron orbit.

Parameters:

B (float) – Magnetic field strength.
hbar (float) – Reduced Planck constant.
e (float) – Particle charge magnitude.

Return type:

Returns:

float

Examples

>>> magnetic_length(B=1.0)
1.0
>>> round(magnetic_length(B=4.0), 4)
0.5

physicskit.condensed.monatomic_chain_dispersion(k, K=1.0, m=1.0, a=1.0)[source]#

Phonon dispersion of a monatomic harmonic chain, \(\omega(k) = 2\sqrt{K/m}\,|\sin(ka/2)|\).

Masses \(m\) joined by springs \(K\) at spacing \(a\) (Ashcroft and Mermin, Solid State Physics, Eq. 22.29). At long wavelengths \(\omega \approx c|k|\) with sound speed \(c = a\sqrt{K/m}\); at the zone edge \(k = \pi/a\) the group velocity vanishes.

Parameters:

k (ndarray) – Wavenumber.
K (float) – Spring constant.
m (float) – Atomic mass.
a (float) – Lattice spacing.

Return type:

Returns:

numpy.ndarray

Examples

>>> import numpy as np
>>> float(monatomic_chain_dispersion(np.pi, K=4.0))
4.0

physicskit.condensed.plot_band_structure(hamiltonian_func, path_points, labels=None, n_per_segment=100, ax=None)[source]#

Plot bands along a piecewise-linear path through the Brillouin zone.

Parameters:

hamiltonian_func (callable) – A function H(k1, k2) returning the Bloch Hamiltonian.
path_points (sequence of tuple(float, float)) – High-symmetry points (e.g. \(\Gamma, K, M, \Gamma\)) in reduced crystal momentum, connected by straight segments.
labels (sequence of str, optional) – Tick labels for each point in path_points.
n_per_segment (int) – Number of sampled k-points per segment.
ax (matplotlib.axes.Axes, optional) – Axes to draw into; a new figure is created if omitted.

Returns:

fig (matplotlib.figure.Figure)
ax (matplotlib.axes.Axes)

Examples

>>> import numpy as np
>>> from physicskit.condensed.models import graphene_hamiltonian
>>> path = [(0, 0), (2 * np.pi / 3, 4 * np.pi / 3), (2 * np.pi / 3, 2 * np.pi / 3), (0, 0)]
>>> fig, ax = plot_band_structure(lambda k1, k2: graphene_hamiltonian(k1, k2), path, labels=["G", "K", "M", "G"])
>>> isinstance(fig, plt.Figure)
True
>>> sum(1 for line in ax.lines if line.get_color() == "C0")
2

physicskit.condensed.plot_berry_curvature(curvature, ax=None)[source]#

Plot a Berry-curvature field over the Brillouin zone as a heatmap.

Parameters:

curvature (ndarray) – Per-plaquette Berry curvature, as returned by physicskit.condensed.topology.compute_berry_curvature().
ax (matplotlib.axes.Axes, optional) – Axes to draw into; a new figure is created if omitted.

Returns:

fig (matplotlib.figure.Figure)
ax (matplotlib.axes.Axes)

Examples

>>> import numpy as np
>>> from physicskit.condensed.models import haldane_model
>>> from physicskit.condensed.topology import compute_berry_curvature
>>> H = lambda k1, k2: haldane_model(k1, k2, phi=np.pi / 2)
>>> F = compute_berry_curvature(H, grid_size=20, band_index=0)
>>> fig, ax = plot_berry_curvature(F)
>>> isinstance(fig, plt.Figure)
True

physicskit.condensed.plot_edge_state_density(ribbon_hamiltonian_func, k_parallel, ax=None)[source]#

Plot the real-space density of the mid-gap eigenstate closest to zero energy.

Useful for visualizing the exponential localization \(|\psi(x)|^2 \sim e^{-2x/\xi}\) of topological edge states exposed by physicskit.condensed.tight_binding.build_ribbon().

Parameters:

ribbon_hamiltonian_func (callable) – A function H(k_parallel) returning the finite ribbon Hamiltonian, as produced by build_ribbon().
k_parallel (array_like) – Momentum along the periodic direction(s) at which to evaluate the ribbon.
ax (matplotlib.axes.Axes, optional) – Axes to draw into; a new figure is created if omitted.

Returns:

fig (matplotlib.figure.Figure)
ax (matplotlib.axes.Axes)
density (ndarray) – Per-site probability density of the closest-to-zero eigenstate.

Examples

>>> import numpy as np
>>> from physicskit.condensed.models import ssh_lattice_hamiltonian
>>> from physicskit.condensed.tight_binding import build_ribbon
>>> H = ssh_lattice_hamiltonian(v=0.5, w=1.0)
>>> H_wire = build_ribbon(H, open_direction=0, n_cells=15)
>>> fig, ax, density = plot_edge_state_density(H_wire, k_parallel=[])
>>> bool(density[0] > density[len(density) // 2])
True

physicskit.condensed.plot_fermi_surface_3d(dispersion_func, mu=0.0, grid_size=40, k_range=3.141592653589793)[source]#

Render a 3D Fermi-surface isosurface \(\varepsilon(\mathbf{k}) = \mu\).

Parameters:

dispersion_func (callable) – A function eps(kx, ky, kz) returning the band energy.
mu (float) – Chemical potential (Fermi energy) defining the isosurface.
grid_size (int) – Number of samples along each of \(k_x, k_y, k_z\).
k_range (float) – Half-width of the cubic sampling box, [-k_range, k_range].

Returns:

plotly.graph_objects.Figure – A figure containing a single Isosurface trace.

Examples

>>> import numpy as np
>>> eps = lambda kx, ky, kz: -2 * (np.cos(kx) + np.cos(ky) + np.cos(kz))
>>> fig = plot_fermi_surface_3d(eps, mu=0.0, grid_size=20)
>>> fig.data[0].type
'isosurface'

physicskit.condensed.spin_chain_basis(N, n_up=None, parity=None)[source]#

Sorted computational basis states of an \(N\)-site spin-1/2 chain, optionally restricted to a sector.

Parameters:

N (int) – Number of sites.
n_up (int | None) – Keep only states with this many up spins, i.e. total \(S^z = n_\uparrow - N/2\).
parity (int | None) – Keep only states with this eigenvalue of \(P = \prod_i \sigma^z_i = (-1)^{n_\downarrow}\).

Return type:

Returns:

numpy.ndarray of int64 – Sorted basis states.

Examples

>>> spin_chain_basis(4, n_up=2)
array([ 3,  5,  6,  9, 10, 12])
>>> len(spin_chain_basis(6, parity=+1))
32

physicskit.condensed.square_lattice_dynamical_matrix(kx, ky, K1=1.0, K2=0.5, m=1.0, a=1.0)[source]#

Dynamical matrix of a square lattice with nearest (\(K_1\)) and next-nearest (\(K_2\)) central springs.

For central-force springs to neighbors \(\mathbf R\) with unit vectors \(\hat{\mathbf e}\), \(D(\mathbf k) = \sum_{\mathbf R} \frac{K_R}{m}(1 - \cos\mathbf k\cdot\mathbf R)\, \hat{\mathbf e}\hat{\mathbf e}^T\) (Ashcroft and Mermin, Eq. 22.58), which here gives

\[\begin{split}D = \frac{2K_1}{m}\begin{pmatrix} 1 - \cos k_xa & 0 \\ 0 & 1 - \cos k_ya \end{pmatrix} + \frac{K_2}{m}\left[(1 - \cos(k_x + k_y)a)\begin{pmatrix}1&1\\1&1\end{pmatrix} + (1 - \cos(k_x - k_y)a)\begin{pmatrix}1&-1\\-1&1\end{pmatrix}\right].\end{split}\]

The diagonal springs give the lattice its shear stiffness: with \(K_2 = 0\) the transverse branch is dispersionless along the axes.

Parameters:

kx (ndarray) – Wavevector components (broadcast together).
ky (ndarray) – Wavevector components (broadcast together).
K1 (float) – Nearest- and next-nearest-neighbor spring constants.
K2 (float) – Nearest- and next-nearest-neighbor spring constants.
m (float) – Atomic mass.
a (float) – Lattice constant.

Return type:

Returns:

numpy.ndarray – Shape broadcast(kx, ky).shape + (2, 2), real symmetric.

Examples

>>> D = square_lattice_dynamical_matrix(0.0, 0.0)
>>> D.tolist()
[[0.0, 0.0], [0.0, 0.0]]

physicskit.condensed.square_lattice_phonon_dispersion(kx, ky, K1=1.0, K2=0.5, m=1.0, a=1.0)[source]#

Phonon frequencies \(\omega_s(\mathbf k) = \sqrt{\lambda_s(D(\mathbf k))}\) of the square lattice.

Along \(\mathbf k \parallel \hat x\) the long-wavelength sound speeds are \(c_L = a\sqrt{(K_1 + K_2)/m}\) (longitudinal) and \(c_T = a\sqrt{K_2/m}\) (transverse).

Parameters:

kx (ndarray) – Wavevector components.
ky (ndarray) – Wavevector components.
K1 (float) – Spring constants.
K2 (float) – Spring constants.
m (float) – Atomic mass.
a (float) – Lattice constant.

Return type:

Returns:

numpy.ndarray – Shape broadcast(kx, ky).shape + (2,), ascending per wavevector.

Examples

At \((\pi, \pi)/a\) both branches have \(\omega^2 = 4K_1/m\):

>>> import numpy as np
>>> bool(np.allclose(square_lattice_phonon_dispersion(np.pi, np.pi), 2.0))
True

physicskit.condensed.ssh_hamiltonian(k, v=1.0, w=1.0)[source]#

Bloch Hamiltonian of the Su-Schrieffer-Heeger (SSH) dimerized chain.

Parameters:

k (float) – Reduced crystal momentum, periodic on \([0, 2\pi)\).
v (float) – Intracell hopping amplitude.
w (float) – Intercell hopping amplitude.

Return type:

Returns:

ndarray, shape (2, 2) – Bloch Hamiltonian \(H(k) = \begin{pmatrix}0 & v+we^{-ik}\\ v+we^{ik} & 0\end{pmatrix}\).

Notes

The chain is topologically nontrivial (hosts protected zero-energy edge states under open boundary conditions) when \(v < w\), and trivial when \(v > w\). See Also ssh_lattice_hamiltonian() for the real-space builder used to expose those edge states.

Examples

>>> import numpy as np
>>> H = ssh_hamiltonian(k=0.0, v=0.5, w=1.0)
>>> np.round(np.linalg.eigvalsh(H), 8)
array([-1.5,  1.5])

physicskit.condensed.ssh_lattice_hamiltonian(v=1.0, w=1.0)[source]#

Build the SSH chain as a Hamiltonian.

Parameters:

v (float) – Intracell hopping amplitude (A-B bond within a unit cell).
w (float) – Intercell hopping amplitude (B-A bond across the cell boundary).

Return type:

Returns:

Hamiltonian – A two-orbital-per-cell chain, suitable for build_ribbon() to expose edge zero modes.

See also

ssh_hamiltonian: Closed-form 2x2 Bloch Hamiltonian equivalent to this builder.

Examples

>>> import numpy as np
>>> H = ssh_lattice_hamiltonian(v=0.5, w=1.0)
>>> np.round(H.bands([0.0]), 8)
array([-1.5,  1.5])

physicskit.condensed.surface_dirac_hamiltonian(kx, ky, v_f=1.0)[source]#

Low-energy effective Hamiltonian of a single topological-insulator surface Dirac cone.

\[H_{\text{surf}}(\mathbf{k}) = \hbar v_F\left(k_x\sigma_y - k_y\sigma_x\right),\]

the massless, spin-momentum-locked Dirac fermion measured directly by ARPES on the Bi₂Se₃ and Bi₂Te₃ surfaces in 2008-2009: spin polarization locked perpendicular to momentum, with no Kramers-degenerate partner at the same energy and momentum (backscattering off nonmagnetic disorder is forbidden).

Parameters:

kx (float) – Momentum measured from the surface Dirac point.
ky (float) – Momentum measured from the surface Dirac point.
v_f (float) – Surface Fermi velocity.

Return type:

Returns:

ndarray, shape (2, 2) – Hermitian Hamiltonian with linear (Dirac) spectrum \(E_\pm(\mathbf{k}) = \pm v_F|\mathbf{k}|\).

Examples

>>> import numpy as np
>>> H = surface_dirac_hamiltonian(0.3, -0.4, v_f=2.0)
>>> np.round(np.linalg.eigvalsh(H), 8)
array([-1.,  1.])

physicskit.condensed.tfim_free_fermion_spectrum(N, J=1.0, h=1.0, parity=None)[source]#

Exact spectrum of the periodic TFIM from its Jordan-Wigner free-fermion solution.

The Jordan-Wigner transformation maps tfim_hamiltonian() onto free fermions, \(H = \sum_k \varepsilon_k (n_k - \tfrac12)\) with

\[\varepsilon_k = 2\sqrt{J^2 + h^2 - 2Jh\cos k}\]

(Pfeuty 1970, Eq. 2.15; Lieb, Schultz and Mattis, Ann. Phys. 16, 407 (1961)). The parity \(P = +1\) sector has antiperiodic fermions, \(k = \pm\pi(2n+1)/N\), and an even number of excitations; \(P = -1\) has periodic fermions, \(k = 2\pi n/N\), and an odd total fermion number. The unpaired modes \(k = 0, \pi\) carry the signed energy \(2(h - J\cos k)\) and count as bare fermions.

Parameters:

N (int) – Number of sites (keep \(N \lesssim 16\): all \(2^N\) levels are enumerated).
J (float) – Ising coupling.
h (float) – Transverse field.
parity (int | None) – Restrict to one sector (None for both).

Return type:

Returns:

numpy.ndarray – Sorted energies.

Examples

>>> import numpy as np
>>> H, _ = tfim_hamiltonian(5, J=1.0, h=0.6)
>>> bool(np.allclose(tfim_free_fermion_spectrum(5, 1.0, 0.6), np.linalg.eigvalsh(H.toarray())))
True

physicskit.condensed.tfim_ground_state_energy(N, J=1.0, h=1.0)[source]#

Exact ground-state energy of the periodic TFIM, the even-sector Bogoliubov vacuum.

\[E_0 = -\frac12 \sum_{k = \pm\pi(2n+1)/N} \varepsilon_k,\]

which tends to \(-\tfrac{N}{2\pi}\int_{-\pi}^{\pi}\varepsilon_k\,dk/2\) as \(N \to \infty\) (Pfeuty 1970, Eq. 3.2).

Parameters:

N (int) – Number of sites.
J (float) – Ising coupling.
h (float) – Transverse field.

Return type:

Returns:

float

Examples

>>> import numpy as np
>>> bool(np.isclose(tfim_ground_state_energy(8, 1.0, 1.0), -2 / np.sin(np.pi / 16)))
True

physicskit.condensed.tfim_hamiltonian(N, J=1.0, h=1.0, periodic=True, parity=None)[source]#

Sparse transverse-field Ising chain Hamiltonian.

\[H = -J \sum_{\langle ij \rangle} \sigma^x_i \sigma^x_j - h \sum_i \sigma^z_i\]

(Pfeuty, Ann. Phys. 57, 79 (1970), in Pauli matrices). The Ising coupling acts along \(x\) so that the conserved \(\mathbb{Z}_2\) parity \(P = \prod_i \sigma^z_i\) is diagonal in the computational basis; parity selects its sector. The chain is critical at \(h = J\).

Parameters:

N (int) – Number of sites.
J (float) – Ising coupling.
h (float) – Transverse field.
periodic (bool) – Periodic (ring) or open boundary conditions.
parity (int | None) – Parity sector (None for the full space).

Return type:

tuple[csr_matrix, ndarray]

Returns:

H (scipy.sparse.csr_matrix) – The Hamiltonian in the sector basis.
basis (numpy.ndarray) – The sector’s basis states.

Examples

At zero field the ground state is the doubly degenerate ferromagnet with energy \(-NJ\) on a ring:

>>> import numpy as np
>>> H, _ = tfim_hamiltonian(4, J=1.0, h=0.0)
>>> np.round(np.linalg.eigvalsh(H.toarray())[:2], 12)
array([-4., -4.])

physicskit.condensed.topological_insulator_3d_hamiltonian(kx, ky, kz, m, t=1.0)[source]#

Bloch Hamiltonian of the minimal cubic-lattice 3D topological insulator model.

Parameters:

kx (float) – Reduced crystal momenta, each periodic on \([0, 2\pi)\).
ky (float) – Reduced crystal momenta, each periodic on \([0, 2\pi)\).
kz (float) – Reduced crystal momenta, each periodic on \([0, 2\pi)\).
m (float) – Band mass. The model is a strong topological insulator for \(1 < |m/t| < 3\), and trivial for \(|m/t| < 1\) or \(|m/t| > 3\).
t (float) – Hopping amplitude setting the bulk bandwidth.

Return type:

Returns:

ndarray, shape (4, 4) – Hermitian Bloch Hamiltonian.

Examples

>>> import numpy as np
>>> H = topological_insulator_3d_hamiltonian(0.3, -0.7, 0.5, m=1.0)
>>> np.allclose(H, H.conj().T)
True
>>> H0 = topological_insulator_3d_hamiltonian(0.0, 0.0, 0.0, m=1.0, t=1.0)
>>> np.round(np.linalg.eigvalsh(H0), 8)
array([-4., -4.,  4.,  4.])

physicskit.condensed.topological_insulator_3d_slab_hamiltonian(kx, ky, n_layers, m, t=1.0)[source]#

Slab Hamiltonian: periodic in x, y, open (finite) along z.

Truncating the z direction exposes the top and bottom (001) surfaces. Diagonalizing at each (kx, ky) and looking for mid-gap states localized at the outer layers reveals the surface Dirac cone directly, with no separate topological-invariant calculation needed.

Parameters:

kx (float) – Reduced crystal momenta, periodic on \([0, 2\pi)\).
ky (float) – Reduced crystal momenta, periodic on \([0, 2\pi)\).
n_layers (int) – Number of atomic layers stacked along z.
m (float) – Band mass (see topological_insulator_3d_hamiltonian()).
t (float) – Hopping amplitude.

Return type:

Returns:

ndarray, shape (4 * n_layers, 4 * n_layers) – Hermitian, open-boundary slab Hamiltonian.

Examples

In the strong-topological-insulator regime (\(1 < |m/t| < 3\)), a thick slab has a state pinned to exactly zero energy at the surface Brillouin zone center, localized entirely on one outer layer – the surface Dirac point – that a topologically trivial slab (\(|m/t| < 1\)) lacks:

>>> import numpy as np
>>> gap = np.linalg.eigvalsh(topological_insulator_3d_slab_hamiltonian(0.0, 0.0, n_layers=40, m=-2.0))
>>> bool(np.min(np.abs(gap)) < 1e-8)
True
>>> gap_trivial = np.linalg.eigvalsh(topological_insulator_3d_slab_hamiltonian(0.0, 0.0, n_layers=40, m=0.0))
>>> bool(np.min(np.abs(gap_trivial)) > 0.5)
True

physicskit.condensed.weyl_node_locations(m, t=1.0)[source]#

Exact \(k_z\) positions of the two Weyl nodes on the \(k_z\) axis.

Solves \(m - t(2 + \cos k_z) = 0\) – the mass term of weyl_semimetal_hamiltonian() at \(k_x=k_y=0\) – for \(k_z\).

Parameters:

m (float) – Band mass.
t (float) – Hopping amplitude.

Return type:

ndarray | None

Returns:

tuple of float or None – (-k0, +k0) with \(k_0 \in (0, \pi)\), or None if \(m/t\) lies outside \((1, 3)\) (no nodes on this axis).

Examples

>>> import numpy as np
>>> k0 = weyl_node_locations(m=2.0, t=1.0)
>>> np.round(k0, 8)
array([-1.57079633,  1.57079633])
>>> weyl_node_locations(m=4.0, t=1.0) is None
True

physicskit.condensed.weyl_semimetal_hamiltonian(kx, ky, kz, m, t=1.0)[source]#

Bloch Hamiltonian of the minimal two-band cubic-lattice Weyl semimetal.

Parameters:

kx (float) – Reduced crystal momenta, each periodic on \([0, 2\pi)\).
ky (float) – Reduced crystal momenta, each periodic on \([0, 2\pi)\).
kz (float) – Reduced crystal momenta, each periodic on \([0, 2\pi)\).
m (float) – Band mass. Two Weyl nodes exist on the \(k_z\) axis for \(1 < m/t < 3\) (see weyl_node_locations()); outside that window the model is fully gapped.
t (float) – Hopping amplitude setting the bulk bandwidth.

Return type:

Returns:

ndarray, shape (2, 2) – Hermitian Bloch Hamiltonian.

See also

weyl_node_locations: Exact \(k_z\) positions of the two Weyl nodes.
weyl_semimetal_slab_hamiltonian: Open-boundary slab exposing Fermi-arc surface states.

Examples

At m/t = 2 the two nodes sit at \(k_z = \pm\pi/2\), where the bulk gap closes exactly:

>>> import numpy as np
>>> H = weyl_semimetal_hamiltonian(0.0, 0.0, np.pi / 2, m=2.0, t=1.0)
>>> np.round(np.linalg.eigvalsh(H), 8)
array([0., 0.])
>>> H_away = weyl_semimetal_hamiltonian(np.pi, np.pi, np.pi, m=2.0, t=1.0)
>>> np.round(np.linalg.eigvalsh(H_away), 8)
array([-5.,  5.])

physicskit.condensed.weyl_semimetal_slab_hamiltonian(ky, kz, n_layers, m, t=1.0)[source]#

Slab Hamiltonian: periodic in y and z, open (finite) along x.

Truncating the x direction exposes two open surfaces whose surface Brillouin zone is the \((k_y, k_z)\) plane the two Weyl nodes project onto at \((0, \pm k_0)\). Diagonalizing at fixed \(k_z\) and scanning \(k_y\) reveals, for \(|k_z| < k_0\), a single chiral mode crossing zero energy (the Fermi arc), and a full gap for \(|k_z| > k_0\).

Parameters:

ky (float) – Reduced crystal momenta, periodic on \([0, 2\pi)\).
kz (float) – Reduced crystal momenta, periodic on \([0, 2\pi)\).
n_layers (int) – Number of layers stacked along x.
m (float) – Band mass (see weyl_semimetal_hamiltonian()).
t (float) – Hopping amplitude.

Return type:

Returns:

ndarray, shape (2 * n_layers, 2 * n_layers) – Hermitian, open-boundary slab Hamiltonian.

Examples

Between the two Weyl nodes (\(k_z=0\)), a thick slab has a state pinned near zero energy at some \(k_y\) – the Fermi-arc surface mode – absent outside the node range (\(k_z=\pi\)):

>>> import numpy as np
>>> ky_grid = np.linspace(0, 2 * np.pi, 200, endpoint=False)
>>> gaps_inside = [np.min(np.abs(np.linalg.eigvalsh(weyl_semimetal_slab_hamiltonian(ky, 0.0, 40, m=2.0)))) for ky in ky_grid]
>>> bool(min(gaps_inside) < 1e-6)
True
>>> gaps_outside = [np.min(np.abs(np.linalg.eigvalsh(weyl_semimetal_slab_hamiltonian(ky, np.pi, 40, m=2.0)))) for ky in ky_grid]
>>> bool(min(gaps_outside) > 0.5)
True

physicskit.condensed.xxz_hamiltonian(N, J=1.0, delta=1.0, h=0.0, periodic=True, n_up=None)[source]#

Sparse spin-1/2 Heisenberg XXZ chain Hamiltonian.

\[H = J \sum_{\langle ij \rangle} \left( S^x_i S^x_j + S^y_i S^y_j + \Delta\, S^z_i S^z_j \right) - h \sum_i S^z_i, \qquad \mathbf{S} = \tfrac12 \boldsymbol{\sigma}.\]

The flip-flop term \(S^x S^x + S^y S^y = \tfrac12(S^+S^- + S^-S^+)\) conserves total \(S^z\), so passing n_up builds the Hamiltonian directly in that magnetization sector. \(\Delta = 1\) is the isotropic Heisenberg (XXX) antiferromagnet solved by Bethe (Z. Phys. 71, 205 (1931)).

Parameters:

N (int) – Number of sites.
J (float) – Exchange coupling (\(J > 0\) antiferromagnetic).
delta (float) – Ising anisotropy \(\Delta\).
h (float) – Longitudinal magnetic field.
periodic (bool) – Periodic (ring) or open boundary conditions.
n_up (int | None) – Number of up spins of the sector (None for the full space).

Return type:

tuple[csr_matrix, ndarray]

Returns:

H (scipy.sparse.csr_matrix) – The Hamiltonian in the sector basis.
basis (numpy.ndarray) – The sector’s basis states, from spin_chain_basis().

Examples

The two-site singlet has energy \(-3J/4\):

>>> import numpy as np
>>> H, basis = xxz_hamiltonian(2, periodic=False, n_up=1)
>>> np.linalg.eigvalsh(H.toarray())
array([-0.75,  0.25])

physicskit.condensed.z2_invariant(hamiltonian_func, grid_size=30, spin_block=(0, 2))[source]#

Compute the \(\mathbb{Z}_2\) invariant of an \(s_z\)-conserving spinful model.

Valid whenever hamiltonian_func is block-diagonal in spin (e.g. the Kane-Mele or BHZ models with Rashba coupling set to zero), in which case the \(\mathbb{Z}_2\) invariant reduces to the spin-up sector’s Chern number modulo 2 (Sheng et al. / spin-Chern-number approach).

Parameters:

hamiltonian_func (callable) – A function H(k1, k2) returning the full (2N, 2N) Bloch Hamiltonian with the spin-up block occupying spin_block.
grid_size (int) – Discretization resolution for the Chern-number calculation.
spin_block (tuple) – (start, stop) slice selecting the spin-up block’s rows/columns.

Return type:

int

Returns:

int – \(\mathbb{Z}_2\) invariant, 0 (trivial) or 1 (topological).

Examples

>>> from physicskit.condensed.models import kane_mele_hamiltonian
>>> trivial = lambda k1, k2: kane_mele_hamiltonian(k1, k2, lambda_so=0.0, lambda_v=0.3)
>>> topological = lambda k1, k2: kane_mele_hamiltonian(k1, k2, lambda_so=0.06)
>>> z2_invariant(trivial, grid_size=20)
0
>>> z2_invariant(topological, grid_size=20)
1

physicskit.condensed.zak_phase(hamiltonian_func_1d, grid_size=200, band_index=0)[source]#

Compute the Zak phase of a 1D band: the Berry phase accumulated across the BZ.

Parameters:

hamiltonian_func_1d (callable) – A function H(k) returning the (N, N) complex Bloch Hamiltonian at reduced momentum k, periodic on \([0, 2\pi)\).
grid_size (int) – Number of discretization steps around the 1D Brillouin zone.
band_index (int) – Band index (0 = lowest energy).

Return type:

Returns:

float – The Zak phase, wrapped to \((-\pi, \pi]\). For a chiral-symmetric model such as SSH, this is quantized to \(0\) (trivial) or \(\pi\) (topological).

Examples

>>> import numpy as np
>>> from physicskit.condensed.models import ssh_hamiltonian
>>> topological = lambda k: ssh_hamiltonian(k, v=0.5, w=1.0)
>>> trivial = lambda k: ssh_hamiltonian(k, v=1.0, w=0.5)
>>> round(abs(zak_phase(topological)), 4)
3.1416
>>> round(abs(zak_phase(trivial)), 4)
0.0

Generic tight-binding lattice and Bloch Hamiltonian construction.

This module provides the building blocks that physicskit.condensed.models uses to assemble specific band models (graphene, the Haldane model, …): a Lattice describing the Bravais lattice and orbital basis, and a Hamiltonian that accumulates real-space hoppings and Bloch-sums them into \(H(\mathbf{k})\).

Convention#

Bloch Hamiltonians in this package use the periodic gauge: a hopping term with integer cell offset \(\mathbf{R} = (n_1, n_2, \dots)\) (expressed in units of the primitive lattice vectors) contributes a phase \(e^{i \mathbf{k} \cdot \mathbf{R}}\) where \(\mathbf{k}\) is the reduced crystal momentum conjugate to the primitive-vector indices, each component periodic on \([0, 2\pi)\). This intracell-position-free convention keeps \(H(\mathbf{k})\) manifestly periodic on the Brillouin zone torus, which is exactly what the Fukui-Hatsugai-Suzuki Chern-number algorithm in physicskit.condensed.topology needs.

class physicskit.condensed.tight_binding.Hamiltonian(lattice, onsite=None)[source]#

Bases: object

Real-space tight-binding Hamiltonian, Bloch-summed into \(H(\mathbf{k})\).

Parameters:

lattice (Lattice) – The underlying Bravais lattice and orbital basis.
onsite (array_like, shape (n_orbitals,), optional) – Onsite energies. Defaults to zero for every orbital.

Examples

Nearest-neighbor graphene, reproducing the linear Dirac dispersion:

>>> import numpy as np
>>> lat = Lattice.honeycomb()
>>> H = Hamiltonian(lat)
>>> t = 1.0
>>> H.add_hopping(0, 1, (0, 0), t)
>>> H.add_hopping(0, 1, (-1, 0), t)
>>> H.add_hopping(0, 1, (0, -1), t)
>>> K = np.array([2 * np.pi / 3, 4 * np.pi / 3])
>>> np.round(np.linalg.eigvalsh(H.bloch(K)), 8)
array([-0.,  0.])

add_hopping(i, j, cell_offset, amplitude)[source]#

Add a hopping term \(t\, c_i^\dagger(0) c_j(\mathbf{R})\) (+ h.c.).

Parameters:

i (int) – Orbital indices within the unit cell.
j (int) – Orbital indices within the unit cell.
cell_offset (tuple of int) – Integer offset \(\mathbf{R}\) (in primitive-vector units) of orbital j’s cell relative to orbital i’s cell.
amplitude (complex) – Hopping amplitude. The Hermitian conjugate term is added automatically; do not add both a bond and its reverse.

Raises:

ValueError – If (i, j, cell_offset) describes a diagonal onsite term (i == j and cell_offset all zero); use onsite instead.

Return type:

None

bands(k)[source]#

Eigenvalues of bloch() at k, sorted ascending.

Return type:: ndarray

Examples

>>> lat = Lattice.chain()
>>> H = Hamiltonian(lat)
>>> H.add_hopping(0, 0, (1,), 1.0)
>>> import numpy as np
>>> np.round(H.bands([np.pi]), 8)
array([-2.])

bloch(k)[source]#

Evaluate the Bloch Hamiltonian at reduced crystal momentum k.

Parameters:: k (array_like, shape (dim,)) – Reduced crystal momentum, each component periodic on \([0, 2\pi)\).
Return type:: ndarray
Returns:: ndarray, shape (n_orbitals, n_orbitals) – Hermitian Bloch Hamiltonian matrix.

class physicskit.condensed.tight_binding.Lattice(lattice_vectors, orbitals, labels=<factory>)[source]#

Bases: object

A Bravais lattice with a basis of orbitals.

Parameters:

lattice_vectors (ndarray) – Primitive lattice vectors as rows, in Cartesian coordinates.
orbitals (ndarray) – Orbital positions in fractional coordinates of the primitive cell.
labels (list) – Human-readable names for each orbital (defaults to orb0, orb1, …).

Examples

>>> lat = Lattice.honeycomb()
>>> lat.n_orbitals
2
>>> lat.dim
2

cartesian_orbitals()[source]#

Return orbital positions in Cartesian coordinates.

Return type:: ndarray
Returns:: ndarray, shape (n_orbitals, dim)

Examples

>>> lat = Lattice.square()
>>> lat.cartesian_orbitals()
array([[0., 0.]])

classmethod chain(a=1.0)[source]#

1D monatomic chain with lattice constant a.

Return type:: Lattice
Parameters:: a (float)

classmethod cubic(a=1.0)[source]#

3D simple cubic lattice, one orbital per cell.

Return type:: Lattice
Parameters:: a (float)

property dim: int#

Spatial dimension (1, 2, or 3).

Type:: int

classmethod honeycomb(a=1.0)[source]#

2D honeycomb lattice (graphene structure), sublattices A and B.

a is the lattice constant (nearest-neighbor bond length is a/sqrt(3)).

Return type:: Lattice
Parameters:: a (float)

classmethod kagome(a=1.0)[source]#

2D kagome lattice, three orbitals per cell (edge midpoints of a triangular lattice).

Return type:: Lattice
Parameters:: a (float)

labels: list#

lattice_vectors: ndarray#

property n_orbitals: int#

Number of orbitals per unit cell.

Type:: int

orbitals: ndarray#

reciprocal_vectors()[source]#

Return primitive reciprocal lattice vectors \(\mathbf{b}_i\).

Satisfies \(\mathbf{a}_i \cdot \mathbf{b}_j = 2\pi\delta_{ij}\).

Return type:: ndarray
Returns:: ndarray, shape (dim, dim)

Examples

>>> lat = Lattice.square(a=1.0)
>>> np.allclose(lat.reciprocal_vectors(), 2 * np.pi * np.eye(2))
True

classmethod square(a=1.0)[source]#

2D square lattice, one orbital per cell.

Return type:: Lattice
Parameters:: a (float)

classmethod triangular(a=1.0)[source]#

2D triangular lattice, one orbital per cell.

Return type:: Lattice
Parameters:: a (float)

physicskit.condensed.tight_binding.apply_peierls_phase(positions, hoppings, flux_quanta_per_plaquette, area_per_plaquette)[source]#

Apply a Peierls substitution phase to real-space hoppings for a uniform magnetic field.

In the Landau gauge \(\mathbf{A} = (-By, 0)\), a hopping from site at r_i to site at r_j acquires the phase \(\exp\!\big(i\frac{2\pi\Phi}{\Phi_0}\,\bar{y}\,(x_j-x_i)/a^2\big)\) where \(\Phi/\Phi_0\) is the flux per plaquette in units of the flux quantum, evaluated at the bond midpoint \(\bar y\).

Parameters:

positions (ndarray) – Real-space Cartesian coordinates of every site in a finite lattice.
hoppings (list of tuple(int, int, complex)) – (i, j, amplitude) real-space bonds (indices into positions).
flux_quanta_per_plaquette (float) – Magnetic flux per unit-cell plaquette, in units of the flux quantum \(\Phi_0 = h/e\).
area_per_plaquette (float) – Real-space area of one plaquette, used to convert flux density to the vector potential prefactor.

Returns:

list of tuple(int, int, complex) – The same bonds with Peierls phases multiplied into the amplitude.

Examples

>>> import numpy as np
>>> positions = np.array([[0.0, 0.0], [1.0, 0.0]])
>>> bonds = [(0, 1, 1.0)]
>>> out = apply_peierls_phase(positions, bonds, flux_quanta_per_plaquette=0.25, area_per_plaquette=1.0)
>>> bool(abs(out[0][2] - 1.0) < 1e-12)
True

physicskit.condensed.tight_binding.build_finite_cluster(hamiltonian, n_cells, keep=None)[source]#

Build a fully finite (open in every direction) real-space cluster.

Unlike build_ribbon(), which stays periodic along all but one lattice direction, this truncates every direction, producing a genuinely finite flake with real Cartesian site positions – what a boundary-sensitive real-space plot (e.g. plot_lattice_structure()) needs.

Parameters:

hamiltonian (Hamiltonian) – A Hamiltonian built on a Lattice of any dimension.
n_cells (int or sequence of int) – Number of unit cells spanned along each lattice direction (the bounding-box shape). A single int is broadcast to every direction.
keep (callable, optional) – keep(cell, orbital, position) -> bool predicate selecting which sites within the bounding box survive, for carving non-rectangular shapes (e.g. a disk) out of it. cell is the integer cell-index tuple, orbital the orbital index within the cell, and position its Cartesian coordinate. Defaults to keeping every site in the bounding box.

Returns:

H (ndarray, shape (n_sites, n_sites)) – Dense, open-boundary real-space Hamiltonian.
positions (ndarray, shape (n_sites, dim)) – Cartesian coordinates of the surviving sites, in the same order as H’s rows/columns.
bonds (list of tuple(int, int, complex)) – Surviving real-space hoppings (i, j, amplitude), indexing into positions/H.

Examples

A 6-site open SSH chain (3 unit cells of 2 orbitals each) has exactly 5 nearest-neighbor bonds, one fewer than a periodic ring would:

>>> import numpy as np
>>> from physicskit.condensed.models import ssh_lattice_hamiltonian
>>> H_bulk = ssh_lattice_hamiltonian(v=0.5, w=1.0)
>>> H, positions, bonds = build_finite_cluster(H_bulk, n_cells=3)
>>> H.shape
(6, 6)
>>> len(bonds)
5
>>> positions.shape
(6, 1)

physicskit.condensed.tight_binding.build_ribbon(hamiltonian, open_direction, n_cells)[source]#

Build a ribbon/slab: periodic in-plane, open (finite) along open_direction.

Truncates the periodic boundary condition along one primitive-lattice direction, producing a quasi-1D (2D lattice) or quasi-2D (3D lattice) strip Hamiltonian as a function of the remaining reduced momenta. This exposes edge/surface states localized at the two open boundaries.

Parameters:

hamiltonian (Hamiltonian) – A Hamiltonian built on a 2D (or higher) Lattice.
open_direction (int) – Index of the primitive-lattice direction to truncate.
n_cells (int) – Number of unit cells stacked along open_direction.

Returns:

callable – A function H_ribbon(k_parallel) returning the (n_cells * n_orbitals, n_cells * n_orbitals) ribbon Hamiltonian, where k_parallel is a reduced-momentum vector with one fewer component than the bulk lattice (the remaining periodic directions).

Examples

SSH chain cut open into a finite 20-site wire; the topological phase (v < w) hosts a mid-gap zero mode:

>>> import numpy as np
>>> from physicskit.condensed.models import ssh_lattice_hamiltonian
>>> H = ssh_lattice_hamiltonian(v=0.5, w=1.0)
>>> H_wire = build_ribbon(H, open_direction=0, n_cells=30)
>>> spectrum = np.linalg.eigvalsh(H_wire(np.array([])))
>>> bool(np.any(np.abs(spectrum) < 1e-6))
True

Berry curvature, Chern numbers, Zak phase, and Z2 invariants.

The workhorse here is the Fukui-Hatsugai-Suzuki (FHS) lattice algorithm, which computes an exactly quantized, gauge-invariant Chern number from a finite grid of Bloch eigenvectors – no smooth gauge choice required. The per-plaquette flux accumulation is numba-jitted since it is the inner loop of a double sum over the whole Brillouin-zone grid.

physicskit.condensed.topology.compute_berry_curvature(hamiltonian_func, grid_size=50, band_index=0)[source]#

Compute the discretized Berry curvature of one band over the Brillouin zone.

Uses the Fukui-Hatsugai-Suzuki (FHS) link-variable formula, which is gauge invariant plaquette by plaquette and needs no smooth choice of eigenvector phase.

Parameters:

hamiltonian_func (callable) – A function H(k1, k2) returning the (N, N) complex Bloch Hamiltonian at reduced crystal momentum (k1, k2), each periodic on \([0, 2\pi)\) (see physicskit.condensed.tight_binding).
grid_size (int) – Number of plaquettes along each of \(k_1, k_2\).
band_index (int) – Band index (0 = lowest energy) to compute the curvature for.

Return type:

Returns:

ndarray, shape (grid_size, grid_size) – Berry flux through each plaquette, in \((-\pi, \pi]\). Summing and dividing by \(2\pi\) gives the band’s Chern number.

See also

compute_chern_number: Integrates this curvature over the full zone for every band.

Examples

>>> import numpy as np
>>> from physicskit.condensed.models import haldane_model
>>> H_func = lambda k1, k2: haldane_model(k1, k2, t=1.0, t2=0.2, phi=np.pi / 2, M=0.0)
>>> F = compute_berry_curvature(H_func, grid_size=30, band_index=0)
>>> F.shape
(30, 30)
>>> round(float(F.sum() / (2 * np.pi)))
1

physicskit.condensed.topology.compute_chern_number(hamiltonian_func, grid_size=50)[source]#

Compute the Chern number of every band using the Fukui-Hatsugai-Suzuki method.

Parameters:

hamiltonian_func (callable) – A function H(k1, k2) returning an (N, N) complex Bloch Hamiltonian, using the reduced-momentum convention of physicskit.condensed.tight_binding.
grid_size (int) – Discretization resolution for \(k_1\) and \(k_2\) across \([0, 2\pi)\).

Return type:

list

Returns:

list of int – Chern integer for each energy band, ordered from lowest to highest energy.

Examples

>>> import numpy as np
>>> from physicskit.condensed.models import haldane_model
>>> H_func = lambda k1, k2: haldane_model(k1, k2, t=1.0, t2=0.2, phi=np.pi / 2, M=0.0)
>>> chern_numbers = compute_chern_number(H_func, grid_size=30)
>>> print(chern_numbers)
[1, -1]

physicskit.condensed.topology.z2_invariant(hamiltonian_func, grid_size=30, spin_block=(0, 2))[source]#

Compute the \(\mathbb{Z}_2\) invariant of an \(s_z\)-conserving spinful model.

Valid whenever hamiltonian_func is block-diagonal in spin (e.g. the Kane-Mele or BHZ models with Rashba coupling set to zero), in which case the \(\mathbb{Z}_2\) invariant reduces to the spin-up sector’s Chern number modulo 2 (Sheng et al. / spin-Chern-number approach).

Parameters:

hamiltonian_func (callable) – A function H(k1, k2) returning the full (2N, 2N) Bloch Hamiltonian with the spin-up block occupying spin_block.
grid_size (int) – Discretization resolution for the Chern-number calculation.
spin_block (tuple) – (start, stop) slice selecting the spin-up block’s rows/columns.

Return type:

int

Returns:

int – \(\mathbb{Z}_2\) invariant, 0 (trivial) or 1 (topological).

Examples

>>> from physicskit.condensed.models import kane_mele_hamiltonian
>>> trivial = lambda k1, k2: kane_mele_hamiltonian(k1, k2, lambda_so=0.0, lambda_v=0.3)
>>> topological = lambda k1, k2: kane_mele_hamiltonian(k1, k2, lambda_so=0.06)
>>> z2_invariant(trivial, grid_size=20)
0
>>> z2_invariant(topological, grid_size=20)
1

physicskit.condensed.topology.zak_phase(hamiltonian_func_1d, grid_size=200, band_index=0)[source]#

Compute the Zak phase of a 1D band: the Berry phase accumulated across the BZ.

Parameters:

hamiltonian_func_1d (callable) – A function H(k) returning the (N, N) complex Bloch Hamiltonian at reduced momentum k, periodic on \([0, 2\pi)\).
grid_size (int) – Number of discretization steps around the 1D Brillouin zone.
band_index (int) – Band index (0 = lowest energy).

Return type:

Returns:

float – The Zak phase, wrapped to \((-\pi, \pi]\). For a chiral-symmetric model such as SSH, this is quantized to \(0\) (trivial) or \(\pi\) (topological).

Examples

>>> import numpy as np
>>> from physicskit.condensed.models import ssh_hamiltonian
>>> topological = lambda k: ssh_hamiltonian(k, v=0.5, w=1.0)
>>> trivial = lambda k: ssh_hamiltonian(k, v=1.0, w=0.5)
>>> round(abs(zak_phase(topological)), 4)
3.1416
>>> round(abs(zak_phase(trivial)), 4)
0.0

Landau levels: the 1930 quantization of a charged particle in a uniform magnetic field.

Lev Landau showed that the classically continuous cyclotron motion of a charged particle in a uniform magnetic field \(B\) quantizes into discrete, macroscopically degenerate levels

\[E_n = \hbar\omega_c\left(n+\tfrac12\right), \qquad \omega_c = \frac{eB}{m}, \qquad n = 0, 1, 2, \dots,\]

exactly the harmonic-oscillator spectrum, with \(\omega_c\) set by the field rather than a spring constant. Each level’s macroscopic degeneracy per unit area, \(n_B = 1/(2\pi \ell_B^2)\) with magnetic length \(\ell_B = \sqrt{\hbar/(eB)}\), is the microscopic origin of orbital (Landau) diamagnetism, and – once a 2D electron gas and disorder-broadened levels are added – the direct ancestor of the integer quantum Hall effect (physicskit.condensed.topology) fifty years later.

Uses the same natural-unit convention (\(\hbar=e=1\)) as the rest of physicskit.condensed; pass m, e, and hbar explicitly to work in physical units instead. See physicskit.condensed.tight_binding.apply_peierls_phase() for the lattice (Peierls-substitution) route to the same physics.

physicskit.condensed.landau_levels.cyclotron_frequency(B, m=1.0, e=1.0)[source]#

Cyclotron frequency \(\omega_c = eB/m\).

Parameters:

B (float) – Magnetic field strength.
m (float) – Particle mass.
e (float) – Particle charge magnitude.

Return type:

Returns:

float

Examples

>>> cyclotron_frequency(B=2.0, m=0.5)
4.0

physicskit.condensed.landau_levels.filling_factor(density, B, hbar=1.0, e=1.0)[source]#

Landau-level filling factor \(\nu = n_e / n_B = n_e h/(eB)\).

The number of filled Landau levels (generally non-integer) for a 2D electron density density. Integer \(\nu\) is the condition for an incompressible quantum Hall plateau in physicskit.condensed.topology.

Parameters:

density (float) – 2D electron number density (particles per unit area).
B (float) – Magnetic field strength.
hbar (float) – Reduced Planck constant.
e (float) – Particle charge magnitude.

Return type:

Returns:

float

Examples

>>> round(filling_factor(density=2.0, B=1.0), 6)
12.566371

physicskit.condensed.landau_levels.landau_degeneracy(area, B, hbar=1.0, e=1.0)[source]#

Number of degenerate single-particle states per Landau level.

Each Landau level holds \(n_B \cdot \text{area}\) states, where \(n_B = 1/(2\pi\ell_B^2) = eB/(2\pi\hbar) = B/\Phi_0\) is one state per flux quantum \(\Phi_0 = 2\pi\hbar/e\) threading the sample – the origin of the macroscopic degeneracy behind Landau diamagnetism.

Parameters:

area (float) – Real-space sample area.
B (float) – Magnetic field strength.
hbar (float) – Reduced Planck constant.
e (float) – Particle charge magnitude.

Return type:

Returns:

float

Examples

>>> landau_degeneracy(area=2 * np.pi, B=1.0)
1.0

physicskit.condensed.landau_levels.landau_density_of_states(energies, B, m=1.0, e=1.0, hbar=1.0, n_max=20, broadening=0.05)[source]#

Disorder-broadened density of states per unit area.

Replaces each infinitely sharp level \(n_B\,\delta(E - E_n)\) with a Gaussian of width broadening, the standard phenomenological model for the disorder- or finite-lifetime-broadened Landau levels seen in a real 2D electron gas (and needed to have a finite conductivity between the exactly quantized Hall plateaus).

Parameters:

energies (array_like) – Energies at which to evaluate the density of states.
B (float) – Magnetic field strength.
m (float) – Particle mass.
e (float) – Particle charge magnitude.
hbar (float) – Reduced Planck constant.
n_max (int) – Highest Landau index included in the sum.
broadening (float) – Gaussian standard deviation of each broadened level.

Return type:

Returns:

ndarray – Density of states per unit area, same shape as energies.

Examples

>>> dos = landau_density_of_states([0.5], B=1.0, n_max=0, broadening=0.1)
>>> round(float(dos[0]), 4)
0.6349

physicskit.condensed.landau_levels.landau_level_energies(n_max, B, m=1.0, e=1.0, hbar=1.0)[source]#

Landau energy levels \(E_n = \hbar\omega_c(n+\tfrac12)\) for \(n=0,\dots,n_{max}\).

Parameters:

n_max (int) – Highest Landau index to include.
B (float) – Magnetic field strength.
m (float) – Particle mass.
e (float) – Particle charge magnitude.
hbar (float) – Reduced Planck constant.

Return type:

Returns:

ndarray, shape (n_max + 1,) – Equally spaced energies \(E_0 < E_1 < \dots < E_{n_{max}}\), spacing \(\hbar\omega_c\).

See also

landau_degeneracy: Number of degenerate single-particle states per level.

Examples

>>> landau_level_energies(n_max=2, B=1.0)
array([0.5, 1.5, 2.5])

physicskit.condensed.landau_levels.magnetic_length(B, hbar=1.0, e=1.0)[source]#

Magnetic length \(\ell_B = \sqrt{\hbar/(eB)}\).

The natural length scale of the lowest Landau level’s cyclotron orbit.

Parameters:

B (float) – Magnetic field strength.
hbar (float) – Reduced Planck constant.
e (float) – Particle charge magnitude.

Return type:

Returns:

float

Examples

>>> magnetic_length(B=1.0)
1.0
>>> round(magnetic_length(B=4.0), 4)
0.5

Ginzburg-Landau theory: the 1950 phenomenological theory of continuous phase transitions.

Landau and Ginzburg proposed that a superconductor (or any system undergoing a continuous symmetry-breaking transition) is described entirely by a complex order parameter \(\psi(\mathbf{r})\) – here, the Cooper-pair condensate wavefunction – through a free-energy functional expanded in powers of \(\psi\) and its gradient,

\[f[\psi] = a|\psi|^2 + \frac{b}{2}|\psi|^4 + \frac{\hbar^2}{2m}|\nabla\psi|^2,\]

built from symmetry alone, with no reference to the microscopic pairing mechanism BCS would supply seven years later. Minimizing \(f\) gives a uniform condensate density \(|\psi_0|^2 = -a/b\) below the transition (\(a<0\)) and two emergent length scales – the coherence length \(\xi\) over which \(\psi\) heals back from a boundary, and the magnetic penetration depth \(\lambda\) – whose ratio \(\kappa=\lambda/\xi\) (the Ginzburg-Landau parameter) alone decides whether a superconductor is Type I or Type II.

Uses the same natural-unit convention (\(\hbar=e=1\), and here also \(\mu_0=1\)) as the rest of physicskit.condensed.

physicskit.condensed.ginzburg_landau.ginzburg_landau_parameter(coherence_length, penetration_depth)[source]#

Ginzburg-Landau parameter \(\kappa = \lambda/\xi\).

The single dimensionless number that decides a superconductor’s response to a magnetic field: \(\kappa < 1/\sqrt2\) is Type I (the normal-superconducting interface has positive surface energy, and the field is excluded entirely below \(H_c\)); \(\kappa > 1/\sqrt2\) is Type II (negative surface energy favors flux penetrating as an Abrikosov vortex lattice between \(H_{c1}\) and \(H_{c2}\)).

Parameters:

coherence_length (float) – \(\xi\), from gl_coherence_length().
penetration_depth (float) – \(\lambda\), from gl_penetration_depth().

Return type:

Returns:

float

Examples

>>> round(ginzburg_landau_parameter(coherence_length=1.0, penetration_depth=1.0), 4)
1.0

physicskit.condensed.ginzburg_landau.gl_coherence_length(a, hbar=1.0, m=1.0)[source]#

Ginzburg-Landau coherence length \(\xi = \hbar/\sqrt{2m|a|}\).

The length scale over which the order parameter heals back to its bulk value after being suppressed at a boundary or a vortex core.

Parameters:

a (float) – Quadratic coefficient (only its magnitude matters).
hbar (float) – Reduced Planck constant.
m (float) – Effective mass of the condensate.

Return type:

Returns:

float

Examples

>>> gl_coherence_length(a=-0.5)
1.0

physicskit.condensed.ginzburg_landau.gl_equilibrium_order_parameter(a, b)[source]#

Equilibrium (uniform, field-free) order parameter magnitude \(|\psi_0|=\sqrt{-a/b}\).

Parameters:

a (float) – Quadratic coefficient.
b (float) – Quartic coefficient, b > 0.

Return type:

Returns:

float – \(\sqrt{-a/b}\) if a < 0 (ordered phase), else 0.0 (disordered phase, only the normal state minimizes \(f\)).

Examples

>>> gl_equilibrium_order_parameter(a=-2.0, b=2.0)
1.0
>>> gl_equilibrium_order_parameter(a=1.0, b=2.0)
0.0

physicskit.condensed.ginzburg_landau.gl_free_energy_density(psi, a, b, grad_psi=0.0, hbar=1.0, m=1.0)[source]#

Ginzburg-Landau free energy density \(f = a|\psi|^2 + \tfrac{b}{2}|\psi|^4 + \tfrac{\hbar^2}{2m}|\nabla\psi|^2\).

Parameters:

psi (complex or array_like) – Order parameter value(s).
a (float) – Quadratic coefficient. Changes sign at the transition (a < 0 in the ordered phase, a > 0 in the disordered phase).
b (float) – Quartic coefficient, b > 0 for stability.
grad_psi (complex or array_like, default=0.0) – Gradient \(\nabla\psi\), same shape as psi.
hbar (float) – Reduced Planck constant.
m (float) – Effective mass of the condensate.

Returns:

float or ndarray

Examples

>>> gl_free_energy_density(psi=0.0, a=-1.0, b=1.0)
0.0
>>> psi0 = gl_equilibrium_order_parameter(a=-1.0, b=1.0)
>>> round(gl_free_energy_density(psi0, a=-1.0, b=1.0), 6)
-0.5

physicskit.condensed.ginzburg_landau.gl_order_parameter_profile(x, xi)[source]#

Order parameter healing profile \(\psi(x)/\psi_0 = \tanh(x/(\sqrt2\,\xi))\).

The exact solution of the dimensionless Ginzburg-Landau equation \(\xi^2\psi'' = \psi^3-\psi\) for a condensate pinned to zero at a boundary (x = 0, e.g. a normal-superconducting interface) and recovering its bulk value far from it – the direct, textbook illustration of the coherence length as a healing length.

Parameters:

x (array_like) – Distance from the boundary.
xi (float) – Coherence length, from gl_coherence_length().

Return type:

Returns:

ndarray – \(\psi(x)/\psi_0\), in \([0, 1)\) for \(x \geq 0\).

Examples

>>> import numpy as np
>>> round(float(gl_order_parameter_profile(0.0, xi=1.0)), 8)
0.0
>>> bool(gl_order_parameter_profile(np.array([10.0]), xi=1.0)[0] > 0.999)
True

physicskit.condensed.ginzburg_landau.gl_penetration_depth(psi0, e=1.0, m=1.0)[source]#

London penetration depth \(\lambda = \sqrt{m/(e^2|\psi_0|^2)}\).

The length scale over which an external magnetic field is screened from a superconductor’s interior by the supercurrent it induces in the condensate of density \(|\psi_0|^2\).

Parameters:

psi0 (float) – Equilibrium order parameter magnitude (see gl_equilibrium_order_parameter()).
e (float) – Cooper-pair charge magnitude.
m (float) – Effective mass of the condensate.

Return type:

Returns:

float

Examples

>>> gl_penetration_depth(psi0=1.0)
1.0

Anderson localization: the 1958 discovery that disorder can halt diffusion entirely.

Philip Anderson showed that a quantum particle hopping on a lattice with sufficiently strong, random (quenched) onsite disorder does not diffuse at all: interference between all the scattering paths off the random potential can exponentially localize every eigenstate, at any disorder strength, in one and two dimensions. This was the first demonstration that disorder is not merely a small perturbative correction to a metal’s conductivity, but can drive a genuine transition – absent in three dimensions only above a critical disorder (the “mobility edge”) – from extended, conducting states to localized, insulating ones. It underlies the disorder-driven physics that makes real quantum Hall plateaus finite in width (physicskit.condensed.topology) and the transport (or lack of it) in every real, imperfect crystal.

Uses the same natural-unit convention (\(\hbar=e=1\)) as the rest of physicskit.condensed.

physicskit.condensed.anderson_localization.anderson_chain_hamiltonian(n_sites, disorder_strength, t=1.0, seed=None)[source]#

Real-space 1D Anderson model: a tight-binding chain with random onsite disorder.

\[H = -t\sum_i \left(c_i^\dagger c_{i+1} + \text{h.c.}\right) + \sum_i \epsilon_i\, c_i^\dagger c_i, \qquad \epsilon_i \sim \text{Uniform}\!\left(-\tfrac{W}{2}, \tfrac{W}{2}\right),\]

with open boundary conditions.

Parameters:

n_sites (int) – Number of lattice sites.
disorder_strength (float) – Disorder width \(W\). 0 recovers the clean chain.
t (float) – Nearest-neighbor hopping amplitude.
seed (int | None) – Seed for the random onsite disorder, for reproducibility.

Return type:

Returns:

ndarray, shape (n_sites, n_sites) – Real, symmetric Hamiltonian matrix.

Examples

>>> H = anderson_chain_hamiltonian(n_sites=4, disorder_strength=0.0)
>>> H
array([[ 0., -1.,  0.,  0.],
       [-1.,  0., -1.,  0.],
       [ 0., -1.,  0., -1.],
       [ 0.,  0., -1.,  0.]])

physicskit.condensed.anderson_localization.inverse_participation_ratio(psi)[source]#

Inverse participation ratio \(\text{IPR} = \sum_i|\psi_i|^4 / \left(\sum_i|\psi_i|^2\right)^2\).

A dimensionless measure of how many sites an eigenstate is spread over: \(\text{IPR}\sim 1/N\) for a state extended over all \(N\) sites, and \(\text{IPR} = O(1)\), independent of \(N\), for a state localized on a handful of sites – the diagnostic Anderson used to distinguish the two regimes.

Parameters:: psi (array_like) – Eigenvector components (need not be pre-normalized).
Return type:: float
Returns:: float

Examples

A state spread equally over N sites has IPR = 1/N:

>>> import numpy as np
>>> psi = np.ones(10) / np.sqrt(10)
>>> round(inverse_participation_ratio(psi), 6)
0.1
>>> psi_localized = np.zeros(10); psi_localized[0] = 1.0
>>> inverse_participation_ratio(psi_localized)
1.0

physicskit.condensed.anderson_localization.localization_length(psi, positions=None)[source]#

Estimate an eigenstate’s localization length from the decay of its envelope.

Fits \(\log|\psi_i|^2\) linearly against position on either side of the state’s peak and returns \(\xi = -2/\text{slope}\), averaged over both sides – the length scale over which a localized state’s probability density decays as \(|\psi(x)|^2 \sim e^{-2|x-x_0|/\xi}\). Returns inf for a state too extended (flat/non-monotonic envelope) for the fit to detect exponential decay.

Parameters:

psi (array_like) – Eigenvector components, indexed by site.
positions (array_like, optional) – Site coordinates. Defaults to integer indices 0, 1, ..., N-1.

Return type:

Returns:

float

Examples

An exactly exponential envelope recovers its localization length:

>>> import numpy as np
>>> x = np.arange(200)
>>> xi_true = 5.0
>>> psi = np.exp(-np.abs(x - 100) / xi_true)
>>> round(localization_length(psi), 4)
5.0

Correlated-electron models: Bogoliubov-de Gennes superconductivity and the Hubbard model.

Provides a mean-field BCS/Bogoliubov-de Gennes (BdG) solver for s-wave superconductivity, and small-cluster exact diagonalization (ED) of the 1D Fermi-Hubbard model.

physicskit.condensed.correlated.bdg_bcs_hamiltonian(k, mu=0.0, t=1.0, delta=0.5)[source]#

Mean-field Bogoliubov-de Gennes Hamiltonian for a 1D s-wave BCS superconductor.

Parameters:

k (float) – Reduced crystal momentum, periodic on \([0, 2\pi)\).
mu (float) – Chemical potential.
t (float) – Nearest-neighbor hopping amplitude (sets the normal-state band \(\xi(k) = -2t\cos k - \mu\)).
delta (complex) – s-wave (momentum-independent) pairing amplitude.

Return type:

Returns:

ndarray, shape (2, 2) – BdG Hamiltonian \(H(k) = \xi(k)\tau_z + \mathrm{Re}(\Delta)\tau_x - \mathrm{Im}(\Delta)\tau_y\) in the Nambu basis \((c_k, c_{-k}^\dagger)\).

See also

bdg_spectrum: Quasiparticle energies over a grid of k, exposing the gap.

Examples

>>> import numpy as np
>>> H = bdg_bcs_hamiltonian(k=np.pi / 2, mu=0.0, t=1.0, delta=0.5)
>>> np.round(np.linalg.eigvalsh(H), 8)
array([-0.5,  0.5])

physicskit.condensed.correlated.bdg_spectrum(mu=0.0, t=1.0, delta=0.5, n_k=200)[source]#

Quasiparticle (positive-energy) BdG spectrum over the 1D Brillouin zone.

Parameters:

mu (float) – Chemical potential.
t (float) – Nearest-neighbor hopping amplitude.
delta (complex) – s-wave pairing amplitude.
n_k (int) – Number of momentum points.

Return type:

Returns:

k_grid (ndarray, shape (n_k,)) – Momentum grid over \([0, 2\pi)\).
energies (ndarray, shape (n_k,)) – Positive quasiparticle branch \(E(k) = \sqrt{\xi(k)^2 + |\Delta|^2}\).

Examples

The minimum quasiparticle energy equals the pairing gap \(|\Delta|\):

>>> k, E = bdg_spectrum(mu=0.0, t=1.0, delta=0.5, n_k=400)
>>> round(float(E.min()), 3)
0.5

physicskit.condensed.correlated.hubbard_1d_exact_diagonalization(n_sites, n_up, n_dn, t=1.0, U=4.0, pbc=False, return_eigenvectors=False)[source]#

Exact diagonalization of the 1D Fermi-Hubbard model on a small cluster.

\[H = -t \sum_{\langle ij\rangle, \sigma} c_{i\sigma}^\dagger c_{j\sigma} + U \sum_i n_{i\uparrow} n_{i\downarrow}\]

Builds the Fock-space Hamiltonian in the fixed \((n_\uparrow, n_\downarrow)\) sector via full (dense) diagonalization; suitable for small clusters (\(n_{sites} \lesssim 8\)).

Parameters:

n_sites (int) – Number of lattice sites.
n_up (int) – Number of spin-up and spin-down electrons.
n_dn (int) – Number of spin-up and spin-down electrons.
t (float) – Nearest-neighbor hopping amplitude.
U (float) – Onsite Coulomb repulsion.
pbc (bool) – Periodic (True) or open (False) boundary conditions.
return_eigenvectors (bool) – If True, also return the ground-state vector and the occupation- number bases it is expanded in, e.g. for hubbard_spin_correlations().

Return type:

dict

Returns:

dict – {"ground_state_energy": float, "eigenvalues": ndarray, "dimension": int}, plus {"ground_state_vector": ndarray, "states_up": list, "states_dn": list} when return_eigenvectors=True. ground_state_vector[a * len(states_dn) + b] is the amplitude of the product basis state (states_up[a], states_dn[b]), each an integer bitmask with bit i set if site i is occupied.

Notes

In the large-\(U\) limit at half filling, double occupancy is suppressed and the ground state approaches the Mott-insulating, singly-occupied regime characteristic of spin-charge separation in 1D.

Examples

The 2-site Hubbard dimer at half filling has the exact ground-state energy \(E_0 = \tfrac{1}{2}\big(U - \sqrt{U^2 + 16t^2}\big)\):

>>> import numpy as np
>>> result = hubbard_1d_exact_diagonalization(n_sites=2, n_up=1, n_dn=1, t=1.0, U=4.0)
>>> E_exact = 0.5 * (4.0 - np.sqrt(4.0**2 + 16 * 1.0**2))
>>> bool(round(result["ground_state_energy"], 8) == round(E_exact, 8))
True

physicskit.condensed.correlated.hubbard_spin_correlations(n_sites, ground_state_vector, states_up, states_dn)[source]#

Equal-time z-spin correlations \(\langle S_i^z S_j^z\rangle\) from a Hubbard ground state.

\(S_i^z = \tfrac{1}{2}(n_{i\uparrow} - n_{i\downarrow})\) is diagonal in the occupation-number basis hubbard_1d_exact_diagonalization() builds its Hamiltonian in, so \(\langle S_i^z S_j^z\rangle\) reduces to a weighted sum over basis-state probabilities, \(\sum_\alpha |\psi_\alpha|^2\, s_i(\alpha)\, s_j(\alpha)\), with no off-diagonal matrix elements to track.

Parameters:

n_sites (int) – Number of lattice sites.
ground_state_vector (ndarray) – Ground-state amplitudes, as returned by hubbard_1d_exact_diagonalization() with return_eigenvectors=True.
states_up (list) – Occupation-number bases for each spin species, as returned alongside ground_state_vector.
states_dn (list) – Occupation-number bases for each spin species, as returned alongside ground_state_vector.

Return type:

Returns:

ndarray, shape (n_sites,) – correlations[r] is \(\langle S_i^z S_j^z\rangle\) averaged over every site pair with \(|i-j|=r\); correlations[0] is the onsite \(\langle (S_i^z)^2\rangle\), suppressed toward 0 by double/empty occupancy and toward its \(1/4\) ceiling by strong onsite repulsion.

Notes

A Mott-insulating antiferromagnet’s short-range Neel order shows up as correlations alternating in sign with r – negative (antialigned) at odd separations, positive (aligned) at even ones – the same short-range magnetic correlations believed to survive doping into the cuprate superconductors’ metallic phase.

Examples

At half filling and strong coupling, nearest-neighbor spins are antialigned:

>>> import numpy as np
>>> result = hubbard_1d_exact_diagonalization(n_sites=6, n_up=3, n_dn=3, t=1.0, U=8.0, pbc=True, return_eigenvectors=True)
>>> corr = hubbard_spin_correlations(6, result["ground_state_vector"], result["states_up"], result["states_dn"])
>>> bool(corr[1] < 0)
True

Canonical condensed-matter lattice models: SSH, graphene, Haldane, Kane-Mele, BHZ, Kitaev chain.

All 2D Bloch Hamiltonians in this module are functions of the reduced crystal momentum (k1, k2) described in physicskit.condensed.tight_binding (periodic gauge, each component periodic on \([0, 2\pi)\)), which is the convention expected by physicskit.condensed.topology.compute_chern_number().

physicskit.condensed.models.bhz_hamiltonian(kx, ky, A=1.0, B=1.0, M=1.0, D=0.0)[source]#

Bloch Hamiltonian of the Bernevig-Hughes-Zhang (BHZ) model for HgTe quantum wells.

A minimal 4-band \(\mathbb{Z}_2\) topological insulator model, block-diagonal in a time-reversed pair of 2x2 Dirac-like blocks. Basis order is (E up, H up, E down, H down).

Parameters:

kx (float) – Reduced crystal momenta, each periodic on \([0, 2\pi)\).
ky (float) – Reduced crystal momenta, each periodic on \([0, 2\pi)\).
A (float) – Model parameters controlling the Dirac velocity and quadratic band curvature.
B (float) – Model parameters controlling the Dirac velocity and quadratic band curvature.
M (float) – Band inversion mass. With \(d_z = M - 2B(2-\cos k_x-\cos k_y)\) the model is topological (band-inverted, spin Chern number \(\pm1\)) for \(0 < M/B < 8\) and trivial otherwise (the gap closes at \(\Gamma\), \(X\)/\(Y\), and \(M\) for \(M/B = 0, 4, 8\)).
D (float) – Particle-hole asymmetry parameter, entering as \(\epsilon(k) = -2D(2-\cos k_x-\cos k_y)\), the lattice regularization of the continuum \(-Dk^2\) (Qi, Hughes & Zhang 2008), matching \(B\)’s.

Return type:

Returns:

ndarray, shape (4, 4)

Examples

>>> import numpy as np
>>> H = bhz_hamiltonian(0.1, -0.2, M=1.0, B=1.0)
>>> np.allclose(H, H.conj().T)
True
>>> np.allclose(H[:2, :2], H[2:, 2:].conj())
True

physicskit.condensed.models.bhz_ribbon_hamiltonian(kx, n_cells, A=1.0, B=1.0, M=1.0, D=0.0)[source]#

Real-space BHZ ribbon: periodic along x, open (finite) along y.

Truncating the y direction exposes the pair of helical, spin-locked edge states that make the BHZ model a quantum spin Hall insulator – the effect Konig et al. (2007) measured directly in HgTe/CdTe quantum wells, by two-terminal conductance quantized at \(2e^2/h\).

Parameters:

kx (float) – Reduced crystal momentum along the periodic (x) direction.
n_cells (int) – Number of unit cells stacked along the open (y) direction.
A (float) – Model parameters controlling the Dirac velocity and quadratic band curvature (see bhz_hamiltonian()).
B (float) – Model parameters controlling the Dirac velocity and quadratic band curvature (see bhz_hamiltonian()).
M (float) – Band inversion mass. Topological (edge-state-carrying) for \(0 < M/B < 8\).
D (float) – Particle-hole asymmetry parameter.

Return type:

Returns:

ndarray, shape (4 * n_cells, 4 * n_cells) – Hermitian, open-boundary ribbon Hamiltonian, block-diagonal in the (E up, H up) / (E down, H down) time-reversed sectors.

Examples

In the topological regime (\(0 < M/B < 8\)), the ribbon has a pair of near-zero-energy states crossing at \(k_x = 0\) – the helical edge modes – absent in the trivial regime (e.g. \(M/B < 0\)):

>>> import numpy as np
>>> spectrum = np.linalg.eigvalsh(bhz_ribbon_hamiltonian(kx=0.0, n_cells=40, M=1.0, B=1.0))
>>> bool(np.any(np.abs(spectrum) < 1e-6))
True
>>> spectrum_trivial = np.linalg.eigvalsh(bhz_ribbon_hamiltonian(kx=0.0, n_cells=40, M=-1.0, B=1.0))
>>> bool(np.any(np.abs(spectrum_trivial) < 1e-6))
False

physicskit.condensed.models.graphene_hamiltonian(k1, k2, t=1.0)[source]#

Bloch Hamiltonian of nearest-neighbor graphene (honeycomb lattice).

Parameters:

k1 (float) – Reduced crystal momenta conjugate to the two honeycomb primitive vectors, each periodic on \([0, 2\pi)\).
k2 (float) – Reduced crystal momenta conjugate to the two honeycomb primitive vectors, each periodic on \([0, 2\pi)\).
t (float) – Nearest-neighbor hopping amplitude.

Return type:

Returns:

ndarray, shape (2, 2) – Bloch Hamiltonian on the A/B sublattice basis.

Notes

The Dirac points sit at \((k_1, k_2) = (2\pi/3, 4\pi/3)\) and its time-reversed partner \((4\pi/3, 2\pi/3)\), where the gap closes and the dispersion is linear (massless Dirac cone), as verified in the Examples.

Examples

>>> import numpy as np
>>> K = np.array([2 * np.pi / 3, 4 * np.pi / 3])
>>> np.round(np.linalg.eigvalsh(graphene_hamiltonian(*K)), 8)
array([-0.,  0.])
>>> eigs = np.linalg.eigvalsh(graphene_hamiltonian(*(K + [1e-4, 0])))
>>> round(float(eigs[1] / 1e-4), 4)
1.0

physicskit.condensed.models.graphene_lattice_hamiltonian(t=1.0)[source]#

Nearest-neighbor honeycomb (graphene) tight-binding Hamiltonian.

Parameters:: t (float) – Nearest-neighbor hopping amplitude.
Return type:: Hamiltonian
Returns:: Hamiltonian

See also

graphene_hamiltonian: Closed-form 2x2 Bloch Hamiltonian equivalent to this builder.

physicskit.condensed.models.haldane_lattice_hamiltonian(t=1.0, t2=0.2, phi=1.5707963267948966, M=0.0)[source]#

Haldane model (quantum anomalous Hall) as a Hamiltonian.

Parameters:

t (float) – Nearest-neighbor hopping amplitude.
t2 (float) – Next-nearest-neighbor hopping magnitude.
phi (float) – Next-nearest-neighbor hopping phase, breaking time-reversal symmetry.
M (float) – Sublattice (Semenoff) mass, breaking inversion symmetry.

Return type:

Returns:

Hamiltonian

See also

haldane_model: Closed-form 2x2 Bloch Hamiltonian equivalent to this builder.

physicskit.condensed.models.haldane_model(kx, ky, t=1.0, t2=0.2, phi=1.5707963267948966, M=0.0)[source]#

Bloch Hamiltonian of the Haldane model on the honeycomb lattice.

The Haldane model realizes the quantum anomalous Hall effect: a Chern insulator with zero net magnetic flux per unit cell, arising from a complex next-nearest-neighbor hopping \(t_2 e^{i\phi}\) that breaks time-reversal symmetry while preserving the lattice translational symmetry.

Parameters:

kx (float) – Reduced crystal momenta (see physicskit.condensed.tight_binding), each periodic on \([0, 2\pi)\).
ky (float) – Reduced crystal momenta (see physicskit.condensed.tight_binding), each periodic on \([0, 2\pi)\).
t (float) – Nearest-neighbor hopping amplitude.
t2 (float) – Next-nearest-neighbor hopping magnitude.
phi (float) – Next-nearest-neighbor hopping phase.
M (float) – Sublattice mass. The model is a Chern insulator (\(C = \mathrm{sgn}(\sin\phi)\)) for \(|M| < 3\sqrt{3}\,t_2|\sin\phi|\), and a trivial insulator otherwise.

Return type:

Returns:

ndarray, shape (2, 2)

See also

physicskit.condensed.topology.compute_chern_number: Computes the Chern number of this model.

Examples

>>> import numpy as np
>>> from physicskit.condensed.topology import compute_chern_number
>>> H_func = lambda k1, k2: haldane_model(k1, k2, t=1.0, t2=0.2, phi=np.pi / 2, M=0.0)
>>> compute_chern_number(H_func, grid_size=30)
[1, -1]
>>> H_trivial = lambda k1, k2: haldane_model(k1, k2, t=1.0, t2=0.2, phi=np.pi / 2, M=2.0)
>>> compute_chern_number(H_trivial, grid_size=30)
[0, 0]

physicskit.condensed.models.harper_hofstadter_hamiltonian(k1, k2, p, q, t=1.0)[source]#

Bloch Hamiltonian of the square-lattice Harper-Hofstadter model at flux p/q.

A charged particle on a square lattice threaded by a uniform magnetic flux \(p/q\) (in units of the flux quantum) per plaquette, in Landau gauge \(\mathbf{A} = (0, Bx)\). Translational symmetry along x survives only in steps of q lattice constants, so the q inequivalent sublattice sites m = 0, ..., q-1 within one magnetic unit cell become an internal band index, coupled by

\[H_{mm}(k_1, k_2) = -2t\cos\!\left(k_2 + 2\pi \frac{p}{q} m\right), \qquad H_{m, m+1} = -t,\]

plus the boundary hopping \(H_{q-1, 0} = -t\,e^{ik_1}\) that closes the magnetic unit cell, with \(k_1\) the reduced momentum conjugate to translation by one magnetic cell (period \(q\) sites) and \(k_2\) the ordinary reduced momentum along the unbroken direction. This is the lattice (Bloch) route to the same physics as physicskit.condensed.landau_levels, and, fed one band at a time into compute_chern_number(), the standard numerical verification of the TKNN integer quantum Hall formula: summing the Chern numbers of the lowest \(r\) (non-touching) bands gives the exactly quantized Hall conductance \(\sigma_{xy} = C\,e^2/h\) at the filling between band \(r\) and \(r+1\).

Parameters:

k1 (float) – Reduced crystal momenta, each periodic on \([0, 2\pi)\). k1 is conjugate to the enlarged (q-site) magnetic unit cell.
k2 (float) – Reduced crystal momenta, each periodic on \([0, 2\pi)\). k1 is conjugate to the enlarged (q-site) magnetic unit cell.
p (int) – Flux per plaquette \(p/q\) (in lowest terms), \(0 < p < q\). Odd q avoids the exact band touchings that occur for some even q, where individual-band Chern numbers become ill-defined.
q (int) – Flux per plaquette \(p/q\) (in lowest terms), \(0 < p < q\). Odd q avoids the exact band touchings that occur for some even q, where individual-band Chern numbers become ill-defined.
t (float) – Nearest-neighbor hopping amplitude.

Return type:

Returns:

ndarray, shape (q, q) – Hermitian Bloch Hamiltonian.

See also

physicskit.condensed.topology.compute_chern_number: Chern number of each resulting band.

Examples

At flux 1/3 the three bands carry Chern numbers -1, 2, -1 (summing to zero, as any complete set of bands of a lattice Hamiltonian must), so the Hall conductance is quantized to \(-1\) and \(+1\) (in units of \(e^2/h\)) at the two gaps:

>>> import numpy as np
>>> from physicskit.condensed.topology import compute_chern_number
>>> H_func = lambda k1, k2: harper_hofstadter_hamiltonian(k1, k2, p=1, q=3)
>>> chern_numbers = compute_chern_number(H_func, grid_size=30)
>>> chern_numbers
[-1, 2, -1]
>>> sum(chern_numbers)
0

physicskit.condensed.models.kane_mele_hamiltonian(kx, ky, t=1.0, lambda_so=0.06, lambda_v=0.0, lambda_r=0.0)[source]#

Bloch Hamiltonian of the Kane-Mele quantum spin Hall model.

Two time-reversed Haldane copies (intrinsic spin-orbit coupling with opposite Chern-number-generating phase for each spin), realizing a \(\mathbb{Z}_2\) topological insulator with helical edge states. Basis order is (A up, B up, A down, B down).

Parameters:

kx (float) – Reduced crystal momenta, each periodic on \([0, 2\pi)\).
ky (float) – Reduced crystal momenta, each periodic on \([0, 2\pi)\).
t (float) – Nearest-neighbor hopping amplitude.
lambda_so (float) – Intrinsic spin-orbit coupling strength.
lambda_v (float) – Sublattice (staggered) potential.
lambda_r (float) – Rashba spin-orbit coupling strength. A nonzero value mixes the spin blocks (breaking \(s_z\) conservation) but preserves overall time-reversal symmetry.

Return type:

Returns:

ndarray, shape (4, 4)

Notes

When lambda_r == 0, \(s_z\) is conserved and the Hamiltonian is block-diagonal in spin; each block is a Haldane model with \(t_2 = \lambda_{so}\), \(\phi = \pm\pi/2\). In that case the \(\mathbb{Z}_2\) invariant equals the spin-up Chern number modulo 2 (see physicskit.condensed.topology.z2_invariant()).

Examples

>>> import numpy as np
>>> H = kane_mele_hamiltonian(0.3, 0.7, lambda_so=0.06)
>>> np.allclose(H, H.conj().T)
True
>>> up = H[:2, :2]
>>> down_at_minus_k = kane_mele_hamiltonian(-0.3, -0.7, lambda_so=0.06)[2:, 2:]
>>> np.allclose(down_at_minus_k, up.conj())  # time reversal: H_dn(-k) = H_up(k)*
True

physicskit.condensed.models.kitaev_chain_bdg_real_space(n_sites, mu=0.0, t=1.0, delta=1.0)[source]#

Real-space open-boundary BdG Hamiltonian of the Kitaev chain.

Parameters:

n_sites (int) – Number of lattice sites in the finite, open chain.
mu (float) – Chemical potential.
t (float) – Nearest-neighbor hopping amplitude.
delta (float) – p-wave pairing amplitude.

Return type:

Returns:

ndarray, shape (2*n_sites, 2*n_sites) – BdG Hamiltonian in the Nambu basis \((c_1, \dots, c_N, c_1^\dagger, \dots, c_N^\dagger)\).

Examples

In the topological phase (\(|\mu| < 2t\)), diagonalizing this matrix yields a pair of near-zero-energy eigenvalues (Majorana end modes), exponentially localized at the two ends of the chain:

>>> import numpy as np
>>> H = kitaev_chain_bdg_real_space(n_sites=60, mu=0.0, t=1.0, delta=1.0)
>>> eigs = np.linalg.eigvalsh(H)
>>> bool(np.abs(eigs[np.argmin(np.abs(eigs))]) < 1e-6)
True

physicskit.condensed.models.kitaev_chain_hamiltonian(k, mu=0.0, t=1.0, delta=1.0)[source]#

Bogoliubov-de Gennes Bloch Hamiltonian of the Kitaev p-wave superconducting chain.

Parameters:

k (float) – Reduced crystal momentum, periodic on \([0, 2\pi)\).
mu (float) – Chemical potential.
t (float) – Nearest-neighbor hopping amplitude.
delta (float) – p-wave pairing amplitude.

Return type:

Returns:

ndarray, shape (2, 2) – BdG Hamiltonian \(H(k) = \xi(k)\tau_z + \Delta(k)\tau_y\) in the Nambu basis, with \(\xi(k) = -2t\cos k - \mu\) and \(\Delta(k) = 2\Delta\sin k\).

Notes

The chain is a topological superconductor, hosting unpaired Majorana zero modes at the ends of an open chain, for \(|\mu| < 2t\); it is topologically trivial for \(|\mu| > 2t\).

Examples

>>> import numpy as np
>>> H = kitaev_chain_hamiltonian(k=np.pi / 2, mu=0.0, t=1.0, delta=1.0)
>>> np.round(np.linalg.eigvalsh(H), 8)
array([-2.,  2.])

physicskit.condensed.models.ssh_hamiltonian(k, v=1.0, w=1.0)[source]#

Bloch Hamiltonian of the Su-Schrieffer-Heeger (SSH) dimerized chain.

Parameters:

k (float) – Reduced crystal momentum, periodic on \([0, 2\pi)\).
v (float) – Intracell hopping amplitude.
w (float) – Intercell hopping amplitude.

Return type:

Returns:

ndarray, shape (2, 2) – Bloch Hamiltonian \(H(k) = \begin{pmatrix}0 & v+we^{-ik}\\ v+we^{ik} & 0\end{pmatrix}\).

Notes

The chain is topologically nontrivial (hosts protected zero-energy edge states under open boundary conditions) when \(v < w\), and trivial when \(v > w\). See Also ssh_lattice_hamiltonian() for the real-space builder used to expose those edge states.

Examples

>>> import numpy as np
>>> H = ssh_hamiltonian(k=0.0, v=0.5, w=1.0)
>>> np.round(np.linalg.eigvalsh(H), 8)
array([-1.5,  1.5])

physicskit.condensed.models.ssh_lattice_hamiltonian(v=1.0, w=1.0)[source]#

Build the SSH chain as a Hamiltonian.

Parameters:

v (float) – Intracell hopping amplitude (A-B bond within a unit cell).
w (float) – Intercell hopping amplitude (B-A bond across the cell boundary).

Return type:

Returns:

Hamiltonian – A two-orbital-per-cell chain, suitable for build_ribbon() to expose edge zero modes.

See also

ssh_hamiltonian: Closed-form 2x2 Bloch Hamiltonian equivalent to this builder.

Examples

>>> import numpy as np
>>> H = ssh_lattice_hamiltonian(v=0.5, w=1.0)
>>> np.round(H.bands([0.0]), 8)
array([-1.5,  1.5])

The 3D strong topological insulator: from theory (2007) to Bi2Se3/Bi2Te3 (2008-2009).

Fu, Kane, and Mele extended the 2D \(\mathbb{Z}_2\) classification of physicskit.condensed.models (kane_mele_hamiltonian, bhz_hamiltonian) to three dimensions, predicting a “strong” topological insulator: a bulk band insulator whose every surface, regardless of orientation, hosts a single, gapless, spin-momentum-locked Dirac cone protected by time-reversal symmetry alone. Hasan and Cava’s groups and collaborators confirmed this directly in 2008-2009 by ARPES on Bi₂Se₃ and Bi₂Te₃, resolving exactly one Dirac cone at the surface Brillouin zone center – the 3D generalization of the helical edge states Konig et al. had just observed at the 1D edge of a 2D quantum spin Hall system.

This module implements the minimal 4-band lattice model with this physics (Qi & Zhang, Rev. Mod. Phys. 83, 1057 (2011)):

\[H(\mathbf{k}) = \sin k_x\,\Gamma_1 + \sin k_y\,\Gamma_2 + \sin k_z\,\Gamma_3 + \left(m + t\sum_{i=x,y,z}\cos k_i\right)\Gamma_0,\]

built from mutually anticommuting 4x4 Dirac matrices \(\Gamma_0=\tau_z\otimes\sigma_0\), \(\Gamma_{1,2,3}=\tau_x\otimes\sigma_{x,y,z}\) (\(\tau\) an orbital pseudospin, \(\sigma\) the physical electron spin) – a cubic-lattice generalization of the BHZ model whose bulk gap closes, and topological character changes, at every time-reversal-invariant momentum where the mass term vanishes, \(m/t \in \{-3, -1, 1, 3\}\). Of the resulting windows, \(1 < |m/t| < 3\) is a strong topological insulator (an open slab hosts a single, gapless surface Dirac cone at \(\bar\Gamma=(0,0)\) on each of its two open surfaces); \(|m/t| < 1\) and \(|m/t| > 3\) are trivial.

Uses the same natural-unit convention (\(\hbar=e=1\)) as the rest of physicskit.condensed.

physicskit.condensed.topological_insulator_3d.surface_dirac_hamiltonian(kx, ky, v_f=1.0)[source]#

Low-energy effective Hamiltonian of a single topological-insulator surface Dirac cone.

\[H_{\text{surf}}(\mathbf{k}) = \hbar v_F\left(k_x\sigma_y - k_y\sigma_x\right),\]

the massless, spin-momentum-locked Dirac fermion measured directly by ARPES on the Bi₂Se₃ and Bi₂Te₃ surfaces in 2008-2009: spin polarization locked perpendicular to momentum, with no Kramers-degenerate partner at the same energy and momentum (backscattering off nonmagnetic disorder is forbidden).

Parameters:

kx (float) – Momentum measured from the surface Dirac point.
ky (float) – Momentum measured from the surface Dirac point.
v_f (float) – Surface Fermi velocity.

Return type:

Returns:

ndarray, shape (2, 2) – Hermitian Hamiltonian with linear (Dirac) spectrum \(E_\pm(\mathbf{k}) = \pm v_F|\mathbf{k}|\).

Examples

>>> import numpy as np
>>> H = surface_dirac_hamiltonian(0.3, -0.4, v_f=2.0)
>>> np.round(np.linalg.eigvalsh(H), 8)
array([-1.,  1.])

physicskit.condensed.topological_insulator_3d.topological_insulator_3d_hamiltonian(kx, ky, kz, m, t=1.0)[source]#

Bloch Hamiltonian of the minimal cubic-lattice 3D topological insulator model.

Parameters:

kx (float) – Reduced crystal momenta, each periodic on \([0, 2\pi)\).
ky (float) – Reduced crystal momenta, each periodic on \([0, 2\pi)\).
kz (float) – Reduced crystal momenta, each periodic on \([0, 2\pi)\).
m (float) – Band mass. The model is a strong topological insulator for \(1 < |m/t| < 3\), and trivial for \(|m/t| < 1\) or \(|m/t| > 3\).
t (float) – Hopping amplitude setting the bulk bandwidth.

Return type:

Returns:

ndarray, shape (4, 4) – Hermitian Bloch Hamiltonian.

Examples

>>> import numpy as np
>>> H = topological_insulator_3d_hamiltonian(0.3, -0.7, 0.5, m=1.0)
>>> np.allclose(H, H.conj().T)
True
>>> H0 = topological_insulator_3d_hamiltonian(0.0, 0.0, 0.0, m=1.0, t=1.0)
>>> np.round(np.linalg.eigvalsh(H0), 8)
array([-4., -4.,  4.,  4.])

physicskit.condensed.topological_insulator_3d.topological_insulator_3d_slab_hamiltonian(kx, ky, n_layers, m, t=1.0)[source]#

Slab Hamiltonian: periodic in x, y, open (finite) along z.

Truncating the z direction exposes the top and bottom (001) surfaces. Diagonalizing at each (kx, ky) and looking for mid-gap states localized at the outer layers reveals the surface Dirac cone directly, with no separate topological-invariant calculation needed.

Parameters:

kx (float) – Reduced crystal momenta, periodic on \([0, 2\pi)\).
ky (float) – Reduced crystal momenta, periodic on \([0, 2\pi)\).
n_layers (int) – Number of atomic layers stacked along z.
m (float) – Band mass (see topological_insulator_3d_hamiltonian()).
t (float) – Hopping amplitude.

Return type: