Dynamical Effective Hamiltonian Approach to Second-Harmonic Generation in Quantum Magnets: Application to NiI$_2$

Although second harmonic generation (SHG) is a promising and widely used method recently for studying 2D magnetic materials, the quantitative analysis of the full SHG tensor is currently challenging. In this letter, we describe a first-principles-based approach towards quantitative analysis of SHG in insulating magnets through formulation in terms of dynamical effective operators. These operators are computed by solving local many-body cluster models. We benchmark this method on NiI$_2$, a multiferroic 2D van der Waals antiferromagnet, demonstrating quantitative analysis of reported Rotational Anisotropy (RA)-SHG data. SHG is demonstrated to probe local ring-current susceptibilities, which provide sensitivity to short-range chiral spin-spin correlations. The described methods may be easily extended to other non-linear optical responses and materials.

💡 Research Summary

This paper introduces a first‑principles‑based methodology for quantitatively analyzing second‑harmonic generation (SHG) in insulating quantum magnets. The authors formulate SHG in terms of dynamical effective operators that encapsulate the low‑energy impact of optical coupling to high‑energy electronic excitations. These operators are obtained by exact diagonalization of local many‑body cluster models that include both transition‑metal d‑orbitals and ligand orbitals represented by a minimal set of natural transition orbitals (NTOs).

The theoretical framework starts from the electric‑dipole light‑matter interaction and expresses the SHG susceptibility χμξν(2ω; ω, ω) as a thermal average of material‑specific operators Ōμξν(ω). For magnetic insulators, Ō can be expanded in spin operators; symmetry analysis shows that in a centrosymmetric crystal such as NiI₂ the spin‑independent term vanishes, and the linear spin term is forbidden by inversion symmetry. Consequently, the leading contribution is quadratic in spins and takes the form of a vector product (Si × Sj) multiplied by a complex, frequency‑dependent coefficient Cijμξν(ω), which is the dynamical analogue of a Dzyaloshinskii‑Moriya (DM) vector.

Computationally, the authors first perform DFT (FLEUR) calculations and construct Wannier functions that mix Ni d‑ and I p‑character. Transition matrices between metal‑centered Wannier functions and ligand‑centered NTOs are singular‑value‑decomposed to obtain a compact basis that captures the full optical weight of ligand‑to‑metal charge‑transfer (LMCT) excitations. A nearest‑neighbor Ni–Ni cluster is then built, incorporating on‑site Hubbard U and double‑counting corrections (AMF). The many‑body Hamiltonian is diagonalized in the truncated d + NTO Fock space (≈5 × 10⁵ states) using Krylov subspace techniques. Low‑energy eigenstates are projected onto pure spin states, establishing a mapping from electronic operators to spin operators. Finally, the double‑resonant expression for Ō (Eq. 3) is evaluated for all polarization combinations, yielding the full set of C‑vectors.

Applying this workflow to NiI₂, a van‑der‑Waals antiferromagnet with a centrosymmetric (\bar{R}3m) structure, the authors find that only nearest‑neighbor bonds contribute appreciably to SHG. The computed C‑vectors obey the C₂h symmetry of each edge‑sharing Ni–I–Ni bond: X and Z components are finite for polarization triples (xxx, xyy, yxy, yyx), while Y components are finite for (yyy, yxx, xyx, xxy). The magnitude and phase of each component are dictated by the interplay of intersite d‑d excitations (≈ 3 eV) and ligand‑to‑metal charge‑transfer excitations (≈ 1.5–4 eV). At the experimental photon energy ℏω ≈ 1.25 eV (λ ≈ 991 nm), the condition 2ℏω ≈ ΔE_LMCT is satisfied, and the dominant SHG pathway involves a non‑resonant d‑d transition followed by a resonant LMCT transition, with spin‑orbit coupling on the iodine ligands providing the necessary angular‑momentum transfer. Detailed Feynman‑type diagrams (Fig. 3 i–l) illustrate how hole transfer, ligand SOC (Lz Sz), and photon emission combine to generate specific C‑vector components.

Using the calculated operators, the authors fit rotational anisotropy SHG (RA‑SHG) data from Ref.