Rigorous Coupled-Wave Analysis (RCWA)

• RCWA is a semi-analytical frequency-domain technique that employs Fourier expansion of material parameters and fields to solve Maxwell’s equations in periodic, layered structures.
• It transforms the boundary-value problem into coupled ordinary differential systems or matrix eigenproblems, enabling accurate prediction of diffraction, resonances, and optical responses.
• RCWA supports advanced applications such as photovoltaic optimization, plasmonic resonance analysis, and laser threshold determination through robust numerical and adaptive strategies.

Rigorous coupled-wave analysis (RCWA), also called the rigorous coupled-wave approach and often identified with the Fourier modal method, is a semi-analytical frequency-domain method for solving Maxwell boundary-value problems in structures that are periodic in one or two in-plane directions and layered along a stratification direction. Its central operation is the Fourier expansion of both material parameters and electromagnetic fields, which converts Maxwell’s equations into coupled ordinary differential systems or matrix eigenproblems along the layer-normal coordinate. In this form, RCWA is used to compute diffraction, reflection, transmission, absorptance, near fields, guided and leaky resonances, optical forces, and even lasing thresholds in periodic photonic structures (Civiletti et al., 2019, Gao et al., 2023, Song et al., 2018).

1. Fourier–modal formulation

RCWA begins from the frequency-domain Maxwell curl equations, with either an $e^{-i\omega t}$ or an $e^{+i\omega t}$ convention depending on the formulation. For a structure periodic in $x$ or in $x,y$, the relative permittivity and the electromagnetic fields are expanded in Fourier harmonics indexed by reciprocal-lattice vectors. In a 2D-periodic photonic crystal slab, for example, one writes

$$\varepsilon(\mathbf r)=\sum_G \varepsilon_G(z)e^{iG\cdot \rho},\qquad \mathbf E(\mathbf r)=\sum_G \mathbf E_G(z)e^{i(\mathbf k+G)\cdot \rho},\qquad \mathbf H(\mathbf r)=\sum_G \mathbf H_G(z)e^{i(\mathbf k+G)\cdot \rho},$$

with $\rho=(x,y)$ and $\mathbf k$ the conserved in-plane Bloch wavevector. In one-dimensional gratings, the same idea specializes to quasi-periodic Fourier series in the periodic coordinate, with Floquet shifts $\alpha_n=\alpha+2\pi n/L$ or $K_n=K_0+2\pi n/\Lambda$ (Song et al., 2018, Civiletti et al., 2019, Civiletti et al., 2020).

After substitution into Maxwell’s equations and projection onto each Fourier harmonic, the full vector PDE becomes a coupled system in the stratification variable. In different notational conventions this appears as

$\frac{d\Psi}{dz}=Q(\omega)\Psi$

or

$e^{+i\omega t}$0

where the state vector collects the retained Fourier coefficients of tangential electric and magnetic fields. The matrices are built from diagonal wavevector matrices such as $e^{+i\omega t}$1 and $e^{+i\omega t}$2 together with Toeplitz or convolution matrices of $e^{+i\omega t}$3, $e^{+i\omega t}$4, $e^{+i\omega t}$5, or $e^{+i\omega t}$6 (Song et al., 2018, Gao et al., 2023, Civiletti et al., 2018).

For one-dimensional gratings, polarization reductions are available. For s-polarization, the electric field component $e^{+i\omega t}$7 satisfies the scalar Helmholtz equation

$e^{+i\omega t}$8

whereas for p-polarization the formulation becomes

$e^{+i\omega t}$9

These scalar reductions support rigorous variational analyses of RCWA as a Galerkin method for the corresponding boundary-value problems (Civiletti et al., 2019, Civiletti et al., 2020).

2. Layer eigenproblems and scattering matrices

Within each uniform layer, the Fourier coefficients of the material profile are independent of $x$0, so the coupled-wave system has constant coefficients and can be diagonalized. This yields a modal basis of forward and backward eigenwaves, with field evolution expressed as exponentials of the corresponding propagation constants. In 2D formulations, the layer solution is commonly written in terms of eigendecompositions of products such as $x$1 or of a first-order system matrix; in 1D and 2D alike, the field in a layer is a superposition of modal amplitudes multiplied by $x$2 or $x$3 factors (Gao et al., 2023, Robertson et al., 2018).

Continuity of tangential $x$4 and $x$5 at every interface produces the global coupling between layers. Although transfer matrices can be written directly, stable RCWA implementations usually assemble scattering matrices for individual layers and then cascade them. In the standard notation,

$x$6

where $x$7 and $x$8 are incoming and outgoing plane-wave amplitude vectors over all retained diffraction channels. For multilayer stacks, layer matrices are combined by the Redheffer star product, a construction used both in conventional RCWA and in hybrid variants such as dipole-augmented Fourier-modal schemes (Song et al., 2018, Fradkin et al., 2018).

This scattering-matrix architecture is important numerically. Several formulations explicitly note that direct transfer-matrix multiplication becomes unstable when many layers or strongly evanescent channels are present, whereas stable S-matrix cascading avoids numerical overflow and permits treatment of thick gratings, metallic layers, and high-contrast structures (Ahmad et al., 2018, Hugonin et al., 2021). The same framework underlies software implementations such as "RETICOLO" (Hugonin et al., 2021).

3. Resonances, poles, and spectral singularities

In RCWA, resonant states are naturally encoded in the analytic structure of the scattering matrix. For passive periodic slabs, the total $x$9-matrix is a meromorphic function of complex frequency, and its poles $x,y$0 satisfy

$x,y$1

These poles correspond to the structure’s natural resonances. In a lossless passive system, the cited formulation states that $x,y$2 by causality (Song et al., 2018).

This pole picture extends RCWA from passive scattering to active devices. In photonic crystal surface-emitting lasers, gain can be introduced by replacing

$x,y$3

or, more realistically, by a Lorentz model analytically continued to complex $x,y$4. As the gain parameter increases, the relevant $x,y$5-matrix pole moves downward in the complex-frequency plane, and lasing threshold is reached when a pole crosses the real axis: $x,y$6 This gives a first-principles threshold criterion without a round-trip cavity argument and is particularly useful for surface-emitting laser systems with complex in-plane structures (Song et al., 2018).

Related eigenvalue viewpoints appear in other RCWA extensions. In the modified treatment of twisted bilayer photonic slabs, guided eigenmodes are identified from a round-trip operator in the gap-modal space,

$x,y$7

with resonances occurring as $x,y$8. In plasmonic grating structures with a $x,y$9-thin conductive channel, resonant absorption is resolved through modified boundary conditions incorporating a surface-current admittance, and RCWA recovers the corresponding far-field and near-field resonance structure (Xie et al., 15 May 2025, Lyaschuk et al., 2021).

4. Convergence theory, numerical errors, and near-field reconstruction

A major development in the mathematical analysis of RCWA is the identification of the method as a Galerkin scheme. For s-polarized gratings, Civiletti, Lakhtakia, and Monk show that RCWA with $$\varepsilon(\mathbf r)=\sum_G \varepsilon_G(z)e^{iG\cdot \rho},\qquad \mathbf E(\mathbf r)=\sum_G \mathbf E_G(z)e^{i(\mathbf k+G)\cdot \rho},\qquad \mathbf H(\mathbf r)=\sum_G \mathbf H_G(z)e^{i(\mathbf k+G)\cdot \rho},$$0 retained Fourier modes and a stairstep approximation in $$\varepsilon(\mathbf r)=\sum_G \varepsilon_G(z)e^{iG\cdot \rho},\qquad \mathbf E(\mathbf r)=\sum_G \mathbf E_G(z)e^{i(\mathbf k+G)\cdot \rho},\qquad \mathbf H(\mathbf r)=\sum_G \mathbf H_G(z)e^{i(\mathbf k+G)\cdot \rho},$$1 is exactly the variational discretization of a perturbed problem. Under non-trapping hypotheses, the analysis yields

$$\varepsilon(\mathbf r)=\sum_G \varepsilon_G(z)e^{iG\cdot \rho},\qquad \mathbf E(\mathbf r)=\sum_G \mathbf E_G(z)e^{i(\mathbf k+G)\cdot \rho},\qquad \mathbf H(\mathbf r)=\sum_G \mathbf H_G(z)e^{i(\mathbf k+G)\cdot \rho},$$2

that is, $$\varepsilon(\mathbf r)=\sum_G \varepsilon_G(z)e^{iG\cdot \rho},\qquad \mathbf E(\mathbf r)=\sum_G \mathbf E_G(z)e^{i(\mathbf k+G)\cdot \rho},\qquad \mathbf H(\mathbf r)=\sum_G \mathbf H_G(z)e^{i(\mathbf k+G)\cdot \rho},$$3 in $$\varepsilon(\mathbf r)=\sum_G \varepsilon_G(z)e^{iG\cdot \rho},\qquad \mathbf E(\mathbf r)=\sum_G \mathbf E_G(z)e^{i(\mathbf k+G)\cdot \rho},\qquad \mathbf H(\mathbf r)=\sum_G \mathbf H_G(z)e^{i(\mathbf k+G)\cdot \rho},$$4 and $$\varepsilon(\mathbf r)=\sum_G \varepsilon_G(z)e^{iG\cdot \rho},\qquad \mathbf E(\mathbf r)=\sum_G \mathbf E_G(z)e^{i(\mathbf k+G)\cdot \rho},\qquad \mathbf H(\mathbf r)=\sum_G \mathbf H_G(z)e^{i(\mathbf k+G)\cdot \rho},$$5 in $$\varepsilon(\mathbf r)=\sum_G \varepsilon_G(z)e^{iG\cdot \rho},\qquad \mathbf E(\mathbf r)=\sum_G \mathbf E_G(z)e^{i(\mathbf k+G)\cdot \rho},\qquad \mathbf H(\mathbf r)=\sum_G \mathbf H_G(z)e^{i(\mathbf k+G)\cdot \rho},$$6. For the stairstep approximation, the general bound is $$\varepsilon(\mathbf r)=\sum_G \varepsilon_G(z)e^{iG\cdot \rho},\qquad \mathbf E(\mathbf r)=\sum_G \mathbf E_G(z)e^{i(\mathbf k+G)\cdot \rho},\qquad \mathbf H(\mathbf r)=\sum_G \mathbf H_G(z)e^{i(\mathbf k+G)\cdot \rho},$$7 in $$\varepsilon(\mathbf r)=\sum_G \varepsilon_G(z)e^{iG\cdot \rho},\qquad \mathbf E(\mathbf r)=\sum_G \mathbf E_G(z)e^{i(\mathbf k+G)\cdot \rho},\qquad \mathbf H(\mathbf r)=\sum_G \mathbf H_G(z)e^{i(\mathbf k+G)\cdot \rho},$$8, improving to $$\varepsilon(\mathbf r)=\sum_G \varepsilon_G(z)e^{iG\cdot \rho},\qquad \mathbf E(\mathbf r)=\sum_G \mathbf E_G(z)e^{i(\mathbf k+G)\cdot \rho},\qquad \mathbf H(\mathbf r)=\sum_G \mathbf H_G(z)e^{i(\mathbf k+G)\cdot \rho},$$9 when $\rho=(x,y)$0 is smooth in $\rho=(x,y)$1. Their numerical experiments show measured $\rho=(x,y)$2-error $\rho=(x,y)$3 and slice-thickness decay approximately $\rho=(x,y)$4 (Civiletti et al., 2019).

For p-polarized gratings, the same authors derive a Galerkin interpretation for a perturbed problem with approximate permittivity and obtain a combined estimate

$\rho=(x,y)$5

Representative numerical examples show $\rho=(x,y)$6-convergence roughly $\rho=(x,y)$7 and $\rho=(x,y)$8-convergence $\rho=(x,y)$9, while the authors state that the numerical results suggest further work (Civiletti et al., 2020).

A recurring misconception is that RCWA’s known near-field difficulties imply comparable far-field inaccuracy. The cited studies separate these effects. Far-field quantities such as reflectance, transmittance, and absorptance typically converge rapidly with Fourier order, whereas direct Fourier reconstruction of discontinuous field components suffers Gibbs oscillations at material boundaries (Weismann et al., 2015, Robertson et al., 2018). The principal pathology is the reconstruction of the discontinuous normal electric field. Schuster et al. therefore propose a “reconstruct $\mathbf k$0, then divide” formulation,

$\mathbf k$1

and its 2D analogue based on a normal-vector field, so that only continuous combinations such as $\mathbf k$2 and $\mathbf k$3 are reconstructed in Fourier space (Weismann et al., 2015).

For nanowire solar cells, Robertson et al. introduce two additional remedies: the Continuous-Variable Formulation (CVF), which reconstructs continuous quantities $\mathbf k$4 and $\mathbf k$5, and a far-field-based rescaling

$\mathbf k$6

used to force agreement between layer-integrated near-field and far-field absorptance. In their test cell, far-field absorptance converges very rapidly, with $\mathbf k$7 giving $\mathbf k$8, whereas unmodified near fields at short wavelengths remain poorly converged even at $\mathbf k$9. With CVF and rescaling, $\alpha_n=\alpha+2\pi n/L$0 yields $\alpha_n=\alpha+2\pi n/L$1 total-generation error under AM1.5G and enables computational speedups between $\alpha_n=\alpha+2\pi n/L$2 and $\alpha_n=\alpha+2\pi n/L$3 times for spectrally integrated calculations (Robertson et al., 2018). The same separation between robust far-field momentum fluxes and unreliable boundary-local stress integrals motivates optical-force calculations based on momentum conservation rather than direct Maxwell-stress integration (Gao et al., 2023).

5. Extensions, hybridizations, and generalized formulations

Several extensions modify the standard Fourier-modal workflow to address geometries for which bare RCWA converges slowly or requires excessive slicing. One limitation arises in nanoparticle arrays: small, high-contrast inclusions generate local fields with large spatial gradients, so the Fourier series converges very slowly and may require $\alpha_n=\alpha+2\pi n/L$4 harmonics. Fradkin et al. therefore combine RCWA with the discrete dipole approximation (DDA), replacing the problematic particle layer by an effective dipole-sheet scattering matrix built from the lattice-corrected polarizability

$\alpha_n=\alpha+2\pi n/L$5

In the cited examples, the DDA-augmented scheme typically converges for $\alpha_n=\alpha+2\pi n/L$6 and captures localized surface plasmon resonances, lattice plasmon resonances, and their hybridization with guided modes (Fradkin et al., 2018).

A different limitation is longitudinal variation of the cross section. Standard RCWA treats $\alpha_n=\alpha+2\pi n/L$7-variation by staircasing into many thin layers. "VarRCWA: An Adaptive High-Order Rigorous Coupled Wave Analysis Method" (Zhu et al., 2022) reinterprets constant-cross-section RCWA as a zeroth-order approximation and derives a Born-series-style perturbative expansion in the variation of $\alpha_n=\alpha+2\pi n/L$8 and $\alpha_n=\alpha+2\pi n/L$9. Using a reference cross section and adaptive subdivision, VarRCWA reduces the number of required sections while meeting a user-specified tolerance. In the reported benchmarks, speedups for $K_n=K_0+2\pi n/\Lambda$0 range from $K_n=K_0+2\pi n/\Lambda$1 up to $K_n=K_0+2\pi n/\Lambda$2, and a linear taper is solved with four adaptive sections and first-order expansion in under $K_n=K_0+2\pi n/\Lambda$3 of the time needed by conventional RCWA with $K_n=K_0+2\pi n/\Lambda$4 slices (Zhu et al., 2022).

A more radical approach appears in "A rigorous coupled-wave analysis of birefringent holographic gratings with periodically-modulated dielectric tensor along an in-plane direction and tensor variations in the thickness direction" (Mansuripur et al., 19 Sep 2025). Mansuripur and Jakobsen avoid slicing entirely by expanding the dielectric tensor in both $K_n=K_0+2\pi n/\Lambda$5 and $K_n=K_0+2\pi n/\Lambda$6, producing one augmented $K_n=K_0+2\pi n/\Lambda$7 eigenproblem with $K_n=K_0+2\pi n/\Lambda$8. This unsliced formulation introduces a $K_n=K_0+2\pi n/\Lambda$9-fold degeneracy for each physical diffraction order, but the cited numerical results report stable solutions, energy conservation $\frac{d\Psi}{dz}=Q(\omega)\Psi$0 to better than $\frac{d\Psi}{dz}=Q(\omega)\Psi$1, and matrix sizes up to $\frac{d\Psi}{dz}=Q(\omega)\Psi$2 (Mansuripur et al., 19 Sep 2025).

Other generalized formulations incorporate singular sheets or modified basis functions. Lyaschuk et al. embed a $\frac{d\Psi}{dz}=Q(\omega)\Psi$3-thin two-dimensional electron gas in RCWA by replacing the usual tangential-field continuity with jump conditions involving a boundary admittance $\frac{d\Psi}{dz}=Q(\omega)\Psi$4, which enables percent-level accurate modeling of THz plasmonic grating structures in less than $\frac{d\Psi}{dz}=Q(\omega)\Psi$5 s on a desktop for $\frac{d\Psi}{dz}=Q(\omega)\Psi$6 (Lyaschuk et al., 2021). Xie et al. replace evanescent gap harmonics in twisted bilayer photonic slabs by flux-carrying basis functions and obtain a modified eigenmode analysis, together with a five-layer uniform-slab approximation of accuracy around $\frac{d\Psi}{dz}=Q(\omega)\Psi$7 (Xie et al., 15 May 2025). These developments suggest that RCWA is best understood not as a fixed algorithmic recipe but as a family of Fourier-domain scattering and eigenvalue formulations.

6. Applications, software, and inverse design

RCWA is widely used wherever periodicity and layered structure coincide. In thin-film photovoltaics, it has been coupled to differential evolution for geometry and bandgap optimization of three-dimensional corrugated tandem solar cells. In the cited optimization study, the population size is $\frac{d\Psi}{dz}=Q(\omega)\Psi$8, the crossover is $\frac{d\Psi}{dz}=Q(\omega)\Psi$9, the differential weight is $e^{+i\omega t}$00, and $e^{+i\omega t}$01 generations were found sufficient for near-convergence (Civiletti et al., 2018). In a related study of a two-dimensional periodically corrugated metallic backreflector, RCWA is used to separate total absorptance from useful semiconductor absorptance and to correlate spectral peaks with surface-plasmon-polariton waves and waveguide modes, emphasizing that not every resonance enhancing total absorption is useful for the photovoltaic junctions (Ahmad et al., 2018).

In subwavelength coronagraphy, RCWA is used to design and reverse-engineer annular groove phase masks. The cited AGPM study reports that coronagraphic performance is very sensitive to small errors in etch depth and grating profile; after post-fabrication re-etch optimization, starlight rejection up to $e^{+i\omega t}$02 is demonstrated in a broadband L filter, corresponding to a raw contrast of about $e^{+i\omega t}$03 at two resolution elements from the star for a perfect input wave front on a circular, unobstructed aperture (Catalan et al., 2016). In graded chiral TiO$e^{+i\omega t}$04 thin films, RCWA combined with Bruggeman homogenization predicts filtering frequency and polarization selectivity, with a reported theoretical Bragg wavelength of $e^{+i\omega t}$05 nm versus an experimental $e^{+i\omega t}$06 nm and maximum selectivity of about $e^{+i\omega t}$07 for the stated deposition geometry (Babaei et al., 2010).

For optical-force calculations on periodic nanophotonic actuators, RCWA is presented as a practical alternative to FDTD. In the benchmark asymmetric Si–SiO$e^{+i\omega t}$08 dimer, RCWA converges at $e^{+i\omega t}$09 Fourier harmonics, gives $e^{+i\omega t}$10 difference from FDTD, and costs about $e^{+i\omega t}$11 s per single-wavelength evaluation on an Nvidia P100, versus about $e^{+i\omega t}$12 s for FDTD on the cited HPC configuration. Across a geometry sweep of $e^{+i\omega t}$13 gap positions, RCWA requires about $e^{+i\omega t}$14 min versus about $e^{+i\omega t}$15 h for FDTD, and on consumer laptops with CUDA GPUs the single-precision RCWA runtime is about $e^{+i\omega t}$16 s per simulation while remaining $e^{+i\omega t}$17 faster than FDTD (Gao et al., 2023).

The software ecosystem around RCWA includes "RETICOLO software for grating analysis" (Hugonin et al., 2021), a MATLAB implementation for 1D classical and conical diffraction and 2D crossed gratings, with Bloch-mode and field-visualization toolboxes and, in version V9 and later, a toolbox for arbitrarily anisotropic multilayer thin films. More recently, surrogate models have been built directly on RCWA outputs. In "Physics-constrained neural networks for surrogate modeling of lossless periodic structures" (Prehn et al., 26 Jun 2026), Prehn and Jung treat the block Jones matrix predicted by RCWA as a point on the complex Stiefel manifold, enforce energy conservation by differentiable symmetric orthogonalization, and report test-set mean-squared errors on the order of $e^{+i\omega t}$18 to $e^{+i\omega t}$19, energy-conservation error of order $e^{+i\omega t}$20, and about $e^{+i\omega t}$21 Jones-matrix predictions per iteration in approximately $e^{+i\omega t}$22 s on an NVIDIA A100 GPU. A plausible implication is that RCWA now functions both as a direct solver and as a generator of physically constrained reduced models for large-scale inverse design (Prehn et al., 26 Jun 2026).