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.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.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.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.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.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.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.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.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.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.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.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:

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.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)
>>> 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)
>>> 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: