1 Introduction↩︎

Photonics engineering leverages the geometric design of materials to control electromagnetic (EM) responses. Periodic structures are widely employed in photonic systems, with applications such as beam splitters [1], beam steering [2], [3], guided-mode resonance filters [4], couplers in integrated optics [5], [6], imaging [7] and waveguide combiners [8]. The design of these structures relies on accurate and efficient EM simulations, which are essential for both forward modeling and inverse design workflows.

Rigorous coupled-wave analysis (RCWA) [9]–[11] provides a semi-analytic solution for layered periodic structures and is widely used in photonic device design. However, its computational cost can become prohibitive for large-scale optimization and inverse design over high-dimensional design spaces. Neural networks (NNs) offer a promising alternative as surrogate EM solvers, learning from data to accelerate simulations by orders of magnitude [12], [13]. Once trained, they enable rapid inference and flexible inverse design through readily accessible gradients [14], [15].

To enforce consistency with the underlying physical laws, several physics-informed and physics-augmented neural-network approaches have been proposed for EM modeling [16], [17]. Physics-informed neural networks typically incorporate governing equations through residual terms in the loss function [18], while physics-augmented surrogate models embed domain knowledge through tailored architectures or hybrid solvers such as WavYNet [19]. These methods improve physical fidelity but often rely on soft constraints and weighted loss terms.

Here, we introduce a PCNN that enforces energy conservation as a hard output constraint. For lossless periodic structures, energy conservation restricts Jones matrices describing the amplitude, phase, and polarization of scattered far fields to a Stiefel manifold [20]. By projecting the network output directly onto the manifold, energy conservation is satisfied by construction. We evaluate the resulting surrogate against an unconstrained NN baseline and demonstrate its use in a differentiable waveguide-combiner simulation. The PCNN approach provides a scalable framework for waveguide engineering and offers a practical route toward high-performance augmented reality (AR) display systems. More broadly, it establishes a general methodology for the design of lossless layered periodic structures.

2 Energy-conserving Jones matrix representation↩︎

For a lossless periodic structure, the sum of diffraction efficiencies over all \(M\) propagating reflected and transmitted orders satisfies \[\label{eq:DE} \sum_{m=1}^M \eta_m = 1.\tag{1}\] Here, \(\eta_m\) is the diffraction efficiency of the non-evanescent order \(m\), defined as the ratio between the flux of the diffracted Poynting vector and the flux of the incident Poynting vector [21]. In direct RCWA calculations, this condition must be satisfied, with Moharam et al. recommending a numerical accuracy on the order of \(10^{-10}\) [10]. Strict energy conservation is essential, as even a single RCWA prediction in which the total diffraction efficiency exceeds unity can compromise an entire optical simulation.

For an incident electric field with Jones vector \(e\), the field diffracted into the \(m\)th order is \[e_m = J_m e = \begin{pmatrix} j^{\mathrm{TE},\mathrm{TE}}_m & j^{\mathrm{TM},\mathrm{TE}}_m \\ j^{\mathrm{TE},\mathrm{TM}}_m & j^{\mathrm{TM},\mathrm{TM}}_m \end{pmatrix} \begin{pmatrix} e_{\mathrm{TE}} \\ e_{\mathrm{TM}} \end{pmatrix} = \begin{pmatrix} a_m & b_m \\ c_m & d_m \end{pmatrix}e.\] The Jones matrix \(J_m \in \mathbb{C}^{2 \times 2}\) is expressed in transverse electric (TE) and transverse magnetic (TM) polarization bases, defined with respect to planes of incidence and diffraction [21]. For clarity, we adopt the shorthand notation \(a_m, b_m, c_m, d_m\) for matrix elements. The condition in Eq. (1 ) is equivalent to \[\sum_{m=1}^M \|J_m e\|^2 = e^\dagger\left(\sum_m J_m^\dagger J_m\right)e \overset{\displaystyle !}{=} \|e\|^2,\] which must hold for arbitrary \(e\), and hence \[\label{eq:constraints} \sum_m J_m^\dagger J_m = \begin{pmatrix} \sum_m |a_m|^2 + |c_m|^2 & \sum_m a_m^* b_m + c_m^* d_m \\ \sum_m b_m^* a_m + d_m^* c_m & \sum_m |b_m|^2 + |d_m|^2 \end{pmatrix} = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix}=\mathbb{I}_2.\tag{2}\] A convenient construction satisfying these constraints is obtained by forming vectors of Jones matrix elements, \[\label{eq:uv95def} \begin{align} u &= (a_1,a_2,\ldots,a_M,c_1,c_2,\ldots,c_M), \\ v &= (b_1,b_2,\ldots,b_M,d_1,d_2,\ldots,d_M). \end{align}\tag{3}\] Enforcing orthonormality of \(u\) and \(v\) ensures that Eqs. (1 ) and (2 ) are satisfied: \[\label{eq:constraints-uv} \sum_m J_m^\dagger J_m = \begin{pmatrix} \|u\|^2 & u^\dagger v \\ v^\dagger u & \|v\|^2 \end{pmatrix} = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix}.\tag{4}\]

3 Physics-constrained Neural Network↩︎

The surrogate model aims to learn the nonlinear mapping from RCWA input parameters to the corresponding Jones matrices, \[\label{eq:mapping} (p_1,p_2,\ldots,p_D,k_x,k_y) \mapsto \{J_m(p_1,p_2,\ldots,p_D,k_x,k_y)\}_{m=1}^M.\tag{5}\] Here, \(p_d\) denotes a geometric design parameter, such as grating duty cycle, \(D\) is the dimensionality of the design space, and \((k_x,k_y)\) are the transverse wavevector components of the incident light. Each input configuration is mapped to Jones matrices for \(M\) diffraction orders retained from the RCWA simulation. The layered periodic structure is invariant along the longitudinal \(z\)-direction within each layer, and periodic in the transverse plane. We propose a PCNN architecture that first predicts an unconstrained complex matrix \(V=[u\;v]\) and then projects it onto the Stiefel manifold \(\mathrm{St}(2,2M,\mathbb{C})\), the set of orthonormal 2-frames in \(\mathbb{C}^{2M}\) [20]. At the output of the PCNN, the resulting orthonormal vectors \(\tilde{V}=[\tilde{u}\;\tilde{v}]\) are reshaped into Jones matrices.

The projection is performed using Löwdin symmetric orthogonalization [22], [23]. Writing the thin singular-value decomposition (SVD) [24] of \(V\) as \[V = U\Sigma W^\dagger,\] the orthonormalized output is \[\label{eq:NN95svd95out} \tilde{V} = U W^\dagger.\tag{6}\] This projection treats the two output vectors symmetrically, in contrast to sequential Gram–Schmidt orthogonalization [24]. For further details, a derivation of Eq. (6 ) is provided (Section S1, Supplement 1). In summary, the network learns within the manifold of admissible solutions, guaranteeing that the final outputs satisfy energy conservation.

In this work, the models are implemented as fully connected multilayer perceptrons (Fig. 1) with rectified linear unit activations [25]. The projection is independent of the particular network architecture or input parameterization and provides a general approach for incorporating energy conservation into RCWA surrogate models.

4 Dataset Generation and Model Training↩︎

In this work, we investigate the surrogate modeling of one-dimensional all-dielectric binary gratings composed of single-crystal silicon and silicon dioxide. These materials were modeled as lossless owing to their low optical loss across the visible regime [26], Sell2016_visible? and were selected for their compatibility with standard semiconductor fabrication processes. RCWA simulations were performed using Reticolo [21] to generate datasets for the two design spaces. The first design space, denoted \(D=2\), comprises a grating parameterized by the duty cycle (\(DC\)) and the thickness of a \(\mathrm{SiO}_2\) underlayer (\(h_{\mathrm{ul}}\)). The \(D=4\) design extends the structure to a multilevel configuration, in which an additional binary grating is introduced on top of the first. Its geometric parameters include the duty cycles of the two gratings, \(h_{\mathrm{ul}}\) and the lateral displacement of the upper grating. Datasets comprising 25,000 and 100,000 samples were generated for the \(D=2\) and \(D=4\) design spaces, respectively, with samples spanning a range of incident wavevectors (Section S2 of Supplement S1). Representative gratings from the datasets are shown in Fig. 2, and these datasets are subsequently used in the waveguide optimization demonstration.

The models are trained using supervised learning (training details are provided in Section S3 of Supplement 1), and performance is evaluated by comparing the predicted Jones matrices \(\hat{J}_m\) to the simulated RCWA matrices \(J_m\) using

\[\label{eq:NN95loss} \mathrm{MSE} = \frac{1}{4NM} \sum_{n=1}^{N} \sum_{m=1}^{M} \sum_{i,j \in \{\mathrm{TE},\mathrm{TM}\}} \left| \hat{J}_m^{n,ij}-J_m^{n,ij} \right|^2 .\tag{7}\] The mean-squared error (MSE) captures deviations in both amplitude and phase via the real and imaginary components across all \(N\) samples and \(M\) diffraction orders. Unlike related approaches that rely on multiple weighted loss terms with potentially competing gradients [19], [27], the PCNN employs a single loss term. This avoids gradient-balancing issues and eliminates the need to tune loss-term weights.

To test the impact of the constraint, the PCNN models were compared with unconstrained NN models. Both were trained with the same supervised MSE objective and differ only by the hard projection layer. The NN uses the same input features and hidden-layer form, but omits the projection step and directly outputs Jones matrices. To ensure a fair comparison, the NN and PCNN were independently selected using the same validation-based model-selection protocol. For both models, the number of hidden layers \(S\) and number of nodes \(h_S\) in each layer were swept over an identical search space, and this procedure was repeated for different fractions of the training dataset (Section S3, Supplement 1). The models with the lowest validation error were selected and used for evaluation on test sets.

5 Impact of the Energy Constraint↩︎

Here, we evaluate the models trained on the complete training dataset and report on reduced training-data regimes in Section S3 of Supplement 1. Across all three wavelengths and both design spaces, the PCNN achieved a lower test MSE than the NN, indicating that the energy-conservation constraint acts as an inductive bias. Although both methods achieve similar test MSE values on the order of \(10^{-5}\) to \(10^{-6}\) (Section S3, Supplement 1), substantial differences in physical consistency remain. To quantify this, the energy conservation for a structure is measured by the deviation of the predicted response from Eq. (2 ), \[\label{eq:energy95error} \Delta_E = \left\| \sum_{m=1}^{M} \hat{J}_m^{\dagger} \hat{J}_m - \mathbb{I}_2 \right\|_{\infty} = \; \max_{i,j} \left| \left( \sum_{m=1}^{M} \hat{J}_m^{\dagger} \hat{J}_m - \mathbb{I}_2\right)_{ij} \right| .\tag{8}\] Here, \(\Delta_E=0\) corresponds to exact conservation of total diffracted energy for arbitrary incident polarization. Because the PCNN output is projected onto the Stiefel manifold, Eq. (8 ) is satisfied up to numerical precision, resulting in test-set energy-conservation errors of order \(\Delta_E\sim10^{-6}\). In comparison, the unconstrained NN yields errors of order \(\Delta_E\sim10^{-2}\). This discrepancy highlights the fact that low prediction error does not necessarily imply physical consistency. For lossless structures, even small deviations from the admissible manifold can produce significant violations of energy conservation, limiting the reliability of unconstrained surrogate models.

The accuracy of the surrogate model is governed by several factors, including the density at which the design space is sampled, the surrogate model itself, and the underlying physics within the parameter space. The converged RCWA solution, obtained by solving Maxwell’s equations, can exhibit significantly different sensitivities across regions of the design space and wavevector coordinates \((p_1,\ldots,p_D,k_x,k_y)\). In sensitive regions, small parameter variations can produce large changes in the Jones matrices. For the grating structures considered here, the duty cycle is a sensitive parameter as it directly modifies the Fourier coefficients of the periodic permittivity distribution in the RCWA formulation [9]–[11], reshaping the Fourier spectrum of the layer. We therefore perform one-dimensional duty-cycle sweeps on additional test sets while keeping the remaining geometric and wavevector parameters fixed. With reliability in mind and to assess worst-case behavior beyond MSE, we evaluate the infinity norm error for a structure as \[\label{eq:J95infty95error} \mathcal{L}_{\infty} = \max_{m,i,j} \left| \hat{J}^{ij}_m - J^{ij}_m \right|.\tag{9}\]

The \(\mathcal{L}_\infty\) error corresponds to the maximum absolute element-wise error across all Jones matrix elements and diffraction orders. Fig. 3 shows the corresponding behavior for a sweep in the \(D=2\) space. The upper panels compare the \(\mathcal{L}_\infty\) and energy conservation errors as a function of duty cycle. Although the PCNN and NN achieve comparable Jones-matrix errors, the NN exhibits substantial violations of energy conservation, resulting in nonphysical loss and gain. The lower panels show the evolution of the zeroth reflected-order Jones matrix as a function of duty cycle for \(\lambda=0.53~\mu\mathrm{m}\). Both models accurately reproduce the RCWA trajectory, while the PCNN predictions remain restricted to the manifold.

Data-efficiency comparisons further show that the constraint usually preserves or improves the predictive accuracy once the design spaces in this work are sufficiently sampled (Section S3, Supplement 1). Together, these results indicate that the energy-conservation constraint provides an effective inductive bias and improves the physical reliability of the learned model. Having evaluated worst-case \(\Delta_E\) and \(\mathcal{L}_\infty\) metrics along sensitive trajectories in the design space, while also observing strong overall MSE performance, we conclude that the PCNNs are capable of serving as accurate and reliable surrogates. This enables their application to the inverse-design tasks considered in the next section.

6 Waveguide Optimization Demonstration↩︎

To illustrate the capabilities of the proposed PCNN framework, we demonstrate a ray-tracing-based optimization of the outcoupling (OC) region in a full-color diffractive waveguide combiner for AR glasses. Gratings are located on the upper surface of a planar single-crystal silicon waveguide and form an incoupling (IC) region and an OC region, resulting in one-dimensional exit-pupil expansion [8], [28]. The IC grating design is fixed, whereas the OC region is discretized into a two-dimensional grid of \(40{,}000\) grating subregions. The OC region is modeled using a locally periodic approximation (LPA) [29], assuming slowly varying grating geometries across neighboring subregions. Each OC subregion is assigned a periodic grating with parameter vector \(\mathbf{p}\in\mathbb{R}^D\) and the collection of parameter vectors across all subregions defines a large design space \(P\). To mimic manufacturable designs, selected geometric parameters are constrained to remain uniform across the OC region. Fig. 4 sketches the waveguide layout and eyebox region, where rays are collected to form a viewable image.

The forward model combines geometric ray tracing with PCNN evaluations in a reduced-order treatment of the dominant waveguide propagation paths (Section S5, Supplement 1). It starts by tracing rays originating from a microdisplay, spanning a 20° diagonal field of view (FOV). These rays enter the IC and undergo total internal reflection within the waveguide. When rays reach the discretized OC region, ray–grating interactions are evaluated using PCNNs. Relevant Jones matrices are applied sequentially to accumulate the transmitted and reflected fields, yielding intensity distributions across the eyebox. The polarization-averaged intensities are binned by wavevector direction and form the basis of the optimization objective, denoted \(C(P)\). The objective uses approximate CIE 1931 color-matching functions [30] to convert the simulated three-wavelength intensities into XYZ tristimulus values from which CIE-style \(L^*a^*b^*\) chromatic coordinates are computed relative to a D65 white point. Nonzero \(a^*\) and \(b^*\) values are penalized, together with luminance and chromaticity variations across angular eyebox bins (Section S5 of Supplement 1).

In conventional EM solvers, evaluating the optical response for each ray–grating interaction within an iterative loop can incur substantial computational cost. Here, the PCNNs enable efficient simulation of large ray ensembles while preserving energy conservation (time and memory complexity comparisons are provided in Section S4 of Supplement 1). Forward and gradient computations are implemented within the JAX framework [31], enabling reverse-mode (adjoint) automatic differentiation, just-in-time compilation, and efficient vectorization over rays. The purpose of the demonstration is to illustrate the scale of differentiable optimization enabled by the surrogate framework rather than to provide a production-ready waveguide design. The optimization launches 40,015 rays per wavelength into the waveguide. Each iteration evaluates approximately \(5\times10^6\) Jones matrices, recursively propagates electric fields, computes eyebox intensity distributions, evaluates adjoint gradients, and updates the design parameters. When JIT-compiled on an A100 GPU, the complete iteration requires approximately \(0.5\) seconds. This computational performance enables repeated optimization studies under varying initial conditions and objective functions \(C(P)\).

Fig. 5 shows the results of an optimization for the \(D=4\) design space. The final lower and upper grating duty-cycle maps exhibit spatially varying structures across the OC region, reflecting the field redistribution required to homogenize the eyebox. Polarization-averaged eyebox renderings are shown in display sRGB for visualization purposes (Section S5 of Supplement 1). The optimized design produces a more uniform full-color eyebox, which is consistent with the optimization objective. Visualization 1 (Supplement 2) provides an animation of the optimization in Fig. 5, along with additional optimization runs. The \(D=4\) design space yielded more uniform eyebox distributions than the \(D=2\) design space while maintaining a comparable optimization runtime.

7 Conclusion↩︎

We presented a physics-constrained neural network for the surrogate modeling of RCWA outputs in the form of Jones matrices. By projecting the network output onto the complex Stiefel manifold associated with lossless scattering, the model enforces energy conservation by construction while remaining fully differentiable. Our results demonstrate that unconstrained surrogate models can exhibit substantial violations of energy conservation despite achieving low prediction error, motivating the direct incorporation of the manifold constraint. Across both grating design spaces examined for AR waveguide applications, the PCNN maintained or improved predictive accuracy while enforcing physical consistency.

The surrogate is particularly attractive for applications in which only the far-field scattering response is required. In these cases, Jones matrices provide a compact representation that retains the relevant amplitude, phase, and polarization information without requiring the full electromagnetic field distribution. The PCNN allows the optical response of periodic structures to be evaluated rapidly within larger simulation and optimization frameworks. We demonstrated this through a gradient-based optimization of a full-color diffractive waveguide combiner, involving millions of Jones-matrix evaluations per iteration.

More generally, the proposed framework provides a physically constrained surrogate-modeling approach for RCWA-based design of lossless layered periodic structures.

Acknowledgment↩︎

This project was made possible by the DLR Quantum Computing Initiative and the Federal Ministry for Research, Technology and Space; qci.dlr.de/projects/qcoptsens/

Disclosures↩︎

The authors declare they have no conflicts of interest.

Data availability↩︎

The RCWA-generated datasets used for training and evaluating the surrogate models are not publicly available at this time but may be obtained from the authors upon reasonable request.

Code Availability↩︎

The surrogate-model implementations are not publicly available at this time but may be obtained from the authors upon reasonable request.

Supplemental document↩︎

See Supplement 1 for supporting content and Supplement 2 for Visualization 1, both are provided as ancillary files with this arXiv submission.

References↩︎