🔬 materials science

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

This paper presents a perturbative second-order optical susceptibility calculation method for bulk materials using non-orthogonal pseudoatomic orbitals and symmetry-enforced Slater-Koster-like integrals in the velocity gauge, successfully validated on cubic silicon carbide and gallium arsenide.

Original authors: Angiolo Huaman, Luis Enrique Rosas-Hernandez, Salvador Barraza-Lopez

Published 2026-02-03

📖 5 min read🧠 Deep dive

CC BY 4.0

Original authors: Angiolo Huaman, Luis Enrique Rosas-Hernandez, Salvador Barraza-Lopez

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine you are trying to predict how a material will react when you shine a very specific, powerful light on it. Specifically, you want to know if the material can take two photons (packets of light) and smash them together to create one new photon with double the energy. This is called Second Harmonic Generation (SHG). It's the magic trick behind things like green laser pointers and advanced medical imaging, but to design better materials for these tricks, scientists need to calculate a complex number called the second-order optical susceptibility (let's call it $\chi^{(2)}$ ).

For a long time, scientists had two main ways to do these calculations:

The "Plane Wave" Method: Imagine trying to describe a bumpy landscape by laying down a giant, perfectly flat sheet of graph paper over it. You have to use a massive amount of grid squares (computing power) to capture the tiny bumps and valleys. This is accurate but computationally expensive.
The "Localized" Method: Imagine describing that same landscape by placing small, custom-shaped clay models only where the hills and valleys actually are. This is much more efficient, but for a long time, the math for doing this with light was tricky and often required a messy middle-step called "Wannierization" (translating the clay models into a different language just to do the math).

What This Paper Does
The authors of this paper, Angiolo Huamán and colleagues, have built a new, streamlined tool to calculate this light-matter interaction using the efficient "clay model" approach, but without the messy translation step.

Here is the breakdown of their approach using simple analogies:

1. The Building Blocks: "Pseudo-Atomic Orbitals" (PAOs)

Instead of using giant flat sheets, the team uses PAOs. Think of these as tiny, fuzzy clouds of electron probability centered right on top of each atom in the material (like Silicon or Carbon).

The Old Way: To calculate how light moves between these clouds, previous methods often required converting these clouds into a different mathematical format first.
The New Way: The authors say, "Let's just do the math directly on these clouds." They use a method called perturbation theory, which is like asking, "If I nudge this electron cloud slightly with light, how does it wiggle?"

2. The "Slater-Koster" Shortcut: Using Symmetry as a Cheat Code

The hardest part of the math is calculating how these electron clouds interact with each other across space. It's like trying to calculate the wind resistance between two specific trees in a forest. If you have a forest with 1,000 trees, calculating every single pair is a nightmare.

The authors realized that nature is symmetrical.

The Analogy: Imagine you are in a perfectly symmetrical room. If you know how a ball bounces off the floor in the center, you automatically know how it bounces off the floor in the corners because the room is symmetrical. You don't need to measure every corner; you just measure one and apply the rule.
The Paper's Trick: They use symmetry to identify which interactions are identical and which are zero. They calculate a few "master" interactions (called two-center integrals) and then use symmetry rules to fill in the rest of the map. This saves a massive amount of computer time.

3. The "Non-Orthogonal" Twist

In math, "orthogonal" usually means things are at right angles and don't interfere with each other. In this specific type of chemistry software (called SIESTA), the electron clouds do overlap and interfere (they are "non-orthogonal").

The Challenge: Most standard math tools break when things overlap.
The Solution: The authors developed a specific set of equations that handle this overlap naturally. They treat the overlapping clouds like a team of people passing a ball; they account for the fact that the ball is being held by two people at once, rather than pretending it's only with one.

4. Testing the Tool

To prove their new calculator works, they tested it on two famous materials:

Silicon Carbide (3C-SiC): A very hard, durable material used in electronics.
Gallium Arsenide (GaAs): A common material used in lasers and solar cells.

They ran their new "clay model" calculator and compared the results to:

Older, more expensive "flat sheet" (plane-wave) calculations.
Results from other established scientific papers.

The Result: Their new method matched the expensive, heavy-duty calculations almost perfectly but did it much faster and without needing the extra "translation" steps.

Summary

This paper is essentially a new, highly efficient instruction manual for a specific type of computer simulation. It tells scientists how to predict how materials will bend and twist light using a "local" approach (focusing on individual atoms) rather than a "global" approach (looking at the whole crystal at once).

By using symmetry as a shortcut and handling overlapping electron clouds correctly, they have made it easier and faster to design new materials for:

Telecommunications: Sending data faster with light.
Metrology: Measuring things with extreme precision.
Quantum Information: Creating entangled pairs of photons for future quantum computers.

The paper does not claim to have built a new laser or a new quantum computer; it simply provides a better, faster way to do the math required to design the materials that will eventually power those technologies.

Technical Summary: Perturbative Second-Order Optical Susceptibility of Bulk Materials

Problem Statement
The calculation of the second-order optical susceptibility, $\chi^{(2)}_{ijk}(-2\omega; \omega, \omega)$ , is critical for applications in metrology, spectroscopy, telecommunications, and quantum information. While theoretical formalisms for computing this quantity in semiconductors generally follow perturbation theory (within the independent-particle approximation, IPA) or time-dependent density functional theory (TDDFT), most existing methodologies rely on orthogonal basis sets, such as plane waves or post-processed Wannier functions. Although recent tools have begun recasting these problems using localized orbitals (Wannier functions), there remains a need for robust frameworks that utilize non-orthogonal localized basis sets directly, particularly for simulations involving significant vacuum regions (e.g., molecules, surfaces, 2D materials) where Wannierization post-processing is not required.

Methodology
The authors present a formalism for calculating $\chi^{(2)}_{ijk}$ within the IPA and perturbation theory in the velocity gauge, relying on numerical pseudo-atomic orbitals (PAOs) from the SIESTA DFT package. The approach, termed "Nonlinear Algorithms by Post Processing in SIESTA," proceeds through the following steps:

Construction of Bloch Eigenstates: The method retrieves auxiliary radial functions from SIESTA's output files and reconstructs PAOs in real space. These are combined with spherical harmonics (using the Condon-Shortley phase convention) to form a non-orthogonal basis set. Bloch sums are constructed from these localized orbitals to solve the generalized eigenvalue problem $H_0(k)C^0_n(k) = \epsilon^0_n(k)S_0(k)C^0_n(k)$ , yielding unperturbed eigenvalues and eigenvectors.
Perturbative Expansion: The second-order susceptibility is decomposed into two contributions, $A_{ijk}(2\omega)$ and $B_{ijk}(\omega)$ , derived from a perturbative expansion of the density matrix up to second order in the vector potential. These expressions involve sums over valence and conduction states and integrals over the first Brillouin zone (FBZ).
Symmetry-Enforced Two-Center Integrals: A core technical innovation is the calculation of momentum matrix elements between Bloch states. Instead of relying on Fourier transforms or position operators, the authors compute "Slater-Koster-like" two-center integrals of the momentum operator ( $\partial_i$ $\partial_{i}$ ) between PAOs.
- The method utilizes the spatial symmetries of the crystal to rotate the integration axis such that the vector connecting two orbitals aligns with the $z$ -axis.
- This rotation allows the reduction of general two-center integrals into a small subset of independent, non-zero integrals along the $z$ -axis.
- Symmetry relations (parity and angular momentum conservation) are used to identify which integrals vanish and which are related, significantly reducing the number of numerical integrations required.
Brillouin Zone Reduction: The integration over the FBZ is reduced to the irreducible first Brillouin zone (IFBZ) and further to the time-reversal symmetric irreducible zone ($trIFBZ$) by applying point group operations and time-reversal symmetry constraints.
Kramers-Kronig Transformation: The real part of the susceptibility is obtained from the imaginary part via the Kramers-Kronig relation, utilizing a principal value integral.

Key Contributions

Non-Orthogonal Basis Framework: The work establishes a rigorous perturbative framework for $\chi^{(2)}$ calculations using non-orthogonal PAOs, avoiding the need for Wannierization.
Symmetry-Driven Integration: The paper introduces a specific algorithmic approach to evaluate momentum matrix elements by exploiting crystal symmetries to reduce complex two-center integrals to a minimal set of independent integrals along a symmetry axis.
Implementation in SIESTA: The method is implemented as a post-processing tool for SIESTA, leveraging its existing PAO basis sets and output formats (ion files, HSX/HS files).
Benchmarking: The approach is tested on paradigmatic cubic zinc-blende materials: Silicon Carbide (3C-SiC) and Gallium Arsenide (GaAs).

Results
The methodology was applied to 3C-SiC and GaAs using various basis sets (SZP, DZP, TZP).

Band Structure and DOS: The calculated band structures and density of states for 3C-SiC using the DZP basis set showed fair agreement with plane-wave (VASP) calculations and experimental indirect band gaps.
Convergence: The imaginary part of $\chi^{(2)}_{xyz}$ was found to converge with respect to the number of $k$ -points in the $trIFBZ$, the number of nearest neighbors included in the two-center integral summation, and the basis set size.
Comparison: The calculated $\chi^{(2)}$ spectra for 3C-SiC were compared with previous literature (specifically Ref. [65]), showing consistency in the onset of susceptibility and peak positions.
Linear Response: The framework also successfully computed linear optical responses, serving as a benchmark for the underlying formalism.

Significance and Claims
The authors position this work as a "first step" toward eventual incorporation of many-body effects (such as excitonic interactions via GW-BSE) within a non-orthogonal PAO framework. The primary significance claimed is the ability to perform second-order optical susceptibility calculations for systems with significant vacuum (molecules, nanotubes, surfaces, 2D materials) without the computational overhead or approximation errors associated with Wannierization. The method is particularly suited for symmetry analysis of Second Harmonic Generation (SHG) responses in 2D materials. The paper emphasizes that while it currently operates within the IPA, the formalism provides a necessary foundation for extending these calculations to include electron-hole interactions in the future.

1. The Building Blocks: "Pseudo-Atomic Orbitals" (PAOs)

2. The "Slater-Koster" Shortcut: Using Symmetry as a Cheat Code

3. The "Non-Orthogonal" Twist

4. Testing the Tool

Summary

More like this