Berry Curvature Numerics

Numerical eigenvectors carry arbitrary phases, so differentiating them directly can produce a meaningless Berry curvature. Gauge-invariant overlaps, projectors, or velocity-matrix formulas avoid that instability when their assumptions are satisfied. For the geometric interpretation, see Berry curvature; here we develop its numerical evaluation.

Helpful background. Band-structure numerics explains isolated bands and subspace tracking.

Connection convention and gauge cancellation

Use the connection convention

$$\begin{aligned} \mathcal A_i&=i\langle u|\partial_i u\rangle,\\ \Omega_{xy}&=\partial_x\mathcal A_y-\partial_y\mathcal A_x\\ &=-2\,\operatorname{Im}\langle\partial_xu|\partial_yu\rangle. \end{aligned}$$

Under $u\mapsto e^{i\chi}u$, $\mathcal A_i\mapsto\mathcal A_i-\partial_i\chi$, while curvature is unchanged. Independent eigensolver phases can make finite differences of $u$ large even for a smooth projector.

For a normalized isolated band define link phases

$$U_\mu(\boldsymbol k)= \frac{\langle u(\boldsymbol k)|u(\boldsymbol k+\Delta_\mu)\rangle} {|\langle u(\boldsymbol k)|u(\boldsymbol k+\Delta_\mu)\rangle|}.$$

On a positively oriented $x$-then-$y$ plaquette form

$$P_\square= U_x(\boldsymbol k)U_y(\boldsymbol k+\Delta_x) U_x(\boldsymbol k+\Delta_y)^*U_y(\boldsymbol k)^*.$$

Every vertex phase cancels. Since a forward overlap is approximately $e^{-i\mathcal A_\mu\Delta k_\mu}$,

$$\Omega_{xy}(\text{plaquette center}) \approx-\frac{\operatorname{Arg}P_\square}{\Delta k_x\Delta k_y}.$$

The minus sign is tied to $\mathcal A=+i\langle u|\nabla u\rangle$. Reversing the plaquette orientation also reverses the result.

Projector and spectral formulas

For $P=|u\rangle\langle u|$,

$$\Omega_{xy}=i\,\operatorname{Tr}\!\left(P[\partial_xP,\partial_yP]\right).$$

Finite differences of projectors avoid phase fixing but still have mesh error. For an isolated eigenvalue in a fixed orthonormal basis, differentiation of the eigenproblem gives

$$\Omega_{n,xy}=-2\,\operatorname{Im}\sum_{m\ne n} \frac{\langle u_n|\partial_xH|u_m\rangle \langle u_m|\partial_yH|u_n\rangle} {(E_n-E_m)^2}.$$

Near small gaps the denominator amplifies both physical curvature and numerical error. At a degeneracy an individual-band curvature is not defined; use the isolated multi-band projector instead.

The spectral expression follows by projecting the differentiated eigenvalue equation onto another band: $\langle u_m|\partial_i u_n\rangle=\langle u_m|\partial_iH|u_n\rangle/(E_n-E_m)$. Inserting this into the curvature and canceling the component parallel to $u_n$ gives the displayed result. This provides an independent implementation from overlap plaquettes.

For a group of occupied bands, the projector formula gives the total curvature of that group. An equivalent spectral sum uses $n$ in the group and $m$ outside it; contributions between two occupied bands cancel in pairs. Do not evaluate singular individual-band denominators at an internal degeneracy and hope that their floating-point infinities cancel. Only the gap to excluded states belongs in this grouped calculation.

A momentum-dependent basis requires its connection or correctly constructed overlaps. Differentiating only the coefficient vector omits basis geometry. Keep the orbital convention from tight-binding assembly.

Two-level local benchmark

For $H(\boldsymbol k)=\boldsymbol d(\boldsymbol k)\cdot\boldsymbol\sigma$, the lower-band projector is $(I-\hat{\boldsymbol d}\cdot\boldsymbol\sigma)/2$. Substitution yields

$$\Omega_{-,xy} =\frac12\hat{\boldsymbol d}\cdot (\partial_x\hat{\boldsymbol d}\times\partial_y\hat{\boldsymbol d}).$$

For $\boldsymbol d=(k_x,k_y,m)$,

$$\Omega_{-,xy} =\frac{m}{2(k_x^2+k_y^2+m^2)^{3/2}}.$$

Compare overlap plaquettes centered near the origin with this value as the mesh is refined. The upper-band result has opposite sign. This is a local benchmark; integrating this continuum model over the plane does not by itself define an integer lattice Chern number.

At $m=1$ and the origin, the lower curvature is $+1/2$ and the upper curvature is $-1/2$. Compare three independent routes: the analytic expression, the spectral derivative formula with $\partial_xH=\sigma_x$ and $\partial_yH=\sigma_y$, and centered finite differences of the projector. A square overlap plaquette centered at the same point should converge to these values with second-order local error for this smooth model.

When $|m|$ decreases, the curvature peak narrows on momentum scale $|m|$ and its magnitude grows as $1/(2m^2)$. A mesh adequate at $m=1$ need not resolve $m=0.01$. The relevant refinement parameter is spacing relative to the curvature scale, not spacing alone.

Mesh diagnostics

Record the smallest neighbor-overlap magnitude, spectral gap, and largest plaquette phase magnitude. Nearly zero overlaps make phase extraction unreliable. Plaquette phases close to $\pi$ can alias a larger underlying flux into the principal branch.

Apply random independent phases to all eigenvectors and recompute; curvature should be unchanged to roundoff. Refine the mesh and compare local curvature, not only its integral. An integer Chern number can remain stable while the curvature profile is poorly resolved.

Exercise

Why does normalizing an overlap not solve the problem when its magnitude is nearly zero?

Solution

Normalization divides by a tiny number, so small absolute eigenvector errors create order-one phase errors. The remedy is a finer mesh, a correctly selected isolated subspace, or a different representation—not merely replacing the denominator by a small cutoff.

References

• Fukui, T., Hatsugai, Y., and Suzuki, H., “Chern numbers in discretized Brillouin zone: efficient method of computing (spin) Hall conductances,” Journal of the Physical Society of Japan 74, 1674–1677 (2005), arXiv:cond-mat/0503172.
• Xiao, D., Chang, M.-C., and Niu, Q., “Berry phase effects on electronic properties,” Reviews of Modern Physics 82, 1959–2007 (2010), doi:10.1103/RevModPhys.82.1959.