Numerical Estimation of Bending in Holographic Volume Gratings by Means of RCWA and Deep Learning

Abstract

In this paper, we introduce a novel approach to model bending phenomena on holographic volume gratings based on Rigorous Coupled Wave Analysis (RCWA), in which the bending as a phase in the dielectric permittivity expansion is introduced, and the Shooting Method (SM) is employed to solve the resulting system of equations. Further validation of our model is conducted by comparing its predictions to those obtained from reference Finite-Difference Time-Domain (FDTD) simulations and Coupled Wave Theory (CWT, referring to Kubota’s model that includes the bending phenomenon). Furthermore, we propose a methodology for estimating the bending from the diffraction efficiency curves in transmission volume gratings based on deep learning models, with a subsequent study of their accuracy and applicability.

1. Introduction

Holographic gratings have emerged as versatile and indispensable tools in optics and photonics, finding applications in fields as diverse as spectroscopy, telecommunications, and laser technology. They occupy a central position in various scientific, industrial, and technological domains because they can disperse, reflect, or diffract light in precise and controlled ways [1,2,3,4].

Often, imperfections in the recording or fabrication process may lead to an unwanted bending of the structure of diffraction gratings, which inevitably impacts its optical properties in uncontrolled ways. A non-uniform deformation of its fringes characterises this phenomenon, and it is mainly associated with a pre-stress in the recording material. On the other hand, many researchers have focused on bending as a crucial element for holographic sensors [5,6]. Zhou et al. [7] considered the bending deformation of a hydrogel actuator to be an essential deformation parameter, and Yetisen et al. have corroborated the application of holographic sensors as versatile optical sensing platform to characterise the deformation of the polymer due to its narrow diffraction spectrum, high visualisation, and stability [8,9]

The phenomenon of bending within holographic gratings has been understudied despite its frequent occurrence. The leading theory previously developed on the subject has been carried out by T. Kubota, which is based on Kogelnik’s coupled wave theory with an attenuated grating model [10,11]. This method focuses solely on the substrate itself and fails to account for air–medium discontinuities, thereby presenting limitations in its applicability. Indeed, over the years, Rigorous Coupled Wave Theory (RCWT) has been proven to be a powerful and versatile numerical method in the context of diffractive optics and photonics. This technique, initially developed by Moharam and Gaylord [12,13,14,15], has significantly advanced our understanding of light diffraction and its interaction with periodic structures. RCWT is founded on the rigorous solution of Maxwell’s equations and has been widely adopted in the study of diffraction gratings, photonic crystals, and other periodic optical structures [16]. Over the years, RCWT has evolved and been extended to address complex and practical problems, including incorporating anisotropic and inhomogeneous materials [17]. Its continued development and application have cemented RCWT as a reference of numerical techniques for diffractive optical analysis and design.

Accurately determining the bending of holographic grating fringes measured in laboratory settings implies using cost-consuming advanced techniques such as Scanning Electron Microscopy, Ellipsometry, or Interferometry. Due to the complex nature and the non-derivability of the formalism developed for the analysis of bending, an analytical approach to characterise these deformations from diffraction efficiency patterns is not feasible. Therefore, this research aims to set up a novel scheme based on deep learning solutions for estimating the parameters that define the bending in holographic gratings from the diffraction efficiency curves obtained numerically, addressing these limitations and providing a more accessible solution. Neural network-based solutions, particularly due to their nonlinear activation functions, are widely known precisely at regression problems in cases where analytical, linear solutions are not feasible. Thus, the most straightforward and viable solution is developing a computational regression model based on these technologies  [18,19,20]. For feeding the deep learning engine, we proposed a novel RCWA-SM scheme that considers bending as a phase in dielectric permittivity expansion and uses the shooting method to solve the system of equations. A numerical SF-FDTD scheme was also considered for validating the RCWA-SM implementation. It is worth noting that RCWA-SM was found to be the best candidate for use in the deep learning engine due to its performance compared to the SF-FDTD scheme.

In the present work, we propose a novel method for predicting bending parameters from theoretical data using deep learning tools. This study continues our previous work on the estimation of characteristic parameters of holographic volume gratings based on machine learning, as presented in  [21], where we first demonstrated the feasibility of using deep learning models to estimate grating parameters. Here, we extend that foundation by investigating optimal architectures for predicting bending parameters and comparing the performance of a Feedforward Neural Network (FNN) and a hybrid Convolutional Neural Network (CNN)–FNN model. Slight improvements in accuracy and robustness were observed in the hybrid model, even when a smaller number of weights was applied.

FNNs and CNNs have played pivotal roles in advancing machine learning capabilities, offering robust solutions to vast real-world problems [22,23,24,25]. FNNs, also known as dense networks, consist of layers where each neuron is connected to every neuron in the subsequent layer, with intertwined nonlinear activation functions. This architecture enables FNNs to learn intricate relationships within data, making them suitable for various tasks such as classification, regression, and pattern recognition [26,27,28]. On the other hand, CNNs leverage convolutional layers to extract spatial features from input data. Although CNNs are not as commonly used for regression as for classification due to their design being intrinsically well suited for pattern extraction and hierarchical feature learning, they have proven effective for regression tasks like image super-resolution  [29,30], depth estimation [31], Image Reconstruction [32], and others.

Regarding the structure of this paper, in the first section, we will set the theoretical basis for developing our RCWA-Shooting method approach and the theoretic basis for the implemented deep learning solutions. The results of our RCWA-SM study on the standard RCWA and FDTD models will be presented. Then, the application and results of the deep learning solutions will be discussed. Finally, and serving as the conclusion, a discussion on the validity and suitability of both approaches will be conducted.

2. Theory

In this section, we present the theoretical models and methodologies applied in this work. These include the development of a framework that incorporates bending into the RCWA for accurately modelling nonplanar grating structures. Furthermore, we outline the theoretical basis for the deep learning approach employed to predict bending behavior based on the RCWA-SM model.

2.1. Rigorous Coupled Wave Analysis—Shooting Method Approach

The exposure of a material to an interference pattern results in a modulation of the medium’s dielectric permittivity that can be expressed as

$$\epsilon (z,x)=\sum _h{\epsilon }_h\left(z\right)exp\left[jhK\left(xsin\phi +zcos\phi \right)\right],$$

(1)

where K is the grating vector, and $\phi$ is the grating slant angle.

We can parametrize a bending in this dielectric modulation pattern as a characteristic polynomial $W\left(z\right)$ as follows:

$$W\left(z\right)=a_0+a_1\left({\displaystyle \frac{z}{d}}\right)+a_2{\left({\displaystyle \frac{z}{d}}\right)}^2+a_3{\left({\displaystyle \frac{z}{d}}\right)}^3,$$

(2)

where d is the optical depth.

Introducing it as an expansion factor in the dielectric permittivity expression embedded in ${\epsilon }_h\left(z\right)$, we have the following:

$${\epsilon }_h\left(z\right)={\epsilon }_h·exp\left[jhW\left(z\right)\right],$$

(3)

and then, the more specific formula for the permittivity can be written as

$$\epsilon (z,x)=\sum _{h=-1}^1{\epsilon }_h·exp\left\{jh\left[W\left(z\right)+K\left(xsin\phi +zcos\phi \right)\right]\right\}.$$

(4)

Figure 1 shows the model of the thick hologram used for the analysis where the bending in the interference fringes can be identified. The terms $U_x^{\left(i\right)}$ and $S_y^{\left(i\right)}$ are related to the amplitudes of the electromagnetic field, which is detailed next.

Following a similar formalism to [33], we can describe the electromagnetic interaction as follows:

$${\displaystyle \frac{d}{dz}}\left[\begin{array}{c}⋮\\ S_y^{\left(i\right)}\\ ⋮\\ U_x^{\left(i\right)}\\ ⋮\end{array}\right]=\left[\begin{array}{cc}j{\mathbf{k}}_z& k_0\mathbb{I}\\ {\mathbf{k}}_x^2/k_0-k_0\mathit{ϵ}& j{\mathbf{k}}_z\end{array}\right]\left[\begin{array}{c}⋮\\ S_y^{\left(i\right)}\\ ⋮\\ U_x^{\left(i\right)}\\ ⋮\end{array}\right]$$

(5)

where $S_y^{\left(i\right)}$ and $U_x^{\left(i\right)}$ represent the amplitudes of the ith space harmonic components of the electric and magnetic fields, respectively, $\mathbb{I}$ is the identity matrix, ${\mathbf{k}}_z={\left[k_{g,z}^{\left(i\right)}{\delta }_{i,j}\right]}_{i,j}={\left[k_{g,z}^{\left(i\right)}\right]}_{i,i}$ and ${\mathbf{k}}_x={\left[k_x^{\left(i\right)}\right]}_{i,i}$ are diagonal matrices, and the wave vector is $k_x^{\left(i\right)}=k_{0}nsin\theta -iKsin\phi$; the z-component of the grating wave vector is defined as follows:

$$k_{g,z}^{\left(i\right)}=\left\{\begin{array}{cc}\sqrt{n^2{k}_0^2-k_x^{\left(i\right)2}}\hfill & \mathrm{if}n{k}_0\ge k_x^{\left(i\right)}\hfill \\ j\sqrt{-n^2{k}_0^2+k_x^{\left(i\right)2}}\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

(6)

where $k_0$ is the free-space wave number, n and $\theta$ are the refractive index values of the recording medium and the angle of incidence in the recording medium, respectively, and $\phi$ is the grating slant angle.

$${\mathbf{k}}_z=\left[\begin{array}{ccccc}\ddots & & & & ⋰\\ & k_{g,z}^{(i-1)}& 0& 0& \\ & 0& k_{g,z}^{\left(i\right)}& 0& \\ & 0& 0& k_{g,z}^{(i+1)}& \\ ⋰& & & & \ddots \end{array}\right],{\mathbf{k}}_x=\left[\begin{array}{ccccc}\ddots & & & & ⋰\\ & k_x^{(i-1)}& 0& 0& \\ & 0& k_x^{\left(i\right)}& 0& \\ & 0& 0& k_x^{(i+1)}& \\ ⋰& & & & \ddots \end{array}\right]$$

and $\mathit{ϵ}={\left[{ϵ}_{i-j}\right]}_{i,j}$ is a Toeplitz matrix, wherein we then have the following:

$$\mathit{ϵ}=\left[\begin{array}{ccccc}\ddots & & & & ⋰\\ & {\epsilon }_0& {\epsilon }_{-1}& {\epsilon }_{-2}& \\ & {\epsilon }_1& {\epsilon }_0& {\epsilon }_{-1}& \\ & {\epsilon }_2& {\epsilon }_1& {\epsilon }_0& \\ ⋰& & & & \ddots \end{array}\right]$$

with respect to the boundary conditions and letting $R_i$ be the normalised electric-field amplitude of the ith reflected wave in region I and $T_i$ be the normalised electric-field amplitude of the transmitted wave in region II:

$$\begin{array}{ccc}{\mathbf{S}}_y\left(0\right)& =& {\left[\begin{array}{ccccc}\dots & R_{-1}& 1+R_0& R_{+1}& \dots \end{array}\right]}^T,\\ \\ {\mathbf{S}}_y\left(d\right)& =& {\left[\begin{array}{ccccc}\dots & T_{-1}{e}^{j{k}_{g,z}^{(-1)}d}& T_0{e}^{j{k}_{g,z}^{\left(0\right)}d}& T_{+1}{e}^{j{k}_{g,z}^{(+1)}d}& \dots \end{array}\right]}^T,\\ \\ {\mathbf{U}}_x\left(0\right)& =& {\displaystyle \frac{-j}{k_0}}{\left[\begin{array}{ccccc}\dots & -k_{I,z}^{(-1)}{R}_{-1}& k_{I,z}^{\left(0\right)}(1-R_0)& -k_{I,z}^{(+1)}{R}_{+1}& \dots \end{array}\right]}^T,\\ \\ {\mathbf{U}}_x\left(d\right)& =& {\displaystyle \frac{-j}{k_0}}{\left[\begin{array}{ccccc}\dots & k_{II,z}^{(-1)}{T}_{-1}{e}^{j{k}_{g,z}^{(-1)}d}& k_{II,z}^{\left(0\right)}{T}_0{e}^{j{k}_{g,z}^{\left(0\right)}d}& k_{II,z}^{(+1)}{T}_{+1}{e}^{j{k}_{g,z}^{(+1)}d}& \dots \end{array}\right]}^T.\end{array}$$

(7)

Indeed, with $\left[k_{\mathcal{L},z}\right]$ being the diagonal matrices of the z-components of the wave vectors in the outer media, we can rewrite the last two in terms of the first ones:

$$\left[\begin{array}{c}{\mathbf{U}}_x\left(0\right)={\displaystyle \frac{j}{k_0}}{\mathbf{k}}_{I,z}\left({\mathbf{S}}_y\left(0\right)-2{\left[{\delta }_{i0}\right]}_i\right)\\ \\ {\mathbf{U}}_x\left(d\right)={\displaystyle \frac{-j}{k_0}}{\mathbf{k}}_{II,z}{\mathbf{S}}_y\left(d\right)\end{array}\right].$$

(8)

Similarly, following [34], if we now define the reflectance and transmittance vectors, $R=\left(R_i\right)i$ and $T=\left(T_j\right)j$, we can express the boundary conditions in an alternative manner:

$$\left[\begin{array}{ccccccc}{\mathbf{S}}_y\left(0\right)& =& \mathbf{R}+{\left[{\delta }_{i0}\right]}_i& & {\mathbf{S}}_y\left(d\right)& =& \mathbf{T}{e}^{j{\mathbf{k}}_{z}d}\\ \\ {\mathbf{U}}_x\left(0\right)& =& {\displaystyle \frac{-j}{k_0}}{\mathbf{k}}_{I,z}\left({\left[{\delta }_{i0}\right]}_i-R\right)& & {\mathbf{U}}_x\left(d\right)& =& {\displaystyle \frac{-j}{k_0}}{\mathbf{k}}_{II,z}T{e}^{j{\mathbf{k}}_{z}d}\end{array}\right].$$

(9)

We will modify the problem, transforming it to an Initial Value Problem (IVP), with the initial conditions given by

$${\widehat{S}}_y(z_0;p)={\left({\delta }_{ip}\right)}_i={\left[\begin{array}{cccccc}\dots & \dots & 0& 1& 0& \dots \end{array}\right]}^T$$

(10a)

$${\displaystyle \frac{d}{dz}}{\widehat{S}}_y(z_0;p)=-j{\left(k_{\mathcal{L},z}^{\left(i\right)}{\delta }_{ip}\right)}_i={\left[\begin{array}{cccccc}\dots & \dots & 0& -j{k}_{\mathcal{L},z}^{\left(p\right)}& 0& \dots \end{array}\right]}^T$$

(10b)

which, in turn, define those of ${\displaystyle {\widehat{U}}_x(z_0;p)}$ and its derivative through Equation (8):

$$\begin{array}{c}{\displaystyle {\widehat{U}}_x(z_0;p)=\frac{-j}{k_0}{\left[\begin{array}{cccccc}\dots & \dots & 0& k_{II,z}^{\left(p\right)}& 0& \dots \end{array}\right]}^T}\end{array}.$$

The general solution will then be a linear combination of the particular solutions, with some constants ${\{a_p,b_p\}}_{p\in \left\{-N,\dots ,N\right\}}$ to be determined:

$$S_y\left(z\right)=\sum _{p=-N}^N{a}_p{\widehat{S}}_y(z;p),$$

(11a)

$$U_x\left(z\right)=\sum _{p=-N}^N{b}_p{\widehat{U}}_x(z;p).$$

(11b)

By construction, Equation (11) at $z=z_0$ becomes

$$\left[\begin{array}{c}{\displaystyle S_y^{\left(i\right)}\left(z_0\right)=\sum _{p=-N}^N{a}_p{\widehat{S}}_y^{\left(i\right)}(z_0;p)=\sum _{p=-N}^N{a}_p{\delta }_{ip}=a_i,}\\ {\displaystyle U_x^{\left(i\right)}\left(z_0\right)=\sum _{p=-N}^N{b}_p{\widehat{U}}_x^{\left(i\right)}(z_0;p)=-j\sum _{p=-N}^N{b}_p\frac{k_{II,z}^{\left(i\right)}}{k_0}{\delta }_{ip}=-j\frac{k_{II,z}^{\left(i\right)}}{k_0}{b}_i,}\end{array}\right].$$

(12)

This helps us to rewrite the boundary conditions by considering $z_0=d$.

$$\begin{array}{cc}\hfill \blacksquare z=0& \blacksquare z=d\hfill \end{array}$$

(13a)

$$\begin{array}{cc}\hfill {\delta }_{i0}+R_i=\sum _p{a}_p{\widehat{S}}_y^{\left(i\right)}(0;p)& T_i=a_i{e}^{-j{k}_{g,z}^{\left(i\right)}d}\hfill \end{array}$$

(13b)

$$\begin{array}{cc}\hfill k_{I,z}^{\left(i\right)}\left({\delta }_{i0}-R_i\right)=j{k}_0\sum _p{b}_p{\widehat{U}}_x^{\left(i\right)}(0;p)& T_i=b_i{e}^{-j{k}_{g,z}^{\left(i\right)}d}\hfill \end{array}$$

(13c)

In the previous relationships, there is a total of $4N$ equations for the four sets of N unknowns: $a_p$, $b_p$, $R_i$, and $T_j$.

Proceeding as we did for the previous formalism, we only need to choose the boundary condition at $z=0$ with $a_i=b_i$:

$$2k_{I,z}^{\left(i\right)}{\delta }_{i0}=\sum _p\left[k_{I,z}^{\left(i\right)}{\widehat{S}}_y^{\left(i\right)}(0;p)+j{k}_0{\widehat{U}}_x^{\left(i\right)}(0;p)\right]a_p.$$

(14)

This is done to construct a system of equations for the constants $a_p$ in the form ${\sum }_p{M}_{ip}{a}_p=n_i$, where

$$\left.\begin{array}{ccc}{M}_{ip}\hfill & =\hfill & k_{I,z}^{\left(i\right)}{\widehat{S}}_y^{\left(i\right)}(0;p)+j{k}_0{\widehat{U}}_x^{\left(i\right)}(0;p)\hfill \\ \\ n_i\hfill & =\hfill & 2k_{I,z}^{\left(i\right)}{\delta }_{i0}\hfill \end{array}\right\}a=M^{-1}n$$

(15)

are known.

In this way, we can compute the field value at the first surface using Equation (15), thus obtaining the reflection and transmission coefficients.

Figure 2 summarises the sequence of the implementation of the RCWA-SM method detailed in this section. The flowchart shows how an initial boundary value problem (BVP) is transformed into an initial value problem (IVP) and finally to a linear system of equations to solve the reflection and transmission coefficients.

2.2. FNNs and CNNs Theory Foundations

Deep learning (DL), a subset of machine learning encompassing all neural network architectures, has proven to be a revolutionary field with applications that impact a wide range of scientific and technological areas  [18,35,36]. These algorithms mimic the structure and function of the human brain and are used to recognize patterns for supervised, semi-supervised, or unsupervised purposes. The most significant distinction of DL algorithms compared to traditional machine learning algorithms is that DL algorithms autonomously learn hierarchical representations of data, eliminating the need for hand-crafted feature engineering to transform the data into a more suitable form.

All artificial neural networks are composed of interconnected nodes, or neurons, organized into layers. FNNs are the simplest form of neural networks, consisting of an input layer, an output layer, and one or more hidden layers in between, as illustrated in Figure 3a. Each hidden layer is essentially a matrix and an activation function that introduces nonlinearity to the system, which is a fundamental trait of neural networks. The components of these matrices are called weights, and they are updated iteratively to minimize the loss function between the model and the target outputs, typically using a gradient descent-based optimization algorithm  [18,19,35,36].

On the other hand, a one-dimensional CNN hidden layer consists of several filters or kernels that slide along the input or preceding layers sequence, as shown in Figure 3b, with its subsequent activation function [35,37]. The input data are typically a one-dimensional array featuring some spatial auto-correlation, such as a time series, a sequence of words in natural language processing or, in our case, a diffraction pattern. This allows kernels to extract the underlying patterns derived from that spatial auto-correlation of the input. In this case, the weights consist of the components of the convolutional kernels, which are updated in a completely analogous fashion to the FNNs mentioned above.

FNNs and CNNs differ fundamentally in how they handle data. FNNs treat individual data points as separate, fixed units, while CNNs are designed to capture local correlations, spatial hierarchies, and patterns within data through their sliding convolutional kernels. This makes CNNs ideal for image and signal processing, in which the spatial and temporal dependencies are present and critical. By optimising their kernel weights through backpropagation and gradient descent, CNNs can improve their ability to generalise spatial or temporal auto-correlated data predictions.

Given our objective of predicting the bending polynomial coefficients and other fringe parameters from spatially auto-correlated experimental diffraction efficiency patterns, it is essential to consider the inclusion of CNN layers in our deep learning solution. This study thoroughly evaluates their suitability, as discussed in the present work.

3. Results

This section presents the results of our study on modelling and predicting bending phenomena in holographic volume gratings. We begin by evaluating RCWA-SM for modelling diffraction efficiency in gratings with varying degrees of bending, comparing with Kubota’s method and Split-Field Finite-Difference Time-Domain (SF-FDTD) simulations. Following this, we explore deep learning approaches to estimate bending parameters. Specifically, we compare two architectures: an FNN and a CNN–FNN hybrid designed to predict bending polynomial coefficients, optical depth, and index modulation based on diffraction patterns.

3.1. Suitability of the RCWA-Shooting Method

This section covers the results of different gratings with different bending degrees. The bending scheme considered is the one defined by Equation (2). The term $a_0$ is not taken into account, since it only implies a displacement in the fringes of the grating, not influencing the diffraction efficiency of the gratings. The results here compare the diffraction efficiency obtained through Kubota’s formalism, RCWA-SM (introduced in Section 2.1), and the results from SF-FDTD method. The diffraction grating simulated is a transmission grating $\varphi ={90}^{\circ }$ with ${\theta }_{\mathrm{ref}}={\theta }_{\mathrm{obj}}=13.{65}^{\circ }$. These angles produce a grating with a spatial frequency of $\mathsf{\Delta }=0.823$ $\mathsf{\mu }$m. The average refractive index was $n_0$ = 1.63, and the modulation index was $n_1$ = 0.025. The thickness of the grating wa $d=12$ $\mathsf{\mu }$m. This setup defines a volume grating where only one diffracted order will be relevant in efficiency since the parameter $Q\gg 1$. Regarding FDTD simulations, the grid density was ${\lambda }_0/40$, and the illumination wavelength wa ${\lambda }_0=633$ nm. The time resolution was varied as a function of the incidence angle ${\theta }_0$, ensuring in all cases the Courant–Friederich–Lewy condition defined in Equation (14) in [38], and the number of time steps in the simulation was chosen as 5000, ensuring the steady state of the electromagnetic waves inside the diffractive element. Due to the spatial symmetry, a 2D simulation with 1238 × 52 grid cells was considered.

Firstly, simulations for a specific value of the $a_1,a_2$, and $a_3$ were considered. More specifically, Figure 4 shows the results for the grating matching the cases summarised in

Table 1. Parameter values considered for $a_i$ in Figure 4.

Case	${\mathit{a}}_1$	${\mathit{a}}_2$	${\mathit{a}}_3$
0	0	0	0
1	0.5	0	0
2	0	0.5	0
3	0	0	0.5
4	0.1	0.1	0.1
5	0.1	0.2	0.2
6	0.1	0.3	0.3

. For Figure 4a, the different methods show accurate results when no bending was present in the diffraction grating. The results in Figure 4b–d show how the bending induced a displacement on the curves for the SM-RCWA and Kubota compared to the FDTD. This displacement when bending is present is mainly because both RCWA-SM and Kubota consider a constant grating vector modulus along z and introduce the bending as a phase term. In contrast, the FDTD computes the diffraction efficiency directly using the interaction of light inside the periodic structure with the bending as they travel along the z direction. This approach is more accurate, and it translated in the results as a slight displacement of the diffraction efficiency curves compared to the RCWA-SM and Kubota. For the cases where $a_2$ and $a_3$ were the only parameters inducing the bending, SM-RCWA and Kubota’s implementations maintained the deviation to some degree (compared to the FDTD). As the modulus of the bending roe ($\left|\mathbf{a}\right|$) and the amplitude of $a_2$ and $a_3$ were prominent, the SM-RCWA and Kubota significantly reduced its performance compared to the FDTD. The influence of the bending can be summarised as an angular displacement of the diffraction curve; the angular diffraction response was not symmetric around Bragg’s angle, and a rise in the curves induced a disappearance of the null diffraction efficiency values.

For analysing more precisely the performance of the different methods, Figure 5 shows the error between the SM-RCWA and Kubota’s against the FDTD for the various scenarios here considered (summarised in Table 1). From Figure 5a, it can be seen that the SM-RCWA showed quite good results under the ones obtained by the traditional Kubota implementation.

Figure 6 shows the results for a more profound analysis when only varying $a_1$ from −0.5 up to 0.5. Figure 6a represents the variation in the amplitude of the maximum of the diffraction efficiency. Figure 6b shows the angle where the maximum value of the diffraction efficiency was located. Figure 6c,d represent the amplitude of the first lobes around the Bragg’s angle. Following the same procedure, Figure 6d,f show the position of these lateral lobes as a function of the $a_1$ parameter. The bottom graphs in Figure 6g,h represent the correlation coefficient of the diffraction efficiency compared to the response of the grating without bending and the deviation of the Bragg’s angle compared to the ideal case (no bending), respectively. The analysis of Figure 6 implies that the amplitudes analysed remained constant for different values of $a_1$.

Figure 7 and Figure 8 show the same analysis for $a_2$ and $a_3$, respectively. The linear dependence of the angular variations for $a_2$ seems to be not as linear as the case shown in Figure 6. The variations in amplitude also became more significant. The performance of the SM-RCWA and Kubota seems to be consistent with the previous results.

Figure 8 shows the more extreme results due to the variation in $a_3$. Regarding amplitude, it can be seen that as $a_3$ rises, the diffraction efficiency significantly reduces. Figure 8g shows that the diffraction efficiency curve is progressively less correlated to the ideal response dealing with the disappearance of the lateral lobes of the curve, as it can be seen in the abrupt discontinuities on the tracking of the lateral lobes in Figure 8c–f. These discontinuities in the curves are produced by the severe degradation of the diffraction efficiency curve that induces the vanishing of lateral lobes.

3.2. Deep Learning-Based Bending Estimator

In this section, we present and show the results of the two deep learning models developed: one strictly based on FNNs architecture and the other one being a hybrid of CNNs and FNNs for the estimation of bending parameters as three-order polynomial coefficients alongside optical depth and index modulation. The models were trained using simulated data generated across a range of standardised index modulations, optical depths, and the three coefficients ($a_1$, $a_2$, $a_3$) characterising the third-degree bending polynomial ($a_0$ leaves the grating’s diffraction efficiency pattern invariant). To ensure a fair comparison between the two models and to investigate the impact of CNN layers on the training process, both architectures were constructed to possess a comparable number of weights.

A five feedforward layer model, with interleaved Rectified Linear Unit (ReLU) activation functions and 20% dropouts (to prevent overfitting), was employed as our first simpler FNN model, with a total 1,459,120 weights, as the code image shows in Appendix A.1. For our hybrid model instead, we used 1- to 64- and 64- to 128-channel seven-sized kernel CNN layers with interleaved ReLU activation functions [39] and a maximum pooling of two kernels before three more straightforward feedforward layers interleaved again with ReLU activation functions and 20% dropouts, with a total number of weights of 1,247,808, as shown in Appendix A.2.

We opted for an L1 mean absolute error loss function and a batch size of 100. For the optimisation algorithm, we chose Root Mean Square Propagation (RMSProp) after we tried various optimisation algorithms and observed that it yielded faster convergence and improved stability compared to others such as Adam, Adagrad, or Stochastic Gradient Descent (SGD). RMSProp is designed to address some of the limitations of traditional gradient descent methods, such as the sensitivity to learning rates and the scaling of gradients, by adjusting the learning rates for each parameter in the neural network based on the magnitudes of their gradients. It also scales down the gradients of parameters with more significant updates and those with more minor updates. This helps prevent gradients from exploding or vanishing [40,41]. A learning rate of $1\times {10}^{-5}$ with a momentum of 0.8 was chosen for the training.

FNN Only Versus Hybrid Model Comparison

In this section, the FNNs architecture and the hybrid of CNNs and FNNs are compared in terms of their accuracy, robustness, and generalisation. A dataset of 18, 000 (17,000 for training and 1000 for testing) Shooting method solved RCWA bending simulations of diffraction patterns with $a_1$, $a_2$, and $a_3$ bending parameters ranging from −0.5 to 0.5, optical depths of 10–20 $\mathsf{\mu }$m, and index modulations of 0.01–0.03. The reason behind keeping an amount of 17,000 of data for the training and 1000 of data for the test is that while the test set mades up only 5.6% of the total dataset, it was indeed large enough to validate the model itself by allowing us to employ a larger amount of data for the training. We present both the bending-only (BO) estimation and bending and n1, d estimation versions of the models. The recording theta and the object theta were symmetrically set to 15 and −15 degrees, respectively, and the recording lambda was set to 633 nm. The simulations were conducted for an air medium air, with a medium average refraction index of 1.63. Both trainings were conducted with 8000 epochs. This choice was empirically founded, as, after approximately this number of epochs, the training and validation losses stabilised, offering no significant improvement with further epochs.

Regarding the simulation time and computational performance, both the RCWA-SM simulations and the deep learning models were developed fully in Python and executed in a computer with the following main specs: Intel^® Core^™ i7 13700, 16 GB RAM DDR4-3200 MHz, and 512 GB PCIe NVMe. A total of 8000 epoch trainings took approximately 1 h and 32 min for the FNN model and 2 h and 21 min for the hybrid one—due to its convolutional layers being more computationally costly.

In

Table 2. Mean absolute error comparison between hybrid and FNN models (standardised values).

Parameter	Error FNN	Error FNN (BO)	Error Hybrid	Error Hybrid (BO)
$a_1$	0.044554	0.043941	0.041073	0.040534
$a_2$	0.048321	0.048787	0.045769	0.045369
$a_3$	0.015857	0.017267	0.020191	0.018063
d	0.007719	-	0.008780	-
$n_1$	0.007626	-	0.010430	-
total bending	0.036244	0.036665	0.035678	0.034655
total	0.024815	0.036665	0.025249	0.034655

are shown the results for the Mean Absolute Errors for both models and for $a_1$, $a_2$, and $a_3$ for the optical depth and index modulation. While both methods achieved good results for the prediction of the aforementioned parameters, slight differences in the overall error were observed between the proposed methods, although the hybrid model achieved slightly better accuracy for the overall bending coefficient parameters, while the FNN achieved better optical depth and index modulation prediction performance. A further standard deviation analysis of the difference between the model and target values was made, shown in

Table 3. Standard deviation comparison between hybrid and FNN models.

Parameter	sd FNN	sd FNN (BO)	sd Hybrid	sd Hybrid (BO)
$a_1$	0.087227	0.088635	0.075095	0.073375
$a_2$	0.087632	0.090564	0.076039	0.075379
$a_3$	0.021851	0.023284	0.028523	0.026634
d	0.010985	-	0.011818	-
$n_1$	0.011211	-	0.013305	-

, showcasing that the hybrid model featured lower standard deviations for the prediction of the $a_1$ and $a_2$ coefficients, suggesting a more robust behaviour with less predictions falling outside a normal range. An outlier count analysis (

Table 4. Outliers comparison between hybrid and FNN models. The threshold for the outlier criteria was set to ±0.3 for $a_1$ and $a_2$, ±0.08 for $a_3$, and ±0.05 for d and $n_1$ standarised values.

Parameter	Outliers FNN	Outliers FNN (BO)	Outliers Hybrid	Outliers Hybrid (BO)
$a_1$	22	25	16	15
$a_2$	22	24	15	15
$a_3$	9	9	12	13
d	3	-	1	-
$n_1$	4	-	2	-

) was also conducted, which confirmed our previous observations.

When we compare the bending phase polynomials that the FNN (BO) model predicted and that the hybrid (BO) model predicted with the base true data in the test subset, an excellent concordance can be observed, as shown in Figure 9a,b. When we did the proper RCWA-SM reconstruction of the diffraction patterns out of the aforementioned predicted and base bending parameters, we observe that the results are consistently closely to being identical, as shown in Figure 9c,d, proving the excellent efficiency of our bending prediction despite the non-negligible (0.04) bending coefficient losses. It can be noted that, unlike variations in optical depth and index modulation, variations of this magnitude (0.04) in bending did not significantly impact the diffraction patterns. This points out that the limit of accuracy in our prediction models are primarily constrained by such variations not translating into noticeable changes in optical properties.

4. Conclusions

This work presents the results of a novel RCWA approach for transmission holographic gratings with bending. These results were compared with those obtained from a set of rigorous SF-FDTD and Kubota’s simulations, providing a benchmark for the performance of our approach. The analysis covers a range of different bendings and a parametric analysis of the various parameters that define the bending. The results show that our method is consistent and quite accurate for low and moderate amplitudes of the bending of the fringes ($\sum |a_i|<0.5$). However, as the amplitude of the bending increased, our approach began to deviate from the numerical results, leading to differences in the angles at which the maximum diffraction efficiency could be found, as well as in the amplitudes of this parameter.

A deep learning optimum model for the estimation of bending polynomial coefficients was developed, as we compared the performance of two different architecture-based models. The results indicate that while both models achieved good predictive performance, the hybrid model slightly outperformed the FNN-only model in predicting bending coefficients, featuring lower standard deviations and fewer outliers, while the former model excelled in predicting optical depth and index modulation. Despite small bending coefficient losses (0.04), the predicted and true bending phase polynomials revealed high concordance, with the RCWA-SM reconstructions of diffraction patterns being nearly identical. This indicates that the limits of the accuracy of the estimators are due to the insensitivity of the optical system to small variations in bending.

In future works, the authors plan to extend the current formalism for considering different type of gratings, e.g., slanted and reflection gratings, and to compare the performance of the scheme with experimental results.