Perturbative second-order optical susceptibility of bulk materials: a symmetry-enforced return to non-orthogonal localized basis sets

The second-order optical susceptibility of semiconductors $χ_{ijk}^{(2)}(-2ω;ω,ω)$ finds application in metrology, spectroscopy, telecommunications, material characterization, and quantum information. Pioneering calculations of $χ_{ijk}^{(2)}(-2ω;ω,ω)$ utilized non-orthogonal Gaussian orbitals centered at atoms. That formulation transitioned into plane-wave-based algorithms as time went by. As of late, nevertheless, multiple tools for calculating optical susceptibilities have recast the problem using Wannier ({\em i.e.}, {\em localized}) orbitals, making a comeback onto frameworks based on localized basis sets. Here, we present an approach for calculating $χ_{ijk}^{(2)}(-2ω;ω,ω)$ reliant on numerical pseudoatomic orbitals (PAOs) within perturbation theory in the velocity gauge. Its salient feature is a calculation of `Slater-Koster-like’ two-center integrals of the momentum operator in between PAOs identified by symmetry. The approach was successfully tested on paradigmatic cubic silicon carbide (3C-SiC) and gallium arsenide, for which linear responses are contributed as well.

💡 Research Summary

The manuscript introduces a novel first‑principles framework for computing the second‑order optical susceptibility χ^(2)_{ijk}(−2ω; ω, ω) of bulk semiconductors using non‑orthogonal numerical pseudo‑atomic orbitals (PAOs) within the velocity‑gauge perturbation theory. Historically, χ^(2) calculations have relied on either orthogonal plane‑wave bases or on post‑processed Wannier functions, both of which become cumbersome for systems with large vacuum regions (e.g., molecules, nanowires, 2D layers). By staying within the PAO representation generated by the SIESTA code, the authors avoid Wannierization and exploit the locality of the basis to achieve substantial computational savings.

The core of the method consists of three tightly coupled components: (1) construction of Bloch states from a non‑orthogonal PAO set, (2) evaluation of momentum‑operator matrix elements ⟨φ_α|p̂_i|φ_β⟩ as two‑center integrals, and (3) exploitation of crystal point‑group symmetry to reduce the number of distinct integrals to a Slater‑Koster‑like set of parameters. Each PAO is expressed as a product of a radial function (read from SIESTA ion files) and a real spherical harmonic with the Condon‑Shortley phase. The Bloch sum φ_{αk}(r) = N^{-1/2}∑_R e^{ik·R} φ_α(r−R−δ_α) is built, and the generalized eigenvalue problem H_0(k)C_n(k)=ε_n(k)S_0(k)C_n(k) yields the band energies ε_n(k) and expansion coefficients C_n(k). Because the basis is non‑orthogonal, the overlap matrix S_0(k) must be retained throughout the derivation.

In the velocity gauge, the density matrix is expanded to second order in the monochromatic vector potential A(t). The resulting χ^(2) tensor separates into an “A‑term” (two‑photon resonant processes) and a “B‑term” (one‑photon resonant processes). Both terms are expressed as Brillouin‑zone integrals over products of momentum matrix elements and energy denominators, with explicit δ‑functions enforcing energy conservation. The authors provide compact expressions for Im