Berry Flux Diagonalization

• Berry Flux Diagonalization is a computational framework that computes electric polarization changes and topological invariants through the diagonalization of Berry flux operators on k-space plaquettes.
• The method decomposes global phase shifts into locally gauge-invariant increments, bypassing the need for dense interpolation and enhancing efficiency in ab initio calculations.
• Benchmark studies show that BFD achieves high accuracy in predicting polarization and topological properties, with results agreeing within 1–2% of traditional methods across various materials.

Berry Flux Diagonalization (BFD) is a computational framework for evaluating the change in electric polarization and quantized topological invariants in crystalline systems using the geometric properties of the many-body wavefunction, specifically through direct diagonalization of Berry flux operators on discrete k-space plaquettes. The method addresses inherent branch ambiguities in Berry-phase polarization formulas by decomposing global phase changes into a sum of small, gauge-invariant increments. BFD bypasses the need for intermediate path sampling between two endpoint wavefunctions, yielding computational advantages in efficiency and robustness for both ab initio and model-derived band structures (Bonini et al., 2020, Poteshman et al., 23 Nov 2025, Tyner et al., 2021).

1. Conceptual Foundation and Motivation

The modern theory of polarization relates the change in electronic polarization, $\Delta P$, between two insulating states A and B to the integrated Berry flux, $\Phi$, through a surface in (k, λ) space parameterizing the evolution between A and B:

$\Delta P_{A\to B} = -\frac{e}{2\pi}\Phi$

where

$\Phi = \iint_{S} \Omega(k,\lambda)\,dk\,d\lambda$

with $\Omega(k,\lambda)$ the Berry curvature. In practice, only the endpoints' ground-state wavefunctions are typically accessible, and $\Phi$ can only be computed modulo $2\pi$ unless the gauge is continuous along the path.

Conventional Berry-phase polarization calculations resolve this modulo $2\pi$ ambiguity by constructing a densely sampled interpolation path and tracking the Berry phase so that increments are always $\ll\pi$. BFD circumvents the need for this interpolation by decomposing the global phase difference into a sum of locally gauge-invariant phases across many small k-space plaquettes, under the minimal-evolution assumption that the physical path corresponds to the smallest possible gauge winding (Bonini et al., 2020).

2. Mathematical Structure of BFD

2.1. Discretization and Plaquette Construction

For a uniform k-point mesh along a direction $\alpha$, BFD forms closed four-vertex plaquettes straddling state A (reference) and state B (polarized):

• Vertices: $$(k_j, A) \rightarrow (k_j, B) \rightarrow (k_{j+1}, B) \rightarrow (k_{j+1}, A)$$

Overlap matrices between these vertices are defined as: $$M^{\langle ab \rangle}_{mn} = \langle u_m(k_a)|u_n(k_b) \rangle$$ where $|u_{n}(k)|$ are the occupied Bloch eigenstates at $k$.

Each $M^{\langle ab \rangle}$ is decomposed by singular value decomposition (SVD) as $M=V\Sigma W^\dagger$; the corresponding "unitary approximant" is $\mathcal{M}=VW^\dagger$.

2.2. Wilson Loop and Diagonalization

The plaquette Berry-flux evolution operator is constructed as: $$U_{p} = \mathcal{M}^{\langle12\rangle}\mathcal{M}^{\langle23\rangle}\mathcal{M}^{\langle34\rangle}\mathcal{M}^{\langle41\rangle}$$ $U_p$ is a unitary $N_{\textrm{occ}}\times N_{\textrm{occ}}$ matrix whose eigenvalues are $\{e^{i\varphi_n^p}\}$. Each eigenphase provides a gauge-invariant, branch-consistent phase increment for that plaquette.

The total phase difference, and thus the electronic polarization change, is then: $$\Phi = \sum_p \sum_{n=1}^{N_{\textrm{occ}}} \varphi_n^p \qquad \Delta P_{\textrm{elec}} = -\frac{e}{2\pi} \Phi$$

For non-Abelian (e.g., Kramers-degenerate) bands, the same logic applies, and the non-Abelian Wilson loop diagonalization yields quantized spin-Chern numbers or other topological invariants in terms of the eigenphases (Tyner et al., 2021).

3. Algorithmic Workflow

The canonical BFD workflow comprises the following steps (Bonini et al., 2020, Poteshman et al., 23 Nov 2025):

1. Gauge Alignment: At each $k$, overlap the A and B states, SVD the overlap matrix, and rotate B to maximize overlap (set $\det[M^{AB}(k)] \gt 0$).
2. Plaquette Processing: For each adjacent k-point pair along $\alpha$,
   • Construct all four overlap matrices on the plaquette,
   • SVD to form nearest unitary $\mathcal{M}$ for each edge,
   • Build $U_p$ by multiplying the four $\mathcal{M}$,
   • Diagonalize $U_p$ to extract eigenphases $\{\varphi_n^p\}$.
3. Validation: Ensure all $\Sigma$ singular values $\gt$ threshold (typically $0.15$); all $|\varphi_n^p| \ll \pi$ for branch consistency.
4. Summation: Compute total phase $\Phi$ across all plaquettes and bands.
5. Ionic Contribution: Add

$$\Delta P_\textrm{ion} = (e/\Omega) \sum_i Z_i \Delta r_i$$

for total $$\Delta P = \Delta P_\textrm{elec} + \Delta P_\textrm{ion}$$.

1. Branch Heuristics: If maximal per-atom displacement $\Delta r_\textrm{max} \leq 0.3$Å, no intermediate interpolated structures needed; otherwise, interpolate such that adjacent structures have $\Delta r_\textrm{max}^{adj}\leq 0.3$ Å for reliable gauge alignment (Poteshman et al., 23 Nov 2025). If failure is detected (e.g., $\Sigma_\textrm{min} < 0.15$), augment the path via extra interpolations.

4. Physical Interpretation and Scope

BFD imposes "minimal gauge winding" between endpoints, mirroring the shortest-ionic-path mapping, ensuring that each gauge-invariant increment $|\varphi_n^p|<\pi$. This decomposition matches the experimental switching polarization for typical ferroelectrics under the empirical minimal-evolution assumption. BFD is manifestly applicable to multiband systems, pyroelectrics, antiferroelectrics, heterostructures, and is generalizable to Born–Oppenheimer paths in finite electric fields (Bonini et al., 2020).

The Wilson loop diagonalization extracts full topological content, and the method extends to non-Abelian scenarios, such as computation of $SU(2)$ spin-Chern numbers for Kramers-degenerate bands in PT-symmetric materials (Tyner et al., 2021). The core principle is that the gauge-invariant content of the Berry flux is always encoded in the eigenphases of the associated Wilson loops, regardless of Abelian or non-Abelian structure.

5. Practical Implementation and High-Throughput Screening

Automated BFD pipelines have been deployed for high-throughput computation of electric polarization across databases of materials (Poteshman et al., 23 Nov 2025). The workflow leverages:

• A prior real-space assessment of atomic displacement to estimate required interpolation count.
• DFT evaluation of wavefunctions at endpoints (and minimal intermediates if needed).
• Post-processing to compute BFD contributions, with automated singular value and branch monitoring.
• Compatibility with standard plane-wave DFT codes (VASP, Quantum ESPRESSO, ABINIT).

BFD outperforms traditional dense-interpolation Berry-phase workflows, particularly in high-throughput settings. For example, in a benchmarking of 176 candidate ferroelectrics, BFD required $\leq$4 interpolations in $174/176$ cases, compared with standard workflows using $10-20$, and delivered agreement within $1\%$ for $154/176$ materials; failures due to band-gap closure or poor overlap are automatically flagged and can be remediated by mesh refinement or additional interpolations (Poteshman et al., 23 Nov 2025).

6. Benchmark Data and Limitations

BFD has been benchmarked on prototypical materials such as BaTiO$_3$, KNbO$_3$, PbTiO$_3$, LiNbO$_3$, Bi$_2$MoO$_6$, CrO$_3$, and BiFeO$_3$. Results for the polarization change $\Delta P$ closely match established values from dense-path Berry-phase integration (within $1.5\%-2\%$ across codes and materials) (Bonini et al., 2020, Poteshman et al., 23 Nov 2025).

The key limitation is the requirement of high wavefunction overlap (quantified by $\Sigma_\textrm{min}>0.15$) between adjacent structures; failures arise for paths involving metallic intermediates or large atomic displacements. These are detected by automated branch and overlap diagnostics. The method is memory-intensive due to the need for all k-point wavefunctions but is under optimization for irreducible-zone reading (Poteshman et al., 23 Nov 2025).

7. Generalizations and Future Developments

BFD provides a robust, gauge-fixed, branch-continuous tool for Berry-phase and topological-invariant computations, replacing costly intermediate sampling with endpoint-only evaluation wherever minimal-evolution assumptions are justified. Future extensions include direct calculation in higher-dimensional parameter spaces (e.g., strain, electric/magnetic field), computation of Chern numbers from 2D Brillouin zone plaquettes, and tracking of Weyl node charges by Berry flux diagonalization in 3D momentum-space cubes (Bonini et al., 2020). The non-Abelian generalization, exploiting Wilson loop diagonalization within $SU(N)$ gauge bundles, underpins unified approaches to bulk-boundary correspondence and higher-order topology (Tyner et al., 2021).