# Analysis of Faceted Gratings Using C-Method and Polynomial Expansion

Gérard Granet, Kofi Sényo Edee

► To cite this version:

Gérard Granet, Kofi Sényo Edee. Analysis of Faceted Gratings Using C-Method and Polynomial Expansion. Photonics, 2024, 11 (3), pp.215. 10.3390/photonics11030215. [hal-04827027](https://hal.science/hal-04827027v1)

HAL Id: hal-04827027

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

Submitted on 6 Jan 2025

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

# Analysis of Faceted Gratings Using C-Method and Polynomial Expansion

Gérard Granet *,† ID and Kofi Edee †,D

Université Clermont Auvergne, Clermont Auvergne INP, CNRS, Institut Pascal, F-63000 Clermont-Ferrand, France; kofi.edee@uca.fr

* Correspondence: gerard.granet@uca.fr

† These authors contributed equally to this work.

**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 us 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. One very efficient and fast method for the simulation of the diffraction by optical gratings is the integral equation method which has many variants and is capable of tackling profile curves with corners in its recent versions [5,6]. In return, the theory behind it involves sophisticated mathematics. 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 [7,8] (the C method) and the coupled-wave method [9,10] (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

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

check for
updates

**Citation:** Granet, G.; Edee, K. Analysis of Faceted Gratings Using C-Method and Polynomial Expansion. *Photonics* **2024**, *11*, 215. https://doi.org/10.3390/photonics11030215

Received: 21 December 2023

Revised: 17 February 2024

Accepted: 22 February 2024

Published: 27 February 2024

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

**Copyright:** © 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).
slowly for profiles with sharp edges, while FMM is not well-suited for metal-dielectric
structures in the optical region. Fortunately, improvements achieved for one method can
generally be applied to the other. Factorization rules initially discovered [11,12] and later
rigorously established for lamellar gratings in the case of the FMM [13] have been suc-
cessfully applied to gratings with sharp edges [14]. 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 and solutions of Maxwell’s equations, it was shown that poly-
nomial expansion allows for the enforcement of boundary conditions at jump locations,
resulting in exponential convergence of the spectral operator [15–17]. 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 continuous at the jumps’ location, as outlined in [18]. 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 grat-
ings, including the case of overhanging gratings within the framework of the C method.
Plumey et al. [19] and Preist et al. [20] were the first to extend the C method to encompass
isotropic gratings with overhanging profiles or vertical facets. Plumey et al.’s approach
exploits the fact that only the derivative of the grating profile function, not the function
itself, appears in the eigenvalue problem. Similarly, Preist et al., as well as Li [21], utilize an
oblique coordinate system where the grating profile can be represented by a single-valued
periodic function. In a recent paper, Ming and Sun [22] adopted a similar approach to
Plumey et al. but formulated the C-method based on a simple change in variables instead of
tensor calculus as conducted by Li et al. [23] The approaches of Plumey et al. and Preist et al.
are fully equivalent, and in this paper, we utilize the latter. Therefore, the main difference
from Preist et al.’s work arises from the use of polynomial expansion.

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 polynomial 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 *λ*, angular
frequency *ω*, and vacuum wavenumber *k*. The wave vector is inclined at an angle *θ* to
the *Oz* axis. The time dependence is expressed by the factor *exp*(i*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_4.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 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 the 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 that enables such a kind of solution a “translation coordinate system”. For a cylindrical surface invariant along the $Oy$ axis, Plumey et al. [24]
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, see Figure 2. This results in an inclined (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_5.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^i J_j, \qquad (6) $$

where the exponent $i$ 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 & 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 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 the 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}
\xi^{ijk} \partial_j E_k &= -ik \sqrt{g} g^{ij}(Z_0 H_j) \\
\xi^{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{g}g^{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} \quad \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 with 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 operator [25]. 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 ith testing and jth 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 depends 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 the 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 the one hand and $G_1$ on the other hand should not share the same basis [18]. 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 us 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 us examine the components $F_2$ and $G_3$. In accordance with Edee and Plumey [18], 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}_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,m}$ 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:

$$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 to the $q$th sub-domain and is a matrix with $Q$ rows and $M_q$ columns. Its non-zero elements are given below:

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

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

The newly defined basis being constructed as a linear combination of Legendre polynomials it is convenient to introduce the change in 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 a manner similar to $\tilde{B}_n$ but with coefficients $a_{M_{k,q}}$ being equal to zero. Numerical experiments have shown that the latter basis is also convenient to serve as the basis of the testing functions. As conducted for the $\tilde{B}_m$, we define the change in basis matrix $K$ whose elements are:

$$K_{mn} = \begin{cases} \delta_{mn} & \text{if } m = M_q \\ 0 & \text{if } m \neq M_q \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}} \quad \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 be 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}(\mathbf{d}) \tilde{K} \qquad (27)$$

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

### 4.3. Algebraic Eigenequation

We are now able to transform into an algebraic eigensystem of equations. Since they are continuous, $F_2$ and $G_3$ field components are expanded onto a $\tilde{B}_m$ basis whereas $G_1$ is expanded onto a $B_m$ basis which serves also as a test basis. Using a bracket 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 eigensystem:
$$ \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}. \qquad (32) $$

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

$$ \mathbf{D} = \frac{i}{k} (\mathbf{K}^t \mathbf{K})^{-1} \tilde{\mathbf{D}}. \qquad (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}^{2\text{dim}} 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}. \qquad (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_6.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. (1) and (2) refer to medium 1 and medium 2 respectively</figcaption></figure>

Solving the matrix eigen equation (31) for each medium gives us two sets of eigenvectors from which we easily reconstruct 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}^{2\text{dim}} 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}. \qquad (35) $$

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

$$ \sum_{q=1}^{2\text{dim}} A_q^{(1)} \begin{bmatrix} \frac{1}{\nu} |F_{2q}^{(1)m}\rangle \\ |G_{1q}^{(1)m}\rangle \end{bmatrix} = \sum_{q=1}^{2\text{dim}} A_q^{(2)} \begin{bmatrix} \frac{1}{\nu} |F_{2q}^{(2)m}\rangle \\ |G_{1q}^{(2)m}\rangle \end{bmatrix}. \qquad (36) $$
To go further, we need to sort the eigensolution of the 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} = S \begin{bmatrix} A_q^{(1)-} \\ A_q^{(2)+} \end{bmatrix}, \qquad (37) $$

with

$$ 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 and backward waves are those with a computed eigenvalue with a negative and positive imaginary part, respectively. 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 a 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 us to distinguish among the computed real eigenvalues those that correspond to forward waves or backward waves. For the diffracted waves such as $-\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 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_q^{(1)+}$ and $A_q^{(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{aligned} R_n^{te} &= \frac{|A_n|^2 \rho_n^{(1)}}{\gamma_0^{(1)}} & \text{for reflected TE polarized orders} \\ T_n^{te} &= \frac{|A_n|^2 \rho_n^{(2)}}{\gamma_0^{(1)}} & \text{for transmitted TE polarized orders} \\ R_m^{tm} &= \frac{|A_n|^2 \rho_n^{(1)}}{\gamma_0^{(1)}} & \text{for reflected TM polarized orders} \\ 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} \end{aligned} \quad (45) $$

# 5. Results

## 5.1. Validation by Comparison with Published Data

In this section, we first validate our code by comparing our results with previously published data by Plumey et al. [19], which were obtained through another implementation of the C-method and the Fourier Modal Method (FMM). The investigation focuses on two gratings displayed in Figure 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 the polynomial expansion, we incorporate only 15 basis elements per facet. Results are shown in Table 1. 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. Next, for comparison purposes, we consider a more recent paper by Ming et al. [22]. The profile shape corresponds to case (c) in Figure 4. The grating is illuminated from the vacuum under normal incidence with a wavelength $\lambda = 1/1.5$. Results are shown in Table 2. The optical index of the substrate is $\nu_2 = 1.5$. It is observed that the agreement between the two variants of the C-method is excellent. It should be noted that, regardless of polarization, our results were obtained using 60 polynomials in the field component expansion, while Ming et al. mention 229 and 449 Fourier harmonics for TE and TM polarization, respectively.
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Overhanging gratings. The endpoint coordinates of the facets of the grating (a) are (0,0), (1.1,1), and (1,1). Gratings (b) and (c) have similar profile shapes with different parameters. The endpoint coordinates of the facets of the grating (b) are (0,0), 0.75, 1.25, (1, 1.25), (0.72,0), (1,0) whereas the endpoint coordinates of the facets of the grating (c) are (0,0), (3/8,1), (5/8,1), (4/8,0), (1,0). All units are arbitrary.</figcaption></figure>

<table><caption>**Table 1.** Comparison between FMM and C-method by Plumey et al. [19] on one hand and our implementation of C-method in an oblique coordinate system on the other hand. The exponent (a), is for Plumey's results, taken from Table 3 in [19]. The profile shape is that of grating (a) in Figure 4. The optical indices are such that $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 <sup>(a)</sup></th><th colspan="2">FMM <sup>(a)</sup></th><th colspan="2">CPE</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>

<table><caption>**Table 2.** Comparison between C-method by Ming et al [22] on the one hand and our implementation of C-method (CPE, for C-method by polynomial expansion) in an oblique coordinate system on the other hand. The exponent (a), is for Mings's results, taken from Table 3 in [22] The profile shape is that of grating (c) in Figure 4. The optical index is $v_2 = 1.5$. The incident wave parameters are $\theta = 0^\circ$ and $\lambda = 1/1.5$.</caption><thead><tr><th>Orders</th><th>C-Method <sup>(a)</sup></th><th>CPE</th></tr></thead><tbody><tr><td colspan="3"><strong>TE polarization</strong></td></tr><tr><td>R<sub>1</sub></td><td>0.1156 (-1)</td><td>0.01162</td></tr><tr><td>R<sub>0</sub></td><td>0.3981 (-2)</td><td>0.00397</td></tr><tr><td>R<sub>1</sub></td><td>0.2673 (-3)</td><td>0.00270</td></tr><tr><td>T<sub>-2</sub></td><td>0.9877 (-1)</td><td>0.98399</td></tr><tr><td>T<sub>-1</sub></td><td>0.4771</td><td>0.47659</td></tr><tr><td>T<sub>0</sub></td><td>0.9975 (-1)</td><td>0.09944</td></tr><tr><td>T<sub>1</sub></td><td>0.2685</td><td>0.26885</td></tr><tr><td>T<sub>2</sub></td><td>0.3831 (-1)</td><td>0.03428</td></tr></tbody></table>
<table><caption>**Table 2.** *Cont.*</caption><thead><tr><th>Orders</th><th>C-Method <sup>(a)</sup></th><th>CPE</th></tr></thead><tbody><tr><td>TM polarization</td><td></td><td></td></tr><tr><td>R<sub>-1</sub></td><td>0.5258 (-2)</td><td>0.00523</td></tr><tr><td>R<sub>0</sub></td><td>0.3197 (-2)</td><td>0.00320</td></tr><tr><td>R<sub>1</sub></td><td>0.5832 (-3)</td><td>0.00058</td></tr><tr><td>T<sub>-2</sub></td><td>0.3419 (-1)</td><td>0.03405</td></tr><tr><td>T<sub>-1</sub></td><td>0.7470</td><td>0.74813</td></tr><tr><td>T<sub>0</sub></td><td>0.4760 (-1)</td><td>0.04749</td></tr><tr><td>T<sub>1</sub></td><td>0.1502</td><td>0.15022</td></tr><tr><td>T<sub>2</sub></td><td>0.1240 (-1)</td><td>0.01239</td></tr></tbody></table>

## 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 a vacuum, and the waveguide is made of glass with a refractive index $v = 1.53$. The coupler's operation requires a single order to be reflected in the 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 6, illustrating the computation time as a function of the error function. The code was implemented in the Matlab language (version 2021b) 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<sup>-19</sup>. Taking this criterion into account, the reference value for the Te<sub>+1</sub> order was found to be Te<sub>+1</sub> = 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.

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Error in the calculation of the $Te_{+1}$ efficiency of grating (b) in Figure 4. N is the number of Floquet harmonics, M is the number of slices. The grating is enlightened from the 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 6.** Error in the calculation of the $Te_{+1}$ efficiency of grating (b) in Figure 4. 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_10.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>

## 6. Conclusions

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.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data are contained within the article.

**Conflicts of Interest:** The authors declare no conflicts of interest.

## Appendix A

In this section we derive the Legendre polynomial expansion of a plane wave on the $q$th facet of a faceted grating whose end-points coordinates of the $q$th 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:

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

with $\alpha_m = v \sin(\theta) + m \frac{\lambda}{d}$. The exponents + and - are for forward and backward waves respectively. The real diffracted waves are such that $v^{(p)} \le \alpha_m \le v^{(p)}$. For an exp($i\omega t$ time dependence, $\rho_m^{(p)}$ is defined as follows:
$$ \rho_m = \begin{cases} \sqrt{\nu^2 - \alpha_m^2} & \text{if } (\nu^2 - \alpha_m^2) \in \mathbb{R}^+ \\ -i\sqrt{\alpha_m^2 - \nu^2} & \text{if } (\nu^2 - \alpha_m^2) \in \mathbb{R}^- \end{cases} \qquad (A2) $$

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

$$ \begin{aligned} f(x^1) &= 0.5(x_q - x_{q-1})x^1 + 0.5(x_q + x_{q-1}) \\ h(x^1) &= 0.5(z_q - z_{q-1})x^1 + 0.5(z_q + z_{q-1}) \end{aligned} \qquad (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 $q$th facet, the $x^1$ function to consider is:

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

with:

$$ \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} \qquad (A5) $$

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

$$ \exp(-i\alpha x) = |e^n(\alpha)L_n\rangle \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} \qquad (A6) $$

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

$$ \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} \qquad (A7) $$

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

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

$$ \phi_{m,q,l}^{l\pm} = e^n(Y_q) \exp(-i\Theta_q) \exp(-i\Gamma_q) \qquad (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. [CrossRef]

2. Song, J.; Ding, J.F. Echelle diffraction grating demultiplexers with a single diffraction passband. Opt. Commun. 2010, 283, 537–541. [CrossRef]

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. [CrossRef] [PubMed]

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. [CrossRef] [PubMed]

5. Popov, E. (Ed.) Gratings: Theory and Numeric Applications, 2nd ed.; AMU, CNRS, Institut Fresnel: Marseille, France, 2014.
6. Bernd, H.; Kleemann, A.M.; Wyrowski, F. Integral equation method with parametrization of grating profile theory and experiments. *J. Mod. Opt.* **1996**, *43*, 1323–1349. [CrossRef]

7. Chandezon, J.; Raoult, G.; Maystre, D. A new theoretical method for diffraction gratings and its numerical application. *J. Opt.* **1980**, *11*, 235. [CrossRef]

8. 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. [CrossRef]

9. Moharam, M.G.; Gaylord, T.K. Rigorous coupled-wave analysis of metallic surface-relief gratings. *J. Opt. Soc. Am. A* **1986**, *3*, 1780–1787. [CrossRef]

10. 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. [CrossRef]

11. Lalanne, P.; Morris, G.M. Highly improved convergence of the coupled-wave method for TM polarization. *J. Opt. Soc. Am. A* **1996**, *13*, 779–784. [CrossRef]

12. 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. [CrossRef]

13. Li, L. Use of Fourier series in the analysis of discontinuous periodic structures. *J. Opt. Soc. Am. A* **1996**, *13*, 1870–1876. [CrossRef]

14. 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. [CrossRef]

15. 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. [CrossRef]

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. [CrossRef]

17. Granet, G.; Randriamihaja, M.H.; Raniriharinosy, K. Polynomial modal analysis of slanted lamellar gratings. *J. Opt. Soc. Am. A* **2017**, *34*, 975–982. [CrossRef]

18. 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. [CrossRef]

19. 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. [CrossRef]

20. Preist, T.W.; Harris, J.B.; Wanstall, N.P.; Sambles, J.R. Optical response of blazed and overhanging gratings using oblique chandezon transformations. *J. Mod. Opt.* **1997**, *44*, 1073–1080. [CrossRef]

21. Li, L. Oblique-coordinate-system-based Chandezon method for modeling one-dimensionally periodic, multilayer, inhomogeneous, anisotropic gratings. *J. Opt. Soc. Am. A* **1999**, *16*, 2521–2531. [CrossRef]

22. Ming, X.; Sun, L. Simple reformulation of the coordinate transformation method for gratings with a vertical facet or overhanging profile. *Appl. Opt.* **2021**, *60*, 4305–4314. [CrossRef] [PubMed]

23. Li, L.; Chandezon, J.; Granet, G.; Plumey, J.P. Rigorous and efficient grating-analysis method made easy for optical engineers. *Appl. Opt.* **1999**, *38*, 304–313. [CrossRef] [PubMed]

24. 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. [CrossRef]

25. Harrington, R.; Antennas, I.; Society, P. *Field Computation by Moment Methods*; IEEE Press Series on Electromagnetic Waves; IEEE: New York, NY, USA, 1996.

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.