<figure><img src="image_1.png" /></figure>

HAL
open science

# Analysis of Faceted Gratings using C-method and Polynomial Expansion

Gérard Granet, Kofi Edee

► To cite this version:

Gérard Granet, Kofi Edee. Analysis of Faceted Gratings using C-method and Polynomial Expansion. 2023. hal-04362636

HAL Id: hal-04362636

https://hal.science/hal-04362636v1

Preprint submitted on 23 Dec 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Article

# Analysis of Faceted Gratings using C-method and Polynomial Expansion

Gérard Granet <fnref n="1" />,†,‡ ID, Kofi Edee <fnref n="2" />,‡ ID

**Abstract:** The coordinate-transformation-based differential method developed by Chandezon et al is recognized as one of the simplest and most versatile approaches for modeling surface-relief gratings. In this study, we present a novel numerical solution using Legendre polynomial expansion, enabling to deal efficiently with faceted gratings. Additionally, we introduce an oblique coordinate transformation to analyze overhanging faceted gratings. Notably, the C-method with polynomial expansion (CPE) demonstrates a dramatic improvement in convergence speed compared to the Fourier Modal Method (FMM).

**Keywords:** Diffraction gratings; Computational electromagnetic methods; Spectral method; Curvilinear coordinates; Method of moments; Legendre polynomials; Oblique coordinates.

## 1. introduction

Faceted gratings are employed in various applications involving diffractive optical elements. These applications span security features [1], optical communication [2] [3], and laser systems. Currently, they are frequently utilized in augmented and mixed reality applications. In the realm of head-mounted displays [4], slanted gratings are commonly used to couple light into optical light guides due to their high efficiency in a specific diffraction order, enabling compact designs. Regardless of the application, the grating's role is to manipulate light in terms of direction, intensity, or polarization. While the choice of a faceted grating may be influenced by superior performance compared to other solutions, ease and cost of manufacturing are also crucial criteria, particularly for general public applications. The desired optical response of a faceted grating is achieved by manipulating its period, which is relatively straightforward, and, more importantly, by configuring the arrangement of its facets. In simpler terms, an optimization process is required. While there are various optimization algorithms available, they all entail solving the same problem multiple times with different adjustable parameters. Since the diffraction effect of the grating can only be simulated accurately through rigorous methods, the overall process is time-consuming. In recent years, the advancement of computers has facilitated the application of the finite element method (FEM) to solve Maxwell's differential equations. This method is particularly well-suited for addressing periodic structures or models with arbitrary shapes. Among the existing grating methods, the curvilinear coordinate method [5], [6] (the C method) and the coupled-wave method [7], [8] (also known as the Fourier modal method or FMM) are considered the simplest. These methods share similarities in their respective domains of application, formulating the diffraction problem as a spectral one solved in discrete Fourier space. Consequently, their strengths and weaknesses are mainly related to the use of a Fourier basis. For instance, the C method converges slowly for profiles with sharp edges, while FMM is not well-suited for metal-dielectric structures

**Citation:** Lastname, F.; Lastname, F.; Lastname, F. Title. *Sensors* **2023**, *1*, 0.
https://doi.org/

Received:

Revised:

Accepted:

Published:

**Copyright:** © 2023 by the authors.
Submitted to *Sensors* for possible open
access publication under the terms and
conditions of the Creative Commons
Attribution (CC BY) license (https://
creativecommons.org/licenses/by/
4.0/).

<footnote n="1">¹ Affiliation 1; gerard.granet@uca.fr</footnote>

<footnote n="2">² Affiliation 2; kofi.edee@uca.fr</footnote>

<footnote>* Correspondence: gerard.granet@uca.fr;</footnote>

<footnote>† Current address: Université Clermont Auvergne, Clermont Auvergne INP, CNRS, Institut Pascal, F-63000 Clermont-Ferrand, France;</footnote>

<footnote>‡ These authors contributed equally to this work.</footnote>
in the optical region. Fortunately, improvements achieved for one method can generally be applied to the other. Factorization rules initially discovered [9],[10] and later rigorously established for lamellar gratings in the case of the FMM [11] have been successfully applied to gratings with sharp edges [12],[13]. In the case of smooth surface-relief gratings, the C method is considered more suitable than FMM, especially for avoiding the slicing of sidewalls when they are not aligned with a coordinate axis. However, this advantage becomes questionable for faceted gratings because the derivative of the profile function, which is part of the operator derived from Maxwell's equations, exhibits jumps. These jumps slow down the convergence of the method when relying on Fourier expansion. Returning to FMM and lamellar gratings, which also introduce jumps in the coefficients of some field components, solutions of Maxwell's equations, it was shown that polynomial expansion allows the enforcement of boundary conditions at jump locations, resulting in exponential convergence of the spectral operator [14],[15],[16],[17],[18]. In the initial stages of promoting this approach, the required boundary conditions were explicitly incorporated into the operator through additional equations. While this formulation proved to be highly efficient, its implementation was overly complex, particularly when dealing with crossed binary gratings. Consequently, an alternative approach utilizing polynomial expansion was developed. This approach involves expanding only the continuous field components onto a basis whose elements are themselves continuous at the location of the jump, as outlined in [19]. The latter formulation also demonstrated high efficiency but, in contrast to the former, resulted in a much simpler implementation.

The goal of this paper is to establish a parallel scheme for analyzing faceted gratings, including the case of overhanging gratings within the framework of the C method. Section 2 provides a concise review of the essential aspects of the C method. In Section 3, a matrix eigenvalue equation independent of polarization is derived for oblique coordinates. Section 4 introduces, for the first time, the solution of the spectral problem using Legendre polynomials expansion. The final section is dedicated to evaluating the method through numerical results and comparing them with existing literature.

## 2. Statement of the problem

### 2.1. *Geometry of the profile*

We consider an electromagnetic problem involving two homogeneous, non-magnetic media separated by a cylindrical periodic surface with a period of $d$, which remains invariant along the $y$ axis. This surface is illuminated from either above or below by a linearly polarized monochromatic plane wave with a vacuum wavelength $\lambda$, angular frequency $\omega$, and vacuum wavenumber $k$. The wave vector is inclined at an angle $\theta$ with respect to the $Oz$ axis. The time dependence is expressed by the factor $\exp(i\omega t)$. The unknowns in this context are the complex amplitudes of the reflected and transmitted plane waves. This paper confines its focus to profiles composed of juxtaposed line segments. Figure (1) illustrates some shapes that can be achieved using this approach.
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Illustration of faceted gratings. (a): blazed grating, (b) trapezoidal grating, (c) overhanging grating, (d) Omega-like grating.</figcaption></figure>

Triangular gratings and trapezoidal gratings are evidently part of this family. The most comprehensive method to depict a profile $p(x, z)$ is to parameterize it using two functions, $x = f(x^1)$ and $z = h(x^1)$, where $x^1$ can be chosen, for instance, along the arclength of the profile. This type of description becomes essential for overhanging gratings, such as the one illustrated in case (c) in the Figure (1). The tangent at the profile is expressed as:

$$ \dot{p} = \frac{d_1 h}{d_1 f} \qquad (1) $$

Here, $d_1$ denotes $d/dx^1$. In the case of a faceted grating, one period of the profile is divided into as many sub-domains as there are facets:

$$ \Omega = \bigcup_{q=1}^{Q} \Omega_{q} \qquad (2) $$

For faceted gratings, each facet, bounded by two points with Cartesian coordinates $(x_q, z_q)$ and $(x_{q-1}, z_{q-1})$, is parameterized by a two-parameter functions $x = f_q(x^1)$ and $z = \dot{p}_q f_q(x^1)$, having a local slope $\dot{p}_q$:

$$ \dot{p}_q = \frac{d_1 h_q}{d_1 f_q} = \frac{z_q - z_{q-1}}{x_q - x_{q-1}} \qquad (3) $$

However, for the derivation of a general formulation, we assume that the grating profile is parameterized by two independent parametric functions $f$ and $h$.

## 2.2. Overview of C-method

The solution to the electromagnetic problem requires enforcing boundary conditions at the profile surface. The key idea of C-method is to make this step easy by mapping the profile of the grating to a plane using a new coordinate system $(x^1, x^2, x^3)$ deduced from the Cartesian coordinate system $(x, y, z)$. Maxwell's equations are then expressed under the covariant form, leading to a first-order differential system of equations applicable to both polarizations. The unknowns of these equations are independent of $x^3$, allowing the search for solutions in the form:

$$ \psi(x^1, x^3) = \exp(-ik\gamma x^3)\phi(x^1). \qquad (4) $$

Chandezon called a coordinate system which enables such a kind of solution a "translation coordinate system". For a cylindrical surface invariant along the $Oy$ axis, Plumey et al ([20])
demonstrated that the most general transformation leading to a translation coordinate system takes the form:

$$ 
\begin{aligned}
x &= f(x^1) + c_x x^3 \\
y &= x^2 \\
z &= h(x^1) + c_y x^3
\end{aligned} 
\qquad (5) $$

By choosing $c_x = \cos(\varphi)$ and $c_y = \sin(\varphi)$, the Ox$^3$-axis forms an angle $\varphi$ with respect to the Ox-axis. This results in an inlined (oblique) translation coordinate system. Such a coordinate system becomes essential for solving diffraction problems by profiles similar to profile c) in Figure(1).

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Illustration of an inclined translation coordinate system. Two coordinate surfaces $x^3 = cte$ are translated from each other parallel to $x^3$-axis</figcaption></figure>

A curvilinear coordinate system is characterized by a metric tensor $[g_{ij}]$ where covariant components are obtained from the Jacobian matrix. The $(m, n)$th entry of the Jacobian is denoted as $J_{mn} = \partial_n(x'_m)$, where $x'_m$ represents the coordinates $(x, y, z)$ and $\partial_n$ denotes the partial derivative with respect to $x^n$:

$$ [g_{ij}] = J^t J, \qquad (6) $$

where the exponent $t$ indicates the transposition. For the coordinate system given by Equation (5), the matrix composed of the covariant components of the metric tensor is:

$$ [g_{ij}] = \begin{bmatrix} (\dot{f}^2 + \dot{h}^2) & 0 & \cos(\varphi)\dot{f} + \sin(\varphi)\dot{h} \\ 0 & 1 & 0 \\ \cos(\varphi)\dot{f} + \sin(\varphi)\dot{h} & 0 & 1 \end{bmatrix}, \qquad (7) $$

The determinant $g_q$ is expressed as:

$$ g = (\dot{f} \sin(\varphi) - \dot{h} \cos(\varphi))^2. \qquad (8) $$

Here $\dot{f}$ and $\dot{h}$ denote the derivative with respect to $x^1$ of functions $f$ and $h$ respectively. It has to be emphasized that not all surface relief gratings can be analyzed with the C-method associated with an oblique translation coordinate system. The profile labeled by letter (d) in the figure (1) serves as an example. However, the discussion regarding the existence or non-existence of a translation coordinate system compatible with a profile shape goes beyond the scope of this paper. We are now ready to derive the system of differential equations resulting from the application of C-method

### 3. Spectral formulation of the problem with C-method

#### 3.1. Maxwell's equations under the covariant form

Given a non-orthogonal coordinate system $x^1, x^2, x^3$, it is appropriate to use covariant Maxwell's equations. Using the Einstein convention, in a homogeneous medium with an optical index $v$, free of electrical charges and currents, these equations can be expressed as follows:

$$ 
\begin{aligned}
\bar{\zeta}^{ijk} \partial_j E_k &= -ik \sqrt{g} g^{ij}(Z_0 H_j) \\
\bar{\zeta}^{ijk} \partial_j (Z_0 H_k) &= ikv^2 \sqrt{g} g^{ij} E_j.
\end{aligned} 
\qquad (9) $$
Here the contravariant components $g^{ij}$ of the metric tensor associated with coordinates $x^1, x^2, x^3$ are calculated by inverting the matrix formed by the $g_{ij}$, $g$ being the determinant of the latter. $\partial_i$ represents the partial derivative versus $x^i$ ($\partial/\partial x^i$) and $\xi^{ijk}$ is the Levi-Civita indicator with only non-zero elements being $\xi^{123} = \xi^{231} = \xi^{312} = 1$ and $\xi^{132} = \xi^{213} = \xi^{321} = -1$. Lastly, $Z_0$ is the vacuum wave impedance. For convenience, let $\Lambda^{ij}$ denote the non-null $\sqrt{gg^{ij}}$. They are expressed as follows:

$$ 
\begin{aligned}
\Lambda^{11} &= \frac{1}{\sin(\varphi)\dot{f} - \cos(\varphi)\dot{h}} \\
\Lambda^{13} &= -\frac{\cos(\varphi)\dot{f} + \sin(\varphi)\dot{h}}{\sin(\varphi)\dot{f} - \cos(\varphi)\dot{h}} \\
\Lambda^{31} &= \Lambda^{13} \\
\Lambda^{22} &= (\sin(\varphi)\dot{f} - \cos(\varphi)\dot{h}) \\
\Lambda^{33} &= \frac{\dot{f}^2 + \dot{h}^2}{\sin(\varphi)\dot{f} - \cos(\varphi)\dot{h}}
\end{aligned} 
\quad (10) $$

### 3.2. 2D operator

As the problem is independent of $x^2$ and the medium is homogeneous, the general system represented by Equation (9) separates into two independent subsystems corresponding to TM and TE polarizations. The TE and TM systems involve components $E_2, H_1, H_3$ and components $H_2, E_1, E_3$ respectively. While each component depends on $x^3$ and $x^1$, Equation (4) indicates that only the $x^1$ dependence is an unknown function. Consequently, the partial derivative with respect to $x^3$ is replaced by $-ik\gamma$, and the partial derivative with respect to $x^1$ becomes an ordinary derivative (the symbol $\partial_1$ becomes $d_1$). Introducing $F_2, G_1$ and $G_3$ such that:

$$ \begin{pmatrix} F_2 = \nu E_2 \\ G_1 = Z_0 H_1 \\ G_3 = Z_0 H_3 \end{pmatrix} \text{ for TE polarization } \begin{pmatrix} G_1 = -\nu E_1 \\ G_3 = -\nu E_3 \\ F_2 = Z_0 H_2 \end{pmatrix} \text{ for TM polarization , } \quad (11) $$

we obtain the following equations valid for the two polarizations:

$$ \gamma F_2 = -\nu (\Lambda^{11} G_1 + \Lambda^{13} G_3) \quad (12a) $$

$$ d_1 F_2 = -ik\nu (\Lambda^{31} G_1 + \Lambda^{33} G_3) \quad (12b) $$

$$ \gamma G_1 - \frac{i}{k} d_1 G_3 = -\nu \Lambda^{22} F_2 \quad (12c) $$

## 4. Numerical solution

In essence, the numerical solution of a problem involves the transformation of operators into matrices. During this process, it is crucial to carefully consider the electromagnetic boundary conditions of the field across the discontinuity lines of the tensor $\Lambda^{ij}$. In the problem at hand, the discontinuity comes from the derivative of the profile along $x^1$. From tensor theory, we know that $E_2, E_3, H_2$, and $H_3$ are continuous at nodes $x_q$ whereas $E_1$ and $H_1$ are discontinuous at the same places. We will take that property into account when deriving the matrices associated to the $\Lambda^{ij}$. In this section, we explain the method of solving the eigenvalue equation (12) using the method of moments (MoM). The MoM will allow us to associate matrices with the derivative operator $d_1$ and with the parameters $\Lambda^{ij}$ also considered as operators.
## 4.1. Method of Moments

The MoM is a general method suitable for solving differential equations of the type $\gamma f = \mathcal{L}f$ where $\mathcal{L}$ is a differential operators [21]. Discretizing the equation via the method of MoM yields the matrix equation $\gamma \mathbf{f} = \mathbf{L}\mathbf{f}$ in which the elements of matrices $\mathbf{L}$ and vector $\mathbf{f}$ are such that:

$$L_{ij} = \langle T_i, W_j \rangle, f_j = \langle f, W_j \rangle. \qquad (13)$$

Here $T_i$ and $W_j$ refer to the $i$th testing and $j$th expansion functions respectively. The inner product of any functions $T$ and $W$, $\langle T, W \rangle$, defined on a domain $\Omega$ is defined, for all $x \in \Omega$ as:

$$\langle T, W \rangle = \int_{\Omega} T^{*}(x)W(x)dx, \qquad (14)$$

where the exponent $*$ denotes complex conjugation. From a practical standpoint, the efficiency of numerical computation and memory requirements significantly depend on the choice of the testing basis and the expansion basis, with smaller matrices generally resulting in faster computations and lower memory usage. Of course, from a practical point of view, the smaller the matrix, the faster the numerical computation and the less memory requirement. The efficiency of numerical computation and memory requirements is highly influenced by the choice of the testing basis and the expansion basis. The subsequent sections are dedicated to matrices associated with the derivation operators.

## 4.2. polynomial basis

### 4.2.1. Legendre polynomials

In this section, we solve the eigenvalue problem defined by Equation (12) using the method of moments and Legendre polynomials although more general Gegenbauer polynomials could be used as well. The Legendre polynomials are defined on the unit interval $[-1, 1]$ on which they are orthogonal and complete:

$$\int_{-1}^{1} L_{m}(x)L_{n}(x)dx = \frac{2}{2m+1}\delta_{mn}. \qquad (15)$$

where $L_m$ is the Legendre polynomial of degree $m$ and $\delta_{mn}$ designates the Kronecker symbol. The Legendre polynomials and their derivative satisfy:

$$\int_{-1}^{1} \frac{dL_{n}(x)}{dx} L_{m}(x)dx = \begin{cases} 2 & \text{for } n > 0 \text{ and } m - n \text{ odd} \\ 0 & \text{otherwise} \end{cases}. \qquad (16)$$

The above relation may also be considered as the definition of the entry of the $m$th row and $n$th column of a matrix $\mathbf{d}$

### 4.2.2. Expansion basis and test basis

The initial step in the method of moments involves expanding the unknown functions onto a carefully chosen basis. These basis functions must efficiently consider the continuity and discontinuity properties of the unknowns in the problem. Various approaches can be employed to account for these discontinuities or, more precisely, the physical properties they impose on the unknown functions $F_i$ and $G_j$. One approach is to incorporate them into the basis itself. In our particular problem, the functions $F_2$, $G_1$, and $G_3$ are all pseudo-periodic. Additionally, $F_2$ and $G_3$ exhibit continuity at the nodes $x = x_p$, where $p = 2, \dots, P$. Consequently, when expanding onto a basis, $F_2$ and $G_3$ on one hand and $G_1$ on the other hand should not share the same basis [19]. Our objective is to construct expansion and test bases defined over the entire period of the grating, starting from Legendre polynomials defined on the interval $[-1, 1]$. To achieve this, let's define a set of polynomials $P_n(x)$:

$$P_{n,q}(x) = \begin{cases} L_n(x^1) & \text{if } x \in \Omega_q \\ 0 & \text{if } x \notin \Omega_q \end{cases}, \qquad (17)$$
where $x$ and $x^1$ are linked by the mapping

$$x = 0.5(x_q - x_{q-1})x^1 + 0.5(x_q + x_{q-1}). \qquad (18)$$

Now, let's examine the components $F_2$ and $G_3$. In accordance with Edee and Plumey [19], we introduce a new basis whose elements $\tilde{B}_n$ are derived from $P_{n,q}$. These elements should exhibit continuity at the nodes $x_q$ for $q = 1, \dots, Q-1$ and be linked by the pseudo-periodicity coefficient $\tau$ at nodes $x_0$ and $x_Q$. The pseudo-periodicity coefficient is a phase coefficient, such that the field components satisfy $U(x^1+d) = \tau U(x^1)$, implying that we need to enforce $\tilde{B}_m(x_Q) = \tau \tilde{B}_m(x_0)$. Given that there are $M_q$ basis elements per sub-domain $\Omega_q$, we seek the $\tilde{B}_n$ in the following form: $\tilde{B}_n$ under the form

$$\tilde{B}_m(x) = P_{n,q}(x) + \sum_{q=1}^{Q} a_{M_q,n} P_{M_q,q}, \qquad (19)$$

with

$$m = \begin{cases} n & \text{if } q = 1 \\ n + \sum_{k=1}^{q} M_k & \text{if } 2 \le q < Q-1 \end{cases}. \qquad (20)$$

Hence, there are $Q \times \sum_q M_q$ weighting coefficients $a_{M_q,n}$ which have to be computed by imposing the continuity of $\tilde{B}_m(x)$ at nodes $x_q$, $q = 1, \dots, Q-1$ and its pseudo periodicity at nodes $x_0$ and $x_Q$. In passing, let us denote by $dim$ the total number of basis vectors:

$$dim = \sum_{q=1}^{Q} M_q. \qquad (21)$$

The coefficients $a_{M_q,n}$ are obtained by solving, for every sub-domain $\Omega_q$, $M_k$ linear systems of the form $RA_{M_k,q} = L_{(q)}$ where the matrices $R$ and $L_{(q)}$ are built by using the values of the $P_{n,q}$ at $x = x_{q-1}$ and $x = x_q$, that is $P_{n,q}(x_{q-1}) = (-1)^n$ and $P_{n,q}(x_q) = 1$. $A_{M_q,m}$ is the matrix whose $q$th row of length $Q$ is formed by the $a_{M_q,m}$. Matrix $R$ is as follows:

$$\mathbf{R} = \begin{bmatrix} 1 & -(-1)^{M_1} & 0 & 0 \\ 0 & 1 & -(-1)^{M_2} & 0 \\ 0 & 0 & \cdots & \cdots \\ -\tau(-1)^{M_Q} & 0 & 0 & 1 \end{bmatrix}. \qquad (22)$$

$L_{(q)}$ where the index $q$ refers the $q$th sub-domain is a matrix with $Q$ rows and $M_q$ columns. Its non-zero elements are given below:

$$\begin{aligned} L_{(1),(1,n)} &= -1 \\ L_{(1),(Q,n)} &= \tau(-1)^{n-1} \quad \text{if} \quad n \in [1, M_1] \end{aligned} \qquad (23a)$$

$$\begin{aligned} L_{(q),(q-1,n)} &= (-1)^{n-1} \\ L_{(q),(q,n)} &= 1 \end{aligned} \quad \text{if} \quad 2 \le q \le Q, n \in [1, M_q]. \qquad (23b)$$

The newly defined basis being constructed as linear combination of Legendre polynomials it is convenient to introduce the change of basis matrix $\tilde{K}$ whose columns are the expansion coefficients of the $\tilde{B}_n$ in terms of the $P_{n,q}$ its elements are:

$$\tilde{K}_{mn} = \begin{cases} \delta_{mn} & \text{if } m = M_q \\ a_{M_q,n} & \text{if } m = M_q \end{cases} \qquad (24)$$

Contrary to $F_2$ and $G_3$, $G_1$ is not continuous at nodes $x_q$. That is why we may expand it on a basis $B_m(x)$ defined in manner similar to $\tilde{B}_n$ but with coefficients $a_{M_k,q}$ being equal to zero. Numerical experiment have shown that the latter basis is also convenient to serve
as the basis of the testing functions. As we did for the $\tilde{B}_m$ we define the change of basis matrix $K$ whose elements are:

$$K_{mn} = \begin{cases} \delta_{mn} & \text{if } m = M_q. \\ 0 & \end{cases} \qquad (25)$$

Properties of the new basis

* It has to be emphasized that the test functions $B_n(x)$ are orthogonal to themselves but also to the $\tilde{B}_m(x)$. Indeed we have:

$$\begin{align} \langle B_n, B_m \rangle &= \langle B_n, \tilde{B}_m \rangle = \delta_{mn} \frac{2}{2q+1} \frac{2}{x_q - x_{q-1}} & \text{with} \nonumber \\ M_q &\le q < M_{q+1} \nonumber \\ m &= q + \sum_{k=1}^{q} M_{k-1}, \quad M_0 = 0 \tag{26} \end{align}$$

* The relation (16) between the Legendre polynomials and their derivative may generalized to the new basis. The inner product $\langle B_m(x), d\tilde{B}_n/dx \rangle$ generates the following matrix:

$$\tilde{D} = K^t \mathrm{diag}(d) \tilde{K} \qquad (27)$$

where the exponent $t$ is for transposition and where $\mathrm{diag}(d)$ is a bloc diagonal matrix constituted of the matrices $d$ the dimension of which are $(M_q + 1) \times (M_q + 1)$

### 4.3. Algebraic eigen equation

We are now able to transform into an algebraic eigen system of equations. Since they are continuous $F_2$ and $G_3$ field components are expanded onto $\tilde{B}_m$ basis whereas $G_1$ is expanded onto $B_m$ basis which serves also as test basis. Using a bra-ket notation and Einstein's convention we have:

$$|F_2\rangle = F_2^m |\tilde{B}_m\rangle, \quad |G_3\rangle G_3^m |\tilde{B}_m\rangle, \quad |G_1\rangle = G_1^m B_m\rangle. \qquad (28)$$

The projection of each line of (12) onto the $B_n$ yields:

$$\gamma \mathbf{K}^t \tilde{\mathbf{K}} |F_2\rangle = -\nu (\Lambda^{11} \mathbf{K}^t \tilde{\mathbf{K}} |G_1\rangle + \Lambda^{13} \mathbf{K}^t \tilde{\mathbf{K}} |G_3\rangle), \qquad (29a)$$

$$\tilde{D}|F_2\rangle = -ik\nu (\Lambda^{31} \mathbf{K}^t \tilde{\mathbf{K}} |G_1\rangle + \Lambda^{33} \mathbf{K}^t \tilde{\mathbf{K}} |G_3\rangle), \qquad (29b)$$

$$\gamma \mathbf{K}^t \tilde{\mathbf{K}} |G_1\rangle = -\nu \Lambda^{22} \mathbf{K}^t \tilde{\mathbf{K}} |F_2\rangle + \frac{i}{k} \tilde{D} |G_3\rangle. \qquad (29c)$$

where the $\Lambda^{ij}$ are diagonal matrices whose elements are $\Lambda_{mn}^{ij}$:

$$\Lambda_{mn}^{ij} = \begin{cases} \delta_{mn} \Lambda_1^{ij} & \text{if } 1 \le m \le M_1 \\ \delta_{mn} \Lambda_q^{ij} & \text{if } q + \sum_{k=1}^{q-1} M_k \le m \le q + \sum_{k=1}^{q} M_k. \end{cases} \qquad (30)$$

The fact that the $\Lambda^{ij}$ gives rise to diagonal matrices in polynomial space is an extremely interesting property. It avoids the intricate computation of the convolution matrices of Fourier-based methods which should respect the so-called Fourier factorization rules. The $\mathbf{K}^t \tilde{\mathbf{K}}$ is also a diagonal matrix, hence its product with the $\Lambda^{ij}$ is commutative. Eliminating $|G_3\rangle$ thanks to (29b) yields the following algebraic eigen system:

$$\gamma \begin{bmatrix} |F_2\rangle \\ |G_1\rangle \end{bmatrix} = \mathbf{C} \begin{bmatrix} |F_2\rangle \\ |G_1\rangle \end{bmatrix}, \qquad (31)$$
with

$$ \mathbf{C} = \begin{bmatrix} -\mathbf{\Lambda}^{13} (\mathbf{\Lambda}^{33})^{-1} \mathbf{D} & -\nu (\mathbf{\Lambda}^{33})^{-1} \\ -\nu \mathbf{\Lambda}^{22} + \mathbf{D} (\mathbf{\Lambda}^{33})^{-1} \mathbf{D} & -\mathbf{D} \mathbf{\Lambda}^{31} (\mathbf{\Lambda}^{33})^{-1} \end{bmatrix}. \quad (32) $$

and where the matrix D, which is associated with the derivative operator $(i/k)d_1$ is:

$$ \mathbf{D} = (\mathbf{K}^t \mathbf{K})^{-1} \tilde{\mathbf{D}}. \quad (33) $$

Finally, the field in a translation oblique coordinate system may be expressed as a linear combination of vectors deduced from the eigenvectors of the matrix C. For example, in TE polarization, $E_2(x^1, x^3)$ and $H_1(x^1, x^3)$ write:

$$ \begin{bmatrix} E_2(x^1, x^3) \\ H_1(x^1, x^3) \end{bmatrix} = \sum_{q=1}^{2dim} A_q \exp(-ik\gamma_q x^3) \begin{bmatrix} \frac{1}{\nu} F_{2q}^{m} |\tilde{B}_m\rangle \\ G_{1q}^{m} |B_m\rangle \end{bmatrix}. \quad (34) $$

### 4.4. Application to diffraction gratings

We consider a surface relief grating separating two materials with optical indices $v_1$ and $v_2$. The latter can be either real or complex. The grating is illuminated from above by a unit amplitude plane wave whose wave vector $k_1$ is inclined at $\theta$ to the z axis, see Figure (3).

<figure><img src="image_5.png" /><figcaption>Figure 3. Illustration of the incident and diffracted waves on a surface relief faceted grating. The incident medium has a real optical index.</figcaption></figure>

Solving the matrix eigen equation 31 for each medium gives us two sets of eigenvectors from which we reconstruct easily the vector formed by the components that are continuous at the grating surface. Let us restrict ourselves to TE polarization. We have:

$$ \begin{bmatrix} E_2^{(p)}(x^1, x^3) \\ H_1^{(p)}(x^1, x^3) \end{bmatrix} = \sum_{q=1}^{2dim} A_q^{(p)} \exp(-ik\gamma_q^{(p)} x^3) \begin{bmatrix} \frac{1}{\nu} F_{2q}^{(p)m} |\tilde{B}_m\rangle \\ G_{1q}^{(p)m} |B_m\rangle \end{bmatrix}. \quad (35) $$

At the grating surface, $E_2$ and $H_1$ are continuous, which yields the linear system of equations

$$ \sum_{q=1}^{2dim} A_q^{(1)} \begin{bmatrix} \frac{1}{\nu} |F_{2q}^{(1)m}\rangle \\ |G_{1q}^{(1)m}\rangle \end{bmatrix} = \sum_{q=1}^{2dim} A_q^{(2)} \begin{bmatrix} \frac{1}{\nu} |F_{2q}^{(2)m}\rangle \\ |G_{1q}^{(2)m}\rangle \end{bmatrix}. \quad (36) $$

To go further, we need to sort the eigensolution of C matrix according to whether they correspond to forward waves or backward waves. The forward and the backward, waves correspond to waves that propagate or decay in the positive or the negative direction of
the z-axis respectively. This step is mandatory for enforcing boundary conditions on the grating surface with the S matrix algorithm which results in a matrix equation of the form:

$$ \begin{bmatrix} A_q^{(1)+} \\ A_q^{(2)-} \end{bmatrix} = \mathbf{S} \begin{bmatrix} A_q^{(1)-} \\ A_q^{(2)+} \end{bmatrix}, \qquad (37) $$

with

$$ \mathbf{S} = \begin{bmatrix} \frac{1}{\nu^{(1)}} |F_{2q}^{(1)m}\rangle^+ & -\frac{1}{\nu^{(2)}} |F_{2q}^{(2)m}\rangle^- \\ |G_{1q}^{(1)m}\rangle^+ & -|G_{1q}^{(2)m}\rangle^- \end{bmatrix}^{-1} \begin{bmatrix} -\frac{1}{\nu^{(1)}} |F_{2q}^{(1)m}\rangle^- & \frac{1}{\nu^{(2)}} |F_{2q}^{(2)m}\rangle^+ \\ -|G_{1q}^{(1)m}\rangle^- & |G_{1q}^{(2)m}\rangle^+ \end{bmatrix}. \quad (38) $$

Where the exponents + and - denote forward and backward waves respectively. With a time dependence $\exp(i\omega t)$ the decaying forward, respectively backward waves, are those with a computed eigenvalue with a negative, respectively positive imaginary part. Regarding the computed real eigenvalues, it is numerically observed that they do not depend on the profile function. Providing that the truncation order is large enough, the real eigenvalues obtained in the inclined translation coordinate can be identified with those calculated in the Cartesian coordinate system according to the grating formula. In medium (p), the real diffracted orders are such that:

$$ E_2^{(p)\pm}(x,z) = \sum_m A_m^{(p)\pm} \exp(\mp ik\rho_m^{(p)}z) \exp(-ik\alpha_m x), \quad (39) $$

with $\alpha_m = \nu_1 \sin(\theta) + m \frac{\lambda}{d}$, and $\nu^{(p)} \le \alpha_m \le \nu^{(p)}$ and where $\rho_m^{(p)}$ is defined as follows:

$$ \rho_m^{(p)} = \begin{cases} \sqrt{\nu_p^2 - \alpha_m^2} & \text{if } (\nu_p^2 - \alpha_m^2) \in \mathbb{R}^+ \\ -i\sqrt{\alpha_m^2 - \nu_p^2} & \text{if } (\nu_p^2 - \alpha_m^2) \in \mathbb{R}^- \end{cases}. \quad (40) $$

at a point M, the field can be computed either in terms of variables z and x or in terms of x³ and x¹. Hence, we may write:

$$ E_2^{(p)\pm}(x,z) = E_2^{(p)\pm}(x^1, x^3) \sum_m A_m^{(p)\pm} \exp(-ik\gamma_m^{(p)\pm}x^3)\phi(x^1). \quad (41) $$

then substituting $z(x^3, x^1)$ and $x(x^1, x^3)$ for z and x, in Equation (39), we have:

$$ E_2^{(p)\pm}(x^1, x^3) = \sum_m A_m^{(p)\pm} \exp(-ik(\pm\rho_m^{(p)}\sin\varphi + \alpha_m\cos\varphi))x^3 \exp(-ik(\pm\rho_m^{(p)}h(x^1) + \alpha_m\sin\varphi f(x^1))). \quad (42) $$

The importance of the latter expression is twofold. First, it allows to distinguish among the computed real eigenvalues those that corresponds to forward waves or backward waves. For the diffracted waves such that $-\nu < \alpha_m < \nu$ we have:

$$ \begin{aligned} \gamma_m^{(p)+} &= \rho_m^{(p)} \sin(\varphi) + \alpha_m \cos(\varphi) \\ \gamma_m^{(p)-} &= -\rho_m^{(p)} \sin(\varphi) + \alpha_m \cos(\varphi) \end{aligned} \quad (43) $$

second the dependence on $x^1$ of the diffracted waves, is clearly stated as:

$$ \phi_m^{(p)\pm}(x^1) = \exp(-ik(\pm\rho_m^{(p)}h(x^1) + \alpha_m f(x^1))). \quad (44) $$

In the particular case of faceted gratings, the above expression may be expanded analytically onto the Legendre polynomials and, in a second step, onto the $\tilde{B}_m$ in terms of the Bessel functions of the first kind of order $n + 1/2$. It remains to calculate the other component $H_1$ or necessary to write boundary conditions at $x^3 = 0$. The details of all these calculations
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Overhanging gratings. The end point coordinates of the facets of grating (a) are (0,0), (1.1,1) and 1,1). The end point coordinates of the facets of grating (b) are (0,0), 0.75, 1.25, (1, 1.25),(0.72,0),1,0). All units are arbitrary.</figcaption></figure>

are reported in Appendix A. In conclusion, we can separate the eigenvectors into two sets corresponding to forward waves and backward waves. Each set contains eigenvectors with real eigenvalues and complex eigenvalues. The eigenvectors with real eigenvalues are replaced by their Cartesian plane wave counterpart. The amplitude coefficient of the incident plane wave $A_0^{(p)}$ is set to one and (37) is solved for $A_0^{(1)+}$ and $A_0^{(2)-}$. Since the eigenvectors with real eigenvalues are nothing more than plane waves expressed in the new coordinate system the reflected and transmitted efficiencies in the $n$th order are computed with the same formulas as in the Cartesian coordinate system.

$$ 
\begin{align} 
R_n^{te} &= \frac{|A_n|^2 \rho_n^{(1)}}{\gamma_0^{(1)}} & \text{for reflected TE polarized orders} \nonumber \\ 
T_n^{te} &= \frac{|A_n|^2 \rho_n^{(2)}}{\gamma_0^{(1)}} & \text{for transmitted TE polarized orders} \tag{45} \\ 
R_m^{tm} &= \frac{|A_n|^2 \rho_n^{(1)}}{\gamma_0^{(1)}} & \text{for reflected TM polarized orders} \nonumber \\ 
T_m^{tm} &= \frac{\nu_2^2 |A_n|^2 \rho_n^{(2)}}{\nu_1^2 \gamma_0^{(1)}} & \text{for reflected TM polarized orders} \nonumber 
\end{align} 
$$

# 5. Results

## 5.1. Validation by comparison with published data

In this section, we validate our code by comparing our results with previously published data by Plumey et al [21], which was obtained through another implementation of the C-method and the Fourier Modal Method (FMM). The investigation focuses on two gratings displayed in Fig. 4. Plumey et al’s approach differs from ours in philosophy for the C-method. They exploit the fact that only the grating profile function, not the function itself, appears in the eigenvalue equation. Additionally, their numerical solution is obtained in Fourier space. Plumey et al mention the use of 131 Floquet harmonics for the C-method, while 111 Floquet harmonics and 100 layers were utilized for FMM and the lossy dielectric. With polynomial expansion, we incorporate only 15 basis elements per facet. As expected, excellent agreement is observed between the two variants of the C-method, showing good results for TE polarization and the dielectric grating, but only moderately acceptable results for the lossy dielectric.

## 5.2. Comparison of convergence between FMM and CPE

In this section, we examine a slanted surface relief grating designed for coupling light inside a waveguide, as illustrated by grating (b) in Figure (4). The incident medium is vacuum, and the waveguide is made of glass with a refractive index $v = 1.53$. The
<table><caption>**Table 1.** Comparison between FMM and C-method by Plumey et al [21] on the one hand and our own implementation of C-method in an oblique coordinate system on the other hand. The exponent (a), are for Plumey's results, taken from table 3 in [] The grating is grating (a) in Figure(4). The optical indices are $v^2 = 2.25(D)$, and $v^2 = -21 - i60.4(LD)$. The incident wave parameters are $\theta = 25^\circ$ and $\lambda = 0.7$.</caption><thead><tr><th rowspan="2">orders</th><th colspan="2">C-method (a)</th><th colspan="2">FMM (a)</th><th colspan="2">C-method</th></tr><tr><th>D</th><th>LD</th><th>D</th><th>LD</th><th>D</th><th>LD</th></tr></thead><tbody><tr><td colspan="7"><strong>TE polarization</strong></td></tr><tr><td>R<sub>-2</sub></td><td>0.04506</td><td>0.23791</td><td>0.04502</td><td>0.23602</td><td>0.45071</td><td>0.23869</td></tr><tr><td>R<sub>-1</sub></td><td>0.00316</td><td>0.31502</td><td>0.00315</td><td>0.30576</td><td>0.00316</td><td>0.31410</td></tr><tr><td>R<sub>0</sub></td><td>0.00019</td><td>0.10992</td><td>0.00019</td><td>0.10593</td><td>0.00020</td><td>0.10922</td></tr><tr><td>T<sub>-2</sub></td><td>0.35193</td><td></td><td>0.35197</td><td></td><td>0.35199</td><td></td></tr><tr><td>T<sub>-1</sub></td><td>0.02459</td><td></td><td>0.02459</td><td></td><td>0.02459</td><td></td></tr><tr><td>T<sub>0</sub></td><td>0.5659</td><td></td><td>0.56594</td><td></td><td>0.56591</td><td></td></tr><tr><td>T<sub>1</sub></td><td>0.00913</td><td></td><td>0.00912</td><td></td><td>0.00913</td><td></td></tr><tr><td colspan="7"><strong>TM polarization</strong></td></tr><tr><td>R<sub>-2</sub></td><td>0.03438</td><td>0.45797</td><td>0.03418</td><td>0.42275</td><td>0.03445</td><td>0.45382</td></tr><tr><td>R<sub>-1</sub></td><td>0.00114</td><td>0.14879</td><td>0.00111</td><td>0.14854</td><td>0.00115</td><td>0.15215</td></tr><tr><td>R<sub>0</sub></td><td>0.00004</td><td>0.00841</td><td>0.00004</td><td>0.00322</td><td>0.00005</td><td>0.00616</td></tr><tr><td>T<sub>-2</sub></td><td>0.11085</td><td></td><td>0.11061</td><td></td><td>0.11087</td><td></td></tr><tr><td>T<sub>-1</sub></td><td>0.16566</td><td></td><td>0.16515</td><td></td><td>0.16569</td><td></td></tr><tr><td>T<sub>0</sub></td><td>0.68405</td><td></td><td>0.68501</td><td></td><td>0.68396</td><td></td></tr><tr><td>T<sub>1</sub></td><td>0.00384</td><td></td><td>0.00389</td><td></td><td>0.00386</td><td></td></tr></tbody></table>

coupler's operation requires a single order to be reflected in vacuum and three orders to be transmitted in the guide, with a constraint on the direction of propagation in the guide. The parameters governing this, such as the pitch of the grating, depend on the optical index of the guide and the chosen wavelength. For the coupler to be ideal, all light should be transmitted in a single order capable of propagating in the guide through total internal reflection. The distribution of energy among different orders can only be determined through rigorous electromagnetic simulation. Achieving this objective involves optimizing the depth and inclination angles of the grating through an optimization process. This requires both maximizing a specific order while simultaneously minimizing all others. Therefore, the numerical code must converge rapidly, ensuring that computed efficiencies, which should be as small as possible, are greater than the unavoidable calculation error. Considering grating (b) operating at a wavelength of $\lambda = 1.031$, our FMM results with 101 Floquet harmonics and 800 slices provide the following TE efficiencies: $Re = 9.3822(-4)$, $Te_{-1} = 9.524e(-1)$, $Te_0 = 1.7000(-3)$, and $Te_{+1} = 4.4521e(-2)$. Our focus is on examining the convergence of the $Te_{+1}$ efficiency concerning the number of slices for different numbers of Floquet harmonics. We define an error function as follows:

$$err(N, M) = \log(|Te_{+1}^{N,M} - Te_{+1}^{ref}|) \qquad (46)$$

where N is for the number of Floquet harmonics, M is the number of slices and $Te_{+1}^{ref}$ is the value obtained with N = 101 and M = 800 slices. Figure (5) depicts the error functions $err(10, M)$, $err(15, M)$, and $err(20, M)$ for M ranging from 8 to 200. It is observed that the error function reaches a minimum value whose location depends on N, suggesting a link between horizontal and vertical spatial resolution. Surprisingly, the error function does not exhibit significant dependence on N for a small number of slices. Regardless of the number of Floquet harmonics, M = 42 slices are sufficient to achieve an error function below -3. Specifically, with the couple (N, M) = (21, 42), the error function is $err = -3.81$. In a subsequent numerical experiment, we held N constant and varied the number of slices M. For each couple (N, M), we executed the code 100 times and measured the computation time. The results are presented in Figure (7), illustrating the computation time as a function
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Error in the calculation of the $Te_{+1}$ efficiency of grating (b) in figure. *N* is the number of Floquet harmonics, *M* is the number of slices. The grating is enlightened from vacuum under normal incidence. The parameters are $\lambda = 1.031, v_1 = 1, v_2 = 1.5$ and its end facets coordinates in arbitrary units are: $(x, z)$ are 0, 0), 0.75, 1.25, (1, 1.25), (0.72, 0), 1, 0)</figcaption></figure>

of the error function. The code was implemented in the Matlab language on a Dell XPS 13-9305 with an Intel i5 processor. The same analysis was conducted with our C-method code utilizing polynomial expansion and an oblique coordinate system. The inclination angle was chosen as the mean value of the inclination of the first and third facets, with an equal number of basis functions chosen for each facet. The error function, defined similarly to the FMM computation, now depends on the single parameter *dim*. It is observed that the error function also reaches a minimum, but the reason for this differs from the previous simulation with FMM. In the case of the C-method, it has been demonstrated that a too-large truncation number leads to ill-conditioned systems due to round-off errors when writing boundary conditions at the grating surface. This issue is illustrated by plotting an estimate of the reciprocal of the condition of the **S** matrix in the 1-norm, given by the function *rcond* in Matlab. If **S** is well-conditioned, *rcond(S)* is near 1.0; if **S** is badly conditioned, *rcond(S)* is close to 0. In the present simulation, reliable computed efficiencies are obtained only if *rcond(S)* is greater than $10^{-19}$. Taking this criterion into account, the reference value for the $Te_{+1}$ order was found to be $Te_{+1} = 4.4523e(-2)$. Choosing 9 basis elements per facet with the C-method resulted in an error function of *err* = -4.18. Using the same rule as for the FMM, we investigated the computation time as a function of the error, with the results also reported in Figure (7). The C-method with polynomial expansion outperforms the FMM. Additionally, it is important to note that the chosen test case is very favorable to the FMM, being the TE polarization with low index contrast. As the performance of the C-method is independent of materials or polarization, one may conclude that it will consistently outperform the FMM whenever applicable.

## 6. Conclusion

We formulated the C-method using an oblique coordinate system and implemented it for the first time with Legendre polynomial expansion (CPE). We compared it to the Fourier Modal Method (FMM) in a scenario that was initially favorable to the FMM. Remarkably, for the same level of precision, CPE significantly outperformed FMM in terms of computational speed. While no single method can efficiently address all diffraction problems related to gratings, CPE stands out as an excellent candidate, particularly for faceted profiles.

**Author Contributions:** Conceptualization, G.G. and K.E.; methodology, G.G.; software, G.G.; validation, G.G. and K.E.; formal analysis, G.G. and K.E.; writing—original draft preparation, G.G.; writing—review and editing, G.G. and K.E.; visualization, K.E.; All authors have read and agreed to the published version of the manuscript.

**Funding:** "This research received no external funding".

**Data Availability Statement:** Data sharing is not applicable to this article.
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Error in the calculation of the $Te_{+1}$ efficiency of grating (b) in figure. The grating is enlightened from vacuum under normal incidence. The parameters are $\lambda = 1.031, v_1 = 1, v_2 = 1.5$ and its end facets coordinates in arbitrary units are: (x, z) are 0,0), 0.75, 1.25, (1, 1.25), (0.72, 0), 1, 0)</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Computation time as a function of error with FMM and CPE. The computation concerns the $Te_{+1}$ transmitted order. The graphs plot 100 times the true computation time. For Fmm, the number of Floquet harmonics being fixed at 21, the error diminishes with the number of slices. For CPE the error diminishes with the number of basis elements.</figcaption></figure>
**Conflicts of Interest:** "The authors declare no conflict of interest." 369

**Appendix A** 370

In this section we derive the Legendre polynomial expansion of a plane wave on the qth facet of a faceted grating whose end-points coordinates of the qth facet are $(x_{q-1}, z_{q-1})$ and $(x_q, z_q)$ respectively. In a medium with a real optical index $v$, the y-component of a real diffracted order is such that: 371 372 373 374

$$ \phi^{\pm}(x, z) = \sum_{m} A_{m}^{\pm} \exp(\mp ik\rho_{m}z) \exp(-ika_{m}x), \quad (A1) $$

with $\alpha_m = v \sin(\theta) + m \frac{\lambda}{d}$. The exposants + and - are for forward and backward waves respectively. The real diffracted waves are such that $v^{(p)} \le \alpha_m \le v^{(p)}$. For a $\exp(i\omega t)$ time dependence, $\rho_m^{(p)}$ is defined as follows: 375 376 377

$$ \rho_m = \begin{cases} \sqrt{v^2 - \alpha_m^2} & \text{if } (v^2 - \alpha_m^2) \in \mathbb{R}^+ \\ -i\sqrt{\alpha_m^2 - v^2} & \text{if } (v^2 - \alpha_m^2) \in \mathbb{R}^- \end{cases} \quad (A2) $$

The facet is parameterized by functions $f(x^1)$ and $h(x^1)$ 378

$$ \begin{aligned} f(x^1) &= .5(x_q - x_{q-1})x^1 + .5(x_q + x_{q-1}) \\ h(x^1) &= .5(z_q - z_{q-1})x^1 + .5(z_q + z_{q-1}) \end{aligned} \quad (A3) $$

In the new coordinates $x^1, x^3$, the surface of the grating is the coordinate surface $x^3 = 0$. So $\phi^{\pm}$ is a function of the only $x^1$ coordinate which is obtained by setting $x = f(x^1)$ and $z = h(x^1)$ in (A1). So, adding an extra index $q$ to refer to the qth facet, the $x^1$ function to consider is: 379 380 381 382

$$ \phi_{m,q}^{\pm}(x^1) = \exp(-iY_q x^1) \exp(-i\Theta_q) \exp(-i\Gamma_q) \quad (A4) $$

with: 383

$$ \begin{aligned} Y_q &= 0.5k(\pm\rho_m(z_q - z_{q-1}) + \alpha_m(x_q - x_{q-1})) \\ \Theta_q &= 0.5k(\pm\rho_m(z_q + z_{q-1})) \\ \Gamma_q &= 0.5k\alpha_m(x_q + x_{q-1}) \end{aligned} \quad (A5) $$

Lets denote by $e^n(\alpha)$ the expansion coefficient of $\exp(-i\alpha x)$ on to the Legendre polynomials: 384

$$ \exp(-i\alpha x) = |e^n(\alpha)L_n(x)| \quad \text{where} \quad e^n(\alpha) = \frac{\langle\exp(-i\alpha x)L_m(x)\rangle}{\langle L_n(x)L_m(x)\rangle} \quad (A6) $$

The inner product $\langle\exp(-i\alpha x), L_n(x)\rangle$ has the following analytical expression 385

$$ \int_{-1}^{1} L_n(x) \exp(-i\alpha x^1) dx^1 = \begin{cases} (-i)^n \sqrt{\frac{2\pi}{\alpha}} J_{n+\frac{1}{2}}(\alpha) & \text{if } \alpha > 0 \\ (i)^n \sqrt{\frac{2\pi}{|\alpha|}} J_{n+\frac{1}{2}}(\alpha) & \text{if } \alpha < 0 \end{cases} \quad (A7) $$

where $J_{n+\frac{1}{2}}$ is a Bessel function of the first kind of order $n+1/2$. From Eq(A4) and Eq(A7), we may derive the sought expansion coefficients 386 387

$$ \phi_{m,q}^{\pm}(x^1) = |\phi_{m,q,l}^{l\pm} L_l\rangle \quad (A8) $$

$$ \phi_{m,q,l}^{l\pm} = e^n(Y_q) \exp(-i\Theta_q) \exp(-i\Gamma_q) \quad (A9) $$
Let us assume that $\phi$ is $E_2$. Then, in Cartesian coordinates, it is trivial to express $Z_0H_x$ and $Z_0H_z$ in terms of $E_2$:

$$ \begin{aligned} Z_0H_x &= \mp\rho_m E_2 \\ Z_0H_z &= \alpha_m E_2 \end{aligned} \qquad (A10) $$

from which we get $G_1$

$$ G_1 = \frac{df}{dx^1}(Z_0H_x) + \frac{dh}{dx^1}(Z_0H_z) \qquad (A11) $$

References

1. Song, Q.; Pigeon, Y.E.; Heggarty, K. Faceted gratings for an optical security feature. Appl. Opt. 2020, 59, 910–917. https://doi.org/10.1364/AO.378122.

2. Song, J.; Ding, J.F. Echelle diffraction grating demultiplexers with a single diffraction passband. Optics Communications 2010, 283, 537–541. https://doi.org/10.1016/j.optcom.2009.10.095.

3. Miller, J.M.; de Beaucoudrey, N.; Chavel, P.; Turunen, J.; Cambril, E. Design and fabrication of binary slanted surface-relief gratings for a planar optical interconnection. Appl. Opt. 1997, 36, 5717–5727. https://doi.org/10.1364/AO.36.005717.

4. Mattelin, M.A.; Radosavljevic, A.; Missinne, J.; Cuypers, D.; Steenberge, G.V. Design and fabrication of blazed gratings for a waveguide-type head mounted display. Opt. Express 2020, 28, 11175–11190. https://doi.org/10.1364/OE.384806.

5. Chandezon, J.; Raoult, G.; Maystre, D. A new theoretical method for diffraction gratings and its numerical application. Journal of Optics 1980, 11, 235.

6. Chandezon, J.; Dupuis, M.T.; Cornet, G.; Maystre, D. Multicoated gratings: a differential formalism applicable in the entire optical region. J. Opt. Soc. Am. 1982, 72, 839–846. https://doi.org/10.1364/JOSA.72.000839.

7. Moharam, M.G.; Gaylord, T.K. Rigorous coupled-wave analysis of metallic surface-relief gratings. J. Opt. Soc. Am. A 1986, 3, 1780–1787. https://doi.org/10.1364/JOSAA.3.001780.

8. Moharam, M.G.; Grann, E.B.; Pommet, D.A.; Gaylord, T.K. Formulation for stable and efficient implementation of the rigorous coupled-wave analysis of binary gratings. J. Opt. Soc. Am. A 1995, 12, 1068–1076. https://doi.org/10.1364/JOSAA.12.001068.

9. Lalanne, P.; Morris, G.M. Highly improved convergence of the coupled-wave method for TM polarization. JOSA A 1996, 13, 779–784.

10. Granet, G.; Guizal, B. Efficient implementation of the coupled-wave method for metallic lamellar gratings in TM polarization. J. Opt. Soc. Am. A 1996, 13, 1019–1023. https://doi.org/10.1364/JOSAA.13.001019.

11. Li, L. Use of Fourier series in the analysis of discontinuous periodic structures. J. Opt. Soc. Am. A 1996, 13, 1870–1876. https://doi.org/10.1364/JOSAA.13.001870.

12. Li, L.; Chandezon, J. Improvement of the coordinate transformation method for surface-relief gratings with sharp edges. J. Opt. Soc. Am. A 1996, 13, 2247–2255. https://doi.org/10.1364/JOSAA.13.002247.

13. Granet, G.; Chandezon, J.; Plumey, J.P.; Raniriharinosy, K. Reformulation of the coordinate transformation method through the concept of adaptive spatial resolution. Application to trapezoidal gratings. J. Opt. Soc. Am. A 2001, 18, 2102–2108. https://doi.org/10.1364/JOSAA.18.002102.

14. Morf, R.H. Exponentially convergent and numerically efficient solution of Maxwell’s equations for lamellar gratings. J. Opt. Soc. Am. A 1995, 12, 1043–1056. https://doi.org/10.1364/JOSAA.12.001043.

15. Edee, K. Modal method based on subsectional Gegenbauer polynomial expansion for lamellar gratings. J. Opt. Soc. Am. A 2011, 28, 2006–2013. https://doi.org/10.1364/JOSAA.28.002006.

16. Edee, K.; Fenniche, I.; Granet, G.; Guizal, B. Modal method based on subsectional Gegenbauer polynomial expansion for lamellar gratings: weighting function, convergence and stability. Prog. Electromagn. Res. 2013, 133, 13–35.

17. Randriamihaja, M.H.; Granet, G.; Edee, K.; Raniriharinosy, K. Polynomial modal analysis of lamellar diffraction gratings in conical mounting. J. Opt. Soc. Am. A 2016, 33, 1679–1686. https://doi.org/10.1364/JOSAA.33.001679.

18. Granet, G.; Randriamihaja, M.H.; Raniriharinosy, K. Polynomial modal analysis of slanted lamellar gratings. J. Opt. Soc. Am. A 2017, 34, 975–982. https://doi.org/10.1364/JOSAA.34.000975.

19. Edee, K.; Plumey, J.P. Numerical scheme for the modal method based on subsectional Gegenbauer polynomial expansion: application to biperiodic binary grating. J. Opt. Soc. Am. A 2015, 32, 402–410. https://doi.org/10.1364/JOSAA.32.000402.

20. Plumey, J.P.; Granet, G. Generalization of the coordinate transformation method with application to surface-relief gratings. J. Opt. Soc. Am. A 1999, 16, 508–516. https://doi.org/10.1364/JOSAA.16.000508.

21. Plumey, J.P.; Guizal, B.; Chandezon, J. Coordinate transformation method as applied to asymmetric gratings with vertical facets. J. Opt. Soc. Am. A 1997, 14, 610–617. https://doi.org/10.1364/JOSAA.14.000610.