A generalized perturbative approach for the computation of nonlinear scattering problems

1. Problem Statement

The simulation of light scattering by nonlinear optical materials is computationally challenging because the constitutive relations are nonlinear, leading to a system of fully coupled partial differential equations (PDEs).

• Current State of the Art: The standard "rigorous" approach solves these coupled nonlinear systems using iterative algorithms (e.g., Picard or Newton-Raphson) combined with linear solvers like the Finite Element Method (FEM).
• Limitations: This iterative process is time-consuming, often requiring ten or more iterations to converge. In computationally demanding scenarios, such as brute-force optimization or 3D simulations, this approach becomes prohibitively expensive.
• Existing Approximations: Previous perturbative methods often rely on the "undepleted-pump approximation," which retains only the first two equations of the expansion. While fast, this neglects pump depletion and higher-order effects, limiting accuracy for strong nonlinearities or long interaction lengths.

2. Methodology

The authors propose a generalized perturbative approach that decouples the nonlinear system into a sequence of linear equations solvable sequentially, eliminating the need for iterative loops.

Theoretical Framework

• Field Expansion: The total electric field $\mathbf{E}(\mathbf{r}, t)$ is expanded as a power series in a perturbation parameter $\eta$:
  $$\mathbf{E}_p = \sum_{j=1}^{n} \eta^j \mathbf{E}^{(j)}_p$$
  where $p$ denotes the harmonic order and $j$ denotes the perturbative order.
• Decoupling: By substituting this expansion into the governing Maxwell equations (considering up to third-order nonlinearities), the authors derive a set of polynomial equations. Setting coefficients of each power of $\eta$ to zero yields a triangular (cascade) system:
  • Order $j=1$: A standard linear Maxwell equation driven by the source.
  • Order $j > 1$: Linear equations where the source term is composed of previously computed lower-order fields ($\mathbf{E}^{(l)}$ where $l < j$).
• Key Properties: The authors prove two critical properties that drastically reduce the number of equations:
  1. Harmonic Limit: The field component $\mathbf{E}^{(j)}_p$ is zero if the harmonic order $|p| > j$.
  2. Parity Rule: The field $\mathbf{E}^{(j)}_p$ is zero if $p$ and $j$ have different parity (i.e., $p - j$ is odd).
     These properties ensure that at any perturbative order $j$, only specific harmonics need to be calculated, creating a sparse, sequential system.

Numerical Implementation

• Solver: The linear equations at each step are solved using the Finite Element Method (FEM) via the open-source software Gmsh and GetDP.
• Geometry: The method is applied to 2D structures (infinite in $y$ and $z$) illuminated by plane waves, reducing the problem to 1D formulations with appropriate boundary conditions.
• Validation: The results are compared against a rigorous iterative FEM solution (previously published by the same authors) and validated using an energy conservation rule (sum of reflection, transmission, and losses).

3. Key Contributions

1. Generalization: The method extends perturbative analysis beyond the undepleted-pump approximation to arbitrary orders ($n$-th order), allowing for the inclusion of pump depletion and higher-order harmonic generation.
2. Algorithmic Efficiency: By transforming a coupled nonlinear problem into a cascade of linear problems, the method removes the need for iterative convergence loops.
3. Analytical Proofs: The paper provides rigorous proofs for the independence of the solution from the perturbation parameter $\eta$ (allowing $\eta=1$) and the vanishing conditions for specific field components, which optimize the computational load.
4. Broad Applicability: The framework handles arbitrary polarization, conical incidence, general susceptibility tensors ($\chi^{(2)}$ and $\chi^{(3)}$), and mixed nonlinear effects (e.g., simultaneous second- and third-order processes).

4. Results

The authors validated the method using two nonlinear materials: Potassium Titanyl Phosphate (KTP) and Lithium Niobate (LiNbO$_3$).

Case Study 1: Second-Harmonic Generation (SHG) in KTP

• Accuracy:
  • 2nd Order (Undepleted Pump): Showed significant errors (up to 25% energy balance error) due to neglecting pump depletion.
  • 4th Order: Reduced error to ~1.9%.
  • 8th Order: Achieved relative errors as low as 0.31% for a 30 $\mu$m crystal with 38% conversion efficiency.
• Convergence: The method converges as long as the incident amplitude is below a critical radius defined by $E_{0,max} \propto \lambda_0 / (L \chi^{(2)})$. Beyond this radius, the series diverges, similar to the rigorous method's failure to converge.
• Computation Time:
  • The perturbative approach showed constant computation time regardless of incident amplitude (unlike the rigorous method, which scales exponentially with amplitude due to iteration count).
  • At high amplitudes ($3 \times 10^9$ V/m), the 4th-order perturbative method was 4 times faster than the rigorous method while maintaining ~0.1% error.

Case Study 2: Coupled 2nd and 3rd Order Effects in LiNbO$_3$

• Complexity: Simulated a scenario involving Second-Harmonic Generation (SHG), Third-Harmonic Generation (THG), and sum-frequency generation simultaneously.
• Performance: The 5th-order perturbative method achieved an energy balance error of $6 \times 10^{-3}\%$ with a computation time of 5 seconds, compared to 16 seconds for the rigorous method.
• Conclusion: Even for complex, multi-harmonic interactions, the perturbative method maintains high accuracy with a substantial speedup.

5. Significance and Outlook

• Computational Efficiency: The primary significance is the drastic reduction in computational cost. Since the most expensive part of FEM is solving the matrix system, avoiding repeated iterations makes this method ideal for large-scale 2D and 3D problems.
• Optimization Enabler: The speedup makes brute-force optimization of nonlinear photonic devices (e.g., designing gratings or metasurfaces for maximum harmonic conversion) feasible, which was previously too costly with rigorous iterative solvers.
• Versatility: The approach is not limited to 1D slabs; the authors argue it is directly applicable to complex geometries where the gain in efficiency would be even more pronounced.
• Limitations: The method is valid within the radius of convergence. For extremely high intensities or very long interaction lengths where the rigorous method fails to converge, this perturbative approach will also diverge. However, for the vast majority of practical nonlinear optical regimes, it offers a superior balance of speed and accuracy.

In summary, this paper presents a robust, mathematically rigorous, and highly efficient alternative to iterative nonlinear solvers, enabling the rapid simulation and optimization of complex nonlinear scattering problems.