Complex-Scaled Discrete Dipole Approximation

• CSDDA is a computational framework that extends DDA by incorporating complex scaling to impose outgoing-wave conditions and damp nonphysical divergences.
• It uses iterative solvers with optimized mixing parameters to achieve numerically stable solutions even for materials with strong anisotropy and high optical losses.
• The method accurately simulates localized surface plasmon resonances and hyperbolic modes, validated by close agreement between simulated and experimental spectral features.

The complex-scaled discrete dipole approximation (CSDDA) is a computational framework for modeling electromagnetic response of nanoscale materials with strong anisotropy, high optical losses, or supporting hyperbolic regimes. Developed in the context of anisotropic klockmannite CuSe nanocrystals, CSDDA extends the conventional discrete dipole approximation by incorporating complex scaling of spatial and material parameters. This approach enforces outgoing-wave boundary conditions, damps nonphysical field divergences, and enables numerically stable solutions even for metallic and hyperbolic domains, facilitating accurate simulation of localized surface plasmon resonances (LSPRs) and complex modal structures in anisotropic nanomaterials (Parekh et al., 26 Dec 2025).

1. Theoretical Basis

The foundation of CSDDA is the frequency-domain vector wave equation for non-magnetic media,

$$\nabla\times\nabla\times \mathbf E(\mathbf r)\;-\;k^{2}\,\varepsilon(\mathbf r)\,\mathbf E(\mathbf r)\;=\;\mathbf 0,$$

where $k = \omega/c$ and $\varepsilon(\mathbf r)$ is the generally tensorial permittivity. In DDA, the scatterer is discretized as $N$ point-dipoles at positions $\mathbf r_n$, each characterized by a polarizability tensor $\boldsymbol\alpha_n$. The induced dipole moments $\{\mathbf p_n\}$ satisfy

$$\mathbf p_n = \boldsymbol\alpha_n\,\mathbf E_{\rm loc}(\mathbf r_n),$$

with the local field

$$\mathbf E_{\rm loc}(\mathbf r_n) = \mathbf E_{\rm inc}(\mathbf r_n) + \sum_{m \neq n} \mathbf G_{nm}\, \mathbf p_m,$$

where $\mathbf G_{nm}$ is the free-space or substrate-modified Green's tensor. Complex scaling is introduced via the transformation $\mathbf r \to \mathbf r e^{i\theta}$, which imposes outgoing-wave boundary conditions and damps far-field components. This modifies each dipole's response such that $\varepsilon \to \varepsilon e^{-2i\theta}$ and the Green’s tensor acquires a factor $e^{ik|\mathbf r_n-\mathbf r_m|(1 - e^{i\theta})}$. The choice of $\theta$ (typically $\sim20^\circ$) is critical for stability but may remain unstated.

2. Iterative Solvers: CSDDA and CSDDA++

The CSDDA algorithm employs an iterative successive-over-relaxation scheme to solve the $3N \times 3N$ dipole system. For isotropic cases,

$$\mathbf E_{\rm loc,n}^{(i+1)} = g^{(i)} \Bigl( \mathbf E_{\rm inc}(\mathbf r_n) + \mathbf E_{\rm scat,sum,n}^{(i)} \Bigr) + (1-g^{(i)})\,\mathbf E_0,$$

where the iterative mixing parameter $g^{(i)}$ is optimized at each step by minimizing the global relative error

$$R^{(i)} = \frac{\sum_n \| \mathbf E_{\rm inc}(\mathbf r_n) + \mathbf E_{\rm scat,sum,n}^{(i)} - \mathbf E_{\rm loc,n}^{(i+1)} \|^2}{\sum_n \|\mathbf E_{\rm loc,n}^{(i+1)}\|^2 }.$$

For highly anisotropic or hyperbolic materials (where propagation and evanescent contributions are comparable), CSDDA++ generalizes the update to four complex mixing parameters $\{a^{(i)}, b^{(i)}, y^{(i)}, 8^{(i)}\}$: $$\mathbf E_{\rm loc,n}^{(i+1)} = a^{(i)}\Bigl[ \mathbf E_{\rm inc}(\mathbf r_n) + \sum_m \mathbf G_{nm}\, \mathbf E_{\rm loc,m}^{(i)} \Bigr] + b^{(i)} \mathbf E_0 + y^{(i)} \sum_m \mathbf G_{nm}\, \mathbf E_{\rm loc,m}^{(i)} + 8^{(i)} \sum_m \mathbf G_{nm}\, \mathbf E_{\rm loc,m}^{(i)},$$ with all parameters chosen to minimize the global error at each iteration.

3. Numerical Implementation and Parameters

CSDDA simulations are performed using a high-resolution cubic dipole mesh:

• Dipole side: $d = 1.5$ Å.
• Mesh dimensions: For a 20 nm triangular prism (edge 19 nm, thickness 5 nm), up to $1.15 \times 10^6$ dipoles are deployed, corresponding to a lattice of $80 \times 80 \times 20$ sites.
• Convergence criterion: Iterations proceed until $R < 10^{-6}$, or to $\sim10^{-4}$ in the metallic regime, where diverging tip fields hinder convergence.
• Anisotropic dielectric response: Each $\boldsymbol\alpha_n$ is assigned from QSGW+RPA-calculated $\varepsilon_{xx}(\omega)$ (in-plane) and $\varepsilon_{zz}(\omega)$ (out-of-plane).
• Substrate effects: The presence of a silica substrate (SiO$_2$) is incorporated via the image-dipole (plate) Green’s function for self-interaction: $$\mathbf G_{nn}^{\rm plate}(\omega) = \frac{1}{32\pi\varepsilon_0 k \zeta} \begin{pmatrix} 1/(\varepsilon(\omega)+1) & 0 & 0\ 0 & 1/(\varepsilon(\omega)+1) & 0\ 0 & 0 & 1 \end{pmatrix}$$ with $\zeta$ the distance to the image plane.
• Boundary conditions: Simulations use $n=1.45$ (hexane) above a semi-infinite SiO$_2$ substrate.

4. Application to Anisotropic Klockmannite CuSe

CSDDA was applied to triangular prism, hexagon, truncated triangle, and disk-shaped CuSe nanocrystals of $\sim19$ nm diameter and thickness $\sim5$ nm, with up to $1.1 \times 10^6$ dipoles per model. The spectral window covered $400$–$1600$ nm (photon energies $0.8$–$3.1$ eV), encompassing visible, hyperbolic, and metallic regimes. Key features:

• In-plane and out-of-plane dielectric constants were derived from QSGW+RPA calculations.
• Complex scaling angle $\theta \sim 20^\circ$ (unstated exact value) was employed to suppress backscattering of evanescent fields.
• Mesh refinement to sub-nanometer scale ($d = 1.5$ Å) was maintained throughout.

5. Validation, Mode Analysis, and Physical Interpretation

Validation was performed via spectral, modal, and field-mapping comparisons:

• Absorption cross-section matching: Simulated orientation-averaged $C_{\rm abs}$ peaks at $1045$ nm (simulated) versus $1143$ nm (experimental) with deviation $<0.1$ eV and no empirical fit.
• LSPR mode analysis: A high-energy shoulder near $773$ nm, sensitive to plate shape, was consistently observed in experiment and simulation, verifying quadrupolar contributions in LSPR response.
• Hyperbolic regime mapping: Only simulations with full dielectric anisotropy ($\varepsilon_{xx} \neq \varepsilon_{zz}$) reproduced the "envelope" plasmonic mode in the $550$–$880$ nm window, characterized by $\Re[\varepsilon_{xx}(\omega)] < 0$, $\Re[\varepsilon_{zz}(\omega)] > 0$. Isotropic models yielded only corner-confined or leaky dielectric modes.

6. Computation of Observable Spectral Quantities and Identification of Hyperbolic Behavior

After convergence, the dipole moments $\{\mathbf p_n\}$ are used to compute physically meaningful cross-sections:

• Absorption: $$C_{\rm abs} = \frac{4\pi k}{|E_0|^2} \sum_n \Im[ \mathbf p_n \cdot \mathbf E_{\rm loc}(\mathbf r_n)^* ]$$
• Scattering: $$C_{\rm sca} = \frac{ k^4 }{ 6\pi |E_0|^2 } \sum_n |\mathbf p_n|^2$$
• Extinction: $C_{\rm ext}=C_{\rm abs}+C_{\rm sca}$

Observed spectra (absorption, extinction) are obtained via orientation-averaged $C_{\rm abs}(\lambda)$ or normalized absorbance $A(\lambda) \propto C_{\rm abs}(\lambda)$. In the hyperbolic window, CSDDA++ uniquely resolves plasmonic "envelope" modes uniformly spanning the nanocrystal's surface, resulting from the interplay of propagating and evanescent fields in the regime where $\varepsilon_{xx}$ and $\varepsilon_{zz}$ have opposite signs (Parekh et al., 26 Dec 2025).

---

The CSDDA and its CSDDA++ extension enable rigorous, numerically stable simulation of anisotropic and hyperbolic nanoplasmonics, provided accurate dielectric tensors, sub-nanometer meshing, and specialized iterative solvers are available. These developments underlie key advances in understanding the optical and plasmonic behavior of anisotropic nanomaterials such as klockmannite CuSe.