Density of States

The density of states counts independent one-particle states per unit energy. For crystalline bands it compresses the full dispersion $\varepsilon_n(\mathbf k)$ into an energy distribution, retaining how much state space lies near each energy while discarding where those states occur in momentum space.

That compression is powerful but incomplete. The density of states controls filling integrals, low-temperature thermodynamics, phase space for transitions, and many spectroscopic line shapes. It does not by itself determine velocities, scattering times, optical matrix elements, topology, or whether states are extended.

Required background. Brillouin Zones supplies the full-zone measure and state-counting domain; Density of States: First Encounter supplies the continuum counting and delta-function definition generalized here to bands.

Counting electronic states per energy interval

For a finite one-particle spectrum,

$$D(E) = \sum_{\lambda} \delta(E-E_{\lambda}).$$

$D(E)$ is a distribution with units of inverse energy. The label $\lambda$ must enumerate every independent state exactly once, including any spin, valley, orbital, or other degeneracy intended by the convention.

The integrated density of states is

$$\mathcal N(E) = \int_{-\infty}^{E} dE'\, D(E').$$

For a finite system it is a staircase that jumps by the degeneracy of each level. Distributionally,

$$D(E) = \frac{d\mathcal N}{dE}.$$

A smooth curve appears only after a thermodynamic limit, a coarse graining, or an explicitly stated broadening.

Total, per-cell and per-volume densities

Three normalizations are common and must not be mixed.

Total density

For $N_c$ cells,

$$D_{\mathrm{tot}}(E) = \sum_{n,\mathbf k} \delta \left[ E-\varepsilon_n(\mathbf k) \right].$$

It scales extensively with system size.

Density per primitive cell

Define

$$\rho_c(E) = \frac{D_{\mathrm{tot}}(E)}{N_c}.$$

For a model with $N_b$ bands, where $n$ already includes all counted internal states,

$$\int_{-\infty}^{\infty} dE\, \rho_c(E) = N_b.$$

Thus $\rho_c$ has units of states per cell per energy.

Density per physical volume

With $V=N_c\Omega_c$,

$$\rho_V(E) = \frac{D_{\mathrm{tot}}(E)}{V} = \frac{\rho_c(E)}{\Omega_c}.$$

Its units are states per physical volume per energy. In two-dimensional materials one often quotes density per area; in one dimension, per length.

If a displayed band structure omits spin and every band is spin degenerate, multiply by $g_s=2$ only once. When spin–orbit coupling splits the bands and spin is already included in $n$, no additional factor is allowed.

Counting States in the Brillouin Zone

For a finite periodic crystal,

$$\rho_c(E) = \frac{1}{N_c} \sum_{n,\mathbf k} \delta \left[ E-\varepsilon_n(\mathbf k) \right].$$

The thermodynamic rule

$$\frac{1}{N_c} \sum_{\mathbf k} \longrightarrow \frac{\Omega_c}{(2\pi)^d} \int_{\mathrm{BZ}} d^dk$$

gives

$$\rho_c(E) = \frac{\Omega_c}{(2\pi)^d} \sum_n \int_{\mathrm{BZ}} d^dk\, \delta \left[ E-\varepsilon_n(\mathbf k) \right].$$

Equivalently,

$$\rho_V(E) = \frac{1}{(2\pi)^d} \sum_n \int_{\mathrm{BZ}} d^dk\, \delta \left[ E-\varepsilon_n(\mathbf k) \right].$$

The Brillouin-zone volume

$$\operatorname{Vol}(\mathrm{BZ}) = \frac{(2\pi)^d}{\Omega_c}$$

immediately proves the band-count sum rule.

Integrated form

The number of states per cell below $E$ is

$$\mathcal N_c(E) = \frac{\Omega_c}{(2\pi)^d} \sum_n \int_{\mathrm{BZ}} d^dk\, \Theta \left[ E-\varepsilon_n(\mathbf k) \right].$$

This is the Brillouin-zone volume occupied by all sublevel sets $\varepsilon_n(\mathbf k)\leq E$, measured in units of one state per band per cell.

Constant-Energy Geometry

Away from critical points, the delta-function integral can be converted into an integral over the constant-energy surface

$$\mathcal S_n(E) = \left\{ \mathbf k: \varepsilon_n(\mathbf k)=E \right\}.$$

The coarea formula gives

$$\rho_V(E) = \frac{1}{(2\pi)^d} \sum_n \int_{\mathcal S_n(E)} \frac{ dS_{\mathbf k} }{ \left| \boldsymbol\nabla_{\mathbf k} \varepsilon_n(\mathbf k) \right| }.$$

Using

$$\mathbf v_n(\mathbf k) = \frac{1}{\hbar} \boldsymbol\nabla_{\mathbf k} \varepsilon_n(\mathbf k),$$

this becomes

$$\rho_V(E) = \frac{1}{(2\pi)^d} \sum_n \int_{\mathcal S_n(E)} \frac{ dS_{\mathbf k} }{ \hbar|\mathbf v_n(\mathbf k)| }.$$

Large constant-energy surfaces and small group velocities both enhance the density of states. The formula is not regular at a critical point where $\boldsymbol\nabla_{\mathbf k}\varepsilon_n=\mathbf0$; those points generate van Hove nonanalyticities.

In one dimension, a constant-energy surface is a collection of roots $k_i$:

$$\rho_c(E) = \frac{a}{2\pi} \sum_{n,k_i} \frac{1}{ \left| d\varepsilon_n/dk \right|_{k_i} },$$

for a primitive length $a$. Each root must be counted, including symmetry-related left- and right-moving states.

Parabolic Band Edges

Consider an isotropic extremum

$$\varepsilon(\mathbf k) = E_0 + \frac{\hbar^2k^2}{2m^\ast}, \qquad m^\ast>0,$$

with total internal degeneracy $g$. Define

$$\xi = E-E_0.$$

The number of states per volume below $E$ is

$$\frac{\mathcal N(E)}{V} = g \frac{S_{d-1}}{ d(2\pi)^d } \left( \frac{2m^\ast\xi}{\hbar^2} \right)^{d/2} \Theta(\xi),$$

where

$$S_{d-1} = \frac{ 2\pi^{d/2} }{ \Gamma(d/2) }$$

is the area of the unit $(d-1)$-sphere. Differentiation gives

$$\rho_d(E) = g \frac{S_{d-1}}{(2\pi)^d} \frac{m^\ast}{\hbar^2} \left( \frac{2m^\ast\xi}{\hbar^2} \right)^{(d-2)/2} \Theta(\xi).$$

The threshold law is therefore

$$\rho_d(E) \propto \xi^{d/2-1} \Theta(\xi).$$

One dimension

$$\rho_1(E) = \frac{g}{\pi\hbar} \sqrt{ \frac{m^\ast}{2\xi} } \Theta(\xi).$$

The inverse-square-root divergence reflects the vanishing velocity at the band edge.

Two dimensions

$$\rho_2(E) = \frac{g m^\ast}{2\pi\hbar^2} \Theta(\xi).$$

An ideal parabolic two-dimensional band has a constant density of states with a step at its edge.

Three dimensions

$$\rho_3(E) = \frac{g}{4\pi^2} \left( \frac{2m^\ast}{\hbar^2} \right)^{3/2} \sqrt{\xi} \, \Theta(\xi).$$

The density rises continuously with a square-root onset.

Dimensionality fixes the ideal threshold law near a nondegenerate parabolic minimum: divergence in $d=1$, a step in $d=2$, and a square-root onset in $d=3$. A nondegenerate saddle in two dimensions instead produces a logarithmic van Hove singularity. Disorder, interaction self-energies, and numerical broadening can round these ideal DOS features. Thermal occupations and thermally convolved measurements smear their observed signatures; temperature alone does not change the DOS of a fixed noninteracting Hamiltonian.

Van Hove Singularities

A van Hove critical point satisfies

$$\boldsymbol\nabla_{\mathbf k} \varepsilon_n(\mathbf k_c) = \mathbf0.$$

Near an isolated nondegenerate critical point,

$$\varepsilon_n(\mathbf k_c+\mathbf q) = E_c + \frac12 \mathbf q^{\mathsf T} \mathcal H \mathbf q + O(q^3),$$

where $\mathcal H$ is the Hessian matrix of second derivatives. The numbers of positive and negative Hessian eigenvalues distinguish minima, maxima, and saddles.

For ordinary Morse critical points:

• one-dimensional extrema produce inverse-square-root divergences;
• two-dimensional extrema produce step discontinuities;
• two-dimensional saddles produce logarithmic divergences;
• three-dimensional extrema produce square-root onsets;
• three-dimensional saddles generally produce finite square-root cusps with a singular derivative.

At a two-dimensional saddle,

$$\varepsilon-E_c \approx \alpha q_x^2 - \beta q_y^2, \qquad \alpha,\beta>0.$$

With a finite momentum cutoff $\Lambda$ set by the range of the local expansion,

$$\rho_{\mathrm{sing}}(E) \sim C \ln \left( \frac{\Lambda_E}{ |E-E_c| } \right),$$

where $C>0$ depends on the curvature and normalization. The logarithm is universal; its additive constant and ultraviolet scale are not.

Beyond ordinary critical points

A flat band, a line of stationary points, or a critical point with a singular Hessian can produce stronger or nonstandard behavior. Such higher-order van Hove singularities require their own local expansion rather than the quadratic Morse classification.

Symmetry can enforce several equivalent critical points at the same energy, multiplying the singular contribution. Conversely, spin–orbit coupling, strain, or symmetry breaking can split one feature into several nearby peaks.

Example: One-Dimensional Tight-Binding Band

For

$$\varepsilon(k) = \epsilon_0 - 2t\cos(ka),$$

the spinless density per cell is

$$\rho_c(E) = \frac{ \Theta \left( 2|t|-|E-\epsilon_0| \right) }{ \pi \sqrt{ 4t^2-(E-\epsilon_0)^2 } }.$$

It diverges at both band edges,

$$E = \epsilon_0\pm2|t|,$$

and obeys

$$\int dE\, \rho_c(E) = 1.$$

The singularities do not represent extra states. They describe the compression of a fixed number of $k$ states into a narrow energy interval where $d\varepsilon/dk$ vanishes.

Example: Square-Lattice Saddle

For nearest-neighbor hopping,

$$\varepsilon(\mathbf k) = -2t \left[ \cos(k_xa) + \cos(k_ya) \right].$$

The saddle points at $(\pi/a,0)$ and $(0,\pi/a)$ have energy $E_c=0$. The exact spinless density per cell can be written

$$\rho_c(E) = \frac{ \Theta(4|t|-|E|) }{ 2\pi^2|t| } K \left[ 1- \left( \frac{E}{4t} \right)^2 \right],$$

where the complete elliptic integral is defined by

$$K(m) = \int_0^{\pi/2} \frac{d\theta}{ \sqrt{ 1-m\sin^2\theta } }.$$

$K(m)$ diverges logarithmically as $m\to1$, giving the van Hove singularity at $E=0$. At the band edges $E=\pm4|t|$, the density approaches the finite two-dimensional parabolic-edge value.

Example: Two-Dimensional Dirac Cone

For one isotropic cone,

$$\varepsilon_{\pm}(\mathbf q) = \pm\hbar v|\mathbf q|,$$

the density per area, including a declared degeneracy $g$, is

$$\rho_A(E) = \frac{ g|E| }{ 2\pi\hbar^2v^2 }.$$

It vanishes linearly at the node. A zero density at one energy does not imply a finite band gap: a Dirac semimetal has states arbitrarily close to the nodal energy.

Flat and Nearly Flat Bands

An exactly flat band with $w$ states per cell at $E_f$ contributes

$$\rho_c^{\mathrm{flat}}(E) = w\delta(E-E_f).$$

A small bandwidth $W$ spreads this weight over a narrow interval, with typical scale $w/W$. The resulting curve can still contain divergent van Hove singularities; finite bandwidth does not guarantee a finite peak height. Interactions and disorder can qualitatively reorganize a partially filled flat band, so the noninteracting density alone cannot determine the resulting phase.

Projected and Local Densities

For normalized Bloch eigenvectors in an orthonormal orbital basis, the orbital-projected density per cell is

$$\rho_{\alpha}(E) = \frac{1}{N_c} \sum_{n,\mathbf k} \left| u_{n\alpha}(\mathbf k) \right|^2 \delta \left[ E-\varepsilon_n(\mathbf k) \right].$$

Completeness gives

$$\sum_{\alpha} \rho_{\alpha}(E) = \rho_c(E).$$

Individual orbital projections depend on the retained basis and on unitary rotations among orbitals. They are useful diagnostics, not basis-independent observables.

The real-space local density of states is

$$\rho(\mathbf r,E) = \sum_{\lambda} \left| \psi_{\lambda}(\mathbf r) \right|^2 \delta(E-E_{\lambda}).$$

It resolves surfaces, defects, sublattices, and orbital texture. Integrating over the sample returns the total density of states.

Green Functions and Interactions

For an energy-resolved retarded electron Green function normalized so that a noninteracting pole is $(E-\varepsilon+i0^+)^{-1}$,

$$\rho(E) = -\frac{1}{\pi} \operatorname{Im} \operatorname{Tr} G^R(E),$$

with the trace and normalization chosen for total, local, or per-cell density. In a noninteracting band model, poles reproduce delta functions. Interactions can shift and broaden peaks, reduce pole residues and create incoherent continua. For a fixed noninteracting disordered Hamiltonian the exact finite-system spectrum still consists of delta functions; averaging over samples or finite resolution can smooth it.

The resulting single-particle spectral density is not the many-body level density

$$\sum_m \delta(E-E_m^{\mathrm{many-body}}),$$

which counts entire many-body eigenstates and has very different scaling. Spectral Functions and Green Functions and Density of States develop that distinction.

Filling and Thermodynamics

For noninteracting fermions, the particle number per cell is

$$\nu(T,\mu) = \int_{-\infty}^{\infty} dE\, \rho_c(E) f(E-\mu),$$

where

$$f(E-\mu) = \frac{1}{ e^{(E-\mu)/(k_{\mathrm B}T)} +1 }.$$

At zero temperature,

$$\nu(0,\mu) = \mathcal N_c(\mu).$$

If there is a delta-function contribution exactly at $\mu$, its occupation must be specified separately; the zero-temperature step-function convention alone does not fix an arbitrary partial filling.

The internal energy per volume of a noninteracting band system is

$$\frac{U}{V} = \int dE\, E \rho_V(E) f(E-\mu).$$

If $\rho_V(E)$ is smooth on the scale $k_{\mathrm B}T$ around the Fermi energy,

$$\frac{C_V}{V} = \frac{\pi^2}{3} k_{\mathrm B}^2 T \rho_V(E_{\mathrm F}) + O(T^3).$$

For a Landau Fermi liquid the analogous leading coefficient uses the quasiparticle level density, $\rho_{\mathrm{qp}}(0)=\sum_{n,\mathbf k}\delta(\varepsilon^*_{n\mathbf k})/V$, with quasiparticle energies measured from $\mu$. It is generally not the electron spectral density $-\operatorname{Im}\operatorname{Tr}G^R/\pi$: the coherent part of that density contains an extra residue $Z_{n\mathbf k}$. For example, a momentum-independent self-energy giving $\varepsilon^*=Z\xi$ yields $\rho_{\mathrm{qp}}(0)=\rho_0(0)/Z$, whereas the coherent spectral density at zero energy is $Z\rho_{\mathrm{qp}}(0)=\rho_0(0)$. Specific heat is enhanced even though that spectral height is not. Coleman develops this Fermi-liquid distinction.

A nearby van Hove singularity or narrow band can invalidate the smooth-density approximation. In an interacting equilibrium calculation, particle number can still be obtained from the full electron spectral density times the Fermi function, but the simple $E\rho f$ formula above is not a general interacting-energy formula.

The density of states is not an occupation function. Empty and filled states contribute equally to $\rho$; occupation enters through $f$, a Bose factor, or a nonequilibrium distribution.

Spectroscopy and Response

Tunneling

Under restrictive assumptions of weak tunneling, slowly varying matrix elements, and a featureless counterelectrode,

$$\frac{dI}{dV}$$

tracks a thermally broadened local density of states near energy $eV$. Tip orbital symmetry, barrier transmission, temperature, and many-body tunneling effects can reshape the signal.

Photoemission

Angle-resolved photoemission measures an occupied, matrix-element-weighted spectral function,

$$I(\mathbf k,\omega) \propto |M(\mathbf k,\omega)|^2 f(\omega) A(\mathbf k,\omega),$$

not the bare density of states. Momentum integration can resemble an occupied spectral density only after accounting for matrix elements and resolution.

Optical absorption

The joint density of states contains pairs of bands separated by photon energy:

$$J(\hbar\omega) = \sum_{n,m} \int_{\mathrm{BZ}} \frac{d^dk}{(2\pi)^d} \left[ f_n(\mathbf k)-f_m(\mathbf k) \right] \delta \left[ \varepsilon_m(\mathbf k) - \varepsilon_n(\mathbf k) - \hbar\omega \right].$$

It supplies phase space, but optical intensity also requires polarization-dependent transition matrix elements and selection rules. A peak in $J$ need not appear in a forbidden optical channel.

Transport

Conductivity is not determined by $\rho(E_{\mathrm F})$ alone. It weights states by velocity, scattering, and sometimes geometric factors. A large density caused by very flat states can coexist with poor charge transport.

Numerical Evaluation

Histogram

Sample $\varepsilon_n(\mathbf k)$ on a mesh and count states in energy bins. The result depends on bin width and can miss narrow or singular features unless the mesh is converged.

Explicit broadening

Replace each delta function by a normalized kernel. A Gaussian choice is

$$\delta_{\eta}^{\mathrm G}(x) = \frac{1}{ \eta\sqrt{\pi} } e^{-x^2/\eta^2},$$

while a Lorentzian is

$$\delta_{\eta}^{\mathrm L}(x) = \frac{1}{\pi} \frac{\eta}{ x^2+\eta^2 }.$$

The width $\eta$ is part of the result and must be reported. Gaussian instrumental resolution and Lorentzian lifetime broadening have different physical interpretations.

Interpolation methods

Tetrahedron and related simplex methods interpolate bands within momentum-space cells and integrate the delta function analytically. Adaptive meshes focus resolution near band extrema, crossings, and critical points. Wannier interpolation can cheaply generate dense meshes when its band subspace is validated.

Band Structure Workflows explains convergence of the underlying bands. For a DOS calculation, vary both the momentum mesh and the energy resolution: a smooth-looking curve alone does not show that a narrow feature is resolved.

Sum-rule checks

The zeroth moment must approach $N_b$ per cell when all bands and kernel tails are included. For the unbroadened finite density, the first moment is

$$\int dE\,E\rho_c(E)=\frac1{N_c}\sum_{\mathbf k}\operatorname{Tr}h(\mathbf k).$$

For a complete orthonormal tight-binding basis, the thermodynamic Brillouin-zone average equals $\operatorname{Tr}t(\mathbf0)$. On a finite mesh, longer-range hopping harmonics can alias into the constant term; the mesh must resolve them before using that equality as a check.

A normalized Gaussian centered on each level preserves both zeroth and first moments. A Lorentzian preserves the zeroth moment but has no ordinary first moment: its positive- and negative-energy contributions diverge separately. A symmetric principal value is a different prescription. Thus a Lorentzian-broadened plot should not be tested against a purported convergent first-moment integral. Finite energy windows can also lose substantial Lorentzian tail weight.

Exercises

Exercise 1: dimensional threshold law

Derive the exponent of the density of states for

$$E-E_0 = A k^s, \qquad A>0,$$

in $d$ dimensions.

Solution

The number of states below energy $E$ is proportional to the volume of a $d$-ball of radius

$$k = \left( \frac{E-E_0}{A} \right)^{1/s}.$$

Therefore

$$\mathcal N(E) \propto (E-E_0)^{d/s} \Theta(E-E_0).$$

Differentiating,

$$\rho(E) \propto (E-E_0)^{d/s-1} \Theta(E-E_0).$$

The parabolic case has $s=2$, giving $d/2-1$. A two-dimensional Dirac cone has $s=1$, giving a linear density.

Exercise 2: logarithm from a saddle

Show qualitatively why

$$E-E_c = \alpha q_x^2-\beta q_y^2$$

produces a logarithmic density in two dimensions.

Solution

Rescale coordinates so the dispersion is

$$\xi = x^2-y^2.$$

Introduce

$$u = x+y, \qquad v = x-y,$$

so $\xi=uv$. Up to a constant Jacobian, the singular integral is

$$\int du\,dv\, \delta(\xi-uv).$$

Integrating over $v$ gives

$$\int \frac{du}{|u|}$$

over limits set by the local momentum cutoff and by $|uv|=|\xi|$. The result is

$$\rho_{\mathrm{sing}}(\xi) \propto \ln \left( \frac{\Lambda_E}{|\xi|} \right).$$

The cutoff fixes the nonsingular additive part, while the logarithmic dependence is universal for a nondegenerate two-dimensional saddle.

Exercise 3: broadening preserves state count

Show that either normalized kernel preserves state count. For one level at zero energy broadened by a Lorentzian, find the fraction of weight retained in the window $[-W,W]$ and the positive-half first moment up to $W$.

Solution

Both kernels satisfy

$$\int_{-\infty}^{\infty} dx\, \delta_\eta(x) = 1.$$

For a broadened finite density,

$$D_\eta(E) = \sum_\lambda \delta_\eta(E-E_\lambda),$$

one has

$$\int dE\, D_\eta(E) = \sum_\lambda \int dE\, \delta_\eta(E-E_\lambda) = \sum_\lambda 1.$$

Broadening redistributes spectral weight in energy but must not create or remove total weight. A finite plotting window can appear to violate the rule by cutting off kernel tails.

For the Lorentzian,

$$\begin{aligned} \int_{-W}^{W}\delta_\eta^{\mathrm L}(E)\,dE &=\frac2\pi\arctan(W/\eta),\\ \int_0^W E\delta_\eta^{\mathrm L}(E)\,dE &=\frac{\eta}{2\pi}\ln[1+(W/\eta)^2]. \end{aligned}$$

The first expression tends to one; the second diverges logarithmically. At $W=10\eta$, the displayed window retains only about $93.65\%$ of the weight.

Exercise 4: Equal DOS, different conductivity

Compare two two-dimensional parabolic bands with the same degeneracy, band edge, particle density and constant relaxation time $\tau$. One has masses $(m,m)$; the other has $(m/r,mr)$, where $r>0$. Compare their DOS and diagonal conductivities.

Solution

For $E=\hbar^2(k_x^2/m_x+k_y^2/m_y)/2$, the substitution $p_i=k_i/\sqrt{m_i}$ makes constant-energy ellipses circular. Its Jacobian is $\sqrt{m_xm_y}$, so

$$\rho_A(E)=\frac{g\sqrt{m_xm_y}}{2\pi\hbar^2}\Theta(E).$$

Both bands have the same DOS and the same Fermi energy at fixed density. In the relaxation-time approximation, however,

$$\sigma_{xx}=\frac{ne^2\tau}{m_x},\qquad \sigma_{yy}=\frac{ne^2\tau}{m_y}.$$

Relative to the isotropic band, the anisotropic conductivities are multiplied by $r$ and $1/r$. At $r=9$ their ratio is $81$, although the DOS curves coincide exactly. The missing information is the directional velocity distribution. Effective Mass develops these mass distinctions.

References

• N. W. Ashcroft and N. D. Mermin, Solid State Physics (Holt, Rinehart and Winston, 1976), Chapters 2, 8, and 12.
• P. E. Blöchl, O. Jepsen, and O. K. Andersen, “Improved tetrahedron method for Brillouin-zone integrations,” Physical Review B 49, 16223–16233 (1994), doi:10.1103/PhysRevB.49.16223.
• P. Coleman, Introduction to Many-Body Physics (Cambridge University Press, 2015), doi:10.1017/CBO9781139020916.
• E. N. Economou, Green’s Functions in Quantum Physics, 3rd ed. (Springer, 2006), Chapters 1–3.
• P. A. Lee, “Tight Binding and van Hove Singularity,” lecture notes for Theory of Solids I, MIT OpenCourseWare (2004), course materials.
• L. Van Hove, “The Occurrence of Singularities in the Elastic Frequency Distribution of a Crystal,” Physical Review 89, 1189–1193 (1953), doi:10.1103/PhysRev.89.1189.