## Abstract

The main goal of the method proposed in this paper is the numerical study of various kinds of anisotropic gratings deposited on isotropic substrates, without any constraint upon the diffractive pattern geometry or electromagnetic properties. To that end we propose a new FEM (Finite Element Method) formulation which rigorously deals with each infinite issue inherent to grating problems. As an example, 2D numerical experiments are presented in the cases of the diffraction of a plane wave by an anisotropic aragonite grating on silica substrate (for the two polarization cases and at normal or oblique incidence). We emphasize the interesting property that the diffracted field is non symmetric in a geometrically symmetric configuration.

©2007 Optical Society of America

## 1. Introduction

The field of application of diffraction gratings has been extended over the past decades. They are commonly used in various industrial and scientific domains, such as integrated optics, acoustooptics or spectroscopy [1]. For instance, they proved to be particularly efficient as anti-reflective surfaces [2], beam shapers, waveguide couplers [3] or spectral filters [4]. More than a dozen of numerical methods (coupled-wave approach [5, 6, 7], C method [8, 9], differential method [10], integral method [11], OAGSM method [12], FDTD method [13, 14], FEM [15, 16, 17], MFS [18], …) have already been developed as an answer to this strong need of conceiving and optimizing gratings. Among these abundant numerical approaches, a few of them were adapted to the modeling of anisotropic gratings [6, 7, 9, 10, 12, 16].

The main difficulty in modeling diffraction of a plane wave by a grating arises from the fact that we are dealing with infinite issues while the computational domain has to be bounded. Most of the field calculation methods (FEM [17] and FDTD [14]) use Perfectly Matched Layers -PML[19]- to truncate the free space along one direction and pseudo-periodic conditions along the others, but still remains the problem of the localization of the plane wave sources.

For instance Bao *et al*. raised this source localization problem in the case of an isotropic grating [17]. They pointed out that many other FEM formulations used transparent boundaries to get around this difficulty which is resulting in the truncation of infinite series due to the use of a nonlocal operator. In this context Bao *et al*. confined the sources in the top PML as an alternative solution. Furthermore, Ohkawa *et al*. adapted a FEM to the study of anisotropic gratings [16], but without using any PML which leads to a rather complicated formulation.

In this paper, we present a new approach treating rigorously this issue resulting in a radiation problem with sources localized in the diffractive element itself. We mathematically split the whole problem into two parts. The first one consists in the classical calculation of the *total field* solution of a simple interface. The second one amounts to looking for a *radiated field* with sources confined within the diffractive obstacles and deduced from the first elementary problem. From this viewpoint, the latest *radiated field* can be interpreted as an e*xact perturbation of the total field*.

Traditionally, the diffracted field is the response of the obstacle *i.e*. the region of space where *ε *and *µ* are different from their vacuum values. The difficulty in the present problem is that the *diffractive obstacle* is of finite extent in the basis cell as the substratum is different from the superstratum as it is often the case in optics. In this case, the classical formulation based on the *diffracted field* involves equivalent sources whose domain extents to infinity under the groove region. In our approach, the sources are confined in the groove region and this eases the design of the PML since our *unknown radiated field* satisfies the OutgoingWave Condition (O.W.C.). Moreover, it is straightforward to generalize our formulation to obtain a simple computation for obstacles embedded in a multilayer structure. Finally, its implementation is independent of the geometry of the diffractive object.

Our method can be set up for anisotropic materials such as CaCO_{3} [20], LiNbO_{3} [16], Ni:YIG [21] and even in lossy cases (CoPt, CoPd [22]) without increasing the computational time or ressource. Since the dielectric tensor of magneto-optic (electro-optic) materials can be modified once placed into a constant magnetic (electric) field, it strongly raises the possibility to conceive controlable optical devices.

## 2. Set up of the problem and notations

We denote by **x, y** and **z**, the unit vectors of the axes of an orthogonal co-ordinate system *Oxyz*. We deal only with time-harmonic fields; consequently, the electric and magnetic fields are represented by the complex vector fields **E** and **H**, with a time dependance in exp(-*iωt*).

Besides, in this paper, we assume that the tensor fields of relative permittivity ε͇ and relative permeability µ ͇can be written as follows:

where *ε _{xx},εa,…µ_{zz}* are possibly complex valued functions of the two variables

*x*and

*y*and where $\overline{\epsilon}$

*(resp. µ*

_{a}*̄*) represents the conjugate complex of

_{a}*ε*(resp.

_{a}*µ*).

_{a}*These kinds of materials are said to be z-anisotropic*. It is of importance to note that with such tensor fields, lossy materials can be studied (the lossless materials correspond to tensors with real diagonal terms represented by Hermitian matrices) and that the problem is invariant along the

*z*-axis but the tensor fields can vary continuously (gradient index gratings) or discontinuously (step index gratings). Moreover we define

*k*

_{0}:=

*ω/c*.

The gratings that we are dealing with are made of three regions (See Fig. 1).

• *The superstratum* (*y>H*) which is supposed to be homogeneous, isotropic and lossless and characterized solely by its relative permittivity ε^{+} and its relative permeability *µ*
^{+} and we denote ${k}^{+}\u2254{k}_{0}\sqrt{{\epsilon}^{+}{\mu}^{+}}$

• *The substratum* (*y*<0) which is supposed to be homogeneous and isotropic and therefore characterized by its relative permittivity *ε-* and its relative permeability *µ-* and we denote ${k}^{-}\u2254{k}_{0}\sqrt{{\epsilon}^{-}{\mu}^{-}}$

• *The groove region* (0<*y*<*H*) which can be heterogeneous and *z*-anisotropic and thus characterized by the two tensor fields *ε͇ ^{g}*(

*x,y*) and

*µ͇*(

^{g}*x,y*). It is worth noting that the method presented in this paper does work irrespective of whether the tensor fields are piecewise constant. The groove periodicity along

*x*-axis will be denoted

*d*.

This grating is illuminated by an incident plane wave of wave vector **k**
^{+}=α**x**+*β*+y=*k*
^{+ }(sinθ0**x**-cosθ_{0}y), whose electric field (TM or s-polarization case) (resp. magnetic field (TE or p-polarization case)) is linearly polarized along the *z*-axis:

where **A**
^{0}
_{e} (resp. **A**
^{0}
* _{m}*) is an arbitrary complex number. The magnetic (resp. electric) field derived from

**E**

^{0}

*(resp.*

_{e}**H**

^{0}

*) is denoted*

_{m}**H**

^{0}

*(resp.*

_{e}**E**

^{0}

*) and the electromagnetic field associated with the incident field is therefore denoted (*

_{m}**E**

^{0},

**H**

^{0}) which is equal to (

**E**

^{0}

*,*

_{e}**H**

*0*

*) (resp. (*

_{e}**E**

^{0}

*,*

_{m}**H**

^{0}

*)).*

_{m}The problem of diffraction that we address in this paper is therefore to find Maxwell equation solutions in harmonic regime *i.e*. the unique solution (**E,H**) of:

such that the diffracted field satisfies an *OutgoingWaves Condition* (O.W.C. [24]) and where **E** and **H** are quasi-periodic functions with respect to the *x* co-ordinate.

## 3. Theoretical developments of the method

#### 3.1. Decoupling of fields and z-anisotropy

We assume that δ͇(*x,y*) is a *z*-anisotropic tensor field (*δ _{xz}*=

*δ*=

_{yz}*δ*=

_{zx}*δ*=0). Moreover, the left upper matrix extracted from δ͇ is denoted $\underset{\xaf}{\underset{\xaf}{\widehat{\delta}}}$ , namely:

_{zy}For *z*-anisotropic materials, in a non-conical case, the problem of diffraction can be split into two fundamental cases (TE case and TM case). This property results from the following equality which can be easily derived:

where *u* is a function which does not depend on the *z* variable. Relying on the previous equality, it appears that the problem of diffraction in a non conical mounting amounts to looking for an electric (resp. magnetic) field which is polarized along the *z*-axis ; **E**=*e(x,y)z* (resp. **H**=*h(x,y)z)*. The functions *e* and *h* are therefore solutions of similar differential equations:

with

in the TM case and

in the TE case.

#### 3.2. How to reduce a problem of diffraction to a problem of radiation with localized sources?

In its initial form, the problem of diffraction summed up by Eq. (6) is not well suited to the Finite Element Method. In order to overcome this difficulty, we propose to split the unknown function *u* into a sum of two functions *u*
_{1} and *u*
* ^{d}_{2}*, the first term being known as a closed form and the latter being a solution of a problem of radiation whose sources are localized within the obstacles.

We have assumed that outside the groove region (cf. Fig. 1), the tensor field ξ and the function *χ* are constant and equal respectively to ξ͇- and χ- in the substratum (*y*<0) and equal respectively to ξ͇+ and χ+ in the superstratum (*y>H*). Besides, for the sake of clarity, the superstratum is supposed to be made of an isotropic and lossless material and is therefore solely defined by its relative permittivity *ε*
^{+} and its relative permeability *µ*
^{+}, which leads to:

or

where Id_{2} is the 2×2 identity matrix. With such notations, ξ͇ and *χ* are therefore defined as follows:

It is now apropos to introduce an auxiliary tensor field ξ͇1 and an auxiliary function *χ*1:

these quantities corresponding, of course, to a simple plane interface. Besides, we introduce the constant tensor field ξ͇0 which is equal to *ξ*͇^{+ }everywhere and a constant scalar field *χ*0 which is equal to *χ*+ everywhere. Finally, we denote *u*0 the function which equals the incident field *u*
^{inc} in the superstratum and vanishes elsewhere

We are now in a position to define more precisely the problem of diffraction that we are dealing with. The function *u* is the unique solution of

In order to reduce this problem of diffraction to a radiation problem, an intermediate function is necessary. This function, called *u*
_{1}, is defined as the unique solution of the equation:

The function *u*
_{1} corresponds thus to *an annex problem* associated to a simple interface and can be solved in closed form and *from now on is considered as a known function*. As written above, we need the function *u ^{d}_{2}* which is simply defined as the difference between

*u*and

*u*

_{1}:

The presence of the superscript *d* is, of course, not irrelevant : As the difference of two diffracted fields, the O.W.C. of *u ^{d}_{2}* is guaranteed (which is of prime importance when dealing with PML cf. 3.4). As a result, the Eq. (14) becomes:

where the right hand member is a scalar function which may be interpreted as a *known source term -𝒮1(x,y) and the support of this source is localized only within the groove region*. To prove it, all we have to do is to use Eq. (15):

Now, let us point out that the tensor fields ξ͇ and ξ͇1 are identical outside the groove region and the same holds for *χ* and *χ*1. The support of 𝒮_{1} is thus localized within the groove region as expected. It remains to compute more explicitly the source term 𝒮_{1}. Making use of the linearity of the operator 𝓛 and the equality *u*
_{1}=*u*
^{d}
_{1}+*u*0, the source term can be splitted into two terms

where

and

Now, bearing in mind that u0 is nothing but a plane wave *u*
_{0}=exp(*i*
**k**+·**r**) (with **k**
^{+}=α**x**+β+y), it is sufficient to give ∇* _{u0}*=

*i*

**k**+

*u*for the weak formulation associated with Eq. (17):

_{0}The same holds for the term associated with the diffracted field (*u ^{d}*

_{1}=

*ρ*exp(

*i*

**k**-·

**r**)).

where *ρ* is nothing but the complex reflection coefficient associated with the simple interface :

#### 3.3. Quasi-periodicity and weak formulation

The weak formulation follows the classical lines and is based on the construction of a weighted residual of Eq. (6), which is multiplied by the complex conjugate of a weight function *u*′ and integrated by part to obtain :

The solution *u* of the weak formulation can therefore be defined as the element of the space *L*
^{2}(**curl**,*d,k*) of quasiperiodic functions (i.e. such that *u(x,y)*=*u#(x,y)e ^{ikx}* with

*u#(x,y)*=

*u#(x+d,y)*, a

*d*-periodic function) of

*L*

^{2}(

**curl**) on Ω such that:

As for the boundary term introduced by the integration by part, it can be classically set to zero by imposing Dirichlet conditions on a part of the boundary (the value of *u* is imposed and the weight function *u*′ can be chosen equal to zero on this part of the boundary) or by imposing homogeneous Neumann conditions (ξ͇∇*u*)·**n**=0 on another part of the boundary (and *u* is therefore an unknown to be determined on the boundary). A third possibility are the so-called quasiperiodicity conditions of particular importance in the modelling of gratings. Denote by Γ* _{l}* and Γ

*the lines parallel to the*

_{r}*y*-axis delimiting a cell of the grating respectively from its left and right neighbor cell. Considering that both

*u*and

*u*′ are in

*L*

^{2}(

**curl**,

*d,k*), the boundary term for Γ

*∪Γ*

_{l}*is*

_{r}${\int}_{{\Gamma}_{l}\cup {\Gamma}_{r}}\overline{u\prime}(\underset{\xaf}{\underset{\xaf}{\xi}}\nabla u)\xb7\mathbf{n}dS={\int}_{{\Gamma}_{l}\cup {\Gamma}_{r}}\overline{{u}_{\#}^{\text{'}}}{e}^{-\mathrm{ikx}}(\underset{\xaf}{\underset{\xaf}{\xi}}\nabla \left({u}_{\#}{e}^{+\mathrm{ikx}}\right))\xb7\mathbf{n}dS=$${\int}_{{\Gamma}_{l}\cup {\Gamma}_{r}}\overline{{u}_{\#}^{\text{'}}}\left(\underset{\xaf}{\underset{\xaf}{\xi}}\left(\nabla {u}_{\#}+\mathrm{ik}{u}_{\#}\mathbf{x}\right)\right)\xb7\mathbf{n}dS=0$

because the integrand $\overline{{u}_{\#}^{\text{'}}}\left(\underset{\xaf}{\underset{\xaf}{\xi}}\left(\nabla {u}_{\#}+\mathit{ik}{u}_{\#}\mathbf{x}\right)\right)\xb7\mathbf{n}$ is periodic along *x* and the normal **n** has opposite directions on Γ* _{l}* and Γ

*so that the contributions of these two boundaries have the same absolute value with opposite signs. The contribution of the boundary terms vanishes therefore naturally in the case of quasiperiodicity.*

_{r}The finite element method is based on this weak formulation and both the solution and the weight functions are classically chosen in a discrete space made of linear or quadratic Lagrange elements, i.e. piecewise first or second order two variable polynomial interpolation built on a triangular mesh of the domain Ω (cf. Fig.2a). Dirichlet and Neumann conditions may be used to truncate the PML domain in a region where the field (transformed by the PML) is negligible. The quasiperiodic boundary conditions are imposed by considering the *u* as unknown on Γ* _{l}* (in a way similar to the homogeneous Neumann condition case) while, on Γ

*,*

_{r}*u*is forced equal to the value of the corresponding point on Γ

*(i.e. shifted by a quantity -*

_{l}*d*along

*x*) up to the factor

*e*. The practical implementation in the finite element method is described in details in [25, 26]

^{ikd}#### 3.4. Perfectly Matched Layer for z–anisotropic materials

The main drawback encountered in electromagnetism when tackling theory of gratings through the finite element method is the non-decreasing behaviour of the propagating modes in superstratum and substratum (if they are made of lossless materials): The PML has been introduced by [19] in order to get round this obstacle. The computation of PML designed for *z*-anisotropic gratings is the topic of what follows. Note that the calculation is made for the general case of a *z*-anisotropic substrate, whereas only the case of isotropic substrate will be dealt in the numerical case of the last part. Besides, from now on, the substratum (resp. superstratum) will refer to the truncated space between the groove region and the upper (lower) bound of the lower (upper) PML, whereas before these regions were considered infinite.

The PML may be introduced as a change of co-ordinate corresponding to a *complex stretch* of the co-ordinate corresponding to the direction along which the field must decay [28, 29, 30, 31]. In our study, rectangular PML have been chosen. For this kind of PML, we first introduce new tensor fields $\underset{\xaf}{\underset{\xaf}{{\epsilon}_{s}}}(x,y)$ and $\underset{\xaf}{\underset{\xaf}{{\mu}_{s}}}(x,y)$ defined as follows:

where *J _{s}* is the so-called stretched Jacobian matrix. In our study, this latest matrix is defined as a diagonal matrix field,

*J*=diag (1,

_{s}*s*(

_{y}*y*),1), where

*s*(

_{y}*y*) is a possibly complex valued function of

*y*. For this kind of Jacobian, the computation of $\underset{\xaf}{\underset{\xaf}{{\delta}_{s}}}$ presents no difficulty:

We are now in a position to introduce a new electromagnetic field (called substituted field in the sequel) *F** _{s}*=(

**E**

_{s},

**H**

*s*)

*which is solution of Eq. (3) except that we have replaced ε͇ and µ by $\underset{\xaf}{\underset{\xaf}{{\epsilon}_{s}}}$ and $\underset{\xaf}{\underset{\xaf}{{\mu}_{s}}}$. The main feature of this latest field is the remarkable correspondence with the first field*

^{T}**; whatever the function sy provided that it equals 1 for**

*F**Y*+∗<

*y*<

*Y*-∗ (cf. Fig. 2a), the two fields

**and**

*F*

*F**are identical in the region*

_{s}*Y*+∗<

*y*<

*Y*-∗ [27]. In other words, the PML is completely reflectionless. In addition, for complex valued functions

*s*(

_{y}*∑m*{

*s*} strictly positive in PML), the field

_{y}*does converge exponentially towards zero (as*

**F**_{s}*y*tends to ±∞, cf. Fig. 2c and 2d) although its physical counterpart

*does not. As a consequence,*

**F***is of finite energy and for this substituted field a weak formulation can be easily derived which is essential when dealing with Finite Element Method.*

**F**_{s}Moreover, when dealing with simple stretched functions it is possible to give a simple criterion which ensures the exponential decreasing of the field * F_{s}*. For instance, take the following example :

${s}_{y}\left(y\prime \right)=\{\begin{array}{ll}{\zeta}^{+}& y>{Y}_{*}^{+}\\ {\zeta}^{-}& y<{Y}_{*}^{-}\\ 1& \mathrm{elsewhere}\end{array}$

where ζ^{+} and ζ- are complex numbers. In that case, the complex function *y*(*yc*) is given by:

$y\left({y}_{c}\right)=\{\begin{array}{ll}{Y}_{*}^{+}+{\zeta}^{+}\left({y}_{c}-{Y}_{*}^{+}\right)& \mathrm{for}\phantom{\rule{.2em}{0ex}}{y}_{c}>{Y}_{*}^{+}\\ {y}_{c}& \mathrm{for}\phantom{\rule{.2em}{0ex}}{Y}_{*}^{-}<{y}_{c}<{Y}_{*}^{+}\\ {Y}_{*}^{-}+{\zeta}^{-}\left({y}_{c}-{Y}_{*}^{-}\right)& \mathrm{for}\phantom{\rule{.2em}{0ex}}{y}_{c}<{Y}_{*}^{-}.\end{array}$.

Now, let us consider a propagating plane wave in the substratum *u _{n}*(

*x,y*) :=exp(

*i*(

*αx*-β-+,

*)) can be reexpressed in stretched coordinates in the PML*

_{n}y*as per*:

The behaviour of this latest function along the *y _{c}* direction is governed by the function

*U*:=

^{sc}(yc)*e*-

*iβ-+,nζ-yc*. Letting

*β*′-+,

*n*:=ℜ

*e*{

*β*-+,

*}, β″-+,*

_{n}*n*:=ℑ

*m*{β-+,

*}, ζ′- :=ℜ*

_{n}*e*{ζ-} and ζ″- :=ℑ

*m*{ζ-}, the non-oscillating part of the function

*U*is given by exp ((

^{sc}(yc)*β*′-+,

*ζ″-+β″-+,*

_{n}*ζ′-)*

_{n}*yc*). Keeping in mind that

*β*′-+,

*and/or*

_{n}*β*″-+,

*are positive numbers, the function*

_{n}*U*decreases exponentially towards zero as

^{sc}*y*tends to -∞ (Fig. 2d) provided that ζ- belongs to ℂ+ :={

_{c}*z*∈ℂ,ℜe{

*z*}>0, and ℑ

*m*{

*z*}>0}. In the same way, it can be shown that ζ

^{+}belongs to ℂ

^{+}.

#### 3.5. Synthesis of the method

In order to give a general view of the method, all information is collected here that is necessary to set up the practical Finite Element Model. First of all, the computation domain Ω (cf. Fig. 2a) corresponds to a truncated cell of the grating which is a finite rectangle divided into five horizontal layers. These layers are respectively from top to bottom upper PML, the superstratum, the groove region, the substratum, and the lower PML. The unknown field is the scalar function *u ^{d}_{2}* defined in Eq. (16). Its finite element approximation is based on the second Lagrange elements built on a triangle meshing of the cell (cf. Fig. 2b). A complex algebraic system of linear equations is constructed via the Galerkin weighted residual method,

*i.e*. the set of weight functions

*u*′ is chosen as the set of shape functions of interpolation on the mesh [25].

• In region 1 (upper PML, see Fig. 2a),

with ξ^{͇}+* _{s}* and

*χ+*depending on the equivalent anisotropic properties of the PML given by Eq. (7), Eq. (8), Eq. (27) and Eq. (28).

_{s}• In region 2 (superstratum),

with *ξ*͇^{+} and *χ*
^{+ }depending on the homogeneous isotropic properties of the superstratum given by Eq. (7), Eq. (8), Eq. (9) and Eq. (10).

• In region 3 (groove region),

with *ξ͇ ^{g}* and

*χ*depending on the heterogeneous possibly anisotropic properties given by Eq. (7), Eq. (8), Eq. (11) and 𝒮1 given by Eq. (19), Eq. (22), Eq. (23) and Eq. (24).

^{g}• In region 4 (substratum),

with *ξ͇*- and *χ*- depending on the homogeneous isotropic properties of the substratum given by Eq. (7), Eq. (8), Eq. (9) and Eq. (10).

• In region 5 (lower PML),

with ξ͇-s and *χ*-s depending on the equivalent anisotropic properties of the PML given by Eq. (7), Eq. (8), Eq. (27) and Eq. (28).

## 4. Numerical experiments

#### 4.1. Diffraction efficiencies calculation

The rough result of the FEM calculation is the total complex field solution of Eq. (6) at each point of the bounded domain. We deduce from *u ^{d}* (cf Eq. (14)) the diffraction efficiencies with the following method. The superscripts

^{+}(resp. -) correspond to quantities defined in the superstratum (resp. substratum) as previously.

On the one hand, since *u ^{d}* is quasi-periodic along the

*x*-axis, it can be expanded as a Rayleigh expansion (see for instance [24]):

where

On the other hand, introducing Eq. (35) into Eq. (6) leads to the Rayleigh coefficients :

For a temporal dependance in *e ^{-iωt}*, the O.W.C. imposes

*s*=

_{n}*u*=0. Combining Eq. (36) and Eq. (37) at a fixed

_{n}*y*

_{0}altitude leads to :

We extract these two coefficients by trapezoidal numerical integration along *x* from a cutting of the previously calculated field map at *y*
_{0}. It is well known that the mere trapezoidal integration method is very efficient for smooth and periodic functions (integration on one period) [32]. Now the restriction on a horizontal straight line crossing the whole cell in homogeneous media (substratum and superstratum) is of*C*∞ class. From a numerical point of view, it appears that the interpolated approximation of the unknown function, namely *u ^{d}_{2}* preserves the good behaviour of the numerical computation of these integrals. From this we immediately deduce the reflected and transmitted diffracted efficiencies of propagative orders (

*T*and

_{n}*R*) defined by :

_{n}This calculation is performed at twenty different *y*
_{0} altitudes in the superstratum and the substratum, and the mean value for each propagative transmitted or reflected diffraction order is presented in the following numerical examples.

#### 4.2. Numerical validation of the method

We can refer to [23] in order to test the accuracy of our method. The studied grating is isotropic, but there is a lack of numerical experiments in the literature treating the very same alignments of the incident plane wave and grating co-ordinate system with the principal axes of z-anisotropic materials. We computed the following problem (cf. Fig.3), as described in [17] and [23]. The wavelength of the incident plane wave is set to 1*µm* and is incoming with an angle of *π*/6 with respect to the normal to the grating. In this example, *ε*
^{+}=1 and *ε*-=*ε ^{g}*=-44.9757+2.9524

*i*. We present the

*R*

_{0}efficiency (cf. Table 1) in the two cases of polarization in function of the mesh refinement. So we have a good agrement to the reference values, and the accuracy reached is independent from the polarization case.

#### 4.3. Experiment set up based on existing materials

Let us now consider a trapezoidal (cf. Fig. 4) anisotropic grating made of aragonite (CaCO_{3}) deposited on an isotropic substratum (SiO_{2}, εS_{i}O_{2}=2.25). Along the anisotropic crystal axis, its dielectric tensor can be written as follows [20] :

Now let’s assume that the natural axis of our aragonite grating are rotated by +45° around the grating infinite dimension. The dielectric tensor becomes :

We shall here remind that our method remains strictly the same whatever the diffractive element geometry is. The 2D computational domain is bounded along the y-axis by the PML and along the x since we consider only one pseudo period. We propose to calculate the diffractive efficiencies at *λ*
_{0}=633*nm* in the two cases of polarization TE and TM, and for different incoming incidences (0°, -20° and -40°). Since both µ͇ and ε͇ are Hermitian, the whole incident energy is diffracted and the sum of theses efficiencies ought to be equal to the incident energy, which will stand for validation of our numerical calculation.

Finally, the resulting bounded domain is meshed with a maximum mesh element side size of *λ*
_{0}/10√*ε*. Efficiencies are still post-processed in accordance with the calculation presented section 4.1.

#### 4.4. Numerical results

At normal incidence, the *h* field in the TE case (cf. Fig. 5b) is non symmetric whereas the *e* field in the TM case is (cf. Fig. 5a). This is illustrated by the obvious non-symmetry of *T*
^{TE}-1 and *T*
^{TE}1 (cf. Table 2: 0.322510 versus 0.124722), whereas *T*
^{TM}-1=*T*
^{TM}1=0.20313. Note that the *material asymmetry ratio*
$\frac{0.251}{2.592}$ between the diagonal term and off-diagonal term of ${\underset{\xaf}{\underset{\xaf}{\epsilon}}}_{\mathrm{Ca}{\mathrm{CO}}_{3}}^{+45\xb0}$ is much smaller than the *optical asymmetry ratio* at normal incidence, *i.e*. $\frac{{T}_{1}-{T}_{-1}}{{T}_{1}}$.

## 5. Concluding remarks

A novel FEM formulation was adapted to the analysis of z-anisotropic gratings relying on a rigorous treatment of the plane wave sources problem through an equivalent radiation problem with localized sources. The developed approach presents the advantage of being very general in the sense that it is applicable to every conceivable grating geometry. Furthermore, it can be easily extended to the case of an anisotropic grating in a multilayered dielectric stack.

Numerical experiments based on existing materials at normal and oblique incidences in both TE and TM cases showed the efficiency and the accuracy of our method. We demonstrated we could generate imbalanced symmetric propagative orders in the TE polarization case and at normal incidence with an aragonite grating on a silica substratum. A deeper research upon magneto-optic (electro-optic) materials could for instance lead to the conception of a controllable polariser since their dielectric tensor can be modified with a external constant magnetic (electric) field. Finally, the 3D case relies on the implementation of 3D PML and bi-pseudoperiodic conditions but the same approach concerning the source localization problem can be implemented. Studies are in progress in this direction.

## Acknowledgments

The authors gratefully acknowledge the Conseil Général des Bouches-du-Rhône for financial support. They also would like to thank the ANR MAGNETO-PHOT (BLAN06-2-135594) who supported this research.

## References and links

**1. **H. Kikuta, H. Toyotai, and W. Yul, “Optical elements with subwavelength structured surfaces,” Opt. Rev. **10**, 63–73 (2003) [CrossRef]

**2. **E. N. Glytsis and T. K. Gaylord, “High-spatial-frequency binary and multilevel stairstep gratings: Polarization-selective mirrors and broadband antireflection surfaces,” Appl. Opt. **31**, 4459–4470 (1992). [CrossRef] [PubMed]

**3. **S. M. Schultz, E. N. Glytsis, and T. K. Gaylord, “Design of a high-efficiency volume grating coupler for line focusing,” Appl. Opt. **37**, 2278–2287 (1998). [CrossRef]

**4. **K. Knop, “Diffraction gratings for color filtering in the zero diffraction order (ET),” Appl. Opt. **17**, 3598–3603 (1978). [CrossRef] [PubMed]

**5. **M. G. Moharam and T. K. Gaylord, “Rigorous coupled-wave analysis of planar-grating diffraction,” J. Opt. Soc. Am. **71**, 811–818 (1981). [CrossRef]

**6. **K. Rokushima and J. Yamanaka, “Analysis of anisotropic dielectric gratings,” J. Opt. Soc. Am. **73**, 901–908 (1983). [CrossRef]

**7. **E. N. Glytsis and T. K. Gaylord, “Rigorous three-dimensionnal coupled-wave diffraction analysis of single and cascaded anisotropic gratings,” J. Opt. Soc. Am. A **4**, 2061–2080 (1987). [CrossRef]

**8. **J. Chandezon, D. Maystre, and G. Raoult, “A new theoretical method for diffraction gratings and its numerical application,” J. Opt. (Paris) **11**, 235241 (1980). [CrossRef]

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

**10. **K. Watanabe, “Numerical integration schemes used on the differential theory for anisotropic gratings,” J. Opt. Soc. Am. A **19**, 2245–2252 (2002). [CrossRef]

**11. **D. Maystre, “A new general integral theory for dielectric coated gratings,” J. Opt. Soc. Am. **68**, 490–495 (1978). [CrossRef]

**12. **V. V. Belyaev, E. M. Kushnir, A. V. Klyshkov, and V. I. Tsoi, “Numerical modeling of the diffraction of light at periodic anisotropic gratings with rectangular surface microrelief,” J. Opt. Technol. **72**, 725–728 (2005). [CrossRef]

**13. **K. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media,” IEEE Trans. Antennas Propag. **AP-14**, 302–307 (1966).

**14. **W. M. Saj, “FDTD simulations of 2D plasmon waveguide on silver nanorods in hexagonal lattice,” Opt. Express **13**, 4818–4827 (2005). [CrossRef] [PubMed]

**15. **T. Delort and D. Maystre, “Finite-element method for gratings,” J. Opt. Soc. Am. A **10**, 2592–2601 (1993). [CrossRef]

**16. **Y. Ohkawa, Y. Tsuji, and M. Koshiba, “Analysis of anisotropic dielectric grating diffraction using the finiteelement method,” J. Opt. Soc. Am. A **13**, 1006–1012 (1996). [CrossRef]

**17. **Gang Bao, Zhiming Chen, and Haijun Wu, “Adaptive finite-element method for diffraction gratings,” J. Opt. Soc. Am. A **22**, 1106–1114 (2005). [CrossRef]

**18. **F. Zolla and R. Petit, “Method of fictitious sources as applied to the elctromagnetic diffraction of a plane wave by a grating in conical mounts,” J. Opt. Soc. Am. A **13**, 796–802 (1996). [CrossRef]

**19. **J.-P. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves,” J. Comput. Phys. **114**, 185–200 (1994). [CrossRef]

**20. **G. Tayeb, “Contribution à létude de la diffraction des ondes électromagnétiques par des réseaux. Réflexions sur les méthodes existantes et sur leur extension aux milieux anisotropes, ‘Université Paul Cézanne, PhD Thesis’ 1990.

**21. **N. Kono and Y. Tsuji, “A novel finite-element method for nonreciprocal magneto-photonic crystal waveguides,” J. Lightwave Technol. **22**, 1741–1747 (2004). [CrossRef]

**22. **A. F. Zhou, J. K. Erwin, C. F. Brucker, and M. Mansuripur, “Dielectric tensor characterization for magnetooptical recording media,” Appl. Opt. **31**, 6280–6286 (2005). [CrossRef]

**23. **G. Granet, “Reformulation of the lamellar grating problem through the concept of adaptive spatial resolution,” J. Opt. Soc. Am. A **16**, 2510–2516 (1999). [CrossRef]

**24. **R. Petit, “Electromagnetic Theory of Gratings,” Springer Verlag (1980).

**25. **F. Zolla, G. Renversez, A. Nicolet, B. Khulmey, S. Guenneau, and D. Felbacq, *Foundations of Photonic Crystal Fibres* (Imperial College Press, London, 2005). [CrossRef]

**26. **A. Nicolet, S. Guenneau, C. Geuzaine, and F. Zolla, “Modeling of electromagnetic waves in periodic media with finite elements,” J. Comput. Appl. Math. **168 (1–2)**, 321–329 (2004). [CrossRef]

**27. **Y. Ould Agha, F. Zolla, A. Nicolet, and S. Guenneau, “On the use of PML for the computation of leaky modes : an application to Microstructured Optcal Fibres,” The International Journal for Computation and Mathematics in Electrical and Electronic Engineering 27, no. 1, 95–109 (Jan.2008)

**28. **A. Nicolet, F. Zolla, Y. Ould Agha, and S. Guenneau, “Leaky modes in twisted microstructured optical fibres,” Waves Random Media, 17:4, 559–570 (2007).

**29. **M. Lassas, J. Liukkonen, and E. Somersalo, “Complex riemannian metric and absorbing boundary condition,” J. Math Pures Appl. **80**, 739–768 (2001).

**30. **M. Lassas and E. Somersalo, “Analysis of the PML equations in general convex geometry,” Proc. R. Soc. Edinburgh **131**, 1183–1207 (2001). [CrossRef]

**31. **G. Bao and H. Wu, “On the convergence of the solutions of PML equations for time harmonic Maxwell’s equations,” SIAM J. Numer. Anal. **43**, 2121–2143 (2005). [CrossRef]

**32. **P. Helluy, S. Maire, and P. Ravel, “Intégrations numériques d’ordre élevé de fonctions régulières ou singulières sur un intervalle,” CR. Acad. Sci. Paris, Sér. I, Math , **327**, 843–848 (1998).