# Topology optimization of blazed gratings under conical incidence

Simon Ans  Aix Marseille Univ, CNRS, CNES, LAM, Marseille, France  Aix
Marseille Univ, CNRS, Centrale Med, Institut Fresnel, Marseille, France
Frédéric Zamkotsian  Aix Marseille Univ, CNRS, CNES, LAM, Marseille, France
Guillaume Demésy  Aix Marseille Univ, CNRS, Centrale Med, Institut Fresnel,
Marseille, France

###### Abstract

A topology optimization method is presented and applied to a blazed
diffraction grating in reflection under conical incidence. This type of
gratings is meant to disperse the incident light on one particular diffraction
order and this property is fundamental in spectroscopy. Conventionally, a
blazed metallic grating is made of a sawtooth profile designed to work with
the $\pm 1^{\text{st}}$ diffraction order in reflection. In this paper, we
question this intuitive triangular pattern and look for optimal opto-geometric
characteristics using topology optimization based on Finite Element modelling
of Maxwell’s equations. In practical contexts, the grating geometry is mono-
periodic but it is enlightened by a 3D plane wave with a wavevector outside of
the plane of invariance. Consequently, this study deals with the resolution of
a direct and inverse problem using the Finite Element Method in this
intermediate state between 2D and 3D: the so-called conical incidence. A
multi-wavelength objective is used in order to obtain a broadband blazed
effect. Finally, several numerical experiments are detailed. The results show
that it is possible to reach a $98\%$ diffraction efficiency on the
$-$1${}^{\text{st}}$ diffraction order if the optimization is performed on a
single wavelength, and that the reflection integrated over the [400,1500] nm
wavelength range can be 29% higher in absolute, 56% in relative, than that of
the sawtooth blazed grating when using a multi-wavelength optimization
criterion (from 52% to 81%).

##### Keywords:

Maxwell’s equations, Conical incidence, Finite Element Method, Topology
optimization, Blazed gratings, Metasurface, Broadband optimization.

##  1 Introduction

Topology optimization [1] is a powerful modelling tool allowing to adapt the
design of devices with respect to a given performance target. Since its
introduction in the late nineties, it has been substantially improved and
applied to many areas of physics and engineering. In this method, a given
_design space_ $\Omega_{d}$ is chosen as a subset of the whole computational
domain $\Omega$ and a _design variable_ $\rho$ (also called density field in
the literature) is defined over $\Omega_{d}$. It takes values in $[0,1]$, 0
representing the presence of a chosen material, 1 the presence of another one.
The regularity of this variable upon $\Omega_{d}$ is then ensured by
supplementary constraints of binarization and connectedness. This method
allows to keep the same mesh throughout the optimization process. It finds
applications the numerous domains of physics and engineering where Partial
Differential Equations (PDEs) are implied. A non exhaustive list includes
solid mechanics and acoustics [2, 3, 4, 5], fluid mechanics [6, 7], electro-
mechanics [8, 9] and of course photonics [10, 11, 12, 13].

Topology optimization is often compared with shape optimization [14], which
has a very similar goal. However, in shape optimization, the Degrees of
Freedom (DoFs) are linked with the boundaries of the geometry which has to
remain consistent with the Finite Element Method (FEM). Thus, a major
difference with topology optimization is that it does not modify the topology
of the structure, which implies some constraints. For instance the boundaries
must always have at least a Lipschitz regularity. From a practical point of
view, a remeshing step at each iteration is therefore necessary [15, 16].
Shape optimization can and has been used to optimize blazed gratings for a
long time [17].

The blazed grating is the fundamental optical component of spectrometers. It
is expected to reflect most of the incoming light in one particular
diffraction order. Parametric optimizations over established geometries
allowed to optimize blazed gratings over a quite large range of frequencies
[18, 19, 20]. Typically, the blazed efficiency reaches either 90% for a
particular wavelength or lies between 20% and 60% in the visible/NIR/SWIR
frequency range (from 0.4 to 2.5$\,\mu$m).

More recently, impressive performance breakthroughs have been made with the
concept of metasurfaces coupled with advanced optimization techniques.
Heuristic optimizations of metasurfaces based on the lateral phase shift
induced by a blazed grating [21, 22] can in principle be outperformed by
topology optimization. A noticeable example is the open-source optimization
repository MetaNet provided by Fan _et al_ in Ref. [23]. This repository is
aimed at designing mono- and bi-periodic metasurface gratings (also called
metagratings) with respect to a deflection objective using topology
optimization (illustrated in Ref. [24] for instance), enabling to create a
metasurface database. The electromagnetic modelling part of this database
relies on the Rigorous Coupled Wave Analysis (RCWA), which restricts the
design to oblique edges. A more general tool is presented by Vial _et al_ in
Ref. [25] using Finite elements for the scalar case and auto-differentiation.
More general inverse design analyses for metasurfaces such as Ref. [26, 27,
28] include topology optimization.

In this paper, we report on the design of a blazed grating in reflection under
conical incidence in the visible and Near-InfraRed (NIR) wavelength range
[400,1500] nm, using the topology optimization on a model based on the
Maxwell’s equations. To our knowledge, this is the best performance ever
numerically obtained for a design of broadband blazed metasurface, reaching a
reflection averaged on this spectral range of 81%. Another goal of this paper
is also to clarify all the derivations leading to the Jacobian of the cost
function, which are rather intricate in the periodic case under conical
incidence. Finally, this paper is aimed at providing to the community an open-
source framework available
online111https://gitlab.onelab.info/doc/models/-/tree/master/DiffractionGratingsTopOpt
based on the FEM open-source suite Gmsh/GetDP [29, 30] and on the Globally
Convergent Method of Moving Asymptotes (GCMMA) optimization method [31].

The paper is organized as follows. The physical problem is described in
section 2, along with a quick overview of metasurfaces. The section 3 deals
with the direct problem which allows to compute the response of a mono-
periodic grating, enlightened by a 3D arbitrary plane wave, _i.e._ under
conical incidence. The latter is often overlooked although it constitutes a
truly valuable intermediate case, being more general than 2D and
computationally way lighter than 3D. The optimization problem for the design
of a blazed metasurface is introduced in section 4. This problem leads to the
need to compute the Jacobian of the target function, using results detailed in
section 5. The section 6 deals with the discretization aspects in order to
implement the optimization problem. Lastly, as the entire optimization process
and implementation are set out and clarified, numerical examples are shown to
illustrate the use of this method in section 7. These examples illustrate
various ways to beat the efficiencies of a classical sawtooth blazed grating
over a wide range of wavelengths. Moreover two supplementary supports are
provided. First, a technical support presents some proofs and details (see the
Appendix of this document). They are not necessary for the understanding of
the work and the results, but provide a more global scope of the calculations
that led to them and can be useful to the reader interested in modifying the
cost function. Secondly, the code that supplied all the results presented in
section 7 is available.

##  2 Problem description

###  2.1 Blazed nanostructured grating

In photonics, a diffraction grating is a periodic structure that deflects an
incident electromagnetic radiation to several directions, called the
_diffraction orders_ which depend on the period, the wavelength and the angle
of incidence. The wavelength dependence of this phenomenon justifies the
expression _dispersion of light_. A diffraction efficiency, ratio of the power
deviated in one order by the total incident power, is associated to each
diffraction order. This property leads to many applications, as shown by the
variety of articles cited in the introduction of this paper.

The purpose of a blazed grating is to maximize the diffraction of the
electromagnetic field on one particular diffraction order $n$ (typically the
$\pm 1$ order). In this article, a blazed grating in reflection on the $-1$
order is studied. However the demonstration to design it for transmission
and/or another order remains the same.

The most common reflective blazed gratings are made of triangles covered by a
thin layer of a reflective metal (often silver or gold, depending on the
targeted wavelength range) [32], as illustrated on Fig. 1a. This kind of
gratings have good performances on one octave of frequencies, but it would be
better to extend the blazing effect of the structure to at least two octaves
(visible and near infrared). Their historical design rules are particularly
simple, based on mere angular considerations as will be detailed later.

Figure 1: Two examples of mono-periodic structures. Their period is denoted by
$d$. The grey color represents a metal (_e.g._ silver Ag) and the blue color
represents a dielectric material (_e.g._ silica SiO2). (a) Sawtooth blazed
grating in reflection. The angle $\alpha$ is called the blaze angle. (b) Mono-
periodic metasurface grating.

We propose in this paper to study nanostructured gratings, the so-called
metasurfaces, as sketched on Fig. 1b. The proposed base structure consists of
an etched dielectric layer deposited on a reflective metal substrate. These
changes with respect to the sawtooth blazed grating presents considerable
theoretical advantages. While the triangles have only 1 DoF (the blaze angle
$\alpha$, see Fig. 1a), the flexibility of the structure is now tremendously
increased. The thickness of the layer, the dielectric materials considered,
the geometry of the pattern are parameters that become additional DoFs in
order to optimize the response of the structure.

In this article, the design of a mono-periodic metasurface is discussed,
however the incident field is considered to be 3D. This intermediate case
between the scalar 2D case and the full 3D vector case is called the _conical
case_. It is indeed useful since it allows to model the response of mono-
periodic gratings without any limitation on the incident light wave vector or
polarization, while keeping the calculation time far lower than for a full 3D
case.

The physical problem describing the response of such a grating is called the
_direct problem_ and its resolution using the FEM is detailed in the next part
of this section.

###  2.2 Computational domain and design space

Let consider a period $d$ of a mono-periodic grating, considered invariant
along the $z$ axis, as shown on Fig. 2. The numerical domain is denoted by
$\Omega$ (the dashed rectangle in Fig. 2). A point in this domain is denoted
by $\mathbf{r}=(x,y)$, whereas a point in the whole 3D space is denoted by
$\mathbf{x}=(x,y,z)$. The grating is enlightened by an incident linearly
polarized plane wave $\mathbf{E}^{\text{inc}}$ with a freespace wavelength
$\lambda$. The angles of incidence $\theta_{i}$, $\varphi_{i}$ define its
wavevector and $\psi_{i}$ denotes the polarization angle as shown on the
right-hand side of Fig. 2.

Figure 2: Description of the numerical domain for the conical case. The
incident electric field $\mathbf{E}^{\text{inc}}$ is defined by its wave
vector, depending on the incidence angles $\theta_{i}$, $\varphi_{i}$ and
$\psi_{i}$. For the conical case, a 2D section is considered, split into the
PMLs (violet regions), the superstrate $\mathcal{S}^{+}$, the substrate
$\mathcal{S}^{-}$ and the design region $\Omega_{d}$. Each region of the
domain $\Omega$ is defined by its permittivity and permeability tensors
fields, respectively $\underline{\underline{\varepsilon_{r}}}$ and
$\underline{\underline{\mu_{r}}}$. The (isotropic) relative permittivity
function of $\Omega_{d}$ is parametrized by a design field $\rho$.

The design region $\Omega_{d}$ (with shades of blue in $\Omega$) is deposited
on the substrate $\mathcal{S}^{-}$ (dark grey region) with isotropic relative
permittivity $\varepsilon_{r}^{-}(\lambda)$ and seats underneath the
superstrate $\mathcal{S}^{+}$ (transparent light grey region) with isotropic
relative permittivity $\varepsilon_{r}^{+}(\lambda)$ (typically 1 for
air/vacuum). The definition of the relative permittivity in the design region
is the keystone of the topology optimization. It can take any real value
between two relative permittivities $\varepsilon_{r,1}$ and
$\varepsilon_{r,2}$ by defining a density field
$\rho:\Omega_{d}\ni\mathbf{r}\mapsto\rho(\mathbf{r})\in[0,1]$. Commonly, the
SIMP method [33, 34] corresponds to the following bijection:

$\varepsilon_{r}^{d}(\lambda,\mathbf{r})=(\varepsilon_{r,2}(\lambda)-\varepsilon_{r,1}(\lambda))\rho(\mathbf{r})+\varepsilon_{r,1}(\lambda).$
(1)

The typical values taken by $\varepsilon_{r,1}$ and $\varepsilon_{r,2}$ are
respectively $\varepsilon_{r}^{+}$ and $\varepsilon_{r}^{\text{diel}}$, the
relative permittivity of a given dielectric. The upper and lower boundaries of
the numerical domain are then completed by Perfectly Matched Layers (PMLs, in
violet on Fig. 2) which permittivities are known tensors noted
$\underline{\underline{\varepsilon_{r}}}^{\text{PML}^{\pm}}$ [35]. The PMLs
have a magnetic permeability tensor denoted by
$\underline{\underline{\mu_{r}}}^{\text{PML}^{\pm}}$, whereas all the other
bulk materials in the domain have a relative permeability of 1. Eventually,
Bloch (quasi-periodic) conditions [36, 37] are imposed on the left and right
boundaries of the numerical domain.

##  3 Direct problem

###  3.1 Mono-periodic grating under conical incidence

Let introduce the incident field
$\mathbf{E}^{\text{inc}}=\mathbf{A_{e}}\exp(i\mathbf{k}^{+}_{\downarrow}\cdot\mathbf{x}),\,\,\forall\mathbf{x}\in\mathbb{R}^{3}$,
where $\mathbf{A}_{e}$ is the polarization of the electric field [38] defined
by

$\mathbf{A}_{e}=A_{e}\begin{pmatrix}\cos\psi_{i}\cos\theta_{i}\sin\varphi_{i}+\sin\psi_{i}\cos\varphi_{i}\\\
-\cos\psi_{i}\sin\theta_{i}\\\
\cos\psi_{i}\cos\theta_{i}\cos\varphi_{i}-\sin\psi_{i}\sin\varphi_{i}\end{pmatrix},$
(2)

with $A_{e}$ the amplitude of the electric field, and
$\mathbf{k}^{+}_{\downarrow}$ the wave vector that depends on the angles of
incidence (cf. conventions detailed on Fig. 2) and the wavelength. More
precisely it is the vector
$\mathbf{k}^{+}_{\downarrow}=(\alpha,\beta^{+},\gamma)^{\top}$ where, defining
$\omega=2\pi c/\lambda$ and $k_{0}=\omega/c$,

$\alpha=-k^{+}\sin\theta_{i}\sin\varphi_{i},\quad\beta^{+}=-k^{+}\cos\theta_{i}\quad\text{and}\quad\gamma=-k^{+}\sin\theta_{i}\cos\varphi_{i},$

denoting $k^{+}=k_{0}\sqrt{\varepsilon_{r}^{+}\mu_{r}^{+}}$. In the same way,
the wave vectors $\mathbf{k}^{+}_{\uparrow}=(\alpha,-\beta^{+},\gamma)^{\top}$
and $\mathbf{k}^{-}_{\downarrow}=(\alpha,\beta^{-},\gamma)^{\top}$ are
defined, with
$\beta^{-}=-\sqrt{k_{0}^{2}\varepsilon_{r}^{-}\mu_{r}^{-}-\alpha^{2}-\gamma^{2}}$.
These wave vectors are useful to define some intermediate electric fields
introduced below for the scattered field formulation.

When the grating presents one axis of invariance $(Oz)$ while enlightened by a
3D vector plane wave, the following traditional ansatz [39] is considered:

$\mathbf{E}_{\text{3D}}(\mathbf{x})=\mathbf{E}_{\text{3D}}(x,y,z)=\mathbf{E}(x,y)e^{i\gamma
z}=\mathbf{E}(\mathbf{r})e^{i\gamma z},$ (3)

which allows to reduce the unknown field to a 2D vector field $\mathbf{E}$.
This 2D field $\mathbf{E}$ can be split into its tangential component
$\mathbf{E}_{t}(\mathbf{r})=E_{x}(\mathbf{r})\hat{\mathbf{x}}+E_{y}(\mathbf{r})\hat{\mathbf{y}}$
and longitudinal component $E_{\ell}(\mathbf{r})\hat{\mathbf{z}}$ which is
continuous by construction. The following 2D (called tangential) operators
[39] are then introduced:

$\displaystyle\operatorname{\mathbf{grad}}_{t}E$ $\displaystyle=\frac{\partial
E}{\partial x}\hat{\mathbf{x}}+\frac{\partial E}{\partial y}\hat{\mathbf{y}},$
(4) $\displaystyle\operatorname{\mathbf{curl}}_{t}\mathbf{E}$
$\displaystyle=\Big{(}\frac{\partial E_{x}}{\partial y}-\frac{\partial
E_{y}}{\partial x}\Big{)}\hat{\mathbf{z}},$

while the $\operatorname{\mathbf{curl}}$ operator can be conveniently replaced
by the $\operatorname{\mathbf{curl}}_{\gamma}$ operator such that:

$\operatorname{\mathbf{curl}}_{\gamma}\mathbf{E}(\mathbf{r})=\operatorname{\mathbf{curl}}\left(\mathbf{E}(\mathbf{r})e^{i\gamma
z}\right)e^{-i\gamma z}.$ (5)

We are looking for solutions of finite energy, which means that the 2D vector
field $\mathbf{E}$ belongs to the functional space
$\bm{\mathcal{H}}(\operatorname{\mathbf{curl}}_{\gamma},\mathbb{R}^{3},e^{i\alpha
d})$ which corresponds to the functions of the Sobolev space
$\bm{\mathcal{H}}(\operatorname{\mathbf{curl}}_{\gamma},\mathbb{R}^{3})$ that
are pseudo-periodic with a factor $\alpha$ and a period $d$ along the $x$ axis
(a property coming from the Floquet-Bloch theorem). Eliminating the magnetic
field in Maxwell’s equations in the harmonic regime leads to the direct
problem on the total field $\mathbf{E}^{\text{tot}}$ [38]:

Find
$\mathbf{E}^{\text{tot}}\in\bm{\mathcal{H}}(\operatorname{\mathbf{curl}}_{\gamma},\mathbb{R}^{3},e^{i\alpha
d})$ such that $\forall\mathbf{x}\in\mathbb{R}^{3}$ with
$\omega\in\mathbb{R}^{+}$ given,

$\left\\{\begin{array}[]{llc}\operatorname{\mathbf{curl}}_{\gamma}(\underline{\underline{\mu_{r}}}^{-1}\operatorname{\mathbf{curl}}_{\gamma}\mathbf{E}^{\text{tot}})-k_{0}^{2}\underline{\underline{\varepsilon_{r}}}\mathbf{E}^{\text{tot}}=\mathbf{0},\\\
\text{such that an associated diffracted field $\mathbf{E}^{d}$ satisfies an
outgoing condition.}\end{array}\right.$ (6)

This problem is called the conical Helmholtz propagation equation and has a
unique solution. Note that it writes exactly the same for a given bi-periodic
grating, replacing the 2D $\operatorname{\mathbf{curl}}_{\gamma}$ operator by
the 3D $\operatorname{\mathbf{curl}}$ operator. The outgoing condition
satisfied by the diffracted field remains to be clarified.

###  3.2 Annex problem for the scattered field formulation

For the direct problem, we choose to work with a scattered field formulation
and introduce to that extent the following _annex problem_ [37] made a single
interface between the superstrate and substrate, which is the same as the
original problem removing the design region. The annex problem is thus
characterized by the following relative permittivity and permeability fields:

$\displaystyle\underline{\underline{\varepsilon_{r,a}}}(\lambda,\mathbf{r})$
$\displaystyle=\underline{\underline{\varepsilon_{r}}}(\lambda,\mathbf{r})-\varepsilon_{r}^{d}(\lambda,\mathbf{r})\mathbb{1}_{\Omega_{d}}(\mathbf{r})$
(7) $\displaystyle\underline{\underline{\mu_{r,a}}}(\mathbf{r})$
$\displaystyle=\underline{\underline{\mu_{r}}}(\mathbf{r}),$

noting $\mathbb{1}_{D}$ the characteristic function of a set $D$. This annex
problem can be easily solved making use of the Fresnel coefficients in the
appropriate orthonormal basis
$(\hat{\mathbf{p}},\hat{\mathbf{s}},\mathbf{k}^{+}_{\downarrow}/k^{+})$ (see
Fig. 2), which basically extends the notion of TE/TM polarization. More
explicitly, the Fresnel coefficients are

$\displaystyle r_{s}$
$\displaystyle=\frac{\beta^{+}-\beta^{-}}{\beta^{+}+\beta^{-}},\quad\quad\quad$
$\displaystyle t_{s}=\frac{2\beta^{+}}{\beta^{+}+\beta^{-}},$ (8)
$\displaystyle r_{p}$
$\displaystyle=\frac{\beta^{+}\varepsilon_{r}^{-}-\beta^{-}\varepsilon_{r}^{+}}{\beta^{+}\varepsilon_{r}^{-}+\beta^{-}\varepsilon_{r}^{+}},\quad\quad\quad$
$\displaystyle
t_{p}=\frac{2\beta^{+}\varepsilon_{r}^{-}}{\beta^{+}\varepsilon_{r}^{-}+\beta^{-}\varepsilon_{r}^{+}}\cdot$

Hence, denoting
$Z^{+}=\sqrt{\dfrac{\mu_{0}\mu_{r}^{+}}{\varepsilon_{0}\varepsilon_{r}^{+}}}$
the impedance of the superstrate, the fully $\hat{\mathbf{s}}$-polarized
electric and $\hat{\mathbf{s}}$-polarized magnetic fields write respectively:

$\left\\{\begin{array}[]{llc}\mathbf{E}^{\text{inc}}_{\hat{\mathbf{s}}}&=&\exp(i\mathbf{k}^{+}_{\downarrow}\cdot\mathbf{x})\,\hat{\mathbf{s}}\\\\[2.84526pt]
\mathbf{E}^{r}_{\hat{\mathbf{s}}}&=r_{s}&\exp(i\mathbf{k}^{+}_{\uparrow}\cdot\mathbf{x})\,\hat{\mathbf{s}}\\\\[2.84526pt]
\mathbf{E}^{t}_{\hat{\mathbf{s}}}&=t_{s}&\exp(i\mathbf{k}^{-}_{\downarrow}\cdot\mathbf{x})\,\hat{\mathbf{s}}\end{array}\right.\quad\text{and}\quad\left\\{\begin{array}[]{llc}\mathbf{H}^{\text{inc}}_{\hat{\mathbf{s}}}&=1/Z^{+}&\exp(i\mathbf{k}^{+}_{\downarrow}\cdot\mathbf{x})\,\hat{\mathbf{s}}\\\\[2.84526pt]
\mathbf{H}^{r}_{\hat{\mathbf{s}}}&=r_{p}/Z^{+}&\exp(i\mathbf{k}^{+}_{\uparrow}\cdot\mathbf{x})\,\hat{\mathbf{s}}\\\\[2.84526pt]
\mathbf{H}^{t}_{\hat{\mathbf{s}}}&=t_{p}/Z^{+}&\exp(i\mathbf{k}^{-}_{\downarrow}\cdot\mathbf{x})\,\hat{\mathbf{s}}\end{array}\right..$
(9)

Then the fully $\hat{\mathbf{p}}$-polarized electric field can be deduced:

$\left\\{\begin{array}[]{llc}\mathbf{E}^{\text{inc}}_{\hat{\mathbf{p}}}&=-\mathbf{k}^{+}_{\downarrow}\times\mathbf{H}^{i}_{\hat{\mathbf{s}}}/(\omega\varepsilon_{0}\varepsilon_{r}^{+})\\\\[2.84526pt]
\mathbf{E}^{r}_{\hat{\mathbf{p}}}&=-\mathbf{k}^{+}_{\uparrow}\times\mathbf{H}^{r}_{\hat{\mathbf{s}}}/(\omega\varepsilon_{0}\varepsilon_{r}^{+})\\\\[2.84526pt]
\mathbf{E}^{t}_{\hat{\mathbf{p}}}&=-\mathbf{k}^{-}_{\downarrow}\times\mathbf{H}^{t}_{\hat{\mathbf{s}}}/(\omega\varepsilon_{0}\varepsilon_{r}^{-})\end{array}\right..$
(10)

Finally the linearly polarized annex electric field $\mathbf{E}_{a}$ solution
of the annex problem writes:

$\mathbf{E}_{a}=A_{e}(\cos\psi_{i}\mathbf{E}_{a,\hat{\mathbf{p}}}-\sin\psi_{i}\mathbf{E}_{a,\hat{\mathbf{s}}}),$
(11)

where

$\mathbf{E}_{a,\\{\hat{\mathbf{s}},\hat{\mathbf{p}}\\}}=(\mathbf{E}^{\text{inc}}_{\\{\hat{\mathbf{s}},\hat{\mathbf{p}}\\}}+\mathbf{E}^{r}_{\\{\hat{\mathbf{s}},\hat{\mathbf{p}}\\}})\mathbb{1}_{\mathcal{S}^{+}}+\mathbf{E}^{t}_{\\{\hat{\mathbf{s}},\hat{\mathbf{p}}\\}}\mathbb{1}_{\mathcal{S}^{-}}.$
(12)

This last expression completes the tools needed to compute the response of a
periodic structure lying upon a substrate through the scattered field
formulation. It allows to write the solution of the Helmholtz equation as an
outgoing field with a source localized into the design region. This outgoing
field condition is modelled numerically by the PMLs that absorb the field
radiating from the design region. Homogeneous Neumann or Dirichlet conditions
on the lower and upper boundaries of $\Omega$ can be chosen to truncate the
PMLs. The former allow to keep track of the field values and the PMLs endings,
while the latter (chosen here) set the field to zero at the PMLs endings while
reducing the number of FEM unknowns. The corresponding functional space
becomes
$\mathcal{V}\coloneqq\mathcal{H}_{0}(\operatorname{\mathbf{curl}}_{\gamma},\Omega,e^{i\alpha
d})$ in order to point out these extra boundary conditions. The scattered
field formulation of direct problem reads:

######  Proposition 1.

The diffracted field
$\mathbf{E}_{*}^{d}\coloneqq\mathbf{E}^{\text{tot}}-\mathbf{E}_{a}$ is the
solution of the following PDE: Find $\mathbf{E}_{*}^{d}\in\bm{\mathcal{V}}$
such that $\forall\mathbf{r}\in\Omega$ with $\omega\in\mathbb{R}^{+}$ fixed,
$\operatorname{\mathbf{curl}}_{\gamma}(\underline{\underline{\mu_{r}}}^{-1}\operatorname{\mathbf{curl}}_{\gamma}\mathbf{E}_{*}^{d})-k_{0}^{2}\underline{\underline{\varepsilon_{r}}}\mathbf{E}_{*}^{d}=k_{0}^{2}(\underline{\underline{\varepsilon_{r}}}-\underline{\underline{\varepsilon_{r,a}}})\mathbf{E}_{a}.$
(13)

Since by construction
$\underline{\underline{\varepsilon_{r}}}-\underline{\underline{\varepsilon_{r,a}}}=0$
out of $\Omega_{d}$, the source (the so-called Right-Hand Side, RHS) is
localized in the design region, which indeed guaranties the outgoing nature of
the chosen diffracted field.

###  3.3 Weak formulation

The decomposition into a tangential and longitudinal part also applies to the
characteristics of the materials. The permittivity and permeability tensors
are decoupled in the following manner:

$\underline{\underline{\varepsilon_{r}}}=\underbrace{\begin{pmatrix}\tilde{\underline{\underline{\varepsilon_{r}}}}&0\\\
0&1\end{pmatrix}}_{\text{\normalsize$\underline{\underline{\varepsilon_{r,t}}}$}}\underbrace{\begin{pmatrix}\underline{\underline{I_{2}}}&0\\\
0&\varepsilon_{r,zz}\end{pmatrix}}_{\text{\normalsize$\underline{\underline{\varepsilon_{r,\ell}}}$}}\quad\text{and}\quad\underline{\underline{\mu_{r}}}=\underbrace{\begin{pmatrix}\tilde{\underline{\underline{\mu_{r}}}}&0\\\
0&1\end{pmatrix}}_{\text{\normalsize$\underline{\underline{\mu_{r,t}}}$}}\underbrace{\begin{pmatrix}\underline{\underline{I_{2}}}&0\\\
0&\mu_{r,zz}\end{pmatrix}}_{\text{\normalsize$\underline{\underline{\mu_{r,\ell}}}$}}$
(14)

so that
$\underline{\underline{\varepsilon}}\mathbf{E}=\underline{\underline{\varepsilon_{t}}}\mathbf{E}_{t}+\underline{\underline{\varepsilon_{\ell}}}E_{\ell}\hat{\mathbf{z}}$.
Gathering all the definitions above and using and using the identity
$\operatorname{\mathbf{curl}}_{\gamma}\mathbf{E}=\operatorname{\mathbf{curl}}_{t}\mathbf{E}_{t}+\left(\operatorname{\mathbf{grad}}_{t}E_{\ell}-i\gamma\mathbf{E}_{t}\right)\times\hat{\mathbf{z}}$,
the following weak formulation in the conical case is obtained:

######  Corollary 1.

The weak formulation of the direct problem in the conical case is: Find
$\mathbf{E}^{d}_{*}\in\bm{\mathcal{V}}$ such that for all
$\mathbf{E}^{\prime}\in\bm{\mathcal{V}}$,
$\displaystyle\int_{\Omega}\Biggl{[}\mu_{r,zz}^{-1}$
$\displaystyle\operatorname{\mathbf{curl}}_{t}\mathbf{E}_{*,t}^{d}\cdot\overline{\operatorname{\mathbf{curl}}_{t}\mathbf{E}_{t}^{\prime}}+\left(\underline{\underline{\mu_{r,t}}}^{-1}(\operatorname{\mathbf{grad}}_{t}E_{*,\ell}^{d}\times\hat{\mathbf{z}})\right)\cdot\overline{\operatorname{\mathbf{grad}}_{t}E_{\ell}^{\prime}\times\hat{\mathbf{z}}}$
(15)
$\displaystyle+i\gamma\left(\underline{\underline{\mu_{r,t}}}^{-1}(\operatorname{\mathbf{grad}}_{t}E_{*,\ell}^{d}\times\hat{\mathbf{z}})\right)\cdot\overline{\mathbf{E}_{t}^{\prime}\times\hat{\mathbf{z}}}-i\gamma\left(\underline{\underline{\mu_{r,t}}}^{-1}(\mathbf{E}_{*,t}^{d}\times\hat{\mathbf{z}})\right)\cdot\overline{\operatorname{\mathbf{grad}}_{t}E_{\ell}^{\prime}\times\hat{\mathbf{z}}}$
$\displaystyle+\gamma^{2}\left(\underline{\underline{\mu_{r,t}}}^{-1}(\mathbf{E}_{*,t}^{d}\times\hat{\mathbf{z}})\right)\cdot\overline{\mathbf{E}_{t}^{\prime}\times\hat{\mathbf{z}}}$
$\displaystyle-
k_{0}^{2}\left(\underline{\underline{\varepsilon_{r,t}}}\mathbf{E}_{*,t}^{d}\cdot\overline{\mathbf{E}_{t}^{\prime}}+\varepsilon_{r,zz}E_{*,\ell}^{d}\overline{E_{\ell}^{\prime}}-(\underline{\underline{\varepsilon_{r,t}^{a}}}-\underline{\underline{\varepsilon_{r,t}}})\mathbf{E}_{a,t}\cdot\overline{\mathbf{E}_{t}^{\prime}}-(\varepsilon_{r,zz}^{a}-\varepsilon_{r,zz})E_{a,\ell}\overline{E_{\ell}^{\prime}}\right)\Biggr{]}\,\mathrm{d}\Omega=0.$

The decomposition of the periodic part of the solution
$\mathbf{E}_{*}^{d}e^{-i\alpha x}$ into Fourier series allows to obtain the
complex amplitudes of each diffraction order [38]:

$\left\\{\begin{array}[]{llc}\displaystyle
r_{n}^{u}=\dfrac{1}{d}\int_{0}^{d}e^{-i\alpha_{n}x}\,\mathbf{E}^{d}(x,y_{0})\cdot\hat{\mathbf{u}}\,\mathrm{d}x&\mbox{for
}y_{0}>y_{\Omega_{d}}\\\ \displaystyle
t_{n}^{u}=\dfrac{1}{d}\int_{0}^{d}e^{-i\alpha_{n}x}\,\mathbf{E}^{\text{tot}}(x,y_{0})\cdot\hat{\mathbf{u}}\,\mathrm{d}x&\mbox{for
}y_{0}<0\end{array}\right.,$ (16)

where $u$ designates $x$, $y$ or $z$ and
$\mathbf{E}^{d}=\mathbf{E}_{*}^{d}+\mathbf{E}_{a}^{d}$ is the total diffracted
field. Making use of the Poynting theorem, these amplitudes lead to the
reflection and transmission efficiencies

$\left\\{\begin{array}[]{llc}R_{n}=\dfrac{\beta_{n}^{+}}{A_{e}^{2}\beta^{+}}(\lvert
r_{n}^{x}\rvert^{2}+\lvert r_{n}^{y}\rvert^{2}+\lvert
r_{n}^{z}\rvert^{2})\\\\[8.53581pt]
T_{n}=\dfrac{\beta_{n}^{-}}{A_{e}^{2}\beta^{+}}(\lvert
t_{n}^{x}\rvert^{2}+\lvert t_{n}^{y}\rvert^{2}+\lvert
t_{n}^{z}\rvert^{2})\end{array}\right..$ (17)

It is possible to access the quantity $\mathbf{E}^{d}\cdot\hat{\mathbf{y}}$,
the normal trace of $\mathbf{E}^{d}$ on the boundary of interest, in the
quantity $r_{n}^{y}$ by using a Lagrange multiplier mapping the normal trace
of the transverse field $\mathbf{E}^{d}_{*,t}$ on the PML/superstrate
boundary. Indeed as this particular component $E^{d}_{*,y}$ is discontinuous
across the superstrate/PML interface, it cannot be simply postprocessed from
$\mathbf{E}^{d}\cdot\hat{\mathbf{y}}$. But, considering the adjoint problem in
section 5, it is more appropriate to use another expression [40] of the
diffraction efficiencies that only involves the components of the field
tangential to the superstrate/PML interface:

$\left\\{\begin{array}[]{llc}R_{n}=\dfrac{1}{\beta_{n}^{+}\beta^{+}A_{e}^{2}}\Bigl{[}((\beta_{n}^{+})^{2}+\alpha_{n}^{2})\lvert
r_{n}^{x}\rvert^{2}+((\beta_{n}^{+})^{2}+\gamma^{2})\lvert
r_{n}^{z}\rvert^{2}+2\alpha_{n}\gamma\text{Re}(r_{n}^{x}\overline{r_{n}^{z}})\Bigr{]}\\\\[5.69054pt]
T_{n}=\dfrac{1}{\beta_{n}^{-}\beta^{+}A_{e}^{2}}\Bigl{[}((\beta_{n}^{-})^{2}+\alpha_{n}^{2})\lvert
t_{n}^{x}\rvert^{2}+((\beta_{n}^{-})^{2}+\gamma^{2})\lvert
t_{n}^{z}\rvert^{2}+2\alpha_{n}\gamma\text{Re}(t_{n}^{x}\overline{t_{n}^{z}})\Bigr{]}\end{array}\right..$
(18)

##  4 Optimization problem

###  4.1 Design variables and constraints

In discrete topology optimization, the permittivity in $\Omega_{d}$ is
considered constant in each element of the mesh and can take any value given
by the SIMP interpolation law defined in Eq. (1) so as to minimize a target
function. This process can be left completely free or guided by additional
constraints. In this paper, two classical constraints are added:

  * •

a connectedness filter in order to avoid unrealistic isolated elements ;

  * •

a binarization filter in order to get only one material or the other in the
optimal pattern.

Let us consider a density (or design) variable $\bm{\rho}$ constant per
element of the mesh. Let thus assume that if the mesh is made of $N$
triangles, $\bm{\rho}$ is a vector of size $N$ and that $0\leq\rho_{i}\leq 1$
for all $i$. Since this density field defines the relative permittivity of the
design space, changing $\bm{\rho}$ modifies the structure and thus the
scattered field $\mathbf{E}_{*}^{d}$ and the diffraction efficiency $R_{n}$.
From now on, the notation $R_{n}$ designates a function
$R_{n}:[0,1]^{N}\ni\bm{\rho}\mapsto R_{n}(\bm{\rho})\in[0,1]$.

The _connectedness constraint_ is imposed using a connectedness map
$\bm{\rho}_{f}$, which is a simple sliding averaging as detailed in Appendix
A. The _binarization_ is applied over this first filter and obtained using the
usual function [25] called $\hat{\bm{\rho}}$ for all $i$:

$\hat{\bm{\rho}}:[0,1]^{N}\ni\bm{\rho}_{f}\mapsto\frac{\tanh(\beta_{f}\nu)+\tanh[\beta_{f}(\bm{\rho}_{f}-\nu)]}{\tanh(\beta_{f}\nu)+\tanh[\beta_{f}(1-\nu)]}\in[0,1]^{N}$
(19)

with $\nu\in[0,1]$ and where $\beta_{f}$ is increased during the optimization
process. A standard configuration is $\nu=1/2$ and $\beta_{f,m}=2^{m}$ with
$m=1,\dots,7$ increasing during the optimization process.

To summarize, the constraints on the design variables are gathered in
functions that coerce the design into having a specific behaviour. It can be
written through a composition of maps:

$\hat{\bm{\rho}}_{f}:[0,1]^{N}\ni\bm{\rho}\mapsto(\hat{\bm{\rho}}\circ\bm{\rho}_{f})(\bm{\rho})\in[0,1]^{N}.$
(20)

Now we are in position to introduce the optimization problem.

###  4.2 Optimization of diffraction efficiencies

The objective of the study is to provide a blazed metasurface which is the
most efficient as possible on one particular diffraction order in reflection.
Note that extending what follows to the case of transmission is
straightforward. Therefore we want to maximize the energy related quantity
$R_{n}$ for a specific value of $n$ on the constrained density field
$\hat{\bm{\rho}}_{f}$. Let then define the composed function
$\mathcal{R}_{n}:[0,1]^{N}\ni\bm{\rho}\mapsto(R_{n}\circ\hat{\bm{\rho}}\circ\bm{\rho}_{f})(\bm{\rho})\in[0,1]$
that returns the diffraction efficiency. We choose to minimize
$\mathcal{F}_{n}\coloneqq 1-\mathcal{R}_{n}$ since $\mathcal{R}_{n}$ has
values between 0 and 1. The DoFs are the densities $\rho_{i}$ into the design
region with $i\in\llbracket 1,N\rrbracket$. The optimization problem writes:

$\displaystyle\underset{\bm{\rho}}{\min}$
$\displaystyle\mathcal{F}_{n}(\bm{\rho})$ (21) such that
$\displaystyle\left\\{\begin{array}[]{llc}\mathcal{L}(\bm{\rho},\mathbf{E}^{\prime})=0&\mbox{$\forall\mathbf{E}^{\prime}\in\bm{\mathcal{V}}_{c}$}\\\
0\leq\rho_{i}\leq 1&\mbox{$\forall i\in\llbracket
1,N\rrbracket$}\end{array}\right.,$

where $\mathcal{L}$ designates the weak formulation in Eq. (15), that is,
written without the tangential/longitudinal decomposition for compactness:

$\begin{array}[]{cccl}\mathcal{L}:&\mathbb{R}^{N}\times\bm{\mathcal{V}}&\longrightarrow&\mathbb{R}\\\
&(\bm{\rho},\mathbf{E}^{\prime})&\longmapsto&\displaystyle\int_{\Omega}\bigg{[}\underline{\underline{\mu_{r}}}^{-1}\operatorname{\mathbf{curl}}_{\gamma}\mathbf{E}_{*}^{d}(\bm{\rho})\cdot\operatorname{\mathbf{curl}}_{\gamma}\overline{\mathbf{E}^{\prime}}-k_{0}^{2}\Big{(}\underline{\underline{\varepsilon_{r}}}(\bm{\rho})\mathbf{E}_{*}^{d}(\bm{\rho})+(\underline{\underline{\varepsilon_{r}}}(\bm{\rho})-\underline{\underline{\varepsilon_{r,a}}})\mathbf{E}_{a}\Big{)}\cdot\overline{\mathbf{E}^{\prime}}\bigg{]}\mathrm{d}\Omega\end{array}.$
(22)

Note that $\mathbf{E}_{*}^{d}$ is the solution of Eq. (13) for a given design
variable $\bm{\rho}_{*}$ called equilibrium point in the framework of
optimization. Given the number of unknowns, this kind of optimization problem
is intractable without the Jacobian of the target $\mathcal{F}_{n}$ with
respect to the design variable $\bm{\rho}$.

###  4.3 Jacobian of the target

As mentioned above, the target function is actually a composition of three
maps, namely:

$\mathcal{F}_{n}(\bm{\rho})=(F_{n}\circ\hat{\bm{\rho}}\circ\bm{\rho}_{f})(\bm{\rho}),$
(23)

where $F_{n}=1-R_{n}$. The chain rule is applied to get the expression of the
derivatives with respect to all $\rho_{i}$:

$\frac{\partial\mathcal{F}_{n}}{\partial\bm{\rho}}(\bm{\rho})=\frac{\partial\bm{\rho}_{f}}{\partial\bm{\rho}}(\bm{\rho})\frac{\partial\hat{\bm{\rho}}}{\partial\bm{\rho}_{f}}(\bm{\rho}_{f})\frac{\partial
F_{n}}{\partial\hat{\bm{\rho}}_{f}}(\hat{\bm{\rho}}_{f}).$ (24)

The notation (bold or not) shows whether the Jacobian is a vector of size $N$
(differentiation of a scalar with respect to a vector) or an $N\times N$
matrix (differentiation of a vector with respect to a vector). The first two
factors are known since the definitions of the constraints on the design
variable are differentiable analytic functions (see Appendix A). Moreover
starting from Eq. (18) and recalling that $\partial\lvert
r_{n}^{u}\rvert^{2}/\partial\hat{\rho}_{i}=2\mathrm{Re}\\{\overline{r_{n}^{u}}\,\partial
r_{n}^{u}/\partial\hat{\rho}_{i}\\}$ [17], we obtain:

$\frac{\partial F_{n}}{\partial\hat{\rho}_{f,i}}=-\frac{\partial
R_{n}}{\partial\hat{\rho}_{f,i}}=-\frac{2}{\beta_{n}^{+}\beta^{+}A_{e}^{2}}\text{Re}\Biggl{[}((\beta_{n}^{+})^{2}+\alpha_{n}^{2})\frac{\partial
r_{n}^{x}}{\partial\hat{\rho}_{f,i}}\overline{r_{n}^{x}}+((\beta_{n}^{+})^{2}+\gamma^{2})\frac{\partial
r_{n}^{z}}{\partial\hat{\rho}_{f,i}}\overline{r_{n}^{z}}+\alpha_{n}\gamma\Bigl{(}\frac{\partial
r_{n}^{x}}{\partial\hat{\rho}_{f,i}}\overline{r_{n}^{z}}+r_{n}^{x}\frac{\partial\overline{r_{n}^{z}}}{\partial\hat{\rho}_{f,i}}\Bigr{)}\Biggr{]}.$
(25)

However, the complex amplitudes $r_{n}^{u}$ depend explicitly on the vector
diffracted field $\mathbf{E}_{*}^{d}$ that is only computed numerically using
the FEM. Consequently the analytic derivatives are not available. An intuitive
but extremely costly solution would consist in computing the derivative
numerically using finite differences for each design variable $\rho_{i}$. This
would require solving the direct problem $N$ times per iteration of the
optimization process, which is clearly prohibitive. This issue led to the
development of the so-called adjoint method.

##  5 Adjoint method

The last hurdle to compute numerically the Jacobian in Eq. (25) is to know the
value of $\partial r_{n}^{u}/\partial\hat{\rho}_{f,i}$ around the current
equilibrium point $\bm{\rho}_{*}$. It is provided by the following
proposition:

######  Proposition 2.

The derivatives of $r_{n}^{u}$ with respect to the $\hat{\rho}_{f,i}$ around
the equilibrium point $\hat{\bm{\rho}}_{f,*}$ are given by $\frac{\partial
r_{n}^{u}}{\partial\hat{\rho}_{f,i}}(\hat{\bm{\rho}}_{f,*})=\int_{\mathcal{T}_{i}}k_{0}^{2}(\varepsilon_{r,\text{diel}}-\varepsilon_{r}^{+})\mathbf{E}^{\text{tot}}(\bm{\rho}_{*})\cdot\bm{\lambda}_{*}^{u}\,\mathrm{d}\Omega,$
(26) where $\mathcal{T}_{i}$ is the triangle of the mesh element with the
density $\hat{\rho}_{f,i}$ and $\bm{\lambda}_{*}^{u}$ is the unique solution
of the _adjoint problem_ , noting
$\bm{\mathcal{V}}^{\text{adj}}=\bm{\mathcal{H}}_{0}(\operatorname{\mathbf{curl}}_{-\gamma},\Omega,e^{-i\alpha
d})$: Find $\bm{\lambda}_{*}^{u}\in\bm{\mathcal{V}}^{\text{adj}}$ such that
for all $\bm{\lambda}^{\prime}\in\bm{\mathcal{V}}^{\text{adj}}$,
$\int_{\Omega}\Biggl{[}\underline{\underline{\mu_{r}}}^{-1}\operatorname{\mathbf{curl}}_{\gamma}\bm{\lambda}_{*}^{u}\cdot\operatorname{\mathbf{curl}}_{\gamma}\overline{\bm{\lambda}^{\prime}}-k_{0}^{2}\underline{\underline{\varepsilon_{r}}}\bm{\lambda}_{*}^{u}\cdot\overline{\bm{\lambda}^{\prime}}\Biggr{]}\,\mathrm{d}\Omega=\frac{1}{d}\int_{0}^{d}e^{-i\alpha_{n}x}\hat{\mathbf{u}}\cdot\overline{\bm{\lambda}^{\prime}}\,\mathrm{d}x.$
(27)

The proof of this proposition is detailed in Appendix B. Note that if the
grating is 3D, then the adjoint problem (27) remains the same, but the source
becomes an integral over a surface $\Gamma^{+}$ in the superstrate. The
following corollary (which proof is also in Appendix B) provides the weak
formulation in order to solve Eq. (27) with the FEM in conical mounting.

######  Corollary 2.

The weak formulation of the adjoint problem in the conical case is: Find
$\bm{\lambda}_{*}^{u}\in\bm{\mathcal{V}}^{\text{adj}}$ such that for all
$\bm{\lambda}^{\prime}\in\bm{\mathcal{V}}^{\text{adj}}$,
$\displaystyle\int_{\Omega}\Biggl{[}\mu_{r,zz}^{-1}$
$\displaystyle\operatorname{\mathbf{curl}}_{t}\bm{\lambda}_{*,t}^{u}\cdot\overline{\operatorname{\mathbf{curl}}_{t}\bm{\lambda}_{t}^{\prime}}+\left(\underline{\underline{\mu_{r,t}}}^{-1}(\operatorname{\mathbf{grad}}_{t}\lambda_{*,\ell}^{u}\times\hat{\mathbf{z}})\right)\cdot\overline{\operatorname{\mathbf{grad}}_{t}\lambda_{\ell}^{\prime}\times\hat{\mathbf{z}}}$
(28)
$\displaystyle\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\boldsymbol{-}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}i\gamma\left(\underline{\underline{\mu_{r,t}}}^{-1}(\operatorname{\mathbf{grad}}_{t}\lambda_{*,\ell}^{u}\times\hat{\mathbf{z}})\right)\cdot\overline{\bm{\lambda}_{t}^{\prime}\times\hat{\mathbf{z}}}\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\boldsymbol{+}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}i\gamma\left(\underline{\underline{\mu_{r,t}}}^{-1}(\bm{\lambda}_{*,t}^{u}\times\hat{\mathbf{z}})\right)\cdot\overline{\operatorname{\mathbf{grad}}_{t}\lambda_{\ell}^{\prime}\times\hat{\mathbf{z}}}$
$\displaystyle+\gamma^{2}\left(\underline{\underline{\mu_{r,t}}}^{-1}(\bm{\lambda}_{*,t}^{u}\times\hat{\mathbf{z}})\right)\cdot\overline{\bm{\lambda}_{t}^{\prime}\times\hat{\mathbf{z}}}\Biggr{]}\,\mathrm{d}\Omega=\underbrace{\frac{1}{d}\int_{0}^{d}e^{-i\alpha_{n}x}\Bigl{(}\overline{\bm{\lambda}_{t}^{\prime}}+\overline{\lambda_{\ell}^{\prime}}\hat{\mathbf{z}}\Bigr{)}\cdot\hat{\mathbf{u}}\,\mathrm{d}x}_{\text{\emph{RHS}}}.$

Note that the only differences with Corollary 1 are the signs before the
$i\gamma$ terms. With these two results, the direct problem and the adjoint
problems on the $x$ and $z$ axes are solved to find the derivatives of
$r_{n}^{u}$, and thus of $\mathcal{F}_{n}$. The numerical resolution for $u=x$
is of the same cost as the direct problem. Then the same Finite Element
factorisation matrix is used for $u=z$ and the resolution on the periodic
pattern only necessitates to generate the new RHS. It allows to get all the
derivatives of the target with a bit more than two Finite Element resolutions,
showing the power of this method.

##  6 Discretization and numerical aspects

The following discrete spaces are used in the simulations:

  * •

The scattered vector field $\mathbf{E}^{d}_{*}\in\bm{\mathcal{V}}$ is by
construction split into its possibly discontinuous transverse components and
longitudinal continuous one:
$\mathbf{E}^{d}_{*}=\mathbf{E}_{*,t}^{d}+E_{*,\ell}^{d}\,\hat{\mathbf{z}}$.
The transverse vector field $\mathbf{E}_{*,t}^{d}$ is discretized using
hierarchical Webb elements [41, 42, 43], with 3 FEM unknowns per edge and 2
per face, that is 11 FEM unknowns per triangle. The longitudinal scalar field
$E_{*,\ell}^{d}$ is discretized using Lagrange elements of the second order
(denoted $P^{2}$), that is 6 FEM unknowns per triangle. Bloch boundary
conditions are applied to both discrete spaces: the unknowns defined on the
right boundary of $\Omega$ ($x=+d/2$) are the same as those defined on the
left boundary ($x=-d/2$) up to a phase shift $e^{+i\alpha d}$. Dirichlet
boundary conditions are applied to the top and bottom boundaries of
$\text{PML}^{+}$ and $\text{PML}^{-}$ respectively.

  * •

The discrete version of the functional space $\bm{\mathcal{V}}^{adj}$ on which
the adjoint vector field $\bm{\lambda}_{*}^{u}$ is defined follows the exact
same steps but one: the phase shift for the Bloch boundary conditions is
$e^{-i\alpha d}$.

  * •

The densities and the Jacobian are constant scalars per mesh triangles. For
the Jacobian, following the method described in Ref. [16], integrals over each
mesh element defined in Eq. (26) are in fact performed by solving a weak
projection of the integrand on the discontinuous constant per element space
($P^{0}$ elements). This projection corresponds to a trivial weak formulation
defined on the design space, which is light and fast to run.

The direct and adjoint problems are solved in parallel with the GetDP software
[30] on a mesh generated by Gmsh [29]. The optimization is led with the GCMMA
[31] using the NLopt package [44]. The accuracy of the numerical scheme is
detailed in Appendix C222The verification of the Jacobian is actually made
with an even more general function: a multi-wavelength target, defined in the
subsection 7.1.2 of this paper..

##  7 Optimized blazed metasurfaces

###  7.1 Patterning a silica slab above a silver substrate

####  7.1.1 Mono-wavelength optimization

Here the study is focused on an application of the method on a concrete
example, using silica (relative permittivity described in Ref. [45]) over a
silver substrate (relative permittivity described in Ref. [46]). The period
$d$ is set up to $3300\text{\,}\mathrm{n}\mathrm{m}$ and for now a single
incident plane wave is considered, with wavelength $\lambda=$
$700\text{\,}\mathrm{n}\mathrm{m}$ and the incident angles
$\theta_{i}=5\degree$, $\varphi_{i}=-66\degree$ and $\psi_{i}=90\degree$.
These parameters can be found in experiments such as in Ref. [47].

Figure 3: Optimization of the pattern on the $-1$ diffraction order for a
period of $3300\text{\,}\mathrm{n}\mathrm{m}$, a target wavelength of
$700\text{\,}\mathrm{n}\mathrm{m}$, the angles $\theta_{i}=5\degree$,
$\varphi_{i}=-66\degree$, $\psi_{i}=90\degree$, an oscillation tolerance of
the target function of $10^{-5}$ and 15726 optimization DoFs in the design
region. Four steps of constraints are displayed: (a) Initial configuration
chosen for this optimization: a graded-indexed blazed "silica to air" layer of
dielectric with a thickness of $650\text{\,}\mathrm{n}\mathrm{m}$ and a linear
spatial variation of the density from 1 (left) to zero (right) ; (b) no
constraint, $R_{-1}=98.3\%$ reached in 18 iterations ; (c) only connectedness
(filter radius of $200\text{\,}\mathrm{n}\mathrm{m}$), $R_{-1}=98\%$ reached
in 100 iterations ; (d) only binarization, $R_{-1}=99\%$ reached in 503
iterations, and (e) with binarization and connectedness, $R_{-1}=97.9\%$
reached in 492 iterations.

The initial configuration of the optimization process remains to be discussed.
Indeed, given the large number of DoFs, the function
$\mathcal{F}_{n}(\bm{\rho})$ exhibits many local minima, and the initial
configuration has a large impact on the resulting local minimum found at the
end of the iterative process. We choose to start from a non-realistic blazed
grating designed to reproduce the same phase shift as the one induced by a
blazed sawtooth silver grating with a blaze angle of $5\degree$ and a period
$d=$ $3300\text{\,}\mathrm{n}\mathrm{m}$. This equivalent layer has a
thickness of $650\text{\,}\mathrm{n}\mathrm{m}$ with a linear graded
permittivity distribution as shown on Fig. 3a and detailed in Appendix D. Such
a configuration already gives an efficiency in the blazed order of 81% with
room for improvement, as shows its spectral response in grey color on Fig. 4a.

The effect of each constraint is illustrated by running four optimizations at
a single wavelength with the same initial configuration. The patterns obtained
on one period as well as the number of iterations and the efficiency reached
at the targeted wavelength 700 nm are shown in Fig. 3. First, neither the
connectedness nor the binarization filter are used (Fig. 3b). Then the
connectedness filter solely (Fig. 3c) or the binarization filter solely (Fig.
3d) is applied. Finally both constraints are taken into account (Fig. 3e).

In this figure, the level of blues designates the density of silica (between
zero and 100%, zero is white) on every triangle of the Finite Element mesh in
one period of the design space. Without binarization, it leads to an
unrealistic blurry (_i.e._ graded-indexed) shape. The sole connectedness
constraint slightly enlarges small details. The effect of the sole
binarization is clear in Fig. 3d since a binary design is obtained. The
combination of both constraints leads to a readable shape. The optimum found
is close to 100% (99%, given that 1% of the total light is absorbed at
$700\text{\,}\mathrm{n}\mathrm{m}$ by Joule effect), therefore an important
and encouraging remark is that neither the binarization nor the connectedness
highly affect the maximal efficiency. Actually, another even higher efficiency
maximum is found with the sole binarization. Indeed, the function to optimize
and its Jacobian change for each constraint. It explains why there is no
intuitive continuity between these results.

The spectral response on the $-1^{\text{st}}$ order as well as the evolution
of the efficiency at the targeted wavelength of
$700\text{\,}\mathrm{n}\mathrm{m}$ during the optimization process with
binarization and connectedness are shown in Fig. 4.

Figure 4: Optimization process leading to the pattern of Fig. 3e. (a) Spectral
response of the $-1^{\text{st}}$ diffraction efficiency in reflection on the
[400,1500] nm spectral range under the target incidence angles. The maximum
efficiency is reached exactly at the targeted wavelength
$700\text{\,}\mathrm{n}\mathrm{m}$, highlighted by the green dashed line. The
green arrow illustrates the improvement enabled by the optimization algorithm
on this very wavelength, starting from the initial configuration which
spectral response is represented in grey ; (b) Evolution of the
$-1^{\text{st}}$ diffraction efficiency at $700\text{\,}\mathrm{n}\mathrm{m}$
during the optimization process, zoom on the [95,98]% efficiency interval.
Each drop corresponds to an increment of the binarization process $m$ which
are imposed at the red dashed lines.

A maximum of efficiency is clearly reached (blue curve in Fig. 4a) at the
wavelength $700\text{\,}\mathrm{n}\mathrm{m}$, which shows that the optimized
result corresponds to a resonant mechanism at the targeted wavelength. The
Fig. 4b shows that the optimization is stable, and that the binarization
process does perturb the optimization for a few iteration only (see the red
dashed lines, showing each increment of the integer binarization parameter $m$
introduced in Eq. (19)). Moreover, the threshold where the binarization
prevents the efficiency to go higher is visible after $m=5$ (iteration 400),
when the density $\rho$ can only take values really close to 0 or 1.

The fact that the efficiency nearly reaches 100% at this wavelength can be
highlighted by displaying the corresponding diffracted field (projected on the
$(Oxy)$ plane), as shown on Fig. 5a. On the latter, three periods of the
pattern and 3$\mu$m of the air above are visible. The target incidence angles
are taken as well as the target wavelength. In particular, with
$\theta_{i}=5\degree$, the angle of reflection on the $-1^{\text{st}}$
diffraction order is $\theta_{-1}=-5\degree$. The fact that the diffracted
field is really close to a plane wave with a 5$\degree$ deflection is a direct
illustration of the 98% efficiency of the pattern in this order. To complete
this study, the field inside the design region is also displayed on Fig. 5b,
revealing a resonance at the bottom left of the design space.

Figure 5: (a) Diffracted field with the optimized pattern of Fig. 3e on three
periods, for the target wavelength
$\lambda=$700\text{\,}\mathrm{n}\mathrm{m}$$ with the target incidence angles.
It is close to a plane wave reflected with the angle $\theta_{-1}=-5\degree$.
(b) Diffracted field in the groove region on the same three periods with the
same specifications. With the resonance on the left of each period, the scale
is changed.

With this mono-wavelength optimization process, a resonant topology was
drilled into the design space. This resonant process has a quite low quality
factor given the width of the resonance, but it is not broadband given the
spectral range targeted. The reflection averaged over the spectral range
[400,1500] nm (defined as
$\int_{400}^{1500}R_{-1}(\lambda)\,\mathrm{d}\lambda/(1500-400)$) is only 51%,
whereas the sawtooth grating already reaches 52%.

####  7.1.2 Multi-wavelength optimization

In order to broaden the spectral interval of the blaze effect, a multi-
wavelength objective is now considered on $N_{\lambda}$ wavelengths in
$\Lambda=\\{\lambda_{1},\dots,\lambda_{N_{\lambda}}\\}$, a set of targeted
wavelengths chosen within the spectral range of interest:

$F_{n}(\bm{\rho})=1-\frac{1}{N_{\lambda}}\sum_{\lambda\in\Lambda}R_{n}(\bm{\rho},\lambda).$
(29)

The density of targeted wavelengths is higher in spectral intervals exhibiting
resonances (see green dots in abscissa of Fig. 6b).

Figure 6: Results for the multi-wavelength optimization on the
$-$1${}^{\text{st}}$ diffraction order for a
$650\text{\,}\mathrm{n}\mathrm{m}$-high pattern on 24 target wavelengths
(green dots on Fig. b) between 400 and $1500\text{\,}\mathrm{n}\mathrm{m}$,
for the angles $\theta_{i}=5\degree$, $\varphi_{i}=-66\degree$,
$\psi_{i}=90\degree$, an solver tolerance of the target of $10^{-5}$, with
connectedness (filter radius of $200\text{\,}\mathrm{n}\mathrm{m}$) and
binarization. (a) Optimized pattern. (b) Spectral response in reflection on
the [400,1500] nm spectral range under the target incidence angles. A
comparison is displayed with the classical blazed grating (in red) and the
initial configuration (in grey) with the same angles of incidence. (c)
Spectral response in reflection with 11 evenly spaced polarization angles
$\psi_{i}\in[0,180]\degree$.

Since each Jacobian of $R_{n}(\bm{\rho},\lambda),\,\,\lambda\in\Lambda$ can be
found following the steps described in the previous mono-wavelength case, the
multi-wavelength Jacobian is trivial to compute. Note that one has to solve
$N_{\lambda}$ direct problems along with $2N_{\lambda}$ adjoint problems,
which increases the computational burden. This is why the code has been
parallelized to run all the wavelengths at once. As at this stage of the
resolution, the problems are totally independent from a wavelength to another,
there is no communication requirement for this parallelization.

The multi-wavelength optimization process has a striking effect on the
bandwidth of the blaze effect as shown on Fig. 6b. For this optimization, the
target is composed of 24 different target wavelengths, separated by
$25\text{\,}\mathrm{n}\mathrm{m}$ on the interval $[400,600]$ nm, by
$50\text{\,}\mathrm{n}\mathrm{m}$ on the interval $[600,1200]$ nm and
$100\text{\,}\mathrm{n}\mathrm{m}$ on the interval $[1200,1500]$ nm. More
precisely, the diffraction efficiency averaged on the bandwidth is
significantly increased, reaching 66%, which is an absolute increasing of 14%
in comparison with the sawtooth grating (52%) on the same wavelength interval.

In fact, the blaze response (blue curve in Fig. 6b) obtained is equal or
higher than that of the sawtooth grating (red curve) over _the whole_ interval
of interest. This higher performance result demonstrates the relevance of the
multi-wavelength approach.

Finally, we stress that although the multi-wavelength optimization was carried
out for a particular incident polarization angle, the dependency of the
response with respect to the polarization angle $\psi$ is moderate, as shown
in Fig. 6c for multiple values of the polarization angle $\psi$ ranging in
$[0,180]\degree$. The maximal discrepancy between all the polarization angles
reaches 17% around $1050\text{\,}\mathrm{n}\mathrm{m}$, while on the
[400,1500] nm interval its averaged value is 8.4%.

###  7.2 Larger design space

The longest wavelengths are not diffracted as efficiently as in the visible
range, which is due to the thickness of the design layer $h_{\ell}$. When
$\lambda>h_{\ell}$, the efficiencies keep dropping down because the material
becomes sub-wavelength in the vertical direction, which is not enough to
provide a sufficient phase shift. Therefore, the thickness of the design
region is increased.

Two ways to adapt the initial configuration are chosen and lead to different
optimal patterns, as displayed on Fig. 7. The first possibility is to widen
the linearly decreasing permittivity on the full design space (Fig. 7a), the
second one is to keep the same initial pattern as on Fig. 3a, while the rest
of the design space above is let to $\bm{\rho}=\mathbf{0}$ (Fig. 7d). Note
that the first choice redshifts the blazing effect central wavelength of the
initial configuration while the second does not.

Figure 7: Results for the multi-wavelength optimization on the
$-$1${}^{\text{st}}$ diffraction order for a
$1300\text{\,}\mathrm{n}\mathrm{m}$-high pattern on 24 target wavelengths
(green dots) between 400 and $1500\text{\,}\mathrm{n}\mathrm{m}$, for the
angles $\theta_{i}=5\degree$, $\varphi_{i}=-66\degree$, $\psi_{i}=90\degree$
with a tolerance of the target of $10^{-5}$, with connectedness (filter radius
of $200\text{\,}\mathrm{n}\mathrm{m}$) and binarization, using two different
initial configurations. (a) Initial configuration with the linearly decreasing
permittivity widened to the full design space. (b) The corresponding optimized
pattern ; (c) Spectral response in reflection on the [400,1500] nm spectral
range under the target incidence angles for this optimal grating (blue curve),
as compared to the initial configuration (grey curve) and the sawtooth grating
(red curve). (d) and (e) Same analysis with the initial configuration of Fig.
3a with $650\text{\,}\mathrm{n}\mathrm{m}$ of the design space above let free.
(f) Spectral response in reflection on the [400,1500] nm spectral range under
the target incidence angles for this optimal grating (blue curve), as compared
to the sawtooth grating (red curve). Two drops are spotted in this spectral
response, at $440\text{\,}\mathrm{n}\mathrm{m}$ and
$820\text{\,}\mathrm{n}\mathrm{m}$ (red dots and dashed lines). Another
optimization has been processed, adding these two wavelengths (violet curve).

The first consequence of the initial configuration of Fig. 7a is the redshift
of its maximum efficiency (see the grey curve on Fig. 7c). The challenge of
the optimization now consists in significantly improving the diffraction of
the visible light. Running the same optimization as in the previous section on
this larger design space leads to the final binary pattern shown on Fig. 7b.
It provides a convincing broadband efficiency (blue curve on Fig. 7c) since
the reflection averaged over the spectral range [400,1500] nm is now of 77%,
outperforming the pattern shown in Fig. 6a by 11% in absolute value, and thus
the classical sawtooth grating by 25%. Therefore the relative difference
between the two average reflectivity is 48%, which points out an outstanding
improvement.

For the initial configuration shown in Fig. 7d, leading to the pattern of Fig.
7e, the improvement is even more impressive, providing an averaged diffraction
efficiency of 81% over the considered bandwidth (blue curve on Fig. 7f, as
compared to the red curve of the sawtooth grating, 29% more in absolute terms,
56% in relative terms). However, two main drops appear at
$440\text{\,}\mathrm{n}\mathrm{m}$ and $820\text{\,}\mathrm{n}\mathrm{m}$ (red
dots and dashed lines on Fig. 7f). Including these two wavelengths in the
targeted-wavelength set $\Lambda$ allows to remove these drops, as highlighted
by the red arrows and the new spectral response in violet. Other sharp drops
appear on this new spectral response (for example at
$620\text{\,}\mathrm{n}\mathrm{m}$) and the efficiency averaged on the
wavelength range slightly decreases to 80%. However the global response is
better distributed, which is a non negligible quality in spectroscopy
applications. While the response dropped to 35% on the first optimal pattern,
it is now higher than 65% except at the end of the interval, where the two
responses are the same anyway. This kind of consideration is a way to improve
the existing code, by automating the correction of deep drops in the spectral
response.

###  7.3 Patterning the traditional sawtooth profile

The model of a dielectric metasurface deposited on a metallic substrate allows
to avoid the plasmonic resonances due to the sharp features of the pattern
[48], a phenomenon observed on the sawtooth grating. However, at this point,
one wonders what would be the outcome if the latter was considered as an
initial configuration to see if the classical triangular design is actually
the best metallic blazed grating.

Moreover, this analysis completes the optimization process because it has been
shown that a different kind of interpolation is more efficient and stable for
the metallic gratings. More precisely, Christiansen _et al._ have developed a
non-linear interpolation scheme in Ref. [28], since the complex refractive
index $n$ is linearly interpolated instead of the relative permittivity:

$\displaystyle\varepsilon_{r}^{d}(\lambda,\mathbf{r})=n(\lambda,\mathbf{r})^{2}$
(30) with $\displaystyle
n(\lambda,\mathbf{r})=n_{1}(\lambda)+(n_{2}(\lambda)-n_{1}(\lambda))\bm{\rho}(\mathbf{r}),$

$n_{2}$ being the refractive index of silver and $n_{1}$ being that of air.
This new interpolation method induces a slight modification in the Proposition
2, detailed in Appendix B, and leads to the new patterns presented in Fig. 8.

Figure 8: Results for the optimization on the $-$1${}^{\text{st}}$ diffraction
order for a $650\text{\,}\mathrm{n}\mathrm{m}$-high silver pattern for the
angles $\theta_{i}=5\degree$, $\varphi_{i}=-66\degree$, $\psi_{i}=90\degree$,
a tolerance on the target of $10^{-5}$, with connectedness (filter radius of
$50\text{\,}\mathrm{n}\mathrm{m}$) and binarization. The initial configuration
is a silver sawtooth pattern with a blaze angle $\alpha=5\degree$. (a)
Optimized pattern for a mono-wavelength optimization at
$573\text{\,}\mathrm{n}\mathrm{m}$ ; (b) Optimized pattern for a multi-
wavelength optimization at the wavelengths designated by the green dots on the
Fig. c ; (c) Spectral response in reflection on the [400,1500] nm spectral
range under the target incidence angles for the mono-wavelength (blue curve)
and the multi-wavelength (green curve) optimized patterns, as respectively
compared to the sawtooth gratings with the blaze angles $\alpha=5\degree$
(grey curve) and $\alpha=6.5\degree$ (black curve).

Both a mono- and a multi-wavelength optimizations have been made on the
traditional silver sawtooth grating, providing respectively the patterns of
Fig. 8a and 8b and the spectral responses of the blue and green curves on Fig.
8c. With the mono-wavelength optimization, a gain of 4% on the maximal
efficiency is observed as well as a slight shift of the maximal efficiency due
to the truncation of the triangle displayed on Fig. 8a. With the multi-
wavelength optimization (performed on the wavelengths highlighted by the green
dots on Fig. 8c), this shift is even more pronounced. Therefore, even if the
reflection averaged over the spectral range [400,1500] nm is of 64% (against
52% for the sawtooth grating with a blaze angle of 5$\degree$), the higher
frequencies in the visible range are too much cut. The response obtained is
actually slightly widened with respect to that of a sawtooth grating with a
higher blaze angle (here $\alpha=6.5\degree$, as highlighted by the black
curve). The multi-wavelength optimization has only exhibited the performances
of a traditional triangular blazed grating with a maximum more centered in the
[400,1500] nm wavelength range (here at $700\text{\,}\mathrm{n}\mathrm{m}$).

This final study presents two advantages. On the one hand, it shows that the
optimization of metallic patterns is stable and accurate with this non-linear
interpolation method. However on the other hand, it exhibits that the
performances of the triangular profile are actually almost optimal among the
metallic gratings. The plasmonic resonances are consequently not the only
reason why the optimization is performed on dielectrics for the blazed
gratings: the triangular metallic pattern does not have much room for
improvement.

##  8 Conclusion

In this article, we provide all the theoretical steps allowing to compute the
topology optimization with the adjoint method using the FEM under a conical
incidence (2.5D). By numerically solving only two PDEs under variational form,
an accurate Jacobian of the reflection/transmission diffraction efficiencies
with respect to a relative permittivity density is computed in order to make
the (GCMMA) optimization algorithm functional. The applications shown in the
numerical section are concrete examples demonstrating the high-performance of
the algorithm. The blazed diffraction efficiencies can reach 98% for a single-
wavelength optimization. Even more striking, for a multi-wavelength or
broadband optimization, the traditional sawtooth grating blazed on its
$-1^{\text{st}}$ diffraction order is outperformed by topology optimized
gratings that increases the diffraction efficiency averaged on the [400,1500]
nm spectral range by up to 29% in absolute terms and 56% in relative terms.
The algorithm even points out that the usual sawtooth pattern cannot be truly
improved, hence the choice to use the dielectric-on-metal model.

A lot of complementary studies are possible from now on. In the frame of open
science, we provide the source code allowing to retrieve the numerical
examples of this paper. Other materials, geometries, constraints and targets
can be easily taken into account starting from the provided Gmsh/GetDP
template model.

The preliminary results presented here are an introduction to a more practical
work on the manufacturing of optimized blazed metasurfaces which is under way.
The manufacturability of the device is a very novel field of research since
the gratings already produced are either for larger wavelengths [22], or
really small devices [27, 49]. To our knowledge two nano-blazed grating on a
convex surface have been manufactured and tested [50, 51]. Being able to also
fabricate the optimized grating would enable to oversee all the steps of the
device development: theoretical modelling, manufacturing, characterization,
comparison with the model and eventually integration in an actual instrument.

## Acknowledgements

This work has been partially supported by CNES and Thales Aliena Space with a
PhD grant.

## References

  * [1] M. P. Bendsøe and O. Sigmund.  Topology optimization.  Springer, 2004.
  * [2] S. Dou and J. S. Jensen.  "Optimization of nonlinear structural resonance using the harmonic balance method".  Journal of Sound and Vibration, 334:239–245, 2015.
  * [3] J. S. Jensen.  "A simple method for coupled acoustic-mechanical analysis with application to gradient-based topology optimization".  Structural and Multidisciplinary Optimization, 59(5):1567–1580, 2019.
  * [4] A. Gupta, R. Chowdhury, A. Chakrabarti, and T. Rabczuk.  "A 55-line code for large-scale parallel topology optimization in 2D and 3D".  2020\.
  * [5] D. G. Nielsen, F. T. Agerkvist, and J. S. Jensen.  "Optimization of realistic loudspeaker models with respect to basic response characteristics".  In Proceedings of 23rd International Congress on Acoustics, pages 6219–6226. German Acoustical Society (DEGA), 2019.
  * [6] T. Borrvall and J. Petersson.  "Topology optimization of fluids in Stokes flow".  International Journal for Numerical Methods in Fluids, 41:77–107, 2003.
  * [7] Y. Deng, Z. Liu, P. Zhang, et al.  "Topology optimization of unsteady incompressible Navier–Stokes flows".  Journal of Computanional Physics, 230:6688–6708, 2011.
  * [8] E. Kuci, F. Henrotte, P. Duysinx, and C. Geuzaine.  "Combination of topology optimization and Lie derivative-based shape optimization for electro-mechanical design".  Structural and Multidisciplinary Optimization, 59(5):1723–1731, 2019.
  * [9] X. Qian and O. Sigmund.  "Topological design of electromechanical actuators with robustness toward over- and under-etching".  Computer Methods in Applied Mechanics and Engineering, 253(10):237–251, 2013.
  * [10] J. S. Jensen and O. Sigmund.  "Topology optimization of photonic crystal structures: a high bandwidth low-loss T-junction waveguide".  Journal of the Optical Society of America B, 22(6):1191–1198, 2005.
  * [11] K. Friis and O. Sigmund.  "Robust topology design of periodic grating surfaces".  Journal of the Optical Society of America B, 29(10):2935–2943, 2012.
  * [12] Y. Deng and J. G. Korvink.  "Self-consistent adjoint analysis for topology optimization of electromagnetic waves".  Journal of Computational Physics, 361:353–376, 2018.
  * [13] J. L. Pita Ruiz, A. S. A. Amad, L. H. Gabrielli, and A. A. Novotny.  "Optimization of the electromagnetic scattering problem based on the topological derivative method".  Optics Express, 27(23):33586–33605, 2019.
  * [14] G. Allaire, C. Dapogny, and F. Jouve.  Shape and topology optimization.  Springer, 2004.
  * [15] G. Allaire, C. Dapogny, and P. Frey.  "Shape optimization with a level set based mesh evolution method".  Computer Methods in Applied Mechanics and Engineering, 282:22–53, 2014.
  * [16] E. Kuci, F. Henrotte, P. Duysinx, and C. Geuzaine.  "Design sensitivity analysis for shape optimization based on the Lie derivative".  Computer Methods in Applied Mechanics and Engineering, 317:702–722, 2017.
  * [17] G. Bao, P. Li, and J. Lv.  "Numerical solution of an inverse diffraction grating problem from phaseless data".  Journal of the Optical Society of America A, 30(3):293–299, 2013.
  * [18] Z. Xiong, W. He, Q. Wang, et al.  "Design and optimization method of a convex blazed grating in the Offner imaging spectrometer".  Applied Optics, 60(2):383–391, 2021.
  * [19] O. Sandfuchs, M. Kraus, and R. Brunner.  "Structured metal double-blazed dispersion grating for broadband spectral efficiency achromatization".  Journal of the Optical Society of America A, 37(8):1369–1380, 2020.
  * [20] B. Borguet, V. Moreau, E. Renotte, et al.  "The dual-blazed diffraction grating of the CHIME hyperspectral instrument: design, modelling, and breadboarding".  In Proceedings of SPIE 12777, International Conference on Space Optics – ICSO 2022, page 127772H, 2023.
  * [21] W. T. Chen, J. S. Park, J. Marchioni, et al.  "Dispersion-engineered metasurfaces reaching broadband 90% relative diffraction efficiency".  Nature Communications, 14(2544), 2023.
  * [22] X. Li, M. Memarian, K. Dhwaj, and T. Itoh.  "Blazed metasurface grating: the planar equivalent of a sawtooth grating".  IEEE MTT-S International Microwave Symposium (IMS), pages 1–3, 2016.
  * [23] J. Jiang, E. Lupoiu, E. W. Wang, et al.  "MetaNet: a new paradigm for data sharing in photonics research".  Optics Express, 28(9):13670–13681, 2020.
  * [24] T. Phan, D. Sell, E. W. Wang, et al.  "High-efficiency, large-area, topology-optimized metasurfaces".  Light: Science & Applications, 8(48), 2019.
  * [25] B. Vial and Y. Hao.  "Open-source computational photonics with auto-differentiable topology optimization".  Mathematics, 10(20):3912, 2022.
  * [26] R. E. Christiansen and O. Sigmund.  "Inverse design in photonics by topology optimization: tutorial".  Journal of the Optical Society of America B, 38(2):496–509, 2021.
  * [27] Z. Li, R. Pestourie, Z. Lin, et al.  "Empowering metasurfaces with inverse design: principles and applications".  ACS Photonics, 9(7):2178–2192, 2022.
  * [28] R.E. Christiansen, J. Vester-Petersen, S. P. Madsen, and O. Sigmund.  "A non-linear material interpolation for design of metallic nano-particles using topology optimization".  Computer Methods in Applied Mechanics and Engineering, 343:23–29, 2021.
  * [29] C. Geuzaine and J. F. Remacle.  "Gmsh: a three-dimensional finite element mesh generator with built-in pre- and post-processing facilities".  International Journal for Numerical Methods in Engineering, 79(11):1309–1331, 2009.
  * [30] P. Dular, C. Geuzaine, F. Henrotte, and W. Legros.  "A general environment for the treatment of discrete problems and its application to the finite element method".  IEEE Transactions on Magnetics, 34(5):3395–3398, 1998.
  * [31] K. Svanberg.  "A class of globally convergent optimization methods based on conservative convex separable approximations".  Society for Industrial and Applied Mathematics, 12(2):559–573, 2002.
  * [32] F. Zamkotsian, I. Zhurminsky, P. Lanzoni, et al.  "Blazed gratings on convex substrates for high throughput spectrographs for Earth and Universe Observation".  In Proceedings of SPIE 11852, International Conference on Space Optics – ICSO 2020, page 118520N, 2021.
  * [33] M. P. Bendsøe and O. Sigmund.  "Material interpolation schemes in topology optimization".  Archive of Applied Mechanics, 69(10):635–654, 1999.
  * [34] Z. Houta, F. Messine, and T. Huguet.  "Topology optimization for magnetic circuits with adjoint method in 3D".  Structural and Multidisciplinary Optimization, 2023.
  * [35] J. P. Berenger.  "A perfectly matched layer for the absorption of electromagnetic waves".  Journal of Computational Physics, 114:185–200, 1994.
  * [36] G. Floquet.  "Sur les équations différentielles linéaires à coefficients périodiques".  Annales Scientifiques de l’Ecole Normale Supérieure, 12:47–88, 1883.
  * [37] G. Demésy, A. Nicolet, and F. Zolla.  "Open Source Models for the Parametric Study of Diffraction Gratings in 2D/2.5D/3D with ONELAB".  2022\.
  * [38] G. Demésy, F. Zolla, A. Nicolet, and B. Vial.  Chapter 5: Finite element method.  In E. Popov, editor, Gratings: Theory and Numeric Applications. Presses Universitaires de Provence, 2nd revisited edition, 2014.
  * [39] F. Zolla, G. Renversez, A. Nicolet, et al.  Chapter 3: Electromagnetism – prerequisties.  In Foundations of Photonics Crystal Fibres. Imperial College Press, 2nd edition, 2012.
  * [40] F. Zolla.  Contribution à létude de la diffraction et de l’adsorption des ondes électromagnétiques : Structures bipériodiques minces, structures cylindriques par la méthode des sources fictives, section 3.4.  PhD thesis, Aix-Marseille Université, 1993.
  * [41] C. Geuzaine, B. Meys, P. Dular, and W. Legros.  "Convergence of high order curl-conforming finite elements [for EM field calculations]".  IEEE Transactions on Magnetics, 35(3), 1999.
  * [42] J. P. Webb and B. Forgahani.  "Hierarchal scalar and vector tetrahedra".  IEEE Transactions on Magnetics, 29(2):1495–1498, 1993.
  * [43] J. Jin.  Chapter 8: Vector finite elements.  In The Finite Element Method in Electromagnetics. John Wiley & Sons Inc., 3rd edition, 2014.
  * [44] S. G. Johnson.  The NLopt nonlinear-optimization package.  https://github.com/stevengj/nlopt, 2007.
  * [45] I. H. Malitson.  "Interspecimen comparison of the refractive index of fused silica".  Journal of the Optical Society of America, 55:1205–1208, 1965.
  * [46] P. B. Johnson and R. W. Christy.  "Optical constants of the noble metals".  Physical Review B, 6(4370–4379), 1972.
  * [47] F. Zamkotsian, P. Lanzoni, N. Tchoubaklian, et al.  "BATMAN @ TNG: instrument integration and performance".  In Proceedings of SPIE 10702, Ground-based and Airborne Instrumentation for Astronomy VII, page 107025P, 2018.
  * [48] S. A. Maier.  In Plasmonics: Fundamentals and Applications. Springer, 2007.
  * [49] M. Mansouree, A. McClung, S. Samudrala, and A. Arbabi.  "Large-scale parametrized metasurface design using adjoint optimization".  ACS Photonics, 8(2):455–463, 2021.
  * [50] F. Burmeister, T. Flügel-Paul, U. Zeitner, et al.  "Binary blazed reflection grating for UV/VIS/NIR/SWIR spectral range".  In Proceedings of SPIE 11180, International Conference on Space Optics – ICSO 2018, page 111801J, 2018.
  * [51] M. S. L. Lee, J. Cholet, A. Delboulbé, et al.  "Wide Band UV/Vis/NIR blazed-binary reflective gratings for spectro-imagers: two lithographic technologies investigation".  Journal of the European Optical Society - Rapid Publications, 19(1):7, 2023.
  * [52] R. Potthast.  "Fréchet differentiability of boundary integral operators in inverse acoustic scattering".  Inverse Problems, 10(2):431–447, 1994.

## Appendix

Technical proofs and validations that are not essential for the understanding
of the paper are detailed here. The same notations and definitions as in the
article are used.

First, additional information is given on the connectedness and binarization
filters on the density field and their Jacobian in appendix A. Then the proof
of proposition 2 and its corollary are given in appendix B, which can be
useful to the curious reader who needs to modify the cost function. A
validation of the computation of the Jacobian through this adjoint method is
shown afterwards in appendix C, by direct comparison to a brute element-by-
element finite difference approach. Finally, the choice of the initial
configuration of the optimization is explained in appendix D.

##  Appendix A Filters applied and their Jacobian

The diffusion filter used here is a linear function that homogenizes the
values around every mesh element with a given radius $r_{f}$. Let
$\mathcal{T}_{i}$ be the triangle with the density $\rho_{i}$. Let then
$\text{dist}_{j,r_{f}}$ be a function such that

$\text{dist}_{j,r_{f}}(\rho_{i})=\left\\{\begin{array}[]{llc}\text{dist}(\mathcal{T}_{i},\mathcal{T}_{j})&\mbox{if
the latter is smaller than $r_{f}$}\\\ 0&\mbox{else}\end{array}\right.,$ (31)

where the distance between two triangles is the distance between their
barycenters. This distance function selects the surroundings triangles closer
than the filter radius. Then the diffusion filter function is defined by:

$\rho_{f,i}(\bm{\rho})=\frac{\sum_{j=1}^{N}\text{dist}_{j,r_{f}}(\rho_{i})\rho_{j}}{\sum_{j=1}^{N}\text{dist}_{j,r_{f}}(\rho_{i})}.$
(32)

In essence, a sliding averaging is performed on the densities. One can
actually recognize a linear matrix function. Let $\mathbb{D}$ be the matrix of
the different values of $\text{dist}_{j,r_{f}}$:
$\mathbb{D}_{i,j}=\text{dist}_{j,r_{f}}(\rho_{i})$ and let $\mathfrak{1}$ be
the column vector with ones. Then we can rewrite the filter function

$\rho_{f,i}(\bm{\rho})=\frac{\mathbb{D}\bm{\rho}}{\mathbb{D}\mathfrak{1}}$
(33)

where the division between vectors in the last expression designates the
division element per element. This expression directly leads to the Jacobian
of the filter:

$\frac{\partial\bm{\rho}_{f}}{\partial\bm{\rho}}(\bm{\rho})=\frac{\mathbb{D}}{\mathbb{D}\mathfrak{1}}$
(34)

where the division between a matrix and a vector in the last expression is the
division between each line of the matrix by this vector.

Concerning the binarization filter, a mere analytic differentiation leads to

$\frac{\partial\hat{\rho}_{i}}{\partial\rho_{f,j}}(\bm{\rho}_{f})=\frac{\beta_{f}(1-\tanh^{2}[\beta_{f}(\rho_{f,i}-\nu)])}{\tanh(\beta_{f}\nu)+\tanh[\beta_{f}(1-\nu)]}\delta_{ij},$
(35)

where $\delta_{ij}$ is the Kronecker delta. The matrix multiplication
(34)$\times$(35) (in this order) provides the filtering factor in the
derivative of $\mathcal{F}_{n}$. The last factor for the total derivative is
given by the proposition 2, which proof is detailed in the next section of
this technical support.

##  Appendix B Proof of the proposition 2 and its corollary

###  B.1 Proof of the proposition

This proof is adapted from the general result presented in Ref. [8]. For the
sake of simplicity, we denote only by $\bm{\rho}$ the filtered density field
(instead of $\hat{\bm{\rho}}_{f}$). First it is important to notice that
$r_{n}^{u}$ depends explicitly on $\mathbf{E}_{*}^{d}$ that depends implicitly
on $\bm{\rho}$. The derivative of the composite function writes then for all
$\bm{\rho}\in\mathbb{R}^{N}$:

$\frac{\partial
r_{n}^{u}}{\partial\rho_{i}}(\mathbf{E}_{*}^{d}(\bm{\rho}))=\mathchoice{{\frac{\partial
r_{n}^{u}}{\partial\rho_{i}}\,\smash{\vrule
height=10.07558pt,depth=3.7778pt}}_{\,\frac{\partial\mathbf{E}}{\partial\rho_{i}}=0}}{{\frac{\partial
r_{n}^{u}}{\partial\rho_{i}}\,\smash{\vrule
height=7.06667pt,depth=3.62411pt}}_{\,\frac{\partial\mathbf{E}}{\partial\rho_{i}}=0}}{{\frac{\partial
r_{n}^{u}}{\partial\rho_{i}}\,\smash{\vrule
height=5.38318pt,depth=3.62411pt}}_{\,\frac{\partial\mathbf{E}}{\partial\rho_{i}}=0}}{{\frac{\partial
r_{n}^{u}}{\partial\rho_{i}}\,\smash{\vrule
height=5.38318pt,depth=3.62411pt}}_{\,\frac{\partial\mathbf{E}}{\partial\rho_{i}}=0}}(\mathbf{E}_{*}^{d}(\bm{\rho}))+\Bigl{\\{}\mathrm{D}_{\mathbf{E}}r_{n}^{u}\Bigr{\\}}\biggl{(}\frac{\partial\mathbf{E}_{*}^{d}}{\partial\rho_{i}}(\bm{\rho})\biggr{)},$
(36)

where $\mathrm{D}_{\mathbf{E}}$ designates the Fréchet derivative which is the
linear operator such that [52]:

$r_{n}^{u}(\mathbf{E}_{*}^{d}(\bm{\rho})+\delta\mathbf{E}(\bm{\rho}))=r_{n}^{u}(\mathbf{E}_{*}^{d}(\bm{\rho}))+\Bigl{\\{}\mathrm{D}_{\mathbf{E}}r_{n}^{u}\Bigr{\\}}(\delta\mathbf{E}(\bm{\rho}))+o(\delta\mathbf{E})$
(37)

with

$\lim_{\lVert\delta\mathbf{E}\rVert\to
0}\frac{1}{\lVert\delta\mathbf{E}\rVert}o(\delta\mathbf{E})=0.$

As $r_{n}^{u}$ is linear with respect to $\mathbf{E}_{*}^{d}(\bm{\rho})$, the
first term in (36) vanishes and moreover

$\displaystyle
r_{n}^{u}(\mathbf{E}_{*}^{d}(\bm{\rho})+\delta\mathbf{E}(\bm{\rho}))$
$\displaystyle=\frac{1}{d}\int_{0}^{d}e^{-i\alpha_{n}x}(\mathbf{E}_{*}^{d}(\bm{\rho})+\delta\mathbf{E}(\bm{\rho})+\mathbf{E}_{a}^{d})\cdot\hat{\mathbf{u}}\,\mathrm{d}x$
$\displaystyle=\frac{1}{d}\int_{0}^{d}e^{-i\alpha_{n}x}(\mathbf{E}_{*}^{d}(\bm{\rho})+\mathbf{E}_{a}^{d})\cdot\hat{\mathbf{u}}\,\mathrm{d}x+\frac{1}{d}\int_{0}^{d}e^{-i\alpha_{n}x}\delta\mathbf{E}(\bm{\rho})\cdot\hat{\mathbf{u}}\,\mathrm{d}x.$

Thus we deduce that for all $\bm{\rho}\in\mathbb{R}^{N}$:

$\Bigl{\\{}\mathrm{D}_{\mathbf{E}}r_{n}^{u}\Bigr{\\}}\biggl{(}\frac{\partial\mathbf{E}_{*}^{d}}{\partial\rho_{i}}(\bm{\rho})\biggr{)}=\frac{1}{d}\int_{0}^{d}e^{-i\alpha_{n}x}\frac{\partial\mathbf{E}_{*}^{d}}{\partial\rho_{i}}(\bm{\rho})\cdot\hat{\mathbf{u}}\,\mathrm{d}x=\frac{\partial
r_{n}^{u}}{\partial\rho_{i}}(\mathbf{E}_{*}^{d}(\bm{\rho})).$ (38)

This derivative cannot be calculated numerically since it would rely on a
computation of every derivative of
$\partial\mathbf{E}_{*}^{d}/\partial\bm{\rho}$. The same consideration as for
$r_{n}^{u}$ applies here: one could use numerical finite differences. However
an alternative method is to consider the augmented Lagrangian of $r_{n}^{u}$
regarding the Helmholtz PDE.

For $u=\\{x,z\\}$ let $r_{n}^{u}$ be a performance function and $r_{n,a}^{u}$
its so-called augmented Lagrangian:

$r_{n,a}^{u}(\mathbf{E}_{*}^{d}(\bm{\rho}),\mathbf{E}^{\prime})=r_{n}^{u}(\mathbf{E}_{*}^{d}(\bm{\rho}))-\mathcal{L}(\bm{\rho},\mathbf{E}^{\prime}).$
(39)

Then the same reansoning as equations (36) to (38) is made for $\mathcal{L}$
and the derivatives of this object are for all
$(\bm{\rho},\bm{\lambda})\in\mathbb{R}^{N}\times\bm{\mathcal{V}}$,

$\displaystyle\frac{\partial
r_{n,a}^{u}}{\partial\rho_{i}}(\mathbf{E}_{*}^{d}(\bm{\rho}),\mathbf{E}^{\prime})=$
$\displaystyle\int_{0}^{d}\frac{1}{d}e^{-i\alpha_{n}x}\frac{\partial\mathbf{E}_{*}^{d}}{\partial\rho_{i}}(\bm{\rho})\cdot\hat{\mathbf{u}}\,\mathrm{d}x-\int_{\Omega}\Biggl{[}\underline{\underline{\mu_{r}}}^{-1}\operatorname{\mathbf{curl}}\frac{\partial\mathbf{E}_{*}^{d}}{\partial\rho_{i}}(\bm{\rho})\cdot\operatorname{\mathbf{curl}}\overline{\mathbf{E}^{\prime}}$
(40) $\displaystyle-
k_{0}^{2}\biggl{(}\underline{\underline{\varepsilon_{r}}}\frac{\partial\mathbf{E}_{*}^{d}}{\partial\rho_{i}}(\bm{\rho})+\frac{\partial\underline{\underline{\varepsilon_{r}}}}{\partial\rho_{i}}(\bm{\rho})\underbrace{(\mathbf{E}_{*}^{d}(\bm{\rho})+\mathbf{E}_{a})}_{\mathbf{E}^{\text{tot}}}\biggr{)}\cdot\overline{\mathbf{E}^{\prime}}\Biggr{]}\,\mathrm{d}\Omega.$

We consider now that, for a given $i\in\llbracket 1,N\rrbracket$,
$\bm{\lambda}\coloneqq\overline{\mathbf{E}^{\prime}}$ is a function of an
adjoint functional space $\bm{\mathcal{V}}^{\text{adj}}$ that is to be
characterized and
$\bm{\lambda}^{\prime}\coloneqq\partial\overline{\mathbf{E}_{*}^{d}}/\partial\rho_{i}$
is a test function in this same adjoint space. The crucial remark to
understand the adjoint method is that, if there exists a function
$\bm{\lambda}_{*}^{u}=\overline{\mathbf{E}^{\prime}}$ such that only the last
term in (40) (with the underbrace) remains, then it would be possible to
compute the derivatives of the augmented Lagrangian. Furthermore, for the
equilibrium design variable $\bm{\rho}_{*}$ such that
$\mathcal{L}(\bm{\rho_{*}},\mathbf{E}^{\prime})=0$ for all
$\mathbf{E}^{\prime}\in\bm{\mathcal{V}}$, we can state [16] that
$\frac{\partial\mathcal{L}}{\partial\rho_{i}}(\bm{\rho}_{*})=0$, so that

$\frac{\partial
r_{n,a}^{u}}{\partial\rho_{i}}(\mathbf{E}_{*}^{d}(\bm{\rho}_{*}),\bm{\lambda}_{*}^{u})=\frac{\partial
r_{n}^{u}}{\partial\rho_{i}}(\mathbf{E}_{*}^{d}(\bm{\rho}_{*})).$ (41)

In other terms, finding such an adjoint variable $\bm{\lambda}_{*}^{u}$ would
solve the issue of calculating the derivatives of $r_{n}^{u}$. Let then
consider the so-called adjoint problem:

Find $\bm{\lambda}_{*}^{u}\in\bm{\mathcal{V}}^{\text{adj}}$ such that for all
$\bm{\lambda}^{\prime}\in\bm{\mathcal{V}}^{\text{adj}}$,

$\int_{\Omega}\Biggl{[}\underline{\underline{\mu_{r}}}^{-1}\operatorname{\mathbf{curl}}\bm{\lambda}_{*}^{u}\cdot\operatorname{\mathbf{curl}}\overline{\bm{\lambda}^{\prime}}-k_{0}^{2}\underline{\underline{\varepsilon_{r}}}\bm{\lambda}_{*}^{u}\cdot\overline{\bm{\lambda}^{\prime}}\Biggr{]}\,\mathrm{d}\Omega-\frac{1}{d}\int_{0}^{d}e^{-i\alpha_{n}x}\hat{\mathbf{u}}\cdot\overline{\bm{\lambda}^{\prime}}\,\mathrm{d}x=0.$
(42)

This problem has a unique solution indeed, since it is the weak formulation of
a Maxwell’s equation with a surface current $\mathbf{j}\coloneqq
e^{-i\alpha_{n}x}/d\,\,\hat{\mathbf{u}}$ as a source.

Eventually, combining Eq. (41) with the injection of the solution of Eq. (42)
into Eq. (40) provides all the derivatives of $r_{n}^{u}$ with respect to the
$\rho_{i}$ around the equilibrium point $\bm{\rho}_{*}$:

$\frac{\partial
r_{n}^{u}}{\partial\rho_{i}}(\bm{\rho}_{*})=\int_{\Omega}k_{0}^{2}\frac{\partial\underline{\underline{\varepsilon_{r}}}}{\partial\rho_{i}}(\bm{\rho})\mathbf{E}^{\text{tot}}(\bm{\rho}_{*})\cdot\bm{\lambda}_{*}^{u}\,\mathrm{d}\Omega.$
(43)

The last step is to express the derivatives of
$\underline{\underline{\varepsilon}}_{r}$ with respect to $\rho_{i}$. The
relative permittivity actually varies only in the design region. Moreover, a
variation around a single density $\rho_{i}$ only induces a variation in the
triangle $\mathcal{T}_{i}$. Lastly, since the interpolation methods seen in
the article provide a direct link between $\varepsilon_{r}^{d}$ and
$\bm{\rho}$, the expression of
$\partial\underline{\underline{\varepsilon}}_{r}/\partial\rho_{i}$ is:

  * •

$\varepsilon_{r,\text{diel}}-\varepsilon_{r}^{+}$ in a dielectric, using the
linear SIMP interpolation method [33], which concludes the proof of the
proposition 2 ;

  * •

in a metal, using the non-linear interpolation method [28],

$\frac{\partial\underline{\underline{\varepsilon_{r}}}}{\partial\rho_{i}}(\bm{\rho})=2(n_{2}-n_{1})\cdot\Big{(}n_{1}+\bm{\rho}(n_{2}-n_{1})\Big{)}.$

###  B.2 Proof of the corollary

As it can be seen between the corollaries 1 and 2, the only difference between
both weak formulations is the sign before $\gamma$ (written in red in the
corollary 2). This aspect is detailed here. More precisely, the central
element of this proof is the characterization of
$\bm{\mathcal{V}}^{\text{adj}}$. Let $\mathbf{E}$ be a function of
$\bm{\mathcal{V}}$. In the conical case, it means that it can be decomposed
using its quasi-periodicity with a factor $\alpha$ on the one hand and the
exponential dependency in $z$ on the other hand. Then for all
$\mathbf{x}\coloneqq(x,y,z)\in\mathbb{R}^{3}$,

$\mathbf{E}(\mathbf{x})=\mathbf{E}^{\\#}(x,y)e^{i\alpha x}e^{i\gamma z}$ (44)

with $\mathbf{E}^{\\#}$ a wave periodic along $x$. Also by definition
$\overline{\mathbf{E}}\in\bm{\mathcal{V}}^{\text{adj}}$ and we have that

$\overline{\mathbf{E}}(\mathbf{x})=\overline{\mathbf{E}^{\\#}}(x,y)e^{-i\alpha
x}e^{-i\gamma z}.$ (45)

The reciprocal of this statement is immediate by using the exact same
arguments. Furthermore, using now $\operatorname{\mathbf{curl}}_{-\gamma}$ for
the elements of $\bm{\mathcal{V}}^{\text{adj}}$ leads to the sign change in
the weak formulation of the corollary 2.

##  Appendix C Numerical validation of the adjoint method

The Finite Element Method for the direct problem or the convergence of the
methods are common elements that have to be checked in order to certify that
the algorithm provides a trustable result. An important point here is to check
the computation of the Jacobian with the adjoint method. To do so, we compare
it to the Jacobian found with finite differences.

Let consider a grid with a random density distribution (see Fig. 9a). The
derivative of the target is computed both ways, as illustrated on Fig. 9:
first with the adjoint method and then with the finite differences. The direct
problem is solved $N$ times in order to have the variation of
$\mathcal{F}_{n}$ when adding a small step $h$ to $\rho_{i}$ for every $i$
(for every triangle). Note that the derivatives of the filters described in
section 2 of this document are included in this study. For this test, the step
of the finite differences is equal to $h=10^{-2}$ and the binarization factor
is $\beta_{f}=8$. Also, a conical incidence is chosen, with
$\theta_{i}=5\degree$, $\varphi_{i}=-66\degree$ and $\psi_{i}=90\degree$, just
like in the optimization case of the article. Finally, this validation is made
for a multi-wavelength target function with wavelengths $\\{400,700,900\\}$nm.
The relative differences between the adjoint method and finite differences
displayed on Fig. 9c are defined on each mesh element $\mathcal{T}_{i}$ by:

$d_{a,f}=\left|\frac{\partial_{i}\mathcal{F}_{\text{finiteDiff}}-\partial_{i}\mathcal{F}_{\text{adjoint}}}{\text{mean}|\partial_{i}\mathcal{F}_{\text{finiteDiff}}|}\right|.$
(46)

It is smaller than $4\%$ in all mesh elements, even with this coarse mesh.
This error decreases when the mesh is refined.

Figure 9: Test of the adjoint method on a small grid for an incident field
with angles $\theta_{i}=5\degree$, $\varphi_{i}=-66\degree$,
$\psi_{i}=90\degree$ and with a multi-wavelength target function with
wavelengths $\\{400,700,900\\}$ nm on the specular order. (a) Initial random
distribution of densities. (b) Derivatives obtained for each mesh element,
both lead to the same distribution. (c) Relative differences between both
methods (adjoint and finite differences) for each mesh element.

##  Appendix D Choice of the initial configuration

For any optimization process, an initial configuration is required. Its choice
is not obvious at first sight and it is important since its choice conditions
the minimum of the cost function found by the optimization algorithm. In this
study, we decided to consider as initial configuration blazed grating made of
a graded indexed dielectric lying on a silver substrate. The linear decay of
the graded permittivity (see Fig 10b) in the design region tries to mimic the
phase shift obtained with a geometrical linear ramp of the silver triangle
blazed grating (see Fig 10a). The only unknown in this problem is the height
$y_{\Omega_{d}}$ of the graded layer. It is found by considering optical paths
difference between $x=$ $0\text{\,}\mathrm{n}\mathrm{m}$ and $x=$
$3300\text{\,}\mathrm{n}\mathrm{m}$ as shown in red color in Fig 10b,
similarly to the (approximate) rule of thumb used for traditional blazed
gratings.

Figure 10: Imitation of the phase shift $\Delta\phi$ induced by a triangular
grating providing the initial configuration for the optimization, described
with the optical paths of the light in both configurations. (a) Optical path
for the triangular grating. (b) Approximated optical path for the equivalent
density distributed pattern. (c) (black lines) Phase of the diffracted field
on a period of the triangular pattern ($x\in$[0,3300]nm) with a 2D incoming
field with $\lambda=600$nm, $\theta_{i}=5\degree$ in TE (solid lines) and TM
(dashed lines) polarizations. (red lines) Ideal case where a plane wave is
reflected on the $-$1 diffraction order. (d) Same analysis for the equivalent
density distribution pattern.

On Fig. 10a, the height $y_{\Omega_{T}}$ of the silver triangle is
$d\tan\alpha$. Therefore the difference $\Delta d_{o,T}^{n}(\lambda)$ of
optical paths between $x=$ $0\text{\,}\mathrm{n}\mathrm{m}$ and $x=$
$3300\text{\,}\mathrm{n}\mathrm{m}$ for the triangle is

$\Delta
d_{o,T}^{n}(\lambda)=-d\tan\alpha\Big{(}\frac{1}{\cos\theta_{i}}+\frac{1}{\cos\theta_{n}(\lambda)}\Big{)}.$
(47)

For the graded-indexed grating shown in Fig. 10b, the optical rays of interest
are the ray entering the cell at $x=0$ and the ray leaving the cell at $x=d$.
The horizontal size of the region where those rays and their respective
reflection travel is
$s_{h}(\lambda)=d|\sin\theta_{i}-\sin\theta_{n}(\lambda)|$. Let now consider
that $\theta_{i}\ll\pi/2$ and $\theta_{n}(\lambda)\ll\pi/2$. In that case
$s_{h}(\lambda)\ll d$ and it can be considered that the first ray evolves in a
dielectric with $n=n_{2}(\lambda)$ and the second one with $n=n_{1}(\lambda)$.
This is why approximately, the optical distance for the first light beam is
$n_{2}(\lambda)y_{\Omega_{d}}(\lambda)(1/\cos\theta_{i}+1/\cos\theta_{n}(\lambda))$,
and
$n_{1}(\lambda)y_{\Omega_{d}}(\lambda)(1/\cos\theta_{i}+1/\cos\theta_{n}(\lambda))$
for the second one. The optical path difference $\Delta
d_{o,\ell}^{n}(\lambda)$ between the two ends of the design region is thus:

$\Delta
d_{o,\ell}^{n}(\lambda)=y_{\Omega_{d}}(\lambda)(n_{1}(\lambda)-n_{2}(\lambda))\Big{(}\frac{1}{\cos\theta_{i}}+\frac{1}{\cos\theta_{n}(\lambda)}\Big{)}.$
(48)

Setting this optical path difference $\Delta d_{o,\ell}^{n}(\lambda)$ to that
of the traditional triangular grating $\Delta d_{o,T}^{n}(\lambda)$ leads to:

$y_{\Omega_{d}}(\lambda)=\frac{d\tan\alpha}{n_{2}(\lambda)-n_{1}(\lambda)}.$
(49)

The higher $n_{2}(\lambda)$ is, the lower $y_{\Omega_{d}}(\lambda)$ is. For
the multi-wavelength optimization, the range [400,1500]nm is considered and
the lower refractive index of silica is for $\lambda=$
$1500\text{\,}\mathrm{n}\mathrm{m}$. In this case $n_{2}(\lambda)=1.445$ and
by taking $d=$ $3300\text{\,}\mathrm{n}\mathrm{m}$ and $\alpha=5\degree$, we
obtain $y_{\Omega_{d}}(\lambda)\simeq$ $650\text{\,}\mathrm{n}\mathrm{m}$.

An illustration of the phase shift induced by the graded-indexed structure is
shown in Fig. 10d and can be compared to the phase shift induced by the
triangular grating shown in Fig. 10c. The phase of the diffracted field
computed using the FEM is plotted in black lines for an incoming field with a
small angle of incidence ($5\degree\simeq
0.087\,\text{rad}\ll\pi/2\,\text{rad}$). The wavevector lies in the plane of
incidence ($\varphi_{i}=0$) and the computation is made in TE (solid lines)
and TM (dashed lines) polarizations cases. We compare the phase shift obtained
numerically to a perfect plane wave reflected on the $-$1${}^{\text{st}}$
diffraction order with a $y$-intercept chosen the same as the for real field
(red lines). It shows that the approximate design rule guided by physical
optics provides a satisfying phase shift. Actually, the spectral response of
this initial configuration (in grey on Fig. 4a for instance) indicates that
this grating is quite blazed for a broad wavelengths range, the blazed
efficiency being over 50% in the [450,1100] nm range.