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

TL;DR

This work develops a perturbative framework to compute the second-order optical susceptibility $χ^{(2)}(-2ω;ω,ω)$ using non-orthogonal PAOs in the velocity gauge. By enforcing crystal symmetries, two-center Slater–Koster-like momentum integrals are efficiently evaluated, and the irreducible Brillouin zone is exploited to reduce computational cost. The method is validated on 3C-SiC and GaAs, showing convergence with DZP basis, adequate $k$-point sampling, and agreement with literature, including a scissor correction for GaAs. The approach offers a compact, symmetry-guided route to nonlinear optical properties for systems with vacuum regions and lays groundwork for including many-body effects in future work.

Abstract

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.

Figures (11)

• Figure 1: Typical radial functions: (a) $R^{aux}_{1,\textsf{l}(1),\textsf{ir}(\textsf{l},1)}(r)$ for Si; $a_B=0.529177$ Å is the Bohr radius. (b) $R^{aux}_{2,\textsf{l}(2),\textsf{ir}(\textsf{l},2)}(r)$ for C. Cutoff radii are written, too.
• Figure 2: (a) 3C-SiC conventional (cubic) unit cell (with lattice parameter $a$ shown) and primitive unit cell (parallelepiped). (b) FBZ with high-symmetry points indicated. The irreducible first Brillouin zone including time reversal symmetry is highlighted, too.
• Figure 3: The heart of this development is to calculate linear momentum integrals among PAOs in highly-symmetric configurations, and to use symmetries to determine actual integrals on materials.
• Figure 4: Examples of two-center integrals (equation \ref{['d_integrals']}), as a function of the number of integration points ngp$^3$, which is another input in our code ($a_1=1.911$ Å). Convergence of the order of $10^{-3}$ is achieved for $\textsf{ngp}=60$.
• Figure 5: Flowchart showing the step-by-step procedure used in this paper to obtain the second-order optical susceptibility.